Quantifying the Escape of Lyα at z ≈ 5–6: A Census of Lyα Escape Fraction with Hα-emitting Galaxies Spectroscopically Confirmed by JWST and VLT/MUSE

The James Webb Space Telescope provides an unprecedented opportunity for unbiased surveys of Hα-emitting galaxies at z > 4 with the NIRCam's wide-field slitless spectroscopy (WFSS). In this work, we present a census of Lyα escape fraction (f esc,Lyα ) of 165 star-forming galaxies at z = 4.9–6.3, utilizing their Hα emission directly measured from FRESCO NIRCam/WFSS data. We search for Lyα emission of each Hα-emitting galaxy in the Very Large Telescope/MUSE data. The overall f esc,Lyα measured by stacking is 0.090 ± 0.006. We find that f esc,Lyα displays a strong dependence on the observed UV slope (β obs) and E(B − V), such that the bluest galaxies (β obs ∼ −2.5) have the largest escape fractions (f esc,Lyα ≈ 0.6), indicative of the crucial role of dust and gas in modulating the escape of Lyα photons. f esc,Lyα is less well related to other parameters, including the UV luminosity and stellar mass, and the variation in f esc,Lyα with them can be explained by their underlying coupling with E(B − V) or β obs. Our results suggest a tentative decline in f esc,Lyα at z ≳ 5, implying increasing intergalactic medium attenuation toward higher redshift. Furthermore, the dependence of f esc,Lyα on β obs is proportional to that of the ionizing photon escape fraction (f esc,LyC), indicating that the escape of Lyα and ionizing photon may be regulated by similar physical processes. With f esc,Lyα as a proxy to f esc,LyC, we infer that UV-faint (M UV > −16) galaxies contribute >70% of the total ionizing emissivity at z = 5–6. If these relations hold during the epoch of reionization, UV-faint galaxies can contribute the majority of UV photon budget to reionize the Universe.


Introduction
The epoch of reionization (EoR) represents a phase transition when the intergalactic medium (IGM) evolved rapidly from mostly neutral to highly ionized (Fan et al. 2006;Robertson 2022).Although it is generally acknowledged that EoR is completed by z ∼ 5-6 (e.g., Kulkarni et al. 2019;Zhu et al. 2021;Jin et al. 2023), its onset and major contributors remain elusive.Recent measurements on quasar luminosity functions (LFs) at z ∼ 6 suggest that the quasar population only contributes a negligible fraction (7%-10%) of total photons to keep the Universe ionized (Matsuoka et al. 2018;Jiang et al. 2022).Young star-forming galaxies are thought to be the main producers of ionizing photons for reionization.However, one of the biggest uncertainties is whether a vast number of faint galaxies (e.g., Finkelstein et al. 2019) or a small population of bright oligarchs (e.g., Naidu et al. 2020) dominate the cosmic ionizing emissivity during the EoR.To address these questions, extensive studies in the past decades focused on three key quantities observationally: the ionizing photon production efficiency ξ ion , the UV luminosity density ρ UV , and the escape fraction of Lyman-continuum (LyC) photons f esc,LyC (Robertson et al. 2013;Bouwens et al. 2016).
Among the above three elements, the biggest uncertainty is from the estimates of f esc,LyC at high redshifts (Vanzella et al. 2010;Bian & Fan 2020;Naidu et al. 2022).Directly observing the leakage of LyC photons at z  4 (Steidel et al. 2018) is difficult due to the high opacity of the IGM.Therefore, many works studied low-redshift galaxies with similar properties (e.g., Leitherer et al. 2016;Flury et al. 2022a;Saldana-Lopez et al. 2022) or indirect tracers of f esc,LyC (e.g., Flury et al. 2022b;Xu et al. 2022;Choustikov et al. 2023), although it is still unclear whether low-redshift analogs can be good representatives of the main contributors to reionization (e.g., Katz et al. 2020;Schaerer et al. 2022), and whether those indirect indicators of f esc,LyC derived from very small sample could be generalized (e.g., Katz et al. 2022;Saxena et al. 2022).The typically required f esc,LyC for reionization to end at z ∼ 5-6 is approximately 5% (Robertson et al. 2015;Finkelstein et al. 2019), and recent work have constrained the average f esc,LyC to be 10% through various indirect methods (Steidel et al. 2018;Pahl et al. 2021;Saldana-Lopez et al. 2023).
To directly constrain the f esc,LyC for galaxies at z  4, one of the most promising ways is to study the escape of Lyα.Numerous studies have shown a strong connection between the escape of ionizing photons and Lyα photon (Chisholm et al. 2018;Marchi et al. 2018;Maji et al. 2022).Their escapes are both determined by the distribution and geometry of neutral gas in the interstellar medium (ISM), and facilitated by optically thin channels within the gas reservoir (Dijkstra et al. 2016;Gazagnes et al. 2020;Lin et al. 2023).These ionized channels could be carved by feedback from massive young stars, which also modulate the production of LyC and Lyα photons (Vanzella et al. 2022;Kim et al. 2023).
A number of previous efforts have been made to probe the escape of Lyα ( f esc,Lyα ), which could be derived by comparing Lyα and Hα emissions.Before the launch of James Webb Space Telescope (JWST), Hα of galaxies at z  5 could be estimated only through the Spitzer broadband color (e.g., Faisst et al. 2019;Lam et al. 2019;Stefanon et al. 2022).Hence, the measurements of f esc,Lyα in the prelaunch era suffer from large modeling and photometric uncertainties.Improved estimates of the Hα flux excess have been done recently using JWST narrowband observations (Ning et al. 2023;Simmonds et al. 2023).Moreover, the spectroscopic capability of JWST in the infrared regime enables us to detect and accurately measure Hα emission lines in faint, high-redshift galaxies.The JWST Near Infrared Spectrograph (NIRSpec) can capture both the Lyα and Hα emission of high-redshift galaxies simultaneously (Chen et al. 2024;Saxena et al. 2024) and thus quantify f esc,Lyα , but a census of f esc,Lyα by NIRSpec is difficult due to the limited sample size it can observe and the trade-off between high sensitivity and high spectral resolution.The wide-field slitless spectroscopy (WFSS) of JWST/NIRCam (Greene et al. 2016;Rieke et al. 2023), instead, offers a unique capability to map the spectra of all sources within its field of view (FoV; Sun et al. 2022Sun et al. , 2023)), ideal to build a large sample of emission line galaxies through blind spectroscopic surveys.Combining the Hα emission lines captured by WFSS with the corresponding Lyα lines observed by ground-based facilities, we are able to, for the first time, conduct a direct census of f esc,Lyα with an unbiased sample limited by Hα flux.
In this work, we utilize JWST/NIRCam WFSS around 4 μm to detect Hα emission over the redshift range of z ≈ 4.9-6.3,when the Universe was experiencing rapid transition, and galaxies were undergoing substantial mass, dust, and chemical content assembly (Madau & Dickinson 2014;Davidzon et al. 2017;Faisst et al. 2019).The ground-based Very Large Telescope (VLT)/Multi Unit Spectroscopic Explorer (MUSE; Bacon et al. 2010), working in the full optical domain, perfectly covers their Lyα so that we can obtain more precise measurements of f esc,Lyα .
The paper is organized as follows.In Section 2, we outline the photometric and spectroscopic data used and sample selection based on Hα emission lines.We describe the measurements of galaxy properties, search for Lyα lines, and methods to calculate f esc,Lyα in Section 3. In Section 4, we present the distribution of f esc,Lyα , its dependence on galaxy properties and redshift evolution.We then investigate the connection between LyC and Lyα escape in Section 5, and discuss the implication for reionization.The final summary is presented in Section 6.We show more details about the analysis and address potential systematics in Appendix A. Throughout this work, a flat ΛCDM cosmology is assumed, with H 0 = 70 km s −1 Mpc −1 , Ω Λ,0 = 0.7, and Ω m,0 = 0.3.

