Rest-frame UV Colors for Faint Galaxies at z ∼ 9–16 with the JWST NGDEEP Survey

NGDEEP is a deep imaging and slitless spectroscopic JWST Cycle 1 treasury program (PID 2079, PIs: S. Finkelstein, C. Papovich, N. Pirzkal) targeting constraints in feedback processes in galaxies across cosmic time (Bagley et al. 2023). NGDEEP's parallel ultradeep NIRCam imaging in the Hubble Ultra-Deep Field "Parallel-2" (HUDF-Par2) field should ultimately detect ≳100 galaxies at redshifts z = 9–16 over $5\ {\mathrm{arcmin}}^{2}$ at 1–5 μm on the deepest HST F814W imaging in the sky (m ${}_{\mathrm{lim},{\rm{F}}814{\rm{W}}}\sim$ 30) from the HUDF12 campaign (Ellis et al. 2013; Koekemoer et al. 2013). NGDEEP should ultimately reach m ∼ 30.8, up to ∼0.5 mag deeper than Guaranteed Time Observation (GTO) NIRCam surveys in some filters, and will provide robust measurements for β and galaxy morphologies for UV-faint, low-mass galaxies. While the full program was planned for early 2023, only half of the program was obtained due to a temporary suspension of operations for NIRISS, causing the NGDEEP observations to be pushed to the edge of the visibility window. The next visibility window satisfying the PA requirement of the parallel observations will occur in early 2024 when the remaining observations are expected to be taken. In this study, we report results using NIRCam data from the first half of the NGDEEP program, which is among the deepest data yet obtained by JWST.

We use the NIRCam and HST/ACS F814W imaging mosaics and photometric catalogs presented in Leung et al. (2023). We refer the reader there for full details but summarize the reduction and cataloging process here. The NIRCam imaging data are reduced using the JWST pipeline ¹⁴ (Bushouse et al. 2022) with custom modifications. We process the exposures through Stage 1 of the pipeline, applying custom procedures to remove snowballs, wisps, and 1/f noise. The images are then processed through Stage 2. Since only preflight flats are available for the short-wavelength filters through the CRDS, we apply custom sky flats to these filters, and we use the recently updated CRDS flats for the long-wavelength filters. Before combining the exposures into mosaics through Stage 3, a custom version of the TweakReg routine is used to align the NIRCam images with the ACS mosaics (Koekemoer et al. 2011, 2013), whose astrometry is tied to Gaia DR3. ¹⁵ Finally, we perform a custom background subtraction procedure on all the mosaics.

Aperture photometry was performed using Source Extractor (Bertin & Arnouts 1996). Colors were measured using small elliptical Kron apertures, accounting for the varying point-spread function (PSF) across 0.8–5 μm via matching the PSF of bands bluer than F277W to the F277W PSF. For F356W and F444W, source-by-source correction factors were derived to account for missing flux in the F277W-defined aperture. Total fluxes were derived first by scaling to the flux ratio in F277W between the flux in the small elliptical aperture to the default Source Extractor "MAG_AUTO" Kron parameters, with a residual aperture correction (typically 5%–10%) measured and applied following source injection simulations. Flux uncertainties were measured empirically from the data and assigned to each object for each filter based on the size of its elliptical aperture.

We make use of the galaxy sample from Leung et al. (2023), summarizing their sample selection and completeness here. They incorporate a combination of flux detection signal-to-noise ratios (S/Ns) with quantities derived from photometric redshift fitting with EAZY (Brammer et al. 2010), using the full probability density functions of the photometric redshift, P(z), and the best-fit redshift, z_a. Their final sample was then visually inspected to reject any spurious sources and to ensure a real dropout in all filter bands blueward of the Lyman break corresponding to the estimated redshift. We here make use of 36 of the full sample of 38 galaxies from Leung et al. (2023), which span photometric redshifts z ∼ 8.5–15.5 and UV magnitudes −21 ≤ M_UV ≤ −16. We exclude the two galaxies that Leung et al. (2023) found had red colors in the long-wavelength bands. As discussed in Leung et al. (2023), this may indicate that these sources have an AGN component contributing to their emission, similar to other sources identified at brighter magnitudes and lower redshifts (e.g., Kocevski et al. 2023; Labbe et al. 2023; Matthee et al. 2023; Barro et al. 2024).

To explore any potential sample bias against UV-red objects, we make use of the same source injection simulations described in Leung et al. (2023). They measured the completeness of their photometric selection by injecting and attempting to recover 50,000 mock objects with varied F277W magnitudes, colors, and surface brightness profiles across images, measuring completeness via the proportion recovered in both photometric and sample selection criteria within specified magnitude bins. Galaxy size, a potential influence on sample completeness, is incorporated as an estimation parameter. Any bias in size measurements, identified by comparing input and recovered half-light radii, is corrected during the completeness calculation to ensure accurate data evaluation and selection robustness. We explore the completeness as a function of β and find that our sample selection remains sensitive to significantly redder UV slopes than our reddest object, and we observe that the completeness is only ∼15% lower at β = −1.5 than at β = −2.5. We conclude that our photometric selection criteria do not bias the observed β distribution.

For each source in the sample, we run Bagpipes (Carnall et al. 2021), a Bayesian SED fitting code where photometric data is provided as input along with a suite of user-defined priors, and the posterior distributions of galaxy properties and corresponding model spectra are returned. Our priors span a wide parameter space to ensure that we are not omitting valuable information (see Table 1) and account for a variable star formation history (SFH), by fitting a delayed-tau SFH model. We aim to measure β from these model spectra directly, following Finkelstein et al. (2012).

Table 1.  Bagpipes Priors Defined for This Work

Parameter	Range	Description
SFH Component
sfh["age"]	(0.002, 13.0)	Age of the galaxy in Gyr
sfh["tau"]	(0.1, 14.0)	Delayed decay time in Gyr
sfh["metallicity"]	(0.0, 2.5)	Metallicity in units of Z⊙
sfh["massformed"]	(4.0, 13.0)	${\mathrm{log}}_{10}$ of total mass formed in units of M⊙
Dust Component
dust["type"]	Calzetti	Shape of the attenuation curve
dust["Av"]	(0.0, 4.0)	Dust attenuation in units of magnitude
dust["eta"]	1.0	Multiplicative factor on Av for stars in birth clouds
Nebular Component
nebular["logU"]	(−4.0, −1.0)	${\mathrm{log}}_{10}$ of the ionization parameter
nebular["fesc"]	(0.0, 1.0)	Lyman continuum escape fraction
Additional Fit Instructions
fit_instructions["redshift"]	(0.0, 20.0)	Redshift range tested
fit_instructions["redshift_prior"]	CDF(z) function	CDF of the EAZY P(z) as prior on redshift

Download table as:  ASCIITypeset image

We made two modifications to the Bagpipes source code: (1) We ensure that Bagpipes is consistent with the EAZY results used to select the galaxy sample by implementing an additional redshift prior that is the cumulative distribution function (CDF) of our EAZY P(z) while allowing the redshift range tested to be free from z = 0 to 20. To ignore any potential low-redshift solutions with Bagpipes, we set the CDF(z < 5) to zero. (2) We also incorporate a new free parameter added to Bagpipes, the Lyman continuum escape fraction, f_esc, to allow for the flexibility of measuring bluer colors within our fitting (f_esc was otherwise by default set to zero, which can lead to nebular continuum dominating at young ages, and thus β > −3 even for dust-free, metal-poor populations). By incorporating this f_esc component, alongside other tunable priors, we are able to allow Bagpipes to generate a model reaching a minimum blue "floor" of β_SED = −3.25 when f_esc = 1 (compared to β_SED = −2.44 when f_esc = 0). Ideally, this β floor would be bluer than the bluest observed value (which, as we have shown, is not possible if f_esc is fixed to zero). This should result in the recovery of the correct median β without any systematic bias. Although we incorporate this new parameter, we do not interpret the f_esc posteriors below, as they are quite broad based on the photometric constraints.

