Towards unveiling the Cosmic Reionization: the ionizing photon production efficiency ( ξ ion ) of Low-mass H α emitters at z ∼ 2 . 3

We investigate the galaxy properties of ∼ 400 low-mass ( < 10 9 M ⊙ ) H α emitters (HAEs) at z ∼ 2.3 in the ZFOURGE survey. The selection of these HAEs is based on the excess in the observed K s broad-band flux compared to the stellar continuum estimated from the best-fit SED. These low-mass HAEs have elevated SFR(H α ) above the star formation main sequence (SFMS), making them potential analogs of the galaxies that reionized the universe during the epoch of reionization. The ionizing photon production efficiencies ( ξ ion ) of the low-mass HAEs have a median value of log( ξ ion /erg − 1 Hz ) = 25 . 24 +0 . 10 − 0 . 13 (25 . 35 +0 . 12 − 0 . 15 ), assuming the Calzetti (SMC) curve for the stellar continuum dust correction. This value is higher than that of main sequence galaxies by ∼ 0.2 dex at similar redshift, indicating that the low-mass HAEs are more efficient in producing ionizing photons. Our results also consolidate the trend of increasing ξ ion with redshift, but reveal a ”downsizing” relationship between ξ ion and stellar mass ( M ⊙ ) with increasing redshift. We further explore the dependence of ξ ion on other galaxy properties, such as the UV spectral slope ( β UV ), the UV magnitude ( M UV ), the equivalent widths ( EWs ) of H α and [O iii ] emission lines. Galaxies with the bluer UV slopes, fainter UV luminosities and higher equivalent widths exhibit elevated ξ ion by a factor of ∼ 2 compared to the median ξ ion of our sample. JWST data will provide an opportunity to extend our method and further investigate the properties of low-mass galaxies at high redshifts.


INTRODUCTION
The Epoch of Reionization (EoR) is a significant period in the cosmic history, when the ubiquitous neutral intergalactic medium (IGM) was ionized by the enormous amounts of ionizing photons (also called LyC or Lyman continuum) from the first luminous sources at that epoch.Current understanding suggests that reionization happens at z > 6 (e.g., Pentericci et al. 2014;Schenker et al. 2014).Several important aspects of the EoR still remain unclear, such as identifying the specific sources drive reionization (Madau & Haardt 2015;Robertson et al. 2015).Some studies propose that lowmass galaxies are expected to be major contributors to the cosmic reionization as the UV radiation budget is dominated by faint galaxies from the UV (∼ 1500 Å) luminosity function at z > 6 (e.g., Atek et al. 2015;Bouwens et al. 2015;Finkelstein et al. 2019).On the other hand, massive and bright galaxies are also proposed as major LyC leakers because of their larger escape fraction in model estimations (e.g., Gnedin et al. 2008;Naidu et al. 2020).
In order to fully understand how the universe was reionized, we must constrain the following, (i) the ionizing photon production efficiency, ξ ion , defined by the number of Lyman continuum (LyC) photons, i.e, ionizing photons, produced per UV luminosity (Robertson et al. 2013) and (ii) the escape fraction of ionizing photons, f esc , defined as a ratio of transmitted ionizing photons (λ < 912 Å) to input ionizing photons.Rest-frame optical emission lines serve as valuable tools for determining these parameters.The primary method for estimating ξ ion involves using nebular recombination lines, such as Hα, to infer the ionizing photon production rate.f esc , though derived from LyC flux emitted, exhibits correlations with various observational features of emission lines.However, accurately measuring these parameters at z > 6 presents considerable challenges.Directly measuring LyC flux is impossible due to the attenuation from the IGM along the line of sight.Meanwhile, galaxies at that epoch are too faint to construct a large and precise sample of emission line information.
Nevertheless, by conducting a detailed study of lowerredshift analogs, we can gain valuable insights into the physical properties of these early, faint galaxies at z > 6.Previous studies have already identified several lowredshift analogs.Individual galaxies with strong LyC leakage (Izotov et al. 2016a(Izotov et al. ,b, 2018a,b) ,b) are characterized by high emission line ratios of O32 = [OIII] / [OII] ≳ 5.This high O32 ratio has been suggested to be positively correlated with f esc (Nakajima & Ouchi 2014;Izotov et al. 2018a).Also, the inter-dependence of strong [OIII] emission line and f esc is pointed out by Nakajima & Ouchi (2014) through photoionization models.In another aspect, at lower redshifts, the estimation of ξ ion for individual galaxies at lower redshift has already been achieved using nebular recombination lines (e.g., Bouwens et al. 2016;Nakajima et al. 2016;Shivaei et al. 2018;Nanayakkara et al. 2020;Emami et al. 2020;Atek et al. 2022;Stefanon et al. 2022;Matthee et al. 2023).
Previous studies on ξ ion have used various samples: (i) Large sample of relatively massive galaxies (> 10 9.5 M ⊙ ) (e.g., Shivaei et al. 2018), (ii) Small sample of low-mass galaxies (e.g., Hayashi et al. 2016;Emami et al. 2020), (iii) Large sample of low-mass galaxies but at relatively low redshift (z ∼ 1; e.g., Atek et al. 2022).To gain a comprehensive understanding of ξ ion closer to EOR, pushing forward to large sample of low-mass galaxies at high redshift would be an important and necessary task.A large number of low-mass Hα emitters (HAEs) have been revealed at z ∼ 2.3 by Terao et al. (2022).These galaxies exhibit strong Hα emission lines primarily produced by massive stars (O and B stars) with lifetimes of less than 10 Myr.The possible presence of high-ionization states in these galaxies makes it reasonable to assume these low-mass HAEs as analogs of LyC leakages during EoR.Exploring the properties of them is therefore valuable in enhancing our understanding of cosmic reionization.In this work, we focus on a substantial sample of 401 low-mass Hα emitters (HAEs) at 2.05 < z < 2.5, which are identified using the photometric catalog from the FourStar galaxy evolution sur-vey (ZFOURGE).The sample reaches down to a stellar mass of 10 8 M ⊙ .For a comprehensive understanding of the selection criteria and the parent HAEs catalog, we refer readers to the detailed description provided in Chen et al. (2024) (hereafter C24).
The outline of this paper is as follows.We introduce the data in Section 2 and measurements of physical parameters in Section 3. We exhibit the relationship between ξ ion and various galaxy properties of the selected low-mass HAEs in Section 4. We further explore the observation and modeled results of ξ ion and discuss the implications of low-mass galaxies for their contribution to the cosmic reionization in Section 5. Finally we summarize our result in Section 6 and discuss future observations to further investigate the properties of HAEs at z ∼ 2.3.

