The ultraviolet continuum slopes of high-redshift galaxies: evidence for the emergence of dust-free stellar populations at z > 10

We present an analysis of the ultraviolet (UV) continuum slopes ( 𝛽 ) for a sample of 172 galaxy candidates at 8 < 𝑧 phot < 16 selected from a combination of JWST NIRCam imaging and COSMOS/UltraVISTA ground-based near-infrared imaging. Focusing primarily on a new sample of 121 galaxies at ⟨ 𝑧 ⟩ ≃ 11 selected from ≃ 320 arcmin 2 of public JWST imaging data across 15 independent data sets, we investigate the evolution of 𝛽 in the galaxy population at 𝑧 ≥ 9. We find a significant trend between 𝛽 and redshift, with the inverse-variance weighted mean UV slope evolving from ⟨ 𝛽 ⟩ = − 2 . 17 ± 0 . 06 at 𝑧 = 9 . 5 to ⟨ 𝛽 ⟩ = − 2 . 59 ± 0 . 06 at 𝑧 = 11 . 5. Based on a comparison with stellar population models including nebular continuum emission, we find that at 𝑧 > 10 . 5 the average UV continuum slope is consistent with the intrinsic blue limit of dust-free stellar populations ( 𝛽 int ≃ − 2 . 6 ) . These results suggest that the moderately dust-reddened galaxy population at 𝑧 < 10 was essentially unattenuated at 𝑧 ≃ 11. The extremely blue galaxies being uncovered at 𝑧 > 10 place important constraints on dust attenuation in galaxies in the early Universe, and imply that the already observed galaxy population is likely supplying an ionising photon budget capable of maintaining ionised IGM fractions of ≳ 5 per cent at 𝑧 ≃ 11.


INTRODUCTION
The capability of JWST to provide very deep, near-/mid-infrared imaging and spectroscopy at  > 2m is revolutionising our understanding of the earliest galaxies.Before the launch of JWST, a small number of  > 10 galaxy candidates had been discovered using the Hubble Space Telescope (HST) (e.g.Ellis et al. 2013;Oesch et al. 2016).Now, deep NIRCam imaging surveys are revealing hundreds of these systems (e.g.Finkelstein et al. 2022;Adams et al. 2023a,b;Austin et al. 2023;Bouwens et al. 2023;Casey et al. 2023;Donnan et al. 2023a,b;Finkelstein et al. 2023b;Harikane et al. 2024;Hainline et al. 2024;Leung et al. 2023;McLeod et al. 2024;Robertson et al. 2023;Donnan et al. 2024).Moreover, increasing numbers of these galaxies have been confirmed spectroscopically (e.g.Arrabal Haro et al. 2023a,b;Bunker et al. 2023;Curtis-Lake et al. 2023;Fujimoto et al. 2023a;Harikane et al. 2023;Hsiao et al. 2023;Roberts-Borsani et al. 2023;Tang et al. 2023).The discovery of a large abundance of galaxies at  > 10 has been one of the key early JWST results.
The fact that we are finding large numbers of  > 10 galaxies, including a number of relatively ultraviolet (UV) bright objects ( UV ≤ −20; e.g.Naidu et al. 2022;Castellano et al. 2023;Casey et al. 2023), represents a challenge for models of early galaxy for-★ E-mail: fergus.cullen@ed.ac.uk mation.At  ≲ 8 the evolution of the UV luminosity function (LF) can be readily explained assuming no redshift evolution in the star formation efficiency (e.g.Bouwens et al. 2015;Mason et al. 2015;Mashian et al. 2016;Tacchella et al. 2018;Harikane et al. 2022).If this model is extrapolated to  > 8, a continued evolution of the UV LF is predicted, but this is in tension with the JWST observations (e.g.Finkelstein et al. 2023a;Harikane et al. 2024).In fact, current constraints suggest that the UV LF does not evolve at all between  = 9 and  = 11 (at least at  UV ≲ −20; McLeod et al. 2024).
A number of explanations for the non-evolving UV LF have been proposed including (i) an evolving star-formation efficiency, for example via a reduction in the stellar feedback efficiency at  > 10, including the potential for 'feedback-free' star formation (Dekel et al. 2023;Yung et al. 2024); (ii) a bias towards young, highly star-forming galaxies up-scattered with respect to the median UV magnitude versus halo mass relation (Mason et al. 2023;Shen et al. 2023); (iii) inherent uncertainties related to the spectral energy distributions (SEDs) of low-metallicity stellar populations (Inayoshi et al. 2022); and (iv) a transition to essentially dust-free star formation at  > 10 (Ferrara et al. 2023;Ferrara 2023;Ziparo et al. 2023).In this paper, we focus on the last of these scenarios, presenting an analysis of the UV continuum SEDs of galaxy candidates in the redshift range 8 <  phot < 16 to investigate whether a transition to zero/negligible dust attenuation at  > 10 is consistent with the current data.
It is well established that the stellar UV continuum, which can be described by a power-law (   ∝   for  rest ≃ 0.12 − 0.3 m), is an excellent probe of dust obscuration in star-forming galaxies (e.g.Calzetti et al. 1994;Meurer et al. 1999).In general, the UV continuum slope, , is sensitive to the light-weighted age, metallicity and dust attenuation of the population of massive stars in a galaxy (i.e.O-and B-type stars, with ages ≲ 100 Myr).However, at the highest redshifts -and more broadly for all young star-forming galaxiesage and metallicity effects become subdominant, and  is especially sensitive to dust (e.g.Tacchella et al. 2022).Indeed, this connection between  and dust attenuation has been demonstrated directly in studies sensitive to infrared dust emission up to  ≃ 8 (see e.g.Bowler et al. 2024).By providing deep infrared imaging out to  = 5 m, JWST/NIRCam now enables robust estimates of  for galaxies at  > 10 (e.g.Cullen et al. 2023;Topping et al. 2024), probing dust in galaxies at the earliest cosmic epochs.
In an earlier work, we presented an examination of the UV continuum slopes of galaxies at  > 8 selected from a combination of early JWST/NIRCam imaging and ground-based near-infrared imaging from COSMOS/UltraVISTA (Cullen et al. 2023).This initial sample comprised 61 galaxies at a mean redshift of ⟨⟩ = 10 spanning a factor of ≃ 80 in UV luminosity (−22.6 <  UV < −17.9).We found  slopes that were, on average, bluer than their lower-redshift counterparts at fixed  UV .These results had been tentatively anticipated with HST (e.g.Dunlop et al. 2013;Wilkins et al. 2016;Bhatawdekar & Conselice 2021;Tacchella et al. 2022) and were in agreement with a number of theoretical model predictions (e.g.Yung et al. 2019;Vijayan et al. 2021;Kannan et al. 2022).Independent JWST analyses painted a similar picture (e.g.Topping et al. 2022;Nanayakkara et al. 2023).These early studies of  ≃ 10 objects suggested stellar populations consistent with the young, low-metallicity, and moderately dust-reddened stellar populations anticipated by theoretical models.
Here, we expand upon our initial analysis using a new sample of galaxy candidates at  > 9 selected from a number of public JWST imaging surveys.Our new set of candidates builds on the sample presented in McLeod et al. (2024), which was used to robustly estimate the  ≃ 11 luminosity function.These new  > 9 galaxy candidates are selected across ≃ 320 arcmin 2 of deep JWST/NIRCam imaging, which represents a ≃ 7× increase in area compared to Cullen et al. (2023).As a result, we can now significantly improve upon our previous analysis in terms of the total number of  > 9 candidates.Crucially, our new sample now contains a significant number of galaxy candidates at  > 10, with an average redshift of ⟨⟩ = 10.7, and enables us to place robust constraints on the evolution of  up to  ≃ 12.
The new analysis presented here complements the recent study of Topping et al. (2024), who present an investigation of the UV continuum slopes of galaxies up to  ≃ 12 in the JADES survey (Eisenstein et al. 2023), finding that by  = 12 the average value of  is extremely blue ( ≃ −2.5).In this paper, we explore evidence for a similar trend in an independent sample of galaxies selected across a wide range of current JWST imaging data sets.Our sample is selected across an area ≃ 2× larger than the full JADES area (320 versus 175 arcmin 2 ) and, crucially, across 11 independent, noncontiguous fields, thereby mitigating the effect of cosmic variance.Although the JADES imaging used in Topping et al. (2024) is deeper, the increase in area means that our sample is on average brighter, and the resulting samples sizes at  > 9 are comparable.
The paper is structured as follows.In Section 2 we describe the data and sample selection, and provide details of the sample properties and our method for determining .In Section 3 we present the results of our analysis, focusing on the evolution of  with  and absolute UV magnitude ( UV ).Our primary new finding is that, by  ≃ 11, the typical UV continuum slope of the galaxy population is extremely blue ( ≃ −2.6).In Section 4 we discuss the implications of this result and demonstrate, based on an analysis of theoretical stellar population models, that the UV continuum slopes of galaxies at  ≳ 10 are consistent with dust-free star formation at this epoch.Finally, we summarise our main conclusions in Section 5. Throughout we use the AB magnitude system (Oke 1974;Oke & Gunn 1983), and assume a standard cosmological model with  0 = 70 km s −1 Mpc −1 , Ω  = 0.3 and Ω Λ = 0.7.

SAMPLE SELECTION AND PROPERTIES
Our primary high-redshift galaxy sample is drawn from the widearea JWST search for  > 8.5 galaxies presented in McLeod et al. (2024).We additionally incorporated three other ultra deep datasets not included in the original McLeod et al. (2024) sample: NGDEEP (Bagley et al. 2023), the first data release of the JWST Advanced Deep Extragalactic Survey (JADES, Eisenstein et al. 2023), and the additional ≃ 25 arcmin 2 region taken in parallel to the recent UNCOVER NIRSpec observations (hereafter "UNCOVER-South").This primary sample was augmented by additional galaxies at  > 7.5 drawn from COSMOS/UltraVISTA and JWST that were presented in our earlier work (Cullen et al. 2023).In the following, we discuss the sample selection and properties of the various data sets used in this work.