Figure 1 shows the prior on β_SED with the escape fraction parameter now included in the parameter space, derived by fitting a mock source whose S/N is zero in all bands, allowing the redshift to be free from z = 0 to 20. Our β_SED prior spans a range of −3.25 < β_SED < 10.0. We also test the ability of Bagpipes to accurately recover a mock spectrum whose "true" UV slope is −3.0 with an S/N of ∼5 and ∼10 in all bands. With the same priors, we show that Bagpipes can recover β_SED values close to the "true" UV slope value (there is a slight bias toward redder values, which we further discuss in Section 4.2). For this work, we make use of the default BC03 stellar models (Bruzual & Charlot 2003) and Calzetti et al. (1994) dust models Bagpipes provides, though we note that future work using Population III star models may be able to reach bluer values.

Zoom In Zoom Out Reset image size

We utilize the median (and the difference between the median and the 68% confidence bounds) for redshift, stellar mass, UV magnitude, dust attenuation, SFR, and mass-weighted age as our final values (and error bars) in this work. For each source, we set Bagpipes to fit model SEDs to the photometric data points and their errors, return the best-fit SED model, and return 1000 draws from the resulting posteriors for various galaxy properties as a function of these SEDs.

We measured the UV slope, β_SED, from each of the 1000 posterior model spectra via a power-law fit to all model flux density data points from λ_rest = 1500 to 3000 Å, while excluding a 10 Å region around wavelengths where we expect high-redshift emission lines to lie (λ_rest = 1549, 1909, 2796, 2803 Å), resulting in an average of 131 data points per model SED for the UV slope calculation. While Calzetti's work primarily focused on lower-redshift sources that often contained significant emission lines, and thus the need for wavelength windows, our study concentrates on galaxies at much higher redshifts, where the prevalence of emission lines is notably reduced. This distinction is crucial since the presence of strong emission lines can significantly affect the interpretation of UV spectral slopes. Therefore, our approach prioritizes capturing the bluer colors, particularly in cases of dust-free, metal-poor populations, which are more common at high redshifts. We note that we explicitly exclude the bluest windows from the Calzetti et al. (1994) definition of β to exclude any potential reddening due to the observed strong Lyα damping wings at high redshift (e.g., Arrabal Haro et al. 2023; Umeda et al. 2023; Heintz et al. 2023). We adopt the median and the difference from the 1σ bounds from these measurements as our final β_SED value and uncertainties, where $\mathrm{log}({f}_{\lambda })=\beta \mathrm{log}(\lambda )+\mathrm{log}({y}_{\mathrm{int}})$. This approach, while deviating from the wavelength windows of Calzetti et al. (1994), similarly avoids strong emission lines. For our SED figures (see Section 4.1 and Appendix Figure A1), we display the lowest χ² model SED, the central 68%, and the full range covered by the full model posterior.

Here we describe our process of measuring the UV spectral slope via photometric power-law fitting to the observed photometry. This is a more commonly used method to measure the UV spectral slope that is not reliant on stellar population models, though it is very sensitive to the number of photometric data points available. When measuring the UV spectral slope with photometric power-law fitting to the observed photometry, β_PL, we redshift the rest-frame λ = 1500–3000 Å regime to the corresponding median redshift estimated from EAZY, z_a, for each galaxy in the sample and keep the filters whose central wavelength within each filter curve falls fully within the redshifted wavelength range. The redshift range of our sample, z ∼ 9–16, results in typically two to three photometric data points being available to fit a line (namely, F200W, F277W, F356W, and, at the highest redshifts, F444W). Once the photometric data points within the wavelength range are determined, we fit the data points to a line (as defined at the end of Section 3.2) and run this process through Emcee (Foreman-Mackey et al. 2013). This procedure maximizes the likelihood that the model described by three free parameters matches the observed photometry for a given source: β, y_int, and fractional error factor log(f), where it is assumed that the likelihood function and its uncertainties are simply a Gaussian where the variance is underestimated by some fractional amount. The assumed priors on all three free parameters span a wide parameter space, allowing for walkers to converge, and are defined as follows: −10 < β < 2, −60 < y_int < 30, and −60 < log(f) < 1. Results are derived from the median and 68th percentile of the posterior distribution on these three parameters from a chain consisting of 5000 steps and a burn-in of 100 steps (i.e., the Markov Chain Monte Carlo process is iterated over 5000 times to estimate these parameters).

UV spectral slopes for our sample from both measurement methods are shown in Figure 2. This plot shows both the original and corrected β values (corrections applied are discussed in Section 4.2). We show that the SED fitting method yields smaller uncertainties on average (see Figure 3 for examples of this)—this is simply due to the use of more data points along the specified wavelength range mentioned in Section 3.2, given that the SED is a good fit to the data (future works with varying stellar models may help alleviate any limitations our current model may have; see Section 5). After the correction is applied, we find that β_SED is in fairly good agreement with β_PL for some objects, while the β_PL measurements show significantly larger scatter, extending to both much bluer and much redder values than β_SED. Specifically, while no objects reach β_SED < −3, there are many whose β_PL < −3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In comparing our work to other observations, we wanted to ensure that we used data sets with selection strategies similar to the ones we implemented and that the redshift range overlapped with this work. The observational data sets we used to compare values of β versus redshift, stellar mass, and UV magnitude are as follows: UV slope values and corresponding galaxy parameters at z = 9–11 by Tacchella et al. (2022), whose sample consists of 11 bright (H-band magnitude < 26.6) galaxy candidates selected in the HST CANDELS fields from Finkelstein et al. (2022) with photometry spanning the UV/optical with HST and near-infrared with Spitzer/IRAC. The work of Topping et al. (2022) utilizes the initial NIRCam data release at z = 7–11 for 123 galaxies from the JWST Cosmic Evolution Early Release Science (CEERS; Finkelstein et al. 2017) survey; this survey overlaps with that of the HST EGS field. At z = 8–16, we use the 61 galaxy observations in Cullen et al. (2023b) from three JWST fields, SMACS J0723, GLASS, and CEERS (Pontoppidan et al. 2022), and the ground-based COSMOS/UltraVISTA survey (McCracken et al. 2012). Finally, we include a sample of 18 galaxies in the NGDEEP field at z = 8–16 from the works of Austin et al. (2023).

In Figure 3, we show six galaxies from our sample of 36 galaxies in the NGDEEP field with their corresponding photometry and model SEDs. The figure shows SEDs for the highest-redshift galaxy in our sample at ${z}_{\mathrm{SED}}={15.5}_{-0.3}^{+0.3}$, whose photometric power-law and SED-fitted UV spectral slopes are both very blue, ${\beta }_{\mathrm{SED}}=-{2.4}_{-0.1}^{+0.1}$ (−2.3 before correction) and ${\beta }_{\mathrm{PL}}=-{2.4}_{-0.2}^{+0.2}$. This is consistent with a low-metallicity, young stellar population, though spectroscopic follow-up is needed to confirm both the redshift and the metallicity. We also obtained photometric data for a $z={8.5}_{-0.4}^{+0.3}$ galaxy whose UV magnitude is quite faint, ${M}_{\mathrm{UV}}=-{16.2}_{-0.2}^{+0.2}$, and still exhibits very blue colors with both photometric power-law and SED fitting, ${\beta }_{\mathrm{SED}}=-{2.3}_{-0.2}^{+0.2}$ (−2.2 before correction) and ${\beta }_{\mathrm{PL}}=-{2.7}_{-0.5}^{+0.5}$, highlighting the large scatter in β_PL at faint measured fluxes.

