The production of ionizing photons in UV-faint z ∼ 3 − 7 galaxies

,


Introduction
In recent years, we have obtained increasing evidence that the reionization of hydrogen happened fairly late, approximately one billion years after the Big Bang (z ∼ 5.5 − 10), with a midpoint around z ∼ 7 − 8 (e.g., Fan et al. 2006;Stark et al. 2010;McGreer et al. 2015;Mason et al. 2018;Davies et al. 2018;Qin et al. 2021;Planck Collaboration et al. 2020;Bolan et al. 2022).However, there is evidence for significant star formation before E-mail: gonzalo.prieto@nbi.ku.dk this time (e.g., Oesch et al. 2018;Hashimoto et al. 2018;McLeod et al. 2021) thus it appears that reionization lags behind galaxy formation.The reason for this lag is unknown: we are still lacking a full physical understanding of the reionization process.In particular, we still do not know which types of galaxies drive the process, i.e. which physical mechanisms mediate the production and escape of ionizing photons from galaxies.In order to produce such a late and fairly rapid reionization, the ionizing population could have been dominated by low mass, UV-faint galaxies with a low (∼ 5%) average escape fraction (e.g., Mason et al. 2019;Qin et al. 2021).Alternatively, rarer more massive galaxies with higher escape fractions could have been responsible (e.g., Sharma et al. 2017;Naidu et al. 2020).With only measurements of the timing of reionization, these scenarios are degenerate, thus physical priors on the ionizing properties of galaxies across cosmic time are necessary to pinpoint the sources of reionization.
The total ionizing output of galaxies can be simply parameterized (e.g., Madau et al. 1999;Robertson et al. 2010) as the product of the production rate of ionizing photons relative to non-ionizing UV photons, ξ ion (determined by the stellar populations, e.g., Stanway et al. 2016), and the fraction of ionizing photons which escape the interstellar medium (ISM) into the intergalactic medium (IGM), f esc (determined by the structure and ionization state of the ISM, likely shaped by star formation and feedback, e.g., Trebitsch et al. 2017;Ma et al. 2020).Both of these quantities are also expected to vary with time in an individual galaxy, for example due to the lifetime and properties of young stellar populations, and the effects of feedback and bursty star formation on the ISM.
While we can easily observe the non-ionizing UV photons from galaxies, the high optical depth of the IGM to ionizing photons makes direct measurements of the escaping ionizing spectrum statistically unlikely at z ∼ > 3 (Inoue et al. 2014;Becker et al. 2021;Vanzella et al. 2018).Alternatively, fluxes of nonresonant recombination lines, emitted by gas which was ionized in HII regions around massive stars, can crucially measure the flux of ionizing photons which do not escape galaxies.In particular, Hα emission can be used to directly estimate (1 − f esc )ξ ion (e.g., Leitherer & Heckman 1995;Bouwens et al. 2016;Shivaei et al. 2018;Emami et al. 2020).As f esc is inferred to be low ( ∼ < 10%) on average for Lyman-break galaxies at z ∼ 2 − 4 (Steidel et al. 2018;Begley et al. 2022;Pahl et al. 2022) measurements of Hα should trace the intrinsic production of ionizing photons reasonably well.ξ ion can also be inferred from the strength of [OIII]+Hβ emission, though due to the dependence of [OIII] emission on metallicity and ionization parameter, the correlation is not as tight as with Hα (e.g., Chevallard et al. 2018).
However, due to the limited spatial resolution and sensitivity of Spitzer, previous works had been limited to studying ξ ion in isolated, bright (> L * ) galaxies where deblending IRAC photometry was possible (e.g,.Bouwens et al. 2016) and stacks for fainter galaxies (e.g., Lam et al. 2019;Maseda et al. 2020).With JWST it is finally possible to extend these studies to individual UV-faint galaxies (Endsley et al. 2022), and obtain rest-frame optical spectroscopy at z > 3 (e.g., Sun et al. 2022;Williams et al. 2022).
Results from previous analyses have been intriguing but require further investigation.Using stacked IRAC photometry Lam et al. (2019) found no significant evidence for a strong correlation of ξ ion with M uv .However, Maseda et al. (2020) found a population of extremely UV-faint (M uv > −16) galaxies selected as Lyα emitters in deep MUSE observations, which have very elevated ξ ion compared to higher luminosity galaxies and at fixed Hα EW, implying these efficient ionizing galaxies are particu-larly young and low metallicity.It is thus important to examine the distribution of ξ ion at low UV luminosities, and to compare galaxies with and without Lyα emission to better understand the demographics of the ionizing population.
Furthermore, using early JWST NIRCam data Endsley et al. (2022) discovered a population of UV-faint (M uv ∼ −19) galaxies at z ∼ 6.5 − 8 with high sSFR but low EW [OIII]+Hβ inferred from photometry.The high sSFR would imply high ξ ion due to the increased abundance of O and B stars.To explain the low [OIII]+Hβ EW, Endsley et al. (2022) suggest either these galaxies have extremely low metallicities (reducing oxygen abundance), or alternatively, all nebular lines are reduced.A reduction in all nebular lines could be due to either them being produced in density-bounded HII regions with very high ionizing escape fraction (e.g.Zackrisson et al. 2013;Marques-Chaves et al. 2022), or recent cessation of star formation.At z ∼ 3 − 7 both [OIII]+Hβ and Hα are visible in NIRCam photometry, enabling us to test these scenarios.
In this paper we make use of deep multi-band HST/ACS, WFC3 and JWST/NIRCam imaging with overlapping MUSE observations, enabling us to blindly detect a spectroscopic sample with precision rest-frame ultra-violet to optical photometry.We measure the distribution of ξ ion over a broader luminosity range (−23 ∼ < M uv ∼ < − 15.5) than previously possible in individual galaxies, due to the excellent resolution and sensitivity of NIRCam at rest-optical wavelength compared to Spitzer/IRAC, and the power of gravitational lensing.We explore correlations of ξ ion with empirical galaxy properties.We find significant trends of increasing ξ ion with decreasing UV luminosity, UV β slope and with increasing Hα EW, all implying the strongest ionizers are young sources with expected low metallicities.We also explore whether our sample shows evidence for very low metallicities or extremely high escape fraction.
The paper is structured as follows.In Section 2 we describe the photometric and spectroscopic data for our study.In Section 3 we describe how we infer the ionizing production rate ξ ion and we describe our results of correlations between ξ ion and other galaxy properties, and comparison to the literature in Section 4. We discuss our results and state our conclusions in Section 6.

