MNRAS 000, 1–10 (2022) Preprint 9 January 2023 Compiled using MNRAS LATEX style file v3.

The ultraviolet continuum slopes (𝛽) of galaxies at z ' 8 − 16 from JWST

and ground-based near-infrared imaging
Fergus Cullen 1★ , R. J. McLure1 , D. J. McLeod1 , J. S. Dunlop1 , C. T. Donnan1 , A. C. Carnall1 ,
R. A. A. Bowler2 , R. Begley1 , M. L. Hamadouche1 , T. M. Stanton1
1 Institute for Astronomy, University of Edinburgh, Royal Observatory, Edinburgh, EH9 3HJ, UK
2 Jodrell Bank Centre for Astrophysics, University of Manchester, Oxford Road, Manchester, M13 9PL, UK

Accepted XXX. Received YYY; in original form ZZZ

arXiv:2208.04914v4 [astro-ph.GA] 6 Jan 2023

ABSTRACT
We study the rest-frame ultraviolet (UV) continuum slopes (𝛽) of galaxies at redshifts 8 < 𝑧 < 16 (h𝑧i = 10), using a
combination of James Webb Space Telescope (JWST) ERO and ERS NIRCam imaging and ground-based near-infrared imaging
of the COSMOS field. The combination of JWST and ground-based imaging provides a wide baseline in both redshift and
absolute UV magnitude (−22.6 < 𝑀UV < −17.9), sufficient to allow a meaningful comparison to previous results at lower
redshift. Using a power-law fitting technique, we find that our full sample (median 𝑀UV = −19.3 ± 1.3) returns an inverse-
variance weighted mean value of h𝛽i = −2.10 ± 0.05, with a corresponding median value of 𝛽 = −2.29 ± 0.09. These values
imply that the UV colours of galaxies at 𝑧 > 8 are, on average, no bluer than the bluest galaxies in the local Universe (e.g.,
NGC 1705; 𝛽 = −2.46). We find evidence for a 𝛽 − 𝑀UV relation, such that brighter UV galaxies display redder UV slopes
(d𝛽/dMUV = −0.17 ± 0.05). Comparing to results at lower redshift, we find that the slope of our 𝛽 − 𝑀UV relation is consistent
with the slope observed at 𝑧 ' 5 and that, at a given 𝑀UV , our 8 < 𝑧 < 16 galaxies are bluer than their 𝑧 ' 5 counterparts,
with an inverse-variance weighted mean offset of hΔ𝛽i = −0.38 ± 0.09. We do not find strong evidence that any objects in
our sample display ultra-blue UV continuum slopes (i.e., 𝛽 . −3) that would require their UV emission to be dominated by
ultra-young, dust-free stellar populations with high Lyman-continuum escape fractions. Comparing our results to the predictions
of theoretical galaxy formation models, we find that the galaxies in our sample are consistent with the young, metal-poor and
moderately dust-reddened galaxies expected at 𝑧 > 8.
Key words: galaxies: evolution - galaxies: formation - galaxies: high-redshift - galaxies: starburst - dark ages, reionization, first
stars

1 INTRODUCTION would therefore have important implications for our understanding

of the first galaxies and the process of cosmic hydrogen reionization.
Constraining the physical properties of the first galaxies is a key goal
of the James Webb Space Telescope (JWST). Within the first month In the pre-JWST era, no strong evidence for such primordial
of the data being released, JWST has already revealed a substantial 𝛽 = −3 populations was found. Studies of faint galaxies up to 𝑧 ' 7−8
population of previously unseen galaxies at 𝑧 > 10 (e.g., Adams with the Hubble Space Telescope (HST) revealed an average power-
et al. 2023; Atek et al. 2022; Castellano et al. 2022; Naidu et al. law index of h𝛽i ' −2, indicating moderately young and metal-poor,
2022; Donnan et al. 2022; Finkelstein et al. 2022; Harikane et al. but in no sense extreme, stellar populations (e.g., Dunlop et al. 2013;
2022), as well as one galaxy candidate at 𝑧 ' 16 (within ≈ 200 Finkelstein et al. 2012; Bouwens et al. 2014). Indeed, 𝛽 ' −2 to
Myr of the Big Bang; Donnan et al. 2022). These galaxies provide 𝛽 ' −2.5 is typical of the bluest galaxies observed at 𝑧 = 2 − 4
an unprecedented opportunity to study the properties of primordial (e.g., McLure et al. 2018), and even in the local Universe (e.g., NGC
stellar populations in the early Universe. 1705, 𝛽 = −2.46 ± 0.01, 𝑀UV = −18; Calzetti et al. 1994; Vázquez
One potential indicator of ultra-young, ultra-low metallicity, stel- et al. 2004). Early claims of extremely blue (i.e., 𝛽 ≤ −3) galaxies
lar populations is the power-law index of the rest-frame ultraviolet from HST imaging were later shown to be the result of an observa-
(UV) continuum, 𝛽, where 𝑓𝜆 ∝ 𝜆 𝛽 . At very young ages and low tional bias, pushing measurements towards artificially blue 𝛽 values
metallicities (e.g., 𝑡 < 30 Myr and 𝑍★ . 10−3 ), and in the absence for faint sources near the detection threshold (Bouwens et al. 2010;
of dust extinction and nebular continuum emission, a very low (i.e., Dunlop et al. 2012; Rogers et al. 2013).
blue) value of 𝛽 ' −3 is expected (e.g., Schaerer 2002; Bouwens et al. By providing unprecedentedly deep infrared imaging up to
2010; Chisholm et al. 2022). A robust determination of 𝛽 = −3 would 𝜆 = 5 𝜇m, JWST/NIRCam now enables the first robust estimates of
unequivocally indicate a stellar population that has recently formed 𝛽 for galaxies at 𝑧 > 8. Recently, Topping et al. (2022) have pro-
from pristine (or near-pristine) gas with a large ionizing photon es- vided the first measurements of 𝛽 at 𝑧 ' 7 − 11 from JWST/NIRCam
cape fraction (e.g., Robertson et al. 2010). Finding such galaxies imaging in the EGS field. They find a median value of 𝛽 = −2.0,
consistent with the blue, but otherwise unremarkable, populations
found with HST. Interestingly, however, Topping et al. (2022) also
report two galaxies with seemingly secure 𝛽 ' −3 measurements.
★ E-mail: [email protected] The ultra-blue UV slopes inferred for these sources are bolstered by a

© 2022 The Authors