THE GALAXY SAMPLE
In this study, we utilize the HAEs catalog from C24, which comprises 1318 HAEs at 2.05 < z < 2.5 from the ZFOURGE survey (Straatman et al. 2016).These HAEs are obtained from the flux excess in the observed ZFOURGE-K s broad-band to the stellar continuum estimated from the best-fit SED.The catalog includes information such as redshift, observed emission line fluxes, and SED-derived properties (e.g., stellar mass, dust attenuation) for the HAEs, with IDs corresponding to the ZFOURGE catalog.C24 refined the redshift information by incorporating spectroscopic redshifts (z spec ) from the MOSDEF survey (Kriek et al. 2015) and grism redshifts (z grism ) from the 3D-HST Emission-Line Catalogs (Momcheva et al. 2016).The median redshift of the HAEs is at z med = 2.25.
C24 conducted rigorous SED fitting to obtain accurate galaxy properties by using the CIGALE code (Boquien et al. 2019).The SED fitting, coupled with the reliable stellar continuum, allows for the derivation of emission line fluxes for Hα and [Oiii] based on the flux excesses (F excess ) in the deep ZFOURGE near-IR broadand medium-bandwidth filters.The derived emission line fluxes exhibit excellent agreement with spectroscopic and grism measurements from the MOSDEF and 3D-HST Emission-Line Catalogs.Rest-frame equivalent widths (EW s) of Hα and [Oiii] emission lines are calculated from the F excess and the continuum flux density from the best-fit SED model.In Table 1, we present the emission lines and the information of corresponding filters used in this analysis.

MEASUREMENTS
In this work, we quantify the star formation rates (SFRs) of galaxies by two indicators, Hα and FUV (1500 Å, L 1500 ) luminosities.For the UV luminosity measurement, galaxies are required to be detected in ZFOURGE B and V filter, which yields a cut of 2.5% of total HAEs.We compute the observed UV continuum (L 1500,uncor ) and the UV slope (β UV ) over a rest-frame wavelength range of 1400 − 2800 Å by performing a multi-band fitting to broad-band photometry with the relation of f λ ∝ λ βUV .For ZFOURGE-COSMOS field, the fitting includes broad-band photometry of B, G, V, R, R p , I, Z, Z p .For ZFOURGE-UDS and ZFOURGE-CDFS field, the fitting includes broad-band photometry of B, V, R, I, Z.The amount of dust attenuation is measured through the Bayesian SED fitting, providing the dust extinction E(B −V ) star , E(B − V ) neb .We parameterize the difference between E(B −V ) star and E(B −V ) neb by a factor f = 0.8 (Saito et al. 2020).Once we obtain the intrinsic UV luminosity (L 1500,cor ) and intrinsic Hα luminosity (L Hα,cor ), the star formation rates are converted by using the calibration of Kennicutt & Evans (2012) with a correction to the Chabrier (2003) IMF.
The ionizing photon production efficiency (ξ ion ) is defined as the ratio of the production rate of ionizing photons, N H 0 , in units of s −1 to the intrinsic UV luminosity (L UV,cor ) in units of erg s −1 Hz −1 : In this work, L UV,cor is derived at the rest-frame 1500 Å, and N H 0 is calculated as follows.For an ionizationbounded nebula, we could derive N H 0 from the intrinsic Hα luminosity, L Hα,cor , through the relation of Leitherer & Heckman (1995): × 10 12 L Hα,cor erg s −1 . (2) One important assumption in calculating ξ ion above is the zero LyC escape fraction, i.e, f esc = 0, since we assume the rate of recombinations balances the production rate of ionizing photons.The combination of Equation 1 and 2 actually gives us ξ ion,0 , while the intrinsic ξ ion should be ξ ion = ξ ion,0 /(1 − f esc ).
It is estimated that if the universe is totally reionized by galaxies, a large f esc ∼ 20% at z > 7 is needed (Ouchi et al. 2009).In contrast, after the universe is fully reionized, based on direct LyC imagings at lower redshift (2 < z < 4, e.g., Vanzella et al. 2010;Grazian et al. 2016;Matthee et al. 2017), an upper limit of the escape fraction is f esc < 6%.Therefore, ignoring the mentioned f esc correction would not significantly alter the trends we find in this analysis.We will assume ξ ion ≃ ξ ion,0 when investigating the relationship between ξ ion and galaxy properties in the following sections.
Typically, either the Calzetti et al. (2000) curve or the SMC curve is assumed for the reddening of the stellar continuum in high-redshift galaxies.Because the SMC curve has a steeper intrinsic UV slope, these two dust attenuation curves lead to different ξ ion , varying from less than 0.1 dex (e.g., Bouwens et al. 2016;Tang et al. 2019) to more than 0.2 dex (e.g., Shivaei et al. 2018;Atek et al. 2022).Thus, we also perform SED fitting using the SMC curve for the stellar continuum, while still adopting the Milky Way curve (Cardelli et al. 1989) for fitting the nebular emission in this study.