The wide-area JWST sample
Our wide-area JWST sample was drawn from 15 independent publicly-available JWST imaging datasets within 11 extragalactic fields, covering an on-sky area of ≃ 320 arcmin 2 .A list of these fields, including the proposal ID and the PI name of each dataset is given in Appendix A (table A1).We also include references to the relevant survey paper, ancillary data, and lensing maps where applicable.With the exception of NGDEEP, JADES, and UNCOVER-South, all of the datasets presented here are the same as those described in McLeod et al. (2024); we refer interested readers to Section 2 of McLeod et al. (2024) for a detailed description of the datasets in common.
Our inclusion of three new fields yields a further 50 sources compared to the McLeod et al. (2024) sample.We have identified 9 galaxy candidates at  ≥ 9 in our latest reduction of the NGDEEP observations (Bagley et al. 2023) (see below for candidate selection details).The NGDEEP dataset covers the HST UDF parallel 2 field (Oesch et al. 2007;Bouwens et al. 2011;Ellis et al. 2013), with another NIRCam module covering an adjacent area contained within the GOODS-South CANDELS footprint; therefore, we supplemented the NGDEEP NIRCam imaging with the ACS F435W, F606W, F775W, F814W and F850LP imaging released as part of the Hubble Legacy Fields (Illingworth et al. 2013;Whitaker et al. 2019).We also included data from the first data release of the JWST Advanced Deep Extragalactic Survey (JADES, Eisenstein et al. 2023) 1 .The first data release (Rieke et al. 2023) covers the 'deep' portion of the imaging across an area of ≃ 25 arcmin 2 , from which we have selected a further 27  ≥ 9 galaxy candidates.Finally, we include 14  > 9 candidates that were found in the ≃ 25 arcmin 2 of UNCOVER-South.The coordinates,  phot and  UV for these new candidates are given in Table A3.
All of the JWST/NIRCam data in each field was reduced using the Primer Enhanced NIRCam Image Processing Library (pencil) software, a customised version of the JWST pipeline version 1.6.2(Magee et al. in prep).The CRDS pmap varied slightly between each dataset due to differences in the date on which the data were taken, reduced and analysed, although each pmap was sufficiently up-to-date (≥ 0989) to take into account the most recent zero-point calibrations.Again, we refer interested readers to McLeod et al. (2024) for a detailed description of the various pencil pipeline steps and the full data reduction routine.In Table A2 we report the global 5 limiting magnitudes of the HST and JWST imaging in each field, calculated using 0.35 ′′ −diameter apertures and corrected to total assuming a point-source correction (see McLeod et al. 2024).

Construction of galaxy catalogues
Prior to galaxy candidate selection, we first homogenised the point spread function (PSF) of each of our images to match the PSF full width half maximum (FWHM) of the F444W imaging.This enabled us to measure consistent photometry across all of the photometric filters, removing any systematics that arise as a result of differences in the PSF curve of growth.
For each data set, we constructed multiwavelength photometry catalogues running Source Extractor (Bertin & Arnouts 1996) in dual image mode and performing aperture photometry in 0.35 ′′ −diameter apertures.Catalogues were created primarily with the F200W imaging as the detection band, although we also created additional F277W-, F356W-and F444W-detected catalogues in order to include any sources that had been missed due to a low signal-to-noise ratio (SNR) in the F200W imaging.
To calculate photometric uncertainties, we adopted the method described in McLeod et al. (2024).We first generated a grid of nonoverlapping apertures of diameter 0.35 ′′ spanning the full field of view and identified apertures corresponding to blank sky using a Source Extractor segmentation map.For a given object, we measured the median absolute deviation (MAD) of the nearest 150 − 200 blank sky 0.35 ′′ − diameter apertures on an object-by-object basis and scaled to the standard deviation using  ≃ 1.4826 MAD.At this initial stage, we retained all sources detected at ≥ 5 in our detection catalogues.Finally, to reduce the possibility of artefact contamination, we required a ≥ 3 detection in any one of the other detection bands.
We performed spectral energy distribution (SED) fitting on all sources in our initial catalogues using the SED fitting code LePhare (Arnouts & Ilbert 2011), adopting the Bruzual & Charlot (2003) stellar population synthesis (SPS) library, a Chabrier (2003) IMF and the Calzetti et al. (2000) dust attenuation law.We modelled all sources assuming a declining star formation history with  ranging from 0.1 to 15 Gyr and stellar metallicities of 0.2 Z ⊙ and Z ⊙ .We allowed  V to vary between 0 and 6 and accounted for intergalactic medium (IGM) absorption using the prescription of Madau (1995).
We measured  UV for each candidate from the best-fitting SED using a top-hat filter centred on  rest = 1500 Å with Δ = 100 Å.To correct our  UV to total magnitudes, we scaled to the FLUX_AUTO measurements and then applied an additional correction of ≃ 10 per cent to account for light outside the Kron aperture (McLeod et al. 2024).For the cluster fields, we corrected  UV for gravitational lensing by determining the magnification factor () using publicly available lensing maps (see Table A1).For fields where multiple lensing maps are available, we adopted the median value of  across all maps.For UNCOVER-South, where there are no public lensing maps, we adopted a flat  = 1.2 ± 0.2 magnification, motivated by the typical  values around the southern boundary of the Furtak et al. (2023) UNCOVER lensing map.While  is unaffected by changes in , we include an additional systematic uncertainty in  UV to account for this assumed  floor.

High-redshift galaxy candidate selection
To select our final sample, we first retained galaxies with a photometric redshift  phot > 8.5 , a goodness of fit of  2  ≤ 10, and at least Δ 2 ≥ 4 between the goodness of fit of the primary redshift solution and secondary lower-redshift solution.To ensure a more robust final sample, we then required a detection of ≥ 8 in at least one of the F150W, F200W, or F277W filters.We also required a non-detection (≤ 2) in any of the F090W, F115W and HST ACS imaging filters that were available for a given data set.In practice, the ≤ 2 in F115W criterion restricts our sample to galaxies with  phot ≥ 9.Where there was a lack of F115W imaging, we added a further criterion, whereby we required either a (≤ 2) non-detection in F150W, or at least a factor of two increase in flux between F150W and F200W.This last criterion confined us to redshifts  ≳ 11.5 for those datasets, but mitigated any potential contamination arising from F090W-F150W dropouts that have highly uncertain photometric redshifts.For each candidate, additional checks were performed using the photometric redshift code EAZY (Brammer et al. 2008) as described in McLeod et al. (2024).
Across the 11 fields, we selected a total of 121 galaxy candidates at  phot ≥ 9. Throughout this paper, we refer to these galaxies as our primary wide-area JWST sample.We use this sample to study the relationship between  and redshift over the redshift range 9 ≤  ≤ 12 in Section 3.1.The distribution of  UV and  phot is shown by the blue data points in Fig. 1.The median values are  UV = −19.2and  phot = 10.6.

Additional galaxy candidates: the combined sample
The primary wide-area JWST sample was augmented using the sample analysed in our previous study of the UV continuum slope (Cullen et al. 2023).This sample was initially selected and presented in the UV luminosity function study of Donnan et al. (2023a) and combines a ground-based sample of bright galaxies at  > 7.5 selected from the COSMOS/UltraVISTA field, and a JWST-selected sample at  > 8 from early NIRCam imaging in SMACS J0723, GLASS, and CEERS.
The selection criteria for this sample were less stringent (i.e.requiring only a 5 detection in bands red-ward of the Lyman break; for full details see Donnan et al. 2023a) and therefore contained a number of candidates not selected in the new wide-area JWST sample.We retained the 35 galaxies from the Cullen et al. (2023) JWST sample that were not in our wide-area selection.We included all 16 galaxies from the ground-based COSMOS/UltraVISTA sample; these galaxies, although at a lower median redshift ( = 8.3), provide the valuable dynamic range in  UV needed for studying the relationship between UV continuum slope and galaxy luminosity (which we explore in Section 3.2).
In total, we included 51 galaxy candidates from the sample presented in Donnan et al. (2023a) and Cullen et al. (2023).The distributions of  UV and  phot for the COSMOS/UltraVISTA and JWST galaxies are shown by the red and yellow data points in Fig. 1, respectively.The median values for the COSMOS/UltraVISTA sample are  UV = −21.8and  phot = 8.3; for the JWST sample, the median values are  UV = −19.0 and  phot = 9.9.We note that, because of the less stringent S/N constraints, the  and  UV values of these galaxies have larger uncertainties compared to our primary wide-area sample.
Combined with the primary wide-area sample, the total sample contains 172 candidates at  > 7.5 (Fig. 1); we refer to it as the 'combined' sample throughout this paper.This sample is used primarily in our analysis of the  −  UV relation in Section 3.2.The median values of the absolute UV magnitude and redshift for the combined sample are  UV = −19.4and  phot = 10.3.

Measuring the UV continuum slope
The primary goal of this paper is to study the UV continuum slopes of our high-redshift galaxy candidates.The slope of the UV continuum in the rest frame is characterised by a power-law index, , where   ∝   .To determine , we modelled the galaxy photometry covering rest-frame wavelengths  rest ≤ 3000 Å as a pure power law that includes both the effects of intergalactic medium (IGM) absorption and the Ly  damping wing.At  rest > 1216 Å, we modelled the effect of the Ly  damping wing using the prescription of Miralda-Escudé (1998) (equation 2).At wavelengths below Ly  ( rest ≤ 1216 Å), we assumed complete IGM attenuation (i.e. rest ≤1216 = 0).Our model consisted of three free parameters: (i) , the UV continuum slope; (ii) , the redshift of the galaxy; and (iii)  HI , the neutral hydrogen fraction of the surrounding IGM.To sample the full posterior distribution, we used the nested sampling code dynesty (Speagle 2020) assuming the following flat priors: −10 ≤  ≤ 10; 5 ≤  ≤ 20 and 0 ≤  HI ≤ 1.0.Examples of fits to nine galaxies in our sample are shown in Fig. 2 (blue dashed lines).
We found that for a small minority of objects (≃ 10 per cent of the full sample), our fits returned unrealistically small uncertainties (  < 0.1) that significantly biased the population average estimates.Based on our -recovery simulations (Section 3.1.2),we found that for the typical global imaging depths of our dataset,  could be recovered for the brightest sources ( UV ≲ −20) to within Δ = ±0.25.We therefore adopted this value as a minimum floor on the  uncertainty.Throughout this paper, unless otherwise stated, we adopt the best-fitting redshifts and corresponding uncertainties derived from these fits in our analysis.We use these redshifts to calculate the bestfit values of  UV and their corresponding uncertainties.The derived values of ,  phot , and  UV for our primary wide-area JWST sample are given in Table 1.
Our method differs slightly from other similar studies in the literature which predominantly fix the redshift to  phot and fit only to the photometric data redward of the Lyman break (e.g.Topping et al. 2024;Morales et al. 2024).However, we have verified that both methods ultimately yield comparable results, with offsets of Δ < 0.1.For example, if we adopt the methodology described in Topping et al. (2024) for our sample, we find a median offset of Δ < +0.04.Although we advocate for our adopted method, which has the benefit of incorporating all photometric data and fully accounting for redshift uncertainties, our results remain essentially unchanged if we follow other methodologies common in the literature.