Data
For this work we select fields with multi-band HST/ACS and JWST/NIRCam imaging and overlapping MUSE spectroscopy.We select sources detected with Lyα emission (z ∼ 2.9 − 6.7 in MUSE) and sources with high probability of being in the same redshift range based on photometric redshift, and use the HST + JWST photometry to extract optical emission line fluxes.Below we describe the datasets and the selection of our sample.

Imaging
We use JWST NIRCam imaging in parallel to and of the cluster Abell 2744 from the GLASS-JWST program ERS-1324 (PI Treu Treu et al. 2022) and UNCOVER 1 program GO-2561 (co-PIs Labbé, Bezanson).
In our analysis, we also include new and archival HST imaging, from which ACS imaging is particularly important for constraining photometric redshifts.This includes new HST/ACS data in F606W (59530 s.), F775W (23550 s.), and F814W (123920 s) from HST-GO/DD program 172312 (PI Treu), as well as archival data acquired under the Hubble Frontier Fields program (HST-GO/DD-13495, PI Lotz, Lotz et al. 2017), BUF-FALO (HST-GO-15117, PI Steinhardt, Steinhardt et al. 2020) and programs HST-GO-11689 (PI Dupke), HST-GO-11386 (PI Rodney), HST-GO-13389 (PI Siana), HST-GO-15940 (PI Ribeiro) and HST-SNAP-16729 (PI Kelly).Not all HST bands cover every object in our sample, but we only keep objects in our sample that have a well-constrained photometric redshift, usually meaning that there is ACS coverage (see Section 2.4).We also include HST/WFC3 imaging for completeness, but these are generally not as constraining as the NIRCam fluxes.
The image reduction and calibration, and the methods used to detect sources and measure multi-band photometry in both fields closely follow that of Brammer et al., (in prep).Briefly, we pull calibrated images from the MAST archive3 and process them with the grizli pipeline (Brammer et al. 2022).The pipeline first aligns the exposures to external catalogs and to each other and corrects for any distortion within the image.Following this, we subtract a sky-level background, divide out flatfield structure using custom flat-field images, and correct for 1/ f noise.We also correct for NIRCam image anomalies, which include persistence, any remaining cosmic rays, and 'snowballs' (see Rigby et al. 2022).Finally, we apply zeropoint corrections calculated by G. Brammer4 and drizzle all exposures to a common pixel grid.
For source detection, we use SEP (Barbary 2018) to perform aperture photometry on the F444W detection image in each field.

VLT/MUSE spectroscopy
MUSE spectroscopy in the Abell 2744 cluster were obtained through ESO program 094.A-0115 (PI Richard) and are described by Mahler et al. (2018) and Richard et al. (2021).We use their publicly available catalog to select Lyα emitting galaxies.The data comprise a 4 sq.arcmin mosaic centered on the cluster core.Four 1 sq.arcmin quadrants were observed for a total of 3.5, 4, 4 and 5 hours respectively, and the center of the cluster was observed for an additional 2 hours.The median line flux 1σ uncertainty in the MUSE data is 3.6 × 10 −19 erg s −1 cm −2 .This corresponds to an 5σ EW limit of ∼ 4 − 30 Å over z ∼ 3 − 7 for a galaxy with M uv = −19 (the median for our sample, before accounting for magnification, as EW is invariant under magnification).
VLT/MUSE spectroscopy in the GLASS-JWST NIRCam fields were obtained through a new ESO DDT program 109.24EZ.001(co-PIs Mason, Vanzella) on the nights of July 28 and August 20 2022.The data comprise 5 pointings (4 pointings -4 sq.arcmin -overlapping with NIRCam imaging) each with 1 hour exposure time.The raw data are publicly available on the ESO archive5 .The reduction, calibration and source detection methods used for this work are identical to techniques described in previous works (Caminha et al. 2017;Caminha et al. 2019).A full assessment of the depth is ongoing but based on the ∼ 4 hour depth of the Mahler et al. (2018) observations described above, we estimate a 5σ EW limit of ∼ 8 − 60 Å in these shallower data.
In this work we use 102 spectroscopic confirmations at z ∼ 2.9−6.7:42 from the GLASS-JWST NIRCam fields and 60 from the Abell 2744 cluster field.