Star formation rates of low-mass Hα emitters
In Figure 1, we present the the star formation main sequence (SFMS) of Hα, UV for the HAEs from the C24 catalog.The low-mass HAEs are denoted as blue circles, while the median SFR in 6 mass bins is represented by large open circles.The mass bins are defined as follows: log(M * /M ⊙ ) < 8.5 for the first bin, log(M * /M ⊙ ) > 10.5 for the last bin, and the rest are divided into 0.5 dex widths.We apply the linear least squares regression to the SFR and stellar mass data points, above the mass completeness log(M * /M ⊙ ) = 9.0 (Straatman et al. 2016).The best-fit linear correlation between SFRs and stellar masses, with the 68% confidence interval on slope and intercept, is given by, log SFR(Hα) =( 0 It is obvious that those low-mass HAEs (< 10 9 M⊙) tend to scatter above the SFMS(Hα) but not obvious in SFMS(UV).Left: The star formation main sequence (SFMS) of 1318 HAEs at z med = 2.25 in the ZFOURGE fields based on Hα emission line.Blue circles are 401 low-mass HAEs with log(M * /M⊙) < 9, while grey circles show the other HAEs in our catalog.Open squares are median stacks in 6 mass bins, while the error bars on them represent the scatter in each mass bin.Blue solid line is the best linear fit to the galaxies with log(M * /M⊙) > 9.0, which is extrapolated to lower mass with blue dashed line.The best-fit SFMS from Whitaker et al. (2014), Speagle et al. (2014) and Shivaei et al. (2015) are also shown with black dotted, dashed, and dot-dashed lines, respectively.The error bars on the bottom-right corner represent the median uncertainty of low-mass HAEs (blue) and normal HAEs (black).Right: Same as the left panel, but SFR of each galaxy are calculated from the UV continuum.SFR(Hα) and SFR(UV) are corrected for dust attenuation using the Cardelli/Calzetti curve.We extrapolate the SFR − M * relation from Equation 3into low-mass domain (< 10 9 M ⊙ ) in Figure 1.In the right panel of Figure 1, the low-mass HAEs are found to be closer to the SFMS of UV, showing only a slight elevation in SFR(UV) compared to the SFMS by an average of 0.05 dex.However, a significant fraction of lowmass HAEs lie above the SFMS of Hα, exhibiting an average SFR(Hα) higher than the SFMS by 0.25 dex.To further illustrate this trend, we plot the SFR(Hα) vs. SFR(UV) ratio for our sample in 6 mass bins in Figure 2. In the high-mass domain (> 10 9 M ⊙ ), the ratio of these two SFR indicators is mainly close or below unity.In contrast, for galaxies with log(M * /M ⊙ ) < 9.0, the ratio is clearly above unity.We divide the HAEs in our study into two populations: 401 low-mass HAEs with log(M * /M ⊙ ) < 9.0 and 917 normal HAEs with log(M * /M ⊙ ) > 9.0.
The comparison of SFR(Hα) and SFR(UV) directly visualizes the burstiness of star formation activity.Previous studies (e.g., Weisz et al. 2012;Domínguez et al. 2015;Emami et al. 2019) have explored the time evolution of the Hα-to-UV ratios in variety types of star formation history (SFH) models.Generally, SFR indicators have different timescales, with SFR(UV) having a timescale of ∼100 Myr and SFR(Hα) having a timescale of ∼10 Myr.Emami et al. (2019) has found that in galaxies undergoing rising star formation, the luminosity ratio between Hα and UV is expected to be higher than unity, and SFR(Hα) will reside above the SFMS by up to an order of magnitude.High-resolution hydrodynamical simulations by Sparre et al. (2017) have also shown that the scatter on the SFMS is larger for the Hα-derived SFR because it is more sensitive to short bursts compared to the UV-based indicator.
We further explore the SFH of the HAEs in our sample in Figure 3.In this study, we apply the delayed-τ models to model the SFH of galaxies as follows: The distribution of the t/τ ratio of the HAEs reveals that the low-mass HAEs have a much smaller t/τ , indicating that they are in their early stage of star-forming with rising SFRs.On the other hand, the distribution of normal HAEs shows more variability, indicating a range of SFHs in these galaxies.The rising star formation in low-mass galaxies indicates the abundance of young stellar population (≤ 100 Myr) in these systems, which can contribute to intense Hα emission lines.In Figure 4, we plot the ratio of SFR(Hα) and SFR(UV) as a function of stellar mass, with an additional dimension of the equivalent width of Hα (EW Hα ).It is not surprising that galaxies with high EW Hα exceeding 1000 Å exhibit the highest SFR(Hα)/SFR(UV)ratios.In comparison, galaxies with EW Hα ∼ 100 Å would follow the SFMS.As shown in Figure 4, galaxy with the largest EW Hα (> 1000 Å) have the lowest stellar masses in our sample with a mean value of log(M * /M ⊙ ) = 8.09.Stefanon et al. ( 2022) conducted a stacking analysis of z ∼ 8 galaxies and found that the stacked galaxy exhibits extremely strong emission line, reaching up to EW Hα ∼ 2000 Å. Besides, the stellar mass estimated from their stacked photometry is log(M * /M ⊙ ) = 8.12, which is very close to our result.Therefore, it is reasonable to assume that the low-mass HAEs in our study are possible analogs of galaxies during EoR.

ξ ion evolution with galaxies properties
Here, we further investigate the ionizing photon production efficiency, ξ ion , of the low-mass HAEs and their relationship to various galaxy properties.Since f esc = 0 is assumed in our study, ξ ion is equivalent to the SFR(Hα)/SFR(UV) ratio by amplifying 1.31 × 10 25 .Thus,, ξ ion serves as a direct probe for the star formation activity in galaxies.ξ ion is also an important parameter in modeling the cosmic reionization and determining whether star-forming galaxies are capable of reionizing the universe (e.g., Robertson et al. 2013Robertson et al. , 2015;;Bouwens et al. 2016).In section 3, we have applied two different recipes, Cardelli/Calzetti and Cardelli/SMC, for the nebular/continuum dust correction in our SED fitting analysis.We present the results from both approaches in the following sections, with the individual data points in the figures corresponding to the Cardelli/Calzetti curve.

ξion and stellar mass
We display the relationship between ξ ion and the stellar mass in our sample in the left panel of Figure 5.
We also include open circles to represent the median ξ ion in different mass bins.To balance the number of samples in each bin, we combine the two most massive mass bins in Figure 1.The median ξ ion of the low-mass HAEs is log(ξ ion /erg −1 Hz) = 25.24+0.10 −0.13 (25.35 +0.12 −0.15 ), assuming the Cardelli/Calzetti (Cardelli/SMC) curve.This result is very close to the O3Es subsample in Tang et al. (2019) with EW [OIII] ≃ 300 − 600 Å (25.22).On the other hand, the median ξ ion of the normal HAEs is log(ξ ion /erg −1 Hz) = 25.05 +0.08 −0.10 (25.19 +0.10 −0.13 ), which is quite similar to that of the galaxies from the MOSDEF survey (25.06;Shivaei et al. 2018).When using an SMC extinction curve for the continuum, the ξ ion values are higher by ∼0.15 dex compared to the Calzetti curve.
The ξ ion distributions of low-mass HAEs and normal HAEs are exhibited in the right panel of Figure 5.The intrinsic scatter, represented by the standard deviation of the distribution, is 0.16 (0.19) dex for low-mass HAEs assuming the Cardelli/Calzetti (Cardelli/SMC) curve.For normal HAEs, the intrinsic scatter is 0.16 (0.26) dex.
We find that, on average, the low-mass HAEs have ∼0.2 dex higher ξ ion compared to the normal HAEs, indicating a higher efficiency of producing ionizing photons in the low-mass galaxies.In the high-mass domain, we observe no significant evolution of ξ ion with stellar mass.This trend is consistent with Shivaei et al. (2018) at z ∼ 2 using MOSDEF spectroscopic data.The ξ ion values from MOSDEF remain relatively constant down to 10 9.5 M ⊙ but show an increase in the lowest mass bin of 10 9 M ⊙ , suggesting a possible mass dependence of ξ ion .Our method successfully extends the mass range and provides evidence for the existence of a mass dependence of ξ ion below 10 9 M ⊙ .In contrast, at similar redshift, Emami et al. (2020) has found that ξ ion is generally independent of the stellar mass in their sample of 28 lensed dwarf galaxies, which span the range of 8.0 < log(M * /M ⊙ ) < 10.0.Their sample exhibits a nearly twice as large intrinsic scatter compared to our study.It is worth noting that high lensing magnification can also introduce significant magnification differences across the galaxy sample, leading to uncertainties when measuring ξ ion .Differences in sample selection and observational techniques may contribute to the discrepancy found in these studies.

ξion and UV Properties
It has been suggested that UV-faint galaxies could be significant contributors to cosmic reionization, based on several observations.One reason is the higher number density, indicated by the steep slope at the faint-end of the UV luminosity function (e.g., Atek et al. 2015;Bouwens et al. 2015).Another factor is the potentially higher LyC escape fraction for faint galaxies (Grazian et al. 2017).Investigating the UV properties of the lowmass HAEs can provide valuable insights into EoR.We report our results in Figure 6 and Table 2.
The left panel of Figure 6 shows the relationship between ξ ion and UV slope (β UV ) for our sample.We observe a gradual increase in ξ ion for galaxies with β UV < −2.0, with the bluest β UV objects showing an elevation of more than 0.2 dex.The UV slope is sensitive to the age of stellar population of massive stars (Oand B-type stars) in a galaxy.In ionization-bounded photoionization models (e.g., Topping et al. 2022), the bluest expected UV slope is β UV ≃ −2.6 in a dust-free case with stellar age of < 30 Myr.An older stellar population would lead to a redder UV slope up to β UV ∼ −2.0 and less production of ionizing photons, resulting in a lower ξ ion .Otherwise, ξ ion remains nearly unchanged for galaxies with redder β UV values (> −1.5).The effect of dust attenuation might become dominant in this region.Applying the SMC correction for the UV continuum leads to a similar result, albeit with a slightly lower elevation of ∼0.15 dex.
Our results are consistent with the literature results from Shivaei et al. (2018) 2, but represent different curve for the UV dust correction, while the error bars on them are the scatter in each bin.The error bars on the corner represent the median uncertainty of low-mass HAEs (blue) and normal HAEs (black), added with the number value of median UV slope and UV magnitude.It is clear that galaxies with bluer βUV and fainter UV luminosity are likely to hold larger ξion than others, encouraging that these galaxies would play an important role in reionizing the universe.
On the other hand, studies by Emami et al. (2020) and Onodera et al. (2020) have found no correlation between ξ ion and β UV .As mentioned above, Emami et al. (2020) used a sample of 28 lensed dwarf galaxies, with rest-frame equivalent widths up to 1500 Å.Also, the results from Onodera et al. (2020) are based on ∼ 20 extreme O3Es with rest-frame equivalent widths up to 2000 Å.We speculate that the discrepancy between our results  and previous studies may be attributed to the selection biases, given that their samples encompass a small number of extreme emitters, whereas our sample is more inclusive.
The low-mass HAEs in our study have a median value of β UV = −1.90,much higher than normal HAEs by ∼0.7.Comparing with the dust-free value of β 0 = −2.23 from the Meurer et al. (1999) calibration, we still have 72 low-mass HAEs (more than 1/6) with the bluest β UV < −2.23.Such blue β UV values suggest a very young stellar population with minimal dust content in the system.
Next, we explore the relationship between ξ ion and UV absolute magnitude (M UV ) in the right panels of Figure 6.In our sample, we observe an increase in ξ ion for the faintest galaxies compared to the brighter ones, with a difference of more than 0.3 dex.This trend suggests a dependence between ξ ion and M UV for the HAEs in our study.These results are similar to those observed in Lyman-alpha emitters (LAEs) at z ∼ 3 from Nakajima et al. (2018Nakajima et al. ( , 2020) ) with the faint end of UV magnitude to M UV ≃ −19.5 mag.On the other hand, Shivaei et al. (2018) and Emami et al. (2020) did not find a strong correlation between these two parameters in their studies.While, it should be noted that the galaxy sample in Shivaei et al. (2018) has been selected through spectroscopy, which may introduce a bias towards brighter objects M UV ≤ −21 mag.It can be also inferred from our results that brighter galaxies exhibit a weaker dependence between ξ ion and M UV .

ξion and nebular emission lines
The optical nebular emission lines in galaxies provide a wealth of information on the physical parameters, including the stellar population, chemical abundance, and ionization parameter, which may also be related to ξ ion .
As ξ ion serves as a useful indicator of the stellar populations in galaxies, with younger populations contributing more to the Hα emission lines, we expect to find a universal relationship between ξ ion and the equivalent width of Hα.This relationship is depicted in the left panel of Figure 7. Assuming the Cardelli/Calzetti curve, we fit a linear relationship between these two attributes: log ξ ion =(0.54 ± 0.03) × log(EW Hα ) + (23.76 ± 0.09). (5) For the Cardelli/SMC curve, we find a relationship with a slope of 0.51 ± 0.04 and an intercept of 23.96 ± 0.11.
Since large [Oiii] equivalent widths are typically produced by massive stellar populations with sub-solar metallicities, which meanwhile produce large amounts of ionizing photons, it is suggested that there also exists a correlation between ξ ion and EW [OIII] .This correlation has been observed in both local star-forming galaxies (Chevallard et al. 2018) and high-redshift emitters (Tang et al. 2019;Emami et al. 2020;Nakajima et al. 2020), indicating that systems with higher EW [OIII] are more efficient in producing ionizing photons.Our results also indicate a similar relationship between these two attributes, assuming the Cardelli/Calzetti curve: log ξ ion =(0.27 ± 0.04) × log(EW [OIII] ) + (24.43 ± 0.11). (6) Similarly, for the Cardelli/SMC curve, we find a relationship with a slope of 0.24 ± 0.04 and an intercept of 24.62 ± 0.10.For reference, we also overlay the best-fitting relations from Tang et al. (2019), Emami et al. (2020)  similar trend to literature, we find a clear discrepancy in terms of the slopes and intercepts.One possible reason for this discrepancy could be the treatment of dust extinction.Tang et al. (2019) has derived the relationship for their [Oiii] emitters (O3Es) using both the Calzetti law and the SMC law.Similar to our sample, applying the SMC law results in a similar slope but an elevation in the intercept by ∼ 0.15 dex.The difference in galaxy samples may also contribute to the discrepancy in slopes.Emami et al. (2020) has suggested that the slope in log ξ ion and log EW Hα (log EW [OIII] ) becomes shallower at lower equivalent widths.To explore this further, we derived the best-fitting result only for HAEs with EW Hα (EW [OIII] ) > 250 Å.Both curves shows a steeper slope of 0.61 ± 0.03 (0.32 ± 0.05), supporting the idea in Emami et al. (2020).While, even when considering larger equivalent widths and fitting the galaxies accordingly, the discrepancy with Tang et al. (2019) still remains.The difference in sample selection, with Tang et al. (2019) focusing on the most intense O3Es, may contribute to the discrepancy in the best-fitting relationship.
In any case, it is worth noting that the relationship between ξ ion and the equivalent widths of nebular emission lines, such as Hα and [Oiii], is observed globally at z ∼ 2. These strong nebular emission lines can serve as proxies for measuring ξ ion .

ξ ion evolution with redshift
Constraining ξ ion during EoR is an important task for cosmic reionization models.Measuring ξ ion for a large number of low-mass galaxies at z > 6 remains very challenging.Shivaei et al. (2018) has suggested a possible evolution of ξ ion with redshift, which could provide insights for extrapolating ξ ion to higher redshift.We compare our results with previous studies of ξ ion at various redshifts in Figure 8.Note that the galaxy samples included here are compiled in a heterogeneous manner, encompassing continuum-selected, stacked, and emission line-selected objects, which might introduce additional scatter.Atek et al. (2022) investigated ξ ion of a large sample of 1167 HAEs at 0.7 < z < 1.5 by using 3D-HST grism and imaging data.More than half of their 3D-HST sample are low-mass galaxies below 10 9 M ⊙ , which have a similar stellar mass range to our sample.Since the low-mass 3D-HST galaxies are also the highest-EW galaxies, we consider them as analog sample at lower redshift.Shivaei et al. (2018) used the spectroscopic and UV imaging data from the MOSDEF survey to measure ξ ion at z ∼ 1.5 and z ∼ 2.3 of 673 galaxies.The lowest stellar mass bin of Shivaei et al. (2018) is 10 9.3 M ⊙ , lacking galaxies with log(M * /M ⊙ ) < 9.0.The large number of low-mass HAEs found in our study are essential to complete the whole picture at z ∼ 2. Nakajima et al. (2016) used the spectroscopic measurement of Hβ emission line and SED-derived UV continuum to infer ξ ion of 13 LAEs as well as 2 Lyman Break   (Robertson et al. 2013).The blue solid line and the shaded area show the best-fitting linear regression results to all data with 95% confidence interval.We claim a evolution of ξion with lookback time, with the largest sample at z ∼ 2 so far.This evolution can be explained by more bursty star formation on average at higher redshifts.
Galaxies (LBGs) at z ≃ 3.1 − 3.7.Their LAE sample consists of less dusty system with color excess E(B − V ) almost consistent with zero.Also, the large [Oiii]/[Oii] ratio indicates the higher ionization properties of their LAE sample.Bouwens et al. (2016) measured the Hα emission line from the observed IRAC fluxes in ∼ 300 galaxies at z ≃ 3.8−5.0and 22 galaxies at z ≃ 5.1−5.7.They calculated Hα flux through the flux excesses in IRAC data, i.e, the [3.6]-[4.5]color.Besides, their result showed that assuming the SMC extinction law led to ∼ 0.1 dex higher ξ ion than that with the Calzetti one.Lam et al. (2019) performed a similar estimation of ξ ion from the flux excesses in IRAC as Bouwens et al. (2016) but by using the stacked IRAC colors in a sample of galaxies.Comparing to Bouwens et al. (2016), Lam et al. (2019) extended ξ ion measurements to even lower UV luminosities.The distribution of stellar mass and the large number sample in Lam et al. (2019) was com-parable to the low-mass HAEs in our study, indicating the similarity of the galaxy samples.Stark et al. (2015Stark et al. ( , 2017) ) discussed 4 individual, bright, lensed galaxies at z ∼ 7. ξ ion is inferred from stellar population and photo-ionization models.Significant detection of FUV emission lines of these galaxies, including Lyα, Similar to the sample from Nakajima et al. (2016) Stefanon et al. ( 2022) measured the median-stacked galaxy properties of 102 LBGs at z ∼ 8. Hα flux was inferred from the flux excess in stacked IRAC 5.8 µm band to IRAC 3.6 µm band, and ξ ion was further derived from the UV continuum computed from the stacked SED with the assumption of no dust correction.Their result constituted one of the largest ξ ion estimation.Castellano et al. (2023) calculated the ξ ion of more than 1000 VANDELS galaxies at 2.5 < z < 5 from a multi-band SED fitting approach with BEAGLE.They found no clear evolution of ξ ion with redshift within their probed range, but clear correlations with respect to stellar mass.Also Saldana-Lopez et al. (2023) directly measured ξ ion of the LAEs and non-LAEs from their rest-frame UV spectra.Their LAEs have ∼ 0.3 dex higher ξ ion respect to normal non-LAEs.Matthee et al. (2023) obtained emission line fluxes and physical properties for a sample of 117 O3Es at z ≃ 5.3 − 6.9, using the deep JWST/NIRCam wide field slitless spectroscopic observations.Measurements of physical properties in their study were based on the median stack spectra of O3Es and following SED fiting.ξ ion was obtained from Hβ emission line and SED derived UV luminosity.Also, a subset of 58 O3Es had spectral coverage of Hγ with E(B − V ) neb = 0.14 estimated from the observed Hγ/Hβ ratio.Whitler et al. (2024) selected a candidates of 27 galaxies from the large-scale galaxy overdensities surrounding UV luminous LAEs in the CEERS JWST/NIRCam imaging (Bagley et al. 2023).They modelled the SEDs of these galaxies and obtained the SED-inferred ξ ion by using the BEAGLE code.
Giménez-Arteaga et al. ( 2024) performed resolved SED modelling on a highly-lensed galaxy at z = 6.072 from JWST/NIRCam imaging.The SED results retrieve consistent measurement of emission lines to the IFU spectra.Based on the resolved SED fitting results, they provided 2D distribution of ξ ion of their target.Saxena et al. (2023) used JWST/NIRCam spectroscopic data from the JADES survey (Eisenstein et al. 2023) and directly measured ξ ion of 17 faint LAEs.Similarly, Boyett et al. (2024) used the same dataset and measured ξ ion of 28 [Oiii] or Hα EELGs.The stellar mass of each EELG is fitted by BEAGLE.
We perform a fit to all the listed data points, and the results indicate a clear evolution of ξ ion with redshift, which is consistent with previous studies in the literature (e.g., Shivaei et al. 2018;Atek et al. 2022).The best-fitting relation between ξ ion and redshift yields the following result, log ξ ion = (0.10 ± 0.02) × z + (24.92 ± 0.10). (7) The intrinsic differences among the sample selection in different studies should be taken into consideration.The current limitations in data availability make it challenging to identify a large number of individual galaxies at higher redshifts.Stacking methods, while useful, can introduce systematic differences among the samples.
Constructing large samples of individual low-mass galaxies is a challenging task, and currently, it has been undertaken in Atek et al. (2022) and our study.To obtain a clearer understanding of ξ ion at higher redshifts, it is crucial to construct large samples of individual lowmass galaxies.In section 4.2.1, we found an enhancement of ξ ion at low-masses.We compared our findings with several studies that also examined ξ ion and M * relationships below 10 9 M ⊙ , but at different redshifts.Atek et al. (2022) has analyzed 3D-HST galaxies at z ∼ 1 and also observed an enhancement of ξ ion at lower mass, around 0.5 dex, which is larger than our sample's ∼ 0.2 dex enhancement.Conversely, Lam et al. (2019) has stacked IRAC images of galaxies at z ≃ 4 − 5 and found an almost independent relationship between ξ ion and M * .By combining our result at z ∼ 2.3 in Figure 9, we speculate a possible "downsizing" relationship between ξ ion and M * over cosmic time, suggesting that the correlation between ξ ion and M * weakens from lower redshift to higher redshift.At the low-mass end, the difference in ξ ion among various redshifts is smaller than 0.5 dex.This discrepancy could potentially be mitigated or even canceled out if we extrapolate the M * − ξ ion relationship from our study and Atek et al. (2022) to lower masses, around log(M * /M ⊙ ) ≃ 7.5, as was in Lam et al. (2019).This trend suggests the possibility of an upper limit for ξ ion theoretically, and the low-mass galaxies are gradually approaching the upper limit ξ ion value.
To explore the theoretical limits of ξ ion , we investigate several stellar population synthesis (SPS) models.Notes.a ξ ion,model is derived from the ionizing photon production rate and the luminosity in the FUV band from the spectral synthesis outputs of the models.
6.0 6.5 7.0 7.5 8.0 8.5 9.0 9.5 10.0 log(age/yr)  Whitler et al. 2024).Also, the realistic SFHs of galaxies are more complex than a simple assumption of a constant SFH.If a galaxy's star formation rate is increasing, it is expected to have a higher ξ ion than that of a galaxy with a constant SFH at the same stellar age.Given that low-mass galaxies are more likely to have an increasing star formation rate, their upper limits of ξ ion may be even higher.
At the high-mass end, we find that the ξ ion of massive galaxies at z ∼ 1 is almost one magnitude lower than that at z ≃ 4 − 5. We propose that these galaxies may be at different stages of their star formation.Those galaxies at z ∼ 1 may have decreasing star formation rates, resulting in the lower observed ξ ion .In Figure 10, we also exhibit the evolution of ξ ion for a simple stellar population (SSP) from the BPASS model.If a galaxy undergoes quenching, ξ ion would rapidly decrease to log ξ ion ∼ 24.0 within 10 8 yr.On the other hand, massive galaxies at higher redshift in Figure 9 may still be undergoing continuous star formation.The discrepancy in their ξ ion is likely due to differences in their stellar populations.Referring to Figure 10, ξ ion would gradually decrease but eventually converge after 10 9 yr.For comparison, we include the observed results for massive galaxies at various redshifts in Figure 10.It is suggested that the BPASS models would be more representative for galaxies at z ≃ 4 − 5, while the BC03 models for our sample at z ∼ 2.3.