We also show SEDs for one of the bluest sources in our sample at $z={8.6}_{-0.1}^{+0.1}$ with a ${\beta }_{\mathrm{SED}}=-{2.7}_{-0.1}^{+0.1}$ (−2.5 before correction) and one of the "reddest" sources in our sample with a ${\beta }_{\mathrm{SED}}=-{1.5}_{-0.2}^{+0.3}$ (−1.3 before correction). Although when modeled with Bagpipes the bluest galaxy is considered to be a fairly young galaxy, ∼10⁷ yr old (similar to the bulk of our sample), the lack of dust attenuation and relatively high specific SFR both play a role in the extremely blue slope. The redder galaxy, on the other hand, is slightly fainter, with a bit more dust attenuation contributing to our measurement. Lastly, we show two examples of galaxies whose β_SED fits are either bluer or redder than the measured photometric power-law fits, β_PL. However, in both instances, we show that our SED fitting method yields smaller errors than photometric power-law fitting overall. As discussed, depending on the filter wavelength coverage and the redshift range, the power-law method can be limited to just two data points, yielding not only larger uncertainties but also larger scatter in the measured value.

To test the accuracy and assess any bias in our SED and photometric power-law measurements of the UV spectral slope, we utilize source injection simulations. Full details of this process are described in Leung et al. (2023), but we summarize briefly here. The simulated sources are created with a range of F277W magnitudes, colors, and surface brightness profiles. Once added to a given image, photometry and resulting photometric redshifts are defined in the same manner as what is done on real images. Utilizing the same selection criteria as the observed galaxies, we randomly draw a sample of 1000 simulated galaxies that are defined to have "true" input UV spectral slopes, stellar mass, and UV magnitudes. We ensure that these galaxies span a wide enough parameter space for each property to fully align with the dynamic range of the observations. With this sample, we run Bagpipes and Emcee in the same fashion as the observed sample in order to obtain β, UV magnitude, and stellar mass. We then explore whether there exist any measurement biases by taking the difference between β_SED, β_PL, and the input β as functions of input β. As shown in the left and middle panels of Figure 4, we can see that SED fitting yields a closer fit to the "true" UV slope than with photometric power-law fitting, where anything past Δβ ± 0.5 is to be considered a strong outlier (for SED fitting, this fraction is 75 out of 1000 galaxies, and for photometric power-law fitting, this fraction is 274 out of 1000 galaxies). In particular, the recovered photometric power-law β values preferentially scatter toward bluer slopes.

Zoom In Zoom Out Reset image size

Although our SED fitting method yields tighter agreement between the input and recovered UV spectral slope, we still hit a "minimum β" plateau in what is returned from Bagpipes, as evidenced by the good agreement at β > −2, which begins to flatten at bluer input colors. Photometric power-law fitting, with its larger spread, essentially shows that more filters are needed past 2–3 at these high redshifts in order to tighten measurements. This, however, would require a series of medium-band filters in the rest-frame UV regime that could probe just as deep as the current filters being used, which would be very costly observationally. Alternatively, other approaches to photometric power-law fitting (such as the recent works by Cullen et al. 2023b, 2023a) fit through the Lyman break, wherein redshift is jointly modeled, and this can yield at least one extra filter. We leave a direct comparison between this and SED fitting to future works.

Knowing that SED fitting and photometric power-law fitting have limitations in their ability to precisely measure the UV spectral slope, we use these simulations to derive a bias correction factor. As we cannot know the true value of β for real sources, we derive this correction as a function of observed rest-UV color. For each source, both simulated and observed, we measure a rest-UV color by taking the difference, in magnitude, of the first two filters redward of the filter where the Lyman break is expected to lie.

To adjust for discrepancies in measured versus true UV slopes, we calculate a correction factor, Δβ, which is the difference between measured UV slopes (from SED or photometric fitting) and known input values. We calculate this quantity in UV color bins of 0.25, using only those with over 20 data points. To apply this correction to a given real object, we interpolate the Δβ curves to a given observed UV color value, assigning edge values to outliers. As illustrated in Figure 2, applying this correction even just to β_SED aligns original and adjusted β_SED values closer to β_PL across various redshifts (detailed in Table 2). The Pearson R correlation coefficient analysis shows an improvement from ρ = 0.38 to ρ = 0.48 post-correction, indicating alleviation of bias in the SED-derived UV slopes, evidenced by a shift toward bluer β_SED within the error margins of the sources. These corrected β_SED values will be used as the baseline for all subsequent analyses within this study.

Table 2. Redshift and UV Spectral Slopes for This Sample

ID	R.A.	Decl.	z_a a	zSED b	βPL c	No. Filters	βSED	Adjusted
						for βPL		βSED d
