Ionizing properties of galaxies in JADES for a stellar mass complete sample: resolving the cosmic ionizing photon budget crisis at the Epoch of Reionization

ABSTRACT

We use NIRCam imaging from the JWST Advanced Deep Extragalactic Survey (JADES) to study the ionizing properties of a sample of 14 652 galaxies at |$3 \le z_{\rm {phot}} \le 9$|⁠, 90 per cent complete in stellar mass down to log(M|$_{\star }$|/[M|$_{\odot }])\approx 7.5$|⁠. Out of the full sample, 1620 of the galaxies have spectroscopic redshift measurements from the literature. We use the spectral energy distribution fitting code Prospector to fit all available photometry and infer galaxy properties. We find a significantly milder evolution of the ionizing photon production efficiency (⁠|$\xi _{\rm {ion}}$|⁠) with redshift and UV magnitude than previously reported. Interestingly, we observe two distinct populations in |$\xi _{\rm {ion}}$|⁠, distinguished by their burstiness (given by SFR|$_{10}$|/SFR|$_{100}$|⁠). Both populations show the same evolution with z and M|$_{\rm {UV}}$|⁠, but have a different |$\xi _{\rm {ion}}$| normalization. We convolve the more representative |$\log (\xi _{\rm {ion}} (z,\text{M}_{\rm {UV}}))$| relations (accounting for |$\sim 97$| % of the sample), with luminosity functions from literature, to place constraints on the cosmic ionizing photon budget. By combining our results, we find that one of our models can match the observational constraints from the Ly |$\alpha$| forest at |$z\lesssim 6$|⁠. We conclude that galaxies with M|$_{\rm {UV}}$| between |$-16$| and |$-20$|⁠, adopting a reasonable escape fraction, can produce enough ionizing photons to ionize the Universe, without exceeding the required ionizing photon budget.

1 INTRODUCTION

The Epoch of Reionization (EoR) is one of the major phase transitions of the Universe, when it went from being dark and neutral to highly ionized, allowing Lyman Continuum (Ly C; with energies above 13.6 eV) radiation to travel through the intergalactic medium (IGM). Observations place the end of this epoch at |$z\sim 6$| (Becker et al. 2001; Fan et al. 2006; Yang et al. 2020), with some favouring a later reionization at |$z\sim 5$| (Keating et al. 2020; Bosman et al. 2022; Zhu et al. 2024). There is a debate regarding which sources are the main responsible agents in ionizing the Universe. The community widely agrees that young massive stars within galaxies are key, since they produce copious amounts of ionizing photons, which might be able to escape the interstellar medium (ISM), and eventually ionize the IGM (Hassan et al. 2018; Rosdahl et al. 2018; Trebitsch, Volonteri & Dubois 2020). However, the nature of the galaxies that drive reionization: bright and massive, faint and low-mass, or a combination of them, is still uncertain (Finkelstein et al. 2019; Naidu et al. 2020; Robertson 2022; Yeh et al. 2023). Moreover, it is unclear how much active galactic nuclei (AGNs) contribute to reionization (Dayal et al. 2020; Maiolino et al. 2024a; Madau et al. 2024).

The stellar mass of galaxies has been seen to correlate with how efficiently ionizing photons are produced (Simmonds et al. 2024a). Simulations indicate that it also relates how these ionizing photons escape (Paardekooper, Khochfar & Dalla Vecchia 2015). The latter is measured through their Lyman Continuum escape fractions (f|$_{\rm {esc}}$|(Ly C)), defined as the ratio between the H-ionizing radiation that is emitted intrinsically, and that which reaches the IGM. In order for galaxies to account for the reionization of the Universe, either a significant average f|$_{\rm {esc}}$| value is required (f|$_{\rm {esc}}$| = 10–20 per cent; Ouchi et al. 2009; Robertson et al. 2013, 2015; Finkelstein et al. 2019; Naidu et al. 2020), or a high ionizing photon production efficiency. These ranges of f|$_{\rm {esc}}$| have been observed for individual star-forming galaxies at |$z \lesssim 4$| (e.g. Borthakur et al. 2014; Bian et al. 2017; Vanzella et al. 2018; Izotov et al. 2021), but not usually in large numbers (Leitet et al. 2013; Leitherer et al. 2016; Steidel et al. 2018; Flury et al. 2022). An alternative to high escape fractions is a high ionizing photon production efficiency (⁠|$\xi _{\rm {ion}}$|⁠), given by the ratio between the rate of ionizing photons being emitted (⁠|$\dot{n}_{\rm {ion}}$|⁠), and the monochromatic non-ionizing ultra-violet (UV) luminosity density. Indeed, observational studies up to |$z\sim 9$| have found |$\xi _{\rm {ion}}$| to increase as a function of redshift (e.g. Endsley et al. 2021; Stefanon et al. 2022; Atek et al. 2024; Simmonds et al. 2023; Tang et al. 2023; Harshan et al. 2024; Pahl et al. 2024; Saxena et al. 2024; Simmonds et al. 2024a).

The behaviour of |$\xi _{\rm {ion}}$| as a function of redshift has important consequences on the cosmic budget of reionization (e.g. Muñoz et al. 2024), defined as the number of ionizing photons produced per comoving volume unit of the Universe (⁠|$\dot{N}_{\rm {ion}}$|⁠). Three ingredients must be provided in order to study |$\dot{N}_{\rm {ion}}$|⁠: (1) a prescription for f|$_{\rm {esc}}$|⁠, (2) a UV luminosity density function, |$\rho _{\rm {UV}}$|⁠, describing the number of galaxies per unit volume that have a given UV luminosity, as a function of redshift (e.g. Bouwens et al. 2021; Adams et al. 2024; Donnan et al. 2024; Robertson et al. 2024, Whitler et al., in preparation), and (3) constraints on |$\xi _{\rm {ion}}$|⁠. In addition, the IGM clumping factor of the Universe has to be considered (e.g. Madau, Haardt & Rees 1999; Kaurov & Gnedin 2014; So et al. 2014). This factor is a measure of the uniformity of the matter distribution in the Universe, and has crucial implications on reionization since it relates to the amount of atomic recombinations taking place in the IGM. Briefly, a higher clumping factor implies that more ionizing photons need to be emitted per unit volume at a given redshift in order to sustain hydrogen ionization.

Before the launch of the JWST (Gardner et al. 2023), it was common practice to set f|$_{\rm {esc}}$| and |$\xi _{\rm {ion}}$| as constants. Fortunately, the JWST has given us an unprecedented view of the early Universe at rest-frame optical wavelengths, which has allowed us to place better constraints on |$\xi _{\rm {ion}}$|⁠. In particular, Simmonds et al. (2024a) studied a sample of emission-line galaxies (ELGs) at |$z\sim 4\!-\!9$| using photometry obtained with the Near Infrared Camera (NIRCam; Rieke et al. 2023), on board the JWST. The mean redshift of reionization is at |$z = 7.68\pm 0.79$| (Planck Collaboration VI 2020), meaning that this study (and the many others enabled by JWST) probe deep into the EoR. Through H |$\alpha$| and [O iii] emission line fluxes, Simmonds et al. (2024a) estimated |$\xi _{\rm {ion}}$| for a sample of 677 galaxies. In parallel, they inferred the same quantity by using the spectral energy distribution (SED) fitting code Prospector (Johnson et al. 2019, 2021). They find that the |$\xi _{\rm {ion}}$| measurements estimated by emission line fluxes agree with the values obtained by Prospector. Additionally, they conclude that |$\xi _{\rm {ion}}$| increases with redshift, and that this increase is due to low-mass faint galaxies having more bursty star formation histories (SFHs). Here, burstiness is quantified by the ratio between recent (averaged over 10 Myr) and past (averaged over 100 Myr) star formation rates (SFRs), which is associated with low stellar masses (Weisz et al. 2012; Guo et al. 2016; Looser et al. 2023a), mainly due to the increased importance of stellar feedback. At high redshifts, however, burstiness can also be explained by the imbalance between gas accretion and supernovae (SNe) feedback time-scales, which prevent star formation equilibrium in the ISM (Faucher-Giguère 2018; Tacchella, Forbes & Caplar 2020). We note that the SFR|$_{10}$|/SFR|$_{100}$| ratio is a direct measure of the recent SFH, and its variance for an ensemble of galaxies measures short-term star formation variability (‘burstiness’; Caplar & Tacchella 2019).

The sample constructed in Simmonds et al. (2024a) suffered from one main limitation: since emission lines were required, and with a sufficient strength so that they were measurable from photometry, this sample was biased towards star-forming galaxies with significant H |$\alpha$| and/or [O iii] emission. In fact, Laseter et al. (in preparation) demonstrate |$\xi _{\rm {ion}}$| is consistently high (log(⁠|$\xi _{\rm {ion}}$|/[Hz erg|$^{-1}$|]) |$\approx 25.5$|⁠) down to [O iii] equivalent widths of 200Å due to low metallicities (⁠|$Z \lesssim 1/10 Z_{\odot }$|⁠), further demonstrating this bias from the past results. Given the agreement found between Prospector and the emission line measurements of |$\xi _{\rm {ion}}$|⁠, in this work, we use Prospector to fit the full JWST Advanced Deep Extragalactic Survey (JADES; Bunker et al. 2024; Eisenstein et al. 2023a) photometry set for a sample of JADES galaxies in the Great Observatories Origins Deeps Survey South (GOODS-S; Giavalisco et al. 2004). Our sample is 90 per cent complete in stellar mass down to masses of log(M|$_{\star }$|/[M|$_{\odot }])\sim 7.5$|⁠, providing us a deep statistical view of the ionizing properties of galaxies.