Implications for Reionization
Frameworks on cosmic reionization (e.g., Robertson et al. 2013Robertson et al. , 2015) ) highlighted the importance of constraining two key parameters: the escape fraction f esc and the ionizing photon production efficiency ξ ion , at z > 6. Bouwens et al. (2015) has established a well-matched evolution of the cosmic ionizing emissivity and the galaxy UV-luminosity density, suggest-ing that galaxies are the major sources responsible for reionization.They also constrained a lower limit of log(ξ ion f esc ) = 24.50 ± 0.1, if only galaxies contributed to the reionization.
In our study, we assume that the low-mass HAEs at z ∼ 2.3 can serve as the analog population of the galaxies that reionized the universe at z > 6.The median ξ ion of these low-mass HAEs is comparable to the canonical value of log ξ ion = 25.2 with an escape fraction f esc = 0.2 (Robertson et al. 2013;Bouwens et al. 2015).Note that our measurement of ξ ion assumes no escape fraction.Therefore, we can infer that the required escape fraction to ionize the universe at z > 6 is likely no larger than 0.2.
However, the assumption that low-mass galaxies at z ∼ 2 have similar physical properties to those at higher redshift may not be entirely accurate.For instance, low-mass galaxies at higher redshift are likely to exhibit higher ξ ion compared to our sample.Hydrodynamical simulations show that successive starburst activities are prevalent in low-mass galaxies (e.g., Domínguez et al. 2015;Sparre et al. 2017;Emami et al. 2019).While, variations in the underlying physical conditions of the interstellar medium (ISM) over cosmic time, such as changing metallicities, can lead to different levels of starbursts, resulting in varying equivalent widths of emission lines and different observed ξ ion at different redshifts.In section 5.1, we inferred a possible upper limit of ξ ion at log ξ ion ≳ 25.5.If we accept this value and assume log(ξ ion f esc ) = 24.60,f esc should be 0.13 at maximum for low-mass galaxies during reionization to reionize the universe.
Observations of LyC signal and the corresponding f esc at z ≃ 2−4 have been conducted by several studies in the past decade.Some surveys have reported very few individual detection of LyC emission, but inferred esacpe fraction f esc ∼ 0.05 − 0.1 based on stacking imaging (e.g., Marchi et al. 2017;Naidu et al. 2018;Steidel et al. 2018;Pahl et al. 2021).Still, there are also several individual galaxies that have been observed with very high escape fraction (e.g., Marques-Chaves et al. 2022).The above results, comparable to our estimation, indicate that galaxies could be the dominant source responsible for reionizing the universe.