The effect of Ly 𝛼 damping wings and proximate DLAs
Evidence for Ly  damping wings and even strong proximate damped Ly  systems (DLAs) has been revealed from early JWST NIR-Spec/PRISM spectroscopy of galaxies at  > 8 (Curtis-Lake et al. 2023; Heintz et al. 2023;Umeda et al. 2023).Both effects act to soften the Ly  break and are a potential source of systematic bias in our derived value of .
To account for the effect of the Ly  damping wing, we included  HI as a free parameter in our fits as described above.Although it is not possible to constrain  HI from photometric data alone, our method allowed us to marginalise over the uncertainties related to the unknown Ly  damping wing strength.We found that marginalising over  HI had a minor effect on our results compared to assuming a fixed value of  HI = 1.0; the median value of  became marginally redder (Δ = 0.02) and the median uncertainties were unchanged.Object-to-object variations of up to Δ = ±0.3 were observed, but these were completely consistent with scatter due to the typical photometric uncertainties.In general, we find that the average value of  derived from broadband photometry is not strongly affected by the unknown strength of the Ly  damping wing.
Recently, Heintz et al. (2023) presented evidence for excess DLA absorption in three galaxies at  > 8. Heintz et al. (2023) argue that the strong absorption is likely due to large H i gas reservoirs in the interstellar and circumgalactic medium along the line-of-sight; this leads to an attenuation of flux redward of Ly  in excess of the neutral IGM damping wing (see, for example, Figure 1 of Heintz et al. 2023 and the discussion in Keating et al. 2023).We investigated the effect of DLA absorption by including a DLA at the source redshift with a neutral hydrogen column density of  HI = 10 22.5 cm −2 in our model and re-fitting the sample.This column density is the upper limit of the Heintz et al. (2023) measurements, and also at the upper end of the  HI distribution measured from the afterglow spectra of gamma-ray bursts (GRBs) (Tanvir et al. 2019).Again, we find that the overall effect on  is negligible (median Δ = −0.01 with respect A comparison of UV continuum slope fits to photometry ( phot ) and spectroscopy ( spec ) for the nine galaxies in our sample with JWST/NIRCam and NIRSpec/PRISM observations.In the first nine panels, the black curve is the NIRSpec/PRISM spectrum, with the grey shading representing the ±1 error spectrum.The square data points show the galaxy photometry, with all photometry at  ≥ 0.9 m coming from JWST/NIRCam.The solid red and dashed blue curves are the best-fitting models for the spectrum and photometry, respectively.In the majority of cases, the photometric and spectroscopic fits are fully consistent; however, for some galaxies, small differences in the best-fitting redshift and  are evident.The tenth panel (lower left) shows the comparison between  spec and  phot .The grey circular points show the nine individual galaxies and the red square point is the inverse variance-weighted mean.In general, we find very good agreement between the two estimates with ⟨ spec ⟩ = −2.33 ± 0.02 and ⟨ phot ⟩ = −2.39 ± 0.07.
to the default assumptions).The offset becomes more significant at the lowest redshifts where the difference between a DLA and the Ly  damping wing is more pronounced.For sources at  < 9.5, the median offset is Δ = −0.05,however, this is still not large enough to affect the results presented here.

Spectroscopic Sample
For nine galaxies in our primary sample, fully reduced JWST NIR-Spec/PRISM observations were available through the DAWN JWST Archive (DJA)2 .We used these spectroscopic data to assess the accuracy of our photometric  estimates.
The DJA reductions were performed using the software packages grizli (Brammer 2023) and msaexp (Brammer 2022).The basic details of the data reduction are presented in Valentino et al. (2023) and Heintz et al. (2023).The NIRSpec spectra were affected by wavelength-dependent slit losses due to objects not being fully encompassed within the small MSA slitlets (0.20 ′′ × 0.46 ′′ ) and/or slit misalignment.To calibrate the 1D spectra, we integrated them through the available NIRCam filters and scaled the values to match the observed photometry, employing a linear interpolation between each band.The nine flux-calibrated spectra are shown in Fig. 2.
The method used to fit  from spectroscopic data is essentially the same as the method described above (Section 2.3).The only difference is that, in the spectroscopic case, the model is smoothed to match the resolution of the NIRSpec/PRISM data rather than integrated through the relevant photometric filters.To smooth the models we first resampled to four pixels per full width half maximum (FWHM) element and then convolved with a Gaussian with fixed pixel width of  pix = 1.7.Examples of these model fits are shown in Fig. 2.

Comparing 𝛽 spec and 𝛽 phot
In Fig. 2 we show the fits to the photometry (blue) and to the fluxcalibrated spectroscopy (red) for the nine galaxies in our spectroscopic sample.In general, we find excellent agreement between the two estimates.The lower right panel in Fig. 2 shows the comparison of  spec and  phot , where it can be seen that, for N = 8/9 of the galaxies, the two values are consistent within 1.The inverse-variance weighted mean values of the two estimates are fully consistent with ⟨ spec ⟩ = −2.33 ± 0.02 and ⟨ phot ⟩ = −2.39 ± 0.07.
We find tentative evidence to suggest that, on average,  phot is systematically lower (bluer) than  spec .For the cases in which  phot <  spec , we find that the redshift estimated from the photometry is always larger than that estimated from the spectra, with ⟨ phot −  spec ⟩ ≃ 0.5.Our sample size is clearly small, but other authors (e.g.Arrabal Haro et al. 2023a) have observed a tendency for  phot to be systematically biased high at  > 8 which, if confirmed, might suggest a general bias towards bluer  phot .However, our results suggest that the average effect is likely to be small (Δ ≃ 0.05).Larger samples of deep NIRSpec/PRISM spectra should soon be available to help accurately assess potential biases in ⟨ phot ⟩ due to  phot systematics (e.g.CAPERS; GO 6368; PI: M. Dickinson).

RESULTS
In Fig. 3 we plot the  values for our combined sample as a function of redshift and  UV and highlight the diverse selection of datasets used in this study.The overall range of  values in our new wide-area JWST sample is consistent with that presented in Cullen et al. (2023), with measured values between  ≃ −5 and  ≃ −1, and we observe a number of similar features in the data.
First we again see a large scatter in , with individual estimates as blue as  ≃ −5.These extremely blue values ( < −3) are driven, at least in part, by the well-known blue bias in the estimates of  at faint magnitudes (e.g.Dunlop et al. 2012;Rogers et al. 2014;Cullen et al. 2023) and by the significant uncertainties for individual objects (the median uncertainty for the combined sample is ⟨  ⟩ = 0.5).Unsurprisingly, the scatter is noticeably reduced in our higher S/N wide-area sample, which contains fewer extremely blue outliers; for example, 22 per cent of the Cullen et al. ( 2023) sample have  < −3 compared to 10 per cent of the wide-area sample.Nevertheless, it is important to note that the existence of objects with  ≲ −4 − even in robust, high S/N, samples − emphasises the necessity of adopting empirical approaches to measuring  that will not artificially truncate the full observed distribution (as can happen when measuring  via fitting stellar population models; see e.g.Rogers et al. 2013).
In Table 2 we report the inverse-variance weighted mean  values for the primary sample and for the combined sample.For our primary sample, we find ⟨⟩ = −2.37 ± 0.03.As in Cullen et al. (2023), we prefer the weighted mean over a simple median to mitigate against the blue bias in the  scatter at faint magnitudes (this blue scatter is clearly visible in the bottom panel of Fig. 2).For the combined sample (i.e.including the Cullen et al. 2023 sources) we find ⟨⟩ = −2.32 ± 0.03.This marginally redder value is primarily due to the inclusion of bright objects at  < 10 from the ground-based COSMOS/UltraVISTA survey in the combined sample, as well as the slightly lower median redshift (see Fig. 1, and the results presented in Sections 3.1 and 3.2 below).
The ⟨⟩ values reported in Redshift (z) Figure 3. Plots of UV continuum slope  versus redshift (top) and versus absolute UV magnitude  UV (bottom) for the galaxies in our combined sample sample.In each panel, the points are colour-and symbol-coded by dataset (detailed in the legend in the top panel).The blue circular and brown square data points are taken from our previous sample presented in Cullen et al. (2023).All other data points come from our new wide-area JWST sample (see Section 2 for full sample details and Table A1 for an overview of the wide-area sample).
value of ⟨⟩ = −2.10 ± 0.05 we reported in Cullen et al. (2023).This is in part due to the fact that the new sample is marginally fainter, but the main reason is the higher proportion of candidates at  > 10.
In our new wide-area JWST sample ≃ 70 per cent (84/121) of the candidates have  phot > 10, compared to 43 per cent (26/61) in Cullen et al. (2023).The implication here is that galaxies at  > 10 are significantly bluer, and this can be seen in the top panel of Fig. 3.This clear evolution to bluer  is the main result of this study  Calzetti et al. 1994;Vázquez et al. 2004).As will be discussed in the following, we still find that up to  ≃ 10.5 the typical value of ⟨⟩ is consistent with local blue starburst galaxies.However, we now find that at  > 10.5 a significant fraction of galaxies have extremely blue UV continuum slopes with ⟨⟩ ≃ −2.6.While  = −2.6 represents the extreme blue-end of objects observed in the local Universe (e.g.Chisholm et al. 2022), we find that this is the typical value of  for galaxies at  ≃ 11 with absolute UV magnitudes in the range −20 ≲  UV ≲ −17.
In the following, we present a detailed analysis of the UV continuum slopes of our galaxy candidates and their dependence on redshift and UV luminosity ( UV ).We begin in Section 3.1 by investigating the evolution of ⟨⟩ with redshift and present the main result of this work.Then, in Section 3.2, we derive the  −  UV relation as a function of redshift for our sample, providing an update on the relation derived in Cullen et al. (2023).

The evolution of 𝛽 with redshift
An evolution to bluer UV colours at higher redshifts is expected as galaxies become progressively metal and dust-poor at earlier cosmic times.Before JWST, this trend had been observed up to  ≃ 8 (e.g.Finkelstein et al. 2012;Bouwens et al. 2012;Dunlop et al. 2013;Bouwens et al. 2014;Wilkins et al. 2016).Recently, Topping et al. (2024) have presented evidence for continued evolution beyond  = 8, finding extremely blue average UV slopes of ⟨⟩ ≃ −2.5 at  = 12 (we discuss the comparison between our results and Topping et al.

below).
In Fig. 4 we show the  −  relation for our primary sample (lefthand panel) and combined sample (right-hand panel).Focusing first on the primary sample, it can be seen that we observe a clear  −  trend in our data, with  becoming progressively bluer at higher redshifts.Fitting a linear relation to the individual data points (and accounting for the uncertainties in both  and ) yields where we have restricted the fit to  ≤ 12 due to the lack of candidates at higher redshifts.In Fig. 4 we also show the inverse-variance weighted mean  values in bins of redshift.These binned values (reported in Table 3) are fully consistent with the fit to the individual objects.We see that the typical  changes from ⟨⟩ = −2.59 at  ≃ 11.5 to ⟨⟩ = −2.17 at  ≃ 9.5, a relatively rapid evolution of Δ⟨⟩ ≃ 0.4 over ≃ 100 Myr of cosmic time.
We find that the redshift evolution observed in Fig. 4 is unlikely to ).In the highest redshift bin (11 <  < 12), the population mean value approaches this limit, indicating galaxies unattenuated by dust.In the right-hand panel, we show the same  −  relation for our combined sample (i.e., including the galaxy candidates from Cullen et al. 2023, as indicated in the legend).For clarity, we omit error bars in this panel.At  > 9, the JWST-selected galaxies of Cullen et al. (2023) are consistent with the same trend.Our data show evidence for a steep decline in  with redshift at  > 9, indicating that galaxies are growing progressively more dust and metal poor at early cosmic times.By  ≃ 11, the population-average  is consistent with negligible, even zero, dust attenuation.be strongly affected by differences in the typical  UV as a function of redshift.Significant differences in  UV could result in biases due to the known relation between  and  UV that has been demonstrated in numerous studies up to  ≃ 10 (e.g.Rogers et al. 2014;Bouwens et al. 2014;Cullen et al. 2023; the  −  UV relation for our new sample will be derived in Section 3.2).However, from Fig. 1 it can be seen that the  UV distribution of the primary sample does not evolve strongly between  = 9 and  = 12.From Table 3, the difference in ⟨ UV ⟩ between adjacent redshift bins is approximately |Δ⟨ UV ⟩| ≃ 0.5.If we assume d/d UV = −0.15(see Section 3.2), this difference translates to |Δ⟨⟩| ≃ 0.08, which cannot account for the observed changes in ⟨⟩.In Section 3.1.2and Appendix C, we provide further detailed discussions of possible selection effects and measurement biases present in our analysis.In a variety of tests, we conclude that the  −  trend is a real feature of our data.
Including the additional candidates from the combined sample, we find that they are fully consistent with the  −  relation derived from the primary wide-area JWST sample (right-hand panel of Fig. 4).At  > 9, the JWST-selected galaxies from Cullen et al. (2023) follow the same trend.In particular, it can be seen that these candidates also show extremely blue UV slopes at  ≃ 11.At  < 9, the addition of ground-based COSMOS/UltraVISTA candidates from the Cullen et al. (2023) sample provides an additional redshift bin at ⟨⟩ = 8.5.The value of ⟨⟩ = −1.94± 0.07 in this bin is consistent with an extrapolation of Equation 1.However, we note that the galaxies in the 7.5 <  < 9.0 redshift bin are significantly brighter by Δ UV ≃ −2, meaning that the average is likely biased with respect to the higher redshift bins (by up to Δ ≃ 0.25).Taking into account this bias would suggest that the evolution between  = 8 and  = 9 is likely to be shallower if compared to a sample at  = 8 with the same magnitude.The focus of this paper is on the evolution in  at  > 9, but future JWST studies will enable a robust determination of the redshift-evolution of  for samples matched in magnitude across a wider redshift range.

