An empirical reionization history model inferred from the low-redshift Lyman continuum survey and the star-forming galaxies at z > 8

We present a new analysis of the rest-frame ultraviolet (UV) and optical spectra of a sample of three z > 8 galaxies disco v ered behind the gravitational lensing cluster RX J2129.4 + 0009. We combine these observations with z > 7.5 galaxies from the literature, for which similar measurements are available. As already pointed out in other studies, the high [O III ] λ 5007/[O II ] λ 3727 ratios ( O 32 ) and steep UV continuum slopes ( β ) are consistent with the values observed for low-redshift Lyman continuum emitters, suggesting that such galaxies contribute to the ionizing budget of the intergalactic medium. We construct a logistic regression model to estimate the probability of a galaxy being a Lyman continuum emitter based on the measured M UV , β , and O 32 . Using this probability and the UV luminosity function, we construct an empirical model that estimates the contribution of high-redshift galaxies to reionization. The preferred scenario in our analysis shows that at z ∼ 8, the average escape fraction of the galaxy population [i.e. including both LyC emitters (LCEs) and non-emitters] varies with M UV , with intermediate UV luminosity ( − 19 < M UV < − 16) galaxies having larger escape fraction. Galaxies with faint UV luminosity ( − 16 < M UV < − 13.5) contribute most of the ionizing photons. The relative contribution of faint versus bright galaxies depends on redshift, with the intermediate UV galaxies becoming more important o v er time. UV bright galaxies, although more likely to be LCEs at a given log( O 32 ) and β , contribute the least of the total ionizing photon budget.