250	53.249451	−27.883313	${11.62}_{-0.33}^{+0.21}$	${11.38}_{-0.21}^{+0.18}$	$-{2.85}_{-0.26}^{+0.27}$	3	$-{2.29}_{-0.13}^{+0.15}$	−2.44
1191	53.266583	−27.876581	${12.31}_{-0.30}^{+1.23}$	${12.10}_{-0.21}^{+0.21}$	$-{2.39}_{-0.47}^{+0.47}$	2	$-{2.11}_{-0.16}^{+0.16}$	−2.26
1369	53.249467	−27.875710	${15.82}_{-0.81}^{+0.12}$	${15.46}_{-0.27}^{+0.30}$	$-{2.42}_{-0.19}^{+0.18}$	3	$-{2.29}_{-0.14}^{+0.14}$	−2.44
1716	53.251057	−27.796992	${11.29}_{-1.47}^{+0.30}$	${10.27}_{-0.75}^{+0.69}$	$-{3.70}_{-0.68}^{+0.69}$	3	$-{1.91}_{-0.27}^{+0.31}$	−2.17
2067	53.239797	−27.800244	${10.75}_{-0.66}^{+0.42}$	${10.57}_{-0.36}^{+0.36}$	$-{1.39}_{-0.91}^{+0.89}$	2	$-{2.08}_{-0.20}^{+0.20}$	−2.20
2470	53.248546	−27.802743	${8.56}_{-0.63}^{+0.15}$	${8.17}_{-0.27}^{+0.24}$	$-{2.73}_{-0.27}^{+0.27}$	3	$-{2.40}_{-0.14}^{+0.19}$	−2.64
2497	53.248775	−27.803090	${9.19}_{-0.21}^{+0.60}$	${9.16}_{-0.18}^{+0.18}$	$-{2.27}_{-0.63}^{+0.62}$	2	$-{2.39}_{-0.15}^{+0.14}$	−2.56
3514	53.256875	−27.807957	${11.23}_{-0.90}^{+0.27}$	${10.72}_{-0.60}^{+0.45}$	$-{3.66}_{-0.44}^{+0.44}$	3	$-{2.28}_{-0.22}^{+0.19}$	−2.53
4134	53.245546	−27.814372	${10.69}_{-0.30}^{+0.18}$	${10.57}_{-0.15}^{+0.15}$	$-{2.53}_{-0.22}^{+0.22}$	2	$-{2.37}_{-0.13}^{+0.11}$	−2.50
4330	53.264834	−27.816024	${9.16}_{-0.18}^{+0.96}$	${9.13}_{-0.15}^{+0.18}$	$-{2.79}_{-0.84}^{+0.82}$	2	$-{2.27}_{-0.15}^{+0.15}$	−2.37
4674	53.245601	−27.817588	${8.62}_{-0.21}^{+0.18}$	${8.50}_{-0.12}^{+0.15}$	$-{2.76}_{-0.21}^{+0.21}$	3	$-{2.49}_{-0.26}^{+0.14}$	−2.67
4740	53.257473	−27.817828	${9.19}_{-0.21}^{+0.72}$	${9.19}_{-0.18}^{+0.24}$	$-{2.80}_{-0.79}^{+0.78}$	2	$-{2.47}_{-0.21}^{+0.16}$	−2.66
4919	53.262137	−27.818774	${9.34}_{-0.24}^{+1.02}$	${9.58}_{-0.39}^{+0.45}$	$-{1.82}_{-0.86}^{+0.85}$	2	$-{2.28}_{-0.15}^{+0.16}$	−2.46
5118	53.261009	−27.819895	${9.25}_{-0.06}^{+0.09}$	${9.22}_{-0.03}^{+0.03}$	$-{2.66}_{-0.11}^{+0.11}$	2	$-{2.30}_{-0.04}^{+0.04}$	−2.37
5947	53.248954	−27.822988	${8.89}_{-0.57}^{+1.50}$	${8.80}_{-0.36}^{+0.33}$	$-{2.04}_{-0.46}^{+0.46}$	3	$-{1.83}_{-0.23}^{+0.24}$	−1.94
6134	53.248902	−27.823695	${9.07}_{-0.12}^{+0.09}$	${9.01}_{-0.09}^{+0.06}$	$-{2.32}_{-0.26}^{+0.26}$	2	$-{2.36}_{-0.10}^{+0.09}$	−2.48
6477	53.250469	−27.825099	${8.74}_{-0.30}^{+0.33}$	${8.71}_{-0.21}^{+0.24}$	$-{2.46}_{-0.31}^{+0.31}$	3	$-{2.26}_{-0.16}^{+0.14}$	−2.33
6952	53.239303	−27.827168	${8.68}_{-0.75}^{+0.45}$	${8.53}_{-0.42}^{+0.36}$	$-{2.66}_{-0.49}^{+0.50}$	3	$-{2.16}_{-0.20}^{+0.24}$	−2.28
6980	53.256828	−27.827244	${10.06}_{-0.60}^{+0.45}$	${9.76}_{-0.30}^{+0.33}$	$-{3.33}_{-0.73}^{+0.73}$	2	$-{2.42}_{-0.24}^{+0.19}$	−2.64
7530	53.237017	−27.829892	${8.71}_{-0.45}^{+0.57}$	${8.71}_{-0.30}^{+0.30}$	$-{2.09}_{-0.35}^{+0.36}$	3	$-{1.82}_{-0.24}^{+0.27}$	−2.00
7722	53.242580	−27.830637	${10.93}_{-0.39}^{+0.30}$	${10.66}_{-0.36}^{+0.30}$	$-{2.84}_{-0.46}^{+0.45}$	3	$-{2.53}_{-0.25}^{+0.20}$	−2.69
8042	53.237915	−27.832320	${8.53}_{-0.99}^{+0.18}$	${8.17}_{-0.51}^{+0.36}$	$-{2.92}_{-0.47}^{+0.47}$	3	$-{2.19}_{-0.21}^{+0.20}$	−2.45
8165	53.234410	−27.833172	${10.48}_{-0.96}^{+0.36}$	${9.82}_{-0.48}^{+0.48}$	$-{5.15}_{-0.94}^{+0.95}$	2	$-{2.42}_{-0.27}^{+0.20}$	−2.67
8427	53.235623	−27.834282	${8.98}_{-0.30}^{+0.18}$	${8.89}_{-0.18}^{+0.12}$	$-{2.99}_{-0.28}^{+0.28}$	3	$-{2.38}_{-0.15}^{+0.16}$	−2.61
8461	53.235797	−27.834499	${10.36}_{-1.26}^{+0.27}$	${9.53}_{-0.43}^{+0.64}$	$-{2.88}_{-0.98}^{+0.98}$	2	$-{2.32}_{-0.23}^{+0.17}$	−2.49
8894	53.241008	−27.828822	${9.34}_{-0.06}^{+0.87}$	${9.43}_{-0.18}^{+0.39}$	$-{2.64}_{-0.54}^{+0.55}$	2	$-{2.40}_{-0.14}^{+0.12}$	−2.51
9261	53.270902	−27.841204	${8.74}_{-0.24}^{+0.15}$	${8.62}_{-0.15}^{+0.15}$	$-{3.27}_{-0.33}^{+0.34}$	3	$-{2.44}_{-0.11}^{+0.13}$	−2.55
9555	53.258685	−27.847324	${10.45}_{-0.84}^{+0.39}$	${10.00}_{-0.51}^{+0.45}$	$-{2.65}_{-0.96}^{+0.95}$	2	$-{2.58}_{-0.19}^{+0.17}$	−2.72
10296	53.276910	−27.850568	${10.48}_{-0.30}^{+0.33}$	${9.46}_{-0.03}^{+0.09}$	$-{3.08}_{-0.37}^{+0.37}$	2	$-{2.33}_{-0.11}^{+0.12}$	−2.52
11522	53.242062	−27.855079	${10.84}_{-0.30}^{+0.36}$	${10.75}_{-0.21}^{+0.18}$	$-{2.39}_{-0.37}^{+0.37}$	2	$-{2.15}_{-0.14}^{+0.14}$	−2.26
12453	53.280001	−27.858265	${9.88}_{-0.54}^{+0.78}$	${9.82}_{-0.42}^{+0.42}$	$-{1.10}_{-1.81}^{+1.65}$	2	$-{2.13}_{-0.18}^{+0.24}$	−2.27
13290	53.258417	−27.861651	${8.95}_{-0.18}^{+1.59}$	${9.10}_{-0.30}^{+0.63}$	$-{1.58}_{-0.38}^{+0.38}$	3	$-{1.30}_{-0.24}^{+0.27}$	−1.48
13406	53.240619	−27.862122	${9.31}_{-0.15}^{+1.05}$	${9.52}_{-0.33}^{+0.51}$	$-{2.16}_{-0.81}^{+0.80}$	2	$-{2.32}_{-0.16}^{+0.15}$	−2.46
13782	53.244589	−27.863587	${9.22}_{-0.18}^{+1.53}$	${9.13}_{-0.30}^{+0.51}$	$-{3.26}_{-0.78}^{+0.78}$	2	$-{1.88}_{-0.22}^{+0.23}$	−1.98
17672	53.249058	−27.815575	${9.40}_{-0.09}^{+0.03}$	${9.32}_{-0.04}^{+0.04}$	$-{2.45}_{-0.10}^{+0.10}$	2	$-{2.14}_{-0.05}^{+0.06}$	−2.23
17674	53.233517	−27.816677	${8.56}_{-0.15}^{+0.15}$	${8.56}_{-0.09}^{+0.09}$	$-{2.62}_{-0.10}^{+0.10}$	3	$-{2.54}_{-0.06}^{+0.06}$	−2.67

Notes.