Gravitational lensing magnification
For the galaxies detected in the core of the Abell 2744 cluster we correct for gravitational lensing magnification using the model by Bergamini et al. (2022).The median magnification of the sample is µ = 3.54 with 90% of the galaxies having µ = 2 − 20.We remove sources with a magnification with µ > 50 (12 sources) due to high uncertainties in the model near the critical curves.The galaxies in the parallel fields are ∼ 3 − 10 away from the cluster core where the magnification is expected to be modest (µ ≈ 1), we do not account for magnification of those sources.

Sample selection
For this work, we focus on selecting a sample of galaxies at z ∼ 3 − 7 with high purity.We select 102 MUSE Lyα-detected galaxies with overlapping HST/ACS and JWST/NIRCam data as described above.We also select a comparison sample of galaxies based on peak photometric redshift, within the same footprint as the MUSE observations, which we expect to have slightly lower Hα EW than the Lyα selected sample.
We find the photometric redshift distribution of all sources detected as described in Section 2.1 using EAZY (Brammer et al. 2008), given all available photometric bands.To build a photometric sample with high purity, following Bouwens et al. (2016), we select sources with the peak of their photometric redshift between 2.9 < z < 6.7, and keep only sources which have 90% of the redshift probability density between ∆z ∼ 1 of the peak of their distribution.The resulting high purity photometric sample consists of 268 galaxies.
The redshift and UV magnitude distribution of our sample is shown in Figure 1.The median redshift of the full sample is 4.02, with the Lyα-selected sample having median redshift 3.95.The median M uv is -18.1, with a Kolmogorov-Smirnov (KS) test showing no significant difference between the Lyα and photometrically selected samples.

Inferring nebular emission line strengths from photometry
To estimate nebular emission line fluxes from broad-band photometry, we follow approaches in the literature and fit the spec- .Galaxies studied in this work; in purple Lyα detected galaxies, and in gray galaxies photometrically selected (with no Lyα detected).
Top: Distribution of redshifts for the spectroscopic and photometric samples.We show the spectroscopic redshift where available, or the peak photometric redshift.Bottom: UV magnitude distribution for our sample, we find a median value of −18.14 ± 1.58, with no statistically significant difference between both samples.
tral energy distribution (SED) to the full photometry, excluding bands we expect to contain strong nebular emission lines (e.g., Shim et al. 2011;Stark et al. 2013;Mármol-Queraltó et al. 2016;Bouwens et al. 2016).This provides us with a model for the continuum flux in those bands which we can subtract from the observed photometry to infer the line flux.We use BAGPIPES to fit SEDs (Carnall et al. 2018).We adopt BC03 (Bruzual & Charlot 2003) templates, and exclude any nebular emission contribution.We do not consider any broadbands where Hα or [OIII]+Hβ are observed according to each galaxy's redshift.For ease of comparison to the literature (e.g., Maseda et al. 2020;Lam et al. 2019), we assume a Chabrier (2003) Initial Mass Function and an Small Magellenic Cloud (SMC, Prevot et al. 1984) dust attenuation law, allowing A V to vary from 0−3 mag.Because metallicity is not well known at the range of redshifts we explore, we allow metallicity to vary from 0 − 2 Z .Because star formation histories are notoriously diffi-cult to constrain at high-redshift (Strait et al. 2021), we assume an exponentially rising delayed τ star formation history, allowing τ to vary freely.For the spectroscopically confirmed Lyman α emitters, we fix the redshift at the Lyα redshift.For our photometric sample, we use the photometric redshift obtained from EAZY with uniform prior with ∆z = 1.(see Section 2.4).
We then compare the SED model of the galaxy's continuum to the broadbands where Hα or [OIII]+Hβ fall.We multiply the non-nebular SED posteriors with the transmission of the aforementioned broadbands to obtain the contribution of the galaxy's continuum to the observed flux.By subtracting this continuum flux contribution to the observed photometry, we then are able to recover the flux distributions of Hα and [OIII]+Hβ emission lines for each galaxy.We compared our measurements with a sample of 6 galaxies with [OIII]+Hβ EW measurements from the GLASS-ERS program using JWST/NIRISS (Boyett et al. 2022a), finding that our method recovers the EW of these sources within ∼ 20 − 40%.A full comparison of these photometric inference methods is left to future work.
There are some limitations to our method for obtaining line fluxes, such as contamination from the 4000-Å break in the broadband containing [OIII]+Hβ, the chance that the line falls outside the effective width of any of our broadband filters, or Hα and [OIII]+Hβ falling on the same band.We consider a contribution of 6.8% from [NII] to the calculated Hα flux, and 9.5% from [SII] according to Anders & Fritze-v. Alvensleben (2003).We remove galaxies with a poor χ 2 (>50) score on their SED fit, we choose these value by ignoring all galaxies on the high-end of the χ 2 distribution.
The advantage of this approach, unlike estimating line fluxes directly from the SED fitting, is that it does not depend strongly on star formation history assumptions and allows us to make a mostly empirical measurement of the line fluxes.We obtain comparable results using the flux in the band red-ward of Hα as the continuum flux, assuming a flat optical continuum (see also e.g., Maseda et al. 2020).Estimating other physical parameters such as star-formation rate and stellar mass from the SED fitting did not give reliable results, since it was too dependent on the initial assumptions, and needed extremely young ages (<10Myrs) and an instantaneous burst of star-formation to recreate the observed nebular emissions..
The following results will consist of: 83 and 64 Lyα-emitting galaxies with Hα and [OIII]+Hβ emission line measurements respectively, a photometric sample of 220 and 203 galaxies with Hα, and [OIII]+Hβ emission line measurements respectively.We see both lines in 62 Lyα galaxies and 177 photometrically selected galaxies.We see no apparent biases in our M uv distribution after narrowing down the sample.Nebular emission flux errors are derived from the 68% confidence interval of the resulting distributions.