Data and Sample
We describe the data set that we use in this work and the construction of our sample in this section.We introduce all the photometric and spectroscopic data, including JWST, Hubble Space Telescope (HST), and VLT/MUSE, in Section 2.1.The selection of Hα emitters (HAEs) is presented in Section 2.2.

Photometric and Spectroscopic Data
2.1.1.JWST/NIRCAM Imaging and Slitless Grism in the GOODS-S Field The JWST "First Reionization Epoch Spectroscopically Complete Observations" (FRESCO) Survey (Oesch et al. 2023, GO-1895;PI Oesch), covers 62 arcmin 2 in each of the extragalactic legacy GOODS/CANDELS fields with deep NIRCam imaging and WFSS observations.In GOODS-South field, FRESCO obtains 4 × 2 pointings with 8 × 7043 s exposures using the F444W grism, accompanied by 8 × 3522 s exposures of F210M imaging and 11 × 4456 s of F182M imaging.Additional 8 × 934 s exposures of F444W imaging are taken to associate the slitless spectra with individual objects.
The imaging data were first reduced up to stage2 using the standard JWST pipeline v1.7.29 and the calibration reference files "jwst_1014.pmap."The calibrated single exposures ( * _cal.fits) were further processed by GRIZLI.10GRIZLI mitigated 1/f noise, alleviated "snowball" artifacts from cosmic rays, and converted the world coordinate system (WCS) information in the headers of each exposure to a simple imaging polynomial format so that exposures can be drizzled and combined with ASTRODRIZZLE.The WCS of each exposure image was registered using GAIA DR3 (Gaia Collaboration et al. 2023), and the images were finally drizzled with a pixel scale of 0 031 using pixfrac=0.8.An additional background subtraction using PHOTUTILS.BACK- GROUND2D was performed on the final mosaics after masking bright sources.
For the WFSS data, after the standard stage1 steps, flatfielding was performed using the flat images of direct imaging since the WFSS flat reference files are not available yet in CRDS.1/f noise is removed only along columns as the spectra are dispersed along rows.The WCS for each individual exposure was assigned using CALWEBB.assign_wcs.There were astrometric offsets between this assigned WCS and that of the fully reduced F444W images.To associate the WFSS spectra and individual galaxies accurately, we calculated the offsets using accompanied short-wavelength exposures (i.e., the F182M and F210M imaging taken simultaneously with the WFSS) and corrected them for each WFSS exposure.These preprocessed files, as well as a source catalog based on the final F444W mosaic, were provided as inputs for GRIZLI for spectral extraction using the updated sensitivity functions and tracing models as described in Sun et al. (2023).11To facilitate emission line searching, we performed running median filtering along the dispersion direction on single calibrated grism exposures as described in Kashino et al. (2023), and extracted continuumremoved 1D spectra aiming for emitter searching only.

HST Imaging in GOODS-S Field
The high-level products of the Hubble Legacy Fields12 (HLF; Whitaker et al. 2019) in GOODS-S include images from ultraviolet to infrared bands, observed over the past 18 yr, with 7491 individual exposures and 6.4 million s in total.All the HLF image mosaics are tied to an absolute GAIA DR2 (Gaia Collaboration et al. 2018) reference frame.In this work, we use four optical images in the Advanced Camera for Surveys (ACS)/WFC F435W, F606W, F775W, F814W bands at 30 mas pixel −1 , and four infrared images in the WFC3/IR F105W, F125W, F140W, and F160W bands at 60 mas pixel −1 .

Source Extraction and Photometry
We first resampled all HST images onto a fixed grid with a pixel scale of 0 031.For HST mosaics, we constructed empirical point-spread functions (PSFs) by stacking stars with good signal-to-noise ratios (S/Ns) in the fields.For JWST mosaics, the centers of most bright stars are saturated especially in the F182M and F210M images, so it is hard to build empirical PSFs that can accurately reproduce the complicated PSF shapes with these limited numbers of stars in the final mosaics.Therefore, we adopted the modeled PSFs provided by WEBBPSF. 13We generate PSF matching kernels to match the PSFs of HST WFC3/IR F160W or JWST/NIRCam F444W, using PHOTUTILS, 14 with high-frequency noise filters in Aniano et al. (2011).All HST ACS/WFC mosaics (F435W, F606W, F775W, F814W), JWST/NIRCam F182M, F210M mosaics were convolved to match the PSF of JWST/NIRCam F444W with corresponding kernels; all HST WFC3/IR mosaics (F105W, F125W, F140W) were matched to the PSF of F160W.
We ran SOURCE EXTRACTOR (Bertin & Arnouts 1996) in dual mode for each band, using the F444W mosaics as the detection image.We started with Kron aperture photometry with Kron factors k = 1.2.We measured the flux uncertainty for each object by randomly positioning the same Kron aperture in empty regions nearby.The measured flux and uncertainties were further corrected to total values by multiplying the ratio of flux inside the k = 2.5 versus k = 1.2 Kron apertures on the detection image.Finally, we correct for flux outside the k = 2.5 radii using the PSF of F444W.Note that for HST WFCS/IR bands we homogenized the images to the PSF of F160W, slightly larger than that of JWST/NIRCam F444W; thus, we enlarged the aperture size on these HST IR mosaics by a factor of 1.2, which ensures that all the flux from objects would be included as tested on the F160W PSF.

VLT/MUSE in GOODS-S Field
MUSE is an integral-field spectrograph in the optical wavelength range on the VLT.MUSE has an FoV of 1 arcmin 2 , a spatial sampling of 0 2, and a spectral range from 4750 to 9350 Å with a mean resolution of 3000.Its high sensitivity to emission lines in the full optical domain makes it an ideal instrument to detect high-redshift Lyα emitters (LAEs) over z = 2.9-6.7.
The GOODS-S field is covered by the legacy MUSE-Wide survey (Urrutia et al. 2019) and MUSE Hubble Ultra Deep Field (UDF) survey (Bacon et al. 2023).MUSE-Wide comprises 60 × 1 arcmin 2 pointings, each with a depth of 1 hr.The MUSE UDF survey includes 3 × 3 arcmin 2 mosaic of nine 10 hr depth fields (MOSAIC), a 1 × 1 arcmin 2 31 hr depth field (UDF-10), and the circular MUSE eXtremly Deep Field (MXDF).The MXDF covers a radius of 41″ and 31″ for respectively 10+ and 100+ hr of depth, reaching a final maximum depth of 141 hr, and is the deepest optical spectroscopic survey so far.In Figure 1, we show the footprints of FRESCO, MUSE-Wide, and MUSE UDF.The overlap region between FRESCO and the two MUSE surveys is about 50 arcmin 2 , including ∼43 arcmin 2 of the MUSE-Wide, ∼9 arcmin 2 of the 10 hr MOSAIC, and the whole UDF-10 and MXDF.
We adopt the publicly released DR115 data for MUSE-Wide and AMUSED16 data for MUSE UDF.For MUSE-Wide, DR1  Helton et al. (2023).The crosses denote Hα emitters not covered by MUSE (41 in total); dots are within MUSE coverage (197 in total), and dots with black edges are selected sources we use for analysis (157 in total).Right: the redshift distribution of our Hα emitters.The distribution of selected Hα sources is wrapped with black edges.
only releases 44 of the 60 data cubes, so we reduced the remaining 16 data cubes with our custom pipeline.This includes the standard pipeline17 v2.8.4 (Weilbacher et al. 2020) and improved the sky subtraction using ZAP18 (Soto et al. 2016).We aligned all the data cubes to the Hubble ACS astrometry as done by the MUSE-Wide team.The offset between the HST astrometry and the GAIA DR2, adopted in WFSS and imaging, is ΔR. A. = + 0 094 ± 0 042, and Δdecl.= − 0 26 ± 0 10 (Bacon et al. 2023).
The second criterion requires a complete 912 Å dropout for z > 4.9 galaxies at F435W.The third criterion describes the Lyα break features; we require at least 1 mag drop in F606W that is blueward of the wavelength of Lyα.The fourth criterion is the F444W excess due to Hα emission lines when considering possible Balmer jumps.This is the minimum excess for z = 4.9-6.7 mock galaxies in JAGUAR (Williams et al. 2018) with Hα flux larger than 1.6 × 10 −18 erg s −1 cm −2 (4σ depth of the FRESCO F444W grism data according to Oesch et al. 2023).We also obtained photometric redshifts for each source using EAZY (Brammer et al. 2008) with the corr_sfhz_13 templates, which contain redshift-dependent star formation histories.We excluded objects with robust photo-z, i.e., P(z p > 4.5) = 0, as low-redshift objects.We finally constructed a z > 4.9 emitter catalog by careful visual inspection, during which we removed artifacts on imaging and grism data, false detection by the search algorithms, and misidentified lines from the contamination of other sources.