2 F. Cullen et al.
lack of strong nebular emission-line signatures in the rest-frame opti- Table 1. The best-fitting UV continuum slopes (𝛽) for the full sample of
cal photometry, indicating the large ionizing photon escape fractions galaxies at 𝑧 > 7.5 in the combined JWST and COSMOS/UltraVISTA sam-
expected for such a population. ples. The first column gives the source ID taken from Donnan et al. (2022).
In this paper we use the new galaxy sample described in Donnan Column two gives the sample (COSMOS/UltraVISTA or JWST). Columns
et al. (2022) to present a complementary study of UV continuum three and four give the photometric redshift (𝑧phot ) and absolute UV magni-
slopes for N = 61 galaxies in the redshift range 𝑧 ' 8 − 16 (with tude (𝑀UV ) taken from Donnan et al. (2022). Column four gives the derived
mean h𝑧i ' 10). The Donnan et al. (2022) sample combines galax- UV continuum slope 𝛽.
ies drawn from the early JWST deep fields (at 𝑧 ' 9 − 16) with
an additional sample selected from wide-area ground-based near-IR ID Sample 𝑧phot 𝑀UV 𝛽
imaging in the COSMOS/UltraVISTA field (at 𝑧 ' 8 − 10). Cru-
334330 COSMOS 7.58 −21.30 +0.62
−1.01−0.64
cially, with the inclusion of this ground-based sample, we can also
−21.57 +0.74
−1.87−0.79
probe the brightest galaxies at these redshifts, which evade current 733875 COSMOS 7.58
JWST surveys. This extended baseline in UV luminosity enables us 812867 COSMOS 7.58 −21.02 +0.62
−2.19−0.72
to investigate 𝛽 across a factor of ' 80 in UV luminosity, placing +0.38
688541 COSMOS 7.66 −22.15 −2.88−0.41
early constraints on the relationship between 𝛽 and UV magnitude at +0.27
𝑧 > 8 (i.e., the 𝛽 − 𝑀UV , or the colour-magnitude, relation; Rogers 765906 COSMOS 7.66 −22.61 −1.26−0.25
626972 COSMOS 7.75 −21.49 +0.71
−2.65−1.03
et al. 2014).
Our aim is to provide an exploratory study of the constraints on 𝛽 536767 COSMOS 8.02 −21.40 +0.84
−2.16−0.88
at 𝑧 > 8 enabled by deep JWST multi-band imaging, and to critically 861605 COSMOS 8.02 −21.33 +1.10
−3.91−1.60
assess any early evidence for an evolution in the typical 𝛽 values, +0.75
978389 COSMOS 8.02 −21.68 −1.95−1.12
as well as the relation between 𝛽 and 𝑀UV at these redshifts. We
−22.05 +0.48
−2.07−0.50
also examine evidence for any robust 𝛽 ' −3 sources in our sample, 484075 COSMOS 8.11
and discuss the possibility of spurious 𝛽 ≤ −3 detections for faint 578163 COSMOS 8.20 −22.35 +0.31
−1.04−0.29
sources in the new JWST imaging. 458445 COSMOS 8.38 −21.65 +0.88
−2.89−1.11
The paper is structured as follows. In Section 2 we describe the +0.67
448864 COSMOS 8.57 −21.15 −2.29−0.81
data and galaxy sample constructed by Donnan et al. (2022), and
306122 COSMOS 8.76 −21.76 +0.49
−2.25−0.53
provide the details of our method for determining 𝛽. In Section 3
892014 COSMOS 8.96 −22.16 +0.39
−1.97−0.45
we present our 𝛽 measurements and outline the main results of our
analysis. In Section 4 we discuss the implications for our results 817482 COSMOS 9.89 −22.57 +0.56
−1.15−0.59
before summarising our main conclusions in Section 5. Throughout 43031 JWST 8.57 −18.43 +0.39
−2.17−0.41
we use the AB magnitude system (Oke 1974; Oke & Gunn 1983), +1.28
−18.41 −2.55−1.15
and assume a standard cosmological model with 𝐻0 = 70 km s−1
29274_4 JWST 8.86
+0.45
Mpc−1 , Ω𝑚 = 0.3 and ΩΛ = 0.7. 1434_2 JWST 9.16 −18.82 −2.28−0.45
44085 JWST 9.26 −18.25 +0.41
−1.41−0.38
38697 JWST 9.36 −18.86 +0.50
−1.77−0.46
5071 JWST 9.47 −18.02 +1.00
−2.55−1.01
2 DATA AND UV CONTINUUM SLOPE FITTING
44711 JWST 9.47 −20.14 +0.16
−2.12−0.15
2.1 JWST NIRCam imaging 43866 JWST 9.47 −18.14 +0.32
−2.78−0.31
Our JWST sample was initially presented in Donnan et al. (2022). 34086 JWST 9.47 −17.87 +0.25
−1.92−0.23
The sample is drawn from public NIRCam imaging of three fields 14391 JWST 9.47 −18.81 +0.51
−1.80−0.44
(SMACS J0723, GLASS and CEERS) released as part of the Early +0.50
12682 JWST 9.57 −18.95 −1.38−0.56
Release Observations (ERO, see Pontoppidan et al. 2022) and Early
44566 JWST 9.68 −20.68 +0.11
−1.64−0.12
Release Science (ERS) programmes (Treu et al. 2022). Each of the
22480 JWST 9.68 −18.50 +0.58
−1.60−0.52
three JWST fields were imaged in a combination of the F090W,
F115W, F150W, F200W, F277W, F356W, F410M and F444W fil- 15019 JWST 9.68 −18.67 +1.16
−4.61−1.96
ters, with the specific combination of filters varying slightly from 12218 JWST 9.68 −19.28 +0.27
−2.47−0.28
field to field (see Table 2 of Donnan et al. 2022). This JWST NIRCam +1.37
3398 JWST 9.68 −18.21 −5.72−2.18
imaging was reduced using PENCIL (PRIMER enhanced NIRCam
6200 JWST 9.79 −18.52 +0.88
−3.03−1.27
Image Processing Library) which is a custom version of the JWST
pipeline (1.6.2) with additional steps for background subtraction and 7606 JWST 9.89 −18.08 +2.21
−4.28−2.90
the removal of ‘snowball’ artefacts and including up-to-date calibra- 3763 JWST 9.89 −18.99 +0.39
−3.41−0.38
tions and zero-point corrections (see Donnan et al. 2022). The final
1698 JWST 10.45 −20.62 −2.00+0.14
combined JWST NIRCam imaging area totalled ' 45 arcmin2 (with −0.15

some variation between filters). For this work, prior to catalogue 20976_4 JWST 10.45 −18.80 −1.53+0.58
−0.61
construction, all of the NIRCam imaging was homogenized to the 6647 JWST 10.45 −18.88 −0.23+1.12
−0.92
point-spread-function (PSF) of the F444W filter. 3710 JWST 10.45 −19.06 +0.54
−2.05−0.53
The JWST catalogues were created by running Source Extrac- +0.74
4063 JWST 10.45 −18.03 −3.13−0.77
tor (Bertin & Arnouts 1996) in dual-image mode. The F200W image +0.38
was used as the detection image to optimise for the selection 𝑧 ≥ 8 30585 JWST 10.56 −19.35 −2.98−0.39
galaxies. The photometry for each JWST target was computed in both
0.5-arcsec and 0.36-arcsec diameter apertures. For the purposes of

MNRAS 000, 1–10 (2022)

 UV continuum slopes at 𝑧 ' 8 − 16 3

Table 1. Continued. EAZY and, after applying a robust selection criteria, we retained a final
sample of 16 galaxies at 𝑧 > 7.5. Absolute rest-frame UV magni-
ID Sample 𝑧phot 𝑀UV 𝛽 tudes were calculated for each object from the best-fitting EAZY SED.
Combined, our JWST and COSMOS/UltraVISTA samples yielded a
73150 JWST 10.56 −19.07 −3.57+0.95 total of 61 galaxies at 𝑧 ' 8 − 16.
−1.00
21071_2 JWST 10.68 −19.27 −2.71+0.64
−0.58
20757 JWST 10.68 −17.88 −0.74+1.02
−1.05 2.3 Measuring the UV continuum slope
6415 JWST 10.79 −19.13 −2.02+1.02
−1.05 A number approaches to determining the UV continuum slope from
120880 JWST 10.79 −19.43 −2.73+0.58
−0.59 broadband photometry have been presented in the literature, includ-
26598 JWST 10.79 −18.47 −3.31+0.80
−0.87 ing single colour measurements (e.g., McLure et al. 2011; Dunlop
61486 JWST 11.15 −19.61 −2.61+0.41
−0.52
et al. 2012, 2013) and SED template fitting (e.g., Finkelstein et al.
2012; Tacchella et al. 2022). Here we have adopted the power-law
622_4 JWST 11.27 −18.92 −3.38+0.56
−0.58 fitting method advocated by Rogers et al. (2014) in their study of
33593_2 JWST 11.27 −19.58 −2.07+0.28
−0.30 the 𝛽 − 𝑀UV relation at 𝑧 ' 5. For each source, the redshift was
77241 JWST 11.27 −19.60 −2.51+0.38
−0.42 fixed to the best-fitting photometric redshift estimated by Donnan
5268_2 JWST 11.40 −19.16 −2.41+0.52
−0.55
et al. (2022) and the photometry covering rest-frame wavelengths
127682 JWST 11.40 −19.07 −2.73+0.60 𝜆rest ≤ 3000 Å was modelled as a pure power law ( 𝑓𝜆 ∝ 𝜆 𝛽 ), with
−0.64
IGM absorption at 𝜆 ≤ 1216 Å included using the Inoue et al. (2014)
26409_4 JWST 11.90 −18.84 −3.25+1.02
−1.47 prescription. The only free parameter in this approach is 𝛽, the power-
8347 JWST 11.90 −19.09 −2.93+0.33
−0.38 law spectral index of the UV continuum red-ward of 𝜆 = 1216 Å. We
10566 JWST 12.03 −19.70 −3.44+0.43
−0.43
allowed 𝛽 to vary over the range −10 ≤ 𝛽 ≤ 10 and used the nested
32395_2 JWST 12.29 −19.89 −3.30+0.25 sampling code dynesty (Speagle 2020) to sample the full posterior
−0.30
distribution assuming a uniform prior. The derived values of 𝛽 for
1566 JWST 12.29 −18.77 −2.51+0.51
−0.55 our full sample are given in Table 1.
17487 JWST 12.42 −20.89 −2.64+0.26
−0.27 We investigated the effect of redshift uncertainties by running
27535_4 JWST 12.56 −19.42 −1.70+0.44
−0.46 an additional set of fits in which redshift (𝑧) was included as a
93316 JWST 16.39 −21.66 −1.89+0.15
−0.15
free parameter. For the prior on redshift we assumed a Gaussian
centered on the best-fitting photometric redshift from Donnan et al.
(2022) (𝑧phot ) with 𝜎𝑧 = 1. We find that the effect of fitting for
the present paper, we adopt the 0.5-arcsec apertures to prevent biases redshift on the derived values of 𝛽 is negligible, with a median
in 𝛽 measurements in more extended sources (Rogers et al. 2014). difference across the sample of Δ𝛽 = 0.04, corresponding to a median
However, we have confirmed that adopting the 0.36-arcsec diameter difference in redshift of (𝑧 − 𝑧 phot ) = −0.03. However, we find that
apertures would not change our main results. Redshifts for each ob- marginalizing over a plausible range of redshifts in this way increases
ject were estimated using the photometric redshift fitting code EAZY the typical error on 𝛽 by ' 13%. For the the purposes of this paper,
(Brammer et al. 2008). A thorough selection procedure, described we decided to fix our redshifts to the more-accurate 𝑧phot values
in Donnan et al. (2022), resulted in a final sample of 45 galaxies at presented in Donnan et al. (2022) (i.e., which are derived from fitting
𝑧 > 8.5 across the three fields. to the full rest-frame UV to optical photometry) but increased the
Absolute rest-frame UV magnitudes (𝑀UV ) were calculated for corresponding 𝛽 uncertainties by a factor 1.13.
each object by integrating the best-fitting EAZY spectral energy dis- We note that our approach is similar to the power-law fitting
tribution (SED) through a tophat filter centered on 𝜆rest = 1500 Å method used by Topping et al. (2022), with the main difference being
(Donnan et al. 2022). that our IGM model enables us to include filters encompassing the
Lyman break. However, our results are unchanged if we restrict the
fitting to filters red-ward of 1216 Å. Finally, it is also worth noting
2.2 COSMOS UltraVISTA that we have explicitly assumed that any emission lines present in the
UV spectrum - in particular Ly𝛼 - have a negligible impact on the
Our COSMOS sample was also initially presented in Donnan et al. observed photometry. This assumption appears to be justified based
(2022). The sample was drawn from the UltraVISTA survey (Mc- on the early, low-resolution, JWST/NIRSpec spectroscopy of 𝑧 > 9
Cracken et al. 2012) which provides deep 𝑌 𝐽𝐻𝐾 𝑠 near-IR imaging sources (e.g. Curtis-Lake et al. 2022; Roberts-Borsani et al. 2022).
across 1.8 deg2 in the COSMOS field. The deep near-IR imaging is
supplemented with optical imaging in 𝑢 ∗ 𝑔𝑟𝑖𝑧 from the CFHT Legacy
Survey (Hudelot et al. 2012), and the 𝐺 𝑅𝐼 𝑍 𝑦+NB816+NB921 fil-
3 RESULTS
ters from the Hyper Suprime-Cam Subaru Strategic Program (HSC-
SSP) DR2 (Aihara et al. 2019). All of the near-infrared and optical In Fig. 1 we plot the 𝛽 values for our full JWST and COS-
imaging in COSMOS was aligned to the GAIA EDR3 reference MOS/UltraVISTA sample versus redshift, 𝑧, and absolute UV mag-
and PSF-homogenised to the UltraVISTA 𝑌 −band. Additionally, the nitude, 𝑀UV . As well as illustrating the typical 𝛽 values in our
COSMOS/UltraVISTA dataset was further augmented by 3.6𝜇m and sample, the plots clearly demonstrate the power of combining JWST
4.5𝜇m photometry from Spitzer/IRAC imaging provided by the Cos- with ground-based surveys to probe a large dynamic range in both
mic Dawn Survey (Euclid Collaboration et al. 2022). 𝑧 and, particularly, 𝑀UV . The first point to note is the large scatter
The COSMOS catalogue was produced from inverse variance in observed 𝛽 values, which increases towards the faint luminosity
weighted stacks of the data in the 𝑌 , 𝐽, 𝐻, and 𝐾 𝑠 −bands as described limit in both samples. This effect is seen most dramatically for the
in Donnan et al. (2022). Photometric redshifts were estimated with JWST sample, where values as extreme as 𝛽 < −4 are recovered at

MNRAS 000, 1–10 (2022)

4 F. Cullen et al.

Table 2. Average 𝛽 values and standard errors derived for our full sample 1
and in two bins of absolute UV magnitude. The first column defines each JWST
sample in terms of 𝑀UV . In the second column we report the inverse-variance COSMOS UltraVISTA
weighted mean and standard error of the individual 𝛽 values. In the third 0
column we report the median and 𝜎MAD of the individual 𝑀UV values, where
𝜎MAD = 1.483 × MAD and MAD refers to the median absolute deviation.
−1
Sample h𝛽 i h𝑀UV i
−2

β
Full sample (all 𝑀UV ) −2.10 ± 0.05 −19.3 ± 1.3
𝑀UV ≤ −20.5 −1.80 ± 0.08 −21.6 ± 0.6
−3
𝑀UV > −20.5 −2.32 ± 0.07 −18.9 ± 0.6

−4
𝑀UV & −19. However, the large error bars at these faint luminosi-
ties (𝜎𝛽 ' 1 at 𝑀UV > −19) suggests that this is predominantly
a result of observational uncertainties. Indeed, the preference for −5
low-luminosity galaxies to be scattered blue is a well-known effect,
8 10 12 14 16
caused by the fact that if a galaxy’s flux is boosted into the detec-
Redshift
tion band it will always be biased towards bluer UV slopes (Dunlop
et al. 2013; Rogers et al. 2013). At the brightest UV luminosities in
the JWST sample, where the constraints on individual 𝛽 estimates
improve significantly, the scatter noticeably reduces and fewer ultra- 1
blue (𝛽 < −3) objects are seen. We will discuss implications for the
detection reliability of ultra-blue objects at faint luminosities in more 0
detail in Sections 3.2 and 4.1.
In the 𝛽 versus 𝑧 plot an increase in the scatter at 𝑧 ' 9 − 11
is apparent for the JWST sample. This is caused by a combination −1
of (i) a larger number of intrinsically faint galaxies being detected
in this redshift range, and (ii) a minimum in the number of filters
−2
covering rest-frame wavelengths 𝜆rest ≤ 3000 Å (typically Nfilt = 3
β

at 𝑧 . 10.5 versus Nfilt = 4 at 𝑧 & 10.5, depending on the field).

Interestingly, one of our most robust 𝛽 estimates is the putative 𝑧 ' −3
16.4 galaxy candidate (CEERS 93316) reported in Donnan et al.
(2022), which has 𝛽 = −1.9±0.15. This tight constraint is due, in part,
to the excellent sampling of the rest-frame UV slope for this galaxy −4
(it is covered by the F277W, F356W, F410M and F444W filters).
Promisingly, if these extremely high-redshift objects are confirmed - −5
and if more are uncovered - JWST will be able to accurately constrain
their UV continuum slopes thanks to the excellent coverage of the −23 −22 −21 −20 −19 −18
rest-frame UV continuum at 𝑧 > 11. MUV
It can be seen from Fig. 1 that the typical values of 𝛽 displayed
by the galaxies in our 𝑧 ' 8 − 16 sample are somewhat bluer, but Figure 1. Plots of UV continuum slope 𝛽 versus redshift (top) and versus
not obviously more extreme, than the typical values found at 𝑧 ≤ 8 absolute UV magnitude 𝑀UV (bottom) for the galaxies in our JWST (blue)
with HST (i.e., 𝛽 ' −2; Dunlop et al. 2013). The distribution of and COSMOS/UltraVISTA (red) sub-samples.
points in Fig. 1 is consistent with being drawn from an underlying
population with a relatively narrow intrinsic distribution of 𝛽, with
some evidence for a shallow trend towards bluer 𝛽 values at fainter in the typical values of 𝛽 at 𝑧 > 8. In fact, these early results imply
𝑀UV . In Table 2 we report the inverse-variance weighted mean 𝛽 that even the faintest galaxies that JWST has so far uncovered at
value for our full sample, which we find to be h𝛽i = −2.10 ± 0.05. 𝑧 ' 8 − 16 have, on average, UV colours no more extreme than the
In this instance we preferred the weighted mean over the median so bluest galaxies in the local Universe (e.g., NGC 1705; 𝛽 = −2.46,
as not to be biased by the blue-scatter effect at faint luminosities 𝑀UV = −18).
(i.e., the blue-scattered galaxies are not down-weighted by their large
uncertainties when taking the median). Indeed, the median of the full
3.1 The 𝛽 − MUV relation at z > 8
sample is 𝛽 = −2.29 ± 0.09, where the uncertainty on the median is
estimated using the median absolution deviation estimator (𝜎MAD = The 𝛽 − 𝑀UV relation, often referred to as the colour-magnitude re-
1.483 × MAD). As expected, the median estimate is bluer, although lation, encodes information on the dust and stellar population prop-
the formal difference is only at the ' 2𝜎 level. Our sample average is erties of galaxies as a function of their absolute UV magnitude. A
in decent agreement with the median values reported at 𝑧 ' 7 − 11 in number of studies at 𝑧 . 8 have found strong evidence for a 𝛽 − 𝑀UV
Topping et al. (2022) (𝛽 = −2.0 at 𝑧 ' 7 and 𝛽 = −1.9 at 𝑧 ' 8 − 11). relation in which the UV continuum slopes of galaxies are bluer at
Adopting either the inverse-variance weighted mean or median, it fainter luminosities. This has been used to argue that UV-faint galax-
is clear that our sample shows no evidence for significant evolution ies are typically younger, less metal-enriched, and less dust-obscured

MNRAS 000, 1–10 (2022)

 UV continuum slopes at 𝑧 ' 8 − 16 5
To obtain a rough estimate of the intrinsic scatter in the relation, we
assumed the total variance in the data was a combination of the mea-
1 2 = 𝜎2 + 𝜎2 .
surement error (𝜎m ) and the intrinsic scatter (𝜎sc ): 𝜎tot m sc
2
We found that the value of 𝜎sc that yielded 𝜒 /𝜈 ' 1 was 𝜎sc ' 0.35.
0 Interestingly, this value is again in good agreement with the result of
Rogers et al. (2014), who estimated that the intrinsic scatter of the
𝛽 − 𝑀UV relation at 𝑧 ' 5 increases from 𝜎sc ' 0.1 at 𝑀UV = −18
−1 to 𝜎sc ' 0.4 at 𝑀UV = −21.
Our data also suggest that the normalisation of the relation has
−2 evolved such that, at higher redshifts, the typical 𝛽 values are bluer
β

across the full 𝑀UV range. At bright magnitudes (𝑀UV ≤ −20.5)

the inverse-variance weighted mean offset is hΔ𝛽i = −0.24 ± 0.15,
−3 increasing to hΔ𝛽i = −0.44 ± 0.11 at the faint end (𝑀UV > −20.5).
The offset averaged across the full sample is hΔ𝛽i = −0.38 ± 0.09.
Evidence for a signal in these early datasets is encouraging, and
−4
future larger-area JWST surveys such as PRIMER (GO 1837) will
clarify this situation in the near future. These upcoming surveys will
Rogers et al. 2014 (z = 5)
−5 Bouwens et al. 2014 (z = 5)
serve to both increase the sample size and fill the current magnitude
gap at −21 < 𝑀UV < −20 where JWST can deliver excellent 𝛽
−23 −22 −21 −20 −19 −18 constraints. Overall, the analysis of the 𝑀UV − 𝛽 relation further
MUV emphasises the main result of our analysis: although the galaxies at
𝑧 > 8 are generally bluer than their lower-redshift counterparts, on
Figure 2. A comparison between the 𝛽 versus 𝑀UV relation at 𝑧 > 8 and average, the UV colours of our 𝑧 ' 8 − 16 galaxy sample are not
previously-determined relations at lower redshift. The black solid line shows dramatically bluer than bluest stellar populations observed at lower
the best-fitting 𝛽 − 𝑀UV relation for our full sample which has a slope of redshift, including sources at 𝑧 = 0.
d𝛽/dMUV = −0.17±0.05. The light-grey shaded region represents the 68 per
cent confidence interval around our best-fitting relation. The large diamond
points are the inverse-variance weighted mean values of 𝛽 in the two bins of 3.2 Evidence for ultra-blue objects (𝛽 ' −3)?
absolute UV magnitude given in Table 2. The orange dashed and dotted lines
show the 𝑧 ' 5 relations from Rogers et al. (2014) and Bouwens et al. (2014) Although the typical UV slopes in our sample appear to be no bluer
which have d𝛽/dMUV = −0.12 ± 0.02 and d𝛽/dMUV = −0.14 ± 0.02, re- than the bluest galaxies observed locally, ultra-blue objects (i.e.,
spectively. 𝛽 ≤ −3) may still exist within the population. Indeed, Topping et al.
(2022) have recently identified two sources at 𝑧 ' 7 with reportedly
secure detections of 𝛽 ' −3 from their investigation of the early
than their brighter counterparts (e.g., Bouwens et al. 2014; Rogers
CEERS NIRCam imaging data. If confirmed, this would represent
et al. 2014).
intriguing evidence for young, low metallicity stellar populations
Given the large dynamic range in 𝑀UV provided by our combined
with ionizing continuum escape factions of ' 100% (e.g. Robertson
JWST and COSMOS/UltraVISTA sample, we can examine early
et al. 2010; Chisholm et al. 2022).
evidence for a 𝛽 − 𝑀UV relation at 𝑧 > 8. In Table 2, we report
Our initial JWST and COSMOS/UltraVISTA sample does not
the average 𝛽 values for our sample split into two magnitude bins
provide convincing evidence for such objects. It can be seen in Fig. 3
divided at 𝑀UV = −20.5. We find h𝛽i = −2.32 ± 0.07 for the faint
that the majority of the galaxies in our sample with formal best
bin (median 𝑀UV = −18.9) and h𝛽i = −1.80 ± 0.08 for the bright
fits of 𝛽 ≤ −3 have large uncertainties in the measurement of 𝛽.
bin (median 𝑀UV = −21.6). Our sample is therefore consistent with
In this case it is more likely that the galaxies have been scattered
previous studies, with an evolution to redder colours in brighter
to blue values due to the known blue-bias in the 𝛽 scatter at faint
galaxies.
luminosities (Dunlop et al. 2012; Rogers et al. 2013). We provide a
Interestingly, we find that the formal best-fitting slope 𝛽 − 𝑀UV re-
detailed discussion of this effect in Section 4.1. In contrast, galaxies
lation for our sample is fully consistent with relations derived at lower
with well-constrained UV slopes, which we here define as those with
redshift (within 1𝜎). In Fig. 2 we plot the 𝛽 − 𝑀UV relations at 𝑧 ' 5
an uncertainty of 𝜎𝛽 ≤ 0.22 show no evidence for slopes bluer than
from Rogers et al. (2014) and Bouwens et al. (2014), who both report
𝛽 ' −2.2. Overall, the results shown in Fig. 3 imply that the ultra-
a modest evolution in 𝛽 as a function of absolute UV magnitude,
blue values we see in our sample are a result of statistical, rather than
finding d𝛽/dMUV = −0.12 ± 0.02 and d𝛽/dMUV = −0.14 ± 0.02,
physical, effects.
respectively. Fitting a similar colour-magnitude relation to our indi-
The most plausible ultra-blue candidate in our sample is ID
vidual sources yields a best-fitting slope of d𝛽/dMUV = −0.17±0.05
32395_2 (𝑧 = 12.29; 𝑀UV = −19.89), which has a formal best-
(black solid line in Fig. 2)1 . The full best-fitting color-magnitude re- +0.25 , a value consistent
fitting UV continuum slope of 𝛽 = −3.30−0.30
lation given by
with the two ultra-blue candidates reported in Topping et al. (2022).
+0.05 +1.18
𝛽 = −0.17−0.05 𝑀UV − 5.40−1.25 . (1) The UV slope for this object is formally a 4.8𝜎 deviation from the

Formally, 𝜒2 /𝜈 = 1.5 for the best-fit model with respect to the data.
2 Although this definition is somewhat arbitrary, a galaxy with
𝛽 ≤ −3.0 ± 0.2 would represent a & 5𝜎 deviation from the sample aver-
1 Fitting to the two inverse-variance weighted mean values given in age of 𝛽 = −2.1, and a ' 3𝜎 deviation from the UV slope expected for dust
Table 2 returns a formally steeper, but fully consistent, value of free galaxy with a low escape fraction (e.g., 𝛽 ' −2.4; Cullen et al. 2017), and
d𝛽/dMUV = −0.20 ± 0.06. would thus be a strong candidate for an exotic ‘ultra-blue’ stellar population.

MNRAS 000, 1–10 (2022)

6 F. Cullen et al.
compared to 𝑧 = 5), but no bluer than the bluest objects uncovered
at lower redshifts. We do not find strong evidence for a significant
1
population of ultra-blue 𝛽 < −3 objects in our sample.
In this section we first provide a short discussion of our results,
0 starting with an exploration of the well-know faint-end blue bias and
its affect on our current sample. We then compare our results to
pre-JWST literature measurements at similar redshifts, as well as to
−1 predictions from a number of theoretical galaxy formation models.

−2 4.1 The blue 𝛽 bias at faint magnitudes

β

A bias towards bluer values of 𝛽 at faint magnitudes is a well-known

−3 phenomenon that has been extensively documented in 𝐻𝑆𝑇 studies
(e.g. Bouwens et al. 2010; Dunlop et al. 2012; Rogers et al. 2013).
The bias occurs due the that fact that high-redshift galaxy candidates
−4
are typically selected using a photometric filter as close as possible
to the Lyman break where the UV spectral energy distribution of
−5 young star-forming galaxies peaks. At faint magnitudes, this favours
the selection of objects whose photometry has been ‘up-scattered’
−23 −22 −21 −20 −19 −18 in the short-wavelength detection band; these objects will naturally
MUV appear bluer than they actually are.
To investigate the magnitude of this effect in our JWST sample we
Figure 3. Plot of 𝛽 versus 𝑀UV for the objects with well constrained UV ran a simple simulation. We first constructed 20,000 simple power-
continuum slopes (𝜎𝛽 < 0.2; black filled circles) and the objects with 𝛽 ≤ −3 law SEDs with a intrinsic UV slope of 𝛽int = −2.1 at the median
(grey open circles). We find no strong evidence for UV slopes as blue as redshift of the JWST galaxies (𝑧 = 10.5). The UV magnitudes were
𝛽 ≤ −3 amongst those galaxies with robust measurements of 𝛽. All of the drawn uniformly within the range −20.0 ≤ 𝑀UV ≤ −18.0 and IGM
objects with formal 𝛽 ≤ −3 solutions are poorly constrained and consistent
attenuation was applied using the Inoue et al. (2014) prescription.
with the known blue bias in the 𝛽 scatter near the source-detection threshold.
Photometry was generated in each of the JWST filters and scattered
according the the typical imaging depths (averaging the depths across
sample average, but on closer inspection we find that this galaxy multiple fields where appropriate). The ‘observed’ UV continuum
suffers from above-average systematic uncertainties. For example, slopes (𝛽obs ) were then recovered for the simulated galaxies using
when fitting for both 𝛽 and 𝑧 (Section 2.3) we find a large redshift the same method applied to the real observations.
offset from the Donnan et al. (2022) estimate, with Δ𝑧 = −0.18 (cf. The results of this simulation are shown in Fig. 4. Our results are
the sample median of Δ𝑧 = −0.03). This shift in redshift results in a consistent with the trends observed in previous works (e.g, Rogers
redder best-fitting UV continuum slope of 𝛽 = −3.1. Moreover, we et al. 2013). We find that, at faint 𝑀UV , the galaxies with the highest
find that ID 33593_2 is more sensitive to the chosen aperture size signal-to-noise ratio in the detection band (F200W) are biased blue. If
than the average galaxy in our sample. Adopting smaller 0.36 arcsec we mimic the selection of our sample (i.e., requiring a ≥ 5𝜎 detection
+0.25 . This 0.36 arcsec aperture value is
apertures yields 𝛽 = −2.99−0.30 in F200W) we find that 𝛽int is accurately recovered - on average -
consistent within 3.5𝜎 of the sample average, and within < 2𝜎 of for galaxies brighter than 𝑀UV ' −19.3, but becomes increasingly
the UV continuum slope expected for a standard stellar populations biased to blue values at fainter magnitudes (black dashed line in Fig.
with 𝑍 ' 0.1𝑍 and 𝑓esc = 0.0 (e.g., 𝛽 ' −2.4; Cullen et al. 2017). 4). Fitting for this average bias (Δ𝛽 = 𝛽obs − 𝛽int ) at 𝑀UV > −19.3
Although this galaxy is almost certainly one of the bluest objects we find that Δ𝛽 = −0.275𝑀UV − 5.304. Our results indicate that
in our sample, the combined statistical and systematic uncertainties galaxies at the faint end of our sample (i.e., 𝑀UV ' −18.5) will have
make it difficult to confirm it as a robust ultra-blue, 𝛽 < −3, object. Δ𝛽 = −0.2, on average.
Despite this, Fig. 3 does demonstrate that the currently-available Applying this average bias correction to the individual galaxies in
JWST imaging can undoubtedly deliver well-constrained 𝛽 measure- our sample does not strongly affect our main results. We find that the
ments (𝜎𝛽 ≤ 0.2) for galaxies as faint as 𝑀UV ' −20 at 𝑧 > 8, and recovered slope of the the 𝛽 − 𝑀UV relation becomes slightly shal-
hence should be able to uncover strong candidate 𝛽 ' −3 objects lower (d𝛽/dMUV = −0.10 ± 0.06; still consistent with the Rogers
at these absolute UV magnitudes should they exist (e.g., Topping et al. (2014) and Bouwens et al. (2014) slopes at 𝑧 = 5 within 1𝜎)
et al. 2022), despite the fact that no convincing candidates are found and the inverse-variance weighed mean of the sample slightly red-
here in our current high-redshift galaxy sample. In a future study, der (h𝛽i = −2.05 ± 0.05), but both remain fully consistent with
exploiting data from upcoming, wider-area, JWST Cycle-1 imaging the original non-corrected values. Nevertheless, as sample sizes in-
surveys, we intend to undertake a detailed analysis of the 𝛽 − 𝑀UV crease, and the statistical uncertainties are reduced, this effect will
relation at 𝑧 ≥ 7 and attempt to robustly quantify the intrinsic scatter clearly become more important, potentially requiring more sophis-
in the 𝛽 distribution. ticated simulations including source injection/recovery, treatment of
individual fields, and a consideration of aperture effects and redshift-
dependent systematics.
Fig. 5 also emphasises the fact that even in the magnitude regime
4 DISCUSSION
where the average properties are accurately recovered, the scatter of
We have presented the first estimate of the 𝛽 − 𝑀UV relation at 𝑧 > 8 any individual 𝛽obs can still be substantial. For example, at the mag-
using early JWST data. We find that, on average, galaxies at these nitude of our most robust ultra-blue candidate (𝑀UV ' −20; Section
redshift are bluer than their lower redshift counterparts (Δh𝛽i ' −0.4 3.2) the bias for individual objects in our simulation can be as large at

MNRAS 000, 1–10 (2022)

 UV continuum slopes at 𝑧 ' 8 − 16 7
is larger (redder) than predicted by our best-fitting 𝛽 − 𝑀UV relation,
Selection band SNR (F200W) but clearly consistent with our data within the substantial uncertainty
0 1 2 3 4 5 (Fig. 5). Some degree of systematic offset is perhaps unsurprising
given the limited rest-frame UV coverage of the Dunlop et al. (2013)
data. Technically, the Dunlop et al. (2013) measurement is an estimate
of 𝛽 in the far-ultraviolet, whereas our measurement (and all others
8 0.0
shown in Fig. 5) span far- and near-ultraviolet wavelengths. Indeed,
if we restrict our data to similar rest-frame wavelengths (. 1600Å)
∆β

−0.2
6 we recover a inverse-variance weighed mean of h𝛽i = −1.90 ± 0.10
−0.4
for our JWST sample. In that sense, the rest-frame FUV colors of
−0.6
4 −20 −19 −18 the Dunlop et al. (2013) sample are fully consistent with our new
MUV JWST sample. However, the addition of longer-wavelength anchors
2 at 𝜆 > 2000Å does favour a bluer h𝛽i.
∆β

Wilkins et al. (2016) presented measurements of 𝛽 for five liter-

0 ature sources at 9.6 < 𝑧 phot < 10.2. Their sample was comprised
of four galaxies drawn from the GOODS-South and GOODS-North
−2 field (Oesch et al. 2014), and one gravitatioanlly lensed source re-
ported in Zheng et al. (2012). The full sample spans the UV mag-
−4 nitude range −21.6 ≤ 𝑀UV ≤ −19.4. UV continuum slopes for
∆β = −0.275 × MUV − 5.304 these sources were determined using a single colour measurement
−6 (𝐻F160W − 3.6𝜇m) probing the rest-frame UV continuum in the
−20 −19 −18 range 1200Å . 𝜆 . 3700Å. The five Wilkins et al. (2016) sources
MUV are shown in the 𝛽 − 𝑀UV plane in the left-hand panel of Fig. 5
and are clearly in excellent agreement with our results. Their sample
Figure 4. Plot of the UV continuum slope bias (Δ𝛽 = 𝛽obs −𝛽int ) as a function average, estimated by stacking the photometry of the five individual
of 𝑀UV for of 20,000 simulated galaxies at 𝑧 = 10.5 with 𝛽int = −2.1 (see sources, is h𝛽i = −2.1±0.3. Restricting our sample to the same 𝑀UV
Section 4.1 for a description of the simulations). In the main panel, each data range as the Wilkins et al. (2016) sample returns an inverse-variance
point represents an individual galaxy and is colour-coded according to the weighted mean of h𝛽i = −2.08 ± 0.05.
SNR in the F200W selection filter (the colour bar is saturated at SNR = 5 Bhatawdekar & Conselice (2021) reported h𝛽i = −2.52+0.32
for clarity). The preference for galaxies with the highest SNR in the selection −0.20
for eight lensed 𝑧 = 9 galaxies in the Frontier Field cluster
band to skew towards bluer observed 𝛽 (i.e., Δ𝛽 < 0) is clearly visible at faint
magnitudes (i.e., the colour asymmetry around the Δ𝛽 = 0 line). The black
MACS J0416.1 − 2403. The median magnification-corrected abso-
dashed line shows the average bias for galaxies with SNR ≥ 5 in F200W (i.e., lute UV magnitude of their sample was 𝑀UV = −19.44. Stellar
mimicking our selection criteria) which tends to Δ𝛽 < 0 at 𝑀UV & −19. population models were fit to multi-wavelength photometry span-
The inset panel shows a zoom-in of this average bias, with the purple line ning 0.4 − 4.5𝜇m (i.e., rest-frame UV to optical wavelengths) and 𝛽
showing a linear fit to the relation. From our simple simulation, we find that was calculated from the best-fitting model using the Calzetti et al.
the systematic bias in 𝛽 occurs for galaxies fainter than 𝑀UV > −19.3, and (1994) windows. The Bhatawdekar & Conselice (2021) measure-
has the functional form Δ𝛽 = −0.275𝑀UV − 5.304. ment is clearly in good agreement with our data; our best-fitting
relation predicts h𝛽i = −2.4 at 𝑀UV = −19.44 (Fig. 5). Interest-
ingly, Bhatawdekar & Conselice (2021) also report a tentative color-
Δ𝛽 = ±0.5. Ultimately, while promising ultra-blue candidates can be
magnitude slope of d𝛽/dMUV = −0.19 ± 0.11 at 𝑧 = 9. While the
identified from data at these depths, robust confirmation will likely
uncertainty on this value clearly significant, the formal best-fitting
require either additional photometric tracers (e.g., a lack of emission
value is consistent our estimate.
line signatures in the rest-frame optical photometry; Topping et al.
2022), or deep spectroscopic follow-up. Finally, Tacchella et al. (2022) presented 𝛽 estimates for eleven
bright (𝑀UV < −20.7) galaxies at 𝑧 = 9 − 11, fitting to deep HST
and Spitzer/IRAC photometry and using a similar stellar population
4.2 Comparison with pre-JWST results at z ' 9 − 11 model-fitting technique to Bhatawdekar & Conselice (2021). Again,
Or results can be compared with a number of earlier works at these measurements are fully consistent with our data, although in this
𝑧 ' 9 − 11 undertaken prior to the launch of JWST. These studies case the comparison is primarily with our bright ground-based COS-
were based on either single-colour measurements with HST (using MOS/UltraVISTA sample (Fig. 5). The inverse variance weighed
the 𝐽140 and 𝐻160 filters), or a combination of HST and Spitzer/IRAC mean of the Tacchella et al. (2022) sample is h𝛽i = −1.63 ± 0.04
3.6𝜇m imaging. Most samples consisted of UV-bright galaxies above compared to h𝛽i = −1.79 ± 0.12 for our COSMOS sample. The un-
the knee of the galaxy luminosity function (i.e., 𝑀UV . −20.5; certainties on the Tacchella et al. (2022) estimates are clearly smaller
Bowler et al. 2020; Donnan et al. 2022). At fainter magnitudes, in this regime, primarily due to their deeper HST photometry. How-
galaxies were drawn from the ultra-deep HST imaging in the Hub- ever, as discussed in detail by Rogers et al. (2013), it is worth bearing
ble Ultra Deep Field (UDF-12; Dunlop et al. 2013) or gravitational in mind that the uncertainties resulting from a template fitting ap-
lensing fields. proach can be artificially reduced by the limited parameter space
Based on the UDF-12 dataset, Dunlop et al. (2013) provided the afforded by population synthesis models (which act as a prior on the
first tentative estimate of h𝛽i at 𝑧 = 9. The Dunlop et al. (2013) allowed values of 𝛽).
sample consisted of faint galaxies with 𝑀UV ' −18 and photometric Overall, the comparison with previous work at 𝑧 ' 9 − 11 is
measurements in the HST 𝐽140 and 𝐻160 filters. Using a single-colour clearly encouraging. Fig. 5 also highlights the power the new JWST
estimate, Dunlop et al. (2013) found h𝛽i = −1.80 ± 0.63. This value datasets, which are able to provide - for the first time - 𝛽 estimates

MNRAS 000, 1–10 (2022)

8 F. Cullen et al.

1 1

0 0

−1 −1

−2 −2
β

−3 −3

−4 Tacchella+(2022)
Wilkins+(2016)
−4
Vijayan+(2021) FLARES
Bhatawdekar+(2022) Kannan+(2022) THESAN
−5 Dunlop+(2013) −5 SC-SAM Ultra Wide

−23 −22 −21 −20 −19 −18 −23 −22 −21 −20 −19 −18
MUV MUV
Figure 5. A comparison of the 𝛽 − 𝑀UV relation at 𝑧 > 8 with pre-JWST literature measurements and theoretical galaxy formation models. In the left-hand
panel we show a comparison of our individual data (light grey) and best-fitting relation (black line) to the four pre-JWST studies at 𝑧 ' 9 − 10 indicated in the
legend (Dunlop et al. 2013; Wilkins et al. 2016; Bhatawdekar & Conselice 2021; Tacchella et al. 2022). In general, there is excellent agreement between our new
measurements and the earlier literature data. This comparison also emphasises the power of these new JWST datasets, which have substantially increased the
sample size at faint magnitude (𝑀UV & −20.5). In the right-hand panel we compare our data to three state-of-the-art simulations of galaxy formation: sc-sam
(Yung et al. 2019a), flares (Vijayan et al. 2021) and thesan (Kannan et al. 2022). As discussed in Sec. 4.3, the agreement is is qualitatively very good. The
simulated 𝛽 − 𝑀UV relations have a similar normalization to our data, as well as similar predicted slopes (d𝛽/dMUV ' 0.10 − 0.15). The purple 2D histogram
shows the distribution of ≈ 150, 000 galaxies from sc-sam to illustrate the intrinsic scatter. In this simulation, the predicted intrinsic scatter varies from 𝜎sc = 0.1
at the faint end to 𝜎sc = 0.3 at the bright end, again in reasonable agreement with our estimate of 𝜎sc = 0.35. The good agreement with theoretical predictions
suggests the galaxies in our sample are consistent with the young, metal-poor and moderately dust-reddened population predicted at 𝑧 > 8.

for individual galaxies at 𝑀UV & −20.5 without the assistance of range 0.01 − 1𝑍 . The median mass and metallicity of the simu-
gravitational lensing. lated sample is 𝑀★ ' 108.5 𝑀 and 𝑍★ ' 0.1𝑍 . The mean stellar
age of the sample is ' 100 Myr. We calculated an observed 𝛽 for
each galaxy (i.e., after the application of dust reddening) by fit-
4.3 Comparison with galaxy-formation model predictions ting a power-law to the noise-free mock photometry. The predicted
𝛽 − 𝑀UV relation (observed) is shown in the right-hand panel of Fig.
It is instructive to compare our results to the predictions of state-of- 5, where it can be seen that the overall normalization of the relation
the-are galaxy formation models. As our main comparison, we use is in good agreement with our data. The slope of the sc-sam relation
the Santa Cruz semi-analytic model (sc-sam) for galaxy formation (d𝛽/dMUV ' −0.15) is also fully consistent with our estimate. The
(Somerville et al. 2015, 2021). The sc-sam includes sophisticated 2D histogram in Fig. 5 shows the intrinsic scatter, which increases
prescriptions for cosmological accretion, gas cooling, star-formation, from 𝜎sc = 0.1 at 𝑀UV = −18 to 𝜎sc = 0.3 at 𝑀UV = −22, again in
chemical enrichment and stellar and AGN feedback, and has been reasonable agreement with our estimate of 𝜎sc = 0.35.
shown to successfully reproduce the global properties and scaling
relations of the high-redshift galaxy population out to 𝑧 = 10 (Yung Both the shape and scatter of the 𝛽 − 𝑀UV relation in the sc-sam
et al. 2019a,b). The star-formation and chemical enrichment histories are driven by UV dust attenuation (𝐴UV ). This is unsurprising, as
of the model galaxies are used to generate mock galaxy photometry the intrinsic low-order shape of the UV continuum is not strongly
based on the Bruzual & Charlot (2003) stellar population synthe- metallicity-dependent (Cullen et al. 2019), and at high redshifts the
sis models. Dust attenuation is applied assuming a Calzetti et al. of effect stellar population age is limited by the young age of the
(2000) attenuation law, with the rest-frame 𝑉-band dust attenuation Universe (Tacchella et al. 2022). Across the full range in 𝑀UV , the
calculated based on the surface density and metallicity of cold gas median UV attenuation increases from 𝐴UV = 0.05 to 𝐴UV = 0.64
(Somerville et al. 2012). We refer interested readers to Yung et al. (i.e., galaxies at the bright end in the sc-sam suffer a factor ' 2
(2022b) for a detailed description of the model, including a flowchart decrease in their intrinsic UV flux). The increase in scatter at the
illustrating the full internal workflow of the sc-sam. bright end is also driven by dust, with a larger range of 𝐴UV at
We obtained a sample of ≈ 150, 000 galaxies at 8 < 𝑧 < 10 bright magnitudes. At 𝑀UV = −22, the standard deviation of the 𝐴UV
with dust-attenuated absolute UV magnitudes in the range distribution is 𝜎𝐴UV = 0.4, compared to 𝜎𝐴UV = 0.1 at 𝑀UV = −18.
−23 < 𝑀UV < −18 from the sc-sam ultra-wide lightcone (cover- In the right-hand panel of Fig. 5 we also show 𝛽 − 𝑀UV predictions
ing 2 deg2 ; Yung et al. 2022a)3 . These galaxies have predicted stellar at 𝑧 = 9 from the flares (Lovell et al. 2021; Vijayan et al. 2021) and
masses in the range 107 − 1010 𝑀 and stellar metallicities in the thesan (Kannan et al. 2022; Smith et al. 2022) hydrodymanical sim-
ulations. In both simulations the ISM metallicity is combined with an
assumed dust-to-metal ratio to determine absolute dust attenuation.
3 http://flathub.flatironinstitute.org/group/sam-forecasts To compute the wavelength-dependent attenuation, flares employ a

MNRAS 000, 1–10 (2022)

 UV continuum slopes at 𝑧 ' 8 − 16 9
Charlot & Fall (2000) dust model with an attenuation curve somewhat ACKNOWLEDGEMENTS
stepper than the Calzetti et al. (2000) law, whereas thesan implement
F. Cullen and T. M. Stanton acknowledge support from a UKRI
full dust radiative transfer using skirt (Camps & Baes 2020). Both
Frontier Research Guarantee Grant (PI Cullen; grant reference
simulations predict slightly larger 𝛽 values compared to the sc-sam
EP/X021025/1). R. J. McLure, D. J. McLeod, J. S. Dunlop, C. Don-
model and slightly shallower 𝛽 − 𝑀UV slopes (d𝛽/dMUV ' −0.1).
nan, R. Begley and M. L. Hamadouche, acknowledge the support of
However, neither are obviously incompatible with our data and paint
the Science and Technology Facilities Council. A. C. Carnall thanks
the same basic picture of a uniformly blue (𝛽 . −2) galaxy popula-
the Leverhulme Trust for their support via a Leverhulme Early Career
tion at 𝑧 > 8. Overall, the reasonable agreement between the different
Fellowship. R. A. A. Bowler acknowledges support from an STFC
simulations, and their consistency with early JWST observations, is
Ernest Rutherford Fellowship (grant number ST/T003596/1).
encouraging.
We would like to thank Aaron Yung for kindly helping with the sc-
Based on this comparison, we can infer that our early 𝑧 > 8 data are
sam lightcone data. We would also like to thank Dan Magee for useful
qualitatively consistent with the young, low-metallicity, moderately
discussions and support in the data reduction process. This work is
dust-reddened stellar populations predicted by theoretical models.
based on observations collected at the European Southern Observa-
tory under ESO programme ID 179.A-2005 and 198.A-2003 and on
data products produced by CALET and the Cambridge Astronomy
Survey Unit on behalf of the UltraVISTA consortium. For the pur-
pose of open access, the author has applied a Creative Commons
Attribution (CC BY) licence to any Author Accepted Manuscript
5 CONCLUSIONS version arising from this submission.
We have measured the rest-frame ultraviolet (UV) continuum slopes
(𝛽) of 61 galaxies in the redshift range 8 < 𝑧 < 16, using a com- DATA AVAILABILITY
bination of James Webb Space Telescope (JWST) ERO and ERS
NIRcam imaging and ground-based near-infrared imaging of the All JWST and HST data products are available via the Mikulski
COSMOS/UltraVISTA field. The primary aim of this analysis is to Archive for Space Telescopes (https://mast.stsci.edu). Ultra-
determine whether there is any evidence for an evolution in the typical VISTA DR5 will shortly be made available through ESO. Additional
UV colours of the new population of 𝑧 > 8 galaxies being uncovered data products are available from the authors upon reasonable request.
by JWST. We present early estimates of the average values of 𝛽 at
these redshifts and, using the large dynamic range in UV luminosity
enabled by the combination of JWST and ground-based imaging, we REFERENCES
investigate evidence for a 𝛽 − 𝑀UV relation in our sample. Our main Adams N. J., et al., 2023, MNRAS, 518, 4755
conclusions can be summarized as follows. Aihara H., et al., 2019, PASJ, 71, 114
Atek H., Shuntov M., Furtak L. J., Richard J., Kneib J.-P., Mahler Adi Zitrin
(i) Using a power-law fitting technique, we find that our full sam- G., McCracken Clotilde Laigle Stéphane Charlot H. J., 2022, arXiv e-
ple displays a weighted mean value of h𝛽i = −2.10 ± 0.05, with a prints, p. arXiv:2207.12338
corresponding median value of 𝛽 = −2.29 ± 0.09. This result implies Bertin E., Arnouts S., 1996, A&AS, 117, 393
Bhatawdekar R., Conselice C. J., 2021, ApJ, 909, 144
that even the faintest galaxies that JWST has so far uncovered at
Bouwens R. J., et al., 2010, ApJ, 708, L69
𝑧 ' 8 − 16 have, on average, UV colours no more extreme than the Bouwens R. J., et al., 2014, ApJ, 793, 115
bluest galaxies in the local Universe (e.g., NGC 1705; 𝛽 = −2.46). Bowler R. A. A., Jarvis M. J., Dunlop J. S., McLure R. J., McLeod D. J.,
(ii) We find evidence for a 𝛽 − 𝑀UV relation in our sample, Adams N. J., Milvang-Jensen B., McCracken H. J., 2020, MNRAS, 493,
such that brighter UV galaxies display somewhat redder UV slopes 2059
(d𝛽/dMUV = −0.17 ± 0.06). This value is consistent with the slope Brammer G. B., van Dokkum P. G., Coppi P., 2008, ApJ, 686, 1503
derived at 𝑧 ' 5 by Rogers et al. (2014) and Bouwens et al. (2014) Bruzual G., Charlot S., 2003, MNRAS, 344, 1000
(d𝛽/dMUV = −0.12 ± 0.02 and d𝛽/dMUV = −0.14 ± 0.02, respec- Calzetti D., Kinney A. L., Storchi-Bergmann T., 1994, ApJ, 429, 582
tively). Compared to the the 𝑧 ' 5 relations, we find that the Calzetti D., Armus L., Bohlin R. C., Kinney A. L., Koornneef J., Storchi-
galaxies in our sample are bluer than their 𝑧 ' 5 counterparts. Bergmann T., 2000, ApJ, 533, 682
Camps P., Baes M., 2020, Astronomy and Computing, 31, 100381
The inverse-variance weighted mean offset for the full sample is
Castellano M., et al., 2022, arXiv e-prints, p. arXiv:2207.09436
hΔ𝛽i = −0.38 ± 0.09.
Charlot S., Fall S. M., 2000, ApJ, 539, 718
(iii) Examining the 𝛽 estimates for individual galaxies in our sam- Chisholm J., et al., 2022, MNRAS, 517, 5104
ple, we find no strong evidence for galaxies with ultra-blue UV slopes Cullen F., McLure R. J., Khochfar S., Dunlop J. S., Dalla Vecchia C., 2017,
that would indicate extreme stellar populations (i.e., 𝛽 ≤ −3). The MNRAS, 470, 3006
majority of galaxies in our sample with 𝛽 ≤ −3 have significant Cullen F., et al., 2019, MNRAS, 487, 2038
uncertainties, and appear to be consistent with the well-known blue Curtis-Lake E., et al., 2022, arXiv e-prints, p. arXiv:2212.04568
bias in the 𝛽 scatter for sources near to the source detection threshold. Donnan C. T., et al., 2022, MNRAS,
Dunlop J. S., McLure R. J., Robertson B. E., Ellis R. S., Stark D. P., Cirasuolo
M., de Ravel L., 2012, MNRAS, 420, 901
Overall, we find that the new population of 8 < 𝑧 < 16 galaxies
Dunlop J. S., et al., 2013, MNRAS, 432, 3520
being uncovered with JWST is on average bluer than the population
Euclid Collaboration et al., 2022, A&A, 658, A126
at lower redshfits, but does not show evidence of a substantial fraction Finkelstein S. L., et al., 2012, ApJ, 756, 164
of ultra-blue 𝛽 < −3 sources. Indeed, our sample is consistent with Finkelstein S. L., et al., 2022, arXiv e-prints, p. arXiv:2207.12474
the young (' 100 Myr), metal-poor (' 0.1𝑍 ) and moderately dust- Harikane Y., et al., 2022, arXiv e-prints, p. arXiv:2208.01612
reddened galaxies predicted by current theoretical galaxy formation Hudelot P., et al., 2012, VizieR Online Data Catalog, p. II/317
models at these redshifts. Inoue A. K., Shimizu I., Iwata I., Tanaka M., 2014, MNRAS, 442, 1805

MNRAS 000, 1–10 (2022)

10 F. Cullen et al.
Kannan R., Smith A., Garaldi E., Shen X., Vogelsberger M., Pakmor R.,
Springel V., Hernquist L., 2022, MNRAS, 514, 3857
Lovell C. C., Vijayan A. P., Thomas P. A., Wilkins S. M., Barnes D. J.,
Irodotou D., Roper W., 2021, MNRAS, 500, 2127
McCracken H. J., et al., 2012, A&A, 544, A156
McLure R. J., et al., 2011, MNRAS, 418, 2074
McLure R. J., et al., 2018, MNRAS, 476, 3991
Naidu R. P., et al., 2022, arXiv e-prints, p. arXiv:2207.09434
Oesch P. A., et al., 2014, ApJ, 786, 108
Oke J. B., 1974, ApJS, 27, 21
Oke J. B., Gunn J. E., 1983, ApJ, 266, 713
Pontoppidan K., et al., 2022, arXiv e-prints, p. arXiv:2207.13067
Roberts-Borsani G., et al., 2022, arXiv e-prints, p. arXiv:2210.15639
Robertson B. E., Ellis R. S., Dunlop J. S., McLure R. J., Stark D. P., 2010,
Nature, 468, 49
Rogers A. B., McLure R. J., Dunlop J. S., 2013, MNRAS, 429, 2456
Rogers A. B., et al., 2014, MNRAS, 440, 3714
Schaerer D., 2002, A&A, 382, 28
Smith A., Kannan R., Garaldi E., Vogelsberger M., Pakmor R., Springel V.,
Hernquist L., 2022, MNRAS, 512, 3243
Somerville R. S., Gilmore R. C., Primack J. R., Domínguez A., 2012, MN-
RAS, 423, 1992
Somerville R. S., Popping G., Trager S. C., 2015, MNRAS, 453, 4337
Somerville R. S., et al., 2021, MNRAS, 502, 4858
Speagle J. S., 2020, MNRAS, 493, 3132
Tacchella S., et al., 2022, ApJ, 927, 170
Topping M. W., Stark D. P., Endsley R., Plat A., Whitler L., Chen Z., Charlot
S., 2022, arXiv e-prints, p. arXiv:2208.01610
Treu T., et al., 2022, arXiv e-prints, p. arXiv:2206.07978
Vázquez G. A., Leitherer C., Heckman T. M., Lennon D. J., de Mello D. F.,
Meurer G. R., Martin C. L., 2004, ApJ, 600, 162
Vijayan A. P., Lovell C. C., Wilkins S. M., Thomas P. A., Barnes D. J.,
Irodotou D., Kuusisto J., Roper W. J., 2021, MNRAS, 501, 3289
Wilkins S. M., Bouwens R. J., Oesch P. A., Labbé I., Sargent M., Caruana J.,
Wardlow J., Clay S., 2016, MNRAS, 455, 659
Yung L. Y. A., Somerville R. S., Finkelstein S. L., Popping G., Davé R.,
2019a, MNRAS, 483, 2983
Yung L. Y. A., Somerville R. S., Popping G., Finkelstein S. L., Ferguson
H. C., Davé R., 2019b, MNRAS, 490, 2855
Yung L. Y. A., et al., 2022a, arXiv e-prints, p. arXiv:2210.04902
Yung L. Y. A., et al., 2022b, MNRAS, 515, 5416
Zheng W., et al., 2012, Nature, 489, 406

This paper has been typeset from a TEX/LATEX file prepared by the author.

MNRAS 000, 1–10 (2022)