The transition to extremely dust-poor stellar populations
The key result of this paper, highlighted in Fig. 4, is that by  ≃ 11 the average UV continuum slope of the galaxy population with −21.5 <  UV < −17.5 is consistent with the dust-free limit expected from standard stellar population models and nebular physics.As we discuss in more detail in Section 4.1, this limit assumes a young, dust-free, stellar population surrounded by ionised gas such that the ionising continuum escape fraction is 0 per cent and the strength of the nebular continuum emission is maximised.Under these assumptions, the bluest UV continuum slopes expected are  ≃ −2.4 to −2.6 (e.g.Robertson et al. 2010;Cullen et al. 2017;Reddy et al. 2018).This theoretical limit appears to be validated by observations of dust-poor local star-forming galaxies (e.g.Chisholm et al. 2022).Any value bluer than this limit implies a non-zero escape fraction of ionising photons.Crucially, since dust acts to redden 3 , this limit also implies that it is extremely difficult to obtain UV slopes of  ≲ −2.6 in the presence of dust attenuation.As can be seen in Fig. 4, our results suggest that this limit is reached by  ≃ 11, implying that galaxies at this redshift are uniformly extremely dust-poor and essentially unattenuated.
The rapid transition to extremely blue  values at the highest redshifts in our sample is highlighted further in Fig. 5. Considering the steep nature of the  −  relation (Fig. 4), we fitted a toy model to our wide-area JWST sample with the following form: where  t is the redshift at which the population-average value of  transitions from  0 to  1 .Although this model is not physically motivated, it gives a useful indication of the approximate redshift above which galaxies are becoming uniformly extremely blue.
Fitting the model in Equation 2to our primary wide-area sample yields: The fit is shown in Fig. 5. Again, we can be sure that this difference in ⟨⟩ is not strongly biased due to differences in  UV by comparing the  UV distribution above and below  t .The inset panel of Fig. 5 (lower right) shows that the  UV distributions of the upper and lower redshift samples are fully consistent; both have ⟨ UV ⟩ = −19.2and a standard KS test returns a significance (−value) of  = 0.34, consistent with the null hypothesis that both samples are drawn from the same  UV distribution.
Interestingly, we find that the step function model (Fig. 5) provides a better fit to our data than the linear model (Fig. 4).However, based on the reduced chi-squared values ( 2  ), neither model provides a formally statistically acceptable fit.For the linear model, we find  2  = 3.16 compared to  2  = 2.01 for the step function model.These large  2  values are most likely due to a large intrinsic scatter in  at a fixed redshift.A given redshift bin will span a range of UV magnitude with a width of Δ UV ≃ 3 (Fig. 1), corresponding to Δ ≃ 0.45 assuming d/d UV = −0.15(see Section 3.2).Taking this into account, the fact that neither fit is formally acceptable is not particularly concerning.Although neither fit can be preferred, it is intriguing that the current data prefer a sharp transition above  ≃ 10.5.Larger samples, ideally with a larger fraction of robust spectroscopic redshifts, are needed to clarify the precise nature of the rate of evolution of ⟨⟩ at  > 9.
Taken together, the model fits in Figs. 4 and 5 suggest that at  = 11 the typical UV continuum slope of galaxies is ⟨ obs ⟩ ≲ −2.5.This implies that the galaxy population, across a relatively wide range of  UV , is consistent with the dust-free limit suggested by both theoretical models and observations of local galaxies.Moreover, our data also suggest that the transition to these extremely blue slopes occurs over a narrow redshift range, with galaxies at  ≃ 9.5 having on average blue, but not particularly extreme, UV continuum slopes of ⟨⟩ ≃ −2.2, consistent with a moderate amount of dust attenuation ( UV ≃ 0.5 − 1, assuming  int = −2.6;see McLure et al. 2018).Taken at face value, this implies a rapid buildup of dust in galaxies within the ≃ 100 Myr time frame between  = 11.5 and  = 9.5, consistent with models that predict an efficient ejection of dust in the earliest phases of galaxy formation (Ferrara et al. 2023;Ziparo et al. 2023), negligible amounts of dust being formed in the first phases of star formation (Jaacks et al. 2018) or efficient dust destruction at the highest redshifts (e.g.Esmerian & Gnedin 2023).

Selection effects and measurement bias
It is worth considering the possibility that the extremely blue UV slopes at  > 10 are a result of selection effects and/or measurement biases.Indeed, the measurement of  from broadband photometry is known to be affected by subtle biases.For example, the bias towards bluer values of  at faint UV magnitudes has been extensively documented (e.g.Bouwens et al. 2010;Dunlop et al. 2012;Rogers et al. 2014).In Cullen et al. (2023), we found that for our JWST sample, the average value of  was biased toward bluer values at  UV,obs ≥ −19.3.If we apply the relation derived in Cullen et al. (2023) to the median  UV of galaxies above and below the  = 10.55 transition in Fig. 5 ( UV,obs = −19.2) we find that these ⟨⟩ estimates could be biased blue by Δ = −0.02.A bias at this level would clearly not affect our results.We have also confirmed that the same steep redshift- trend remains if we restrict our sample to the brightest galaxies in our wide-area JWST sample with  UV ≤ −19.5 (i.e.those which should not be affected by a measurement bias).We therefore conclude that it is unlikely that our results are driven by a blue bias due to faint galaxies.We have also explored potential redshift-dependent biases, which might arise from the fact that different combinations of filters are used and different portions of the UV spectrum are sampled depending on the redshift of the candidate.For example, the UV continuum is sampled by the F150W, F200W, and F277W filters at  < 11, compared to the F200W, F277W, and F356W filters at  > 11.To test for a redshift bias we ran a simple simulation in which we first constructed 20, 000 simple power-law SEDs with an intrinsic slope of  int = −2.2spread uniformly across the redshift range 9 <  < 12 and the  UV range −20.5 <  UV < −17.0.We then applied IGM attenuation using the prescription of Inoue et al. (2014).Photometry was generated in each of the observed filters (see Table A2) and scattered according to the typical imaging depths (averaging the depths across multiple fields where appropriate).The 'observed' UV continuum slopes ( recovered ) were then recovered for the simulated galaxies that met our selection criteria.
The resulting bias as a function of redshift is shown in Fig. 6.It can be seen that we do not observe a strong redshift-dependent effect.The largest bias (Δ = −0.11) is seen in a narrow redshift interval around  ≃ 9.8, which corresponds to the redshift at which the Lyman break passes between the F115W and F150W filters.At this specific redshift, the photometric redshift is almost uniformly biased high, resulting in an underestimate of the true UV slope.In general, the bias across all redshifts is negligible (Δ ≃ −0.01) and crucially, at  > 10.5 there is no evidence for a strong blue bias.
It is also clear from Fig. 6 that the scatter in  is driven primarily by the photometric redshift uncertainties.When the photometric redshift is underestimated,  is overestimated and vice versa.We find that the typical scatter in the recovered  at all  and  UV is   ≃ 0.25, which motivates our decision to set this as a minimum error floor on individual  estimates (see Section 2).
Finally, we find that our results are unlikely to be driven by selection effects.This is partly highlighted by the similarity of the  UV distributions above and below the  = 10.55 transition in Fig. 5, which demonstrates that our sample selection does not favour intrin-  4. In the left-hand panel the black solid line is the best-fitting  −  UV relation (Equation 4) with the light-grey shaded region showing the 95 per cent confidence interval around our best-fitting relation.We find evidence for a  −  UV relation with best-fitting slope of d/d UV = −0.15± 0.03 (i.e.5 significance).In the centre and right-hand panels, the black solid line shows the best-fitting  −  UV relations where we have fixed the slope to be the same as the slope determined at  < 10 (see text for details).In the centre and right-hand panels, the black dashed line shows the best-fitting  −  UV relation at  < 10 to highlight the redshift evolution of the normalisation.The best-fitting  −  UV parameters are given in Table 5.
sically fainter, and hence bluer, objects in the higher redshift bins.Two further selection tests are described in detail in Appendix C. Across all tests, we find that the evolution of  with redshift we observe is not caused by a more efficient selection of bluer objects at higher redshifts.