CONCLUSIONS
In this work, we have investigated the galaxy properties, including SFR, ξ ion , UV properties, equivalent width, of the 401 low-mass (< 10 9 M ⊙ ) HAEs from C24 in three ZFOURGE fields.The selection of HAEs is based on the flux excess in ZFOURGE-K s filters due to the intense Hα emission lines relative to the best-fit stellar continuum from SED fitting.The main results can be summarized as follows: • The SFR − M * relation, i.e., SFMS, derived from Hα and UV luminosity exhibits a similar slope (0.56±0.03 vs. 0.60 ± 0.04) above the stellar mass completeness.The low-mass HAEs (< 10 9 M ⊙ ) scatter above the SFMS(Hα) by ∼0.3 dex, while this scattering is not evident in the SFMS(UV).
• Analysis of the SFR(Hα)/SFR(UV) ratio and star formation history suggests that these low-mass HAEs are undergoing early-stage bursty star formation phases with rising SFR.This is reflected in the abundance of younger stellar population with high EW Hα ∼1000 Å in the galaxies.Our population may be the analogs of the major contributors to cosmic reionization at z > 6.
• The ionizing photon production efficiency, ξ ion , of the low-mass HAEs is found to be log(ξ ion /erg −1 Hz) = 25.24+0.10 −0.13 (25.35 +0.12 −0.15 ), assuming the Calzetti (SMC) curve for UV dust correction.This result is higher by ∼0.2 dex compared to other HAEs in our study, suggesting a possible mass dependence of ξ ion .
• We observe a correlation between ξ ion and both the UV slope (β UV ) and UV absolute magnitude (M UV ).Galaxies with bluest UV slopes and faintest UV luminosities exhibit an enhanced ξ ion by nearly a factor of 2 compared to the median ξ ion of our sample.
• Our results confirm that ξ ion is related to the equivalent widths of Hα and [Oiii] (EW [OIII] and EW Hα ).This indicates that the equivalent width of these strong optical lines can serve as a proxy for estimating ξ ion .However, the relationship appears to be less steep compared to literature results, suggesting differences in the galaxy properties of the studied samples.
By combining a comprehensive analysis of literature results, we have strengthened the evidence for the evolution of ξ ion with redshift, extending our study to a significant number of low-mass galaxies at z ∼ 2.Moreover, our findings suggest a potential "downsizing" relationship between ξ ion and stellar mass as we trace back in cosmic time.We utilize stellar population synthesis (SPS) models to highlight that ξ ion observed in lowmass galaxies may approach the upper limit predicted by these models.This finding suggests that low-mass galaxies could be reaching the maximum efficiency of ionizing photon production according to current SPS models.This insight further supports the notion that low-mass galaxies play a crucial role in cosmic reionization.
To further enhance our understanding, the vast amount of data from the James Webb Space Telescope (JWST) will play a crucial role.This will enable the construction of a large sample of low-mass HAEs at higher redshifts, including even lower mass regimes, allowing for a more comprehensive exploration of cosmic reionization.