INTRODUCTION
Reionization occurs rapidly at redshift  ≳ 6 (Becker et al. 2001;Fan et al. 2006;Bañados et al. 2018;Eilers et al. 2018;Planck Collaboration et al. 2020;Becker et al. 2021), where most of the neutral hydrogen in the intergalactic medium (IGM) was ionized by the first sources of Lyman continuum ( < 912 Å) photons.There is substantial evidence that star-forming galaxies are the dominant source of ionizing radiation, since the number density of quasars significantly decreases at  > 3 (Matsuoka et al. 2018;Kulkarni et al. 2019;Jiang et al. 2022;Schindler et al. 2023).Some uncertainties still remain about the relative role of bright versus faint galaxies (Robertson et al. 2015;Sharma et al. 2016;Finkelstein et al. 2019;Naidu et al. 2020), and on the contribution of low luminosity AGN to the overall ionizing photon budget (Giallongo et al. 2019;Jiang et al. 2022;Lu et al. 2022).
The contribution of star-forming galaxies to the reionization of the intergalactic medium (IGM) is often parameterized as a function of three components   =       , where   is the production rate of ionizing photon emitted to the IGM,   is the escape fraction of Lyman continuum (LyC) photons,   is the ionizing photon production efficiency, and   is the UV luminosity density measured at rest-frame 1500 Å (Robertson 2022).
The UV luminosity density is the integral of the UV luminosity function down to a given magnitude limit.The luminosity function (LF), often described with a Schechter function (Schechter 1976), can be directly measured at all redshifts.Before JWST, Hubble and ground based programs provided measurements of the redshift evolution of the LFs out to  ∼ 8 (Bouwens et al. 2015;Livermore et al. 2017;Atek et al. 2018;Bouwens et al. 2021).Its shape, when observations are sufficiently deep to cover the knee of the function, does not deviate substantially from a Schechter LF, although a double power-law parameterization has been suggested for the highest redshift samples at  ∼ 6 (Bowler et al. 2015;Khusanova et al. 2020;Harikane et al. 2022;Donnan et al. 2023).Results before the first JWST data were obtained indicated that luminous galaxies are progressively less numerous toward high-redshifts.For example, luminous galaxies (  < −19) at  ∼ 8 are 25 times less numerous compared to those at  ∼ 0 (Bouwens et al. 2021).The early   results have discovered more bright galaxies at the early universe at  > 8 than we expected (Finkelstein et al. 2022).A double powerlaw LF (e.g., Harikane et al. 2022) is proposed to include the newly discovered population of bright galaxies, since the double powerlaw (DP) LF are higher at the bright UV end (  < −22), and decreases slower toward high redshift at  > 8 than the LFs proposed pre-  (e.g., Bouwens et al. 2021).
is the ionizing production rate relative to the UV luminosity density at 1500 Å.   depends on the shape of the galaxy's ionizing spectrum, which in turns depends on the initial mass function (IMF), the stellar metallicity, and the fraction of binary stars.During the EoR, galaxies are believed to have a "top-heavy" IMF (Davé 2008;van Dokkum 2008;Sharda & Krumholz 2022), lower metallicity, and low dust extinction, where all these properties lead to a higher   (Chisholm et al. 2019;Atek et al. 2022).Indeed, the average   is observed to increase toward higher redshifts (Matthee et al. 2017;Shivaei et al. 2018;Atek et al. 2022;Matthee et al. 2022a).Models of the reionization history of the universe (e.g., Robertson et al. 2013) typically assume a range of log(  ) between 25.2 and 25.3, depending on the specifics of the stellar population models that are used (e.g., Leitherer & Heckman 1995).This range is often referred to as the "canonical" log(  ) (e.g., Shivaei et al. 2018).Some studies find that the average log(  ) evolves with redshift (Matthee et al. 2017).The cause of this evolution is typically ascribed to the evolution of the intrinsic properties of galaxies (e.g., the lower stellar metallicity).The mode of star-formation at any given time, however, can also be important.This is because the definition of   is linked to the lifetimes of massive stars (Kennicutt & Evans 2012).Being produced by recombination in H ii regions, the H luminosity traces the presence of the massive stars with short lifetimes (tens of Myrs), while the UV continuum, can also be generated from intermediate mass stars, which have longer lifetimes.Accordingly, if galaxies were, on average, characterized by shorter and more frequent bursts of starformation at higher redshifts (for example because their dark-matter halo is still growing and feedback could have a stronger impact than in galaxies at low redshift) we would expect a different   distribution, even for the same physical properties such as metallicity and IMF.
In the calculation of   ,   is the term that is least constrained by observations, and the only one that cannot be measured directly for galaxies seen during the EoR.At redshifts  ≳ 4, even if ionizing photons were escaping from galaxies in a proportion similar to their low redshift counterparts, those photons would be absorbed by the neutral hydrogen in the IGM (e.g., Worseck et al. 2014).Therefore, during the EoR,   can only be estimated using empirical indirect indicators calibrated from low-redshift LyC emitters (Flury et al. 2022a).Direct observation of the intrinsically faint LyC is challenging.Until recently, only a few dozens of galaxies at  < 0.4 were spectroscopically confirmed as LyC emitters (   > 0) (Bergvall et al. 2006;Leitet et al. 2011;Borthakur et al. 2014;Leitherer et al. 2016;Izotov et al. 2016bIzotov et al. ,a, 2018a;;Flury et al. 2022a), and even fewer LyC emitters at 2 <  < 4 have been identified (Vanzella et al. 2016;de Barros et al. 2016;Shapley et al. 2016;Bian et al. 2017;Rivera-Thorsen et al. 2019;Fletcher et al. 2019;Marques-Chaves et al. 2021;Saldana-Lopez et al. 2023).Using indirect probes of LyC escape fraction, such as high values of O 32 , steep UV continuum slopes, and high star formation rate surface densities, various authors have been successful in identifying a population of LyC emitting galaxies (Chisholm et al. 2018;Naidu et al. 2020;Flury et al. 2022a;Saldana-Lopez et al. 2022;Chisholm et al. 2022).These indicators, however, may not directly provide the value of   , as the escape is a complex process that depends on the neutral gas density and covering fraction (Bassett et al. 2019).
In this paper, we present a new analysis of a sample of gravitationally lensed galaxies at  > 8 discovered in the RX J2129 galaxy-cluster field (Langeroodi et al. 2022;Williams et al. 2022).We refer to these three galaxies as the RXJ2129 high- galaxies.Using these galaxies and inference based on the analysis of the low-redshift LyC survey by Flury et al. (2022a), we present a new empirical model for the galaxy contribution to the reionization history.The structure of the paper is as follows.In Section 2, we describe the observations and our measurements.We compare our objects with the low redshift Lyman continuum emitters in Section 3. In Section 4, we present our estimation of reionization history and conclude in Section 5. Throughout the paper, we denote   as the singular and plural absolute Lyman continuum "escape fraction" and "escape fractions", and we assume a ΛCDM cosmology with  0 = 67.66km s −1 Mpc −1 , Ω  = 0.04897, and   =8.5988×10 −30 g cm −3 (Planck Collaboration et al. 2020).