The 𝛽 − M UV relation at z = 8 − 12
The  −  UV relation, or colour-magnitude relation, traces the variation of the dust and stellar population properties of galaxies as a function of their UV luminosity.Observations up to  ≃ 10 have shown that a correlation exists such that fainter galaxies are on average bluer (e.g.Meurer et al. 1999;Bouwens et al. 2009;Rogers et al. 2014;Bouwens et al. 2014;Cullen et al. 2023;Topping et al. 2024).These observations suggest that fainter galaxies are typically younger and less dust and metal enriched, in agreement with predictions of theoretical models (e.g.Vĳayan et al. 2021;Kannan et al. 2022).
In Cullen et al. (2023) we provided a fit to  −  UV for a sample of 61 galaxies at ⟨⟩ = 10.Due to the small sample size in Cullen et al. (2023), in this initial fit we included all galaxies across the full redshift range of the sample from  = 8 to  = 16 4 .We found evidence for a  −  UV relation with d/d UV = −0.17± 0.05.This slope was consistent with the slope measured for large samples at lower redshift, for example, the  = 5 relations of Bouwens et al. (2014) and Rogers et al. (2014) which have d/d UV = −0.12± 0.02 and d/d UV = −0.14 ± 0.02 respectively.However, we found that the normalisation for the relation was lower at  ≃ 10, such that galaxies are on average bluer by Δ = −0.4compared to  ≃ 5.
With our enlarged sample, we are now in a position to investigate the redshift evolution of the  −  UV relation at  ≃ 8 − 12.For the analysis in this section, we make use of the combined sample, which includes the COSMOS/UltraVISTA galaxies at 7.5 <  < 10.These ground-based candidates trace the bright end of the UV luminosity  distribution ( UV < −21) and are therefore crucial in providing a large dynamic range in  UV (see Fig. 1).To investigate the evolution of the  −  UV relation with redsfhit, we split our sample into three bins of redshift; for each redshift bin, we then calculated the inversevariance weighted mean  in three bins of  UV .The resulting values of ⟨⟩, ⟨ UV ⟩, and ⟨⟩ are given in Table 4.
In practice, only the redshift bin 7.5 <  < 10 has a sufficient dynamic range in  UV (from  UV = −22.7 to  UV = −17.2;Fig. 7) to allow an accurate estimate of the slope of the  −  UV relation.In this redshift range, we find evidence for a  UV dependence, with ⟨⟩ evolving from ⟨⟩ = −1.79± 0.11 at the brightest magnitudes (median  UV = −21.4) to ⟨⟩ = −2.27± 0.08 at the faintest magnitudes (median  UV = −18.2).The best fitting colour-magnitude relation to the individual candidates (blue data points in Fig. 7) is where we have accounted for the uncertainties in both  and  UV .
In this formulation, the intercept, ( UV = −19) = −2.19± 0.05, represents the typical value of  at  UV = −19.Our new best-fit relation is slightly shallower than the fit in Cullen et al. (2023) but is fully consistent within the uncertainties.We note that the median redshift of galaxies at the bright end ( UV < −21) is lower than the median redshift in the fainter bins by Δ ≃ 1, which may introduce a bias if there is a significant redshift evolution in ⟨⟩ at fixed  UV between  = 8 and  = 9.The clear lack of galaxies at  ≥ 9 and  UV < −21 with JWST photometry is something that will hopefully be addressed with current and future wide-area surveys (e.g.Casey et al. 2023;Franco et al. 2023).
In the two higher redshift bins we do not attempt to constrain  −  UV due to the lack of galaxies with  UV < −21.0.In both redshift bins, the data do not show any clear evidence for a  −  UV relation between  UV = −20.0 and  UV = −18.0.Due to the lack of sufficient dynamic range  UV at these redshifts, we have assumed a fixed slope of d/d UV = −0.15 in these redshifts bins and fitted only the normalisation.This approach yields ( UV = −19) = −2.45± 0.06 at 10 <  < 11 and ( UV = −19) = −2.64 ± 0.06 at 11 <  < 12 (where the uncertainties do not account for the slope uncertainty).The two fits are shown in Fig. 7.The resulting best-fitting  −  UV parameters are given in Table 5.These fits again highlight our main result: at fixed  UV , the typical value of ⟨⟩ is evolving rapidly with redshift between  = 9 and  = 12, reaching ⟨⟩ ≃ −2.6 at the highest redshifts.
In summary, we find evidence (5) for a  −  UV relation in our sample at ⟨⟩ ≃ 9 with a slope of d/d UV = −0.15± 0.03.This slope is fully consistent with the slope derived at  = 5 (e.g.Bouwens et al. 2014;Rogers et al. 2014).Our results currently suggest that the slope of the  −  UV relation is already established within the first ≃ 500 Myr of cosmic time and does not evolve strongly thereafter.At  > 10 we lack the dynamic range in UV to constrain the slope, however, assuming a fixed d/d UV = −0.15,we find that the normalisation evolves rapidly from ( UV = −19) = −2.19± 0.05 at ⟨⟩ ≃ 9 to ( UV = −19) = −2.64 ± 0.06 at ⟨⟩ ≃ 11.

Evidence for piecewise linear relation?
We note that some studies have presented evidence for a change in the slope of the −  UV relation at faint magnitudes such that d/d UV becomes shallower at the faint end.For example, Bouwens et al. (2014) presented tentative evidence for a piecewise-linear  −  UV relation  ≃ 4 − 6, with the change in slope occurring at  UV ≳ −19.On the other hand, Rogers et al. (2014) find no significant evidence for a nonlinear  −  UV relation  = 5.From Fig. 7 it can be seen that our data appear qualitatively consistent with a flattening at faint magnitudes.However, we currently lack the dynamic range in  UV to reliably confirm this trend.Ultimately, larger sample sizes and an extension to fainter  UV are needed to further explore this possibility at  > 9.In the following discussion, we compare our measurements with those in the literature that assume a simple linear form of the  −  UV relation.

Comparison to the literature
In Fig. 8 2014) is ⟨d/d UV ⟩ = −0.12.We also show the wide-area (≃ 0.8 deg 2 ) analysis at  = 5 of Rogers et al. (2014) which, given the total sample size and area, is the benchmark constraint on  −  UV at this redshift.Rogers et al. (2014) find d/d UV = −0.12± 0.02 and ( UV = −19) = −1.93 ± 0.03, in excellent agreement with the  = 5 constraints of Bouwens et al. (2014).Therefore, we consider the  = 5 constraints a robust anchor with which to compare the data at higher redshifts.We also show recent determinations of the −  UV relation across the redshift range 6 <  < 12 from Topping et al. (2024).Their analysis is based on deep JWST/NIRCam observations (taken as part of the JADES survey) and is therefore directly comparable to ours.In their two bins at  > 9 ( = 9.4 and  = 12.0), Topping et al. (2024) (Leitherer et al. 1999;Götberg et al. 2018; left-hand panel), bpassv2.3single-burst models including -enhanced abundance ratios (Byrne et al. 2022; centre panel) and Flexible Stellar Population Synthesis (fsps) single-burst models (Conroy et al. 2009;Conroy & Gunn 2010;Byler et al. 2017;right-hand panel).In each panel, the dashed line shows the UV slopes for pure stellar continuum models.Each line is colour-coded by stellar metallicity, as indicated in the legend.The solid lines show the UV slopes when accounting for nebular continuum emission assuming an escape fraction of 0 per cent (i.e. the maximum nebular continuum strength).The shaded yellow area shows our constraint of ⟨⟩ = 2.60 ± 0.05 at  > 10.5 from the fit described in Section 3.1.All models suggest that the population average value we measure at  > 10.5 is consistent with pure stellar plus nebular emission in the absence of dust.
an essentially flat slope with d/d UV = −0.06± 0.05 (at  = 12.0 the slope is fixed to the  = 9.4 value).At  = 9.4, this is shallower than our determination of d/d UV = −0.15± 0.03.The difference may be attributable in part to the lack of bright-end constraints in Topping et al. (2024): their brightest bin has an absolute UV magnitude of  UV = −20.0,whereas our sample probes to  UV = −22.7 thanks to the inclusion of the COSMOS/UltraVISTA sample (Fig. 1).However, in general, the tension is clearly not significant, and the two constraints are consistent within 2.At  ≃ 12 our results agree fairly well within the uncertainties.Crucially, Topping et al. (2024) also finds that by  = 12 the galaxy population is uniformly extremely blue.
Although there are tensions between the various studies, these differences are not highly significant (i.e.none of current constraints at fixed redshift are incompatible at the ≥ 3 level).The general picture discernible from Fig. 8 is of a relatively shallow  −  UV relation that is in place from  ≃ 10, with a fixed (or at least slowly evolving) slope, probably in the range −0.2 < d/d UV < −0.1.At  > 10 more data are needed to robustly constrain this slope, although current data are compatible with a slowly-/non-evolving d/d UV scenario.In contrast, the normalisation of the relation is clearly changing, with galaxies at fixed  UV becoming bluer towards higher redshift.Based on our results, the evolution is gradual below  ≃ 10, with ( UV = −19) evolving from −2.2 at  = 10 to −1.9 at  = 5.At  > 10, the evolution becomes more rapid (especially when framed in terms of cosmic time), reaching ( UV = −19) ≃ −2.6 by  = 12.The analysis of Topping et al. (2024) suggests a more gradual evolution, but still suggests extremely blue UV colours at the earliest cosmic epochs.

DISCUSSION
We have presented an investigation of the UV continuum slopes of a sample of 172 galaxy candidates at  ≥ 7.5, with the aim of understanding the dependence of  on redshift and  UV .The primary focus of our analysis has been a new wide-area sample of 121 galaxies at  ≥ 9 selected from 15 public JWST/NIRCam imaging datasets (covering an on-sky area of ≃ 320 arcmin 2 ).The main result of this work is that we find a strong redshift evolution of  in our wide-area sample, with a rapid transition from ⟨⟩ ≃ −2.2 at  ≃ 9 to ⟨⟩ ≃ −2.6 at  ≃ 11 (Figs. 4 and 5).The population average value of  at  > 10.5 is consistent with expectations for extremely dust-poor, even dust-free, stellar populations(see discussion below).We then investigated the  −  UV relation using our full sample of 172 galaxies, which included the bright ( UV < −21) groundbased candidates at ⟨⟩ ≃ 8.5 from COSMOS/UltraVISTA.We find evidence (5) for a  −  UV relation at  ≃ 9 with a slope of d/d UV = −0.15± 0.03, similar to the slope observed at lower redshifts (Figs. 7 and 8).At  > 10, our data are not sufficient to accurately constrain the  −  UV slope due to a lack of bright galaxies; however, our analysis confirms a rapid evolution in the normalisation of the relation across the redshift range 9 <  < 12.In this section, we discuss the dust-free galaxy interpretation and consider the implications for dust formation in the early Universe and cosmic reionisation.