Measuring UV absolute magnitude and slope
To infer the UV absolute magnitude, M uv (magnitude at 1500Å) and β slope we fit a power law (e.g., Rogers et al. 2013) f λ ∝ λ β to the fluxes from the HST and JWST bands.We perform the fit using a Markov Chain Monte Carlo sampling using the python module emcee (Foreman-Mackey et al. 2013).We assume flat priors for β and M uv , with bounds −4 < β < 1 and −25 < M UV < −12, sufficient to explore the common value ranges for galaxies (e.g., Bouwens et al. 2014).
To obtain the photometric bands that are observing the UVrestframe of our galaxies, we exclude any bands that fall blue-wards from Lyman-α that might be affected by the Lyman-Break.For the same reason, we exclude bands redwards from the 4000Å break in the restframe.After these requirements, we are left with between 3 and 4 bands for each source.In the case of galaxies with Lyman-α detected in MUSE, we use the line's redshift.For photometrically-selected galaxies, in each call of the likelihood, we randomly draw a redshift from a Normal distribution N(µ = z phot , σ = 0.5), and select the appropriate photometric bands.For lensed sources, we consider magnifications and apply them following the same random draw method as redshift.We use the corresponding magnification and error obtained from the Bergamini et al. (2022) lensing model.

Determination of ξ ion
We define the production rate of ionizing photons, ξ ion .This is given by the ratio between the luminosity of observed ionizing photons and the intrinsic luminosity of the ionizing UV photons (e.g., Leitherer & Heckman 1995): where L Hα is the unattenuated Hα luminosity in erg s −1 and L ν,UV,intr is the intrinsic UV luminosity density at 1500Å.The models from where the conversion factor is derived, assume a young population of massive stars equivalent to a massive HII region.We assume this type of environment to be similar to what we would find in young galaxies.
Because Hα is produced by the excitation of hydrogen gas from ionizing radiation that does not escape the galaxy, and we cannot measure f esc directly in our sample, we note that we are obtaining the production rate of ionizing photons which did not escape the galaxy, ξ ion (1 − f esc ).
We first calculate L Hα directly from the SED obtained in Section 3.1, after accounting for dust attenuation (Prevot et al. 1984).To obtain the intrinsic value of the UV luminosity, we take into account the dust attenuation following Lam et al. (2019), where the intrinsic UV luminosity is defined as L intr UV,ν = L UV,ν / f esc,UV .f esc,UV is the fraction of escaping UV photons not absorbed by the dust.For this, we use the Small Magellanic Cloud dust law defined by Prevot et al. (1984): Where β is the UV slope obtained in Section 3.2.Galaxies with slopes bluer than β < -2.23 are assumed dust-free and we do not correct for dust.
In the following, uncertainties on ξ ion are the 68% confidence intervals, obtained from propagating the uncertainty in the Hα flux from its resulting distribution as described in Section 3.1, and the posterior distributions for β and M uv as described in Section 3.2.