The structure of this paper is the following. In Section 2, we present the data used in this work, along with the sample selection criteria. In Section 3, we present our Prospector fitting method. Some general properties of the sample are given in Section 4, followed by our constraints on the ionizing properties of galaxies in Section 5. Implications for reionization are discussed in Section 6, while the caveats and limitations of our methods are discussed in Section 7. Finally, brief conclusions are presented in Section 8. Throughout this work we assume |$\Omega _0 = 0.315$| and |$H_0 = 67.4$| km s|$^{-1}$| Mpc|$^{-1}$|⁠, following Planck Collaboration VI (2020).

2 DATA AND SELECTION CRITERIA

In this section, we describe the data and selection criteria applied to build our sample of galaxies, with the goal to construct a stellar mass complete sample between redshifts 3 and 9.

2.1 Data

We use the full JADES (Bunker et al. 2024; Eisenstein et al. 2023a) photometry set in the GOODS-S region, including the publicly available NIRCam Deep imaging (Rieke et al. 2023), and the JADES Origins Field (JOF; Eisenstein et al. 2023b), covering an area of |$\sim 45$| arcmin|$^2$| with an average exposure time of 130 h. When available, we also use photometry from the JWST Extragalactic Medium-band Survey (JEMS; Williams et al. 2023), and from the First Reionization Epoch Spectroscopic Complete Survey (FRESCO; Oesch et al. 2023).

2.1.1 Photometry

The photometric catalogue used in this work has been produced in the same way as the one used in Simmonds et al. (2024a). In brief, the source detection and photometry leverage both the JEMS NIRCam medium band and JADES NIRCam broad and medium band imaging. Detection was performed using the photutils (Bradley et al. 2022) software package, identifying sources with contiguous regions of the SNR mosaic with signal |$\gt 3\sigma$| and five or more contiguous pixels. We also use photutils to perform circular aperture photometry with filter-dependent aperture corrections based on model point-spread-functions following the method of Ji et al. (2024), as described in Robertson et al. (2024). In addition to the NIRCam observations, HST images from the Hubble Legacy Field programme (Illingworth et al. 2016; Whitaker et al. 2019), and the Cosmic Assembly Near-Infrared Deep Extragalactic Legacy Survey (CANDELS; Grogin et al. 2011; Koekemoer et al. 2011), as well as reductions of GOODS-S from Giavalisco et al. (2004), are used. The details of the catalogue generation and photometry will be presented in Robertson et al. (in preparation). In this work, we use a Kron aperture placed on images that have been convolved to a common resolution, and impose an error floor of 5 per cent in each band. Our photometry does not take include the eazy derived photometric offsets, however, we find that the NIRCam photometric offsets are small given the uncertainty of the photometry (see Appendix  A). The throughputs of the filters used in this work are shown in Fig. 1, as well as emission lines of interest and how their observed wavelength evolves with redshift.

2.1.2 Redshifts

Due to the richness of the photometry set, photometric redshifts (⁠|$z_{\rm {phot}}$|⁠) can be obtained with great accuracy (Hainline et al. 2024). When fitting SEDs in this work, we use the |$z_{\rm {phot}}$| inferred by the template-fitting code eazy (Brammer, van Dokkum & Coppi 2008), as described in Hainline et al. (2024) and Rieke et al. (2023). Additionally, when available, we make use of spectroscopic redshifts (⁠|$z_{\rm {spec}}$|⁠) from the Near-Infrared Spectrograph (NIRSpec; Jakobsen et al. 2022), as well as those reported in literature (Puskás et al. in preparation), including the JADES NIRSpec redshifts from Bunker et al. (2024) and D’Eugenio et al. (2024). In particular, they have been compiled from the Atacama Large Millimeter/submillimeter Array (ALMA) Spectroscopic Survey in the Hubble Ultra-Deep Field (ASPECS; Walter et al. 2016), CANDELS (Grogin et al. 2011), the 3D-Hubble Space Telescope (HST) Survey (Momcheva et al. 2016), the Multi Unit Spectroscopic Explorer (MUSE; Bacon et al. 2010) Ultra-Deep Field DR2 (Bacon et al. 2023), and a redshift catalogue (F. Sun private communication) produced from grism data of the First Reionization Epoch Spectroscopic Complete Survey (FRESCO; Oesch et al. 2023).

We selected galaxies that were matched to JADES NIRCam positions within 0.25 arcsec, and that were flagged as reliable by each team. Since galaxies often have multiple |$z_{\rm {spec}}$| measurements, we defined four categories to collate the several redshift measurements into one, which we call |$z_{\rm {best}}$|⁠:

• Category A: Only one redshift labelled as having the highest quality. We use this redshift as |$z_{\rm {best}}$|⁠.

• Category B: Multiple redshifts with the same (highest) quality that agree when rounded to the second decimal. We use this rounded value as |$z_{\rm {best}}$|⁠.

• Category C: Multiple solutions with the same (highest) quality that have a non-dramatic disagreement (with a difference smaller than 1). We define |$z_{\rm {best}}$| as the mean between the solutions and add errors reflecting the difference.

• Category D: Multiple solutions with the same (highest) quality that have a dramatic disagreement (⁠|$\Delta \gt 1$|⁠). In these cases, we keep all the highest quality redshifts and follow up with a visual inspection.

2.2 Sample selection criteria

In order to build a sample as complete as possible, we only impose two conditions that have to be met: (1) an S/N of at least 3 in the F444W band, and (2) redshifts |$3 \le z \le 9$|⁠. Based on the nature of the redshift used (photometric or spectroscopic), we construct two samples, that we now introduce.

The photometric sample is initially composed of 37 272 galaxies. We use the SED fitting code prospector (Johnson et al. 2019, 2021) to fit the full sample (see Section 3), yielding a total sample of 35 442 galaxies. The cases that failed (⁠|$\sim 5$|%) were due to either poor photometric coverage (⁠|$\lt 6$| NIRCam photometric points) or to them being false detections (e.g. a diffraction spike). In order to have reliable inferred galaxy properties, we only use the results which have a reduced |$\chi ^2 \le 1$|⁠, resulting in a final photometric sample with 25 319 galaxies. The completeness in stellar mass and UV magnitude is discussed in Section 4. We note that if we extend the limit up to |$\chi ^2$| = 10, we obtain 29 323 galaxies instead. However, the final stellar mass complete sample is virtually unchanged from the one discussed in this work. The spectroscopic sample is composed of 1620 galaxies, all of which are also part of the photometric sample. Their redshift classification is as follows: 1234 in Category A, 363 in category B, 12 in category C, and 11 in category D.

3 SED FITTING WITH PROSPECTOR

As demonstrated in Simmonds et al. (2024a), the ionizing properties inferred with the SED fitting code prospector (Johnson et al. 2019, 2021) are in good agreement with those obtained by emission line fluxes, when such fluxes are detected. Therefore, in this work, we fit the entirety of our samples with prospector, without introducing an additional selection bias (i.e. by only selecting emission line galaxies).

prospector uses photometry and/or spectroscopy as an input in order to infer stellar population parameters, from UV to IR wavelengths. We use photometry from the HST ACS bands: F435W, F606W, F775W, F814W, F850LP, and from the HST IR bands: F105W, F125W, F140W, and F160W. In addition, we use the JADES NIRCam photometry from: F090W, F115W, F150W, F162M, F200W, F250M, F277W, F300M, F335M, F356W, F410M, and F444W. Finally, when available, we include JEMS photometry: F182M, F210M, F430M, F460M, and F480M. The same Kron convolved aperture is used to extract the HST, JADES and JEMS photometry.

For the photometric redshift sample, we adopt a clipped normal distribution using the eazy|$z_{\rm {phot}}$| as the redshift mean, with the sigma given by the |$z_{\rm {phot}}$| errors. Whereas for the spectroscopic redshift sample we fix the redshift to |$z_{\rm {spec}}$| for galaxies in Categories A, B, and D, and use the same prior as for the photometric sample for galaxies in Category C.