^a Redshift determined from EAZY P(z). ^b Redshift determined from Bagpipes with the EAZY redshift prior applied. ^c Photometric power-law fit of the UV spectral slope as determined with Emcee. ^d Adjusted β_SED values with correction applied as discussed in Section 4.2.

Download table as:  ASCIITypeset image

We note that this process was similarly done with photometric power-law fitting. However, we obtain an even stronger bias, given the small number of data points being fit. As such, we do not apply these corrections to the photometric power-law fitted values, though we do show the derived correction values for both methods in the right panel of Figure 4.

We compare the results of the Bagpipes posteriors for galaxy parameters versus the bias-corrected values of β_SED and discuss any correlations. Previous works have explored correlations mainly between β and stellar mass (Finkelstein et al. 2010, 2012), UV absolute magnitude (Bouwens et al. 2010, 2012), and dust attenuation (Calzetti et al. 1994; Meurer et al. 1999). Here we explore monotonic trends with these parameters and with SFR and mass-weighted age using Spearman R correlation coefficients, ¹⁶ ρ (Spearman 1904). We perform a Monte Carlo resampling to estimate the Spearman correlation coefficient and the corresponding p-value between the corrected β_SED and other galaxy parameters, accounting for their uncertainties. We repeatedly add normally distributed random noise to the data and recalculate the ρ and p-value. The median ρ and p-value, along with their 68% and 95% confidence intervals, are noted in the legends of Figures 5, 6, and 7. We offer the median and difference from the 68th confidence bounds for each parameter in Table 3. For Figures 5 and 6, we plot our corrected β_SED values versus redshift, UV magnitude, and stellar mass and offer comparisons with other works that have been done at high redshift (Tacchella et al. 2022; Topping et al. 2022; Cullen et al. 2023b; Austin et al. 2023).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 3.  Bagpipes Median and 1σ Posteriors for Galaxy Parameters