OBSERVATIONS AND ANALYSIS
The details of the observations and the data reduction are reported in the companion papers Williams et al. (2022) and Langeroodi et al. (2022).Briefly, we obtained imaging of the RXJ2129 cluster field with the JWST NIRCam instrument in the F115W, F150W, F200W, F277W, F356W, and F444W filters as part of the Director's Discretionary program (DD-2767; PI: P. Kelly) to obtain followup spectroscopy of the strongly lensed SN 2022riv (Kelly et al. 2022).We identified three high-redshift galaxy candidates using the EAZY (Brammer et al. 2008) photometric redshift estimation algorithm.The follow-up spectroscopy of the RXJ2129 cluster field was obtained using the NIRSpec instrument on JWST in Multi-Object Spectroscopy (MOS) mode.The spectrum wavelength covers from 0.6m to 5.3m, and the spectral resolution ranges from  ≈ 50 on the blue end to  ≈ 400 on the red end.
The flux calibration for the NIRSpec spectroscopy was performed in the PHOTOM step of the JWST Spec2Pipeline 1 .Aperture corrections were applied to the NIRSpec data in the PATHLOSS step of the Spec2Pipeline.This step calculated the expected slitlosses for a point source in a given position within the shutter.Since our sources are so small with half light radius   ≃ 0.04 ± 0.01 arcsecond (Williams et al. 2022), we did not apply any additional aperture correction.
In Figure 1, we show the rest-frame UV spectra with associated errors of RXJ2129-ID 11027 at redshift =9.51, and RXJ2129-ID 11002 at =8.16.The rest-frame UV spectrum of RXJ2129-ID 11022 falls outside of the spectral range covered by the detector.Note that in the RXJ2129-ID 11002 spectrum, the masked peak at 1125 Å is caused by a cosmic-ray hit rather than Ly emission.

Ultraviolet Properties
In this section we describe how we derive the galaxy properties which will be needed in Section 4 to compute the average ionizing background (i.e., to compute   ), namely the 1500Å absolute UV magnitude (  ), the slope of the spectral continuum, , the escape fraction of ionizing radiation,   , and the ionizing photon production efficiency,   .Other properties, such as measurements of rest-frame optical emission lines, are presented in the companion paper on the mass-metallicity relation (Langeroodi et al. 2022).

𝑀 𝑈𝑉 and the slope of the stellar continuum, 𝛽
To measure the   we use the rest-frame UV spectrum, when available.Specifically, we measure   by averaging the flux density between 1400 Å and 1600 Å.This is a spectral region typically free from strong emission lines.For galaxy ID 11022, we calculate the   using the flux density computed in the F150W filter, which corresponds to a rest-frame wavelength of 1650 Å.The analysis of the spectral energy distributions presented in Williams et al. (2022) and Langeroodi et al. (2022) suggests that they suffer low level of dust attenuation.To calculate the observed UV continuum slopes we do not apply any dust correction to the observed magnitudes.The color excess E(B−V) is then estimated with the observed UV continuum slopes (Chisholm et al. 2022).
Since the spectrum of of galaxy ID 11027 suffers from the contamination (Figure 1), we measure the UV continuum slopes, , defined as   ∝   , by fitting a power-law function to the observed flux densities in the photometric bands F1500W, F200W, and F277W.We used Markov Chain Monte Carlo sampling to sample the posterior on the parameter .
The observed UV continuum slopes  ≃ −2 suggest a certain amount of dust attenuation level, with E(B−V)   ≃ 0.03 (see Chisholm et al. 2022, Section 4 and Appendix A).We adopt the Small Magellanic Cloud-like dust attenuation law from Gordon et al. (2016) with R  ≡ A  /E(B−V) = 2.74 to calculate the dust extinction.

Escape Fraction 𝑓 𝑒𝑠𝑐
The escape fraction of ionizing radiation cannot be directly measured at  ≳ 4 because the optical depth for ionizing photons is high (e.g., Worseck et al. 2014).Accordingly, we need to use indirect estimators of   , calibrated at lower redshifts using galaxies as close as possible to those that we observe at  ≳ 6, during the EoR.Schaerer et al. (2022) demonstrated that local extreme emission line galaxies have similar properties to the galaxies JWST has been uncovering.
A number of studies have measured   in the local universe and found relations between   and galaxy properties.Despite a number of   indicators (e.g., O32, ) have successfully been identified using low-redshift LyC emitters, the scatter between the value of   and the value of its indirect estimator is observed to be substantial (Flury et al. 2022b).These indicators, therefore, may not be sufficient for an accurate estimate of   due the complicated escape process involving gas density and covering fraction (Bassett et al. 2019).Here, we use the relation between   and  1500 , the UV slope between 1300Å to 1800 Å, discussed in Chisholm et al. (2022).The relation is given as follows: (2) We report the UV properties of the RXJ2129 high- galaxies in Table 1.For the three galaxies in the SMACS0723 sample, we adopt the magnification, UV absolute magnitudes, and photometric UV continuum slopes reported in Schaerer et al. (2022).Since the H flux is not reported in Schaerer et al. (2022), we adopt the metallicity, emission line fluxes, and ratios from Curti et al. (2023).For completeness, we report the properties of the SMACS0723 high- galaxies in Table 2.