We vary the dust attenuation and stellar population properties following Tacchella et al. (2022a). In summary, we use a two component dust model described in Charlot & Fall (2000) and Conroy, Gunn & White (2009). This model accounts for the differential effect of dust on young stars (⁠|$\lt 10$| Myr) and nebular emission lines, through different optical depths and a variable dust index (Kriek & Conroy 2013). We adopt a Chabrier (Chabrier 2003) initial mass function (IMF), with mass cutoffs of 0.1 and 100 M|$_{\odot }$|⁠, respectively, allowing the stellar metallicity to explore a range between 0.01–1 Z|$_{\odot }$|⁠, and include nebular emission. The continuum and emission properties of the SEDs are provided by the Flexible Stellar Population Synthesis (FSPS) code (Byler et al. 2017), based on cloudy models (v.13.03; Ferland et al. 2013) using MESA Isochrones & Stellar Tracks (MIST; Choi et al. 2016; Dotter 2016), and the MILES stellar library (Vazdekis et al. 2015). We note that he UV extension of the MILES library is based on the Basel Stellar Library (BaSeL; Lastennet et al. 2002). These cloudy grids introduce an upper limit on the permitted ionization parameters (log|$\langle U \rangle$||$_{\rm {max}} = -1.0$|⁠). We briefly remark here that this upper limit might not be appropriate for high-redshift galaxies (see e.g, Cameron et al. 2023). Due to the stochastic nature of the IGM absorption, we set a flexible IGM model based on a scaling of the Madau model (Madau 1995), with the scaling left as a free parameter with a clipped normal prior (⁠|$\mu = 1.0, \sigma = 0.3$|⁠, in a range [0.0, 2.0]). Last but not least, we use a non-parametric SFH (continuity SFH; Leja et al. 2019). This model describes the SFH as eight different SFR bins, the ratios and amplitudes between them are in turn, controlled by the bursty-continuity prior (Tacchella et al. 2022b). For a general view of the goodness of fits of our stellar mass complete sample (described in Section 4.2), we direct the reader to Appendix  A, where we show a comparison between the modelled and observed photometry, as a function of the bands being used.

Fig. 2 shows example SEDs and best-fitting spectra for our sample. From top to bottom the redshift increases from 3.5 to 8.5, as indicated on the top left of each panel. The markers show HST (triangles) and JWST NIRCam (circles) photometry with their corresponding errors. The purple shaded area corresponds to the rest-frame spectral region at |$\lambda _{\rm {rest-frame}} = 1250\!-\!2500$| Å, used to estimate the rest-frame UV continuum slope (⁠|$\beta$|⁠; Calzetti, Kinney & Storchi-Bergmann 1994), defined as F|$_{\lambda } \propto \lambda ^{\beta }$|⁠. We obtain |$\beta$| by fitting a line to the best-fit SED provided by prospector. We stress that there are important limitations to the measurements of |$\beta$| without spectra (for a full review, see Austin et al. 2024). Finally, the vertical dashed lines show the observed wavelength of Ly |$\alpha$|⁠, [O ii], [O iii], and H |$\alpha$|⁠. We note that not all galaxies have obvious emission lines detected in their photometry, illustrating the advantage of our sample selection (see for example, the second panel: JADES-GS + 53.06045−27.82581).

4 SAMPLE GENERAL PROPERTIES

In this section, we discuss some of the general properties of our samples. We first compare the photometric redshifts inferred by eazy to those obtained with prospector, and to the spectroscopic redshifts (when available). We then describe the stellar mass completeness of our sample.

4.1 Redshift comparisons

As previously mentioned, we use eazy-inferred redshifts as priors when fitting the photometric sample. The redshifts obtained by this template-fitting tool have proven to be reliable when using the full JADES NIRCam photometry set (see fig. 13 of Rieke et al. 2023). In Fig. 3, we compare the eazy and the prospector redshifts for the photometric sample. The prospector redshifts are better constrained, as seen by the median error bars (black point). This is not surprising since prospector is already using the eazy results as a prior on the redshift. We note that prospector is more flexible than eazy, since the latter uses a fixed set of templates. Although, we note that the linear combination of templates used by eazy might be outside of the prospector parameter space, or disfavoured by the SFH or other prospector priors. It can also be seen that a lower-redshift solution is preferred for several sources. However, the distribution overall follows a one-to-one relation with a best-fitting slope of |$0.96 \pm 0.01$| (black dashed line), demonstrating that both methods are in general agreement. We note that the number of sources decreases significantly as a function of photometric redshift: 10 215 at |$3 \le z \lt 4$|⁠, 6005 at |$4 \le 5$|⁠, 5203 at |$5 \lt z \le 6$|⁠, 3037 at |$6 \lt z \le 7$|⁠, 651 at |$7 \lt z \le 8$|⁠, and, 208 at |$8\lt z \le 9$|⁠.

In order to test how well the |$z_{\rm {phot}}$| retrieve real (⁠|$z_{\rm {spec}}$|⁠) redshifts, in Fig. 4, we compare the photometric redshifts (from both eazy and prospector) to the spectroscopic sample in Categories A and B. The spectroscopic sample is biased towards brighter galaxies with stronger emission lines, compared to the full photometric sample, so we would expect more accurate |$z_{\rm {phot}}$| estimations. We find a good agreement between |$z_{\rm {phot}}$| and |$z_{\rm {spec}}$| (i.e. |$|z_{\rm {spec}}-z_{\rm {phot}}|\lt 0.15$|⁠), with only a small fraction of outliers. Specifically, 44 (35) for prospector (eazy) derived redshifts, corresponding to |$\sim 3.1$|% (⁠|$\sim 2.5$|%) of the subsample, where photometry alone makes it hard to distinguish a spectral break from another. From this comparison, we cannot say if one code performs better than the other, but we can conclude that both retrieve the correct redshift in the majority of the cases. We highlight that the SED modelling uncertainties and redshift variations are self-consistent.

Given the agreement between |$z_{\rm {phot}}$| and |$z_{\rm {spec}}$|⁠, in the remainder of this work, we use the photometric sample unless explicitly stated.

4.2 Completeness of sample

To estimate the stellar mass completeness of our sample, we use the redshifts and stellar masses inferred with prospector, noting that the spectroscopic sample overlaps with the photometric one, and that the redshifts are in tight agreement. A corner plot showing the mean shape of the posteriors for these parameters can be found in Appendix  B.

To assess the 90 per cent stellar mass completeness limit of our photometric sample, we follow the procedure described in section 5.2 of Pozzetti et al. (2010). In summary, for every redshift, we define a minimum mass (M|$_{\rm {min}}$|⁠), above which all types of galaxies can potentially be observed. To obtain M|$_{\rm {min}}$|⁠, we first need to calculate the limiting stellar mass (M|$_{\rm {lim}}$|⁠) for each galaxy, given by:

where M|$_{\star }$| is stellar mass in units of solar masses. M|$_{\rm {lim}}$| represents the mass a galaxy would have if its apparent magnitude (m) were equal to the limiting magnitude of the survey (m|$_{\rm {lim}}$|⁠) in the F444W band. This band was chosen since it has the longest effective wavelength, and thus, is the tracer of stellar mass. We divide our data into three depths depending on the exposure time: medium (T|$_{\mathrm{exp}}$| < 25 ks), deep (25 ks |$\le$| T|$_{\mathrm{exp}}$| < 65 ks), and ultra-deep (T|$_{\mathrm{exp}}$| > 65 ks), with 5|$\sigma$| flux depths of 6, 4.5, and 2.65 nJy, respectively. Once M|$_{\rm {lim}}$| has been calculated for every galaxy, we compute M|$_{\rm {min}}$| for each field by selecting the faintest 20 per cent of the galaxies with a given T|$_{\mathrm{exp}}$|⁠, and finding the limit above which 90 per cent of the selected M|$_{\rm {lim}}$| values lie. Fig. 5 shows the stellar mass of our sample as a function of redshift (grey circles), the coloured small circles show the M|$_{\rm {lim}}$| of the faintest 20 per cent of galaxies for each field, and the filled lines denote the 90 per cent completeness as a function of redshift. We find that, in general, our photometric sample is stellar mass complete down to log (M|$_{\star }$|/[M|$_{\odot }])\approx$|7.5, with slight variations depending on the depth of the observations (within |$\sim 0.2$| dex).

In the remainder of this study, we focus only on sources that have a prospector-derived stellar mass with a median above the mass completeness limit (log (M|$_{\star }$|/[M|$_{\odot }])\approx$|7.5), yielding a sample of 14 652 galaxies, which are UV complete down to M|$_{\rm {UV}}\approx -16$|⁠.

4.3 M|$_{\rm {UV}}$|–M|$_{\star }$| relation

Fig. 6 shows the stellar mass of our stellar mass complete sample as a function of M|$_{\rm {UV}}$|⁠, divided by redshift bins. The numbers on the top of each panel indicate the number of galaxies in any given bin. We fit the observational data in the M|$_{\rm {UV}}$|–M|$_{\star }$| plane taking into account the completeness limits discussed in the previous section. Specifically, we fit the following relation

