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 inlineformulatexmath idTM0001 notationLaTeX3 le zrm phot le 9texmathinlineformula 90 per cent complete in stellar mass down to logMinlineformulatexmath idTM0002 notationLaTeXstar texmathinlineformulaMinlineformulatexmath idTM0003 notationLaTeXodot approx 7.5texmathinlineformula. 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 inlineformulatexmath idTM0004 notationLaTeXxi rm iontexmathinlineformula with redshift and UV magnitude than previously reported. Interestingly we observe two distinct populations in inlineformulatexmath idTM0005 notationLaTeXxi rm iontexmathinlineformula distinguished by their burstiness given by SFRinlineformulatexmath idTM0006 notationLaTeX10texmathinlineformulaSFRinlineformulatexmath idTM0007 notationLaTeX100texmathinlineformula. Both populations show the same evolution with z and Minlineformulatexmath idTM0009 notationLaTeXrm UVtexmathinlineformula but have a different inlineformulatexmath idTM0010 notationLaTeXxi rm iontexmathinlineformula normalization. We convolve the more representative inlineformulatexmath idTM0011 notationLaTeXlog xi rm ion ztextMrm UVtexmathinlineformula relations accounting for inlineformulatexmath idTM0012 notationLaTeXsim 97texmathinlineformula  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 inlineformulatexmath idTM0013 notationLaTeXalphatexmathinlineformula forest at inlineformulatexmath idTM0014 notationLaTeXzlesssim 6texmathinlineformula. We conclude that galaxies with Minlineformulatexmath idTM0015 notationLaTeXrm UVtexmathinlineformula between inlineformulatexmath idTM0016 notationLaTeX16texmathinlineformula and inlineformulatexmath idTM0017 notationLaTeX20texmathinlineformula adopting a reasonable escape fraction can produce enough ionizing photons to ionize the Universe without exceeding the required ionizing photon budget.
['Astrophysics - Astrophysics of Galaxies']
Ionizing properties of galaxies in JADES for a stellar mass complete sample resolving the cosmic ionizing photon budget crisis at the Epoch of Reionization
{'Ionising properties of galaxies in JADES for a stellar mass complete sample: resolving the cosmic ionising photon budget crisis at the Epoch of Reionisation': "- C. Simmonds, 1 , 2 ★ S. Tacchella, 1 , 2 K. Hainline, 3 B. D. Johnson, 4 D. Puskás, 1 , 2 B. Robertson, 5\n- W. M. Baker, 1 , 2 R. Bhatawdekar, 6 K. Boyett, 7 , 8 A. J. Bunker, 9 P. A. Cargile, 4 S. Carniani, 10\n- J. Chevallard, 9 M. Curti, 11 E. Curtis-Lake, 12 Z. Ji, 3 G. C. Jones, 9 N. Kumari, 13 I. Laseter, 14\n- R. Maiolino, 1 , 2 , 15 M. V. Maseda, 14 P. Rinaldi, 3 A. Stoffers, 1 , 2 H. Übler, 1 , 2 N. C. Villanueva, 1 , 2\n- C. C. Williams, 16 C. Willott, 17 J. Witstok, 1 , 2 and Y. Zhu 3\n- 1 The Kavli Institute for Cosmology (KICC), University of Cambridge, Madingley Road, Cambridge, CB3 0HA, UK\n- 2 Cavendish Laboratory, University of Cambridge, 19 JJ Thomson Avenue, Cambridge, CB3 0HE, UK\n- 3 Steward Observatory, University of Arizona, 933 N. Cherry Avenue, Tucson, AZ 85721, USA\n- 4 Center for Astrophysics | Harvard & Smithsonian, 60 Garden St., Cambridge, MA 02138, USA\n- 5 Department of Astronomy and Astrophysics University of California, Santa Cruz, 1156 High Street, Santa Cruz CA 96054, USA\n- 6 European Space Agency (ESA), European Space Astronomy Centre (ESAC), Camino Bajo del Castillo s/n, 28692 Villanueva de la Cañada, Madrid, Spain 7 School of Physics, University of Melbourne, Parkville 3010, VIC, Australia\n- 8 ARC Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), Canberra 2611, Australia\n- 9 Department of Physics, University of Oxford, Denys Wilkinson Building, Keble Road, Oxford OX1 3RH, UK\n- 10 Scuola Normale Superiore, Piazza dei Cavalieri 7, I-56126 Pisa, Italy\n- 11 European Southern Observatory, Karl-Schwarzschild-Strasse 2, 85748 Garching, Germany\n- 12 Centre for Astrophysics Research, Department of Physics, Astronomy and Mathematics, University of Hertfordshire, Hatfield AL10 9AB, UK\n- 13 AURA for European Space Agency, Space Telescope Science Institute, 3700 San Martin Drive. Baltimore, MD, 21210\n- 14 Department of Astronomy, University of Wisconsin-Madison, 475 N. Charter St., Madison, WI 53706 USA\n- 15 Department of Physics and Astronomy, University College London, Gower Street, London WC1E 6BT, UK\n- 16 NSF's National Optical-Infrared Astronomy Research Laboratory, 950 North Cherry Avenue, Tucson, AZ 85719, USA\n- 17 NRC Herzberg, 5071 West Saanich Rd, Victoria, BC V9E 2E7, Canada  \nAccepted XXX. Received YYY; in original form ZZZ", 'ABSTRACT': 'We use NIRCam imaging from the JWST Advanced Deep Extragalactic Survey (JADES) to study the ionising properties of a sample of 14652 galaxies at 3 ≤ 𝑧 phot ≤ 9, 90% complete in stellar mass down to log(M ★ /[M ⊙ ]) ≈ 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 ionising photon production efficiency ( 𝜉 ion ) with redshift and UV magnitude than previously reported. Interestingly, we observe two distinct populations in 𝜉 ion , distinguished by their burstiness (given by SFR 10 /SFR 100 ). Both populations show the same evolution with 𝑧 and M UV , but have a different 𝜉 ion normalisation. We convolve the more representative log ( 𝜉 ion ( 𝑧, M UV )) relations (accounting for ∼ 97% of the sample), with luminosity functions from literature, to place constraints on the cosmic ionising photon budget. By combining our results, we find that one of our models can match the observational constraints from the Ly 𝛼 forest at 𝑧 ≲ 6. We conclude that galaxies with M UV between -16 and -20, adopting a reasonable escape fraction, can produce enough ionising photons to ionise the Universe, without exceeding the required ionising photon budget.  \nKey words: Galaxies: high-redshift - Galaxies: evolution - Galaxies: general - Cosmology: dark ages, reionization, first stars', '1 INTRODUCTION': "The Epoch of Reionisation (EoR) is one of the major phase transitions of the Universe, when it went from being dark and neutral to highly ionised, allowing Lyman Continuum (LyC; with energies  \nabove 13.6 eV) radiation to travel through the intergalactic medium (IGM). Observations place the end of this epoch at 𝑧 ∼ 6 (Becker et al. 2001; Fan et al. 2006; Yang et al. 2020), with some favouring a later reionisation at 𝑧 ∼ 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 ionising the Universe. The community widely agrees that young massive stars within galaxies are key, since  \nthey produce copious amounts of ionising photons, which might be able to escape the interstellar medium (ISM), and eventually ionise the IGM (Hassan et al. 2018; Rosdahl et al. 2018; Trebitsch et al. 2020). However, the nature of the galaxies that drive reionisation: 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 (AGN) contribute to reionisation (Dayal et al. 2020; Maiolino et al. 2023; Madau et al. 2024).  \nThe stellar mass of galaxies has been seen to correlate with how efficiently ionising photons are produced (Simmonds et al. 2024a). Simulations indicate that it also relates how these ionising photons escape (Paardekooper et al. 2015). The latter is measured through their Lyman Continuum escape fractions (f esc (LyC)), defined as the ratio between the H-ionising radiation that is emitted intrinsically, and that which reaches the IGM. In order for galaxies to account for the reionisation of the Universe, either a significant average f esc value is required (f esc =10-20%; Ouchi et al. 2009; Robertson et al. 2013, 2015; Finkelstein et al. 2019; Naidu et al. 2020), or a high ionising photon production efficiency. These ranges of f esc have been observed for individual star-forming galaxies at 𝑧 ≲ 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 ionising photon production efficiency ( 𝜉 ion ), given by the ratio between the rate of ionising photons being emitted ( / 𝑛 ion ), and the monochromatic non-ionising ultra-violet (UV) luminosity density. Indeed, observational studies up to 𝑧 ∼ 9 have found 𝜉 ion to increase as a function of redshift (e.g. Endsley et al. 2021; Stefanon et al. 2022; Tang et al. 2023; Atek et al. 2023; Simmonds et al. 2023, 2024a; Harshan et al. 2024; Pahl et al. 2024; Saxena et al. 2024).  \nThe behaviour of 𝜉 ion as a function of redshift has important consequences on the cosmic budget of reionisation (e.g. Muñoz et al. 2024), defined as the number of ionising photons produced per comoving volume unit of the Universe ( / 𝑁 ion ). Three ingredients must be provided in order to study / 𝑁 ion : (1) a prescription for f esc , (2) a UV luminosity density function, 𝜌 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 prep.), and (3) constraints on 𝜉 ion . In addition, the IGM clumping factor of the Universe has to be considered (e.g. Madau et al. 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 reionisation since it relates to the amount of atomic recombinations taking place in the IGM. Briefly, a higher clumping factor implies that more ionising photons need to be emitted per unit volume at a given redshift in order to sustain hydrogen ionisation.  \nBefore the launch of the James Webb Space Telescope (JWST; Gardner et al. 2023), it was common practice to set f esc and 𝜉 ion as constants. Fortunately, the JWST has given us an unprecedented view of the early Universe at restframe optical wavelengths, which has allowed us to place better constraints on 𝜉 ion . In particular, Simmonds et al. (2024a) studied a sample of emission-line galaxies (ELGs) at 𝑧 ∼ 4 -9 using photometry obtained with the Near Infrared Camera (NIRCam; Rieke et al. 2023), on board the JWST. The mean redshift of reionisation is at 𝑧 = 7 . 68 ± 0 . 79 (Planck Collaboration et al. 2020), meaning that this study (and the many others enabled by JWST) probe deep into the EoR. Through H 𝛼 and [O /i.pc/i.pc/i.pc] emission line fluxes, Simmonds et al. (2024a) estimated 𝜉 ion for a sample of 677 galaxies. In parallel, they inferred the same quantity by using the spectral energy distribution (SED) fitting code Prospector  \n(Johnson et al. 2019, 2021). They find that the 𝜉 ion measurements estimated by emission line fluxes agree with the values obtained by Prospector . Additionally, they conclude that 𝜉 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 et al. 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).  \nThe 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 𝛼 and/or [O /i.pc/i.pc/i.pc] emission. In fact, Laseter et al. (in prep.) demonstrate 𝜉 ion is consistently high (log( 𝜉 ion /[Hz erg -1 ]) ≈ 25 . 5) down to [O /i.pc/i.pc/i.pc] equivalent widths of 200Å due to low metallicities ( 𝑍 ≲ 1 / 10 𝑍 ⊙ ), further demonstrating this bias from the past results. Given the agreement found between Prospector and the emission line measurements of 𝜉 ion , in this work, we use Prospector to fit the full JWST Advanced Deep Extragalactic Survey (JADES; Eisenstein et al. 2023a; Bunker et al. 2023) 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% complete in stellar mass down to masses of log(M ★ /[M ⊙ ]) ∼ 7 . 5, providing us a deep statistical view of the ionising properties of galaxies.  \nThe structure of this paper is the following. In § 2 we present the data used in this work, along with the sample selection criteria. In § 3 we present our Prospector fitting method. Some general properties of the sample are given in § 4, followed by our constraints on the ionising properties of galaxies in § 5. Implications for reionisation are discussed in § 6, while the caveats and limitations of our methods are discussed in § 7. Finally, brief conclusions are presented in § 8. Throughout this work we assume Ω 0 = 0 . 315 and 𝐻 0 = 67 . 4 km s -1 Mpc -1 , following Planck Collaboration et al. (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 (Eisenstein et al. 2023a; Bunker et al. 2023) 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 ∼ 45 arcmin 2 with an average exposure time of 130 hours. When available, we also use photometry from the JWST Extragalactic Medium-band Survey (JEMS; Williams et al. 2023), and from the First Reionisation Epoch Spectroscopic Complete Survey (FRESCO; Oesch et al. 2023).  \n<!-- image -->  \nFigure 1. Top panel/colon.up Throughputs of the HST and NIRCam filters used in the SED fitting in this work for galaxies at 3 ≤ 𝑧 phot ≤ 9. It is important to note that not all of these filters are available for every galaxy studied in this work. The hatched areas show the HST ACS bands: F435W, F606W, F775W, F814W, F850LP, and the HST IR bands: F105W, F125W, F140W, and F160W. While the filled regions are labelled and show the JWST NIRCam bands, from left to right: F070W, F090W, F115W, F150W, F162M, F182M, F200W, F210M, F250M, F277W, F300M, F335M, F356W, F410M, F430M, F444W, F460M, and F480M. The lines show the observed wavelengths of selected emission lines (Ly 𝛼 , [O /i.pc/i.pc], [O /i.pc/i.pc/i.pc], and H 𝛼 ) with redshift. Bottom panel/colon.up percentage of sources in our final, stellar mass-complete sample that are covered by each JWST NIRCam filter. The medium bands are shown as dotted areas.  \n<!-- image -->', '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 mediumbandandJADESNIRCambroadandmediumbandimaging. Detection was performed using the photutils (Bradley et al. 2022) software package, identifying sources with contiguous regions of the SNRmosaicwithsignal > 3 𝜎 and five or more contiguous pixels. We also use photutils to perform circular aperture photometry with filter-dependent aperture corrections based on model point-spreadfunctions following the method of Ji et al. (2023), 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 NearInfrared Deep Extragalactic Legacy Survey (CANDELS; Grogin et al. 2011; Koekemoer et al. 2011), as well as reductions of GOODSS from Giavalisco et al. (2004), are used. The details of the catalogue generation and photometry will be presented in Robertson et al., (in prep). 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% 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  \n(see Appendix A). The throughputs of the filters used in this work are shown in Figure 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 ( 𝑧 phot ) can be obtained with great accuracy (Hainline et al. 2023). When fitting SEDs in this work, we use the 𝑧 phot inferred by the template-fitting code EAzY (Brammer et al. 2008), as described in Hainline et al. (2023) and Rieke et al. (2023). Additionally, when available, we make use of spectroscopic redshifts ( 𝑧 spec ) from the Near-Infrared Spectrograph (NIRSpec; Jakobsen et al. 2022), as well as those reported in literature (Puskás et al. in prep), including the JADES NIRSpec redshifts from Bunker et al. (2023) 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. 2022), and a redshift catalogue (F. Sun private communication) produced from grism data of the First Reionisation Epoch Spectroscopic Complete Survey (FRESCO; Oesch et al. 2023).  \nFigure 2. Representative example SEDs and best-fit spectra for galaxies in our sample, assuming a continuity (i.e. non-parametric) SFH. From top to bottom the redshift increases from 𝑧 = 3 . 5 to 8.5, in steps of 1, the galaxy identifiers, along with their 𝑧 phot and reduced 𝜒 2 (for the JWST NIRCam photometry), are shown in the top left corner of each panel. The symbols show the photometric points for HST (open squares), JWST NIRCam (open circles), and model photometry (orange diamonds), respectively. The grey curves show the best-fit spectra obtained by Prospector , with the spectral region used to estimate UV continuum slope ( 𝛽 ; 𝜆 rest -frame = 1250 -2500 Å) shaded in purple. The observed wavelengths of Ly 𝛼 , [O /i.pc/i.pc], [O /i.pc/i.pc/i.pc], and H 𝛼 are shown as vertical dotted lines.  \n<!-- image -->  \nobs  \nWe 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 𝑧 spec measurements, we defined four categories to collate the several redshift measurements into one, which we call 𝑧 best :  \n- · Category A: Only one redshift labelled as having the highest quality. We use this redshift as 𝑧 best .\n- · Category B: Multiple redshifts with the same (highest) quality that agree when rounded to the second decimal. We use this rounded value as 𝑧 best .  \n- · Category C: Multiple solutions with the same (highest) quality that have a non-dramatic disagreement (with a difference smaller than 1). We define 𝑧 best as the mean between the solutions and add errors reflecting the difference.\n- · Category D: Multiple solutions with the same (highest) quality that have a dramatic disagreement ( Δ > 1). In these cases, we keep all the highest quality redshifts and follow up with a visual inspection.  \nz  \nz  \nFigure 3. Comparison of input and output photometric redshifts. The vertical axis shows the photometric redshifts inferred using EAzY . These redshifts were used as priors in the SED fitting with Prospector , the medians of these posteriors are shown in the horizontal axis. The black point shows the redshift median errors, we note that the Prospector -inferred values tend to be better constrained, but that they overall follow a one-to-one trend with those inferred using EAzY , as can be seen by the best-fit line (dashed black line). The points that fall outside of the dotted black lines are considered catastrophic outliers, there are 2036 galaxies that fall into this category, corresponding to ∼ 8%. We note that the outlier fraction is estimated by comparing the median values inferred by EAzY and Prospector . The former has considerably larger error bars in general.  \n<!-- image -->", '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 ≤ 𝑧 ≤ 9. Based on the nature of the redshift used (photometric or spectroscopic), we construct two samples, that we now introduce.  \nThe photometric sample is initially composed of 37272 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 35442 galaxies. The cases that failed ( ∼ 5%) were due to either poor photometric coverage ( < 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 𝜒 2 ≤ 1, resulting in a final photometric sample with 25319 galaxies. The completeness in stellar mass and UV magnitude is discussed in Section 4. We note that if we extend the limit up to 𝜒 2 = 10, we obtain 29323 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.  \nz  \nz  \nFigure 4. Comparison of photometric and spectroscopic redshifts, when spectroscopic redshifts in category A and/or B are available. Both Prospector (larger orange circles) and EAzY (smaller grey circles) inferred values generally follow a one-to-one relation with the spectroscopic redshifts, with some exceptions. We conclude that both codes can retrieve 𝑧 spec successfully in the majority of cases. The points that fall outside of the area delimited by the dotted lines are considered catastrophic outliers. There are 44 (35) of such objects inferred by Prospector ( EAzY ), corresponding to ∼ 3 . 1% ( ∼ 2.5%) of the subsample with spectroscopic redshifts in categories A and B.  \n<!-- image -->', '3 SED FITTING WITH PROSPECTOR': 'As demonstrated in Simmonds et al. (2024a), the ionising 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).  \nProspector 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.  \nFor the photometric redshift sample, we adopt a clipped normal distribution using the EAzY 𝑧 phot as the redshift mean, with the sigma given by the 𝑧 phot errors. Whereas for the spectroscopic redshift sample we fix the redshift to 𝑧 spec for galaxies in Categories A, B, and D, and use the same prior as for the photometric sample for galaxies in Category C.  \nWe vary the dust attenuation and stellar population properties  \nFigure 5. Stellar mass completeness of the sample, divided into three fields depending on exposure time (T exp ). The larger grey circles show the stellar mass of each galaxy as a function of redshift, while the smaller coloured circles show the limiting mass for the faintest 20%. The solid curves denote the 90% completeness at each redshift and for each depth, as indicated in the legend, following the prescription of Pozzetti et al. (2010). We find our sample is stellar mass complete down to log(M ★ /[M ⊙ ]) ≈ 7 . 5, with small variations depending on the depth of the field. The dotted vertical lines mark the redshift bins used in the mass completeness estimation, while the numbers provided in each bin indicate the total number of galaxies in the bin (top row, bold faced), and the galaxies above the limiting mass (bottom row, in brackets).  \n<!-- image -->  \nfollowing Tacchella et al. (2022a). In summary, we use a two component dust model described in Charlot & Fall (2000) and Conroy et al. (2009). This model accounts for the differential effect of dust on young stars ( < 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 ⊙ , respectively, allowing the stellar metallicity to explore a range between 0.01 - 1 Z ⊙ , 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 ionisation parameters (log ⟨ 𝑈 ⟩ 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. 2023b). 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 ( 𝜇 = 1 . 0 , 𝜎 = 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,  \nwhere we show a comparison between the modelled and observed photometry, as a function of the bands being used.  \nFigure 2 shows example SEDs and best-fit spectra for our sample. Fromtop 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 JWSTNIRCam(circles) photometry with their corresponding errors. The purple shaded area corresponds to the rest-frame spectral region at 𝜆 rest -frame = 1250 -2500 Å, used to estimate the rest-frame UV continuum slope ( 𝛽 ; Calzetti et al. 1994), defined as F 𝜆 ∝ 𝜆 𝛽 . We obtain 𝛽 by fitting a line to the best-fit SED provided by Prospector . We stress that there are important limitations to the measurements of 𝛽 without spectra (for a full review, see Austin et al. 2024) . Finally, the vertical dashed lines show the observed wavelength of Ly 𝛼 , [O /i.pc/i.pc], [O /i.pc/i.pc/i.pc], and H 𝛼 . 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 JADESNIRCamphotometry set (see Figure 13 of Rieke et al. 2023). In Figure 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 fit slope of 0 . 96 ± 0 . 01 (black dashed line), demonstrating that both methodsare in general agreement. We note that the number of sources decreases significantly as a function of photometric redshift: 10215 at 3 ≤ 𝑧 < 4, 6005 at 4 ≤ 5, 5203 at 5 < 𝑧 ≤ 6, 3037 at 6 < 𝑧 ≤ 7, 651 at 7 < 𝑧 ≤ 8, and, 208 at 8 < 𝑧 ≤ 9.  \nIn order to test how well the 𝑧 phot retrieve real ( 𝑧 spec ) redshifts, in Figure 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 𝑧 phot estimations. We find a good agreement between 𝑧 phot and 𝑧 spec (i.e. | 𝑧 spec -𝑧 phot | < 0 . 15), with only a small fraction of outliers. Specifically, 44 (35) for Prospector ( EAzY ) derived redshifts, corresponding to ∼ 3 . 1% ( ∼ 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  \nFigure 6. Stellar mass as a function of (observed) M UV for our stellar mass complete sample, divided into redshift bins, and colour-coded by the burstiness of their SFHs (i.e. ratio between star formation in the past 10 Myr and the one averaged over the last 100 Myr). For comparison, we overlay the log(M ★ )-M UV mass-to-light relations from Duncan et al. (2014), Song et al. (2016), and Tacchella et al. (2018). The best-fit relation, using forward modelling to take account for the completeness of the sample, is shown as a purple solid curve and shaded area (see Table 1). The black vertical lines in the bottom right panel show the limits of three M UV bins delimited by M UV = -24 , -20 , -18, and -16, which are used later in Figure 16. As expected, brighter galaxies tend to have higher stellar masses, however, there is a considerable scatter in M UV for a fixed stellar mass.  \n<!-- image -->  \nFigure 7. Compilation of 𝜉 ion , 0 from literature, in comparison to our stellar mass complete sample (this work, pentagons). The best-fit line from Stefanon et al. (2022), with a slope of dlog( 𝜉 ion , 0 )/dz = 0 . 09 ± 0 . 01, is shown as a black dashed line and shaded area. In order of increasing redshift, the markers in grey scale represent: UV-faint galaxies at 𝑧 ∼ 3 -7 from Prieto-Lyon et al. (2023) (squares), Ly 𝛼 emitters from Ning et al. (2023) (triangles) and Simmonds et al. (2023) (stars), the H 𝛼 emitters from Rinaldi et al. (2024) (thin diamonds), and finally, the ELGs from Simmonds et al. (2024a), which have a slope of dlog( 𝜉 ion , 0 )/dz = 0 . 07 ± 0 . 02. The thick diamonds show the stacked results from Kumari et al. (2024), for galaxies above and below 𝑧 = 6 . 3. Due to observational limitations and the emission line methods used to estimate 𝜉 ion , 0 (i.e. through H 𝛼 and/or [O /i.pc/i.pc/i.pc] emission) the previous samples were biased towards star-forming galaxies. Our stellar mass complete sample reveals a population of galaxies with low ionising photon production efficiency (log( 𝜉 ion , 0 /[Hz erg -1 ]) ≲ 24.5), and a considerably less significant evolution with redshift (dlog( 𝜉 ion , 0 )/dz ∼ 0 . 01, see Figure 8).  \n<!-- image -->  \nFigure 8. 𝜉 ion , 0 (top) and / 𝑛 ion (bottom) inferred through SED fitting with Prospector , as a function of redshift, colour-coded by a measure of burstiness (defined here as the ratio between the star formation sustained in the past 10 Myr and the one averaged over the past 100 Myr), for the stellar mass complete sample (log(M ★ /[M ⊙ ]) > 7 . 5). The circles with black edges show the spectroscopic sample, while the circles without edges show the photometric sample. Within the latter, the galaxies with zero recent star formation (SFR 10 = 0 M ⊙ yr -1 ) are shown as crosses. The best-fit relations to both samples are shown in the top left corner of each panel. The best fit to the spectroscopic data (blue-filled line) yields a more positive slope than the one derived for the entire photometric sample (black dashed line). Particularly, the fit to 𝜉 ion , 0 in the spectroscopic sample (dlog( 𝜉 ion , 0 )/dz ∼ 0 . 05 ± 0 . 02) is closer to the 𝜉 ion , 0 versus redshift relations from the literature (within errors of the findings of Simmonds et al. 2024a). This is because the spectroscopic sample likely suffers from the same biases as previous studies (i.e. biased towards star-forming galaxies with detectable emission lines). Finally, the pink squares with error bars show the best-fit relations per redshift bin shown in Figures 9 and 10, for a fixed M UV of -18. When the full stellar mass complete sample is considered, both 𝜉 ion , 0 and / 𝑛 ion show only a slight evolution with redshift, 𝜉 ion , 0 shows a strong correlation with burstiness, while / 𝑛 ion does not.  \n<!-- image -->  \nredshift  \nthat the SED modelling uncertainties and redshift variations are selfconsistent.  \nGiven the agreement between 𝑧 phot and 𝑧 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.  \nTo assess the 90% stellar mass completeness limit of our photometric sample, we follow the procedure described in Section 5.2 of  \nPozzetti et al. (2010). In summary, for every redshift we define a minimum mass (M min ), above which all types of galaxies can potentially be observed. To obtain M min , we first need to calculate the limiting stellar mass (M lim ) for each galaxy, given by:  \nlog ( M lim ) = log ( M ★ ) + 0 . 4 ( m -m lim ) , (1)  \nwhere M ★ is stellar mass in units of solar masses. M lim represents the mass a galaxy would have if its apparent magnitude (m) were equal to the limiting magnitude of the survey (m 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 exp < 25 ks), deep (25 ks ≤ T exp < 65 ks), and ultra-deep (T exp > 65 ks), with 5 𝜎 flux depths of 6 nJy, 4.5 nJy, and 2.65 nJy, respectively. Once M lim has been calculated for every galaxy, we compute M min for  \nFigure 9. Dependence of 𝜉 ion , 0 on UV magnitude, separated in redshift bins, and colour-coded by the burstiness of their SFHs. The individual error bars are shown in grey. The top left of each panel shows the redshift range and the amount of galaxies it contains. The number in parenthesis corresponds to a small secondary population of galaxies with log(SFR 10 /SFR 100 ) < -1, that lie systematically below the general 𝜉 ion , 0 trends with M UV and redshift. If these two populations of galaxies are fit separately (when possible), the slope of 𝜉 ion , 0 with M UV is consistent between them, however, their 𝜉 ion , 0 intercept is different. The values for both fits are shown at the bottom of each panel, the darker coloured shaded area and text correspond to the population with no (or very little) recent star formation. We note that this population accounts for only < 3% of the total sample.  \n<!-- image -->  \nFigure 10. Same as Figure 9 but for / 𝑛 ion instead. Contrary to the case of 𝜉 ion , 0 , / 𝑛 ion has a stronger dependence on M UV , where the faintest (brightest) galaxies produce less (more) ionising photons.  \n<!-- image -->  \neach field by selecting the faintest 20% of the galaxies with a given T exp , and finding the limit above which 90% of the selected M lim values lie. Figure 5 shows the stellar mass of our sample as a function of redshift (grey circles), the coloured small circles show the M lim of the faintest 20% of galaxies for each field, and the filled lines denote the 90% completeness as a function of redshift. We find that, in general, our photometric sample is stellar mass complete down to log(M ★ /[M ⊙ ]) ≈ 7.5, with slight variations depending on the depth of the observations (within ∼ 0 . 2 dex).  \nIn 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 ★ /[M ⊙ ]) ≈ 7.5), yielding a sample of 14652 galaxies, which are UV complete down to M UV ≈ -16.', '4.3 M UV -M ★ relation': 'Figure 6 shows the stellar mass of our stellar mass complete sample as a function of M 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 UV -M ★ plane taking into account the completeness limits discussed in the previous section. Specifically, we fit the following relation  \nlog ( M ★ ) = 𝛼 ( M UV + 19 . 5 ) + log ( M ★, 0 ) , (2)  \nwhere M ★ is stellar mass in units of solar masses, 𝛼 is the slope of the M UV -M ★ relation and log ( 𝑀 ★, 0 ) is its normalisation (stellar mass at 𝑀 UV = -19 . 5 AB mag). In addition to these two parameters, we also fit for the scatter 𝜎 M UV -M ★ in M UV at fixed M ★ . 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 UV -M ★ relation (with scatter) with 1 million galaxies assuming an M UV that follows the UV luminosity function (Bouwens et al. 2021), bin this drawn galaxy sample in the M UV -M ★ plane taking into account the completeness limits, and compare the normalised histogram with the observed one. The best-fit parameters can be found in Table 1, for the highest redshift bin (8 < 𝑧 ≤ 9) we adopt the best-fit results from the previous redshift bin (7 < 𝑧 ≤ 8).  \nFrom the nearly linear relation between log ( M ★ ) and log ( SFR ) (i.e. star-forming main sequence; Brinchmann et al. 2004; Daddi et al. 2007; Salim et al. 2007), we would expect a slope 𝛼 of the M UV -M ★ relation close to -0 . 4, assuming a simple, linear conversion between the UV luminosity and SFR (Kennicutt 1998; Madau et al. 1998; Salim et al. 2007). Our steeper slope of ≈ -0 . 6 is consistent with having a higher mass-to-light ratio ( 𝑀 ★ / 𝐿 UV ) for more massive systems, in agreement with the findings of González et al. (2011). We further find that log ( 𝑀 ★, 0 ) decreases with increasing redshift ( 𝑀 ★, 0 ∝ ( 1 + 𝑧 ) -2 . 7 ), i.e. at fixed M UV = -19 . 5 the typical stellar mass of the galaxies is lower at earlier cosmic times. This implies a lower 𝑀 ★ / 𝐿 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 UV -M ★ relation is substantial with 0 . 6 -0 . 7 dex. This scatter includes both the intrinsic scatter of the M UV -M ★ relation and the observational uncertainty, though the latter is only of the order of 0.2 dex.  \nAscomparison,weshowthelog(M ★ )-M 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  \nTable 1. Best-fit parameters for log(M ★ )-M UV relation as parameterised by Eq. 2. Column /one.up/colon.up redshift bin. Column /two.up/colon.up slope of the relation. Column /three.up/colon.up stellar mass normalisation (log(M ★ ) at M UV = -19 . 5. Column /four.up/colon.up scatter in log(M ★ ) at fixed M UV .  \nthe accretion rate of dark matter halos. With our stellar mass complete sample, we find a somewhat steeper slope than previous observational studies (but mostly < 1 𝜎 discrepant; ∼ -0 . 6 instead of ∼ -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 ∼ -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 IONISING PROPERTIES OF GALAXIES': 'In this section, we provide a brief overview of the methods traditionally used to infer ionising properties of galaxies. Followed by the trends of 𝜉 ion and / 𝑛 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 LyC. An excerpt of the properties used in this section is given in Table 2.', '5.1 Measuring 𝜉 ion from emission line fluxes': 'The ionising photon production efficiency, 𝜉 ion , can be measured through H 𝛼 and/or [O /i.pc/i.pc/i.pc] emission line fluxes. In order to use H 𝛼 as a proxy for ionising photon production efficiency ( 𝜉 ion ), one must first correct the H 𝛼 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 LyC f esc , and that dust has a negligible effect on LyC photons, the dust corrected H 𝛼 luminosity relates directly to the rate of ionising photons ( / 𝑛 ion , in units of s -1 ) that are being emitted ( / 𝑛 ion = 7 . 28 × 10 11 L(H 𝛼 ); Osterbrock & Ferland 2006). 𝜉 ion is then the ratio between / 𝑛 ion and the observed monochromatic UV luminosity density (measured at rest-frame 𝜆 = 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 /i.pc/i.pc/i.pc] 𝜆 5007 is available instead, its equivalent width can give a measure of the ionisation 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 gasphase metallicity and nebular conditions.  \nTable 2. Table excerpt of general properties for a selection of galaxies studied in this work. Column /one.up/colon.up JADES identifier, composed of the coordinates of the centroid rounded to the fifth decimal place, in units of degrees. Column /two.up/colon.up photometric redshift inferred using the SED fitting code Prospector . Column /three.up/colon.up logarithm of stellar mass. Column /four.up/colon.up logarithm of the ionising photon production efficiency. Column /five.up/colon.up logarithm of the rate of ionising photons being emitted. Column /six.up/colon.up logarithm of the burstiness of SFH, defined as the ratio between recent star formation (<10 Myr) and the one averaged over the past 100 Myr (using the median values of SFR 10 and SFR 100 ). Column /seven.up/colon.up observed UV magnitude.  \nIn Simmonds et al. (2024a), the ionising properties of a sample of 677 emission line galaxies at 𝑧 ∼ 4 -9 were analysed. These galaxies had signatures of either H 𝛼 and/or [O /i.pc/i.pc/i.pc] emission in their NIRCam photometry, specifically with a 5 𝜎 detection in the filter pairs F335M-F356W and F410M-F444W. Both methods mentioned above were used to estimate 𝜉 ion , while in parallel, Prospector was used to fit their full NIRCam photometry (assuming f 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 ionising 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 ionising photon production efficiencies ( 𝜉 ion , 0 ) include a correcting factor of 1 -f esc , such that:  \n𝜉 ion , 0 = 𝜉 ion × ( 1 -f esc ) . (3)  \nAppendix C shows how well Prospector can constrain stellar masses and ionising photon production efficiencies, as a function of flux in F444W.', '5.2 Trends of ionising properties with redshift': 'Wefirst compare our results to those from the literature, in particular, for the evolution of 𝜉 ion , 0 with redshift. Figure 7 encodes most results to date. Stefanon et al. (2022) compiled 𝜉 ion , 0 measurements up to 𝑧 ∼ 8 (containing data from Stark et al. 2015, 2017; MármolQueraltó et al. 2016; Nakajima et al. 2016; Bouwens et al. 2016; Matthee et al. 2017; Harikane et al. 2018; Shivaei et al. 2018; De Barros et al. 2019; Lam et al. 2019; Faisst et al. 2019; Tang et al. 2019; Nanayakkara et al. 2020; Emami et al. 2020; Endsley et al.  \n2021; Atek et al. 2022; Naidu et al. 2022). This extensive compilation allowed them to estimate the rate of change of 𝜉 ion , 0 with redshift, finding it to have a slope of dlog( 𝜉 ion , 0 )/d 𝑧 = 0 . 09 ± 0 . 01 (black dashed line). The markers show the UV-faint galaxies at 𝑧 ∼ 3 -7 from Prieto-Lyon et al. (2023), the Lyman𝛼 emitters (LAEs) from Ning et al. (2023) and Simmonds et al. (2023), the H 𝛼 emitters from Rinaldi et al. (2024), and the ELGs from Simmonds et al. (2024a). The latter show a slightly less steep evolution of 𝜉 ion , 0 with redshift than the compilation of Stefanon et al. (2022), but consistent within errors (given by dlog( 𝜉 ion , 0 )/dz = 0 . 07 ± 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( 𝜉 ion /[Hz erg -1 ]) ≲ 24 . 5. The inclusion of these galaxies has the result of flattening of the increase of 𝜉 ion , 0 with redshift.  \nFigure 8 shows 𝜉 ion , 0 and / 𝑛 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 𝜉 ion , 0 ; Simmonds et al. 2024a). The crosses represent the galaxies with no recent star formation (SFR 10 = 0 M ⊙ yr -1 ). As a reminder, both samples overlap and the Prospector -inferred redshifts agree with the 𝑧 spec ones, for a great majority of the overlapping galaxies. We present both samples fit separately to highlight that the increase in 𝜉 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 𝜉 ion , 0 for the spectroscopic sample has a slope of dlog( 𝜉 ion , 0 )/dz ∼ 0 . 04 ± 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( 𝜉 ion , 0 )/dz ∼ -0 . 02 ± 0 . 00). The same effect can be seen in the best-fit lines to / 𝑛 ion .', '5.3 Trends of ionising properties with UV magnitude': '𝜉 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 ionising radiation (e.g. Duncan & Conselice 2015; Maseda et al. 2020; Simmonds et al. 2024a; Endsley et al. 2024). Figure 9 shows 𝜉 ion , 0 as a function of M 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 ) ≥ -1) and "mini-quenched" (log(SFR 10 /SFR 100 ) < -1), we distinguish two populations of galaxies that have consistent 𝜉 ion , 0 slopes with M UV , but populate a different 𝜉 ion , 0 range. The star-forming sample lies above log( 𝜉 ion , 0 /[Hz erg -1 ]) ∼ 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/parenleft.upYY/parenright.up , where XX is the number of star-forming galaxies and YY is the number of mini-quenched galaxies. The best fit to both datasets (when enough points to fit a line reliably) are shown in the bottom of each panel. Wefind that in general, there is a slight increase of 𝜉 ion , 0 towards the fainter galaxies (slope of ∼ 0 . 01 -0 . 04), but this increase is lower than the slope of ∼ 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 𝜉 ion , 0 -𝑧 relation shown in the top panel of Figure 8 (for the photometric sample) becomes:  \nlog ( 𝜉 ion , 0 (z) ) = (-0 . 001 ± 0 . 004 ) 𝑧 + ( 25 . 294 ± 0 . 017 ) (4)  \nfor the star-forming sample. Likewise, the / 𝑛 ion -𝑧 relation shown in the bottom panel of Figure 8, for the star-forming sample becomes:  \nlog ( / 𝑛 ion (z) ) = ( 0 . 032 ± 0 . 016 ) 𝑧 + ( 53 . 106 ± 0 . 079 ) . (5)  \nAs with 𝜉 ion , 0 , the mini-quenched galaxies also appear to populate a distinct region in the / 𝑛 ion versus UV magnitude plane. Contrary to 𝜉 ion , 0 , however, / 𝑛 ion has a steeper dependence with M UV (with a slope of ∼ -0 . 4), as shown in Figure 10 1 . In summary, as galaxies become fainter, 𝜉 ion , 0 marginally increases, but / 𝑛 ion significantly decreases. This is an effect that derives naturally from the definition of 𝜉 ion , 0 : the ratio between the ionising photons being produced ( / 𝑛 ion ) and the non-ionising 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.  \nFinally, we combine our redshift and UV magnitude relations to perform a 2-dimensional fit and find a joint relation for the starforming sample, given by:  \n)  \nlog ( 𝜉 ion ( 𝑧, M UV )) = ( 0 . 003 ± 0 . 003 ) 𝑧 + (-0 . 018 ± 0 . 003 ) M UV + ( 25 . 984 ± 0 . 053', '5.4 Unveiling the silent population': 'Wefindatotal of 350 mini-quenched galaxy candidates, corresponding to ∼ 2 . 4% of our total sample. Their number increases with redshift up to 𝑧 = 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 𝑧 = 4 -6 could have a physical origin (Dome  \net al. 2023): 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 𝑧 ∼ 7 . 9 galaxy with a bursty SFH, with evidence of a mini-quenched episode lasting ∼ 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 miniquenched systems. On the other hand, at 𝑧 ∼ 5 -7, the UV part of the spectrum is shifted into the NIRCam F090W filter, which covers our entire stellar mass complete sample (see Figure 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 ( 𝑧 > 7), it becomes increasingly difficult to detect these almost featureless galaxies with photometry alone (and with a S/N > 3).  \n"Mini-quenched galaxies" have low 𝜉 ion , 0 , well below the observed 𝜉 ion , 0 relations (i.e. with redshift and M UV ), and are likely only minor players in the reionisation 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 Looser et al. (2023a); Strait et al. (2023) and Looser et al. (2023b), which are temporarily quenched due to their extremely bursty SFHs (Dome et al. 2023). Interestingly, using simulations, Dome et al. (2023) predict the number of mini-quenched galaxies increases with cosmic time, reaching ∼ 2 -4% at 𝑧 = 4, in broad agreement with our preliminary findings. A spectroscopic follow-up study is needed to confirm this hypothesis.  \nIf we adopt the 𝛽 -f esc (LyC) relation from Chisholm et al. (2022), given by:  \nf esc (LyC) = ( 1 . 3 ± 0 . 6 ) × 10 -4 × 10 (-1 . 22 ± 0 . 1 ) 𝛽 , (6)  \nwe find that a significant number of the mini-quenched galaxy candidates [51 (141)] indicate leakage of ionising radiation [f esc (LyC) > 10% (5%)]. Corresponding to ∼ 16% (40%) of the mini-quenched sample. By comparison, the star-forming sample only has ∼ 4%(32%) of galaxies that obey the same criteria. Figure 11 shows 𝜉 ion , 0 as a function of 𝛽 for all galaxies with f esc (LyC) > 10%. The miniquenched candidates populate a different parameter space in 𝜉 ion , 0 and 𝛽 than the star-forming galaxies, with the mini-quenched galaxies having on average bluer UV continuum slopes and lower 𝜉 ion , 0 values. We now focus on the strongest leaker candidates within the mini-quenched sample (f esc > 10%).  \nObserving LyC leakage directly is impossible at high redshift, the average IGM transmission of hydrogen ionising photons emitted at 𝜆 rest -frame ∼ 900 Å, at 𝑧 ∼ 6, is virtually zero (Inoue et al. 2014). It is even more complicated to observe hydrogen ionising photons emitted at shorter wavelengths, for example at ∼ 700 Å, where the nebular contribution does not contaminate f esc estimations (Simmonds et al. 2024b). Therefore, indirect methods are required to understand how ionising 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 LyC leakage in the fitting routing (Stoffers et al. in prep). This is work in progress and is in the process of being calibrated. Promisingly, Prospector finds signs of leakage for these LyC leaking candidates (finding f esc > 20% for all of them).  \nFigure 12 shows 𝜉 ion , 0 versus f esc (LyC), colour-coded by burstiness. The shaded region shows the area where f esc > 10%, and log(SFR 10 /SFR 100 ) < -1 (the latter roughly coincides with log( 𝜉 ion , 0 /[Hz erg -1 ]) < 24.5). This parameter space is analogous  \nFigure 11. 𝜉 ion , 0 versus 𝛽 for the galaxies with f esc (LyC) > 10%, adopting the relation from Chisholm et al. (2022). The two colours represent star-forming ("SF"; smaller grey circles) and mini-quenched ("MQ": purple larger circles) populations, as shown in the legend. The 𝜉 ion , 0 and 𝛽 medians of each subsample are shown as filled (mini-quenched) and dotted (star-forming) lines. Although, by definition, every galaxy with f esc (LyC) > 10% has a blue UV continuum slope, the mini-quenched candidates have on average bluer 𝛽 and lower 𝜉 ion , 0 values than the star-forming galaxies.  \n<!-- image -->  \nto the one presented in Katz et al. (2023) for "remnant leakers". In short, Katz et al. (2023) propose two modes of LyC leakage (f esc > 20%), 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 ionisation bounded HII 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 HII regions, where the ISM was disrupted enough to halt star formation.  \nFigure 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 ∼ 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 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 miniquenched category (< 3% of the total sample), and their low 𝜉 ion , 0 , we conclude that they are not dominant agents in the reionisation of the Universe.', '6 IMPLICATIONS FOR THE REIONISATION OF THE UNIVERSE': 'Recent works based on JWST observations have arrived to the conclusion that there is an overestimation of ionising photons at the EoR, leading to a so called "crisis" in the ionising photon budget (see e.g.,  \nFigure 12. 𝜉 ion , 0 versus LyC escape fractions, estimated using the 𝛽 relation from Chisholm et al. (2022), and colour-coded by burstiness. The vertical dashed line shows the mean f esc of the star-forming sample. As in previous figures, the crosses show the galaxies with zero recent star formation (SFR 10 = 0 M ⊙ yr -1 ). The dark shaded area highlights the region where log(SFR 10 /SFR 100 ) < -1, and f esc (LyC) > 10%, which broadly coincides with the region where log( 𝜉 ion , 0 /[Hz erg -1 ]) is below 24.5. This parameter space is analogue to the one populated by the "remnant leakers" (RL) presented in Katz et al. (2023). Following this criteria, we find 51 remnant leaker candidates.  \n<!-- image -->  \nFigure 13. Stacked SFHs of the 51 galaxies populating the "remnant leakers" parameter space shown in Figure 12. The SFHs have been normalised to their maximum, the individual SFHs are shown as thin coloured curves, and the stack as a purple thick curve. The shape of the mean SFH of these galaxies indicates a burst of star formation in their recent past (<50 Myr), but not within their immediate past (<10 Myr).  \n<!-- image -->  \nFigure 14. / 𝑁 ion as a function of M UV , separated by redshift bins and assuming a f esc as indicated in each panel. The 𝜉 ion , 0 relations with M UV shown in Figure 9 are adopted (using the relations given for galaxies with log(SFR 10 /SFR 100 ) ≥ -1, which are representative of the general population). The variable f esc prescriptions come from Anderson et al. (2017) and Chisholm et al. (2022), and the luminosity functions are taken from Bouwens et al. (2021). In the top panel, we show for comparison the results adopting a fixed log( 𝜉 ion , 0 /[Hz erg -1 ]) = 25.2 (grey shaded area and dotted lines).  \n<!-- image -->  \nTrebitsch et al. 2022; Chakraborty & Choudhury 2024; Muñoz et al. 2024). In the following sections, we discuss the implications of our 𝜉 ion , 0 estimations for a stellar mass complete sample on the cosmic ionising budget.', '6.1 Constraints on the cosmic ionising budget': 'As in Simmonds et al. (2024a), we now study the implications of our findings on the ionising cosmic budget through / 𝑁 ion , which represents the number of ionising photons produced per volume unit. We first analyse the contributions of different M UV and redshift bins to / 𝑁 ion , by adopting the UV luminosity functions from Bouwens et al. (2021) and some simple prescriptions for f esc : constant of 10 and 20% (Ouchi et al. 2009; Robertson et al. 2013, 2015), and varying according to M UV (Anderson et al. 2017). The latter was constructed using the uniform volume simulation Vulcan , and states the dependence: log(f esc ) = ( 0 . 51 ± 0 . 4 ) M UV + 7 . 3 ± 0 . 08. For 𝜉 ion we use the 𝜉 ion , 0 relations presented in Figure 9 for the star-forming sample (i.e. log(SFR 10 /SFR 100 ) ≥ -1). This sample was chosen because it is representative of most galaxies in this study, accounting for >97% of the total sample. We remind the reader that our Prospector -inferred ionising photon production efficiencies assume f esc =0, such that 𝜉 ion , 0 = 𝜉 ion × (1-f esc ). For simplicity, in the following calculations we adopt 𝜉 ion = 𝜉 ion , 0 , therefore, the resulting / 𝑁 ion values should be taken as lower limits. For reference, the mean f esc for the star-forming sample is ∼ 3% (see Figure 12), which would result in a difference in / 𝑁 ion of <0.1 dex.  \nThe resulting / 𝑁 ion curves as a function of M UV are shown in Figure 14, each panel displaying a different f esc as indicated in the  \ntop left corner. As expected, the curves with fixed f 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 𝜉 ion evolution with M UV . Once all types of galaxies are potentially included, the importance of faint galaxies in the ionising cosmic budget is reduced. We note that our exclusion of mini-quenched galaxies does not change this conclusion, since they only account for ≲ 4% of the sample. In the case of the variable f esc , the curves are still steep, due to the nature of their f esc prescriptions. On the one hand, Anderson et al. (2017) states that ionising 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 ( 𝛽 ) are likely more efficient LyC leakers. The 𝛽 in Equation 6 can be expressed as a function of M UV as:  \n𝛽 = 𝑏 + 𝑎 × ( M UV + 19 . 5 ) (7)  \nAs described in Chisholm et al. (2022), the coefficients 𝑎 and 𝑏 connect M UV to the UV continuum slope, 𝛽 , through the 𝛽 -M 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( 𝜉 ion /[Hz erg -1 ]) = 25.2 (based on stellar populations, as in Robertson et al. 2013). We find that the shapes of the / 𝑁 ion curves with fixed 𝜉 ion are similar to the ones found in this work, albeit slightly offset towards lower / 𝑁 ion .  \nUsing the / 𝑁 ion curves estimated above, in Figure 15 we investigate how the ionising cosmic budget evolves with redshift, given by:  \n/ 𝑁 ion = M UV , max ∫ M UV , min 𝜌 UV ( M UV , 𝑧 ) × f esc ( M UV , 𝑧 ) × 𝜉 ion ( M UV , 𝑧 ) dM UV ,  \n<!-- image -->  \nFigure 15. Cosmic rate of ionising photons emitted per second and per unit volume ( / 𝑁 ion ), as a function of redshift. The error bars in our estimations have been omitted for readability, and are on average ∼ 0 . 3 dex. The same f esc prescriptions assumed in Figure 14 are adopted, as indicated in the legend. The results obtained by adopting the UV luminosity density from Sun & Furlanetto (2016) are shown as filled circles ("S16"), and those obtained by integrating the UV luminosity density curves from Bouwens et al. (2021) down to M UV = -16 are shown as triangles ("B21", curves shown in Figure 14). The stars have been obtained by convolving the stellar mass functions from Weibel et al. (2024) instead ("W24"). We include the curve from Mason et al. (2015) assuming constant 𝜉 ion and f esc , integrated down to a M UV of -15 (as in Mason et al. 2019), as well as the / 𝑁 ion reported in Rinaldi et al. (2024) for H 𝛼 emitters at 𝑧 ∼ 7 -8 (black square). We also include the estimated / 𝑁 ion needed to maintain Hydrogen ionisation in the IGM (Madau et al. 1999), adopting clumping factors of 1, 3 and 10. Finally, as comparison, we show the results from Simmonds et al. (2024a), as a shaded hatched area. These were calculated under the assumption that the ELGs studied in that work were representative of the general galaxy population, with low faint low-mass bursty galaxies creating a turnover in / 𝑁 ion at 𝑧 > 8. Finally, we include observational constraints obtained by observing the Ly 𝛼 forest from Becker & Bolton (2013), Gaikwad et al. (2023), and Davies et al. (2024). With our stellar mass complete sample, we find that galaxies produce enough ionising radiation to ionise the Universe by 𝑧 ≈ 5 -6, without producing an excess in the cosmic ionising photon budget. Importantly, the points estimated adopting the Chisholm et al. (2022) f esc prescription, flatten at 𝑧 ≲ 6, in general agreement with the Ly 𝛼 forest constraints.  \n<!-- image -->  \n(8)  \nwhere / 𝑁 ion is in units of s -1 Mpc -3 , 𝜉 ion is in units of Hz erg -1 , and 𝜌 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. 𝜉 ion -M UV relation). In this figure we integrate the curves shown in Figure 14 down to M 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 UV ∼ -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 𝜉 ion we use the best-fit relation to our data, as shown in Figure 9. For consistency, we use only the star-forming galaxies (accounting for ∼ 94% of the total sample).  \nWhen adopting the Sun & Furlanetto (2016) luminosity functions, we only assume fixed escape fractions of 10 and 20% (filled circles), motivated by the canonical average f esc needed for galaxies to ionise the Universe (Ouchi et al. 2009; Robertson et al. 2013, 2015). As comparison, we add curves indicating the / 𝑁 ion needed to maintain ionisation 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 ionising photons need to be emitted in order to sustain ionisation. 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  \nimportant 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 ionising photons inferred from JWST observations is identified in the cosmic ionising 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, 2019; Naidu et al. 2020; Rinaldi et al. 2024).  \nIn addition to the measurements described above, we include those found through the stellar mass function (instead of 𝜌 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 𝑧 ∼ 4 -9. In this scenario, / 𝑁 ion can be rewritten as:  \n/ 𝑁 ion = M max ∫ M min Φ ( M , 𝑧 ) × f esc ( M , SFR(M) , R(M) ) × / 𝑛 ion ( M , 𝑧 ) dM , (9)  \nwhere Φ is in units of Mpc -3 , f esc is dimensionless and depends on stellar mass, SFR, and size of the galaxy (through the SFR surface density, Σ SFR ), and / 𝑛 ion is in units of s -1 . In this case, we adopt the f esc prescription from Naidu et al. (2020) due to its dependence on Σ SFR , such that:  \nf esc = 𝑚𝑖𝑛 GLYPH<18> 1 , 1 . 6 + 0 . 3 -0 . 3 × GLYPH<16> Σ SFR 1000M ⊙ yr -1 kpc -2 GLYPH<17> 0 . 4 + 0 . 1 -0 . 1 GLYPH<19> (10)  \nIn order to obtain Σ 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  \nFigure 16. Percentage of cosmic rate of ionising photons being emitted per second and per unit volume, as a function of M UV (four leftmost columns, coloured) and stellar mass (rightmost column, black and white). Three bins are shown per redshift in each panel, as indicated in the legends. The three M UV bins correspond to those shown in Figure 6, where brighter galaxies tend to have higher stellar masses. Each row is normalised to the Madau et al. (1999) models shown in Figure 15, where the / 𝑁 ion needed to maintain hydrogen ionisation at every redshift bin is set to 100%: from top to bottom, C=1, C =3 and C=10. The columns show the different f esc prescriptions adopted earlier: constant (10 and 20%), and varying as a function of: M UV following Anderson et al. (2017), 𝛽 as proposed in Chisholm et al. (2022), and Σ SFR from Naidu et al. (2020). It can be seen that as the clumping factors increase, it becomes more difficult for galaxies to produce enough ionising photons by the end of the EoR ( 𝑧 ∼ 5 -6). Importantly, for the calculations made convolving 𝜌 UV (four leftmost columns), the faintest galaxies (fainter than M UV = -20) dominate the ionising photon budget for every f esc prescription and clumping factor. Analogously, the models on the rightmost column indicate that the intermediate and low mass bins dominate the budget, with a significant contribution from the galaxies in the lowest mass bin (log(M ★ /[M ⊙ ])<8.5). Therefore, despite the uncertainties in the clumping factor, we confirm that faint galaxies with low stellar masses are key agents of reionisation.  \n<!-- image -->  \nover the past 10 Myr), such that Σ SFR = SFR / 2 𝜋 R 2 (Shibuya et al. 2019). After estimating Σ SFR , we compute f esc ×/ 𝑛 ion for our stellar masscomplete sample, for each redshift bin. By combining the best-fit relations with Φ , and integrating them down to log(M ★ /[M ⊙ ])=7.5, we finally find / 𝑁 ion as a function of redshift. The values are shown as stars in Figure 15, and are in general agreement with the other variable f esc prescriptions (i.e. Anderson et al. 2017; Chisholm et al. 2022).  \nAlthough the values of / 𝑁 ion are highly uncertain in the early Universe, strong constraints on / 𝑁 ion have been placed by observations of the Ly 𝛼 forest at lower redshifts (e.g. Becker & Bolton 2013; Gaikwad et al. 2023; Davies et al. 2024, included in Figure 15). The Ly 𝛼 forest refers to a series of absorption lines at wavelengths redder than Ly 𝛼 ( 𝜆 rest -frame ∼ 1216 Å) observed in the spectra of high redshift quasars, produced by the intervening neutral IGM. Interestingly, / 𝑁 ion has been observed to flatten once reionisation has been completed (with log( / 𝑁 ion /[s -1 Mpc -3 ]) ∼ 50 . 8 at 𝑧 ≲ 5 -6). Wefind that, out of all the prescriptions adopted in this work, the f esc relations from Chisholm et al. (2022) can match the shape of the Ly 𝛼 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 ★ /[M ⊙ ]) ∼ 7 . 5). Therefore, our points represent lower limits to the cosmic ionising budget.', '6.2 Which galaxies reionised the Universe?': 'In order to determine which galaxies dominate the budget of reionisation, we use the same f esc prescriptions as before, but integrate their contributions in three M UV bins: -24 ≤ M UV < -20, -20 ≤ M UV ≤ -18, and -18 < M UV ≤ -16. Figure 6 shows that M 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 Figure 16 we adopt these M UV luminosity bins and show the relative contribution these bins have in / 𝑁 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 ionising budget, by adopting the stellar mass functions of Weibel et al. (2024) and the f esc prescription from Naidu et al. (2020) (stars in Figure 15). Each row has been normalised 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 ionising radiation in order to sustain hydrogen ionisation, independent of the f esc prescription. As the clumping factor increases, it becomes more difficult for galaxies to ionise the Universe by redshift 5. However, most importantly, we discover that for every clumping factor, f esc assumption, and redshift bin, the fainter galaxies (with M UV > = -20), and galaxies with low  \nand intermediate stellar masses (log(M ★ /[M ⊙ ])<9.5) dominate the cosmic ionising 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; Vĳayan et al. 2021), and the conclusions reached by forward modelling JWST analogues from the SPHINX simulation, presented in Choustikov et al. (2024).  \nTherefore, in this work, we confirm that faint low-mass galaxies with bursty star formation have in general enhanced 𝜉 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 ionising photon budget of the Universe (e.g., Pahl et al. 2024; Trebitsch et al. 2022; Chakraborty & Choudhury 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 𝛼 forest, especially in the flattening of / 𝑁 ion that has been observed at lower redshifts ( 𝑧 ≲ 6). In summary, we have shown that galaxies produce enough ionising photons to ionise the Universe by 𝑧 ≈ 5 -6, without creating a nonphysical excess of ionising photons in the cosmic budget.', '7 CAVEATS AND LIMITATIONS': 'There are a few important limitations to our method, which we now describe. Firstly, 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 ⊙ . 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. (2023a) have found tentative evidence for a top heavy IMF at 𝑧 ∼ 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 𝜒 2 < 1, but we note that by doing so, might have lost information on sources that cannot be reproduced by our models.  \nIn 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. 2024; Ü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 ∗ ≤ 10 9 . 5 M ⊙ ). 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 alon. 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 AGN that can be mimicked by stellar emission.  \nAnother 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  \nNIRCam 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 perfectagreement, suggesting that our results are in general accurate but likely contain a few outliers where the 𝑧 phot is unreliable. We stress that these outliers would be of the order of at most a few percent. As with the stellar populations, our selection of galaxies (fitted with a reduced 𝜒 2 ≤ 1) ensures that the photometry is well represented by the best fit 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.  \nIn 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 ≤ 𝑧 ≤ 9).', '8 CONCLUSIONS': 'We use JWST NIRCam photometry to build a sample of 14652 galaxies at 3 ≤ 𝑧 ≤ 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 < 3% of the total sample. Within the mini-quenched population, we find an interesting subsample with tentative evidence of LyC leakage (through the UV continuum slope 𝛽 ). 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 summarised as follows.  \nWefindthat 𝜉 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 𝜉 ion , 0 with 𝑧 for the more representative star-forming sample is:  \nlog ( 𝜉 ion , 0 (z) ) = (-0 . 001 ± 0 . 004 ) 𝑧 + ( 25 . 294 ± 0 . 017 )  \nThe 2-dimensional fit that accounts for the change of 𝜉 ion with M UV and redshift, for the same sample is given by:  \nlog ( 𝜉 ion , 0 ( 𝑧, M UV )) = ( 0 . 003 ± 0 . 003 ) 𝑧 + (-0 . 018 ± 0 . 003 ) M UV + ( 25 . 984 ± 0 . 053 )  \nTo study the contribution of the galaxies in this study to reionisation, we convolve the star-forming relations (which represent > 97% of the total sample), with luminosity functions from literature. We find that galaxies, which are detected with JWST, can ionise 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 UV = -20 contribute similar amounts of ionising photons (see Figure 14), and that galaxies in the range of M UV = -20 to -16 dominate the budget of reionisation at every redshift bin studied in this work (see Figure 16). With our stellar mass complete sample, our predictions do not overestimate / 𝑁 ion for galaxies at 𝑧 > 8 (see Muñoz et al. 2024). We note that if we extrapolate our trends to fainter magnitudes (to M UV of -14 or -12), the Universe can be reionised with lower escape fractions. Promisingly, by adopting the relation of f esc with M UV presented in Chisholm et al. (2022), we can conciliate our results regarding the ionising  \ncosmic budget with the constraints obtained through observations of the Ly 𝛼 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-exclusiveaccess period.  \nCS thanks Sarah Bosman, Frederick Davies, James Leftley, and Arjen van der Wel for the insightful discussions. CS, W.B., RM, & 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 & 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.', 'REFERENCES': '2011, ApJ, 735, L34  \nTable A1. Photometric offsets defined as the ratio between the observed and modelled photometry. Column /one.up/colon.up name of filter. Colum /two.up/colon.up median offset and errors, given by the 16 th and 84 th percentiles.  \nZhu Y., et al., 2024, arXiv e-prints, p. arXiv:2405.12275', '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 Figure 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 𝜒 values well outside of the bounds of the figures. 𝜒 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-fit models. Figure 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': 'Figure B1 shows the mean shape of the main properties used to estimate the stellar mass completeness in Section 4.2. The redshift and  \nstellar masses have been divided by the 50 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 ionising photon production efficiencies in our sample, we use the Kullback-Leibler definition of information gain (IG), given by:  \nIG = ∫ Posterior ( 𝑝 ) × log 2 Posterior ( 𝑝 ) Prior ( 𝑝 ) d 𝑝 [bits] , (C1)  \nwhere 𝑝 represents the parameter of choice: either log M ★ or log 𝜉 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):  \n- · IG < 1: little-to-no information was gained\n- · 1 ≤ IG ≤ 2: some information was gained\n- · IG > 2: parameter is constrained  \nFigures C1 and C2 show the IG for stellar mass and 𝜉 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 𝜉 ion , 0 , we find a trend of increasing IG with burstiness, but emphasise that all measurements are highly constrained (with IG ≫ 2). Promisingly, most mini-quenched galaxies discussed in this work have a large IG. This is due to the richness of the photometric dataset, which allows to constrain SED shapes even in cases with no obvious emission lines.', 'APPENDIX D: BEST-FIT PARAMETERS FOR RELATIONS OF IONISING PROPERTIES WITH OBSERVED UV MAGNITUDE': 'For readability, the relations of 𝜉 ion , 0 and / 𝑛 ion , as a function of M UV , presented in Figures 9 and 10, are given in Table D1 and D2, respectively.  \nThis paper has been typeset from a T E X/L A T E X file prepared by the author.  \nFigure A1. Comparison between the modelled and observed HST (squares) and JWST NIRCam (circles) photometry for our stellar mass complete sample. 𝜒 scatters around zero with mostly symmetric error bars.  \n<!-- image -->  \nFigure A2. Same as Figure A1 but at rest-frame wavelengths. The grey dots show the individual measurements, while the white circles and error bars show the median and errors when adopting bins of width 100Å. For comparison, the Ly 𝛼 , [O /i.pc/i.pc], H 𝛽 , [O /i.pc/i.pc/i.pc], and H 𝛼 emission wavelengths are shown as vertical dashed lines. It can be seen that the median values scatter around 𝜒 = 0, with symmetrical error bars. Moreover, we find no evidence of emission line fluxes being over (or under) estimated by Prospector .  \n<!-- image -->  \nTable D1. Best-fit parameters for log( 𝜉 ion , 0 )-M UV relations, with form log( 𝜉 ion , 0 ) = 𝛼 M UV + log( 𝜉 ion , int ) shown in Figure 9, where 𝜉 ion , int is the intercept of the relation. Column /one.up/colon.up redshift bin. Column /two.up,/three.up/colon.up slope and 𝜉 ion , 0 normalisation for the star-forming (SF) sample. Column /four.up,/five.up/colon.up slope and 𝜉 ion , 0 normalisation for the mini-quenched (MQ) sample. For 𝑧 ≥ 7 we only provide relations for the star-forming sample, since there are not enough mini-quenched galaxies to obtain a reliable fit.  \nTable D2. Best-fit parameters for log( / 𝑛 ion )-M UV relations, with form log( / 𝑛 ion ) = 𝛼 M UV + log( / 𝑛 ion , 0 ) shown in Figure 10. Column /one.up/colon.up redshift bin. Column /two.up,/three.up/colon.up slope and / 𝑛 ion normalisation for the star-forming (SF) sample. Column /four.up,/five.up/colon.up slope and / 𝑛 ion normalisation for the mini-quenched (MQ) sample. For 𝑧 ≥ 7 we only provide relations for the star-forming sample, since there are not enough mini-quenched galaxies to obtain a reliable fit.  \nFigure B1. Corner plot showing the main properties of interest in the stellar mass completeness estimation: redshift and stellar mass. The values have been normalised to the 50 th percentile (dashed vertical lines) for each galaxy, in order to understand the general shape of the posteriors in our sample. The grey areas show the 16 th and 84 th ranges. We find no significant signatures of asymmetry.  \n<!-- image -->  \n<!-- image -->  \nFigure C1. Kullback-Leibler information gain as a function of flux in the F444W band. A higher flux leads in general to a higher IG. However, all of the galaxies in our sample have IG > 1, with ∼ 94% above IG = 2. Left panel/colon.up colour-coded by stellar mass. Right panel/colon.up colour-coded by burstiness of their SFH, defined as SFR 10 /SFR 100 . The crosses show galaxies with no recent star formation (SFR 10 = 0).  \n<!-- image -->  \n<!-- image -->  \nFigure C2. Same as Figure C1 but for the ionising photon production efficiency, 𝜉 ion , 0 . The information gain tends to increase towards galaxies with burstier SFHs, however, we note that all galaxies in our stellar mass complete sample have a high information gain (well above IG = 2, shown as a dashed horizontal line).  \n<!-- image -->'}