A transition to dust-free star-formation at z ≳ 10?
The bluest/steepest possible value of  is set by the intrinsic stellar spectrum of young metal-poor galaxies.At the youngest ages (≲ 1 Myr) and the lowest metallicities (≲ 1 per cent solar), the UV slope can reach values as blue at  ≃ −3 (e.g.Bouwens et al. 2010;Robertson et al. 2010;Stanway et al. 2016;Topping et al. 2022).However, the presence of ionised gas in galaxies will result in continuum emission as a consequence of free-free, free-bound, and two-photon processes (i.e. the nebular continuum).This nebular continuum acts to redden the slope of the UV continuum (e.g.Byler et al. 2017).Models suggest that the bluest expected UV slopes when accounting for nebular continuum emission are  ≃ −2.6 (e.g.Stanway et al. 2016;Topping et al. 2022).
We reproduce this result for a selection of stellar population synthesis models in Fig. 9 where we plot  as a function of the stellar population age.The stellar populations considered are starburst99 continuous star formation models including a contribution from stripped stars (Leitherer et al. 1999;Götberg et al. 2018), bpassv2.3singleburst models including -enhanced abundance ratios (Byrne et al. 2022)   and the SMC extinction curve (Gordon et al. 2003).These models indicate that to achieve an observed UV continuum slope of  obs = −2.6 including dust would require a counterintuitive scenario in which a young burst, with a large (≳ 50 per cent) escape fraction is at the same time experiencing ≳ 0.2 magnitudes of UV attenuation (see text for discussion).
Appendix B. It can be seen that although variations exist between each model, the overall picture is consistent: when an escape fraction (  esc ) of 0 per cent is assumed (i.e.all of the ionising photons emitted by the stellar population are converted into nebular emission), the bluest expected UV continuum slope is  ≃ −2.6.Moreover, this occurs only for the youngest ages ≲ 30 Myr.If a nonzero  esc is assumed, then bluer slopes are possible, with a minimum value of  ≃ −3, but this occurs at ages of ≲ 3 Myr and for escape fractions close to 100 per cent (i.e.pure stellar emission).
Another notable feature of Fig. 9 is that the youngest, most metalpoor stellar populations produce the strongest nebular continuum emission due to their harder ionising spectra.In some cases, the strong nebular emission can entirely compensate for their intrinsically bluer stellar UV slopes.For example, considering the FSPS models (right-hand panel of Fig. 9), the UV continuum slopes, including nebular continuum emission, are actually redder at the youngest ages.From this we can conclude that invoking more extreme stellar populations with harder ionising spectra (i.e.top-heavy IMFs, Pop-III stars) is unlikely to result in UV slopes intrinsically bluer than  = −2.6 when the effects of nebular continuum emission are included (see e.g.Cameron et al. 2023).
The model predictions in Fig. 9 imply that a stellar population with an observed UV continuum slope of  obs ≃ −2.6 is almost certainly experiencing negligible (essentially zero) dust attenuation, regardless of age.For a stellar population to be dust-attenuated yet still have  obs = −2.6 would require a counter-intuitive scenario in which a significant fraction of ionising continuum photons are escaping into the IGM in the presence of dust.We highlight this argument in Fig. 10, where we show  versus stellar population age for a bpassv2.3burst model across a range of  esc values.Even for extreme escape fractions of  esc = 50 per cent, the intrinsic UV continuum slope only falls below  ≃ −2.6 at the very youngest ages (i.e.< 5 Myr).To achieve  = −2.6 including dust would therefore require a very young, dominant, burst with a large (≳ 50 per cent) escape fraction that at the same time experiences ≳ 0.2 magnitudes of UV attenuation.Although we cannot completely rule such a scenario out (i.e. for very specific star/gas/dust geometries), the degree of fine-tuning required renders it an unlikely explanation for population-average values.
Indeed, it is arguably more likely that the stellar population ages at  > 10.5 are ≳ 30 Myr, such that, on average, these galaxies will have intrinsic UV continuum slopes of  int > −2.6 (Fig. 9).For example, Cullen et al. (2017) analysed the stellar populations of simulated galaxies at  = 5 and found that the typical intrinsic UV slopes assuming  esc = 0 were  ≃ −2.4 for galaxies with light-weighted stellar ages between 20 Myr and ≃ 100 Myr.Lightweighted ages of up to 100 Myr are clearly plausible at  = 11 as this would correspond to a formation redshift of  = 13 and galaxies at this redshift have already been spectroscopically confirmed by JWST (Curtis-Lake et al. 2023).Therefore, it is possible and even likely that the typical intrinsic slopes of the  > 10.5 galaxies in our sample are redder than  = −2.6 (for  esc = 0).In this case, an observed UV continuum slope of  obs = −2.6 indicates both negligible dust attenuation and and a nonzero escape fraction of ionising photons.We discuss implications for the escape fraction and reionisation in more detail below.
Finally, it is worth acknowledging that our ⟨⟩ constraints are formally consistent with a value as red as ⟨⟩ = −2.35within the uncertainties (5).Assuming an intrinsic slope of  int = −2.6 (i.e. an extremely young burst), this value would correspond to an upper limit of  UV ≲ 0.2 − 0.5 depending on the assumed attenuation curve (McLure et al. 2018).However, for older and/or composite stellar populations (with  int ≃ −2.4), the implied upper limit on UV attenuation would be closer to  UV ≲ 0.1 − 0.2.Nevertheless, our formal best estimates are most consistent with negligible, essentially zero, dust attenuation.Below we consider some physical interpretations of dust-free systems in the young Universe.

The physical motivation for dust-free galaxies
The above comparison with theoretical UV spectra models suggests that the average  value we observe at  ≃ 11 is consistent with the expectation for dust-free stellar populations.It is important to ask whether such a scenario is physically plausible.Dust is formed in the ejecta of core-collapse supernovae (e.g.Todini & Ferrara 2001), and as a result galaxies might be expected to accumulate dust immediately after the onset of star formation.In this case, a mechanism for either ejecting or destroying dust formed in supernovae is required to explain the absence of dust in star-forming galaxies.Ferrara et al. (2023) and Ziparo et al. (2023) propose one potential scenario in which intense UV radiation pressure from ongoing star formation in these early galaxies drives dust out of the interstellar medium (see also Nath et al. 2023;Tsuna et al. 2023).Most recently, Ferrara (2023) have postulated that above  ≃ 10 the specific star-formation rate of galaxies crosses a threshold above which powerful radiation-driven outflows are capable of clearing dust from the ISM.Intriguingly, this is the same redshift threshold above which we begin to see uniformly extremely blue colours in our sample.
An alternative hypothesis, proposed by Jaacks et al. (2018), is that the first generation of star formation (i.e.Pop-III stars) produces negligible amounts of dust.In this scenario the stars at  > 10.5 are either forming from recently enriched but dust-free Pop-III gas, or alternatively we are seeing the effect of a significant fraction of current Pop-III star formation in our  > 10.5 sample.Indeed, Jaacks et al. (2018) estimate  = −2.5 ± 0.07 for Pop-III SEDs.We note, however, that, other than their blue UV continuum slopes, the NIRSpec/PRISM spectra shown in Fig. 2 do not show any obvious Pop-III signatures in their FUV spectra (e.g.strong broad He ii emission).Finally, it is also possible that dust is more efficiently destroyed in the earliest star-forming systems.Current estimates of dust grain destruction rates in the ISM vary by up to an order of magnitude due to different assumptions regarding dust mircrophysics (e.g.grain-grain collisions and shattering; e.g.Kirchschlager et al. 2022).It is possible that specific physical conditions at  > 10 could result in enhanced dust destruction (e.g.enhanced surface densities of SFR), although direct evidence is currently lacking.The dust destruction scenario is supported by the theoretical analysis of Esmerian & Gnedin (2023), who used a cosmological fluid-dynamical simulation of galaxies within the first 1.2 Gyr ( > 5) to investigate the effect of different dust models.A comparison of our UV slope measurements up to  ≃ 12 with their model favours enhanced dust destruction with destruction rates elevated by an order of magnitude or more relative to default assumptions.There are therefore a number of possible physical processes which can explain dust-free galaxies.However, our current data cannot distinguish between them.Future observations, in particular deep NIRSpec/PRISM spectroscopy and improved constraints on the redshift evolution of ⟨⟩, will be crucial in this regard.

Direct constraints on dust emission at 𝑧 > 10
Direct detection of dust emission in the rest-frame far-infrared at  > 10 would represent a clear refutation of a dust-free star formation scenario.However, it is worth noting that to date, all deep ALMA observations aimed at detecting dust continuum emission at these redshifts have been unsuccessful (e.g.Fujimoto et al. 2023b;Bakx et al. 2023;Popping 2023;Yoon et al. 2023;Fudamoto et al. 2024).Therefore, currently these observations are fully consistent with the results presented here.Nevertheless, it is clear that further deep ALMA observations (enabling, for example, stacking of multiple targets) will be crucial for robustly testing the inferences that are made from rest-frame UV studies.

The implications for cosmic reionisation
As discussed above, while a UV continuum slope of  ≃ −2.6 clearly indicates negligible dust attenuation, it also implies that the ionising photon escape fraction (  esc ) is likely to be high.In fact, a direct connection between observed  and  esc has been established up to  ≃ 3 (e.g.Chisholm et al. 2022;Begley et al. 2022;Kim et al. 2023).It is likely that this connection extends to higher redshifts.At  ≃ 7.5, Topping et al. (2024) find that ultra-blue UV continuum  2022) using (primarily) data from the LzLCS survey (Flury et al. 2022).Open circular data points show galaxies with unphysical values of  esc > 1 based on this calibration (i.e. the galaxies in our sample that have been scattered to extremely blue  values).The square data points in red and blue show the population average value of  esc above and below  = 10.5, respectively, using the best-fitting  values from Fig. 5.The red and blue horizontal lines show  esc = 0.05 and  esc = 0.20, respectively.Taking the predictions of the low-redshift calibration at face value implies that a substantial fraction of ionising photons (≃ 20 per cent) are escaping into the IGM at  > 10.5.
slopes are associated with a lack of nebular emission line signatures in the rest-frame optical, again suggesting that bluer UV colours are associated with higher  esc .Detailed cosmological radiative transfer simulations also predict that  is a key indicator of  esc at  > 4 (Choustikov et al. 2024).Using a sample 89 galaxies at  ≃ 0.3 with direct LyC measurements from HST/COS spectroscopy, Chisholm et al. (2022) presented a strong (5.7)correlation between  and  esc .The majority of the sample (66 galaxies) was drawn from the Low-redshift Lyman Continuum Survey (LzLCS; Flury et al. 2022)  Interestingly, if ⟨  esc ⟩ ≃ 0.2, then the galaxy population that has already been observed with JWST is likely able to provide enough ionising photons to drive reionisation at  = 11.In Fig. 12 2024).The dashed horizontal lines show the equilibrium value of  ion needed to maintain an ionised IGM fraction of  H ii (as indicated by the labels).The red curve shows  ion as a function of  UV,lim assuming log(  ion /Hz erg −1 ) = 25.8 (i.e. at the upper end of  ion estimates for young metal-poor galaxies; e.g.Tang et al. 2019).The curve turns from solid to dashed below the limiting magnitude of current LF constraints ( UV ≃ −18).At this efficiency, the already observed population of galaxies (i.e. UV ≲ −18; vertical black line) would be supplying enough ionising photons to maintain  H ii > 0.1.The orange and blue curves represent log(  ion /Hz erg −1 ) = 25.5 and log(  ion /Hz erg −1 ) = 25.2, which may be more representative of the typical value during reionisation (e.g.Bouwens et al. 2016;Matthee et al. 2017).At the lowest estimated ionising efficiencies, the observed population would be sufficient to maintain  H ii ≃ 0.05.Taking into account the faint, unseen, galaxy population at  UV > −18 implies relatively high (likely ≳ 5 per cent) ionised fractions at  = 11.
for three values of  ion (i.e. the production efficiency of LyC photons in units of Hz/erg) that bracket the typical range of values inferred for young, low-metallicity galaxies at high redshifts (i.e.log( ion /Hz erg −1 ) = 25.2 − 25.8; e.g.Bouwens et al. 2016;Atek et al. 2024;Tang et al. 2019Tang et al. , 2023)).The ionising photon output needed to maintain an ionised IGM fraction of  H ii is given by the prescription of Madau et al. (1999): where  H ii is the IGM clumping factor which we set following the prescription of Shull et al. (2012).
Setting  H ii = 1 in equation 5 represents the limiting case of maintaining a fully ionised Universe at redshift .However, at  = 11 the Universe is expected to be only partially ionised.Models assuming an early and slow reionisation typically have  H ii ≃ 0.2 (e.g.Finkelstein et al. 2019) while late and rapid models have  H ii ≃ 0.01 (e.g.Naidu et al. 2020).It can be seen from Fig. 12 that, based on the above assumptions, the currently observed galaxy population (i.e. UV ≲ −18) is likely providing enough ionising photons to main-tain  H ii ≳ 0.05.Accounting for the population of fainter galaxies, examples of which have been uncovered down to  UV ≃ −15 (e.g.Atek et al. 2024), implies that the full galaxy population could in principle be supplying a sufficient number of photons to maintain an ionisation fraction of ≳ 20 per cent at  = 11 (the exact value is strongly dependent on the uncertain faint-end slope of the UV LF).Despite systematic uncertainties, it seems probable that the large abundance of blue galaxies at  > 10 unveiled by early JWST studies disfavours low ionised IGM fractions of ≲ 1 per cent.