Figure 1 .
Figure1.Based on the flux excesses in Ks photometry, we define Hα emitters (HAEs) if they show > 3σ Hα detection in Ks photometry.It is obvious that those low-mass HAEs (< 10 9 M⊙) tend to scatter above the SFMS(Hα) but not obvious in SFMS(UV).Left: The star formation main sequence (SFMS) of 1318 HAEs at z med = 2.25 in the ZFOURGE fields based on Hα emission line.Blue circles are 401 low-mass HAEs with log(M * /M⊙) < 9, while grey circles show the other HAEs in our catalog.Open squares are median stacks in 6 mass bins, while the error bars on them represent the scatter in each mass bin.Blue solid line is the best linear fit to the galaxies with log(M * /M⊙) > 9.0, which is extrapolated to lower mass with blue dashed line.The best-fit SFMS fromWhitaker et al. (2014),Speagle et al. (2014) andShivaei et al. (2015) are also shown with black dotted, dashed, and dot-dashed lines, respectively.The error bars on the bottom-right corner represent the median uncertainty of low-mass HAEs (blue) and normal HAEs (black).Right: Same as the left panel, but SFR of each galaxy are calculated from the UV continuum.SFR(Hα) and SFR(UV) are corrected for dust attenuation using the Cardelli/Calzetti curve.

Figure 2 .
Figure2.Relation between the SFR(Hα) and SFR(UV) of the HAEs in our sample.SFR(Hα) and SFR(UV) are corrected for dust attenuation using the Cardelli/Calzetti curve.The color reflects different mass bins adopted in our study, as in Figure1.More clearly, we find the low-mass HAEs (< 10 9 M⊙) have the SFR(Hα)/SFR(UV) ratio close to a factor of 2, indicating their characteristic galaxy properties.

Figure 3 .Figure 4 .
Figure 3.The star formation history (SFH) of the HAEs in our sample.SFH of the low-mass HAEs (< 10 9 M⊙) and normal HAEs (> 10 9 M⊙) are distributed as a histogram of t/τ in cyan column and green column.For reference, a model delayed-τ SFH is shown as black solid line.The SFH suggests that rising star-forming activities are occurring in these low-mass HAEs.

Figure 5 .
Figure 5.The ionizing photon production efficiency (ξion) of all HAEs in our study.ξion is estimated from the Hα and UV luminosities followed by Equation 1 and 2. ξion of the low-mass HAEs is higher than other HAEs in our study by ∼0.2 dex, indicating a possible mass dependence.Left: Dependence of ξion on the stellar mass.The galaxy sample are separated into two populations, the normal HAEs and low-mass HAEs as in Figure 1.The error bars on the upper-right corner represent the median uncertainty of low-mass HAEs (blue) and normal HAEs (black), added with the median stellar mass of each subsample.Open squares and diamonds are median stacks in 5 mass bins, while the error bars on them represent the scatter in each mass bin.The square-shape stacks assume the Calzetti curve with for the UV dust correction, and the diamond-shape stacks assume an SMC curve.Right: The ξion distribution in the normal HAEs and low-mass HAEs assuming the Cardelli/Calzetti dust correction.The median values and errors for each population are log(ξion) = 25.24+0.10 −0.13 , 25.05 +0.08 −0.10 .The dashed lines in each figure indicate the canonical value of log(ξion)= 25.20 from Robertson et al. (2013).

Figure 6 .
Figure6.ξion as a function of the UV spectral slope (βUV, Left), the absolute UV magnitude at 1500 Å (MUV, Right).Symbols are as in Figure5.Open squares and diamonds are median stacks for these two galaxy properties in Table2, but represent different curve for the UV dust correction, while the error bars on them are the scatter in each bin.The error bars on the corner represent the median uncertainty of low-mass HAEs (blue) and normal HAEs (black), added with the number value of median UV slope and UV magnitude.It is clear that galaxies with bluer βUV and fainter UV luminosity are likely to hold larger ξion than others, encouraging that these galaxies would play an important role in reionizing the universe.
Measurments of the UV properties are only based on the observed photometry mentioned in section 3.3, which have no relationship to the SED fittings.a The uncertainty on each number is the median uncertainty of the subsample, which is different from the scatter shown in Figure10.b AGNs and QGs (the total number is 40) are excluded when deriving ξion.c The low-mass HAEs are those with stellar mass < 10 9 M⊙, while the normal HAEs are others.

Figure 8 .
Figure 8. Evolution of estimated ionizing photon production efficiency (ξion) with redshift.Our measurement at z ∼ 2.25 (blue and red stars) are compared with literature results at a wide range of redshifts, include Stark et al. (2015, 2017); Bouwens et al. (2016); Nakajima et al. (2016); Shivaei et al. (2018); Atek et al. (2022); Stefanon et al. (2022); Castellano et al. (2023); Saldana-Lopez et al. (2023); Matthee et al. (2023); Saxena et al. (2023); Whitler et al. (2024); Giménez-Arteaga et al. (2024); Boyett et al. (2024).The spectroscopic-selected samples are marked as pentagonal symbols, the photometric-selected sample as circular ones and the individual sample as square ones.If multiple values exist in the literature, we take the result based on the Calzetti curve for UV continuum.The numbers next to the markers are the average stellar mass from each sample, and those without stellar mass are marked as open symbols.The HAEs from our sample are marked with different colors.The dashed horizontal line indicates the canonical value of log(ξion)= 25.20(Robertson et al. 2013).The blue solid line and the shaded area show the best-fitting linear regression results to all data with 95% confidence interval.We claim a evolution of ξion with lookback time, with the largest sample at z ∼ 2 so far.This evolution can be explained by more bursty star formation on average at higher redshifts.

Figure 9 .
Figure 9.The relation between ξion and stellar mass at different redshifts.Red squares are stacks in 5 mass bins, same as the left panel in Figure 9.A possible downsizing of the relation between ξion and M⊙ with the increase of redshift is shown here.

Table 1 .
The ZFOURGE data used to derive emission line fluxes Notes.aThe Hα emission line at 2.05 < z < 2.5 redshifts into the Ks filter.bThe [Oiii] emission lines redshift into Hs filters at z < 2.26, and into H l filters at z > 2.26.

Table 2 .
Median ξion for HAEs at z ∼ 2.3 separated into different bins of stellar masses, UV luminosities and UVslopes.
Atek et al. (2022))22)in Figure7.Although our sample shows a Figure 7. ξion as a function of the equivalent width of Hα (EWHα, Left), equivalent width of [OIII] (EW [OIII] , Right).[OIII] are derived from the flux excesses in Hs/H l photometry.Symbols are as in Figure 5. Open squares and diamonds are median stacks for various galaxy properties, but represent different curve for the UV dust correction.The error bars on the corner represent the median uncertainty of low-mass HAEs (blue) and normal HAEs (black), added with the number value of median EWHα and EW [OIII] .For both panels, the best-fit relation between ξion and EWHα (EW [OIII]) is shown with a black dashed line where we use the Calzetti curve for dust correction.Other relationship derived for HAEs at 0.7 < z < 1.5(Atek et al. 2022), [OIII] emitters at 1.3 < z < 2.4(Tang et al. 2019), lensed low-mass galaxies at z ∼ 2(Emami et al. 2020) are shown with purple, blue, brown dashed-dotted line, respectively.Note thatEmami et al. (2020)andAtek et al. (2022)used the SMC curve for dust correction, which results in a higher estimation of ξion.

Table 3 .
Various SPS models and corresponding ξion

and
Figure10, we present the results from these models under the assumption of a constant SFH.If we assume a stellar age of 30 Myr in the constant SFH scenario, the estimation of ξ ion from the BPASS model with an upper mass cutoff of 100 M ⊙ and sub-solar metallicity is consistent with the observed log ξ ion ∼ 25.5 at z ≃ 4 − 5.It should be noted that the actual stellar ages may vary among the samples, and even younger stellar populations with an age of 10 Myr would exhibit ∼ 0.15 dex higher ξ ion compared to those at 30 Myr.This estimation is in closer agreement with recent observations from JWST at z > 8 (