HAE Sample
The high-redshift emitter catalog described above comprises all z > 4.9 candidates, including z > 7 [O III] emitters.Bright [O III] emitters are easily identified by their doublets at restframe 4959 and 5007 Å.However, as the [O III] 4959 Å lines are ∼3 times fainter than the 5007 Å, and Hβ can be even fainter, in many cases, we can only detect a single 5007 Å emission while the 4959 Å lines remain undetected (see Appendix A for more details).Since our science goal is more sensitive to the sample purity than completeness, we exclude all possible [O III] candidates based on both photometry and spectroscopy.We first assumed the brightest single emission line in the 1D spectra to be 5007 Å lines and force-fitted 4959 lines.This would effectively identify [O III] emitters with marginally detected 4959 lines.We then further extend the [O III] sample by photometric cuts as follows: 1. S/N <2 in both F606W and F775W, 2. F606W-F115W>1.7 and F814W-F115W>1.7, 3. F115W-F210M<1.0 and F115W-F210M<1.0, 4. F814W-F115W>(F115W-F210M)+1.5.Endsley et al. (2022) to select z > 6.5 emitters with high [O III] equivalent width (EW).Note that we do not require F814W to be completely undetected, allowing for Lyα transmission spikes within its wide wavelength coverage (Kakiichi et al. 2018).

Similar criteria have been applied in
We find 74 possible [O III] candidates.29 of the 74 [O III] candidates reside in the JADES-DEEP footprint (Eisenstein et al. 2023;Rieke et al. 2023).Among the 29 sources, 21 sources have photo-z > 7 based on multiple bands of JADES imaging, with [O III] emission in F444W; eight of them have photo-z of 5 − 6 with Hα in F444W.We argue that our criteria for [O III] could, in most cases, successfully remove z > 7 [O III] emitters and thus keep an Hα sample with high purity.
Finally, we obtain a sample of 222 HAEs.We identify blended clumpy galaxies with irregular morphology and multiple cores as single objects.These clumpy HAEs have Hα emission lines that cannot be distinguished on the 2D grism spectra, and some of them have unsolved Lyα emission captured by MUSE, even though SEXTRACTOR may label them as different individuals.During the extraction of grism spectra, we use a modified segmentation map as the input of GRIZLI, in which clumpy segments are combined.In our final sample, any two "single" objects are required to be separated from each other by at least 3.5 kpc (i.e., 0 6 at z = 5.5).We measure the total Hα flux on 1D extracted spectra by fitting a "Gaussian + constant" model, where the constant is used to depict the local continuum level around Hα.
Moreover, we label sources having neighbors within 25 kpc (about 4″ at z = 5-6) and ±500 km s −1 .Due to the significant disparity between the spatial resolutions of MUSE and JWST, we will exclude these galaxies and their close neighbors in the following analysis, since (1) Lyα around individual galaxies can be extended to 3″ (Wisotzki et al. 2016) and will boost Lyα in the overlapping regions of close pairs, introducing systematics during the stack of Lyα as described in Section 3.3.2,and (2) for very close pairs, it is difficult to obtain accurate Hα flux due to the overlap issue along the dispersion direction of grism.Among all the 222 HAEs, 21 of them have close neighbors.Three of the 222 HAEs are bright point sources with hints as being active galactic nuclei (AGNs).There are 184 HAEs that reside in the MUSE footprints, with 17 having close neighbors and two as potential AGNs.Finally, 165 HAEs are selected for the following studies.The spatial and redshift distributions of our selected sample are shown in Figure 1.Detailed properties of these HAEs are presented in Appendix A, which demonstrate that they are typical main-sequence galaxies in this redshift range.

Measurements
We describe the methods we use to quantify the properties of our HAE sample in Section 3.1.The search and measurements of Lyα emission using VLT/MUSE data are outlined in Section 3.2.The calculation of f esc,Lyα for individual galaxies and through a stack analysis is presented in Section 3.3.

SED Modeling and Galaxy Properties
We infer the physical properties of the 165 z ≈ 5-6 HAEs by spectral energy distribution (SED) modeling with the Bayesian code BEAGLE (Chevallard & Charlot 2016), using the measured photometry (see Section 2.1.3)and Hα line flux as inputs.We do not use the Lyα flux as inputs because modeling Lyα and Hα simultaneously requires correctly taking the Lyα escape process into account, and current SED modeling codes cannot do it.We assume a constant star formation history and the Chabrier (2003) initial mass function (IMF) with an upper limit of 100M e .We adopt flat priors for metallicity, 0.01 < Z/Z e < 0.2, and ionization parameter for nebular emission (Gutkin et al. 2016) in log-space, U 3 log 1 -< <-, which is the dimensionless ratio of the number density of H-ionizing photons to that of hydrogen.For dust attenuation, we adopt two attenuation curves: the Small Magellanic Cloud (SMC) dust law (Pei 1992, hereafter SMC) and the Calzetti et al. (2000) law (hereafter Calz) with R V (ratio of total to selective extinction in V band) being 2.74 and 4.05 respectively.We set the optical depth in the V band varying in the range of 0-0.5 in log-space.Although a large variety of dust attenuation curves for high-redshift galaxies have been used in the literature (Reddy et al. 2015;Scoville et al. 2015;Reddy et al. 2018), these two are commonly suggested to be applicable and utilized broadly in high-redshift studies (e.g., Bouwens et al. 2016;Shivaei et al. 2018;Faisst et al. 2019).In the following analysis, we adopt the Calz as our fiducial dust attenuation law, but we also consider the SMC law for comparison.
The BEAGLE SED fitting results can provide the stellar mass (M å ), the effective optical depth at the V band (τ V ), the attenuated UV luminosity L UV,att , and the intrinsic UV luminosity L UV,int directly.The UV luminosity is defined at rest-frame 1600 Å, and the optical depth could be converted to the color excess E(B − V ) with corresponding attenuation curve.We measure the observed UV slope (β obs , not corrected for dust) by fitting the stellar continua of the best-fit SED models with a power-law function ( f obs l = l b ), spanning from 1260 to 2500 Å.The SED-derived UV luminosity and UV slopes are statistically consistent with the fitting results using the photometry directly.We estimate the specific star formation rate (sSFR) for each galaxy with their intrinsic Hα luminosity and SED-derived M å .The star formation rate (SFR) of each galaxy is converted based on Kennicutt (1998) specialized for the Chabrier (2003) IMF: Note that the conversion factors between the Hα luminosity and SFR for high-redshift galaxies vary among different stellar population models and IMFs, and are highly dependent on the ionizing photon production efficiency ξ ion (e.g., Theios et al. 2019).The conversion factor of subsolar-metallicity galaxies is about 2.5 lower than the canonical conversion factors of lowredshift galaxies (e.g., Reddy et al. 2018;Theios et al. 2019;Shapley et al. 2023).Nevertheless, as we focus on the relative relation between f esc,Lyα and SFR in the following analysis (see Section 4.2), Equation (1) is qualitatively a good proxy to the SFR of our Hα-emitting galaxies.
To derive the sSFR surface density (Σ sSFR ) , we fit the F444W surface brightness distribution, i.e., the approximate Hα surface brightness, for each galaxy using single or multiple Sérsic models.The effective radius (r eff ), which is the radius of the circular aperture including half of the total flux from the galaxy, is measured in the corresponding compound modeled image after PSF deconvolution. 19We do not adopt r eff in simple single Sérsic modeling due to the clumpy morphology as well as the highly coupled parameters during modeling.Then, we calculate Σ sSFR as We further measure the ionizing efficiency ξ ion , i.e., the production rate of LyC ionizing photons per unit nonionizing UV continuum luminosity.ξ ion can be expressed in terms of N H 0 ( ), the production rate of LyC photons, and L UV,int as f esc,LyC is the escape fraction of LyC photons from galaxies.Our result, as shown in Section 5, shows that f esc,LyC has negligible impact to the estimate of ξ ion , so we assume f esc,LyC as zero in Equation (4).Therefore, ξ ion we adopt is actually the production efficiency of ionizing photons failing to escape from the galaxy.We present the EW distribution of our HAE samples in Figure 2 (left panel).The EWs are determined by The distribution of EW Hα and EW Lyα of the 165 selected HAEs in the overlapping region of FRESCO and MUSE footprints.The EWs are calculated following Equation (5) with f cont from the SED-modeled stellar continuum assuming the Calz dust attenuation law.We label the mean uncertainty of EW Hα (pink in the left panel) and EW Lyα (blue for LAEs, and green for non-LAEs in the right panel) as horizontal lines, respectively.For non-LAEs, we note that the uncertainty of forced measured EW Lyα reaches an average of 163 Å when EW Lyα > 100 Å, although the mean uncertainty for non-LAEs is about 47 Å.We mark the EW Lyα = 0 as the gray dashed line.
where λ 0 , λ 1 define the wavelength range for integration, F line is the emission line flux, and f cont is the rest-frame continuum flux density.For Hα, we obtain f cont by fitting the 6530-6600 Å stellar continuum provided by the best-fit SED using powerlaw models, and extrapolate to the wavelength of Hα.The Hα EW (EW Hα ) of our sample ranges from 93 to 4300 Å, with a median value of 665 ± 67 Å.