Correlation Analysis
For the purpose of studying the correlations between galaxy properties we use the python package linmix 6 to do Bayesian linear regression including intrinsic scatter and accounting for two-dimensional errors (Kelly 2007).We fit for log 10 [(1 − f esc )ξ ion ] = α + βX + , where is the intrinsic scatter which is assumed to be normally distributed with variance σ 2 .We recover the best fit trend line from the posteriors, as well the 68% Fig. 2. Distribution of (1f esc )ξ ion .The purple histogram includes galaxies with; Lyα emission detection, and in grey galaxies without Lyα detected.Overall, Lyα emitting galaxies show stronger ionizing photon production than galaxies with no Lyα emission selected galaxies, with median values 25.39±0.64 and 25.31±0.43respectively.We show the median relation from the literature at z ∼ 2 − 5 as a dashed black line (e.g., Shivaei et al. 2018;Bouwens et al. 2016;Lam et al. 2019) confidence interval on the parameters.We report the results in Table 1 and show the best-fit line on plots.

Results
In the following section we present our results.In Section 4.1 we study the trends between ξ ion , Hα EW, M uv , and β slope, and in Section 4.2 we investigate whether our sample shows evidence for very high ionizing photon escape fraction and/or very low metallicity galaxies.

Behaviour of ξ ion
Figure 2 shows the distribution of (1 − f esc )ξ ion for our Lyαselected and photometric samples.We find median values of log 10 ξ ion 25.39±0.64 and 25.31±0.43respectively; with 25.33±0.47[Hz erg −1 ] for the complete data set.We find an intrinsic scatter of 0.42 dex, obtained by subtracting in quadrature the average uncertainty in log 10 ξ ion (= 0.21 dex) from the standard deviation of the observed distribution.The recovered intrinsic scatter is broader by ∼ 0.1 dex than that found by Bouwens et al. (2016) and Shivaei et al. (2018) in M uv ∼ < − 20 galaxies.The broad distribution of ξ ion is likely an outcome of the broad range of stellar populations in these galaxies, i.e. due to a range of star formation histories and thus ages, and stellar metallicity (see e.g., Shivaei et al. 2018).
We perform a two sample Kolmogorov-Smirnov (KS) test to test whether the Lyα-selected and photometric samples are drawn from the same distribution.We recover p-value of 0.03, meaning it is likely that the underlying distributions are different, consistent with the results from (Saldana-Lopez et al. 2022), where a statistically significant difference is found between the ξ ion distributions of LAEs and non-LAEs at z ∼ 3 − 5. Given that galaxies with strong Lyα emission also likely have high ionizing photon escape fractions (e.g., Verhamme et al. 2015;Dijkstra et al. 2016) it is likely that the intrinsic ionizing photon produc-  2016) as colored boxes for comparison.We find evidence for an increase in log 10 (1 − f esc )ξ ion towards fainter UV magnitude, with a slope of 0.03 ± 0.02, only considering the range where our sample is M uv complete (M uv < −18.1).
tion efficiency of these galaxies is even higher than what we can infer based on Hα emission.Figure 3 shows (1 − f esc )ξ ion versus UV magnitude and demonstrates the revolutionary capabilities of MUSE and JWST/NIRCam: we are able to spectroscopically confirm extremely UV-faint galaxies via their high Lyα EW and we are able to infer Hα, and therefore ξ ion , from much fainter individual galaxies than was previously possible with Spitzer, where stacking was necessary at M uv ∼ > − 20 (e.g., Lam et al. 2019;Maseda et al. 2020).We reach ∼ 1 dex lower than any previous studies at similar redshifts without the need of stacking methods.We can reach individual detections of very faint galaxies, M uv < −17.We also find results consistent with those at z ∼ 2 (Shivaei et al. 2018) and at z ∼ 4 − 5 for > L * galaxies (Bouwens et al. 2016) and < L * galaxies (Lam et al. 2019, where a stacking analysis was used), as shown in Figure 2. We note our observations demonstrate the large scatter in (1 − f esc )ξ ion at fixed M uv which was not possible to observe in previous analyses which used stacking of Spitzer photometry for UV-faint galaxies.
As described in Section 3.4 we perform a linear regression to assess correlations in our data.In contrast to Lam et al. (2019) we find significant evidence for a weak trend between ξ ion and M uv , where the highest ξ ion tends to come from the faintest galaxies.Since our sample is not M uv complete, we only study the correlation up to the peak of our M uv distribution (=-18.14) in Figure 1.We find log 10 [(1 − f esc )ξ ion ] = (0.03 ± 0.02)(M uv + 20) + 25.36 ± 0.03, but with large scatter (see Table 1).
Figure 4 shows that ξ ion follows a strong trend with Hα EW, as found in previous work (Harikane et al. 2018;Lam et al. 2019;Tang et al. 2019).Previous works were limited to only the highest Hα EW values, while we reach ∼ 0.75 dex lower due to the sensitivity of NIRCam.This trend is consistent with a picture where ξ ion is elevated in the youngest, most highly star-forming galaxies (e.g., Tang et al. 2019).We find log 10 [(1 − f esc )ξ ion ] = (0.73 ± 0.04)(log 10 EW Hα − 2.5) + 25.15 ± 0.02.The measurement by Maseda et al. (2020), obtained from a stack of extremely UV-faint galaxies with high Lyα EW, lies significantly above our sample and the rest of the literature, with higher (1 − f esc )ξ ion at fixed Hα EW.As discussed by Maseda et al. (2020) this likely implies their sources have a much lower gas-phase metallicity than other samples.
We also find that Lyα-selected galaxies have a higher Hα EW than the photometrically-selected sample, (median EW=732 ± 187 Å compared to 457 ± 161 Å for the photometric sample).A two-sample Kolmogorov-Smirnov test establishes that the EW distributions of the two samples are different (p-value 0.01).This is likely to be the primary driver of the increased ξ ion distribution for the Lyα-selected sample (Figure 2).
At fixed Hα EW we see a clear tendency for galaxies having very blue β UV slopes to have elevated ξ ion (Figure 4).This trend is also seen in the full sample -Figure 5, where we find high ξ ion is weakly correlated to blue β slope, but with large scatter.We find log 10 (1 − f esc )ξ ion = (−0.20 ± 0.04)(β + 2) + 25.41 ± 0.01 (see Table 1).Similar correlations have been seen at z∼6 (i.e.Ning et al. 2022).Using KS test we find no significant difference in the β distributions for the Lyα and photometric samples.Our sample has a median β = −2.1.