where M|$_{\star }$| is stellar mass in units of solar masses, |$\alpha$| is the slope of the M|$_{\rm {UV}}$|–M|$_{\star }$| relation and |$\log (M_{\star ,0})$| is its normalization (stellar mass at |$M_{\rm {UV}}=-19.5$| AB mag). In addition to these two parameters, we also fit for the scatter |$\sigma _{\rm M_{\rm {UV}}-M_{\star }}$| in M|$_{\rm {UV}}$| at fixed M|$_{\star }$|⁠. We then fit these three parameters to the observed distribution using the dynamic nested sampling code dynesty (Speagle 2020). In each model call, we sample the M|$_{\rm {UV}}$|–M|$_{\star }$| relation (with scatter) with 1 million galaxies assuming an M|$_{\rm {UV}}$| that follows the UV luminosity function (Bouwens et al. 2021), bin this drawn galaxy sample in the M|$_{\rm {UV}}$|–M|$_{\star }$| plane taking into account the completeness limits, and compare the normalized histogram with the observed one. The best-fitting parameters can be found in Table 1, for the highest redshift bin (⁠|$8 \lt z \le 9$|⁠) we adopt the best-fitting results from the previous redshift bin (⁠|$7 \lt z \le 8$|⁠).

From the nearly linear relation between |$\log (\text{M}_{\star })$| and |$\log (\text{SFR})$| (i.e. star-forming main sequence; Brinchmann et al. 2004; Daddi et al. 2007; Salim et al. 2007), we would expect a slope |$\alpha$| of the M|$_{\rm {UV}}$|–M|$_{\star }$| relation close to |$-0.4$|⁠, assuming a simple, linear conversion between the UV luminosity and SFR (Kennicutt 1998; Madau, Pozzetti & Dickinson 1998; Salim et al. 2007). Our steeper slope of |$\approx -0.6$| is consistent with having a higher mass-to-light ratio (⁠|$M_{\star }/L_{\rm UV}$|⁠) for more massive systems, in agreement with the findings of González et al. (2011). We further find that |$\log (M_{\star ,0})$| decreases with increasing redshift (⁠|$M_{\star ,0}\propto (1+z)^{-2.7}$|⁠), i.e. at fixed M|$_{\rm UV}=-19.5$| the typical stellar mass of the galaxies is lower at earlier cosmic times. This implies a lower |$M_{\star }/L_{\rm UV}$| toward higher redshifts, consistent with younger stellar populations toward earlier times. A combination of less dust attenuation and higher specific SFR toward higher redshifts and at fixed stellar mass could lead to this outcome. Finally, the scatter of the M|$_{\rm {UV}}$|–M|$_{\star }$| relation is substantial with |$0.6\!-\!0.7$| dex. This scatter includes both the intrinsic scatter of the M|$_{\rm {UV}}$|–M|$_{\star }$| relation and the observational uncertainty, though the latter is only of the order of 0.2 dex.

As comparison, we show the log(M|$_{\star }$|⁠)–M|$_{\rm {UV}}$| mass-to-light relations from Duncan et al. (2014), Song et al. (2016), and Tacchella et al. (2018). The Duncan et al. (2014) and Song et al. (2016) relations are derived from a combination of observations (i.e. Spitzer/IRAC and HST) and SED fitting, while the one from Tacchella et al. (2018) is based on an empirical model that links star formation in galaxies to the accretion rate of dark matter halos. With our stellar mass complete sample, we find a somewhat steeper slope than previous observational studies (but mostly |$\lt 1\sigma$| discrepant; |$\sim -0.6$| instead of |$\sim -0.4$| to |$-0.5$|⁠). This effect is partly a consequence of our sample including a larger population of fainter sources (with lower stellar masses) than the mentioned studies, but also because we consider explicitly the mass completeness of our sample. Specifically, we forward model the completeness when fitting the relation; if we do not do that, we find a shallower relation of |$\sim -0.4$|⁠. Additionally, the stellar masses in these previous studies were estimated assuming parametric SFHs, which cannot fully describe the complexities of galaxies (Lower et al. 2020), while we assume more flexible SFHs.

5 CONSTRAINTS ON THE IONIZING PROPERTIES OF GALAXIES

In this section, we provide a brief overview of the methods traditionally used to infer ionizing properties of galaxies. Followed by the trends of |$\xi _{\rm {ion}}$| and |$\dot{n}_{\rm {ion}}$| with redshift and UV magnitude, including a comparison with values from the literature. Finally, we discuss the emergence of a secondary (previously unseen) sample of galaxies, that are tentatively leaking Ly C. An excerpt of the properties used in this section is given in Table 2.

5.1 Measuring |$\xi _{\rm {ion}}$| from emission line fluxes

The ionizing photon production efficiency, |$\xi _{\rm {ion}}$|⁠, can be measured through H |$\alpha$| and/or [O iii] emission line fluxes. In order to use H |$\alpha$| as a proxy for ionizing photon production efficiency (⁠|$\xi _{\rm {ion}}$|⁠), one must first correct the H |$\alpha$| line flux for dust (for example, using an SMC attenuation curve; Gordon et al. 2003). If Case B recombination is assumed, along with the assumption of no Ly C f|$_{\rm {esc}}$|⁠, and that dust has a negligible effect on Ly C photons, the dust corrected H |$\alpha$| luminosity relates directly to the rate of ionizing photons (⁠|$\dot{n}_{\rm {ion}}$|⁠, in units of s|$^{-1}$|⁠) that are being emitted (⁠|$\dot{n}_{\rm {ion}}$||$= 7.28 \times 10^{11}$| L(H |$\alpha$|⁠); Osterbrock & Ferland 2006). |$\xi _{\rm {ion}}$| is then the ratio between |$\dot{n}_{\rm {ion}}$| and the observed monochromatic UV luminosity density (measured at rest-frame |$\lambda = 1500$| Å). This method suffers from a number of of assumptions that need to be made: it has a great dependence on the chosen attenuation curve, and Case B recombination cannot always be assumed (McClymont et al. 2024; Scarlata et al. 2024). If [O iii]|$_{\lambda 5007}$| is available instead, its equivalent width can give a measure of the ionization field of the galaxy as shown in Chevallard et al. (2018), and later in Tang et al. (2019). This method also depends on the shape of the attenuation curve assumed, as well as on the gas-phase metallicity and nebular conditions.

In Simmonds et al. (2024a), the ionizing properties of a sample of 677 emission line galaxies at |$z\sim 4-9$| were analysed. These galaxies had signatures of either H |$\alpha$| and/or [O iii] emission in their NIRCam photometry, specifically with a 5|$\sigma$| detection in the filter pairs F335M–F356W and F410M–F444W. Both methods mentioned above were used to estimate |$\xi _{\rm {ion}}$|⁠, while in parallel, prospector was used to fit their full NIRCam photometry (assuming f|$_{\rm {esc}}$| = 0). A tight agreement was found between the values measured by emission line fluxes and those inferred by prospector. Therefore, in this work we rely on this SED fitting code to extract the ionizing properties of a stellar mass complete sample, allowing us to potentially study all types of galaxies, not only those with detectable emission lines. Since our prospector fitting routine assumes an escape fraction of zero, our estimated ionizing photon production efficiencies (⁠|$\xi _{\rm {ion},0}$|⁠) include a correcting factor of |$1-$|f|$_{\rm {esc}}$|⁠, such that:

Appendix  C shows how well prospector can constrain stellar masses and ionizing photon production efficiencies, as a function of flux in F444W.

5.2 Trends of ionizing properties with redshift

We first compare our results to those from the literature, in particular, for the evolution of |$\xi _{\rm {ion},0}$| with redshift. Fig. 7 encodes most results to date. Stefanon et al. (2022) compiled |$\xi _{\rm {ion},0}$| measurements up to |$z\sim 8$| (containing data from Stark et al. 2015; Bouwens et al. 2016; Mármol-Queraltó et al. 2016; Nakajima et al. 2016; Matthee et al. 2017; Stark et al. 2017; Harikane et al. 2018; Shivaei et al. 2018; De Barros et al. 2019; Faisst et al. 2019; Lam et al. 2019; Tang et al. 2019; Emami et al. 2020; Nanayakkara et al. 2020; Endsley et al. 2021; Atek et al. 2022; Naidu et al. 2022). This extensive compilation allowed them to estimate the rate of change of |$\xi _{\rm {ion},0}$| with redshift, finding it to have a slope of dlog(⁠|$\xi _{\rm {ion},0}$|⁠)/d|$z=0.09\pm 0.01$| (black dashed line). The markers show the UV-faint galaxies at |$z\sim 3\!-\!7$| from Prieto-Lyon et al. (2023), the Lyman-|$\alpha$| emitters (LAEs) from Ning et al. (2023) and Simmonds et al. (2023), the H |$\alpha$| emitters from Rinaldi et al. (2024), and the ELGs from Simmonds et al. (2024a). The latter show a slightly less steep evolution of |$\xi _{\rm {ion},0}$| with redshift than the compilation of Stefanon et al. (2022), but consistent within errors [given by dlog(⁠|$\xi _{\rm {ion},0}$|⁠)/dz|$=0.07\pm 0.02$|]. Finally, the stellar mass complete sample of this work is shown as grey pentagons. There is an overlap with previous measurements, however, there is a population of previously unseen galaxies with log(⁠|$\xi _{\rm {ion}}$|/[Hz erg|$^{-1}$|]) |$\lesssim 24.5$|⁠. The inclusion of these galaxies has the result of flattening of the increase of |$\xi _{\rm {ion},0}$| with redshift.