Search for Lyα Emission
We search for Lyα emission for all the HAEs covered by MUSE observations.Instead of the blind line algorithm usually adopted in MUSE 3D data cubes, we now know exactly the position and redshift of our targets and can well predict where the Lyα emission lines could emerge.For each Hα source, we extract 20″ × 20″ subcubes and perform median filtering, with a window size of 101 pixels to remove continuum and nearby sources.For sources in overlap regions of multiple pointings, we combine all the individual subcubes weighted by their exposure time. 20We then search for local surface brightness peaks on 3 pixel width pseudo-narrowband slices and 4 × 4 pixel boxes, centered on the source position and its rest-frame 1216 Å.The box size is determined by the typical size of UV star-forming region of high-redshift galaxies shown in Kusakabe et al. (2022).We perform the search within a radius of 1 5 spatially and a spectral range of 1500 km s −1 .This radius ensures the inclusion of potential Lyα flux given a possible offset between the Lyα and Hα due to gas distribution and radiation transfer (e.g., Cai et al. 2019;Zhang et al. 2023).For each selected flux peak, we extract a 1D spectrum using the 4 × 4 pixel box, and fit a skewed Gaussian profile.We visually inspect the extracted spectra and identify 63 HAEs with Lyα emission lines detected at >2σ.In the following analysis, we refer to the HAEs with Lyα lines detected as "LAEs" and the remaining as "non-LAEs"; but we note that this definition fully depends on the depth of the MUSE observations and is not comparable to the literature definitions (e.g., Lyα EW >10 Å).Among the 63 LAEs, 39 have been reported in Kerutt et al. (2022) or Bacon et al. (2023).
We further extract 1D spectra for all the selected 165 HAEs using r = 1 5 circular apertures, and measure the uncertainties by the standard deviation of 500 spectra extracted using equally sized apertures randomly positioned in each subcube.For LAEs, we fix the skewed Gaussian profile derived from the boxy-aperture extracted spectra to maximize S/Ns, and adjust the amplitude to fit the circular-aperture extracted spectra and obtain the total flux.Lyα of two LAEs are overwhelmed by noises in their circular-aperture extracted spectra; we thus adopt the boxy-aperture measured flux for them.Comparing the flux of LAEs we measure with the values reported in Kerutt et al. (2022) and Bacon et al. (2023), the median flux ratio is 1.04 ± 0.38 implying a good consistency.For non-LAE, we measure the integrated flux and uncertainties over ±1000 km s −1 around the wavelength of Lyα within the circular apertures.All non-LAEs do not show emission line peaks around the Lyα wavelength and have integrated Lyα fluxes with S/Ns <3.
Figure 3 shows six examples of LAEs and non-LAEs among our HAE sample.We show the distribution of Lyα EW (EW Lyα ) of LAEs and non-LAEs in the right panel of Figure 2. EW Lyα are measured following Equation (5), with f cont obtained by fitting the 1250-1350 Å stellar continuum.Our detected LAEs have rest-frame EW (EW Lyα ) spanning from 19 Å to 462 Å, with a median of 76 Å.For non-LAE, EW Lyα derived from the forced measured Lyα within ± 1000 km −1 spans from -175 Å to 441 Å, with a median value of 11 Å.Their 3σ upper limits range from 8 Å to 1123 Å with a median value of 91 Å.

Measurement of Lyα Escape Fraction
In this section, we describe our method to measure f esc,Lyα in detail.We outline the measurement for individual HAEs in Section 3.3.1 and the stack procedure in Section 3.3.2,which is aimed at measuring the overall f esc,Lyα of a specific population.

f esc,Lyα of Individual HAEs
Under the assumption of Case B recombination for T = 10 4 K, and n e ≈ 350 cm −3 gas, the Lyα escape fraction could be expressed as (Osterbrock 1989) where F Lyα,obs is the observed Lyα flux, and F Hα,int is the intrinsic Hα flux where the dust attenuation effect has been corrected.For LAEs, f esc,Lyα could be directly calculated following Equation (6) with the detected Lyα emission flux as described in Section 3.2; for non-LAEs, we use the integrated flux over ±1000 km s −1 within the r = 1 5 circular aperture.We correct the dust attenuation with the Calz and SMC laws and the optical depth in the V band derived from Beagle respectively.The differential dust attenuation factor between stars and nebular regions is set as f nebular = 1, as suggested by a number of previous studies that f nebular approaches unity at z > 2 with Hα EW >100 Å (e.g., Faisst et al. 2019).The distribution of f esc,Lyα of individual HAEs will be discussed in Section 4.1

f esc,Lyα by VLT/MUSE Stack
To calculate the typical f esc,Lyα of a specific population comprising HAEs with both Lyα detected and not detected, we stack the MUSE data cubes and measure the stacked Lyα.We shift the median-filtered subcubes of each HAE to its rest frame and then resample the subcubes into a fixed grid along the spectral axis.The rest-frame subcubes are then normalized based on their intrinsic Hα fluxes, with each voxel in the cube divided by the intrinsic Hα flux.We obtain the median-stacked cube by taking the median value of each voxel element to avoid the effect of outliers and asymmetric distribution of f esc,Lyα .The uncertainty of each voxel is estimated by bootstrapping sampling 1000 times.
We extract 1D Lyα spectra using a r = 2 5 circular aperture on the median-stacked cubes.This aperture is conservative, larger than what we use for individual HAEs, as we randomly align the subcubes during the stack without considering the asymmetry of Lyα, and thus, the stacked Lyα could be more extended than individuals.On the Lyα pseudo-narrowband image of the stack of all LAEs integrated over ±1000 km s −1 , an r = 2 5 aperture includes 95% ± 12% of the total flux.We fit skewed Gaussian profiles on the extracted 1D spectra to estimate the normalized Lyα flux.