ID	zSED	MUV	$\mathrm{log}({M}_{* }/{M}_{\odot })$	Av (mag)	log(Age yr−1)	log(SFR/M⊙ yr−1)
250	${11.38}_{-0.21}^{+0.18}$	$-{18.18}_{-0.09}^{+0.1}$	${7.93}_{-0.33}^{+0.16}$	${0.12}_{-0.81}^{+0.14}$	${7.86}_{-0.42}^{+0.20}$	$-{0.16}_{-0.14}^{+0.12}$
1191	${12.10}_{-0.21}^{+0.21}$	$-{18.63}_{-0.09}^{+0.11}$	${8.17}_{-0.33}^{+0.11}$	${0.26}_{-0.14}^{+0.20}$	${7.83}_{-0.46}^{+0.20}$	${0.13}_{-0.18}^{+0.17}$
1369	${15.46}_{-0.27}^{+0.30}$	$-{18.98}_{-0.08}^{+0.08}$	${8.09}_{-0.31}^{+0.19}$	${0.15}_{-0.11}^{+0.15}$	${7.65}_{-0.48}^{+0.22}$	${0.13}_{-0.30}^{+0.15}$
1716	${10.27}_{-0.75}^{+0.69}$	$-{16.85}_{-0.26}^{+0.29}$	${7.72}_{-0.26}^{+0.21}$	${0.40}_{-0.23}^{+0.30}$	${7.94}_{-0.38}^{+0.18}$	$-{0.40}_{-0.19}^{+0.20}$
2067	${10.57}_{-0.36}^{+0.36}$	$-{18.17}_{-0.16}^{+0.18}$	${8.10}_{-0.23}^{+0.16}$	${0.29}_{-0.17}^{+0.20}$	${7.87}_{-0.38}^{+0.21}$	${0.02}_{-0.13}^{+0.13}$
2470	${8.17}_{-0.27}^{+0.24}$	$-{17.63}_{-0.13}^{+0.16}$	${7.06}_{-0.13}^{+0.20}$	${0.27}_{-0.15}^{+0.19}$	${6.16}_{-0.08}^{+0.19}$	$-{0.94}_{-0.14}^{+0.21}$
2497	${9.16}_{-0.18}^{+0.18}$	$-{18.24}_{-0.11}^{+0.11}$	${7.60}_{-0.22}^{+0.24}$	${0.09}_{-0.06}^{+0.12}$	${7.44}_{-0.44}^{+0.34}$	$-{0.30}_{-0.27}^{+0.13}$
3514	${10.72}_{-0.60}^{+0.45}$	$-{18.90}_{-0.18}^{+0.21}$	${8.17}_{-0.33}^{+0.20}$	${0.11}_{-0.08}^{+0.15}$	${7.79}_{-0.55}^{+0.27}$	${0.11}_{-0.23}^{+0.13}$
4134	${10.57}_{-0.15}^{+0.15}$	$-{19.32}_{-0.09}^{+0.08}$	${8.18}_{-0.22}^{+0.15}$	${0.09}_{-0.07}^{+0.12}$	${7.64}_{-0.42}^{+0.24}$	${0.22}_{-0.22}^{+0.06}$
4330	${9.13}_{-0.15}^{+0.18}$	$-{17.90}_{-0.12}^{+0.12}$	${7.85}_{-0.22}^{+0.14}$	${0.11}_{-0.08}^{+0.14}$	${7.91}_{-0.34}^{+0.22}$	$-{0.25}_{-0.07}^{+0.09}$
4674	${8.50}_{-0.12}^{+0.15}$	$-{18.17}_{-0.13}^{+0.11}$	${7.26}_{-0.18}^{+0.20}$	${0.06}_{-0.04}^{+0.08}$	${7.00}_{-0.42}^{+0.34}$	$-{0.69}_{-0.21}^{+0.22}$
4740	${9.19}_{-0.18}^{+0.24}$	$-{18.12}_{-0.13}^{+0.15}$	${7.28}_{-0.22}^{+0.29}$	${0.07}_{-0.05}^{+0.09}$	${7.05}_{-0.56}^{+0.49}$	$-{0.66}_{-0.26}^{+0.32}$
4919	${9.58}_{-0.39}^{+0.45}$	$-{17.47}_{-0.14}^{+0.15}$	${7.58}_{-0.29}^{+0.20}$	${0.12}_{-0.09}^{+0.15}$	${7.76}_{-0.47}^{+0.29}$	$-{0.44}_{-0.17}^{+0.11}$
5118	${9.22}_{-0.03}^{+0.03}$	$-{20.14}_{-0.03}^{+0.03}$	${8.71}_{-0.08}^{+0.05}$	${0.04}_{-0.03}^{+0.06}$	${7.96}_{-0.15}^{+0.08}$	${0.60}_{-0.02}^{+0.03}$
5947	${8.80}_{-0.36}^{+0.33}$	$-{17.44}_{-0.23}^{+0.22}$	${8.10}_{-0.16}^{+0.16}$	${0.46}_{-0.20}^{+0.24}$	${8.05}_{-0.35}^{+0.16}$	$-{0.05}_{-0.15}^{+0.17}$
6134	${9.01}_{-0.09}^{+0.06}$	$-{18.50}_{-0.06}^{+0.06}$	${7.60}_{-0.20}^{+0.14}$	${0.06}_{-0.04}^{+0.08}$	${7.33}_{-0.37}^{+0.25}$	$-{0.33}_{-0.22}^{+0.17}$
6477	${8.71}_{-0.21}^{+0.24}$	$-{17.74}_{-0.13}^{+0.15}$	${7.70}_{-0.23}^{+0.19}$	${0.10}_{-0.07}^{+0.13}$	${7.82}_{-0.39}^{+0.26}$	$-{0.35}_{-0.11}^{+0.08}$
6952	${8.53}_{-0.42}^{+0.36}$	$-{16.21}_{-0.18}^{+0.25}$	${7.30}_{-0.24}^{+0.16}$	${0.18}_{-0.13}^{+0.21}$	${7.95}_{-0.38}^{+0.25}$	$-{0.85}_{-0.12}^{+0.12}$
6980	${9.76}_{-0.30}^{+0.33}$	$-{17.42}_{-0.14}^{+0.16}$	${7.21}_{-0.29}^{+0.30}$	${0.08}_{-0.06}^{+0.10}$	${7.37}_{-0.56}^{+0.42}$	$-{0.70}_{-0.34}^{+0.21}$
7530	${8.71}_{-0.30}^{+0.30}$	$-{17.36}_{-0.17}^{+0.22}$	${8.15}_{-0.19}^{+0.14}$	${0.50}_{-0.18}^{+0.21}$	${8.11}_{-0.29}^{+0.12}$	$-{0.05}_{-0.15}^{+0.15}$
7722	${10.66}_{-0.36}^{+0.30}$	$-{18.21}_{-0.16}^{+0.15}$	${7.15}_{-0.16}^{+0.29}$	${0.06}_{-0.04}^{+0.09}$	${6.63}_{-0.47}^{+0.65}$	$-{0.84}_{-0.17}^{+0.33}$
8042	${8.17}_{-0.51}^{+0.36}$	$-{17.03}_{-0.20}^{+0.24}$	${7.56}_{-0.25}^{+0.17}$	${0.14}_{-0.10}^{+0.18}$	${7.96}_{-0.43}^{+0.24}$	$-{0.58}_{-0.12}^{+0.12}$
8165	${9.82}_{-0.48}^{+0.48}$	$-{18.03}_{-0.19}^{+0.19}$	${7.22}_{-0.23}^{+0.38}$	${0.07}_{-0.05}^{+0.10}$	${6.99}_{-0.69}^{+0.60}$	$-{0.72}_{-0.27}^{+0.39}$
8427	${8.89}_{-0.18}^{+0.12}$	$-{17.41}_{-0.13}^{+0.12}$	${7.75}_{-0.19}^{+0.12}$	${0.10}_{-0.07}^{+0.12}$	${8.07}_{-0.27}^{+0.14}$	$-{0.46}_{-0.07}^{+0.09}$
8461	${9.53}_{-0.43}^{+0.64}$	$-{17.08}_{-0.20}^{+0.21}$	${6.91}_{-0.25}^{+0.33}$	${0.08}_{-0.06}^{+0.11}$	${7.13}_{-0.72}^{+0.50}$	$-{1.02}_{-0.28}^{+0.33}$
8894	${9.43}_{-0.18}^{+0.39}$	$-{17.82}_{-0.10}^{+0.10}$	${7.28}_{-0.23}^{+0.29}$	${0.08}_{-0.06}^{+0.10}$	${7.20}_{-0.48}^{+0.39}$	$-{0.65}_{-0.28}^{+0.26}$
9261	${8.62}_{-0.15}^{+0.15}$	$-{18.33}_{-0.11}^{+0.12}$	${8.11}_{-0.19}^{+0.12}$	${0.08}_{-0.06}^{+0.10}$	${8.12}_{-0.27}^{+0.12}$	$-{0.12}_{-0.06}^{+0.10}$
9555	${10.00}_{-0.51}^{+0.45}$	$-{17.02}_{-0.19}^{+0.19}$	${6.82}_{-0.22}^{+0.35}$	${0.08}_{-0.06}^{+0.11}$	${6.98}_{-0.69}^{+0.55}$	$-{1.14}_{-0.24}^{+0.37}$
10296	${9.46}_{-0.03}^{+0.09}$	$-{19.06}_{-0.07}^{+0.07}$	${7.42}_{-0.11}^{+0.15}$	${0.16}_{-0.10}^{+0.12}$	${6.25}_{-0.15}^{+0.35}$	$-{0.57}_{-0.10}^{+0.14}$
11522	${10.75}_{-0.21}^{+0.18}$	$-{18.54}_{-0.10}^{+0.11}$	${8.20}_{-0.21}^{+0.13}$	${0.19}_{-0.11}^{+0.14}$	${7.95}_{-0.34}^{+0.15}$	$-{0.10}_{-0.10}^{+0.09}$
12453	${9.82}_{-0.42}^{+0.42}$	$-{18.07}_{-0.19}^{+0.23}$	${7.95}_{-0.29}^{+0.21}$	${0.21}_{-0.15}^{+0.21}$	${7.86}_{-0.44}^{+0.25}$	$-{0.13}_{-0.16}^{+0.15}$
13290	${9.10}_{-0.30}^{+0.63}$	$-{17.69}_{-0.21}^{+0.25}$	${8.85}_{-0.21}^{+0.18}$	${1.03}_{-0.21}^{+0.25}$	${8.11}_{-0.19}^{+0.11}$	${0.63}_{-0.19}^{+0.20}$
13406	${9.52}_{-0.33}^{+0.51}$	$-{17.37}_{-0.14}^{+0.14}$	${7.37}_{-0.28}^{+0.25}$	${0.12}_{-0.09}^{+0.14}$	${7.56}_{-0.47}^{+0.39}$	$-{0.54}_{-0.28}^{+0.12}$
13782	${9.13}_{-0.30}^{+0.51}$	$-{17.50}_{-0.19}^{+0.22}$	${8.20}_{-0.16}^{+0.14}$	${0.50}_{-0.19}^{+0.22}$	${8.11}_{-0.24}^{+0.11}$	${0.01}_{-0.17}^{+0.16}$
17672	${9.32}_{-0.04}^{+0.04}$	$-{20.33}_{-0.03}^{+0.03}$	${8.92}_{-0.07}^{+0.05}$	${0.17}_{-0.06}^{+0.07}$	${7.98}_{-0.15}^{+0.11}$	${0.80}_{-0.04}^{+0.05}$
17674	${8.56}_{-0.09}^{+0.09}$	$-{19.01}_{-0.05}^{+0.05}$	${8.21}_{-0.35}^{+0.21}$	${0.05}_{-0.03}^{+0.07}$	${7.95}_{-0.48}^{+0.26}$	${0.09}_{-0.14}^{+0.08}$

Download table as:  ASCIITypeset image

Here we examine the strength of the monotonic correlations between the UV spectral slope, β, and various galactic parameters as presented in Figures 5–7, ranked from the strongest to the weakest correlations. The relationship between β and stellar mass exhibits the most pronounced positive monotonic correlation, with a median Spearman correlation coefficient, ρ = 0.47, and a 1σ confidence interval ranging from 0.36 to 0.56. These values suggest a moderate-to-weak correlation strength. Similarly, SFR and mass-weighted age both demonstrate moderate-to-weak positive correlations with β, with median ρ values of 0.43 (1σ CI = [0.32 − 0.53]) and 0.41 (1σ CI = [0.28 − 0.53]), respectively; this suggests that older and more actively star-forming galaxies tend to have redder UV colors. Dust attenuation, A_v , also shows a moderate-to-weak relationship with the UV spectral slope (median ρ of 0.41 (1σ CI = [0.28 − 0.54])), where larger amounts of dust attenuation along the line of sight lead to redder UV slopes (and vice versa; Calzetti et al. 1994; Meurer et al. 1999). Finally, redshift and UV magnitude both show negligible correlations with the UV spectral slope, where each median ∣ρ∣ < 0.12 (1σ CI = [ − 0.08 − 0.16] and [−0.02 − 0.20], respectively). While combining diverse data sets can potentially reveal more significant trends in β, variations in sample selection and analytical methods introduce potential systematic uncertainties that are challenging to quantify. We also abstain from calculating a combined Spearman correlation coefficient across all studies due to the absence of uncertainties in some data, which restricts the application of Monte Carlo resampling techniques defined above.