Fig. 8 shows |$\xi _{\rm {ion},0}$| and |$\dot{n}_{\rm {ion}}$| as a function of redshift, for our photometric (no edges) and spectroscopic (black edges) sample, colour-coded by the burstiness of their SFH (which correlates the most with |$\xi _{\rm {ion},0}$|⁠; Simmonds et al. 2024a). The crosses represent the galaxies with no recent star formation (SFR|$_{10} = 0$| M|$_{\odot }$| yr|$^{-1}$|⁠). As a reminder, both samples overlap and the prospector-inferred redshifts agree with the |$z_{\rm {spec}}$| ones, for a great majority of the overlapping galaxies. We present both samples fit separately to highlight that the increase in |$\xi _{\rm {ion},0}$| observed previously is mostly due to a selection effect (by selecting predominantly star-forming galaxies with emission lines). Indeed, the best fit to |$\xi _{\rm {ion},0}$| for the spectroscopic sample has a slope of dlog (⁠|$\xi _{\rm {ion},0}$|⁠)/dz|$\sim 0.04\pm 0.02$|⁠, within errors of the findings of Simmonds et al. (2024a). Importantly, when the full photometric sample is taken into account, the slope is considerably flatter (dlog (⁠|$\xi _{\rm {ion},0}$|⁠)/dz|$\sim -0.02\pm 0.00$|⁠). The same effect can be seen in the best-fit lines to |$\dot{n}_{\rm {ion}}$|⁠.

5.3 Trends of ionizing properties with UV magnitude

|$\xi _{\rm {ion}}$| depends on several galaxy properties, such as metallicity, age, and dust content (Shivaei et al. 2018). Moreover, fainter galaxies have been shown to be more efficient in producing ionizing radiation (e.g. Duncan & Conselice 2015; Maseda et al. 2020; Endsley et al. 2024; Simmonds et al. 2024a). Fig. 9 shows |$\xi _{\rm {ion},0}$| as a function of M|$_{\rm {UV}}$| for our data, per redshift bin. We first note that as redshift increases, there are less galaxies per bin, and they tend to be fainter. We find that by dividing the sample into ‘star-forming’ [log (SFR|$_{10}$|/SFR|$_{100}$|⁠) |$\ge$| −1] and ‘mini-quenched’ [log(SFR|$_{10}$|/SFR|$_{100}$|⁠) < −1], we distinguish two populations of galaxies that have consistent |$\xi _{\rm {ion},0}$| slopes with M|$_{\rm {UV}}$|⁠, but populate a different |$\xi _{\rm {ion},0}$| range. The star-forming sample lies above log (⁠|$\xi _{\rm {ion},0}$|/[Hz erg|$^{-1}])\sim 24.5$|⁠, while the mini-quenched (see e.g. Looser et al. 2023a,b) sample mostly lies below. The number of galaxies in each bin is shown in the top of the panels, with format N = XX (YY), where XX is the number of star-forming galaxies and YY is the number of mini-quenched galaxies. The best fit to both data sets (when enough points to fit a line reliably) are shown in the bottom of each panel. We find that in general, there is a slight increase of |$\xi _{\rm {ion},0}$| towards the fainter galaxies (slope of |$\sim 0.01\!-\!0.04$|⁠), but this increase is lower than the slope of |$\sim 0.1$| found in Simmonds et al. (2024a). This is unsurprising given the nature of our sample (i.e. containing potentially all types of galaxies). When taking these two populations in mind, the |$\xi _{\rm {ion},0}$|–z relation shown in the top panel of Fig. 8 (for the photometric sample) becomes:

for the star-forming sample. Likewise, the |$\dot{n}_{\rm {ion}}$|-z relation shown in the bottom panel of Fig. 8, for the star-forming sample becomes:

As with |$\xi _{\rm {ion},0}$|⁠, the mini-quenched galaxies also appear to populate a distinct region in the |$\dot{n}_{\rm {ion}}$| versus UV magnitude plane. Contrary to |$\xi _{\rm {ion},0}$|⁠, however, |$\dot{n}_{\rm {ion}}$| has a steeper dependence with M|$_{\rm {UV}}$| (with a slope of |$\sim -0.4$|⁠), as shown in Fig. 10.¹ In summary, as galaxies become fainter, |$\xi _{\rm {ion},0}$| marginally increases, but |$\dot{n}_{\rm {ion}}$| significantly decreases. This is an effect that derives naturally from the definition of |$\xi _{\rm {ion},0}$|⁠: the ratio between the ionizing photons being produced (⁠|$\dot{n}_{\rm {ion}}$|⁠) and the non-ionizing UV continuum luminosity. The former is dominated by the contribution of young hot stars, while the latter also includes the contribution of older stellar populations.

Finally, we combine our redshift and UV magnitude relations to perform a 2-dimensional fit and find a joint relation for the star-forming sample, given by:

5.4 Unveiling the silent population

We find a total of 350 mini-quenched galaxy candidates, corresponding to |$\sim 2.4$|% of our total sample. Their number increases with redshift up to |$z = 7$| and then decreases significantly. There are a few possible explanations for this behaviour that most likely are working together. On one hand, the preference of mini-quenched galaxies in the redshift window of |$z=4\!-\!6$| could have a physical origin (Dome et al. 2024): at higher redshifts, galaxies are bursty and undergo quenching attempts, but the replenishment time of gas within galaxies is so short that the mini-quenching phase is not observable (i.e. less than a few Myr). Supporting this scenario, recent work by Witten et al. (2024) shows a |$z\sim 7.9$| galaxy with a bursty SFH, with evidence of a mini-quenched episode lasting |$\sim 20$| Myr followed by a rejuvenation event. Towards lower redshifts, the mini-quenching phase might have a longer duration or SFHs are less bursty, leading to fewer mini-quenched systems. On the other hand, at |$z \sim 5\!-\!7$|⁠, the UV part of the spectrum is shifted into the NIRCam F090W filter, which covers our entire stellar mass complete sample (see Fig. 1). The deep NIRCam observations allow for stronger constraints on the SED shape, and could explain the increased detection of mini-quenched galaxies at these redshifts. Finally, at higher redshifts (⁠|$z \gt 7$|⁠), it becomes increasingly difficult to detect these almost featureless galaxies with photometry alone (and with a S/N > 3).

‘Mini-quenched galaxies’ have low |$\xi _{\rm {ion},0}$|⁠, well below the observed |$\xi _{\rm {ion},0}$| relations (i.e. with redshift and M|$_{\rm {UV}}$|⁠), and are likely only minor players in the reionization of the Universe. Their photometry show little-to-no evidence of emission lines, which explains why they have been neglected in previous studies of this nature. These galaxies are possibly analogues to the galaxies studied in Strait et al. (2023); Looser et al. (2023a) and Looser et al. (2023b), which are temporarily quenched due to their extremely bursty SFHs (Dome et al. 2024). Interestingly, using simulations, Dome et al. (2024) predict the number of mini-quenched galaxies increases with cosmic time, reaching |$\sim 2\!-\!4$|% at |$z=4$|⁠, in broad agreement with our preliminary findings. A spectroscopic follow-up study is needed to confirm this hypothesis.

If we adopt the |$\beta$|–f|$_{\rm {esc}}$|(Ly C) relation from Chisholm et al. (2022), given by:

we find that a significant number of the mini-quenched galaxy candidates [51 (141)] indicate leakage of ionizing radiation [f|$_{\rm {esc}}$|(Ly C) > 10 per cent (5 per cent)]. Corresponding to |$\sim 16$|% (40 per cent) of the mini-quenched sample. By comparison, the star-forming sample only has |$\sim$|4 per cent (32 per cent) of galaxies that obey the same criteria. Fig. 11 shows |$\xi _{\rm {ion},0}$| as a function of |$\beta$| for all galaxies with f|$_{\rm {esc}}$|(Ly C) > 10 per cent. The mini-quenched candidates populate a different parameter space in |$\xi _{\rm {ion},0}$| and |$\beta$| than the star-forming galaxies, with the mini-quenched galaxies having on average bluer UV continuum slopes and lower |$\xi _{\rm {ion},0}$| values. We now focus on the strongest leaker candidates within the mini-quenched sample (f|$_{\rm {esc}}$| > 10 per cent).

Observing Ly C leakage directly is impossible at high redshift, the average IGM transmission of hydrogen ionizing photons emitted at |$\lambda _{\rm {rest-frame}} \sim 900$| Å, at |$z \sim 6$|⁠, is virtually zero (Inoue et al. 2014). It is even more complicated to observe hydrogen ionizing photons emitted at shorter wavelengths, for example at |$\sim 700$| Å, where the nebular contribution does not contaminate f|$_{\rm {esc}}$| estimations (Simmonds et al. 2024b). Therefore, indirect methods are required to understand how ionizing photons escape in the early Universe. In order to confirm if these galaxies are indeed leaking, we run prospector on this subsample, with a modified approach that includes Ly C leakage in the fitting routing (Stoffers et al. in preparation). This is work in progress and is in the process of being calibrated. Promisingly, prospector finds signs of leakage for these Ly C leaking candidates (finding f|$_{\rm {esc}}$| > 20 per cent for all of them).