Results
We present the results of our f esc,Lyα measurements as follows: we describe the distribution of f esc,Lyα in Section 4.1; the dependence of f esc,Lyα on various galaxy properties is discussed in Section 4.2; in Section 4.3, we show the redshift evolution of f esc,Lyα .

The Distribution of f esc,Lyα
We measure the f esc,Lyα of individual HAE following the method in Section 3.3.1.The distributions of f esc,Lyα of LAEs and non-LAEs are shown in Figure 5.For our detected LAEs, f esc,Lyα ranges from 0.04 to 1.For non-LAEs, there are negative f esc,Lyα measurements with large uncertainties; the negative  values are all consistent with zero within 3σ.The mean and median values of the individual f esc,Lyα are 0.13 ± 0.01 and 0.08 ± 0.01 (0.15 ± 0.01 and 0.11 ± 0.01) under the assumption of Calz (SMC) dust attenuation law.We model the distribution of f esc,Lyα using an exponential model, as commonly adopted to depict the distributions of Lyα EWs of Lyα-emitting galaxies (e.g., Dijkstra & Westra 2010;Dijkstra & Wyithe 2012;Kerutt et al. 2022).In fact, f esc,Lyα is proportional to the ratio of Lyα and the Hα fluxes, as shown in Equation (6); meanwhile, Lyα EWs are the ratios between Lyα fluxes and the underlying UV continua.Hα fluxes and UV continua both trace the star formation so that f esc,Lyα and Lyα EWs are correlated (e.g., Yang et al. 2016;Sobral & Matthee 2019).The distribution of f esc,Lyα is expressed as follows:  The overall likelihood for our observed f esc,Lyα distribution should be a product of the f esc,Lyα likelihoods of all galaxies.With a flat prior of f * ranging from 0 to 1, we perform Bayesian inference using the Markov Chain Monte Carlo (MCMC21 ) method to model the observed f esc,Lyα distribution.It yields the best-fit f * = 0.15 ± 0.02 (0.17 0.02 0.03 -+ ) for Calz (SMC) dust law, in agreement with the mean value of f esc,Lyα (0.13 ± 0.01 for Calz, and 0.15 ± 0.01 for SMC).The median f esc,Lyα value predicted by the best-fit exponential distribution is f ln 2 * , i.e., 0.10 ± 0.01 (0.12 0.01 0.02 -+ ) for Calz (SMC) dust law, in line with the median stack of all HAEs (0.090 ± 0.006 for Calz; 0.121 ± 0.009 for SMC).

The Dependence of f esc,Lyα on Galaxy Properties
Lyα photon transportation can be modulated by the dust and neutral gas content in galaxies, as well as the geometry and kinematics of ISM (Atek et al. 2009;Smith et al. 2019;Maji et al. 2022).We focus on nine physical properties that are most commonly studied: stellar mass (M å ), ionizing efficiency (ξ ion,0 ), redshift (z), the attenuated (observed) and the intrinsic absolute UV magnitudes (M UV,att and M UV,int ), specific star formation rate (sSFR), sSFR surface density (Σ sSFR ), the observed UV slope (β obs ), and the dust attenuation (E(B − V )).The measurements of these physical parameters are described in Section 3.1.To investigate how f esc,Lyα correlates with the properties of HAEs, we split the HAE sample into five equal-sized bins based on each physical parameter, with 32 or 33 HAEs in each bin.We then estimate the overall f esc,Lyα in each bin by stacking the Lyα.The dependence of f esc,Lyα on the nine physical parameters is presented in Figures 6 and 7. To test whether there are monotonic relations between f esc,Lyα and these physical parameters, we perform the nonparametric Spearman's correlation analysis22 using a Monte Carlo approach following Curran (2014).We include the forced measured f esc,Lyα of non-LAEs and the corresponding uncertainties in the analysis for unbiased statistics.We also test the correlations with the Kendall rank correlation analysis (Isobe et al. 1986), where the 3σ upper limits of non-LAEs' f esc,Lyα , instead of their forced measurements, are used as censored data in the analysis.As shown in Figure 6, f esc,Lyα shows a strong dependence on β obs and E(B − V ), with a Spearman's correlation coefficient (ρ) of −0.37 ± 0.04 (p-value23 0.17 ).The Kendall correlation analysis shows correlation coefficients (τ) of −0.34 ± 0.04 with p ≈ 10 −10 and −0.37 ± 0.03 with p ≈ 4 × 10 −12 for f esc,Lyα −β obs and f esc,Lyα −E(B − V ) respectively.We quote the median values of ρ, τ, and p, and their uncertainties correspond to the 16 and 84 percentiles.For M å , M UV,int , M UV,att , ξ ion , sSFR, and Σ sSFR , we illustrate their correlations with respect to f esc,Lyα Figure 7, and label the correlation coefficients of Spearman's and Kendall's tests as well as their corresponding p-values.We find that f esc,Lyα also weakly correlates with M å and M UV,int , with p of both the two correlation tests much smaller than 0.05.For M UV,att , the correlation coefficients, ρ and τ, are both close to zero with large p-values, so f esc,Lyα is almost not related to M UV,att .In the cases of ξ ion , sSFR, and Σ sSFR , the p-values of Spearman's correlation coefficients can be comparable to or exceed 0.05 within 1σ ranges.The large uncertainties in p are a result of the large uncertainties of the forced measured f esc,Lyα of non-LAEs.On the other hand, their Kendall's correlation analysis reveals τ close to zero with large p-values, indicating no correlations for these quantities.
In conclusion, our results demonstrate that f esc,Lyα has a strong dependence on β obs and E(B − V ).f esc,Lyα is also monotonically related to M å and M UV,int , but no correlations can be confirmed between f esc,Lyα and M UV,att , ξ ion , logsSFR and Σ sSFR .
Among the 165 HAEs, we note that two LAEs (HAE-60916 and HAE-87611) exhibit double-peaked Lyα profiles (Figure D1; see Appendix D for more details).They have blue UV slopes, β obs ≈ − 2.107 ± 0.007 and −1.871 ± 0.023, and high Lyα escape fractions, f esc,Lyα ≈0.58 ± 0.11 and 0.25 ± 0.04.The predicted f esc,Lyα based on their β obs using Equation ( 9) is 0.20 ± 0.14 and 0.10 ± 0.09 respectively.The observed f esc,Lyα of the two double-peaked LAEs is ∼2σ higher than the predicted values.The presence of blue Lyα peaks is likely to signal a boosted escape of Lyα photons compared to normal star-forming galaxies, consistent with previous findings on double-peaked LAEs at the cosmic noon (Furtak et al. 2022;Marques-Chaves et al. 2022).Studies on statistical samples with high S/N Lyα profiles are expected in the future using deep spectroscopy.