SUMMARY AND CONCLUSIONS
We have measured the rest-frame ultraviolet (UV) continuum slopes () for a sample of 172 galaxy candidates at ⟨⟩ = 10.5.The majority of these candidates (≃ 70 per cent) are drawn from our new wide-area JWST sample with absolute UV magnitudes in the range −20.5 <  UV < −17 (Section 2.1).This new sample has allowed us to improve on our previous work (Cullen et al. 2023) in terms of sample size, on-sky area, and median imaging depth.Using a robust power-law fitting technique validated against available spectroscopy (Section 2.3), we report precise estimates of ⟨⟩ at  > 10.The main aim of this analysis is to use the average values of  for these galaxies to place constraints on the typical dust obscuration experienced at the earliest cosmic epochs.Our main results can be summarised as follows: (i) We measure an inverse-variance weighted mean value of ⟨⟩ = −2.37 ± 0.03 for our primary wide-area sample at ⟨⟩ = 10.6.Incorporating candidates from Cullen et al. (2023) − that include UV-bright ground-based COSMOS/UltraVISTA sources at lower redshift (with ⟨⟩ ≃ 8; Fig. 1) − yields a slightly redder value of ⟨⟩ = −2.32 ± 0.03 and a lower average redshift of ⟨⟩ = 10.3.Uniformly, therefore, our sample displays blue UV continuum slopes indicative of dust-poor galaxies at these redshifts.However, a value of ⟨⟩ ≃ −2.37 is not inconsistent with moderate amounts of dust attenuation, nor more extreme than blue objects observed in the local Universe (e.g.NGC 1705;  = −2.46).
(ii) Although the median value of the UV slopes of the widearea sample is ⟨⟩ ≃ −2.37, we find evidence for a steep  −  trend between  = 9 and  = 12 (Fig. 4).The average value of  evolves from ⟨⟩ = −2.17± 0.06 at  ≃ 9.5 to ⟨⟩ = −2.59± 0.06 at  ≃ 11.5.This represents a rapid evolution of Δ⟨⟩ ≃ 0.4 over ≃ 100 Myr of cosmic time.A slope of  = −2.6 represents the extreme end of the local distribution (e.g.Chisholm et al. 2022) but appears to be typical for galaxies at the earliest cosmic epochs.Such steep UV slopes are consistent with pure stellar plus nebular continuum emission that is not attenuated by dust (Fig. 9).
(iii) Fitting a step function model to the  −  data we find that the galaxy population is becoming uniformly extremely blue at  ≳ 10.5 (Fig. 5).Below this transition redshift, our sample displays ⟨⟩ = −2.22 ± 0.05 while at  ≳ 10.5 the typical value is ⟨⟩ = −2.60 ± 0.05.We have verified that this difference is not caused by measurement biases or selection effects.
(iv) In the redshift range 8 <  < 10 we use the large dynamic range in  UV enabled by the inclusion of COSMOS/UltraVISTA candidates to investigate the relationship between  and absolute UV magnitude ( UV ).We find evidence for a  −  UV relation such that brighter galaxies display redder UV slopes (Fig. 7).This is consistent with a picture in which the brighter (and more massive) galaxies at  ≃ 9 are more obscured by dust.Fitting the  −  UV at these redshifts, we find d/d UV = −0.15± 0.03.This slope is consistent with values derived a lower redshifts (down to  = 2; e.g.Bouwens et al. 2014).Therefore, our results suggest that a  −  UV relation has been in place since  ≃ 10 with a slope that does not strongly evolve with redshift.However, the normalisation of the relation is clearly evolving, with ( UV = −19) ≃ −1.8 at  ≃ 5 compared to ( UV = −19) ≃ −2.2 at  ≃ 9.At fixed  UV , galaxies are less dust-enriched (and metal-enriched) than at earlier cosmic epochs.
(v) At  > 10 our data do not cover a sufficient dynamic range in  UV to robustly constrain the  −  UV relation (Fig. 7).However, the data at these redshifts remain consistent with a relatively shallow  −  UV slope and a continued evolution to bluer values of ( UV = −19).Our results suggest that the normalisation evolves from ( UV = −19) = −2.2 at  ≃ 9 to ( UV = −19) = −2.6 at  ≃ 11.At  > 10, wider-area surveys are required to probe the bright end of the UV LF and determine whether the most luminous galaxies ( UV ≲ −20) show evidence for dust-attenuated UV SEDs.
The primary new result of this work is the appearance of a uniformly extremely dust-poor, perhaps even dust-free, galaxy population at  > 10.5, with average UV continuum slopes of ⟨⟩ ≃ −2.6.Similarly blue UV slopes at these redshifts have also been reported in other independent analyses (Austin et al. 2023;Topping et al. 2022;Morales et al. 2024).Competing theoretical explanations for dust-free galaxies exist, but these different physical processes cannot be discriminated with our current data.Nevertheless, these results place important constraints on the origin of dust and metals in the first galaxies.Furthermore, our results imply that the ISM conditions in galaxies at  > 10 favour a significant escape fraction of ionising photons, and that the already observed population of galaxies at these redshifts (i.e. with  UV ≲ −18) is likely capable of maintaining ionised IGM fractions of ≳ 5 per cent.
Table A1.Overview of the different datasets utilised in this study.The first three columns list the name of the field/survey, the JWST proposal ID and the name of the principal investigator.Columns four and five list references to the survey paper and ancillary HST data (where applicable).Column six lists references for the lensing maps used to correct the photometry for the effect of gravitational lensing (where applicable), with the following key: Z15 = Zitrin et al. ( 2015    > 9 the full rest-frame UV SED up to  rest ≃ 3000Å is sampled by the JWST/NIRCam filters 8 , the moderate increase in selection efficiency towards redder  values will apply at all redshifts.From Fig C1 we can conclude that the evolution in ⟨⟩ we observe cannot be a result of an increasing efficiency in selecting bluer galaxies at higher redshift as the relative selection efficiency as a function of  is the same at  = 9 and  = 11.These results highlight a more general point: galaxy selection up to  ≃ 15 with JWST should not be biased towards the intrinsically bluer sources because deep JWST/NIRCam imaging allows us to select across their full rest-frame UV spectrum 9 .
8 Our selection filters trace the rest-frame UV SED up to  rest ≃ 2700Å at  = 9 and  rest ≃ 2300Å at  = 11. 9We note that the increased efficiency in selecting redder objects discussed here is subtly different to the often-discussed  'blue-bias' (e.g.Bouwens et al. 2010;Dunlop et al. 2012;Rogers et al. 2013;Cullen et al. 2023).The blue bias occurs for faint galaxies that are selected via the filter closest to the Lyman break (which will typically be the case for galaxies with  ≲ −2.0, As a second test, we estimated the source recovery fraction and recovered  when shifting the best-fitting SEDs of our 9 <  < 10.5 sources to higher redshift and reapplying our selection criteria.The inverse-variance weighted mean  of these sources is ⟨⟩ = −2.22 with a mean redshift of ⟨⟩ = 9.8.The aim of the test was to discover whether the recovered galaxies were biased towards the bluer or redder objects, and whether there was any bias in the estimated  for the recovered sources.We redshifted each individual 9 <  < 10.5 galaxy by Δ = 1.2 so that the median redshift of the artificially redshifted sample became ⟨⟩ = 11.We then computed new model i.e. a falling   SED).For low S/N objects close to the detection threshold, this selection will favour objects whose photometry has been upscattered in the short-wavelength detection band, and these objects will appear bluer than they actually are.However, because we have restricted our sample to high S/N detections (8), and allowed for selection in any band redward of the Lyman break, the sample average  values for our current sample is not affected strongly by this bias (see discussion in Section 3.1.2).

ID
RA Dec  phot  UV that ≃ 20 per cent drop below our selection thresholds.It can be seen that, on average, the redder sources are more likely to be recovered, and as a result the average  of the recovered sources is biased slightly red (by Δ = 0.04).This result is consistent with the simulations discussed above (i.e. at any given redshift, redder sources are slightly favoured by our selection; Fig. C1).Crucially, for the sources that are recovered, there is no bias in the individual galaxy  estimates at the new redshift (inset panel in Fig. C2).This test conclusively rules out a scenario in which the blue ⟨⟩ we observe at  > 10.5 is the result of a bias favouring the selection of bluer sources.Indeed, this analysis suggests that if the  ≃ 9.5 galaxy population were transplanted to  ≃ 11, the redder sources would preferably be recovered.The fact that we do not see these sources in our data is therefore indicative of a real evolution in the shape of the UV continuum SED with redshift.We also ran this test in reverse, taking the 10.5 <  < 12.0 sources (with ⟨⟩ = −2.60 and ⟨⟩ = 11.4) and blueshifting them to lower redshift by Δ = −1.9so that the median redshift became ⟨⟩ = 9.5.In this case, we find a high average recovery fraction of 96 per cent across the 500 iterations (as expected).Typically, one or two sources per iteration failed to fulfil the selection criteria due to the random photometric perturbations, and these sources were always amongst the bluest in the sample (with  < −2.7).However, because most of the sources remain within our selection criteria, the mean  of the population was accurately recovered (the inverse variance-weighted mean value for the recovered sources was ⟨⟩ = −2.61± 0.01).From this test, we can further conclude that if the blue population we observe at  > 10.5 also existed at  ≃ 9.5, these blue sources would be present in our sample.

Figure 1 .
Figure 1.The distribution of  UV and photometric redshift for our primary wide-area JWST sample (blue points) and the sample from Cullen et al. (2023) comprised of ground-based COSMOS/UltraVISTA and JWST-selected candidates (red and yellow points respectively).The combined sample covers the redshift range 7.6 <  phot < 15.9 (i.e.spanning the time period ≃ 250 − 700 Myr after the Big Bang) and the  UV range −22.7 ≤  UV ≤ −17.2 (i.e. a factor of ≃ 150 in UV luminosity).
Figure2.A comparison of UV continuum slope fits to photometry ( phot ) and spectroscopy ( spec ) for the nine galaxies in our sample with JWST/NIRCam and NIRSpec/PRISM observations.In the first nine panels, the black curve is the NIRSpec/PRISM spectrum, with the grey shading representing the ±1 error spectrum.The square data points show the galaxy photometry, with all photometry at  ≥ 0.9 m coming from JWST/NIRCam.The solid red and dashed blue curves are the best-fitting models for the spectrum and photometry, respectively.In the majority of cases, the photometric and spectroscopic fits are fully consistent; however, for some galaxies, small differences in the best-fitting redshift and  are evident.The tenth panel (lower left) shows the comparison between  spec and  phot .The grey circular points show the nine individual galaxies and the red square point is the inverse variance-weighted mean.In general, we find very good agreement between the two estimates with ⟨ spec ⟩ = −2.33 ± 0.02 and ⟨ phot ⟩ = −2.39 ± 0.07.

Figure 4 .
Figure 4. Plots of UV continuum slope () versus redshift () at  > 8.The left-hand panel shows the individual measurements for our primary wide-area JWST sample (grey circular data points) and the inverse-variance weighted mean values in three redshift bins (yellow square data points).The solid red line shows the best-fitting  −  relation fitted to the individual measurements, which has a slope of d/d = −0.28 ± 0.05.The light-red shaded region shows the 95 per cent confidence interval.The horizontal blue shading represents  in the range −2.6 to −2.4, approximately the minimum value expected for a standard stellar population assuming no dust and a maximum contribution of the nebular continuum (see Section 4.1 and Cullen et al. 2017;Reddy et al. 2018).In the highest redshift bin (11 <  < 12), the population mean value approaches this limit, indicating galaxies unattenuated by dust.In the right-hand panel, we show the same  −  relation for our combined sample (i.e., including the galaxy candidates fromCullen et al. 2023, as indicated in the legend).For clarity, we omit error bars in this panel.At  > 9, the JWST-selected galaxies ofCullen et al. (2023) are consistent with the same trend.Our data show evidence for a steep decline in  with redshift at  > 9, indicating that galaxies are growing progressively more dust and metal poor at early cosmic times.By  ≃ 11, the population-average  is consistent with negligible, even zero, dust attenuation.

Figure 5 .
Figure5.Plot of  versus redshift for our primary wide-area JWST sample at  ≥ 9 (coloured data points) highlighting the transition to dust-free stellar populations at  ≃ 10.5.Below  ≃ 10.5 (black dotted vertical line), the typical value of the UV slope is ⟨⟩ = −2.22(red solid line and red data points); this value is blue, but not inconsistent with moderate amounts of dust attenuation ( UV ≃ 0.5 − 1, assuming  int = −2.6).Above this redshift, however, we find a typical value of ⟨⟩ = −2.60 (blue solid line and blue data points), which is at the extreme blue end of any galaxy observed locally (e.g.Chisholm et al. 2022).Based on stellar population modelling of young, metal-poor galaxies,  ≃ −2.6 is the intrinsic UV slope expected assuming a maximum nebular continuum contribution (i.e. esc = 0) but an absence of dust (see Section 4.1 and discussions inRobertson et al. 2010;Cullen et al. 2017 andReddy et al. 2018).In the inset panels, we show the distribution of  (left) and  UV (right) above (red) and below (blue)  ≃ 10.5.Due to the similar distributions in  UV (both have ⟨  UV ⟩ = −19.2),we do not expect the difference in the  distribution to be strongly biased by differences in galaxy luminosity.Our results suggest uniformly extremely blue  values at  > 10.5.

Figure 6 .
Figure 6.Plot of the UV continuum slope bias (Δ =  recovered −  input ) as a function of  for of 20,000 simulated galaxies at 9 <  < 12 with  int = −2.2(see Section 3.1.2for a description of the simulations).Each data point represents an individual galaxy from the simulation and is colourcoded according to the difference between the recovered and input redshifts ( recovered −  input ).The solid black line shows the running inverse-variance weighted mean value of Δ as a function redshift.In general, the bias is negligible at all redshifts (Δ ≃ −0.01).

Figure 7 .
Figure 7. Plots of   UV for the combined sample in three bins of redshift.In each panel, the redshift interval is shown in the top left-hand corner, and coloured circular data points are the individual galaxy candidates in our sample.The black square data points are the inverse-variance weighted values of ⟨⟩ in the bins of  UV given in Table4.In the left-hand panel the black solid line is the best-fitting  −  UV relation (Equation4) with the light-grey shaded region showing the 95 per cent confidence interval around our best-fitting relation.We find evidence for a  −  UV relation with best-fitting slope of d/d UV = −0.15± 0.03 (i.e.5 significance).In the centre and right-hand panels, the black solid line shows the best-fitting  −  UV relations where we have fixed the slope to be the same as the slope determined at  < 10 (see text for details).In the centre and right-hand panels, the black dashed line shows the best-fitting  −  UV relation at  < 10 to highlight the redshift evolution of the normalisation.The best-fitting  −  UV parameters are given in Table5.

4
The  = 16 candidate has subsequently been shown to be a  ≃ 5 interloper (Arrabal Haro et al. 2023b), although the inclusion of this object does not affect the  −  UV fit in Cullen et al. (2023).

Figure 8 .
Figure 8. Constraints on the  −  UV relation as a function of redshift from this work (blue circles) and three literature studies: Bouwens et al. (2014) (red squares), Rogers et al. (2014) (black diamond) and Topping et al. (2024) (yellow triangles).The upper panel shows the evolution of the slope of the relation (d/d UV ), which remains relatively constant with redshift, with constraints typically in the range −0.2 < d/d UV < −0.1.The open (i.e.unfilled) data points in the upper panel highlight cases in which the slope has been fixed when fitting (as is the case for the  > 10 relations in this work; see text for discussion).The lower panel shows the evolution of normalisation,  (  UV = −19), which clearly evolves such that galaxies at higher redshifts have bluer UV slopes consistent with being typically younger and less dust and metal enriched.At  ≃ 12, the current constraints (our workand Topping et al. 2024) suggest that the galaxy population (at 20.5 ≲  UV ≲ −18) is on average extremely blue, approaching the dust-free limit of  ≃ −2.6.

Figure 10 .
Figure10.The intrinsic UV continuum slope as a function of stellar population age for a bpassv2.3burst model(Byrne et al. 2022) with 1 percent solar metallicity.Each black line shows the relation for a different assumed value of  esc as indicated in the legend.The shaded yellow area shows our constraint of ⟨⟩ = 2.60 ± 0.05 at  > 10.5 from the model fit described in Section 3.1.The red arrows show the dust vectors (i.e.shift in ) assuming 0.2 magnitudes of UV attenuation for theCalzetti et al. (2000) attenuation curve and the SMC extinction curve(Gordon et al. 2003).These models indicate that to achieve an observed UV continuum slope of  obs = −2.6 including dust would require a counterintuitive scenario in which a young burst, with a large (≳ 50 per cent) escape fraction is at the same time experiencing ≳ 0.2 magnitudes of UV attenuation (see text for discussion).

Figure 11 .
Figure11.The fraction of ionising photons (  esc ) as a function of observed UV continuum slope () for our primary wide-area JWST sample.The values of  esc for individual galaxies (grey circles) were calculated using the  −  esc relation at  = 0.3 derived by Chisholm et al. (2022) using (primarily) data from the LzLCS survey(Flury et al. 2022).Open circular data points show galaxies with unphysical values of  esc > 1 based on this calibration (i.e. the galaxies in our sample that have been scattered to extremely blue  values).The square data points in red and blue show the population average value of  esc above and below  = 10.5, respectively, using the best-fitting  values from Fig.5.The red and blue horizontal lines show  esc = 0.05 and  esc = 0.20, respectively.Taking the predictions of the low-redshift calibration at face value implies that a substantial fraction of ionising photons (≃ 20 per cent) are escaping into the IGM at  > 10.5.
which was designed to investigate LyC escape in analogues of young high redshift galaxies.Interestingly, the bluest objects in this sample have  ≃ −2.5.The result of applying the Chisholm et al. (2022)  −  esc relation (their equation 18) to our data is shown in Fig. 11.While the scatter for individual objects is substantial, the population average values yield robust constraints.Using the fitted values of ⟨⟩ from Fig. 5, the Chisholm et al. (2022) relation yields average values of ⟨  esc ⟩ = 0.07 at  < 10.5 to ⟨  esc ⟩ = 0.19 at  > 10.5.

Figure 12 .
Figure12.The emission rate of ionising photons into the IGM per comoving Mpc 3 (  ion / −1 ) for different assumed values of the ionising production efficiency (  ion /Hz erg −1 ) and different integration limits of the UV LF.All calculations assume  esc = 0.2, which we infer using the relation ofChisholm et al. (2022) based on our measured value of ⟨⟩ ≃ −2.6 (see text and Fig.11) and the  = 11 UV LF ofMcLeod et al. (2024).The dashed horizontal lines show the equilibrium value of  ion needed to maintain an ionised IGM fraction of  H ii (as indicated by the labels).The red curve shows  ion as a function of  UV,lim assuming log(  ion /Hz erg −1 ) = 25.8 (i.e. at the upper end of  ion estimates for young metal-poor galaxies; e.g.Tang et al. 2019).The curve turns from solid to dashed below the limiting magnitude of current LF constraints ( UV ≃ −18).At this efficiency, the already observed population of galaxies (i.e. UV ≲ −18; vertical black line) would be supplying enough ionising photons to maintain  H ii > 0.1.The orange and blue curves represent log(  ion /Hz erg −1 ) = 25.5 and log(  ion /Hz erg −1 ) = 25.2, which may be more representative of the typical value during reionisation (e.g.Bouwens et al. 2016;Matthee et al. 2017).At the lowest estimated ionising efficiencies, the observed population would be sufficient to maintain  H ii ≃ 0.05.Taking into account the faint, unseen, galaxy population at  UV > −18 implies relatively high (likely ≳ 5 per cent) ionised fractions at  = 11.

Table 1 .
The best-fitting UV continuum slopes () for our new wide-area JWST sample of galaxy candidates at  > 9.The first column gives the source ID.Columns two and three give the best-fitting photometric redshift ( phot ) and absolute UV magnitude ( UV ) from our fitting procedure described in Section 2.3.Column four gives the derived UV continuum slope .

Table 2 .
Cullen et al. (2023)standard errors derived for our wide-area JWST and combined samples.The first column gives the sample name (as referred to in the text).In the second column we report the inverse-variance weighted mean and standard error of the individual  values.In the third 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.It is worth noting that − as inCullen et al. (2023)− the full sample average values reported in Table2are not more extreme than the bluest galaxies observed in the local Universe (e.g.NGC 1705;  = −2.46, UV = −18;

Table 3 .
Average values and standard errors derived for our wide-area JWST sample and combined sample as a function of redshift.The first column defines the redshift range.In the second column, we report the inversevariance weighted mean and standard error of the individual  values.In the third column, we report the median redshift () and  MAD .In the fourth column, we report the median  UV and  MAD .The values in columns two and three correspond to the yellow square data points shown in Fig.4.Our sample shows a clear trend in ⟨⟩ with redshift.Crucially, the ⟨⟩ values are not expected to be strongly biased due to the  −  UV relation, as the  UV distributions are similar in each redshift bin (with the exception of the 8 <  < 9 bin for the combined sample; see text for discussion).In the highest redshift bin, the extremely blue value of ⟨⟩ ≃ −2.6 is consistent with dust-free stellar populations (see Section 4).

Table 5 .
The best-fitting values for d/d UV and  (  UV = −19) for the fits shown in Fig.7.These fits assume the following functional form: = d/d UV × (  UV + 19) +  (  UV = −19).The first column gives the redshift range.In the second column we report the best-fitting value of d/d UV ; note that this value is fixed for the two bins at  > 10 (see text for details).In the third column we report the best-fitting value of  (  UV = −19).
report Plots of  versus stellar population age for four different SPS models.The stellar populations models considered are starburst99 continuous star formation models including a contribution from stripped stars and Flexible Stellar Population Synthesis (fsps) single-burst models (Conroy et al. 2009; Conroy & Gunn 2010; Byler et al. 2017).Full details of the nebular continuum modelling are given in