Fig. 12 shows |$\xi _{\rm {ion},0}$| versus f|$_{\rm {esc}}$|(LyC), colour-coded by burstiness. The shaded region shows the area where f|$_{\rm {esc}}$| > 10 per cent, and log(SFR|$_{10}$|/SFR|$_{100}$|⁠) |$\lt -1$| (the latter roughly coincides with log(⁠|$\xi _{\rm {ion},0}$|/[Hz erg|$^{-1}$|]) < 24.5). This parameter space is analogous to the one presented in Katz et al. (2023) for ‘remnant leakers’. In short, Katz et al. (2023) propose two modes of Ly C leakage (f|$_{\rm {esc}}$| > 20 per cent), based on galaxies from the sphinx suite of cosmological simulations (Rosdahl et al. 2018, 2022). ‘Bursty leakers’ are galaxies with a recent burst of star formation (within the last 10 Myr, SFR|$_{10}$| > SFR|$_{100}$|⁠), akin to an ionization bounded H ii region with holes. While ‘remnant leakers’ had a strong burst of star formation in the past, but not recently (SFR|$_{100}$| > SFR|$_{10}$|⁠). Making remnant leakers similar to density bounded H ii regions, where the ISM was disrupted enough to halt star formation.

Fig. 13 shows the stacked SFHs of the 51 remnant leaker candidates found in this work, which indicates the presence of a burst in star formation occurring in the past (within the last |$\sim 50$| Myr), but no recent star formation (within the last 10 Myr), in agreement with the SPHINX remnant leakers from Katz et al. (2023). As expected, the photometry for these candidates suggests little-to-no presence of emission lines, which could support a high f|$_{\rm {esc}}$| scenario (Zackrisson et al. 2017). The evidence of these galaxies being in a remnant leaker mode is compelling. An in depth study of these sources will be performed in the future and is beyond the scope of this work. However, given the small amount of galaxies we find to fall into the mini-quenched category (<3 per cent of the total sample), and their low |$\xi _{\rm {ion},0}$|⁠, we conclude that they are not dominant agents in the reionization of the Universe.

6 IMPLICATIONS FOR THE REIONIZATION OF THE UNIVERSE

Recent works based on JWST observations have arrived to the conclusion that there is an overestimation of ionizing photons at the EoR, leading to a so called ‘crisis’ in the ionizing photon budget (see e.g. Trebitsch et al. 2022; Chakraborty & Choudhury 2024; Muñoz et al. 2024). In the following sections, we discuss the implications of our |$\xi _{\rm {ion},0}$| estimations for a stellar mass complete sample on the cosmic ionizing budget.

6.1 Constraints on the cosmic ionizing budget

As in Simmonds et al. (2024a), we now study the implications of our findings on the ionizing cosmic budget through |$\dot{N}_{\rm {ion}}$|⁠, which represents the number of ionizing photons produced per volume unit. We first analyse the contributions of different M|$_{\rm {UV}}$| and redshift bins to |$\dot{N}_{\rm {ion}}$|⁠, by adopting the UV luminosity functions from Bouwens et al. (2021) and some simple prescriptions for f|$_{\rm {esc}}$|⁠: constant of 10 and 20 per cent (Ouchi et al. 2009; Robertson et al. 2013, 2015), and varying according to M|$_{\rm {UV}}$| (Anderson et al. 2017). The latter was constructed using the uniform volume simulation Vulcan, and states the dependence: log(f|$_{\rm {esc}}$|⁠) = |$(0.51\pm 0.4)$|M|$_{\rm {UV}} +7.3 \pm 0.08$|⁠. For |$\xi _{\rm {ion}}$|⁠, we use the |$\xi _{\rm {ion},0}$| relations presented in Fig. 9 for the star-forming sample (i.e. log (SFR|$_{10}$|/SFR|$_{100}$|⁠) |$\ge$| −1). This sample was chosen because it is representative of most galaxies in this study, accounting for >97 per cent of the total sample. We remind the reader that our prospector-inferred ionizing photon production efficiencies assume f|$_{\rm {esc}}$| = 0, such that |$\xi _{\rm {ion},0}$| = |$\xi _{\rm {ion}}$||$\times$|(1−f|$_{\rm {esc}}$|⁠). For simplicity, in the following calculations, we adopt |$\xi _{\rm {ion}}$| = |$\xi _{\rm {ion},0}$| therefore the resulting |$\dot{N}_{\rm {ion}}$| values should be taken as lower limits. For reference, the mean f|$_{\rm {esc}}$| for the star-forming sample is |$\sim 3$|% (see Fig. 12), which would result in a difference in |$\dot{N}_{\rm {ion}}$| of <0.1 dex.

The resulting |$\dot{N}_{\rm {ion}}$| curves as a function of M|$_{\rm {UV}}$| are shown in Fig. 14, each panel displaying a different f|$_{\rm {esc}}$| as indicated in the top left corner. As expected, the curves with fixed f|$_{\rm {esc}}$| are flatter than those derived in Simmonds et al. (2024a), a consequence of now having a stellar mass complete sample that leads to a milder |$\xi _{\rm {ion}}$| evolution with M|$_{\rm {UV}}$|⁠. Once all types of galaxies are potentially included, the importance of faint galaxies in the ionizing cosmic budget is reduced. We note that our exclusion of mini-quenched galaxies does not change this conclusion, since they only account for |$\lesssim 4$|% of the sample. In the case of the variable f|$_{\rm {esc}}$|⁠, the curves are still steep, due to the nature of their f|$_{\rm {esc}}$| prescriptions. On the one hand, Anderson et al. (2017) state that ionizing photons can escape much more easily from faint galaxies. On the other hand, Chisholm et al. (2022) propose that galaxies with bluer UV continuum slopes (⁠|$\beta$|⁠) are likely more efficient Ly C leakers. The |$\beta$| in equation (6) can be expressed as a function of M|$_{\rm {UV}}$| as:

As described in Chisholm et al. (2022), the coefficients a and b connect M|$_{\rm {UV}}$| to the UV continuum slope, |$\beta$|⁠, through the |$\beta$|–M|$_{\rm {UV}}$| relations provided in Bouwens et al. (2014) for selected redshifts. On the top left panel, we also include the results if we adopt a constant log (⁠|$\xi _{\rm {ion}}$|/[Hz erg|$^{-1}$|]) = 25.2 (based on stellar populations, as in Robertson et al. 2013). We find that the shapes of the |$\dot{N}_{\rm {ion}}$| curves with fixed |$\xi _{\rm {ion}}$| are similar to the ones found in this work, albeit slightly offset towards lower |$\dot{N}_{\rm {ion}}$|⁠.

Using the |$\dot{N}_{\rm {ion}}$| curves estimated above, in Fig. 15, we investigate how the ionizing cosmic budget evolves with redshift, given by

where |$\dot{N}_{\rm {ion}}$| is in units of s|$^{-1}$| Mpc|$^{-3}$|⁠, |$\xi _{\rm {ion}}$| is in units of Hz erg|$^{-1}$|⁠, and |$\rho _{\rm {UV}}$| in units of erg s|$^{-1}$| Hz|$^{-1}$| Mpc|$^{-3}$|⁠, and the escape fraction is dimensionless. We assume that the scatter is negligible and that there is an interdependence of variables (e.g. |$\xi _{\rm {ion}}$|–M|$_{\rm {UV}}$| relation). In this figure, we integrate the curves shown in Fig. 14 down to M|$_{\rm {UV}}$||$= -16$| (open triangles), estimated by adopting the luminosity functions from Bouwens et al. (2021). We remind the reader that our sample is UV complete down to M|$_{\rm {UV}} \sim -16$|⁠, and thus, our derived relations are valid in the integrated range. We also include the luminosity densities from Sun & Furlanetto (2016), who fit a power law to the low-mass end. To estimate |$\xi _{\rm {ion}}$| we use the best-fitting relation to our data, as shown in Fig. 9. For consistency, we use only the star-forming galaxies (accounting for |$\sim 94$|% of the total sample).

When adopting the Sun & Furlanetto (2016) luminosity functions, we only assume fixed escape fractions of 10 and 20 per cent (filled circles), motivated by the canonical average f|$_{\rm {esc}}$| needed for galaxies to ionize the Universe (Ouchi et al. 2009; Robertson et al. 2013, 2015). As comparison, we add curves indicating the |$\dot{N}_{\rm {ion}}$| needed to maintain ionization of hydrogen, according to the models of Madau et al. (1999), for clumping factors of 1, 3, and 10 (although see So et al. 2014, for a discussion of the validity of this approach at the EoR). A clumping factor of unity represents a uniform IGM, while larger clumping factors imply a higher number of recombinations taking place in the IGM, and thus, more ionizing photons need to be emitted in order to sustain ionization. Finally, the values found in Simmonds et al. (2024a) for ELGs, are shown as a shaded hatched area. These were provided as upper limits, since they were estimated based on one important assumption: that the sample of ELGs was representative of the entire galaxy population. Indeed, in a recent work by Muñoz et al. (2024), an excess of ionizing photons inferred from JWST observations is identified in the cosmic ionizing budget. With our updated stellar mass complete sample, we are now able to provide more realistic results. We find that the estimations made with our star-forming sample are consistent with those from literature (e.g. Bouwens et al. 2015; Mason et al. 2015; Mason et al. 2019; Naidu et al. 2020; Rinaldi et al. 2024).