Redshift Evolution of f esc,Lyα
We study the evolution of f esc,Lyα as a function of redshift as shown in Figure 8.We split the HAE sample into two redshift bins, z = 5.24 (z = 4, 9-5.5) and z = 5.84 (z = 5.5-6.2), and then measured the stacked Lyα respectively.The observed f esc,Lyα is 0.104 ± 0.007 at z = 5.24 and becomes 0.086 ± 0.010 at z = 5.84.We compare our measurements with the f esc,Lyα measurements in the literature.The observed f esc,Lyα values are only about 56% and 37% of the estimates in Hayes et al. (2011) at z = 5.4 and z = 5.8 respectively.On the other hand, our measurements are generally consistent with the best-fit function in Konno et al. (2016; green dashed line in Figure 8), although they are still lower than their measurement at this redshift range (green circle in Figure 8).We note that, because Hα was not accessible at z = 5-6 until the launch of JWST, all previous measurements rely on UV LFs to estimate the SFR density and Hα luminosity density at these redshifts, while our results present a direct measurement for the first time.The offset between Hayes et al. (2011) and Konno et al. (2016) is due to the differences of the Lyα and UV luminosity limits for deriving the Lyα and UV luminosity densities from their LFs: Hayes et al.
(2011) integrated Lyα LFs from L Lyα = 10 41.3 erg s −1 to Figure 6.f esc,Lyα varying with the observed UV slope β obs and E(B − V ), under the assumption of Calz dust attenuation law.The light blue circles denote LAE individuals, and the light purple triangles are non-LAEs with forced f esc,Lyα measured from the integrated Lyα flux within ±1000 km −1 .For illustration purposes, we do not show the uncertainties of f esc,Lyα of non-LAEs.The orange stars are the overall f esc,Lyα of all HAEs derived from the Lyα stack results.The dashed lines represent the best-fit relations between f esc,Lyα and E(B − V ) or β obs based on the stack results as described in Equation (9).We show their 1σ uncertainties as shaded regions.
L Lyα = 10 42.9 erg s −1 and UV LFs down to −17.6, while Konno et al. (2016) integrated LyαLFs from L Lyα = 10 41.7 erg s −1 to L Lyα = 10 44.4 erg s −1 and UV LFs down to −17.0 mag.In addition to this, the measurements of faint-end LFs and the conversion between UV and Hα SFR could lead to large systematic uncertainties in these UV-based f esc,Lyα estimates (Konno et al. 2016;Sun et al. 2023).The observed f esc,Lyα are in good agreement with Sun et al. (2023), where f esc,Lyα is derived from Hα luminosity densities of HAEs also captured by JWST/ NIRCam WFSS.Our results, combined with those in Sun et al. (2023), suggest a tentative declining trend for the observed f esc,Lyα at z  5.The average E(B − V ) of our HAE sample is about 0.14 at z = 5.24 and 0.09 at z = 5.84, implying an intrinsically higher f esc,Lyα at the higher redshift due to their less dusty content.The decline of f esc,Lyα might be a result of the increasing neutral fraction of IGM as z increases.To completely understand the redshift evolution of f esc,Lyα , direct observations on Lyα and Hα emission lines across wide redshift spans are required in the future.

Discussion
In this section, we discuss the physical origins of the correlations as outlined in Section 4, connect the escape of Lyα and LyC photons in Section 5.1, and further qualitatively constrain ionizing photon budget for the reionziation in Section 5.2.

Implications for Lyα and LyC Escape
As illustrated in Section 4.2, f esc,Lyα correlates with M å and M UV,int , albeit to a lesser extent compared to the strong dependence of f esc,Lyα on β obs and E(B − V ).However, M å and M UV,int may be highly coupled with β obs and E(B − V ).M UV,int are the products of observed M UV,att corrected by the effect of dust attenuation, so it is highly dependent on the dust attenuation parameter adopted.We explore the coupling of f esc,Lyα −M å and f esc,Lyα −M UV,int relations with β obs and E(B − V ) in Figure 9. Comparing the expected f esc,Lyα by Equation (9) based on the median β obs and E(B − V ) in each M å and M UV,int bin, we find that the variation of f esc,Lyα across M å and M UV,int is consistent with the prediction of Equation( 9).It implies that both the f esc,Lyα −M å and f esc,Lyα −M UV,int correlations are modulated by  the underlying coupling among β obs , E(B − V ) and M å , M UV,int .Likewise, we examine the variations of f esc,Lyα versus M UV,att , sSFR, Σ sSFR , and ξ ion , and find that the trends are determined by the variations of median β obs and E(B − V ) in each stack bin.As illustrated in Figure 9, we can reproduce the trend in f esc,Lyα versus these properties by applying the median values of β obs and E(B − V ) to Equation (9).To disentangle their impacts on f esc,Lyα , it is necessary to track the evolution of f esc,Lyα within a fixed M å or M UV,int bin with varying β obs and E(B − V ).This requires a large sample with direct and accurate Lyα and Hα detection, which cannot be achieved with our current sample.
In fact, the dependencies on β obs and E(B − V ) are also found in f esc,LyC of low-redshift LyC leaking galaxies.With a compilation of LyC leakers at z ∼ 0.3, Chisholm et al. (2022) found an analytical relation between f esc,LyC and β obs :  (2016) attenuation law similar in shape to the Calz law.The modeled β obs is sensitive to the effect of dust attenuation.Furthermore, Chisholm et al. (2022) models the SED of local LyC leakers using STARBURST99 (Leitherer et al. 1999), while we perform SED modeling using BEAGLE, which is based on Bruzual & Charlot (2003); the different stellar continuum models also introduce systematics.Nevertheless, the similarity in correlation trends in f esc,Lyα −β obs and f esc,LyC −β obs is consistent with a connection between the escape processes of Lyα and LyC photons.Their escape processes are sensitive to the opacity and geometry of neutral HI gas in both the vicinity of star-forming regions and the circumgalactic medium (Dijkstra et al. 2016;Byrohl & Gronke 2020;Kimm et al. 2022;Maji et al. 2022) f esc,Lyα .We scale our measured f esc,Lyα and f esc,Lyα −β obs relation by 0.15 in Figure 10.The scaled f esc,Lyα −β obs relation is well consistent with the trend of LyC leakers within 1σ.Likewise, we estimate f esc,Lyα ξ ion as a function of β obs using the median stack f esc,Lyα and the median values ofξ ion in each stack bin, and then scale f esc,Lyα ξ ion by 0.15 as a proxy to f esc,LyC ξ ion , the production rate of ionizing photons that can escape to and ionize the IGM.As shown in Figure 10, the scaled f esc,Lyα ξ ion −β obs relation is also in good agreement with the f esc,LyC ξ ion −β obs relation in Chisholm et al. (2022).It implies that f esc,Lyα ξ ion can be used to qualitatively estimate the total emitted ionizing emissivity that impacts the IGM (see Section 5.2 for more details).Despite systematics, such as bias in sample selection and the slight difference in Lyα transmission at z = 4-5 and z = 5-6, the consistency of the scaled f esc,Lyα and f esc,LyC suggests that f esc,Lyα is likely to be proportional to f esc,LyC .With this scaling relation, the overall f esc,LyC for our HAEs are ∼0.015 ( f esc,Lyα ≈0.1) and thus negligible for the estimation of ξ ion in Equation (4).
The link between Lyα and LyC escape is also in line with the correlations we observe between f esc,Lyα and galaxy properties.If the physical conditions that favor LyC photon leakage also facilitate the escape of Lyα, observables that can indicate f esc,LyC are expected to show strong correlations with f esc,Lyα .From cosmological radiation hydrodynamics simulations, Choustikov et al. (2023) concludes that a good diagnostic of f esc,LyC should simultaneously encapsulate at least two of the following three conditions: (1) high sSFR for the production of ionizing photons, (2) the age of stellar population and timescale for stellar feedback, and (3) a proxy for neutral gas content and the ionization state of ISM.In this context, β obs and E(B − V ) serve as strong indicators for the escape of LyC and Lyα photons.E(B − V ) directly tracks the dust content along the line of sight and may also trace neutral hydrogen (e.g., Katz et al. 2022); the production and destruction of dust are also related to supernovae activities (e.g., Priestley et al. 2021).Galaxies with bluer β obs are biased toward having higher sSFRs and younger stellar populations; β obs is also very sensitive to dust attenuation.E(B − V ) and β obs satisfy at least two of the three criteria and thus are good probes of f esc,Lyα .The three physical conditions discussed above also help us to understand the relatively poor dependence of f esc,Lyα on sSFR, ξ ion , Σ sSFR , and M UV,att .sSFR or ξ ion represents just one facet of the three necessary conditions so that a high sSFR or ξ ion only could not guarantee a high f esc,Lyα .Σ sSFR alone cannot be used to reveal feedback processes or the state of the ISM along the sightline.M UV,att is a complex M å , dust, SFR, stellar population, etc., so the f esc,Lyα -M UV,att may suffer from large scatters introduced by the coupling of M UV,att with respect to these parameters.Also, M UV,att does not trace the ISM state along the sightline.Although Choustikov et al. (2023) found M UV,att can be a weak tracer of sSFR and dust attenuation at fixed M å , in which case M UV,att satisfies two out of three conditions while the stellar age would still introduce some scatters, our sample size is not large enough to disentangle on the M UV,att dependence of f esc,Lyα at a fixed M å .