Zoom In Zoom Out Reset image size

When utilizing Bagpipes to build our model SEDs, we set several priors that allow a large range of star-forming galaxy SED models to be fit. While testing to see how blue our models could go with the set list of priors, stellar grids, and dust laws, we were able to reach a blue floor of β_SED = −3.25. Both our sample of 36 galaxies and the simulated 1000 galaxies were shown to reach an artificial plateau at β_SED ∼ −2.8 regardless of this floor (see Figures 2 and 4). Though we attempt to correct for this bias, future work can improve upon our methodology by adopting bluer stellar grid models, such as the Yggdrasil stellar models (Zackrisson et al. 2011) that vary in initial mass function (IMF), metallicity, and SFH.

The works of Bouwens et al. (2010), Rogers et al. (2013), and Dunlop et al. (2013) discuss the possibility of bias in the UV spectral slope, mainly for fainter galaxies at higher redshifts. This bias is possible because fainter red galaxies may not satisfy sample selection criteria, while fainter blue galaxies would, fostering a blue bias, particularly for UV-faint galaxies.

We also conducted tests to evaluate the impact of introducing a free f_esc parameter, with a flat prior between 0 and 1, to allow SEDs without a strong nebular continuum. For our sample, we run Bagpipes with and without f_esc varying, and measure β in both instances. We find that fixing f_esc = 0 versus allowing it to be a free parameter imparts a minimal change in the recovered UV spectral slope. We quantify this via a nebular reddening estimate, ${\rm{\Delta }}{\beta }_{\mathrm{neb}}={\beta }_{\mathrm{SED},\ \mathrm{corr},\ {{\rm{f}}}_{{esc}}=0}-{\beta }_{\mathrm{SED},\mathrm{corr}}=-0.01\pm 0.1$. This indicates that nebular emission does not significantly redden the slope or alter the overall SED shape for this sample. Furthermore, although this incorporation allows β values to become bluer than when not considered, as previously mentioned, we are still limited by how blue the modeling can go with and without the prior added.

In this work, we directly compare our results to the most recent works within the same redshift range as mentioned in Section 3.4. We first compare the results for β, M_UV, and stellar mass for our sample with those of Austin et al. (2023) given their overlapping scope with respect to the field and redshift range of observations. We find that their median M_UV ∼ −18.9 is slightly brighter than our sample with M_UV ∼ −18.1 and our median stellar masses are comparable at ∼10⁸ M_⊙. Any differences in M_UV could be due to their photometric data reduction (the Leung et al. 2023 reduction is significantly deeper) and/or the subsequent SED fitting to the observational data (for this analysis, we quote their values for their fits with LePhare; see their Table 2). Their average β ∼ −2.61 is bluer than ours and spans a wider range of values than our work, at β_SED ∼ −2.45, but this may be a result of their use of photometric power-law fitting as their primary approach to measuring the UV spectral slope.

The work of Topping et al. (2022) studied 123 galaxies in the JWST CEERS field at z = 7–11. These galaxies, utilizing CEERS data, do not reach depths comparable to NGDEEP, and as such, the average luminosities of these galaxies are slightly brighter than those of our sample. Here, the median UV magnitude for their sample is M_UV ∼ −19.6, corresponding to slightly redder UV slopes averaging β ∼ −2.0 (measured by photometric power-law fitting). Cullen et al. (2023b) use data from z = 8 to 16 with JWST SMACS J0723, GLASS, and CEERS surveys along with the COSMOS/UltraVISTA survey. Their work does not note any stellar mass values for their sample of galaxies. However, the diversified survey inputs indicate an average UV magnitude marginally brighter than our data set, at M_UV ∼ −19.3, and a resulting median UV slope of β ∼ −2.3 (also measured by photometric power-law fitting; SED fitting done with EAZY to obtain M_UV).

The work of Tacchella et al. (2022) aimed primarily to study high-mass galaxies at z > 9, at ∼10⁹–10¹⁰ M_⊙, in CANDELS (Grogin et al. 2011; Koekemoer et al. 2011). These galaxies were also primarily brighter and redder than our sample at M_UV ∼ −21.5 with a median β ∼ −1.8 when measured via SED fitting, and others we compared to, but offered a bright and high-mass end to any trends in the UV spectral slope with M_UV and stellar mass. With this addition to the study, we could primarily note that our lower-mass galaxies have significantly bluer UV slopes than their work. This contrast in the characteristics of the two samples highlights the nuanced nature of the relationship between β, M_UV, and stellar mass.

However, it is essential to note that the presence of higher-mass and redder slopes in the Tacchella et al. (2022) sample does not necessarily imply a strong, straightforward trend between β and M_UV at high redshift. This point is corroborated by the findings of Bouwens et al. (2012), who identified a luminosity-dependent relationship between M_UV and β in the redshift range z = 4–7 when defined at slightly redder wavelengths, i.e., λ_rest ∼ 2000 Å. Their study demonstrated that trends between these parameters could be influenced by the specific wavelength at which M_UV is defined. In our analysis and the comparisons we made, M_UV was consistently defined at the rest wavelength of ∼1500 Å, and at high redshifts, we have yet to observe strong correlations between M_UV and the UV spectral slope.

Although we expect β to decrease as we move earlier in cosmic history, where, on average, for the whole population, younger, smaller galaxies with just a few generations of stellar populations existed, the weak trends with redshift we see here could be due to the sample size for each redshift bin we have in our population. When we only look at our sample, there is no significant correlation between β and redshift, while one emerges when combined with other works. Future analyses with larger sample sizes (such as the recent works of Cullen et al. 2023b) that span both UV-bright and faint galaxies and wider ranges of redshift will be more capable of probing any evolution of the UV spectral slope, as well as minimizing the effects of cosmic variance. We also note that previous works, e.g., Topping et al. (2022), Austin et al. (2023), and Cullen et al. (2023b), estimate the UV spectral slope by applying a power-law fit to the photometry. While this is a reasonable approach, as we have discussed, this method can be less reliable when photometry is limited. While our preferred SED method has reduced scatter, it is not free of its own biases. Improvement can be made via higher spectral resolution data (e.g., medium-band imaging or prism spectroscopy) with JWST in the rest-frame UV regime of these galaxy SEDs.

Here we compare the results of our observations and the other observational data sets utilized for this work alongside three cosmological simulations: the Santa Cruz semianalytic model (SC-SAM; Yung et al. 2024), the First Light And Reionisation Epoch Simulations (FLARES; Lovell et al. 2021; Vijayan et al. 2021; Wilkins et al. 2023), and the SIMBA-EoR simulations (Davé et al. 2019; E. Jones et al. 2024, in preparation).