In addition to the measurements described above, we include those found through the stellar mass function (instead of |$\rho _{\rm {UV}}$|⁠). In particular, we use results from Weibel et al. (2024), which are particularly relevant to this work since they were constructed from NIRCam observations of galaxies at |$z\sim 4\!-\!9$|⁠. In this scenario, |$\dot{N}_{\rm {ion}}$| can be rewritten as:

where |$\Phi$| is in units of Mpc|$^{-3}$|⁠, f|$_{\rm {esc}}$| is dimensionless and depends on stellar mass, SFR, and size of the galaxy (through the SFR surface density, |$\Sigma _{\rm {SFR}}$|⁠), and |$\dot{n}_{\rm {ion}}$| is in units of s|$^{-1}$|⁠. In this case, we adopt the f|$_{\rm {esc}}$| prescription from Naidu et al. (2020) due to its dependence on |$\Sigma _{\rm {SFR}}$|⁠, such that:

In order to obtain |$\Sigma _{\rm {SFR}}$|⁠, we first follow the size–mass relation and coefficients (for the full sample) presented in table 3 of Morishita et al. (2024). We then combine the derived sizes with SFR|$_{10}$| (averaged over the past 10 Myr), such that |$\Sigma _{\rm {SFR}} = \frac{\text{SFR}/2}{\pi \text{R}^2}$| (Shibuya et al. 2019). After estimating |$\Sigma _{\rm {SFR}}$|⁠, we compute f|$_{\rm {esc}}$||$\times$||$\dot{n}_{\rm {ion}}$| for our stellar mass-complete sample, for each redshift bin. By combining the best-fitting relations with |$\Phi$|⁠, and integrating them down to log(M|$_{\star }$|/[M|${_\odot }$|]) = 7.5, we finally find |$\dot{N}_{\rm {ion}}$| as a function of redshift. The values are shown as stars in Fig. 15, and are in general agreement with the other variable f|$_{\rm {esc}}$| prescriptions (i.e. Anderson et al. 2017; Chisholm et al. 2022).

Although the values of |$\dot{N}_{\rm {ion}}$| are highly uncertain in the early Universe, strong constraints on |$\dot{N}_{\rm {ion}}$| have been placed by observations of the Ly |$\alpha$| forest at lower redshifts (e.g. Becker & Bolton 2013; Gaikwad et al. 2023; Davies et al. 2024, included in Fig. 15). The Ly |$\alpha$| forest refers to a series of absorption lines at wavelengths redder than Ly |$\alpha$| (⁠|$\lambda _{\rm {rest-frame}} \sim 1216$| Å) observed in the spectra of high redshift quasars, produced by the intervening neutral IGM. Interestingly, |$\dot{N}_{\rm {ion}}$| has been observed to flatten once reionization has been completed (with log(⁠|$\dot{N}_{\rm {ion}}$|/[s|$^{-1}$| Mpc|$^{-3}])\sim 50.8$| at |$z \lesssim 5\!-\!6$|⁠). We find that, out of all the prescriptions adopted in this work, the f|$_{\rm {esc}}$| relations from Chisholm et al. (2022) can match the shape of the Ly|$\alpha$| forest constraints. We remind the reader that in order to construct a stellar mass-complete sample, we have ignored all galaxies below the completeness limit (log(M|$_{\star }$|/[M|$_{\odot }])\sim 7.5$|⁠). Therefore, our points represent lower limits to the cosmic ionizing budget.

6.2 Which galaxies reionized the Universe?

In order to determine which galaxies dominate the budget of reionization, we use the same f|$_{\rm {esc}}$| prescriptions as before, but integrate their contributions in three M|$_{\rm {UV}}$| bins: |$-24 \le \text{M}_{\rm {UV}} \lt -20$|⁠, |$-20 \le \text{M}_{\rm {UV}} \le -18$|⁠, and |$-18 \lt \text{M}_{\rm {UV}} \le -16$|⁠. Fig. 6 shows that M|$_{\rm {UV}}$| and stellar mass are correlated (albeit with a large scatter), where brighter galaxies are more massive (and vice-versa), such that each UV luminosity bin also loosely describes a mass bin in our sample. In the four leftmost columns of Fig. 16, we adopt these M|$_{\rm {UV}}$| luminosity bins and show the relative contribution these bins have in |$\dot{N}_{\rm {ion}}$| for each redshift bin above 5 (where observational studies agree the EoR has ended, e.g. Keating et al. 2020; Yang et al. 2020; Zhu et al. 2024). In the rightmost column, we show the contributions of different stellar mass bins to the ionizing budget, by adopting the stellar mass functions of Weibel et al. (2024) and the f|$_{\rm {esc}}$| prescription from Naidu et al. (2020) (stars in Fig. 15). Each row has been normalized to a different clumping factor using the models from (Madau et al. 1999). We find that if the IGM is uniform (C = 1), then galaxies produce enough ionizing radiation in order to sustain hydrogen ionization, independent of the f|$_{\rm {esc}}$| prescription. As the clumping factor increases, it becomes more difficult for galaxies to ionize the Universe by redshift 5. However, most importantly, we discover that for every clumping factor, f|$_{\rm {esc}}$| assumption, and redshift bin, the fainter galaxies (with M|$_{\rm {UV}} \gt = -20$|⁠), and galaxies with low and intermediate stellar masses (log (M|$_{\star }$|/[M|$_{\odot }$|]) < 9.5) dominate the cosmic ionizing budget. In agreement with the results presented in Seeyave et al. (2023), based the First Light and Reionization Epoch Simulations (FLARES; Lovell et al. 2021; Vijayan et al. 2021), and the conclusions reached by forward modelling JWST analogues from the SPHINX simulation, presented in Choustikov et al. (2024).

Therefore, in this work, we confirm that faint low-mass galaxies with bursty star formation have in general enhanced |$\xi _{\rm {ion}}$| compared to massive galaxies and/or galaxies without recent star formation, in agreement with Simmonds et al. (2024a). However, when taking into account the full galaxy population, their contribution is less extreme as might have been thought previously, resolving any potential crisis in the ionizing photon budget of the Universe (e.g. Trebitsch et al. 2022; Chakraborty & Choudhury 2024; Pahl et al. 2024). Finally, by adopting the relations from Chisholm et al. (2022), we are able to reconcile our results with the constraints provided by observations of the Ly |$\alpha$| forest, especially in the flattening of |$\dot{N}_{\rm {ion}}$| that has been observed at lower redshifts (⁠|$z \lesssim 6$|⁠). In summary, we have shown that galaxies produce enough ionizing photons to ionize the Universe by |$z \approx 5\!-\!6$|⁠, without creating a nonphysical excess of ionizing photons in the cosmic budget.

7 CAVEATS AND LIMITATIONS

There are a few important limitations to our method, which we now describe. First, since our method relies on SED fitting, and includes the assumption of certain stellar populations, we are incapable of detecting and appropriately fitting extreme populations. In particular, we assumed a Chabrier IMF, with a maximum stellar mass of 100 M|$_{\odot }$|⁠. This choice was motivated by the size of the sample, in an effort to chose a representative IMF and stellar populations. However, if there are extreme objects in our sample, we might be missing them. For example, Cameron et al. (2024) have found tentative evidence for a top heavy IMF at |$z \sim 6$| although, Tacchella et al. (2024) find an alternative explanation that does not require to invoke exotic stellar populations. If the IMF evolves with redshift, then our choice of a Chabrier IMF would affect our results (including the derived stellar masses). We circumvent this issue to some extent by selecting only the galaxies that were fit with a reduced |$\chi ^2 \lt 1$|⁠, but we note that by doing so, might have lost information on sources that cannot be reproduced by our models.

In the same vein, our modelling does not include a prescription for AGN. JWST has recently unveiled a hidden population of AGN at high redshifts (Juodžbalis et al. 2023; Madau et al. 2024; Maiolino et al. 2024b; Übler et al. 2024). Underestimating the AGN contribution can lead to a systematic overestimation of stellar mass and SFRs by SED fitting codes (Buchner et al. 2024). Unfortunately, accurately identifying which galaxies in our sample host AGN is not trivial, as demonstrated by Wasleske & Baldassare (2024), who compare different techniques used to select AGN in dwarf galaxies (M|$_* \le 10^{9.5}$|M|$_{\odot }$|⁠). They find that any single diagnostic can retrieve at most half of the AGN sample, and most importantly for this work, the AGN identification is least effective when considering photometry alone. Therefore, quantifying the contribution of AGN to our sample is far from the scope of this work. As such, we report our results and caution that our sample might contain some AGNs that can be mimicked by stellar emission.