Implications for Cosmic Reionization
There is an increasing body of evidences showing that dwarf galaxies are the main drivers of reionization (Atek et al. 2024;Mascia et al. 2023;Simmonds et al. 2023).However, the relative importance of massive and low-mass galaxies in driving reionization remains uncertain (Robertson et al. 2015;Finkelstein et al. 2019;Naidu et al. 2020).It is necessary to quantify the relative contribution of galaxies with different UV magnitude or masses to the total ionizing emissivity, i.e., the number of ionizing photons emitted per unit time and comoving volume.n ion  contributed by a specific population can be cast as a product of the total observed UV luminosity density ρ UV,obs , the production efficiency ξ ion , and the escape fraction of ionizing photons f esc,LyC :  (Bouwens et al. 2016;Steidel et al. 2018;Bian & Fan 2020).We adopt A UV from the best-fit SED models at the rest-frame 1500 Å.
We emphasize that the definitions of ionizing efficiency and escape fraction of ionizing photons may vary among literature (Bouwens et al. 2016;Steidel et al. 2018;Tang et al. 2019;Simmonds et al. 2023), which may cause confusion.In this section, we follow the definition in Bouwens et al. (2016).We define ξ ion as the ionizing production rate per unit intrinsic UV luminosity as described in Section 3.1.We adopt the relative escape fraction, f To provide a reference for future works that may define the ionizing efficiency and escape fraction differently, we also present the analytical relations of f esc,Lyα ξ ion with respect to M UV,att and M å in Appendix C, which can be proportional to the ionizing emissivity per unit intrinsic UV luminosity f esc,LyC ξ ion .
The contribution to the total ionizing photon budget from galaxies fainter than M UV can be expressed as follows: where M lim is the faint-end limit to compute the luminosity density by integrating over the UV LFs.We adopt the UV LF with a Schechter functional form and redshift-dependent Schechter parameters in Bouwens et al. (2022).At z = 5.5, we have the characteristic luminosity M * = − 21.02 ± 0.04, the normalization Φ * = 0.56 ± 0.04, and the faint-end slope α = − 1.90 ± 0.02.Although there are no direct observations on f esc,Lyα of faint galaxies (M UV,att > − 17 mag), we can give a qualitative constraint on the contribution of UV-faint galaxies by extrapolating the f esc,Lyα ξ ion −M UV,att analytical relation to the faint-end.With M 13 lim = - (Bouwens et al. 2016) and Equation (11), we find that galaxies fainter than −16 mag will contribute about 72% of the total ionizing budget, i.e., the accumulative emissivity fraction 0.76 where n ion, tot  is estimated by integrating from −24 mag to M lim .We present the accumulative emissivity fraction as a function of M UV,att in Figure 12.The fraction becomes 94% if M 10 lim = -.refers to the fraction of ionizing budget contributed by galaxies fainter than a specific M UV,att , as in Equation (13).We set the faint-end limit M 13 lim =and −10 for n ion,tot  , shown as the red and blue lines respectively.We mark M UV,att = − 16 as the vertical dotted line for reference.
f esc,LyC ξ ion −β obs relations measured from LyC leakers at z ∼ 0.3 (Chisholm et al. 2022).The scaling factor is consistent with the monotonic f esc,Lyα −f esc,LyC connection at z = 4-5 (Begley et al. 2023), implying that f esc,Lyα can be proportional to f esc,LyC .This may indicate that the escapes of Lyα and LyC are regulated by similar physical processes.It is also in line with the weak dependence of f esc,Lyα on M UV,att , ξ ion , sSFR, and Σ sSFR .The four parameters are not enough to depict the physical conditions that favor LyC and Lyα escape, as predicted by simulations (Choustikov et al. 2023), so they are all weak tracers of LyC or Lyα in isolation.6.We define the relative escape fraction of Lyα photons f esc,Ly rel a as f esc,Lyα divided by the escape fraction of UV photons.We fit the analytical relations of f esc,Ly rel a ξ ion with respect to M UV,att and M å , which is proportional to the ionizing emissivity per unit observed UV luminosity density.By extrapolating the f esc,Ly rel a ξ ion −M UV,att relation to the faint-end, we calculate the fraction of ionizing emissivity from UV-faint galaxies (M UV,att > − 16 mag).Qualitatively, we infer that galaxies with UV fainter than −16 mag may contribute >70% of the ionizing budget at z = 5.5.If the dependence of f esc,Lyα ξ ion on M UV,att holds until the EoR, UV-faint galaxies are predicted to be the main contributor to the reionization.Note that our HAE samples only reach M UV,att ∼ − 17 mag, so direct constraints on f esc,Lyα at the fainter end are required to better quantify the ionizing budget from faint galaxies.
JWST has shown promising prospects to unveil the mysteries of the EoR.Deeper and wider spectroscopic surveys in the future would reach the fainter end with lower stellar masses, and further probe the escape process of Lyα and LyC photons in newborn infant galaxies.
JAGUAR (Williams et al. 2018), which have certain emission features captured by NIRCam/WFSS F444W grism, and model their lines with Gaussian profiles on the mock spectra.To mimic the observational limits, we cut the mock galaxy sample by a magnitude limit of 29 for F444W and flux limit of 10 −18.3 erg s −1 cm −2 for emission lines.As illustrated in Figure A1, the colors of our sample are generally consistent with those mock HAEs and could be distinguished effectively from emitters at lower redshifts.In the left panel of Figure A1, these HAEs stay at a tightly linear region in the magnitude−Hα flux diagram, possibly suggesting main-sequence galaxies at this redshift range.They perfectly overlap with the mock HAEs, and could be distinguished easily from emitters at lower redshifts.The main possible contamination is [O III] emitters at higher redshifts, most of which could be successfully picked out by the procedures in Section 2.2.2.

Appendix B f esc,Lyα under the Assumption of SMC Dust Attenuation Law
As discussed in Section 3.1, the SMC dust attenuation law is broadly adopted in high-redshift studies.It is suggested by previous submillimeter observations (e.g., Capak et al. 2015) for typical z ∼ 5-6 galaxies.Conventionally, Calz and SMC attenuation laws are both applied and complement each other (e.g., Bouwens et al. 2016;Shivaei et al. 2018).Figure B1 illustrates f esc,Lyα of individual galaxies and from median-stack measurements based on the SMC dust law.generate a series of mock 1D spectra in a spectral sampling of 0.5 Å, with random S/Ns and continua but uniform flux.70% spectra have low S/Ns uniformly distributed from 0.1 to 3, 15% have intermediate S/Ns from three to five, and the remaining have the highest S/Ns from five to 15.We then remove the continua using a median filter of 101 pixels, the same as used in the real MUSE data, and stack all the mock spectra by taking the median value at each wavelength grid.We perform experiments for 500 times using 20, 30, 100, 200 spectra.The median values of the recovered flux for each setting reach 96%-99% of the initially set value, suggesting that a 5% underestimate is possible.For a f esc,Lyα ∼ 0.1, this systematics would simply result in an error of 0.005, smaller than the uncertainty as we report in this work.