We match the results from our observations with those found in other studies, using the predictions given by the Santa Cruz semianalytic model (SC-SAM; Yung et al. 2024) at z = 9 (compared to galaxies in redshift bin z = 8.5–9.5) and z = 11 (compared to galaxies from z = 9.5 to 12, as binned in Leung et al. 2023). These predictions are made with halo merger trees extracted from the suite of dark-matter-only simulations gureft, which provide robust halo merger histories with 170 snapshots stored between 40 ≳ z ≳ 6 (Yung et al. 2023). The SC-SAM tracks the full star formation and chemical enrichment histories of predicted galaxies under the influence of a set of key physical processes (see Somerville et al. 2015; Yung et al. 2019, 2022 and references therein) and couples them with the stellar population synthesis (SPS) model Binary Population and Spectral Synthesis (BPASS; Eldridge et al. 2017; Byrne et al. 2022) to produce high-resolution SEDs for each galaxy.

We also compare observations to the predictions from the First Light And Reionisation Epoch Simulations (FLARES; Lovell et al. 2021; FLARES-II, FLARES-V) at z = 9. FLARES resimulates 40 regions selected from a large (3.2 Gpc)³ dark-matter-only simulation using a variant of the eagle (Schaye et al. 2015; Crain et al. 2015) physics model. By simulating a wide range of environments and statistically combining them, flares is able to probe a larger dynamic range of galaxies than possible with comparable volume periodic simulations.

Finally, we compare all observational data sets to the predictions from the Simba-EoR simulations (Davé et al. 2019; E. Jones et al. 2024, in preparation). Simba-EoR builds on the Simba model (Davé et al. 2019), adding a sophisticated subgrid model for the interstellar medium that directly tracks H₂ and dust coevolution via a fully coupled chemical network modulated by a local interstellar radiation field and run at higher resolution. This results in significantly more early star formation versus Simba, leading to good agreement with the z > 9 UV luminosity function (Finkelstein et al. 2023).

For each galaxy in the simulations, we make use of their UV spectral slope, stellar mass, and UV magnitude. In Figure 8, we compare our work and that of previous works utilized in earlier sections to the gureft, flares, and Simba-EoR models. We find that the measurements of β_SED, M_UV, and M_* as derived from Bagpipes for our sample match well with predictions in both redshift bins in comparison to gureft predictions (the brighter, more massive galaxies in our z = 9 and 11 samples are also in agreement with the flares and Simba-EoR simulations). We note that although with SED fitting estimating stellar masses can yield high uncertainties (especially when photometrically limited in the optical regime; Endsley et al. 2023) at such high redshifts, galaxies are still fairly young when compared to those at lower redshifts and as such do not necessarily have such old stellar populations/high stellar masses (see Section 5.2 of Finkelstein et al. 2010). With the photometric data available, we show in Figure 8 that the estimated stellar masses from SED fitting for each individual galaxy reliably align with simulations, as well as with other observations.

Zoom In Zoom Out Reset image size

Notably, we do not observe any ultrablue UV spectral slopes (β < −3), which is in line with these simulations not predicting such extreme values for the galaxy redshifts and masses we examine. Although they are also in agreement with our sample, the gureft and flares simulations do not consider any Population III stellar grids, and our observations do not drive this inclusion for galaxies at the redshift and mass range considered here. Future works incorporating these Population III stellar grids, such as the Yggdrasil models (Zackrisson et al. 2011), will be necessary to explain any potential future observations of β < −3.

JWST, which, as of now, can detect candidate galaxies out to z ∼ 16 as shown in this work and others, may provide insights into the early universe and the existence of some of the earliest generations of stars. However, it remains significantly challenging to distinguish between Population III stars and low-metallicity Population II stars owing to the absence of distinctive photometric features. While potential signatures might be discovered in future works through spectroscopic analysis, observing at substantial depths will be needed to acquire this detailed information.

Although direct observation of Population III stars remains challenging, the galaxies we observe at high redshift are not without their importance. Many of these galaxies are likely low in metallicity, as inferred from the relatively blue average values we have obtained for β and by the evidence for low levels of dust attenuation and low stellar mass, providing a view of the universe when it was younger and less chemically enriched. For example, the local galaxy NGC 1705 has β ∼ −2.5 (Calzetti et al. 1994), similar to our median value, with a measured gas-phase abundance of ∼20% solar. We thus expect that the typical galaxy in our sample is composed of mainly Population II stars, like galaxies we observe at lower redshifts (Greif & Bromm 2006). However, their low metallicities mean that they can serve as analogs for the conditions under which Population III stars might have formed. This contrasts slightly with findings by Pérez-González et al. (2023). They observed a heightened activity in the universe during its early phases, especially around z ≳ 11. Such elevated activity, as they state, could imply a notable presence of Population III stars, as these primordial stars might contribute significantly to the heightened UV photon production. However, our data indicate a more prevalent role for Population II stars. Together, these studies raise intriguing questions about the early universe's stellar populations and their contributions to its early luminosity.

Building on this, the characteristics of galaxies housing different stellar populations are anticipated to vary across several parameters. Here we discuss the differences in what results are expected from Population III stars and what we are obtaining in this work.

Dust attenuation is expected to be zero in galaxies with Population III stars (Bromm et al. 1999). Given that these stars are thought to be the first generation and formed in a dust-free (and metal-free) environment, the corresponding galaxies would have significantly less dust than those with Population II stars along the line of sight. Although with Bagpipes we are obtaining A_v values on the order of zero, based on our values of β, our galaxies are consistent with having low- but not zero-metallicity stellar populations. In such cases, follow-up spectroscopy is vital to verify such situations (e.g., measurements of the rest-frame UV spectral slope directly from the spectra; Balmer decrements are unlikely to be detected at these high redshifts currently).

In terms of age, galaxies with Population III stars are expected to be among the youngest in the universe (existing at high redshifts), as these stars are theorized to have formed within a few hundred million years following the big bang (Bromm et al. 1999; Greif & Bromm 2006). Deep imaging and spectroscopic surveys like NGDEEP are perhaps our best chance at catching a glimpse of these stellar populations, or at least a direct descendant of them (i.e., Population II stars enriched by Population III remnants). Our galaxy sample has ages that span on the order of 1–100 Myr, with the average age being $\mathrm{log}(\mathrm{Age})={7.8}_{-0.8}^{+0.2}$ yr old. The work of Zackrisson et al. (2011) studies the SED tracks of Population I, II, and III (III.1, III.2, and III (Kroupa IMF); see Section 2 of their work for more information on the different Population III types). Defining stellar populations by Type A, B, and C (where Type A has an f_esc = 0, Type B an 0 < f_esc < 1, and Type C an f_esc = 1), we compare the UV slopes of our galaxy sample versus what they predict different stellar populations have when binned by age at 1, 10, and 100 Myr (see their Table 2 for a full list of UV slopes for instantaneous-burst models). Our galaxy sample is fully consistent with that of Population II tracks for both Type A and Type C, where Type A tracks range from −2.2 < β < −2.0 and Type C tracks range from −3.1 < β < −2.0 for 1–100 Myr. We find some evidence that our sample follows that of the less extreme Population III type cases when binned by age due to the uncertainty in f_esc (these UV slopes span −3.5 < β < −2.3 depending on Population III type, stellar population type, and age).

We thus conclude that the galaxies we have observed with NGDEEP do not show evidence of harboring ultra–low-metallicity and/or Population III stars. While we lack current observational data for galaxies harboring Population III stars, theoretical predictions and simulations offer valuable insight for distinguishing potential differences when comparing them to galaxies inhabited by Population I and II stars. The forthcoming observations from JWST will play a pivotal role in validating these predictions, thus enhancing our comprehension of the early universe. Despite the fact that this study has not confirmed the presence of such exotic stellar populations, there remains a possibility that they exist at higher redshifts and/or perhaps exhibit fainter luminosities than what we can detect at this moment.