A search for high escape fraction and extremely low metallicity galaxies
As well as being a tracer for the ionizing photon production of galaxies, nebular emission lines are also sensitive to the escape fraction.Zackrisson et al. (2013) proposed that in galaxies with very high ionizing escape fraction, one would expect a reduction in nebular emission line strength (Hβ EW ∼ < 30 Å) and extremely blue UV slopes (β < −2.5) due to the lack of nebular continuum.Early JWST observations have discovered potentially very blue galaxies (Topping et al. 2022, though c.f. Cullen et al. 2022) and galaxies with weak nebular line emission yet high sSFR (via [OIII]+Hβ, Endsley et al. 2022), potentially indicating a population with high ionizing escape fraction.However, the observation of low [OIII]+Hβ line strengths could also be caused by very low gas-phase metallicity (decreasing the strength of [OIII] emission) or a recent turnoff in star formation (which would also decrease all nebular emission lines).
Given the redshift range of our sample we can infer both Hα and [OIII]+Hβ line strengths for 241 galaxies, allowing us to test these scenarios and to search for galaxies with high escape fraction.We obtained the [OIII]+Hβ nebular line fluxes as described in Section 3.1).
In Figure 6 we show UV β slopes as a function of intrinsic Hα EW for our sample (where we correct for dust attenuation as described in Section 3.1.We compare our sample to the region proposed by Zackrisson et al. (2017) to have f esc > 0.5.While several sources fall into this region, and also have low [OIII]+Hβ EW ( ∼ < 100 Å) the uncertainties are too large to make these robust candidates.We discuss this further in Section 5.2.
Figure 7 shows Hα EW as a function of [OIII]+Hβ for our sample.We see the expected positive correlation between both nebular emission lines, as these lines are all generated by the effects of stellar ionizing radiation.We see a very large scatter (with a range ∼ 1.5 dex) as expected due to variations in metallicity, temperature, and ionization parameter which will affect the strength of [OIII] individual galaxies (e.g., Maiolino   et al. 2008;Steidel et al. 2014;Sanders et al. 2021).We find log 10 EW(Hα) = 0.97 ± 0.06(log 10 EW([OIII] + Hβ) − 2.5) + 2.52 ± 0.03.
Galaxies with detected Lyα tend to occupy the top right of the plot, with strong nebular emission lines, suggesting these are young star-forming galaxies with low metallicity and large ionization parameters, producing copious ionizing photons needed to power these emission lines (see e.g., Yang et al. 2017;Du et al. 2020;Tang et al. 2021, for more detailed studies).We find the Lyα selected galaxies have stronger [OIII]+Hβ EW compared to the photometric population, following the trend with Hα EW in Figure 4. Though, as discussed by Tang et al. (2021), not all galaxies with strong nebular emission are detected in Lyα, indicating that Lyα transmission is reduced due to a high column density of neutral gas in these systems and/or inclination effects.We compare our data to a z ∼ 2 sample by Tang et al. (2019), which was selected based on strong [OIII] emission.We find a similar correlation, but overall our ratio of Hα EW/[OIII]+Hβ EW is higher by ∼0.1 dex.Given that the Tang et al. (2019) sample has significantly sub-solar gas-phase metallicity Z < 0.3Z (Tang et al. 2021), the decrease we observe in [OIII] at fixed Hα EW would likely imply an overall lower metallicity due to a lower number of metal atoms in our sample.

The profile of a strong ionizer
Thanks to the depth of JWST/NIRCam we have been able to assess trends of ξ ion at z > 3 across the broadest range of galaxy properties to-date.From these results, we corroborate previous work at lower redshift and high luminosities, but push the measurement of ξ ion to a large sample of individual UV-faint galaxies for the first time.
We find that galaxies with strong ionizing photon emission tend to have high Hα EW, low UV luminosity, blue UV β slope and Lyα emission -all implying that these galaxies are young, with likely low dust content and metallicity, and have a high O/B star population, capable of producing hard ionizing photons (e.g., Tang et al. 2019;Boyett et al. 2022b).This picture of the integrated emission from galaxies is complemented by high spatial resolution observations of highly magnified arcs with JWST.These have revealed extremely young star clusters ( ∼ < 10 Myr) with [OIII]+Hβ EW > 1000 Å, which dominate the ionizing photon production in their galaxy (Vanzella et al. 2022a,b), indicating that there can be large variations in ξ ion in individual galaxies, if they contain multiple stellar populations, but that the variation is primarily driven by the age of the stellar populations.Ly-Detected no Ly-Detected Fig. 5. UV β slope vs ξ ion (1-f esc ).In purple we show Lyα detected galaxies.In gray is the photometrically selected sample with no Lyα detected.As above, error bars are only shown for 30% of the sources for clarity.We add the stacked measurements from Lam et al. (2019) for comparison.We find a very weak trend of increasing ξ ion with decreasing β, with a linear slope of −0.10 ± 0.06.We also find that overall, our Lyα galaxy sample has higher ξ ion than the photometrically selected one; the primary reason for this difference is that the former has higher Hα EW (Figure 4).The enhanced prevalence of Lyα emission in strong Hα emitters is likely a combination of increased production of Lyα photons due to the young stellar population implied by strong Hα, and potentially also an increase in the Lyα escape fraction in the interstellar medium (Tang et al. 2021;Naidu et al. 2022).In these rapidly star forming galaxies, the hard ionizing radiation may be ionizing the ISM and/or feedback may disrupt the ISM gas leading to a reduced HI column density and dust cover.We note that the galaxies with the highest (1 − f esc )ξ ion are not necessarily all Lyα-emitters, likely due to variance in the geometry and column density of neutral gas and dust in these sources.(Ning et al. 2022) has shown this same correlation between ξ ion and Lyα for a broad range of luminosities and equivalent widths.

The ionizing photon escape fraction
In Section 4.2 we explored whether our sample shows signs of high ionizing photon escape fraction, f esc , using the low Hβ EW -blue UV β slope region defined by Zackrisson et al. (2017) for f esc > 0.5.While several sources fell into this region, with both low Hα EW and [OIII]+Hβ EW ( ∼ < 100 Å), the uncertainties on the line flux measurements are too large for these to be robust candidates.More precise emission line measurements with JWST spectroscopy will be vital for identifying such candidates and their relative abundance in the galaxy population.
The lack of high f esc candidates amongst the Lyα-selected galaxies is also surprising.As the same conditions (low neutral gas covering fraction) facilitate both Lyα escape and Lyman continuum escape, a correlation between Lyα and Lyman continuum escape fraction is expected (e.g., Verhamme et al. 2015;Dijkstra et al. 2016;Reddy et al. 2016).
As discussed by Topping et al. ( 2022), however, it is possible for galaxies with high f esc but with very young ages to still have high nebular emission due to high ionizing photon production.It is likely that the criteria proposed by Zackrisson et al. (2017) can only find high f esc systems within the bounds of the assumptions made for their model, such as galaxy SFH, ages, metallicities, dust, but also the stellar models used.Our results suggest the low luminosity galaxies with high sSFR but low [OIII]+Hβ EW observed by Endsley et al. (2022) may be more likely due to variation in metallicity than high f esc .

Conclusions
We have inferred the hydrogen ionizing photon production rate, modulo the escape fraction, in the largest sample of individual sub-L * z > 3 galaxies to-date, spanning −23 ∼ < M uv ∼ < − 15.5, with a median M uv = −18.1,thanks to deep JWST/NIRCam imaging, enabling us to track the demographics of the ionizing population.Our conclusions are as follows: 1.The median log 10 (1 − f esc )ξ ion of our sample is 25.33 ± 0.47 with an intrinsic scatter of 0.42 dex.The inferred ξ ion distribution of our sample has values in a range of ∼ 1.5 dex, implying a wide range of galaxy properties and ages.2. We find significant trends of increasing (1 − f esc )ξ ion with increasing Hα EW, decreasing UV luminosity, and with decreasing UV slope, all suggesting galaxies which are most efficient at producing ionizing photons are young, highly star forming, which are normally expected to be low metallicity and dust-poor.3. We find galaxies selected with strong Lyα emission to have higher ξ ion than photometrically-selected galaxies, with median log 10 (1− f esc )ξ ion values of 25.39±0.64 and 25.31±0.43respectively.We find the Lyα-detected galaxies have an elevated Hα EW distribution, thus the increased ξ ion is likely driven by the selection based on Lyα selecting a younger population.As strong Lyα emitters also likely have high ionizing photon escape fractions, this implies the intrinsic production rate of ionizing photons in these galaxies could be  2019), which was selected to have strong [OIII] EW, implying that we might be observing lower metallicity galaxies.
significantly higher than what we can infer from Hα luminosities.4. We examine our sample for signs of very high f esc by comparing the inferred strengths of nebular emission lines ([OIII]+Hβ and Hα) and the strength of the nebular continuum via the UV β slope.We find no significant evidence for sources with high escape fraction galaxies with low nebular emission line strength and very blue UV β slopes.The reduced strength of [OIII]+Hβ EW in our z > 3 sample compared to a sample at z ∼ 2 from Tang et al. (2019) implies our sample has likely lower gas-phase metallicity and/or ionization parameter.
We have demonstrated the power of JWST/NIRCam photometry to more precisely constrain the rest-frame optical emission of UV-faint high redshift galaxies than previously possible with Spitzer/IRAC.These observations allow us to constrain the production rate of ionizing photons from early galaxies, corroborating the picture obtained from previous stacking analyses that ξ ion is elevated in young, highly star forming galaxies, but that there is a broad distribution of ξ ion , likely driven by variation in galaxy properties and ages.
With JWST spectroscopy it is becoming possible to obtain direct measurements of optical emission lines in large samples (e.g., Sun et al. 2022;Williams et al. 2022;Matthee et al. 2022).Deriving a census of the ionizing photon production rate across the full galaxy population will be necessary to fully understand reionization.Here we have shown that ξ ion is elevated in UV-faint galaxies with strong nebular emission lines, likely due to young ages.While a thorough analysis of the implications of our results for reionization are beyond the scope of this work, this becomes more prominent at high redshift (e.g., Boyett et al. 2022b;Endsley et al. 2022), implying that it would be possible to complete reionization with modest f esc .Considering the full distributions of ξ ion and f esc across galaxy properties will be required to assess the primary drivers of reionization.
Fig.1.Galaxies studied in this work; in purple Lyα detected galaxies, and in gray galaxies photometrically selected (with no Lyα detected).Top: Distribution of redshifts for the spectroscopic and photometric samples.We show the spectroscopic redshift where available, or the peak photometric redshift.Bottom: UV magnitude distribution for our sample, we find a median value of −18.14 ± 1.58, with no statistically significant difference between both samples.