Another important point is that our work relies heavily on photometric redshift measurements. Specifically, our prospector-inferred redshifts use eazy redshifts as priors. The latter has been proven to provide accurate results for large samples using JADES NIRCam photometry in GOODS-S (Rieke et al. 2023). Moreover, in this work we compare them (when possible) to spectroscopic redshifts compiled from literature. We find a good – but not perfect – agreement, suggesting that our results are in general accurate but likely contain a few outliers where the |$z_{\rm {phot}}$| is unreliable. We stress that these outliers would be of the order of at most a few per cent. As with the stellar populations, our selection of galaxies (fitted with a reduced |$\chi ^2 \le 1$|⁠) ensures that the photometry is well represented by the best-fitting models. Moreover, the addition of HST observations aids in constraining the photometric redshifts. Finally, we do not expect cosmic variance to play an important factor in our results, but plan in the future to perform a similar study in GOODS-N.

In summary, whereas there are intrinsic limitations to our methods, our results are – as much as possible – accurate and representative of the galaxy population in GOODS-S (at |$3 \le z \le 9$|⁠).

8 CONCLUSIONS

We use JWST NIRCam photometry to build a sample of 14 652 galaxies at |$3 \le z \le 9$|⁠, 1640 of them with spectroscopic redshifts from literature. We infer their properties using the SED fitting code prospector, finding two distinct populations of galaxies which can be separated by their burstiness (delimited by log (SFR|$_{10}$|/SFR|$_{100}$|⁠)  = −1). We call these populations star forming and mini-quenched, and note that the mini-quenched galaxies only account for |$\lt 3~{{\ \rm per\ cent}}$| of the total sample. Within the mini-quenched population, we find an interesting subsample with tentative evidence of Ly C leakage (through the UV continuum slope |$\beta$|⁠). These galaxies populate a similar parameter space as the remnant leakers from Katz et al. (2023). Future spectroscopic follow-ups will be necessary to confirm or refute this hypothesis. Our main findings can be summarized as follows.

We find that |$\xi _{\rm {ion},0}$| increases for fainter galaxies with burstier SFHs, in agreement with previous studies, albeit with a milder evolution with redshift. The latter is explained by the nature of our sample, and that previous studies were biased towards galaxies with strong emission lines and/or LAEs. The evolution of |$\xi _{\rm {ion},0}$| with z for the more representative star-forming sample is:

The 2-dimensional fit that accounts for the change of |$\xi _{\rm {ion}}$| with M|$_{\rm {UV}}$| and redshift, for the same sample is given by:

To study the contribution of the galaxies in this study to reionization, we convolve the star-forming relations (which represent |$\gt 97~{{\ \rm per\ cent}}$| of the total sample), with luminosity functions from literature. We find that galaxies, which are detected with JWST, can ionize the Universe by the end of the EoR, if we assume the AGN contribution is minor. In particular, assuming a fixed escape fraction, we find that galaxies fainter than M|$_{\rm {UV}} = -20$| contribute similar amounts of ionizing photons (see Fig. 14), and that galaxies in the range of M|$_{\rm {UV}} = -20$| to |$-16$| dominate the budget of reionization at every redshift bin studied in this work (see Fig. 16). With our stellar mass complete sample, our predictions do not overestimate |$\dot{N}_{\rm {ion}}$| for galaxies at |$z \gt 8$| (see Muñoz et al. 2024). We note that if we extrapolate our trends to fainter magnitudes (to M|$_{\rm {UV}}$| of −14 or -−12), the Universe can be reionized with lower escape fractions. Promisingly, by adopting the relation of f|$_{\rm {esc}}$| with M|$_{\rm {UV}}$| presented in Chisholm et al. (2022), we can conciliate our results regarding the ionizing cosmic budget with the constraints obtained through observations of the Ly |$\alpha$| forest.

ACKNOWLEDGEMENTS

The JADES Collaboration thanks the Instrument Development Teams and the instrument teams at the European Space Agency and the Space Telescope Science Institute for the support that made this program possible. We also thank our program coordinators at STScI for their help in planning complicated parallel observations. We would like to acknowledge the FRESCO team, lead by Pascal Oesch, for developing their observing program with a zero-exclusive-access period.

CS thanks Sarah Bosman, Frederick Davies, James Leftley, and Arjen van der Wel for the insightful discussions. CS, WB, RM, and JW acknowledge support by the Science and Technology Facilities Council (STFC) and by the ERC through Advanced grant number 695671 ‘QUENCH’, and by the UKRI Frontier Research grant RISEandFALL. RM also acknowledges funding from a research professorship from the Royal Society. ST acknowledges support by the Royal Society Research grant G125142. BER, BDJ, PAC & YZ acknowledge support from the NIRCam Science Team contract to the University of Arizona, NAS5-02015. BER also acknowledges support by the JWST Program 3215. DP acknowledges support by the Huo Family Foundation through a P.C. Ho PhD Studentship. This research is supported in part by the Australian Research Council Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), through project number CE170100013. AJB, JC, and GCJ acknowledge funding from the ‘FirstGalaxies’ Advanced Grant from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement no. 789056). SC acknowledges support by European Union’s HE ERC Starting grant no. 101040227 – WINGS. ECL acknowledges support of an STFC Webb Fellowship (ST/W001438/1). IL is supported by the National Science Foundation Graduate Research Fellowship under grant no. 2137424. HÜ gratefully acknowledges support by the Isaac Newton Trust and by the Kavli Foundation through a Newton–Kavli Junior Fellowship. NCV acknowledges support from the Charles and Julia Henry Fund through the Henry Fellowship. The research of CCW is supported by NOIRLab, which is managed by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation.

DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the corresponding author.

Footnotes

REFERENCES

APPENDIX A: GOODNESS OF FITS

Here, we present a comparison between the modelled and observed photometry, defined as the difference between them, divided by the error in the modelled points. The median and errors are shown in Fig. A1, where we exclude the HST bands F435W, F606W, F775W, F814W, and F850LP. The observations in these filters are highly uncertain for our sample, yielding |$\chi$| values well outside of the bounds of the figures. |$\chi$| scatters around zero, with error bars that are mostly symmetric. As expected, the deep NIRCam photometric set is, in general, better represented by the prospector best-fitting models. Fig. A2 shows the same comparison but at rest-frame wavelengths (where the prospector photometric redshift has been adopted). The measurements have been binned into wavelength bins of width 100 Å, and the median values and errors are shown as white circles and error bars. The latter scatter around zero, and most importantly, do not show indication of emission line fluxes being either over or underestimated by prospector. The photometric offsets for each band, defined as the ratio between the observed and modelled fluxes, are given in Table A1.

APPENDIX B: SHAPE OF REDSHIFT AND STELLAR MASS DISTRIBUTIONS

Fig. B1 shows the mean shape of the main properties used to estimate the stellar mass completeness in Section 4.2. The redshift and stellar masses have been divided by the 50|$^{\rm {th}}$| percentile value for each galaxy. We find no significant signs of asymmetry in the posteriors for our sample, and the values are in general well constrained by the fitting routine.

APPENDIX C: INFORMATION GAINED AFTER SED-FITTING

In order to measure how well prospector can constrain stellar masses and ionizing photon production efficiencies in our sample, we use the Kullback–Leibler definition of information gain (IG), given by:

where p represents the parameter of choice: either |$\log \text{M}_{\star }$| or |$\log \xi _{\rm {ion}}$|⁠. IG is a way of measuring the difference between the prior and posterior, in units of bits: an IG  = 0 means the prior is equal to the posterior and no information was gained, whereas a higher value of IG indicates a larger amount of information was gained in the fitting routine. Following the criteria from Simmonds et al. (2018):

• IG < 1: little-to-no information was gained

• 1 |$\le$| IG |$\le$| 2: some information was gained

• IG > 2: parameter is constrained

Figs C1 and C2 show the IG for stellar mass and |$\xi _{\rm {ion},0}$|⁠, respectively, for our stellar mass complete sample. They are shown as a function of the flux in the F444W filter, colour-coded by stellar mass (left) and by burstiness (right). As expected, IG (stellar mass) is highest for brighter massive galaxies, and is lowest for fainter lower-mass galaxies. Importantly, all values are above 1 (the majority lie above 2) and the information gained in stellar mass does not depend on the existence of recent star formation, or strong emission lines. In the case of |$\xi _{\rm {ion},0}$|⁠, we find a trend of increasing IG with burstiness, but emphasize that all measurements are highly constrained (with IG |$\gg 2$|⁠). Promisingly, most mini-quenched galaxies discussed in this work have a large IG. This is due to the richness of the photometric data set, which allows to constrain SED shapes even in cases with no obvious emission lines.

APPENDIX D: BEST-FIT PARAMETERS FOR RELATIONS OF IONIZING PROPERTIES WITH OBSERVED UV MAGNITUDE

For readability, the relations of |$\xi _{\rm {ion},0}$| and |$\dot{n}_{\rm {ion}}$|⁠, as a function of M|$_{\rm {UV}}$|⁠, presented in Figs 9 and 10, are given in Table D1 and D2, respectively.