COMPARISON WITH LOW REDSHIFT ANALOGS
We compare the properties of the  > 7.5 galaxies, including the RXJ2129 high- galaxies, the SMACS0723 high- galaxies, and the galaxies recently reported in the GLASS-JWST program (Mascia et al. 2023) and the CEERS survey (Tang et al. 2023) to those of the low-redshift galaxies studied in the Low- Lyman Continuum survey  Observed, including lensing (Schaerer et al. 2022). () Corrected for lensing (Schaerer et al. 2022), we adopted a 0.2 magnitude uncertainty throughout the calculation.
(LzLCs, at  ∼ 0.3 Flury et al. 2022a).The LzLCs includes 66 newly observed low redshift ( < 0.4) galaxies and 23 galaxies ( < 0.46) in the literature (Izotov et al. 2016b(Izotov et al. ,a, 2018a(Izotov et al. ,b, 2021;;Wang et al. 2019), for which the LyC emitters (LCE) are defined as galaxies with Lyman Continuum detected with 97.725% confidence.We refer to these 89 galaxies as the LzLCs sample.In this sample, 50 galaxies are confirmed as LyC emitters (LCE), and 39 galaxies are not detected in the LyC (non-LCE).The LzLC objects were selected to span a broad range in physical properties associated with a large probability of high escape fraction of ionizing radiation.Specifically, galaxies were selected to have high O 32 ratio, steep UV continuum slope, and high star formation rate surface density, Σ   .

The probability of a galaxy being a LyC emitter
We combine the O 32 ratio and the slope of the UV continuum in a combined indicator,  32 =log(O 32 )−.We choose , O 32 , and   in our model for their simple accessibility in the observations.In Figure 2, we show how, in a diagram of  32 as a function of the absolute UV magnitude, LyC emitters (LCE, open circles) are efficiently separated from non-LyC emitters (non-LCE, filled black circles).We also explored whether a correlation exists between  32 and   , but we do not observe any simple relationship between these quantities, strengthening the conclusion that indirect estimators are mostly useful to identify LyC emitting galaxies, rather than estimating the value of   .Accordingly, we apply a logistic regression to the LzLCs sample, to estimate the probability that a galaxy is a LCE based on   and  32 .The logistic discriminator can be written as: where   is the probability of the galaxy being a LyC emitter, and the best fit values for ( 0 ,  1 ,  2 ) are (−25.82,−1.09, 1.72).
The discriminator is shown in Figure 2 with the black dashed line, in the range where no extrapolation is necessary.Assuming the discriminator has no evolution on redshift, we find that RXJ2129 ID 11002, SMACS 06355, and SMACS 10612 have P  > 0.8, suggesting that these galaxies are likely LyC emitters.SMACS 04590 has higher O 32 and steeper  1500 than SMACS 06355, but its position in Figure 2 (most right yellow point at   ≈ −18) suggests that SMACS 04590 is less likely to be a LCE due to the fainter UV magnitude, with   = 0.39.We note, however, that the local LzLC sample does not extend to magnitudes fainter than ≈ −18.5, such that it is not clear where is the boundary between LyC and non-LyC emitters at the faint end of the UV magnitude > −19, and our logistic discriminator may not be applicable at the faint UV side.The physical connections between the   indicators and   is likely the changes of metallicity, dusts, and star formation.
The faint star-forming galaxies (associated with the low metallicity galaxies) have noticeably different properties compared with their bright counterpart.For instance, the intrinsic UV continuum slope does not get appreciably bluer below 10% solar metallicity (Bouwens et al. 2010;Topping et al. 2022).The changes in O 32 are also less obvious below 10% solar metallicity (Curti et al. 2017;Sanders et al. 2021).Therefore the  and  32 may become less effective as indicators for LyC leaking at the faint UV end.In this paper, we apply the logistic discriminator to   <-18, where the fit is not constrained.To account for the uncertain behavior at the faint end we consider three scenarios, as illustrated in Figure 2: The "extrapolation scenario", represented by the blue dashed line, involves a simple linear extrapolation of the regression to faint magnitudes.In this scenario, faint galaxies are less likely to be classified as Lyman Continuum Emitters (LCEs), as being an LCE requires a large value for ( 32 ) − .The "restrained scenario", denoted by the green dashed line, exhibits a shallower slope compared to the extrapolation scenario.In this case, the criteria for classifying faint galaxies as LCEs are less strict.Finally, the discriminator plateaus in the "flatten scenario", indicated by the red dashed line, where faint galaxies have a significantly higher chance of being categorized as LCEs.The parameters ( 0 ,  1 ,  2 ) in equation 3 for the restrained scenario are (−15.83,−0.55, 1.72), while for the flatten scenario, the values are (−5.655,0, 1.72).
The choice of the discriminator's turning point is primarily determined by observations extending down to   = −18.5.The slope of the extended discriminator fundamentally governs the contribution of the faint galaxy population to reionization.We will demonstrate in Section 4.2 that observational constraints are more in agreement with the reionization history derived from the restrained scenario.

The dependency of 𝜉 𝑖𝑜𝑛 on 𝑀 𝑈𝑉 and redshift
Many works have studied the dependence of   , and whether   depends on the UV magnitudes is still under debate (Matthee et al. 2017;Shivaei et al. 2018;Emami et al. 2020;Nakajima et al. 2020;Atek et al. 2022).We calculate the  0  of the LzLCs sample using equation 1.We calculate the  0  of the  > 7.5 galaxies in the GLASS-JWST program (Mascia et al. 2023) using equation 1 with H luminosities derived from the star formation rate density.We show the log( 0  ) as a function of the UV magnitude in the lower panel of Figure 3, where we observe a positive correlation between  0  and the UV magnitude for both the LzLCs sample and the high-redshift galaxies.However, this  0  vs. UV correlation may be a secondary dependence and the  0  is in fact dependent on other factors such as metallicity and burstiness.In Figure 4, we show all high-redshift galaxies have log( 0  ) ≥ 25.2, consistent with the extrapolation at  = 8 of the observed redshift evolution (Stark et al. 2015;Bouwens et al. 2016;Nakajima et al. 2016;Stark et al.

REIONIZATION HISTORY
In this section, we use the insights gained in the previous sections to build a new model of the reionization history.Specifically, in the calculation of the redshift evolution of the neutral fraction ( H i ) we consider (1) the empirical constraints on the probability that a galaxy is a LCE, (2) the dependency of   on , and (3) the dependency of log(  ) on   or redshift .2017), and the 1- limit.

The empirical model
We calculate the neutral fraction  H i by solving where   is the ionizing photon production rate,  H is the comoving gas number density, and  rec is the recombination time scale (Madau et al. 1999;Robertson et al. 2013;Ishigaki et al. 2018). H and  rec are defined as where   = 0.76,   = 0.24 are the primordial mass fraction of hydrogen and helium (Planck Collaboration et al. 2020), Ω  is the baryon energy density fraction,   is the critical density,  H ii ≡ ⟨ H ii2 ⟩/⟨ H ii ⟩ 2 is the clumping factor, and   () is the case B recombination coefficient.Here we assume  H ii = 3, and   = 2.6 × 10 −13 cm3 s −1 , for an electron temperature of 10 4 K.We solve equation 4 iteratively, assuming the boundary condition that  H i = 1 at  = 15, i.e., we assume that the first sources of ionizing photons appear at this redshift. 2 The ionizing photon production rate depends on the escape fraction of ionizing radiation, the ionizing photon production efficiency and the volume density of UV luminosity (  =       ).In what follows we explain how we include the empirical constraints derived in the previous sections in the calculation of  at each time step.
First, we derive the population average escape fraction ⟨   ⟩ as a function of the UV absolute magnitude, accounting for the dependency of   on  and on the probability of each galaxy being a LCE.In each bin 3 of   , we simulate 100 galaxies, each with a different value of  and log(O 32 ), drawn from the distributions described below.These values are used, together with   , to compute the probability P  (  ,  32 ), using Equation 3. We then draw a Boolean value based on P  to determine whether or not the galaxy is a LCE.If a galaxy is a non-LCE, we set its    = 0.If a galaxy is a LCE, we draw its   value from the multivariate normal distribution with parameters described by Equation 2. The  values are drawn from observations of  > 8 galaxies, assuming the following normal distribution:  = (−0.17)× −5.4 ± 0.4, truncated at  = −3.5 (Cullen et al. 2022).We maintain a fixed slope, /  , in this relation and sample  assuming an intrinsic scatter of 0.4 dex (Rogers et al. 2014;Cullen et al. 2022).Properly accounting for the scatter in  is crucial in our model because a larger scatter results in galaxies with lower  values, which, on average, increase both the (P  ) and the escape fraction.Paalvast et al. (2018) reported no significant correlation between stellar mass and star formation rate with the O 32 values, therefore, we assume that log(O 32 ) values are independent of   .The log(O 32 ) values are generated from a normal distribution of 0.5±0.1 (Sanders et al. 2023).Finally, we compute ⟨   ⟩= 1 100 100 =1    as the average   of the 100 galaxies at each   bin, and we repeat this calculation 500 times.
The results of our simulations for the three scenarios are presented in Figure 5.In the top panel, we display in red the average P  and in green the average   computed without accounting for P  , both as functions of   .The average   without accounting for P  is identical in all three different scenarios since it depends solely on the sampled  values generated from the same distribution.
All three scenarios share the same discriminator for the brighter end (  < −18.5).Therefore, the results for P  and ⟨   ⟩ for   < −18.5 are consistent: brighter galaxies are more likely to be LCEs than fainter ones.However, since brighter galaxies have redder UV continua (larger  values), the   values are generally lower.These trends explain the results in the bottom panel, where we depict the average escape fraction of the population that now includes P  as a function of   .For brighter galaxies (  < −18.5), ⟨   ⟩ increases with decreasing   .
For fainter galaxies, both P  and ⟨   ⟩ change in response to the discriminator's slope.In the extrapolation scenario, most galaxies are less likely to be classified as LCEs.Consequently, the ⟨   ⟩ peaks at   = −18.5 and decreases as   becomes fainter.This ⟨   ⟩ vs.   pattern in our model behaves similarly to the ⟨   ⟩ inferred from the Ly emitters in Matthee et al. (2022b).In the restrained scenario, the slope of the discriminator is lower, resulting in a larger population of faint galaxies being classified as LCEs.In this case, the ⟨   ⟩ continues to increase until   = −16, and the   evolution is less pronounced compared to the previous scenario.In the flatten scenario, the P  reverses at the faint UV end, allowing a significant fraction of faint galaxies being LCEs.As a result, the ⟨   ⟩ keeps increasing with decreasing UV magnitude.
To demonstrate the redshift evolution of   , we also model the   ,   , and ⟨   ⟩ of the galaxies at redshift  = 4 (dashed lines).The  vs.   dependency is adopted from Bouwens et al. (2014) (see also Chisholm et al. 2022), and includes the intrinsic scatter Δ ∼ 0.4, used to draw samples of  for individual galaxies.In our model, as we assume no evolution on the logistic discriminators, the   does not evolve much with redshift.The   and ⟨   ⟩ decrease at lower redshifts as the UV continuum slope  becomes flatter/redder.Our predicted value of ⟨   ⟩≃ 0.03 at  = 4 aligns with the average of   =0.03±0.02 in Saldana-Lopez et al. (2023), where the   is derived using the depth of absorbtion lines from the low-ionization metals and the UV attenuation.Our result is lower than the average of   =0.06±0.01observed with the Lyman Break galaxies at  ∼ 3 in Pahl et al. (2021) and   =0.07±0.02 at  ∼ 3.5 in Begley et al. (2022).The   in Pahl et al. ( 2021) are also higher than the values calculated with equation 2 using the photometric UV slope  of the sample (see Saldana-Lopez et al. (2023) Figure 15).
The   difference can be caused by the choice of dust attenuation law, the intrinsic  value in the stellar population model, or a potential redshift evolution (see Saldana-Lopez et al. 2023, for more discussion).Another possible explanation of discrepancy with some of the  ∼ 3 observations is that we may have underestimated the intrinsic scatter on the  distribution.A larger intrinsic scatter of the  distribution would result in a greater number of individual galaxies with higher   , leading to a higher ⟨   ⟩.
We note that Sanders et al. (2023) found on average log( 32 )=0.9±0.1 at 6.5 <  < 9.3 significantly higher than log( 32 )=0.5±0.1 at  ∼ 5.6.It is unclear whether such high log( 32 ) are commonly seen at  > 6, while the  ∼ 5.6 sample has similar log( 32 ) to the 2.7 <  < 5 galaxies.Therefore we assume the log( 32 ) is independent with redshift and choose to sample log( 32 ) with the mean of 0.5.The value of log( 32 ) does not affect individual   in our model, and only slightly changes the fraction of LCEs population at   < −18.As a result, higher log( 32 ) values do not significantly increase the ⟨   ⟩ at   < −18.
In both the extrapolation and restrained scenarios, the ⟨   ⟩ peaks up at intermediate UV magnitudes (−18.5 <  < −16) and decreases toward the faint UV magnitude end.Grazian et al. (2017); Griffiths et al. (2022) have established an upper limit   < 0.04 at   < −20 4 .Although there are no sufficient observations of faint galaxies yet to study this trend, some simulations have predicted a similar pattern, with   peaking at intermediate UV magnitudes, rather than at the brightest or faintest end.For example, the medium   of SPHINX galaxies in Rosdahl et al. (2022) has the highest   at   ∼ −16, while Ma et al. (2020) predict that the highest   should be in 10 8 stellar mass galaxies, or   ∼ −19 (Stefanon et al. 2021).In these studies, the decline in   for larger and brighter objects is attributed to dust attenuation, while the decreasing for lower mass and fainter objects is attributed to inefficient star formation and increased susceptibility to stellar feedback.For log(  ), we propose two models: the   (UV) model and the   () model.In the   (UV) model, we assume  0  is linearly dependent on   .The linear relation is obtained with the minimum  2 fitting of all galaxies at  > 7.5 in Figure 3:  0  (UV) = 0.11(  + 20.0) + 25.46, and flatten at  0  = 24.5 and 26.We exclude the brightest object in our fitting, as it may potentially be powered by an AGN (Mainali et al. 2018;Tang et al. 2023).In the second model, we adopt   () as a function of redshift.We use the 1- upper limit redshift dependence in Matthee et al. (2017): log(  )() = 24.493+1.180log(1+), which is consistent with more  > 6 objects in Figure 4.
To explore the impact of the newly discovered bright galaxy populations at  > 8 to reionization, we adopt two luminosity functions at  > 8 into our models: a Schechter LF (Bouwens et al. 2021), and a Double Power-law (DP) LF (Harikane et al. 2022).We find the IGM neutral fractions at  > 8 derived the DP LF are only 1% lower than that of the Schechter LF (see Figure 6).Therefore, throughout this paper, we adopt the Schechter luminosity function described in Bouwens et al. (2021) We integrate the UV magnitude from −23 to −13.5.The UV magnitude is truncated at −13.5 to match previous studies (Livermore et al. 2017;Ishigaki et al. 2018;Atek et al. 2018;Naidu et al. 2020;Trebitsch et al. 2022).Using the UV mass-light ratio: log(M * ) = -0.49( +20.5) +8.8 at  = 8 in Stefanon et al. (2021), we estimate a stellar mass of log(M * ) ≃ 5.4 ⊙ at   = −13.5.

The results of the empirical model
We present the reionization history of the three scenarios (blue: extrapolation; green: restrained; red: flatten) in two models (  (UV) and   () in Figure 4.2, with the constraints derived from observations: Ly equivalent width (EW(Ly)) of galaxies (Mason et al. 2018(Mason et al. , 2019;;Hoag et al. 2019;Bruton et al. 2023 (Davies et al. 2018).As reference, we show two simple models where all galaxies have uniform constant   = 0.2 and  0  = 25.3 (  = 25.4),adopting either Schechter luminosity function or a double power-law luminosity function.
In both models, the restrained scenario (green dashed line) provides the closest match with observations, although the   () shows a later reionization history.For our models, we define the beginning redshift of the reionization as  90 when the IGM neutral fraction  H i = 0.9.For the   (UV) model,  90 ∼ 9.6, and for the   () model  90 ∼ 9.0.We define the ending redshift of the reionization as  10 with  H i = 0.1.For the   (UV) model,  10 ∼ 6.4, and for the   () model  10 ∼ 5.4.
In Figure 7 we focus on the restrained scenario and compute, for both models, the contribution to the reionization from galaxies with different UV luminosities.We group galaxies by their   into UV faint (−16 <   < −13.5),UV intermediate (−19 <   < −16), and UV bright (−23 <   < −19) galaxies.In both models, the UV faint galaxies are the major contributors of the reionization, producing ∼ 50% of the ionizing photons at all redshifts.The contribution of the UV intermediate galaxies increases as time goes by, becoming comparable to that of the UV faint galaxies at end of the reionization.UV bright galaxies, even though more likely to be LCEs at a given log( 32 ) and , contribute only 10% to 20% of the total ionizing photon budget.
The distribution of log( 32 ) and  can affect the reionization history and the relative contribution of galaxies with different UV luminosity.If the log( 32 ) has a higher mean, as reported in Sanders et al. (2023), more UV faint galaxies will be classified as LCEs in the restrained and the flatten scenarios, but not in the extrapolation scenario.The relative contribution of UV faint galaxies will be even higher in this case.If the scatter of the  is larger, individual UV bright galaxies with very blue  can bring up ⟨   ⟩ in all scenarios, while the UV intermediate and faint galaxies do not change as much since the  value is truncated at −3.5.In this case the relative contribution of UV bright galaxies in all scenarios will be higher.

CONCLUSIONS
In this paper we present a new analysis of the rest frame UV and optical spectra of a new a sample of  > 8 galaxies discovered behind the gravitational lensing cluster RX J2129.4+0009(Williams et al. 2022;Langeroodi et al. 2022).We combine these observations with those of the  > 7.5 galaxies for which similar data are available (Pontoppidan et al. 2022;Arellano-Córdova et al. 2022;Schaerer et al. 2022;Trump et al. 2022;Carnall et al. 2023;Curti et al. 2023;Rhoads et al. 2023;Mascia et al. 2023;Tang et al. 2023).
We compare the properties of these galaxies with those observed as part of the low redshift Lyman continuum survey (Flury et al. 2022a).The high [O iii]5007/[O ii]3727 emission line ratios (O 32 ) and steep UV continuum slopes,  < −2, of our sample are consistent with the values observed for low redshift Lyman continuum emitters, suggesting that these galaxies potentially contribute to the ionizing budget for the intergalactic medium.We use the H and UV luminosity to estimate the average ionizing photon production efficiency of our sample.
We apply a logistic regression (equation 3) to estimate the probability of a galaxy with   < −18 being a Lyman continuum emitter based on the measured   and  32 values, and explore three scenarios to account for the uncertain behavior at the faint end.Using this probability, we construct an empirical model that estimates the galaxy contribution to the reionization budget based on the observable quantities (  , ,  32 ).The preferred scenario in our analysis shows that at  = 8, the average escape fraction of the galaxy population (i.e., including both LyC emitters and non-emitters) varies with   .  is approximately 4% for bright galaxies (  < −19), and peaked at intermediate UV luminosity (−19 <   < −16).Galaxies with faint UV luminosity (−16 <   < −13.5) contribute half of the ionizing photons throughout the epoch of reionization.The relative contribution of faint versus bright galaxies depends on redshift, with the intermediate UV galaxies becoming more important over time.UV bright galaxies, even though more likely to be LCEs at a given log( 32 ) and , contribute the least of the total ionizing photon budget.

Figure 1 .
Figure 1.The rest-frame UV spectra of RXJ2129-ID 11027 at =9.51 (top) and RXJ2129-ID 11002 at =8.16 (bottom).The black lines are the spectra and the gray shaded areas are the uncertainty.In the RXJ2129-ID 11002 spectrum, the masked peak at 1125 Å is caused by a cosmic-ray hit.

Figure 2 .
Figure 2. The combined   indicator log(O 32 )− as a function of the absolute UV magnitude.The red, blue, and green hexagons are the RXJ2129 high- galaxies presented in this work.The yellow hexagons are the SMACS high- galaxies reported in the JWST Early Release Observation (Pontoppidan et al. 2022), and the yellow squares are the  > 7.5 galaxies in the GLASS-JWST program (Mascia et al. 2023).The low redshift LyC emitters (LCE, white circles) and non-LyC emitters (non-LCE, black dots) are separated in this Figure.The black dashed line is the logistic discriminator described as equation 3, and the gray dashed lines indicate where the P  =75% and P  =25%.It is unclear how this discriminator will extend at   faint than −18.5.We assume three simple extrapolations to account for this uncertainty: extrapolation (blue dashed line), restrained (green dashed line), and flatten (red dashed line).

Figure 3 .
Figure 3. Same as Figure 2 but  0  as a function of the absolute UV magnitude, where  0  =   (   =0).The gray shaded area marks the "canonical" log(  0  ) = 25.2 −− 25.3.The yellow squares are the  > 7.5 galaxies in the GLASS-JWST program(Mascia et al. 2023), and the yellow triangles are the  > 7.5 galaxies in the CEERS survey program(Tang et al. 2023)

Figure 4 .
Figure 4. Ionizing photon production efficiency measured over a wide range of redshift.The blue line and blue shaded area is the redshift evolution derived in Matthee et al. (2017), and the 1- limit.

Figure 5 .Figure 6 .
Figure 5.   and   () as a function of absolute UV magnitude for extrapolation scenario, the restrained scenario, and the flatten scenario.Top panel:   (red) and   () as a function of absolute UV magnitude.  is the probability of whether a galaxy is a LCE or not.  () is the average   value computed without accounting for   .Bottom panel: The blue line shows ⟨   ⟩ the average of   .The shadow areas are the the 90% and 10 % percentiles of the 500 runs.The observational results of Pahl et al. (2021) and Saldana-Lopez et al. (2023) are displayed with white circles and squares, respectively.The median   from SPHINX galaxies from Rosdahl et al. (2022) at  > 10 and  < 9 are indicated through green and orange lines, respectively.

Figure 7 .
Figure 7.The   of galaxies grouped by   as a function of redshift in the restrained scenario.The light to heavy green lines show the contribution from the UV faint (−16 <  < −13.5), intermediate (−19 <  < −16), and bright (−23 <  < −19) galaxies, respectively.The shaded areas show the 90% and 10 % percentiles of the 500 runs.The vertical dashed lines are the redshifts  90 and  10 , when the IGM neutral fractions  H i are 0.9 and 0.1, respectively.

Table 1 .
The derived properties for the RXJ2129 high- galaxies.()

Table 2 .
The derived properties for the SMACS0723 high- galaxies.