Fig. 3 .
Fig. 3. M uv vs (1f esc )ξ ion .In purple circles we show Lyα detected galaxies.In gray circles are the photometrically selected sample with no Lyα detected.We show data from Maseda et al. (2020); Harikane et al. (2018); Lam et al. (2019) and Bouwens et al. (2016) as colored boxes for comparison.We find evidence for an increase in log 10 (1 − f esc )ξ ion towards fainter UV magnitude, with a slope of 0.03 ± 0.02, only considering the range where our sample is M uv complete (M uv < −18.1).We show literature constraints at similar redshifts as colored shapes(Bouwens et al. 2016;Harikane et al. 2018;Lam et al. 2019;Maseda et al. 2020), noting that all constraints fainter than M uv ∼ > − 20 were obtained by stacking IRAC photometry.

Fig. 4 .
Fig.4.We compare Hα equivalent width with the ionizing photon production that does not escape the galaxy.As stars we show Lyα detected galaxies, as circles are photometrically selected galaxies with no Lyα.As above, error bars are only shown for 30% of the sources for clarity.We color code these two samples by UV β slope.On top we show the distribution of Hα EW for the same two samples, compared to the values found byTang et al. (2019).We add data fromHarikane et al. (2018) andLam et al. (2019) for comparison, which is at the high end of our observed Hα EW distribution.We see that stronger ξ ion very strongly correlates with Hα EW.Galaxies with detected Lyα emission have an Hα EW distribution with higher values, median 732 ± 187 Å compared to 457 ± 161 with a Kolmogorov-Smirnov test p-value 0.01.The sources with the reddest UV slopes lie systematically below the best-fit relation at fixed Hα EW.