Figure 1 .
Figure 1.The spatial and redshift distribution of Hα emitters.Left: the distribution of Hα sample in the footprints of multiband legacy data in GOODS-S.The blueshaded region is covered by MUSE observations, including the MUSE-Wide (Urrutia et al. 2019) and MUSE Hubble UDF survey (Bacon et al. 2023), color-coded by the number of overlapping data cubes (darker color indicates more data cubes).The cyan, blue, and magenta frames represent the footprints of three fields in MUSE Hubble UDF survey: the 10 hr UDF, 31 hr UDF-10, and the >100 hr MXDF; the red rectangle shows the FoV of FRESCO JWST/NIRCam F444W imaging.Colorcoded markers indicate the position of our Hα emitters.The color indicates their redshift in exactly the same scale as that on the right panel.An overdensity at z ≈ 5.3-5.4 is obvious, which has been reported in Helton et al. (2023).The crosses denote Hα emitters not covered by MUSE (41 in total); dots are within MUSE coverage (197 in total), and dots with black edges are selected sources we use for analysis (157 in total).Right: the redshift distribution of our Hα emitters.The distribution of selected Hα sources is wrapped with black edges.

Figure 3 .
Figure 3. HAE examples with their JWST/NIRcam F444W images (top first), pseudo-Lyα narrowband images (top second), and their Lyα (top third), Hα spectra (bottom).The solid box on the NB image is the 4 × 4 pixel boxy aperture used to identify Lyα and fit its profile; the dashed circle is the r = 1 5 aperture for the total flux.The top third is the boxy-aperture extracted spectrum (black) and its uncertainty (gray shaded).We show the best-fit skewed Gaussian profile in red.The redshaded region refers to the spectral window (3 pixel width) for Lyα NB.The bottom panels of each example demonstrate the 2D and 1D grism spectrum for Hα.The yellow cross marks the position of Hα emission in the 2D grism spectrum, and the best-fit Gaussian profile is presented in red.

Figure 4 .
Figure 4. Lyα stacks of MUSE data cubes for all HAEs in our sample, as well as the LAE and non-LAE subsamples.All data cubes are normalized by the intrinsic Hα flux following the Calz dust attenuation law.The psedo-Lyα narrowband (NB) images in the left cover 1213-1219 Å, i.e., slices of ±3 Å centered on 1216 Å; the continuum (Cont.for short) NB images are calculated from 1230 to 1240 Å slices and scale to the wavelength width of Lyα slices.The extracted spectra based on r = 2 5 circular apertures (black solid circle in Lyα NB images) are shown, with the best-fit skewed Gaussian profiles in red.
probability of the escape fraction being measured as f esc,Ly ˆa if the true value is x with the uncertainty f esc,Ly D a :

Figure 5 .
Figure5.The distribution of f esc,Lyα for LAEs and non-LAEs, assuming the Calz dust attenuation law.For non-LAEs, we show the forced measured f esc, Lyα based on the integrated flux over ±1000 km s −1 around the wavelength of Lyα.The mean errorbar of f esc,Lyα of LAEs and non-LAEs are marked as the blue and green horizontal lines.The thick red curve is the best-fit exponential model for the f esc,Lyα distribution; the thin red curves represent 100 exponential models with f * (Equation (7)) randomly drawn from the MCMC samples.f esc, Lyα = 0 is labeled as the gray dotted line.The red dashed line denotes the median f esc,Lyα predicted by the best-fit exponential model, and the orange dashed line is the f esc,Lyα measured from the stack of all HAEs.

Figure 7 .
Figure 7. f esc,Lyα varying with different physical parameters of HAEs assuming the Calz dust attenuation law.The light blue dots denote the f esc,Lyα measured for individual LAE; the light purple triangles are the forced measured f esc,Lyα for non-LAEs.The uncertainties of f esc,Lyα of non-LAEs are not shown for illustration purposes.The orange stars represent the stack f esc,Lyα .We label the correlation coefficients of the Spearman's correlation analysis (ρ) and Kendallʼs correlation analysis (τ) for each property, as well as the corresponding p-values.The Spearman's correlation analysis is based on the measured values and uncertainties of all individual LAEs and non-LAEs, while in the Kendallʼs tests we adopt the 3σ upper limits of f esc,Lyα for non-LAEs.

Figure 8 .
Figure 8.The redshift evolution of global f esc,Lyα from z = 0 to z = 8.Our stacked measurement of f esc,Lyα for all HAE samples is shown as a red diamond based on the assumption of Calz dust law.The orange and white diamonds are f esc,Lyα measured by stacking the HAE sample in two redshift bins (z = 5.24, 5.84).f esc,Lyα measurements compiled by Hayes et al. (2011) based on Lyα and UV (z > 2.2) or Hα luminosity densities (z < 2.2) are shown in blue circles, and the best-fit relation for their compilation is shown as the blue dashed line.The green circles denote the measurements in Konno et al. (2016), also using Lyα and UV LFs, and their best-fit relation is shown as the green dashed line.The purple squares are the f esc,Lyα derived from Lyα and Hα luminosity densities at z > 5 in Sun et al. (2023) with HAEs observed by NIRCam/WFSS.

.
The exponential index (−1.22) is in good agreement with that of our best-fit f esc,Lyα −β obs relation (-1.23) in Equation (9).Note that we only compare the Chisholm et al. (2022) f esc,LyC −β obs relation with our Calz dust law-based f esc,Lyα −β obs relation, because the β obs in Chisholm et al. (2022) are measured with the modeled stellar continuum assuming the Reddy et al.

Figure 9 .
Figure 9.The coupling of β obs and E(B − V ) on the dependence of f esc,Lyα on M å , M UV,att , M UV,att , ξ ion , sSFR, and Σ sSFR .We show the f esc,Lyα of individual LAEs as circles and forced measured f esc,Lyα of non-LAEs as triangles, colorcoded by their β obs .The stack f esc,Lyα are marked as stars color-coded by the median β obs in each stack bin.The solid and dashed lines are predicted f esc,Lyα from Equation (9) based on the median β obs and E(B − V ) in each stack bin, respectively.

Figure 10 .
Figure10.The relation between β obs and f esc,LyC , f esc,LyC ξ ion .The gray circles represent LyC leakers at z ∼ 0.3 fromChisholm et al. (2022), where β obs are measured as observed stellar continuum slope at 1550 Å.The gray downward arrows refer to the 1σ upper limits of f esc,LyC for galaxies without LyC detected, i.e., LyC flux is not significantly detected given the observed backgrounds.The gray dashed line and the shaded region are the analytical relations of f esc,LyC and f esc,LyC ξ ion with respect to β obs , and the 1σ uncertainties of the relations.The orange stars are our stack f esc,Lyα measurements assuming the Calz dust law, scaled by a factor of 0.15 according toBegley et al. (2023), and the orange dashed lines are the best-fit f esc,Lyα −β obs and f esc,Lyα ξ ion −β obs relation as described in Equation (9), also scaled by 0.15.

Figure 11 .
Figure 11.f esc,Ly rel a ξ ion as a function of M UV,att , and M å assuming the Calz dust law.The blue circles indicate f esc,Ly rel a ξ ion for each individual LAEs, and purple triangles are f esc,Ly rel a ξ ion of non-LAEs with forced measured f esc,Ly rel a .The fringed orange stars indicate the stack f esc,Ly rel a ξ ion for all HAEs.We adopt the median stack f esc,Lyα values and the median ξ ion and A UV of HAEs in each stack bin.We show the best-fit relations of the stack f esc,Lyα ξ ion of all HAEs with respect to M UV,att , M å as the orange lines.

Figure B1 .
Figure B1.Same as Figure 6 while assuming the SMC dust attenuation law.The dashed lines and gray-shaded regions represent the best-fit analytical relation as described in Equation (9).
ion on M UV,att and M å in Figure11.We fit f esc,Ly rel a ξ ion as a function of M UV,att and M å , yielding the following: