High Equivalent Width of H{\alpha}+[N II] Emission in z~8 Lyman-break Galaxies from IRAC 5.8{\mu}m Observations: Evidence for Efficient Lyman-continuum Photon production in the Epoch of Re-ionization

We measure, for the first time, the median equivalent width (EW) of H$\alpha$+[N II] in star-forming galaxies at $z\sim8$. Our estimate leverages the unique photometric depth of the Spitzer/IRAC $5.8\mu$m-band mosaics (probing $\approx 5500 - 7100$ A at $z\sim8$) of the GOODS Reionization Era Wide Area Treasury from Spitzer (GREATS) program. We median stacked the stamps of $102$ Lyman-break galaxies in the $3.6, 4.5, 5.8$ and $8.0\mu$m bands, after carefully removing potential contamination from neighbouring sources. We infer an extreme rest-frame EW$_0$(H$\alpha$+[N II])$=2328^{+1326}_{-1127}$ A from the measured red $[3.6]-[5.8]=0.82\pm0.27$ mag, consistent with young ($\lesssim10^7$ yr) average stellar population ages at $z\sim8$. This implies an ionizing photon production efficiency of $\log(\xi_{\mathrm{ion},0}/\mathrm{erg\ Hz}^{-1})=25.97^{+0.18}_{-0.28}$. Such a high value for photo production, similar to the highest values found at $z\lesssim4$, indicates that only modest escape fractions $f_\mathrm{esc}\lesssim0.3$ (at $2\sigma$) are sufficient for galaxies brighter than $M_\mathrm{UV}<-18$ mag to re-ionize the neutral Hydrogen at $z\sim8$. This requirement is relaxed even more to $f_\mathrm{esc}\le 0.1$ when considering galaxies brighter than $M_\mathrm{UV}\approx -13$ mag, consistent with recent luminosity functions and as typically assumed in studies addressing re-ionization. These exceptional results clearly indicate that galaxies can be the dominant source of reionizing photons, and provide us with an exciting glimpse into what we might soon learn about the early universe, and particularly about the Reionization Epoch, from forthcoming JWST/MIRI and NIRCam programs.


INTRODUCTION
The characterization of emission lines is a fundamental tool to study the physical processes governing the formation and evolution of galaxies. Hα constitutes one of the most reliable estimators of galaxy's star-formation rates (SFRs -e.g., Moustakas et al. 2006, Madau & Dickinson 2014) over short timescales ( 10 Myr -e.g., Kennicutt 1998, Kennicutt & Evans 2012 because it tightly corre- Email: stefanon@strw.leidenuniv.nl latates with the production of ionizing photons by OB stars, it does not depend on the metallicity, and it is less affected by dust attenuation compared to rest-UV lines. Moreover, because the rest-frame optical light correlates with the stellar mass (M ⋆ -e.g., Stefanon et al. 2017), the equivalent width (EW) of Hα provides a first estimate (modulo a M ⋆ /L optical factor) of the specific star-formation rate (sSFR -e.g., Fumagalli et al. 2012, Mármol-Queraltó et al. 2016, Faisst et al. 2016.
Optical and NIR spectroscopy have allowed astronomers to probe Hα up to z ∼ 2.5 (e.g., Fumagalli et al. 2012, Sobral et al. 2016, Nanayakkara et al. 2020. Progress at z ∼ 4 − 5 has been enabled by analyzing broad-band photometric data from Spitzer /IRAC and interpreting the observed blue [3.6] − [4.5] < 0 mag colors as the result of Hα emission contributing to the flux density in the 3.6µm band (e.g., Smit et al. 2016, Bouwens et al. 2016c, Rasappu et al. 2016, Faisst et al. 2016, Caputi et al. 2017, Lam et al. 2019, Faisst et al. 2019, Harikane et al. 2018a, Maseda et al. 2020. Constraining Hα at z 6 has proven to be quite challenging for a number or reasons. At these redshifts, the flux densities in both the IRAC 3.6 and 4.5µm bands are enhanced by nebular line emission ([O III]+Hβ, and Hα+[N II], respectively -e.g., Stefanon et al. 2021a), making it difficult to ascertain whether the observed colors are due to the combination of nebular line and continuum emission, or just to the continuum. This situation is exacerbated by the general lack of spectroscopic redshifts, essential for identifying which specific nebular lines could be contributing to the flux density in each band. Finally, further uncertainties are introduced by the still unconstrained line ratios at these early epochs (see e.g., Brinchmann et al. 2008, Steidel et al. 2014, Kewley et al. 2015, Faisst et al. 2016, Harikane et al. 2018b, Stefanon et al. 2022 for discussions on line ratios potentially evolving with cosmic time), and by the significantly shallower (3 − 6× depth) data currently available at 5 − 10µm (e.g., Stefanon et al. 2021b). These challenges have largely prevented us from securing an emission-line free continuum estimate at rest-frame optical wavelengths.
Fortunately, a favourable window exists again for 7.0 z 8.7. In this redshift range, Hα crosses into the IRAC 5.8µm band, [O III]+Hβ contribute exclusively to the 4.5µm-band flux density, and the 3.6µm band is free from significant line emission (e.g., Stefanon et al. 2022). Thus the potential exists to isolate the key lines to individual bands, and particularly separate Hα from significant contamination by other lines.
However, the use of the 5.8 and 8.0µm bands observations is not a panacea. Both bands suffer from lower sensitivities since 5.8 and 8.0µm band images could only be acquired during the Spitzer cryogenic mission, whereas those in the 3.6 and 4.5µm bands continued also during the warm mission. As a result, the sensitivities available in the 5.8 and 8.0µm bands are generally > 5−10× shallower than those available in the 3.6 and 4.5µm bands. For this reason, most of the studies characterizing the physical properties of EoR galaxies have so far focused on estimating the intensity of [O III]+Hβ from the measured [3.6] − [4.5] color (Smit et al. 2014, Castellano et al. 2017, De Barros et al. 2019, Stefanon et al. 2019, Bowler et al. 2020, Strait et al. 2020, 2021, Endsley et al. 2021b. Indeed, no detections of Hα from the color excess in the 5.8µm band have been published so far. The few estimates existing in the literature are indirect, and stem from converting the EW([O III]+Hβ) assuming standard Case B recombination coefficients and metal-line ratios either from best-fit SED analyses or extracted from tabulated values (e.g., Smit et al. 2014, Stefanon et al. 2022, Endsley et al. 2021b). Recently, attempts to detect emission in the 5.8µm band for individual z ∼ 7 − 8 Lyman-break galaxies (LBGs) were performed by Asada & Ohta (2022), leveraging the lensing magnification of the foreground galaxy clusters Abell2744, Abell1063, Abell370, and MACS-J0717 from the Hubble Frontier Fields program (HFF -Lotz et al. 2017) were not successful. The lack of direct detections has resulted in a chronic absence of direct constraints on the Hα intensity at these pivotal redshifts.
One possible approach to compensate for the current lack of deep 5 − 10µm data consists in combining the imaging available for samples of galaxies, and extracting their average properties. Our recently released IRAC mosaics from the GOODS Re-ionization Era wide-Area Treasury from Spitzer program (GREATS, PI: I. Labbé -Stefanon et al. 2021b) include all the relevant IRAC observations acquired in the four bands over the CANDELS (Grogin et al. 2011, Koekemoer et al. 2011 GOODS-N and GOODS-S fields (Giavalisco et al. 2004) across the almost 2 decades of Spitzer operations. Notably for this study, the GREATS 5.8µm-band mosaic provides ≈ 4−140× deeper coverage than what exists for the Abell2744, Abell1063, Abell370, and MACS-J0717 HFF fields 1 and 1.5× deeper coverage over a 4× larger area than the IRAC Dark Field 2 (Krick et al. 2009). The GREATS mosaics therefore constitute the deepest and thus most suitable data set for probing Hα emission in z ∼ 8 galaxies prior to JWST operations.
In this study, we explore for the first time the intensity of the Hα emission in z ∼ 8 galaxies by stacking the image stamps in the IRAC bands centered on 102 candidate LBGs at 7.3 < z phot < 8.7 identified by Bouwens et al. (2015b) in the CANDELS GOODS, UDS and COSMOS fields. This sample was already utilized by Stefanon et al. (2022) to study the rest-frame optical properties of z ∼ 8 galaxies, and benefits from minimal neighbour contamination (see Section 2 and Stefanon et al. 2022).
The layout of this paper follows: in Section 2 we present the data set and the sample adopted in this study; in Section 3 we describe the procedure we followed to estimate the average flux densities; the stacked photometry and the estimate of the rest-frame EW 0 (Hα) are presented in Section 4; we place our results in the context of the evolution of the EW 0 (Hα), sSFR and ξ ion in Section 5. A summary of this study together with our conclusions is presented in Section 6.

DATA AND SAMPLE
The sample adopted for this study consists of the 102 candidate z ∼ 8 Lyman-break galaxies previously discussed in Stefanon et al. (2022). Briefly, this sample is based on the Y −dropout LBGs that Bouwens et al. (2015b) identified over the CANDELS (Grogin et al. 2011, Koekemoer et al. 2011) GOODS-N, GOODS-S (Giavalisco et al. 2004, UDS (Lawrence et al. 2007) and COSMOS (Scoville et al. 2007) fields, the ERS field (Windhorst et al. 2011), and the UDF/XDF (Beckwith et al. 2006, Illingworth et al. 2013, Ellis et al. 2013) with the HUDF09-1 and HUFD09-2 parallels (Bouwens et al. 2011) 3 . In Table 1 we summarize the main properties of the adopted data sets. The mosaics are characterized by 5σ depths of ≈ 27.5 mag in the V 606 and I 814 bands, ≈ 26.7 − 27.5 mag in the Y 105 (GOODS fields) and 26.0 mag in the ground-based Y band (UDS and COSMOS), and ∼ 26.8 − 27.8 mag in the J 125 and H 160 bands.
A crucial aspect for this study is that these fields have excellent coverage in the Spitzer /IRAC 3.6, 4.5, 5.8 and 8.0µm bands. In particular, for the GOODS fields we adopted the mosaics and location-dependent pointspread functions (PSFs) from the GREATS program (PI: Labbé -Stefanon et al. 2021b). These mosaics combine all the useful IRAC data acquired across the full scientific life of Spitzer. As a result they are very deep, with 5σ depths of ∼ 26.0−27.0 mag in the IRAC 3.6 and 4.5µm bands, and ∼ 23.0 − 24.0 mag in the IRAC 5.8 and 8.0µm bands. While deep, we also require accurate PSFs to minimize the contamination from neighbours (see e.g., Stefanon et al. 2021b). Given the asymmetric nature of the instrumental IRAC PSF, particularly in the 3.6 and 4.5µm bands, and the variety of programs included in the mosaics, the PSFs can significantly vary across each field. The PSFs for GREATS are reconstructed by combining a high S/N empirical template rotated according to the position angle and weighted through the coverage depth from each program at the specific location.
We validated the redshift of the individual sources in our sample by running EAzY (Brammer et al. 2008) on the set of HST observations available for each source, while requiring 7.3 ≤ z phot ≤ 8.7. At these redshifts, the IRAC colors are sensitive to the intensity of the main rest-frame optical emission lines. A number of studies have shown that the IRAC colors can be successfully used to significantly reduce the photometric redshift uncertainties (e.g., Smit et al. 2014, Roberts-Borsani et al. 2016). However, this could also potentially bias our sample towards sources with strong emission lines. For this reason, we excluded the flux densities in the IRAC bands when running EAzY. Reassuringly, inclusion of IRAC fluxes in estimating the photometric redshifts has no strong impact on the sources we select (91 sources, corresponding to ∼ 89%, are in common between the two samples).
Because of the broad point-spread function of IRAC (≈ 1. ′′ 5 − 2. ′′ 0 from the 3.6µm to the 8.0µm band - Stefanon et al. 2021b), the extended light profiles of neighbouring objects could systematically affect the measurement of the emission of specific sources. For this reason, in our analysis we subtracted the neighbour emission with Mophongo (Labbé et al. 2006(Labbé et al. , 2010a(Labbé et al. ,b, 2013(Labbé et al. , 2015, and removed from the sample those sources where visual inspection still showed residual contamination (Stefanon et al. 2022).
To allow for a more meaningful comparison of our results with the literature, we also estimated the main stellar population parameters. These were computed by Figure 1. Image stamps (∼ 8. ′′ 0 per side) in the IRAC and HST bands, centered on the median stacks. The red circle marks the 2. ′′ 0 diameter aperture adopted for the photometry of the IRAC stacks. The HST stacks are presented to provide a better visual context of the data involved in our study, as the median flux densities in HST bands were estimated from the photometry of individual sources. Each stamp refers to a different band, as labeled at the top; in particular the HST optical stack combines all data available in the B435, V606, i775 and z850 bands. Remarkably, we find a 4.3σ detection for the flux density in the 5.8µm band, while a ∼ 1.8σ measurement in the similarly-deep 8.0µm stack. The striking visual difference between the detection in the 5.8µm band compared to those in the 3.6 and 4.5µm bands is a direct consequence of the ∼ 40× lower sensitivity available in the 5.8µm band. Our measurement suggests a significant contribution from Hα emission to the flux density in the 5.8µm band, as discussed in Section 4.   The black horizontal bars indicate the effective width of the bands. The blue curve corresponds to the best-fitting EAzY template. The inset presents the redshift probability distribution computed by EAzY. The labels at the top-left corner present the number of objects entering the stack, the median redshift, and the MUV computed by EAzY. Remarkably, the flux density in the 5.8µm band is comparable to that in the 4.5µm band, ∼ 2.2× higher than that in the HST and IRAC 3.6µm bands, indicative of strong nebular line emission from Hα. Also evident again is the lack of a prominent Balmer break between the H160 and 3.6µm bands, a result previously emphasised in Stefanon et al. (2022), the lead-up study to this analysis.

STACKING
Following Stefanon et al. (2022), we adopted distinct stacking procedures for the HST and for the IRAC bands, given the different photometric depths. For the HST bands, stacking consisted in evaluating the median of the extracted photometry normalized by the flux density in the H 160 band of each source, as the generally higher S/N characterizing these data reduces the measurement scatter around the true value. For the IRAC bands, however, the lower S/N compared to the HST data could introduce a larger scatter in the final measurement, possibly even systematically affecting it. We therefore constructed image stacks by taking the median of the image stamps centered on each source after they have been cleaned from neighbours using Mophongo and normalized by the H 160 flux density of each source. The stacked IRAC flux densities were measured in 2. ′′ 0diameter apertures. The smaller aperture adopted here compared to what Stefanon et al. (2022) used is a tradeoff between optimizing the S/N and minimizing potential flux loss introduced by the challenges in aligning the sources before taking the median and removal of neighbour contamination, particularly in the 5.8 and 8.0µm data. The aperture photometry was corrected to total using the median of the PSFs reconstructed at the location specific to each source. The applied correction factors are ∼ 2.2, 2.2, 2.9 and 3.3 for the 3.6, 4.5, 5.8 and 8.0µm bands, respectively. Uncertainties associated with the flux densities were computed by bootstrapping the sample 1000 times. Finally, all values were rescaled by the median of the flux densities in the H 160 band. An analysis adopting larger apertures (2. ′′ 5 and 3. ′′ 6) resulted in measurements consistent at 1σ with those obtained with the smaller aperture, albeit with larger uncertainties.
We further validated our 5.8µm-band measurement through a Monte Carlo simulation, presented in Appendix A. Briefly, we applied the same neighbour removal and stacking procedure we adopted for our main analysis to 102 synthetic sources. We added them to the IRAC 5.8µm mosaics, after normalizing their flux densities to those expected for the LBG in our sample, assuming a flat f ν SED and a rest-frame EW 0 (Hα)= 1900Å. This whole process was repeated 100 times. The resulting distribution of flux density measurements shows that, on average, we can recover the input flux density and that the impact of possible contamination by neighbours is negligible, as discussed and shown in more detail in the Appendix and in Figure 6. The depth of the IRAC 5.8 and 8.0µm mosaics in the COSMOS and UDS fields is ∼ 2 mag shallower (corresponding to ∼ 6× brighter flux limits) than the average depth in the GOODS fields, suggesting we should perhaps exclude or de-weight them in our stacks. On the other hand, these two CANDELS-Wide fields do allow us to incorporate 10 of the brightest z ∼ 8 sources with deep HST imaging from the Bouwens et al. (2015b) catalogs, providing a more comprehensive view of the prop-erties of z ∼ 8 galaxies. Even so, the median brightness of the sources selected in the COSMOS and UDS fields is only a factor ∼ 3 brighter than the median for the sources in the GOODS fields. To evaluate the impact of these sources on our stack results, we repeated our stacking analysis excluding the 10 sources in the COS-MOS and UDS fields. Reassuringly, the flux densities of the new measurements in the 3.6, 4.5 and 5.8µm bands differ by ∼ 5 − 10% from those obtained with full sample, after the change in H 160 normalization is taken into account. The flux density in the 8.0µm band is ∼ 60% fainter but still consistent at the 1σ level with that from the full sample. On balance, we therefore opted for including in our stack the z ∼ 8 galaxies from the UDS and COSMOS fields.

RESULTS
The stacked stamps in the IRAC bands are presented in Figure 1, while the photometry in those bands offering coverage for at least 90% of the sources in our sample is listed in Table 2, and displayed in Figure 2. Our photometric measurements are characterized by 20σ detections in the HST J 125 and H 160 bands, and ∼ 10σ in the 3.6 and 4.5µm bands. Remarkably, the stack in the IRAC 5.8µm band has resulted in a ∼ 4.3σ detection, while the stack in the 8.0µm band is characterized by a 1.8σ significance. Because the 5.8µm-and 8.0µm-band mosaics adopted in our study have similar depths and image quality, the 8.0µm band detection, even though of somewhat lower S/N, actually provides valuable added support for our 4.3σ measurement at 5.8µm as being a genuine detection. Together these detections give added confidence that the detection is real and is not significantly affected by neighbour and/or interloper contamination.
To assist in the interpretation of our measurements and to further validate the consistency of the stacked photometry, we also present the best-fit SED template from EAzY in Figure 2. For this step, we complemented the default set with SEDs of young (age∼ 10 6−8 yr) star-forming galaxies from BPASS v1.1 (Eldridge et al. 2017), whose nebular emission was computed with Cloudy (Ferland et al. 2013(Ferland et al. , 2017. . To our knowledge, this measurement constitutes the first detection of Hα from broad-band photometry in normal star-forming galaxies at z > 6.5. Recent attempts at measuring the intensity of Hα at z ∼ 8 for individual sources in HFF cluster fields were unsuccessful (e.g., Asada & Ohta 2022), likely due to the shallow coverage available in the IRAC 5.8µm band over those fields combined with the relatively low magnification values for the considered sources.
Such an elevated EW 0 (Hα) could originate from active galactic nuclei (AGN). Indeed, indication of nuclear activity in z ∼ 8 galaxies has been recently found (e.g., Laporte et al. 2017, Mainali et al. 2018, Topping et al. 2021. However, the extrapolation to z ∼ 8 of recent results at z 7 (e.g., Harikane et al. 2022) suggests that AGN would be a marginal population in the L < L * galaxies which dominate our sample. A more definitive assessment of the fraction of AGN in sub-L * galaxies at z ∼ 8 requires spectroscopic data, still unavailable for statistically significant samples.
The stacked SED also shows a red [3.6] − [4.5] = 0.54 ± 0.13 mag color. The increased flux density in the 4.5µm band is likely to result from substantial en- . This value implies an EW 0 (Hα)= 697 +160 −153Å . Our more direct measurement of EW 0 (Hα) based on the 5.8µm band excess differs by only ∼ 1.3σ from this estimate, providing further confirmation that Hα is very strong in these z ∼ 8 LBGs. Finally, our stacked SED is characterized by a J 125 − H 160 = −0.10 ± 0.07 mag color, indicating a blue UV slope (β ∼ −2.4), and a flat H 160 − [3.6] = −0.01 ± 0.11 mag suggesting young stellar population ages. Stefanon et al. (2022) already provide an extensive discussion of the interpretation of stack results involving these bands (J 125 , H 160 , 3.6µm, and 4.5µm).

Evolution of the EW 0 (Hα)
The large EW 0 (Hα) we infer requires very young stellar populations, few×10 7 yr (e.g., Inoue 2011, Wilkins et al. 2020, consistent with our age estimate based on multi-band photometry and with recent measurements at similar epochs (e.g., Stefanon et al. 2022, Endsley et al. 2021b, Strait et al. 2020, but see e.g., Roberts-Borsani et al. 2020, Tacchella et al. 2022. In Figure 3 we compare the EW 0 (Hα) estimate from this study to measurements at z < 7 from the literature. Specifically, we considered the estimates from  Nanayakkara et al. (2020). We included only measurements with M UV within ∼ ±0.75 mag of M UV = −19.9 mag, or whose stellar masses M ⋆ lie within ±1 dex of the stellar mass we estimated from our stacked photometry, when a typical M UV was not quoted with a result. Our estimate constitutes one of the highest EW 0 (Hα) measurements across the 0 < z < 8 redshift range. Nonetheless, the present large uncertainty of our measurement make it consistent at ∼ 1.5σ with the average EW 0 (Hα) existing at z ∼ 4 − 6.
Given that the samples in Figure 3 were comparably selected over the full redshift range, we can also address the question of the specific star-formation rate (sSFR) over this wide time baseline. At first approximation the EW(Hα) is proportional to the sSFR, and so we can also use the trends seen in Figure 3 to characterie the evolution with redshift of the sSFR. This is shown in the form sSFR∝ (1+z) 5/2 , as derived by Dekel et al. (2013), applying a constant conversion factor to the analytical expression of the evolution of the specific accretion rate of the dark matter halo under the hypothesys of a nonevolving ratio between the stellar mass and the mass of the host dark matter halo. We applied an overall normalization by fitting the curve to the observations. This simple relation can reproduce well the observations at z 4, but underestimates the expected EW(Hα) at lower redshifts, with larger gaps for lower redshift values.
While the uncertainty on the present EW 0 (Hα) measurement is large, one possible explanation for the differential evolution observed between the sSFR and the EW 0 (Hα) is a M ⋆ /L optical ratio (where L optical refers to the luminosity of the continuum at wavelengths close to that of Hα) evolving with cosmic time. This is indeed expected considering the increasingly larger fractions of evolved stellar populations at later cosmic times. The black solid curve in Figure 3 presents the result of applying a redshift-dependent M ⋆ /L optical ratio to the analytical expression of the sSFR(z) evolution. This factor was estimated from our default template set (Section 2), assuming galaxies started forming stars at z ∼ 20 (e.g., Mawatari et al. 2020, Harikane et al. 2022) 4 . We applied a global normalization factor from fitting the curve to the available observations. Here we do not consider the effects of dust attenuation given the growing indication that at the stellar masses considered here they are not a significant factor at z > 2 (e.g., Bouwens  Nanayakkara et al. (2020). The dashed grey curve marks the EW0(Hα) expected when its evolution follows that of the sSFR under a non-evolving star-formation efficiency scenario (e.g., Dekel et al. 2013) and constant M⋆/L optical ratio, while the solid black curve is for a M⋆/L optical ratio that increases with decreasing redshift according to a constant star-formation history. et al. 2016b, Dunlop et al. 2017, McLure et al. 2018, Bouwens et al. 2020) and marginal at z < 2 (e.g., Garn & Best 2010). The curve matches the observations reasonably well, in particular for z 6. Our measurement at z ∼ 8 is consistent at 1σ with the values expected from the new relation, although this is due, at least in part, to the large uncertainties. Nonetheless, the overall agreement in the recovered EW 0 (Hα) with the sSFR evolution to z ∼ 8 supports a scenario of a marginally evolving star-formation efficiency, as suggested by recent observational studies (e.g., Stefanon et al. 2017, Harikane et al. 2018a, Bouwens et al. 2021b, Stefanon et al. 2021a.

Constraints on ξion
We can now use our new determination of the EW 0 (Hα) to derive the efficiency of production of Hionizing photons (ξ ion ). This enables quantifying a key parameter, the total ionizing power of galaxies, in the heart of the reionization epoch and close to the time of instantaneous reionization (z = 8.8, Planck Collaboration et al. 2016). This is particularly valuable since only few, less direct, measurements exist at z > 7 (.e.g., Stark et al. 2015, De Barros et al. 2019, Endsley et al. 2021a), inferred from either spectral analysis of rest-UV emission lines or from SED fitting to broad-band photometry.
Following Bouwens et al. (2016a), we compute ξ ion from the production rate of Lyman-continuum photons, N (H 0 ). This can be inferred from the Hα luminosity L(Hα), using the relation of Leitherer & Heckman (1995): where L(Hα) has units of erg s −1 andṄ (H 0 ) of s −1 . This relation has a small ( 15% or 0.06 dex) dependence on the metallicity and electron temperature (e.g., Charlot & Longhetti 2001), which we assume as systematic uncertainty. This uncertainty is significantly smaller than the stochastic uncertainties associated with the EW 0 (Hα) measurement. The LyC photon production efficiency ξ ion,0 (where the subscript 0 indicates the assumption of an escape fraction f esc = 0, i.e., this is the actual production rate, in the galaxy, excluding any losses) can then be computed as:  Naidu et al. (2022) and Atek et al. (2022). We only considered measurements corresponding to MUV ∼ −19.9 mag or log(M⋆/M⊙) ∼ 8.1 ± 1.0 when the MUV information was missing. The red open circles flag those results at z < 4 whose sample was explicitly selected to have rest-optical lines with EW0 1000Å (Nakajima et al. 2016, Tang et al. 2019 or with likely a hard ionizing spectra . We arbitrarily shifted by ∆z = −0.06 our measurement for AV = 0.2 mag (open star) to improve the readability. The dashed line marks the result of a linear fit, with the 68% confidence interval encompassed by the grey shaded area. This composite set of measurements suggests a steady increase of ξion with increasing redshift.
where L UV is the UV-continuum luminosity computed from the stacked SED. The application of the above relations to our measurements yields ξ ion,0 = 10 25.97 +0.18 −0.28 Hz/erg, assuming negligible dust attenuation, as expected for L < L * LBGs at z > 4 (e.g., Dunlop et al. 2017, Bouwens et al. 2021b, Casey et al. 2021) and from the extrapolation of the results for LAE at lower redshifts (e.g., Naidu et al. 2022). If instead we consider a case with a small amount of dust attenuation, we obtain ξ ion,0 = 10 25.84 +0.18 −0.28 Hz/erg. In deriving this dust-impacted value we assumed, for simplicity, the Calzetti et al. (2000) attenuation law and the same A V = 0.2 mag value for both the stellar continuum and nebular emission, given the relative contribution of the two components is still quite uncertain (e.g., Buat et al. 2018, Shivaei et al. 2020, Li et al. 2021 and references therein). These high values of ξ ion require young stellar populations (ages 10 7 yr -e.g., Robertson 2021), consistent with the values we find from our SED fitting (see also Stefanon et al. 2022).
In Figure 4 we compare the value of ξ ion,0 from this study to previous estimates at similar redshifts and down to z ∼ 2 (Stark et al. 2015, Mármol-Queraltó et al. 2016, Nakajima et al. 2016, Bouwens et al. 2016c, Matthee et al. 2017, Harikane et al. 2018a, Shivaei et al. 2018, De Barros et al. 2019, Lam et al. 2019, Faisst et al. 2019, Tang et al. 2019, Nanayakkara et al. 2020, Emami et al. 2020, Endsley et al. 2021aand Atek et al. 2022. We only considered ξ ion estimates that refer to samples with either M UV within ±1 mag of the UV luminosity for our stack or log M ⋆ within ±1 dex of the stellar mass we estimated with FAST. Our measurement is consistent with the estimates existing at z ∼ 6.5 − 8 (Stark et al. 2015, De Barros et al. 2019, Endsley et al. 2021a. Their estimates derived, respectively, from modelling the intense C IVλ1548Å identified in the spectrum of a z = 7.045 galaxy (Stark et al. 2015), from the [C III] and Lyα lines in three z ∼ 7 galaxies with evidence for significant [O III] emission as suggested by their IRAC colors (Stark et al. 2017), and from SED fitting to multiwavelength photometry (De Barros et al. 2019, Endsley et al. 2021a). Our measurement is also broadly consistent with the estimates at z ∼ 5 (Bouwens et al. 2016c, Harikane et al. 2018a, Lam et al. 2019, Faisst et al. 2019 -see also Maseda et al. 2020 for exception-ally high ξ ion ≈ 10 26.3 Hz erg −1 in lower mass galaxies at z ∼ 4 − 5). Overall, these results suggest that ξ ion ≈ 10 25.6−25.8 Hz/erg could be typical at these epochs, and that ξ ion ≈ 10 25.7 might represent a reasonable estimate.
The values for ξ ion at z ∼ 2 − 3 are characterized by a large dispersion, ranging ≈ 10 24.7 − 10 26.0 Hz/erg. This distribution could be explained at least in part by selection effects. Remarkably, the values of ξ ion from samples characterized by high EW line emission (low-z analogues - Nakajima et al. 2016, Tang et al. 2019, see also Chevallard et al. 2018 for similar values at z ∼ 0), and from Lyα emitters ) are consistent with those found at z ∼ 7 − 8. Instead, the ξ ion estimated from more inclusive samples are generally lower (Matthee et al. 2017, Shivaei et al. 2018, Nanayakkara et al. 2020, Emami et al. 2020, Atek et al. 2022. A formal fit to the evolution of ξ ion , after excluding those from high EW samples, results in log(ξ ion /[Hz erg −1 ]) = (0.09 ± 0.01)z + (24.82 ± 0.08), whose slope is consistent with predictions from recent models (e.g., Finkelstein et al. 2019 -but see e.g., Matthee et al. 2022 for a nonevolving ξ ion model). The high value for ξ ion resulting from large EW samples at all redshifts z > 2 does indicate that we may be settling on a value broadly appropriate for early times for sub-L* LBGs.

Implications for fesc
The ionizing emissivity is generally expressed aṡ N (H 0 ) = ρ UV ξ ion f esc (e.g., Robertson et al. 2013), where ρ UV is the luminosity density at rest-frame UV (≈ 1500 − 1600 AA), and f esc is the fraction of ionizing photons escaping into the IGM. Because it is still uncertain whether the faint-end of the z ∼ 7 − 8 UV LF presents a turn-over (e.g., Atek et al. 2015, Livermore et al. 2017, Yue et al. 2018, Bhatawdekar et al. 2019, Bouwens et al. 2021b), here we explore the impact that different values of turnover have on the f esc . Our estimates are based on the requirement that all the necessary H-ionizing radiation for reionization is generated by stars. For this, we adopteḋ N (H 0 ) = 10 50.75 s −1 Mpc −3 (e.g., Bouwens et al. 2015a, Finkelstein et al. 2019, Mason et al. 2019, Naidu et al. 2020, and ξ ion,0 from this study. In the computation of ρ UV , we approximated the turnover by truncating the z ∼ 8 UV LF of Bouwens et al. (2021b) at values spanning −18 < M T < −12 mag.
The result of this procedure is presented in Figure  5. The larger value we find for ξ ion translates into f esc 30% for M UV > −18 mag. Values f esc 20% have been inferred at z > 6 by recent studies (e.g., Castellano et al. 2017). In particular, absence of a turnover in the faint-end slope down to M T ∼ −13 mag would only require f esc ∼ 5 − 10% to fully ionize the neutral H at z ∼ 8. These values are consistent with the f esc ∼ 5 − 10% inferred for sub-L * LBGs at z 4 by an increasing number of studies (e.g., Marchi et al. 2017, Naidu et al. 2018, Pahl et al. 2021, Siana et al. 2010, Grazian et al. 2016, Rutkowski et al. 2016, Steidel et al. 2018). Together they suggest only a marginal evolution of f esc with cosmic time for the average galaxy population. Furthermore, and qualitatively, such small escape fraction values can more easily be reconciled with the strong emission lines inferred at rest-frame optical for z 7 galaxies (EW 0 ([O III]+Hβ > 1000Å -e.g., Smit et al. 2014, Castellano et al. 2017, De Barros et al. 2019, Stefanon et al. 2019, Bowler et al. 2020, Strait et al. 2020, 2021, Endsley et al. 2021b. The overall consistency is reinforced by considering that photoionization modelling suggests that the production of such strong emission lines already requires very young (≈< 10 7 yr) stellar population ages (e.g., Inoue 2011, Wilkins et al. 2020).
An increasing number of studies are identifying Lyα emission in z 7 galaxies, with rest-frame EW ranging from ≈ 5 − 20Å to > 100 − 200Å (e.g., Pentericci et al. 2014, Stark et al. 2017, Hoag et al. 2019, Fuller et al. 2020, Endsley et al. 2021b, Pelliccia et al. 2021, Larson et al. 2022. A direct conversion of our EW 0 (Hα) using Case B recombination coefficients (L Lyα /L Hα = 8.7) suggests an intrinsic EW 0 (Lyα) intrinsic = 517 +287 −244Å . Under the assumption that the fraction of escaping LyC photons is approximately similar to that of Lyα emission (e.g., Steidel et al. 2018, Izotov et al. 2020, the implied f esc ranges between 10% and ≈ 40 − 50% ( Figure  5). These estimates are likely upper limits, given the still significant fractions of non-detection particularly for sub-L * galaxies (e.g., Pentericci et al. 2014, Jung et al. 2021. Thus the preliminary indications from Lyα studies are consistent with the more direct estimates, and reinforce the likely ready availability of adequate reionizing photons from star-forming galaxies for reionization. To give a sense of sensitivity to any dust, a thin A V = 0.2 mag Calzetti et al. (2000) dust screen would lower the Lyα flux by ∼ 1.7× (EW 0 (Lyα)= 298 +165 −141Å ), and increase the requirement on f esc by the same factor. Since very low dust absorption is likely, this suggests that any likely levels of dust would not change our conclusions significantly.

CONCLUSIONS
Our analysis of the deepest Spitzer /IRAC data available over extragalactic fields for a large sample of z ∼ 8 LBGs has allowed us to detect and measure for the first time the flux in the Hα line in the Figure 5. The black dashed curve and green shaded area mark the escape fraction and 68% confidence interval estimated as a function of the UV LF turn-over magnitude MT , as set for the needed ionizing emissivity, and based on the result for ξion, 0 from this study. The horizontal lines correspond to the fesc estimated matching the rest-frame EW0(Lyα) inferred from our Hα measurement to a compilation of values from the literature. Specifically, we considered the sample averages from Stark et al. (2017), Hoag et al. (2019) and Endsley et al. (2021a), the single z > 7 galaxy in Fuller et al. (2020) sample (C14215A1), and RELICS-DP7 from Pelliccia et al. 2021. early Universe, and to explore the resulting implications for re-ionization. Specifically, we obtained this measurement through a median stacking of Hubble and IRAC data for 102 LBGs initially identified by Bouwens et al. (2015b) from Hubble imaging over the CANDELS GOODS-N/S, ERS, XDF, CANDELS/UDS and CAN-DELS/COSMOS fields. Stefanon et al. (2022) had previously used a similar median stacking procedure to study the main properties of this sample of z ∼ 8 star-forming galaxies as a function of UV luminosity. These fields have deep coverage in the HST /ACS V 606 and I 814 and HST /WFC3 Y 105 , J 125 , JH 140 and H 160 bands. Key for our current study are that these fields also have deep Spitzer /IRAC mosaics from the GOODS Re-ionization Era Wide Area Treasury from Spitzer (GREATS -PI: Stefanon et al. 2021b). These mosaics combine all the relevant observations acquired with IRAC in the 3.6, 4.5, 5.8 and 8.0µm bands over the GOODS-N/S fields across the full scientific lifetime of Spitzer. In particular, the GREATS 5.8µm imaging is the deepest data available at ≈ 6µm before JWST, and represents a unique opportunity to probe Hα at 6.8 z 8.7.
We extracted median flux densities in the IRAC bands after combining the image stamps cleaned from neighbour contamination through Mophongo (Labbé et al. 2006(Labbé et al. , 2010a(Labbé et al. ,b, 2013(Labbé et al. , 2015, and recovered total flux den-sities using the location-specific PSFs from GREATS. Our main results are the following: • Our stack results for 102 galaxies at z ∼ 8 show a 4.3σ detection in the 5.8µm band, and a red [3.6] − [5.8] = 0.82 ± 0.27 mag color. • Interpreting the excess in the 5.8µm band as due to emission from Hα, we infer a rest-frame EW 0 (Hα)= 1960 +1089 −927Å , corresponding to a luminosity of log(L Hα /[erg s −1 ]) = 42.62 +0.15 −0.23 . Our result represents the first direct determination of the Hα intensity at z > 6.5.
These results allow us to draw the following conclusions: • Comparison of our new EW 0 (Hα) measurement with previous determinations at lower redshifts from the literature suggests that the trend of increasing EW 0 (Hα) with redshift (e.g., Faisst et al. 2016) can be extended up to z ∼ 8.
• After accounting for a M ⋆ /L optical ratio that depends on cosmic time, the observed evolution with redshift of EW 0 (Hα) is consistent with the evolution of the specific accretion rate of the dark matter halos, providing further evidence that the star-formation efficiency is at most marginally evolving with cosmic time in the early Universe.
• Following the formalism of Bouwens et al. (2016a), our new measurement of L Hα implies an efficiency of production of LyC photon ξ ion,0 = 10 25.97 +0.18 −0.28 Hz/erg. This constitutes one of the largest ξ ion estimates at 0 < z < 8 for sub-L * galaxies (M UV ∼ −19.8 mag, M ⋆ ≈ 10 8 M ⊙ ). While the uncertainties are large, our new measurement is very consistent with previous estimates at similar redshifts, at z ∼ 5, and with those values at lower redshift inferred from samples with significant nebular line emission. This consistency is not only reassuring but also points to a surprising uniformity across billions of years for star-forming galaxies.
• The large value of ξ ion we find suggests that escape fractions f esc 10% are sufficient for star-forming galaxies to fully ionize the neutral H at z ∼ 8 through escaping LyC radiation. The small value of f esc is consistent with what is seen at lower redshifts z ∼ 2−6 in star-forming galaxies, reinforcing the likelihood that galaxies alone are responsible for reionization.
It is remarkable to step back and realize that this study was enabled by observations in the 5.8µm band, acquired during the first few years of Spitzer scientific operations, a decade and a half ago. The present results highlight once again how powerful and pivotal a small telescope like Spitzer has been, especially when able to leverage robust selections made possible with HST. Fortunately, JWST /NIRSpec, NIRCam and MIRI combine and enhance the capabilities of HST and Spitzer, providing the potential for absolutely game-changing science in the coming years. We performed a Monte Carlo simulation to ascertain whether the signal detected in the 5.8µm band is genuine emission from the sample of LBGs at z ∼ 8 and not the result of residual contamination from neighbouring sources. We generated a new set of 102 flux densities in the H 160 band by randomly scattering the H 160 measurements of our sample according to their associated uncertainties. We then computed the flux densities in the 5.8µm band by assuming a constant ratio f 5.8 /f 160 = 2.2 between the flux density in the 5.8µm band (f 5.8 ) and that in the H 160 band (f 160 ), consistent with what we measure in our stack. This assumption is equivalent to a z ∼ 8 flat f ν SED, with a rest-frame EW 0 (Hα+[N II])=2300Å line emission contributing to the 5.8µm flux density. Point sources having the previously computed 5.8µm flux densities were then added at random locations across the four mosaics, adopting the location-specific PSFs from GREATS. In doing this, we preserved the relative fraction and luminosity distribution of sources in each field present in our original sample. Our adoption of point sources is supported by the smaller sizes (R e 1 kpc, corresponding to ≈ 0. ′′ 2 at z ∼ 8) of sub-L * galaxies at these redshifts (e.g., Shibuya et al. 2015, Bouwens et al. 2021a, compared to the 5.8µm PSF FWHM (≈ 2. ′′ 0). Following the same procedure we implemented for our main analysis (Section 3), we then constructed neighbour-cleaned stamps using Mophongo, and extracted the photometry of the median stack adopting a 2. ′′ 0-diameter aperture, correcting to total using the reconstructed PSF. All these steps were repeated 100 times. The results of this simulation are presented in the left panel of Figure 6. This clearly shows that our analysis is able to recover the median of the input flux densities. To test the amount Figure 6. Left panel: Comparison between the distribution of flux densities adopted as input for our Monte Carlo simulations (red histogram) and that after recovering the synthetic sources following the same procedure we adopted for our main analysis (yellow histogram -see also Section 3). The vertical black dashed line indicates the median of our measurements. Right panel: Distribution of flux densities measured by stacking neighbour-cleaned stamps centered at locations free from existing sources. These two results clearly indicate that our 5.8µm flux density measurement is genuine and any systematics resulting from the imperfect subtraction of neighbouring sources are negligible.

MS
of systematics that non-optimal removal of neighbouring sources could introduce into our measurements, we also extracted the photometry from neighbour-cleaned stacks centered on locations free from sources, as inferred from the combination of the J 125 , JH 140 -and H 160 -band mosaics (i.e., this is equivalent to adopting f 5.8 /f 160 = 0). The outcome of this second experiment is shown in the right panel of Figure 6. As we might expect if there is negligible contamination from the neighboring sources, the measurements are normally distributed around 0 nJy. The present Monte-Carlo simulation results significantly increase our confidence in the overall robustness of our 5.8µm-band flux measurements.