Fig. 6 .
Fig. 6.Comparison between intrinsic (unattenuated) EW of Hα and UV β slope, color-coded by [OIII]+Hβ EW.Lyα galaxies are shown with star shaped markers, and the photometric sample as circles.Galaxies shown in gray, do not have [OIII]+Hβ EW measurements.We show the region predicted by Zackrisson et al. (2017) to show f esc > 0.5.We rescale from Hβ EW to intrinsic Hα with a case B recombination scenario of factor 2.89, assuming a flat optical continuum in f λ , which we confirm from our SED fitting done in 3.1.

Fig. 7 .
Fig. 7. Comparison between EW of Hα and [OIII]+Hβ.We show Lyα detected galaxies as stars and the photometrically selected sample with no Lyα detected as circles.On top we show the distribution of [OIII]+Hβ EW for both of our samples.We find a very strong correlation between Hα EW and [OIII]+Hβ EW, though with large scatter.The dashed lines are the correlation trends found for this work (red) andTang et al. (2019) (green).We find higher Hα EW/[OIII]+Hβ EW than the z ∼ 2 sample fromTang et al. (2019), which was selected to have strong [OIII] EW, implying that we might be observing lower metallicity galaxies.

Table 1 .
Linear fitting parameters for trends with log 10 (1 − f esc )ξ ion Notes.We fit for log 10 [(1 − f esc )ξ ion ] = α + βX + , where is the intrinsic scatter which is assumed to be normally distributed with variance σ 2 .