Lyα Emission at z = 7–13: Clear Evolution of Lyα Equivalent Width Indicating a Late Cosmic Reionization History

We present the evolution of Lyα emission derived from 53 galaxies at z = 6.6–13.2, which have been identified by multiple JWST/NIRSpec spectroscopy programs of Early Release Science, General Observer, Director's Discretionary Time, and Guaranteed Time Observations. These galaxies fall on the star formation main sequence and are typical star-forming galaxies with UV magnitudes of −22.5 ≤ M UV ≤ −17.0. We find that 15 out of 53 galaxies show Lyα emission at the >3σ level, and we obtain Lyα equivalent width (EW) measurements and stringent 3σ upper limits for the 15 and 38 galaxies, respectively. Confirming that Lyα velocity offsets and line widths of our galaxies are comparable to those of low-redshift Lyα emitters, we investigate the redshift evolution of the Lyα EW. We find that Lyα EWs statistically decrease toward high redshifts on the Lyα EW versus the M UV plane for various probability distributions of the uncertainties. We then evaluate neutral hydrogen fractions x H I with the redshift evolution of the Lyα EW and the cosmic reionization simulation results on the basis of a Bayesian inference framework, and obtain x H I < 0.79, = 0.62 − 0.36 + 0.15 ?> , and 0.93 − 0.07 + 0.04 ?> at z ∼ 7, 8, and 9–13, respectively. These moderately large x H I values are consistent with the Planck cosmic microwave background optical depth measurement and previous x H I constraints from galaxy and QSO Lyα damping wing absorption and strongly indicate a late reionization history. Such a late reionization history suggests that major sources of reionization would emerge late and be hosted by moderately massive halos compared with the widely accepted picture of abundant low-mass objects for the sources of reionization.


Introduction
Cosmic reionization is a crucial phase transition in the early universe.During the epoch of reionization, the neutral hydrogen in the intergalactic medium (IGM) was ionized by UV photons emitted from the first sources.Fan et al. (2006) show that reionization is mainly completed by z ∼ 6 from the measurements of Gunn-Peterson troughs (Gunn & Peterson 1965) in the spectra of quasars (QSOs), while Kulkarni et al. (2019) suggest that a broad scatter in the neutral hydrogen fraction x H I exists at redshift as low as z  5.5.Thomson scattering optical depth to the cosmic microwave background (CMB) indicates that significant reionization occurred at z ∼ 7-8 (Planck Collaboration et al. 2020).
Reionization history is investigated by deriving the redshift evolution of x H I in the IGM.Lyα emission from star-forming galaxies is a good probe of x H I because the Lyα emission line is strongly attenuated by the neutral hydrogen via the Lyα damping wing (Ouchi et al. 2020).While Lyα emitters (LAEs) have been commonly found at z ; 6 (e.g., Stark et al. 2011;De Barros et al. 2017), LAEs at higher redshift z  7 are rare (e.g., Stark et al. 2010;Treu et al. 2013;Hoag et al. 2019;Jung et al. 2020).Although reionization history is constrained by many studies of LAE clustering analysis (e.g., Ouchi et al. 2018;Umeda 2023), Lyα luminosity function (e.g., Ouchi et al. 2010;Konno et al. 2014;Hu et al. 2019;Wold et al. 2022), and Lyα equivalent width (EW) distribution (e.g., Mason et al. 2018;Hoag et al. 2019;Bolan et al. 2022), the uncertainties of reionization history are still large at z ∼ 6-8.At a higher redshift of z  8, almost no observational constraints existed before the arrival of the James Webb Space Telescope (JWST; Ferruit et al. 2022;Jakobsen et al. 2022;Böker et al. 2023;Gardner et al. 2023;Rieke et al. 2023;Rigby et al. 2023).JWST allows a deep survey of faint galaxies during the epoch of reionization.In fact, many LAEs are spectroscopically identified at a high redshift of 7  z  11 (e.g., Bunker et al. 2023b;Jung et al. 2023;Tang et al. 2023;Chen et al. 2024;Jones et al. 2024;Saxena et al. 2024).In addition, recent studies of JWST observations place the constraints on reionization history at z  8 (Hsiao et al. 2023;Curtis-Lake et al. 2023;Umeda et al. 2023), although the constraints are not strong.The open questions regarding reionization are thus when reionization started and at what pace reionization proceeded.
Lyα EW, defined by a Lyα flux relative to a UV continuum flux, is a key observable property that reflects the intensity of the Lyα emission.The fraction of LAEs (EW > 25Å) among Lyman-break galaxies (LBGs), the Lyα fraction, is measured as a way to put a constraint on the Lyα optical depth (e.g., Stark et al. 2010;Ono et al. 2012;Schenker et al. 2012;Treu et al. 2012;Pentericci et al. 2018;Jones et al. 2024).The Lyα fraction dramatically decreases from z ∼ 6 to 7, indicating the rapid increase in Lyα optical depth at z  6. Treu et al. (2012Treu et al. ( , 2013) ) suggest that the full Lyα EW distribution, including galaxies with low Lyα EW and no Lyα detections, contains more information about the Lyα optical depth than the Lyα fraction, which is used to evaluate the fraction of galaxies detected at some fixed threshold of Lyα EW.The x H I values in the IGM are not directly constrained by the Lyα optical depth because Lyα photons are absorbed by the neutral hydrogen not only in the IGM, but also in the interstellar medium (ISM) and circumgalactic medium (CGM).The effects of the ISM and CGM on Lyα transmission have been studied by some models or simulations (e.g., Dijkstra et al. 2011;Mesinger et al. 2015).Mason et al. (2018) introduce a flexible Bayesian framework based on that of Treu et al. (2012Treu et al. ( , 2013) ) and the simulations by Mesinger et al. (2015) to infer x H I from LBGs with detections and no detections of Lyα emission.
The source of reionization is another important feature of reionization.Recent observational studies report a steep faintend slope α of a UV luminosity function (Bouwens et al. 2017;Livermore et al. 2017;Atek et al. 2018;Ishigaki et al. 2018;Finkelstein et al. 2019), indicating that faint sources make dominant contributions to UV luminosity correlated with ionizing photon production.Although the escape fraction of ionizing photons f esc is also crucial to understanding ionizing sources, it is impossible to directly measure f esc at z > 4, due to the high optical depth of ionizing photons (e.g., Vanzella et al. 2018).In addition, recent works of simulations and semianalytical models show conflicting results for f esc correlated with halo mass (e.g., Gnedin et al. 2008;Wise & Cen 2009;Sharma et al. 2016) and anticorrelated with halo mass (e.g., Paardekooper et al. 2015;Faisst 2016;Kimm et al. 2017).The f esc value during the epoch of reionization thus remains uncertain.However, given that low-mass (massive) galaxies contribute to early (late) reionization (e.g., Faisst 2016;Dayal et al. 2024), we may be able to constrain the mass dependence of reionization source from reionization history.
There are three scenarios of the reionization history suggested by Naidu et al. (2020), Ishigaki et al. (2018), and Finkelstein et al. (2019), referred to as the late scenario, medium-late scenario, and early scenario, respectively.The three scenarios differ in the escape fraction of ionizing photons f esc and the faint-end slope α of the UV luminosity function.We refer to Model II described in Naidu et al. (2020) as the late scenario because Model II self-consistently explains observational results of the redshift evolution of f esc (e.g., Siana et al. 2010;Rutkowski et al. 2016;Marchi et al. 2017;Steidel et al. 2018;Fletcher et al. 2019).In the late scenario, the faintend slope should be shallow (α > −2) with the assumption that f esc depends on the star formation rate (SFR) surface density Σ SFR .Bright and massive galaxies (M UV < −18 and stellar mass M * > 10 8 M e ) are major contributors to reionization in this scenario.In the medium-late scenario, Ishigaki et al. (2018) find the steep faint-end slope (α  −2) from observations and the constant escape fraction ( f esc ∼ 0.17) with which the scenario can mostly explain the observational measurements.Galaxies with various luminosities contribute to reionization in this scenario.In the early scenario, the faint-end slope should be steep (α  −2) with the assumption that f esc anticorrelates with halo mass M h .Faint and low-mass galaxies (M UV > −15 and M *  10 6 M e ) are major contributors to reionization in this scenario.
In this study, we present the evolution of Lyα emission with 53 galaxies at z  7 whose redshifts are spectroscopically confirmed by JWST/NIRSpec.Combining these high-redshift galaxies and the full distributions of Lyα EW with detections and upper limits, we infer x H I at z ∼ 7-13.The x H I values provide insights into reionization history and sources.This paper is constructed as follows.Section 2 describes the JWST/NIRSpec spectroscopic data and our sample of galaxies.We conduct spectral fittings to our sample galaxies and investigate the Lyα properties obtained from the best-fit parameters in Section 3. In Section 4, we explain our method to infer x H I and present the estimated x H I values.In Section 5, we discuss reionization history and sources.Section 6 summarizes our findings.Throughout this paper, we use a standard ΛCDM cosmology with Ω Λ = 0.7, Ω m = 0.3, and H 0 = 70 km s −1 Mpc −1 .All magnitudes are in the AB system (Oke & Gunn 1983).

Data and Sample
The spectroscopic data sets used in this study were obtained in multiple public observation programs: the Early Release Science (ERS) observations of GLASS (ERS-1324, PI: T. Treu;Treu et al. 2022)  Lützgendorf; Bunker et al. 2023a;Eisenstein et al. 2023).The GLASS data were taken at one pointing position with high resolution (R ∼ 2700) filter-grating pairs of F100LP-G140H, F170LP-G235H, and F290LP-G395H covering the wavelength ranges of 1.0-1.6,1.7-3.1, and 2.9-5.1 μm, respectively.The total exposure time of the GLASS data is 4.9 hr for each filtergrating pair.The CEERS data were taken at 8 and 6 pointing positions with the prism (R ∼ 100) and medium-resolution grating (R ∼ 1000), respectively.The prism covered the wavelength of 0.6-5.3μm.The filter-grating pairs of medium-resolution grating observations were F100LP-G140M, F170LP-G235M, and F290LP-G395M, whose wavelength coverages were 1.0-1.6,1.7-3.1 and 2.9-5.1 μm, respectively.The total exposure time of the CEERS data is 0.86 hr for each filter-grating pair per pointing.The GO-1433 and DDT-2750 observations, each at a single pointing, were conducted with the prism, and the total exposure times are 3.7 and 5.1 hr, respectively.The JADES data were obtained at 3 pointing positions with the prism and four filter-grating pairs of F070LP-G140M (0.7-1.3 μm), F170LP-G235M, F290LP-G395M, and F290LP-G395H.The total exposure times of the JADES observations are 9.3 hr for the prism and 2.3 hr for each filter-grating pair per pointing.Nakajima et al. (2023) and Harikane et al. (2024) have reduced the GLASS, CEERS, GO-1433, and DDT-2750 data with the JWST pipeline version 1.8.5 with the Calibration Reference Data System context file of jwst_1028.pmapor jwst_1027.pmapwith additional processes improving the flux calibration, noise estimate, and the composition, and have measured systemic redshifts z sys with optical emission lines (Hβ, [O III]λ4959, and [O III] λ5007).See Nakajima et al. (2023) and Harikane et al. (2024) for further details of the data reduction and z sys determination.The combination of the data of Nakajima et al. (2023) and Harikane et al. (2024) from the public observational programs provides a sample of 196 galaxies at 4.0 z sys 11.4.We then select 35 galaxies with z sys 7 from the sample that is referred to as the public data sample.The JADES data are publicly available,7 and have already been reduced with the pipeline developed by the ESA NIRSpec Science Operations Team and the NIRSpec GTO Team.The JADES data consists of 178 galaxies at 0.7 z sys 13.2, including 75 galaxies whose redshifts are undetermined.We select 18 galaxies from the JADES data with z sys 6.5 identified by Bunker et al. (2023a) with emission lines (e.g., Hβ, [O III]λ4959, and [O III]λ5007) or spectral breaks.To these 18 galaxies, we add GN-z11 at z sys = 10.6, which is also observed as part of the JADES program (GTO-1181, PI: D. Eisenstein; Eisenstein et al. 2023).We refer to a total of these 19 galaxies as the JADES data sample.Note that we use the values of the rest-frame Lyα EW, Lyα velocity offset, and line width presented in the literature (Bunker et al. 2023b).Combining the JADES data sample with the public data sample, we obtain a total of 54 galaxies whose z sys fall in the range of 6.6 z sys 13.2.Here, 6, 26, 1, 2, and 19 out of all 54 galaxies are found in the GLASS, CEERS, GO-1433, DDT-2750, and JADES data, respectively.We note that these 54 galaxies in the public data and JADES data samples are observed with the gratings and/or prism.The public data sample (35 galaxies) consists of 17, 26, and 8 galaxies that have the grating, prism, and grating+prism spectra, respectively.The JADES data sample (19 galaxies) has both the grating and prism spectra, while the grating spectra at the Lyα wavelengths are missing in three galaxies.There are a total of 33 (=17 + 19 −3) galaxies whose Lyα emission wavelengths are observed with the gratings.The rest of the objects, 21 (=54-33) galaxies, have only the prism spectra, where the measurements of the Lyα velocity offset and line width are not obtained due to the poor spectral resolution.
In Figures 1 and 2, we show the spectra of 53 galaxies with the exception of GN-z11.We note that because CEERS_01027 is observed at two pointing positions with the mediumresolution grating, two spectra of CEERS_01027 (CEERS_P4M_01027 and CEERS_P10M_01027) are shown in Figure 1.
In Figure 3, we compare the stellar mass M * and SFR of our galaxies with those of the star formation main sequence (MS) to confirm whether our galaxies are typical star-forming galaxies or not.We adopt the stellar mass and SFR reported by Nakajima et al. (2023) and Harikane et al. (2024) for 14 and 5 galaxies of the public data sample, respectively.We use the stellar mass and SFR reported by Curti et al. (2024), Harikane et al. (2024), andBunker et al. (2023b) for 11, 3, and 1 galaxies of the JADES data sample, respectively.For 11 galaxies for which Nakajima et al. (2023) only report the stellar mass estimates, we calculate the SFR from the UV magnitude.We first derive the UV luminosity L ν (UV) from the UV magnitude.We then convert the UV luminosity into the SFR with SFR = 1.15 × 10 −28 L ν (UV) (Madau & Dickinson 2014).We correct the SFR for the dust extinction, using an extinction −UV slope β UV relation (Meurer et al. 1999) and β UV -M UV relation (Bouwens et al. 2014).In Figure 3, we do not present the other nine galaxies whose stellar mass is not determined.The stellar masses and SFRs of our galaxies are comparable with the star formation MSs at z ∼ 6 (Santini et al. 2017), indicating that our galaxies are typical star-forming galaxies.

Spectral Fitting
To obtain the Lyα emission properties of galaxies, we fit a continuum+line model to the Lyα emission line.This model is the linear combination of power-law aλ − β and Gaussian functions, where a, β, A, λ cen , and σ are the free parameters.This model is multiplied by where τ(λ) is the optical depth of the IGM to the Lyα photons at the observed wavelength λ (Inoue et al. 2014).We then convolve the model with the line spread functions (LSF) of JWST/ NIRSpec measured with a spectrum of a planetary nebula (Isobe et al. 2023a) to account for instrumental smoothing.
We determine the free parameters of the model, using emcee (Foreman-Mackey et al. 2013) to conduct Markov Chain Monte Carlo (MCMC) simulations.We apply this fitting method to the prism spectra of galaxies in the public data sample.For both the grating and prism spectra of galaxies in the JADES data sample, we add the full width at halfmaximum (FWHM) for the LSF as an additional free parameter.Because the spectral resolution of the JADES data sample is higher than those of the other spectra due to the compact morphologies of the JADES targets (Jones et al. 2024), differences in the FWHM for the LSF have an influence on the fittings.The signal-to-noise ratio of the continua in the grating spectra of our galaxies is worse than that in the prism spectra, which do not allow us to simultaneously determine the amplitude a and the slope β of the power law.We thus fix the slope of the power law to zero for the grating spectra under the assumption of a flat continuum ( f ν = const.),which is typical for star-forming galaxies (Stark et al. 2010), and determine the rest of the model parameters.We show examples of the posterior distributions of the fitting parameters in Figure 4. We determine the model parameter and 1σ uncertainty by the mode (i.e., a peak of the posterior distribution) and 68th percentile highest posterior density interval (HPDI; i.e., the narrowest interval containing 68%) of the posterior distribution, respectively, which are good representatives even for a skewed posterior distribution.We are not able to conduct spectral fittings to the grating spectra of JADES_10013905 and GLASS_100005 due to the small numbers of spectral pixels on the redder side than the observed Lyα wavelengths.As for JADES_10013905, we use the prism spectrum for the fitting.
In summary, we conducted spectral fittings to the grating spectra of 30 (=33 −2 −1) galaxies (not including GNz-11) and the prism spectra of 22 (=21 + 1) galaxies (see Section 2).Note that we have detected continua for 52 (=30 + 22) galaxies, including tentative detections.We confirm that if we use galaxies with continua detected at the >3σ level for the analysis below, our results do not change significantly.
With the public data and JADES data samples, we find that 11 out of 31 galaxies with the grating spectra and four out of 22 galaxies with the prism spectra have Lyα detections, where we place a threshold of the signal-to-noise ratio of >3σ for Lyα detection.In total, we find 15 (=11 + 4) out of 53 galaxies with Lyα detections in our sample.We note that Lyα emission of CEERS_01029 is not detected beyond our threshold, while Tang et al. (2023) claim the detection of Lyα emission.However, we confirm that the Lyα EW measurements of Tang et al. (2023), EW 0,Lyα = 4.2 ± 1.3 Å, is within the 3σ upper limits of the Lyα EW in this work, 9.8 Å.
Table 1 summarizes the properties of all 54 galaxies in our sample.Table 2 compares the rest-frame Lyα EW and Lyα velocity offset of the galaxies, with Lyα detections measured in this study with those in the literature.For CEERS_01027, we show mean values of the Lyα redshift, EW, velocity offset, and line width measured from two grating spectra (Section 2) because the measurements from the two spectra are consistent within errors.

Lyα Velocity Offset
For the galaxies with Lyα emissions detected from the grating spectra, we calculate Lyα redshifts z Lyα with 1 + z Lyα = λ cen /λ Lyα , where λ Lyα is the rest-frame Lyα wavelength, 1215.67Å.We then derive the Lyα velocity offset Δv Lyα, which is defined by where c is the speed of light.In Figure 5, we plot the Δv Lyα values measured in this study and the literature.The Δv Lyα values of our galaxy are comparable with not only z > 6 galaxies, but also z ∼ 2-3 galaxies.We confirm that there is no significant redshift evolution of the Lyα velocity offset.For the same galaxies, we compare our measurements with those of Saxena et al. (2024) and Tang et al. (2023; see also Table 2).The measurements of Δv Lyα for JADES_00016625, JADES_00004297, JADES_10013682, and JADES_00021842 are consistent with those in Saxena et al. (2024).While the measurements of Δv Lyα for CEERS_00698, CEERS_01027, and CEERS_01019 appear to be lower than those of Tang et al. (2023), in Figure 5, they are comparable, considering the slightly lower systemic redshifts reported by Tang et al. (2023) and the uncertainties of the measurements (Table 2).The range and mean value of the velocity offset in our sample are 168-334 and 234 ± 76 km s −1 , respectively.

Lyα Line Profile
In order to investigate the Lyα line profile of our galaxies, we create composite spectra by conducting mean stacking for four (three) galaxies at z ∼ 7 and 8 that have medium-resolution grating spectra.We determine the redshift bin so that the numbers of the galaxies at each bin are comparable.The redshift range and mean redshift of the galaxies at z ∼ 7 and 8 are 6.6 < z < 7.3 (7.4 < z < 8.0) and 〈z〉 = 7 and 8, respectively.Note that at 〈z〉 = 8, we stack four spectra for three galaxies because we have two spectra of CEERS_01027 (Section 2).We do not include the spectrum of CEERS_01019, which has signatures of active galactic nuclei (AGN; Larson et al. 2023;Isobe et al. 2023b) in the composite spectra at 〈z〉 = 8 because the Lyα line width could be broadened by AGN.We note that we do not include CEERS_01019 in the analysis of the Lyα line profile but include it in the analysis of Lyα EW, fraction, and x H I .For stacking individual spectra at z ∼ 7 and 8, we shift the galaxy spectra into the mean redshifts of 〈z〉 = 7 and 8 based on the systemic redshift, respectively.We then obtain the composite spectra, taking averages of the fluxes normalized with the continuum fluxes of individual spectra.We calculate 1σ errors of the composite spectra by taking averages of the error in spectra of individual galaxies.The composite spectra at 〈z〉 = 7 and 8 are shown in Figure 6.We conduct spectral fittings to the composite spectra in the same manner as in Section 3.1.We calculate the intrinsic FWHMs of the Lyα emission lines from the individual spectra and composite spectra by Note that we have already corrected for instrumental broadening of line profiles by conducting convolution with the LSF of JWST/NIRSpec in the fittings (Section 3.1).The FWHM values of the composite spectra at 〈z〉 = 7 and 8 are FWHM = 149 ± 68 and 219 ± 83 km s −1 , respectively.In Figure 7, we compare the FWHMs measured in this study with those taken from the literature.The FWHMs of galaxies at z ∼ 7 and 8 in this work are comparable with those at z ∼ 7 (Ouchi et al. 2010;Endsley et al. 2022) and lower redshift z ∼ 2-3 (Hashimoto et al. 2017;Kerutt et al. 2022).We find no large evolution of FWHMs from z ∼ 2-8, which is consistent with that between z = 5.7 and 6.6 (Ouchi et al. 2010).In Figure 8, we compare the composite spectra at 〈z〉 = 7 with those at 〈z〉 = 8.While we find no large evolution of Lyα line profiles, we find an evolution of the emission line flux.Because the composite spectra are normalized with their continuum fluxes, the evolution of the emission line flux suggests the evolution of the Lyα EW.

Lyα EW
We obtain Lyα fluxes F Lyα by the equation, Ly where l max (l min ) is the maximum (minimum) value of the observed wavelength for the Lyα emission.We estimate flux errors from the MCMC results on the basis of the error For the galaxies with no F Lyα measurements (i.e., no Lyα detections), we estimate the upper limits of the Lyα fluxes.For each galaxy, we conduct Monte Carlo simulations with 100 mock spectra that are random realizations of the observed spectrum with the error spectrum.We conduct the model fitting to the 100 mock spectra to obtain Lyα flux measurements in the same manner as in Section 3.1, and determine the 1σ error of the Lyα flux with the distribution of the Lyα flux measurements of the 100 mock spectra.From these Lyα flux errors, we calculate the 3σ upper limits of EW 0,Lyα for the galaxies with no Lyα detections.We compare our measurements of EW 0,Lyα with those in the literature (Tang et al. 2023;Chen et al. 2024;Jones et al. 2024;Saxena et al. 2024) in Table 2 and find good agreement except for JADES_10013682, JADES_00004297, CEERS_80374, and CEERS_80239.For JADES_10013682, the EW 0,Lyα value measured in this study (31.5 ± 10.0 Å) is lower than those in  (2.7 ± 0.4) × 10 −18 erg s −1 cm −2 , is consistent with that reported by Saxena et al. (2024), (2.7 ± 0.3) × 10 −18 erg s −1 cm −2 .We thus confirm that the differences in the EW 0,Lyα measurements of JADES_10013682 (JADES_ 00004297) come from the continuum fluxes obtained from the different resolutions of the grating and prism spectra, although our Lyα EW measurement is explained by the 2σ (3σ) error of those in Saxena et al. (2024).The impacts of the choice of the grating and prism spectra on Lyα EW measurements are discussed in Appendix A. For CEERS_80374, we do not detect Lyα emission at >3σ significance level and measure only 3σ upper limits of EW 0,Lyα , <86.7 Å, while Chen et al.In Appendix B, we investigate the effects of the potential flux loss of the Lyα halo on our results.These different EW 0,Lyα measurements from those in the literature for a small number of galaxies may not affect our major results.
In Figure 9, we plot EW 0,Lyα as a function of M UV .To investigate the redshift evolution of EW 0,Lyα , we divide our sample into three subsamples by redshift so that each subsample has at least 10 galaxies and a redshift range of Δz  1.Our subsamples consist of 22, 12, and 17 galaxies whose redshift ranges are 6.5 z < 7.5, 7.5 z < 8.5, and 8.5 z < 13.5, referred to as z ∼ 7, 8, and 9-13, respectively.Note that we removed JADES_00008115 and CEERS_01163 from our subsamples because their M UV values are undetermined.We perform linear function fittings to our subsamples combined with samples of galaxies taken from the literature (Ono et al. 2012;Schenker et al. 2012;Oesch et al. 2015;Song et al. 2016;Shibuya et al. 2018;Matthee et al. 2019;Tilvi et al. 2020;Endsley et al. 2022;Jung et al. 2022;Kerutt et al. 2022; see Appendix C of Jones et al. 2024), which do not overlap with the galaxies in this study.For the 3σ upper limits of EW 0,Lyα , we adopt two different distributions in the fittings.First, we adopt normal distributions with a mean value of 0 and standard deviations of the 1σ limits of EW 0,Lyα .Second, we adopt uniform distributions between the 0 and 3σ upper limits of EW 0,Lyα .We first fit the galaxies at z ∼ 7 to determine the slope and the intercept.We then use the fixed slope at z ∼ 7 and determine the intercept for the galaxies at z ∼ 8 and 9-13.The best-fit linear functions are presented in the top panel of Figure 9.There is no significant difference between the results using normal distributions and uniform distributions.We also conduct fittings based on the Kendall correlation test, using the cenken function in the NADA library from the R-project statistics package built based on the results of Akritas et al. (1995) and Helsel (2005).We use the cenken function because it computes the Kendall rank correlation coefficient and an associated linear function for censored data, such as a set of data, including upper limits.The fitting results using the cenken function are shown in the bottom panel of Figure 9.In both the top and bottom panels of Figure 9, we find that EW 0,Lyα becomes lower at higher redshift, which is a clear signature of the redshift evolution of x H I with observational data alone.The evolution of Lyα EW is prominent between z ∼ 8 and 9-13, suggesting that reionization significantly proceeded between z ∼ 8 and 9-13.The p-values for the hypothesis that there is no correlation between M UV and EW 0,Lyα at z ∼ 7, 8, and 9-13 are p = 0.58, 0.98, and 0.96, respectively.These p-values indicate that there is no clear correlation between M UV and EW 0,Lyα at each redshift bin.

Lyα Fraction
To evaluate the ionizing state of the IGM, we calculate the Lyα fraction.We define the Lyα fraction a X Ly EW th by the ratio of the number of galaxies with strong Lyα emission (EW > EW th ) to the number of total galaxies.The Lyα fraction a X Ly 25 is typically investigated for galaxies with −21.75 < M UV < −20.25 and −20.25 < M UV < −18.75 (e.g., Stark et al. 2011;Ono et al. 2012;Schenker et al. 2012Schenker et al. , 2014;;Pentericci et al. 2018;Jones et al. 2024).We only have a small number of galaxies with −20.25 < M UV < −18.75 at z ∼ 8 and 9-13, and thus investigate the Lyα fraction at a rebinned redshift of z ∼ 7 and    (5, 6, 9) 8-13.For 12 (9) galaxies with −20.25 < M UV < −18.75 at z ∼ 7 (8-13), we find 3 (0) galaxies with EW > 25 Å.We obtain = a -+ X 0.25 Ly 25 0.13 0.19 (3/12) and <0.19 (0/9) at z ∼ 7 and 8-13, respectively.We estimate the 1σ error and upper limit, using the statistics for small numbers of events by Gehrels (1986).As shown in the top panel of Figure 10, our a X Ly 25 value at z ∼ 7 is consistent with those in the literature (Ono et al. 2012;Schenker et al. 2012Schenker et al. , 2014;;Pentericci et al. 2018;Jones et al. 2024).We note that the denominators of X Ly 25 at 6 z 7, which is due to the poor Lyα transmission in the moderately neutral IGM.At z ∼ 8-13, we only obtain an upper limit of a X Ly 25 due to no galaxy with EW > 25 Å.Because the Lyα EW becomes lower at higher redshift, as shown in Figure 9, we need a lower threshold of EW than 25 Å to detect Lyα emitting galaxies at z  8.We thus calculate the Lyα fraction 0.11 0.09 0.21 (1/9) at z ∼ 7 and 8-13, respectively.In the bottom panel of Figure 10, we find the trend of a decrease in the Lyα fraction at higher redshift between z ∼ 7 and 8-13.In summary, combining the results of a X Ly 25 in this work and the literature with those of a X Ly 10 in this work, we find that the Lyα fraction monotonically decreases from z ∼ 6 to 8-13 because of the more neutral hydrogen in the IGM at the earlier epoch of reionization.As for the x H I estimation, the full distribution (i.e., the probability distribution function) of Lyα EW allows us to place stronger constraints on x H I than on X Lyα (Treu et al. 2012(Treu et al. , 2013;;Mason et al. 2018), as described in Section 1.We thus mainly focus on x H I estimated with the probability distribution function (PDF) of Lyα EW in Section 4.

Estimates of the Neutral Hydrogen Fraction
We estimate neutral hydrogen fractions x H I , comparing the rest-frame Lyα EW values (simply referred to as EW in this section) obtained from observations with those of the theoretical EW distribution models developed by Dijkstra et al. (2011).As shown in Figure 11, EW values at the low x H I universe are higher than those at the high x H I universe, due to the Lyα damping wing absorption made by the IGM.In the EW distribution models, Dijkstra et al. (2011) quantify the PDF of the fraction of Lyα photons transmitted through the IGM T IGM , combining galactic outflow models with largescale seminumeric simulations of reionization.There are two parameters, the neutral hydrogen column density N H I and the outflow velocity v wind , in the galactic outflow models.The EW distribution models are constructed under the assumption that the IGM at z = 6 is completely transparent to Lyα photons, and that the PDF of EW changes only by the evolution of the ionization state of the IGM.It is also assumed that the EW distribution at z = 6 is described by an exponential function, ), where EW c is a scale factor of EW at z = 6.The EW distributions during the epoch of reionization are computed with where N is a normalization factor, and p(T IGM ) is a PDF of T IGM , which is computed in Dijkstra et al. (2011).As a fiducial model, Dijkstra et al. (2011) use a set of the models with N H I = 10 20 cm −2 , v wind = 200 km s −1 , and EW c = 50 Å, in which T IGM is computed at z ∼ 8.6.For comparison, we show the models of EW distribution with N H I = 10 20 cm −2 , v wind = 25 and 200 km s −1 , and EW c = 30, 40, and 50 Å in Figure 11.The EW distribution does not largely change for different v wind values, but moderately changes for different EW c values.Because Dijkstra et al. (2011) find the relation of Δv Lyα ∼ 2v wind in the galactic outflow model, the Lyα velocity offset of the models with v wind = 200 km s −1 corresponds to 400 km s −1 .We note that the lower Lyα velocity offsets in our sample of Δv Lyα = 234 ± 76 km s −1 (mean values; Section 3.2) than that of the models are unlikely to affect our estimates because v wind does not largely change the EW distribution, as  a The EW 0,Lyα values and 3σ upper limits of these galaxies are measured with the prism spectra.b These galaxies are used to construct the composite spectra (Section 3.3).c Because the grating spectra have a small number of spectral pixels on the redder side of the observed Lyα wavelengths used to determine the continuum fluxes, we are not able to determine the upper limits of Lyα EW.As for JADES_10013905, we adopt the measurement from the prism spectrum.b The EW 0,Lyα values and 3σ upper limits of these galaxies are measured with the prism spectra.c Although Chen et al. (2024) show high EW 0,Lyα values for these galaxies, we do not detect Lyα emission for CEERS_80374 and obtain a lower EW 0,Lyα for CEERS_80239.The deviations in the measurements may come from the aperture size used to extract the spectra of these galaxies.et al. 2003;Pentericci et al. 2011;Vanzella et al. 2011;Willott et al. 2013;Maiolino et al. 2015;Oesch et al. 2015;Stark et al. 2015;Willott et al. 2015;Furusawa et al. 2016;Knudsen et al. 2016;Pentericci et al. 2016;Carniani et al. 2017;Laporte et al. 2017;Mainali et al. 2017;Stark et al. 2017;Pentericci et al. 2018;Hashimoto et al. 2019;Endsley et al. 2022; see Endsley et al. 2022 and Table 4 therein    UV magnitude and obtain the scale factor of Å, which corresponds to about 30 Å for the M UV values of our sample.We thus use the sets of models with N H I = 10 20 cm −2 , v wind = 200 km s −1 , and EW c = 35 ± 5 Å, which are comparable to the observed values.
Figure 12 presents the M UV distribution as a function of redshift for our sample.We find that there is a trend that faint galaxies are included for our galaxies in the low-redshift bins.To avoid systematics introduced by bias in UV magnitudes, we select galaxies, removing 8 (1) galaxies at z ∼ 7 (8) fainter than characteristic luminosity with the criterion of z < 7 or M UV > −18.We conduct two sample Kolmogorov-Smirnov tests between the UV magnitudes for the galaxies of our subsamples (Section 3.4) before and after the selection.Note that the redshift range of our subsamples for z ∼ 7 after the selection is 7 z < 7.5.We confirm that the UV magnitude distributions are not similar between the galaxies at z ∼ 7 and ∼ 8 before the selection at the 99% significance level, but consistent for all three redshift bins after the selection at the 99% significance level.We estimate x H I for our subsamples of z ∼ 7, 8, and 9-13, which consist of 14 (=22 −8), 11 (=12 −1), and 17 galaxies, respectively (Section 3.4).
Here, we need to be careful about the effects of the slit loss on the Lyα EW measurements because we compare the Lyα EW obtained from JWST/NIRSpec MSA observations with the theoretical models of Lyα EW distributions based on the Very Large Telescope/FORS2 observations (Pentericci et al. 2018) to estimate x H I .Our sample galaxies are observed with an NIRSpec MSA microshutter, which could miss the fluxes of extended Lyα halos due to the small area of the microshutters (0 20 × 0 46) (Jakobsen et al. 2022) compared to the FORS2 slit (∼1″ × 8″) (Pentericci et al. 2018).The slit loss effects on Lyα emission are discussed in recent studies.While Jung et al. (2023)   Rest-frame Lyα EW as a function of absolute UV magnitude.The blue, green, and orange circles (arrows) show the measurements (3σ upper limits) of the galaxies at z ∼ 7, 8, and 9-13 of our sample, respectively.The light blue, light green, and light orange symbols are the same as the blue, green, and orange symbols, but for the galaxies obtained from the literature (Ono et al. 2012;Schenker et al. 2012;Oesch et al. 2015;Song et al. 2016;Shibuya et al. 2018;Matthee et al. 2019;Tilvi et al. 2020;Endsley et al. 2022;Jung et al. 2022;Kerutt et al. 2022; see Appendix C of Jones et al. 2024).The top and bottom panels present the results of the linear and cenken function fittings, respectively.In the top panel, the solid (dashed) lines denote the bestfit linear functions obtained by adopting a normal (uniform) distribution for the EW 0,Lyα upper limits.See Section 3.4 for the details of the fittings.values obtained from the galaxies, including those confirmed via the photometry in the literature (Stark et al. 2011;Ono et al. 2012;Schenker et al. 2012Schenker et al. , 2014;;Pentericci et al. 2018;Mason et al. 2019).Bottom: the red star marks show the a X Ly 10 values obtained in this work.
four sources measured from NIRSpec MSA are consistent with those measured from ground-based observations.Jones et al.
(2024) also show that the Lyα EW measurements from NIRSpec MSA are not likely to affect their results.In addition, our a X Ly 25 value obtained from NIRSpec MSA observations is consistent with that obtained from FORS2 observations (Pentericci et al. 2018), suggesting that slit loss effects on Lyα EW measurements may not be significant.However, we need to quantitatively assess the slit loss effects on Lyα EW measurements to estimate x H I accurately.As done by Tang et al. (2024), we conducted a forward-modeling analysis to compare the Lyα EW obtained from NIRSpec MSA with those from FORS2.The details of the analysis are described in Appendix B. The results show that Lyα EW measured with NIRSpec MSA is 72% ± 8% of those with FORS2.We thus use the corrected Lyα EWs divided by this fraction (72% ± 8%) as fiducial values to estimate x H I .
To estimate the probability distributions of x H I , we perform a Bayesian inference based on that of Mason et al. (2018).We note that our Bayesian method only sets x H I as a free parameter, instead of x H I and M UV in the method of Mason et al. (2018) because our subsamples at each redshift bin have similar distributions of UV magnitude as described above.Under the Bayesian method, the posterior probability distribution p(x H I |EW) is based on the equation, where erfc(x) is the complementary error function for x.  0.93 0.07 0.04 at z ∼ 7, 8, and 9-13, respectively.We present a 3σ upper limit for x H I at z ∼ 7.For comparison, we also estimate x H I with two sets of models and measurements.One is the combination of the EW distribution models of Dijkstra et al. (2011) and the uncorrected EWs.The other is the combination of the EW distribution models of Mason et al. (2018, see their Figure 7) with the corrected EWs.As described above, the z ∼ 6 EW distribution in Mason et al. (2018) is based on an LAE observation with FORS2 (De Barros et al. 2017;Pentericci et al. 2018).We thus estimate x H I by combining the corrected EWs with the models of Mason et al. (2018).In Figure 14, we present the results with the three different sets of models and EWs.The x H I estimates with uncorrected EWs are modestly overestimated compared to those with corrected EWs.In the case where we use corrected EWs, there are little differences between the x H I estimates with EW distribution models of Dijkstra et al. (2011) and those of Mason et al. (2018).We thus adopt the x H I estimates with the EW distribution models of Dijkstra et al. (2011) and corrected EWs as fiducial results.

Discussion
In Figure 15, we present our constraints on the redshift evolution of x H I at z = 7-13.For comparison, we also show the  x H I estimates taken from the literature.Our x H I values are consistent with those estimated with Lyα damping wing absorption of QSOs/gamma-ray bursts (GRBs), Lyα luminosity function, Lyα EW of LBGs, and the Thomson scattering optical depth of the CMB within the errors.While these results constrain x H I up to z ∼ 7.5, which is limited to the redshifts of the identified QSOs and GRBs (Totani et al. 2014;Wang et al. 2019;Matsuoka et al. 2023), a recent study by Umeda et al. (2023) extends the x H I measurements to z ∼ 8-12 with JWST galaxies via UV continuum modeling of the Lyα damping wing absorption features.The sample of Umeda et al. (2023) consists of 26 galaxies, which are used in their Lyα damping wing analysis, while our sample of 53 galaxies includes all 26 of Umeda et al.ʼs galaxies.Although there is an overlap of the samples between our study and Umeda et al.ʼs study, our study focuses on Lyα line EW statistics compared with the numerical simulations, which complement the study of Umeda et al. (2023) in the independent method to estimate x H I .We confirm that the x H I estimates of our study and Umega et al.ʼs study are consistent within the errors.This consistency solidifies the validity of our measurements of Lyα EW, which is used to infer the x H I values.We plot the late, medium-late, and early scenarios suggested by Naidu et al. (2020), Ishigaki et al. (2018), and Finkelstein et al. (2019), respectively (Section 1).Our results are consistent with the late and medium-late scenarios.
As explained in Section 1, the late and medium-late scenarios suggest that sources of reionization form at the late epoch, probably at z ∼ 7-8.If the late scenario is correct, the slopes of UV luminosity functions should be as shallow as α > −2 under the assumption that f esc depends on the SFR surface density Σ SFR (Section 1).However, recent observational results have preferred a steep slope of α < −2 (Bouwens et al. 2017;Livermore et al. 2017;Atek et al. 2018;Ishigaki et al. 2018;Finkelstein et al. 2019).Because the requirement of the shallow UV slope does not agree with the observational results, this discrepancy could possibly suggest that the escape fraction does not depend on Σ SFR but on the other quantities.For example, if the escape fraction positively correlates with halo masses (e.g., Gnedin et al. 2008;Wise & Cen 2009;Sharma et al. 2016), the positive correlation could lead to late reionization with the steep UV slope suggested by the UV luminosity function observations.If we consider the possibility of star-forming galaxies as major sources of reionization, ionizing photons are prevented from escaping to the IGM because of the enrichment of metals and dust in star-forming galaxies.To account for the high escape fractions of starforming galaxies hosted by moderately massive halos, channels in the ISM produced by the starbursts and frequent supernovae (e.g., Paardekooper et al. 2015;Steidel et al. 2018) or unknown physical processes may be needed.In fact, Paardekooper et al. (2015) show that young star-forming galaxies in massive halos would have a high escape fraction if supernova feedback causes either an inhomogeneous density or the removal of the dense gas around the galaxies.Another possibility may be that galaxies in massive halos harbor AGNs that efficiently produce ionizing photons.While QSOs (i.e., bright AGNs) make negligible contributions to reionization due to their low number density (Onoue et al. 2017;McGreer et al. 2018;Parsa et al. 2018;Wang et al. 2019;Jiang et al. 2022), some studies suggest that the contributions from faint AGNs to reionization may be significant (Giallongo et   Salvaterra 2015; Giallongo et al. 2019).In fact, many faint AGNs have been found at z = 4-10 by JWST observations (Harikane et al. 2023;Kocevski et al. 2023;Maiolino et al. 2023aMaiolino et al. , 2024;;Übler et al. 2023;Matthee et al. 2024).Harikane et al. (2023) and Maiolino et al. (2023a) suggest that the contributions from faint AGNs (M UV  −22) to the ionizing photons are not very large but do exist.Such many faint AGNs may have caused the late reionization.In summary, our results, together with recent observational results, suggest that the escape fraction may positively correlate with halo masses and that possible major sources of reionization are starforming galaxies or faint AGNs hosted by moderately massive halos.In fact, the galaxies with high neutral hydrogen column densities, which disturb the escape of ionizing photons, are found at z  8 (D' Eugenio et al. 2023;Heintz et al. 2023;Umeda et al. 2023), indicating that significant reionization is caused by the objects that emerge at the late era of reionization (z ∼ 7-8) and are hosted by moderately massive halos.

Summary
In this paper, we present our constraints on x H I from Lyα EWs of galaxies at z ∼ 7-13.We select 54 galaxies at 6.6 < z < 13.2 observed with JWST/NIRSpec in the multiple programs of ERS, DDT, GO, and GTO.The redshifts of the galaxies of our sample are spectroscopically confirmed by emission lines or spectral breaks.Because we are not able to conduct spectral fitting for one galaxy, we investigate the Lyα properties of 53 galaxies.We obtain the Lyα velocity offset, line width, and EW for 15 galaxies with Lyα detections from the spectral fittings.For 38 galaxies with no Lyα detections, we measure the 3σ upper limits of Lyα EW.Our major findings are summarized below: 1. We find no redshift evolution of the Lyα velocity offset and line width from z ∼ 2-3 to ∼ 7-13.While the composite spectra at mean redshifts of 〈z〉 = 7 and 8 also show no evolution of line profiles, we find a clear attenuation of line flux, suggesting the evolution of the Lyα EW.To investigate the redshift dependence of the Lyα EW, we conduct linear function fittings for galaxies at z ∼ 7, 8, and 9-13 by adopting various probability distributions for the upper limits of Lyα EW.We find that EW 0,Lyα becomes lower at higher redshift (Figure 9), which is a clear signature of the redshift evolution of x H I based on the observational data alone.The Lyα fraction measurements at z ∼ 7 also suggest the evolution of x H I between z ∼ 6 and 7. 2. To estimate x H I independent of UV magnitude, we select 42 out of 53 galaxies with z 7 and M UV −18 that have similar distributions of M UV for each redshift bin.
We estimate x H I by comparing our EW 0,Lyα values with those of the EW distribution models (Dijkstra et al. 2011) based on a Bayesian inference (Mason et al. 2018).We correct our EW 0,Lyα values by forward modeling the slit loss effects for comparison with ground-based results.
The estimates are x H I < 0.79, = -+ 0.62 0.36 0.15 , and -+ 0.93 0.07 0.04 at z ∼ 7, 8, and 9-13, respectively.Our x H I values are consistent with those estimated with the Lyα damping wing absorption of QSOs/GRBs, Lyα luminosity function, Lyα EW of LBGs, and the Thomson scattering optical depth of the CMB within errors.Our x H I values are also consistent with those estimated in a complementary approach with the UV continua of bright galaxies for damping absorption estimates (Umeda et al. 2023).The redshift evolution of our x H I values is consistent with the late or medium-late reionization scenarios.Such a late reionization history suggests that major reionization sources are the objects hosted by moderately massive halos, which is in contrast to the widely accepted picture of abundant faint low-mass objects for reionization sources.measurements with the prism spectra in our analysis thus would not affect our results significantly.

Appendix B Forward Modeling
As described in Section 4, we conduct a forward-modeling analysis to estimate the slit loss effects on Lyα EW measurements.We first obtain the intrinsic UV and Lyα spatial profile of 10 galaxies at z ∼ 5-6 based on measurements with MUSE (Leclercq et al. 2017).For UV continuum emission, we construct a 2D exponential profile with a scale length of a core rs cont and a total UV flux obtained from the UV magnitude.For Lyα emission, we construct a 2D exponential profile for two components, which are a core component and a halo component.The profile of the core (halo) component is obtained from a scale length of the core rs cont (halo rs halo ) and a total core (halo) Lyα flux F cont (F halo ).We then convolve the UV and Lyα profiles with the point-spread function (PSF) of NIRSpec and FORS2.We assume that the PSF is approximated by a 2D Gaussian profile.The PSF of NIRSpec is obtained from the WebbPSF package (Perrin et al. 2014).We adopt a seeing size of ∼0 8 (Pentericci et al. 2018)    aperture size of 0 3 (Section 3.4) along the spatial direction centered on the spatial peak position of the slit (0 20 × 0 46; Jakobsen et al. 2022).For FORS2, we measure the fluxes by summing the fluxes within the slit (1″ × 8″; Pentericci et al. 2018).We vary the offset of the source from the center of the slit along the diagonal line of the NIRSpec MSA slit for both the NIRSpec MSA and FORS2 slits.From the flux measurements with NIRSpec and FORS2, we obtain the Lyα EW of 10 galaxies at z ∼ 5-6.We calculate 1σ Lyα EW errors by conducting Monte Carlo simulations of 100 forward modelings with the errors of rs cont , rs halo , F cont , and F halo .We present the fractions of Lyα EW measured within the NIRSpec MSA (FORS2) slit to the Lyα EW measured with the total Lyα and UV fluxes, and fractions of Lyα EW measurements with NIRSpec to those with FORS2 for different source positions in the left and right panels of Figure B2, respectively.We confirm that the fractions do not depend on the source's position in the slit.We thus adopt the fraction of Lyα EW measurements with NIRSpec to those with FORS2 where the source is in the center of the slit, 72% ± 8%, as a correction factor in Section 4.

Figure 1 .
Figure 1.Spectra of the public data sample.Each panel shows the two-dimensional spectrum (top) and the one-dimensional spectrum (bottom).The black solid lines and shaded regions represent the observed spectra and associated 1σ errors, respectively.The vertical red-dashed lines represent the rest-frame Lyα wavelengths of 1215.67 Å at the systemic redshifts.The yellow regions indicate the detected Lyα lines whose signal-to-noise ratio is larger than 3σ (Section 3.1).

Figure 2 .
Figure 2. Same as Figure 1, but for the JADES data sample.
(2024) report a Lyα detection with = a -+ EW 205 0,Ly 27 48 Å.For CEERS_80239, our measurements (105.3 ± 72.1 Å) are lower than those of Chen et al. (2024) ( - + 334 62 109 Å).These differences are due to the aperture used to extract the spectra.Our spectra are produced via the summation of 3 pixels (0 3) along the spatial direction centered on the spatial peak position to minimize the effects coming from the noisy regions close to the edge (Nakajima et al. 2023), while the spectra of Chen et al. (2024) are extracted with the aperture, which is typically 6 pixels.Because the apertures of our spectra are narrower than those of Chen et al. (2024), we may miss the extended fluxes of Lyα halos.Our measurements thus show no Lyα detection or a smaller EW 0,Lyα value than that of Chen et al. (2024).

Figure 3 .
Figure 3. SFR as a function of stellar mass.The red squares and orange circles present our sample galaxies whose SFRs are measured in this work and the literature (Bunker et al. 2023b; Nakajima et al. 2023; Curti et al. 2024; Harikane et al. 2024), respectively.The blue solid line and shaded region indicate the star formation MS at z ∼ 6 and its uncertainty (Santini et al. 2017), respectively.

Figure 4 .
Figure 4.The posterior probability distribution functions of the fitting parameters.The top (bottom) panel presents the fitting results for a grating (prism) spectrum of CEERS-P10M_01027 (CEERS_00439).The red solid lines and black-dashed lines show the mode and the boundaries of the 68th percentile HPDI, respectively (Section 3.1).
d While Tang et al. (2023) claim Lyα detection for this galaxy, we did not detect Lyα emission beyond our threshold of the signal-to-noise ratio of >3σ.The EW 0,Lyα value of Tang et al. (2023) is within the 3σ upper limits of EW 0,Lyα in this work.shown in the top panel of Figure 11.Recent observational studies show column densities of neutral gas of high-z galaxies (Heintz et al. 2023; Umeda et al. 2023), which are N H I ∼ 10 19 -10 23 cm −2 .The observed values of N H I are comparable to those of the models.Dijkstra et al. (2011) chose a scale factor of EW c = 50 Å because the Lyα fraction of LAEs with EW > 75 Å estimated with the EW c value is about 0.2, which corresponds to the median value observed at z ∼ 6 (Stark et al. 2010).However, Pentericci et al. (2018) claim that a scale factor of a best-fit exponential function that matches the observations at z ∼ 6 is EW c = 32 ± 8 Å. Mason et al. (2018) parameterize the z ∼ 6 EW distribution (De Barros et al. 2017; Pentericci et al. 2018) as an exponential function plus a delta function, which depend on

Figure 5 .
Figure 5. Lyα velocity offset as a function of absolute UV magnitude.The red star marks and diamonds show our galaxies and GN-z11 at z = 10.6 (Bunker et al. 2023b), which are included in our sample, respectively.The white and red triangles indicate galaxies identified with ground-based telescopes at z ∼ 2-3 (Erb et al. 2014) and z > 6 (Cubyet al. 2003;Pentericci et al. 2011;Vanzella et al. 2011;Willott et al. 2013;Maiolino et al. 2015;Oesch et al. 2015;Stark et al. 2015;Willott et al. 2015;Furusawa et al. 2016;Knudsen et al. 2016;Pentericci et al. 2016;Carniani et al. 2017;Laporte et al. 2017;Mainali et al. 2017;Stark et al. 2017;Pentericci et al. 2018;Hashimoto et al. 2019;Endsley et al. 2022; seeEndsley et al. 2022 and Table 4 therein), respectively.The red circles and squares denote the measurements of Tang et al. (2023) for the CEERS galaxies and Saxena et al. (2024) for the JADES galaxies, respectively, all of which are included in our sample.For the same galaxies, our measurements (star marks) are connected to the Tang et al. measurements (circles) or the Saxena et al. (2024) measurements (squares) with the black solid lines.The gray-dashed lines represent the empirical relation at z = 7 suggested by Mason et al. (2018).
Figure 5. Lyα velocity offset as a function of absolute UV magnitude.The red star marks and diamonds show our galaxies and GN-z11 at z = 10.6 (Bunker et al. 2023b), which are included in our sample, respectively.The white and red triangles indicate galaxies identified with ground-based telescopes at z ∼ 2-3 (Erb et al. 2014) and z > 6 (Cubyet al. 2003;Pentericci et al. 2011;Vanzella et al. 2011;Willott et al. 2013;Maiolino et al. 2015;Oesch et al. 2015;Stark et al. 2015;Willott et al. 2015;Furusawa et al. 2016;Knudsen et al. 2016;Pentericci et al. 2016;Carniani et al. 2017;Laporte et al. 2017;Mainali et al. 2017;Stark et al. 2017;Pentericci et al. 2018;Hashimoto et al. 2019;Endsley et al. 2022; seeEndsley et al. 2022 and Table 4 therein), respectively.The red circles and squares denote the measurements of Tang et al. (2023) for the CEERS galaxies and Saxena et al. (2024) for the JADES galaxies, respectively, all of which are included in our sample.For the same galaxies, our measurements (star marks) are connected to the Tang et al. measurements (circles) or the Saxena et al. (2024) measurements (squares) with the black solid lines.The gray-dashed lines represent the empirical relation at z = 7 suggested by Mason et al. (2018).

Figure 6 .
Figure 6.Lyα emission lines in the composite spectra of our sample galaxies.Left: the top and bottom axes present the velocity offset and observed wavelength.The y-axis corresponds to the fluxes normalized with the continuum fluxes.The blue solid line and shaded regions represent the composite spectrum for the 〈z〉 = 7 galaxies and associated 1σ errors, respectively.The gray solid lines show the individual spectra whose Lyα lines are redshifted to the mean redshift.The vertical black-dashed line denotes the Lyα wavelength in the observed frame defined with the mean redshift.Right: same as the left panel, but for the 〈z〉 = 8 galaxies.The red solid line indicates the composite spectrum for the 〈z〉 = 8 galaxies.

Figure 7 .
Figure 7. Lyα line width as a function of absolute UV magnitude.The large blue circles and large red squares represent the FHWMs measured from the composite spectra at 〈z〉 = 7 and 8, respectively.The small blue circles and small red squares denote the FHWMs of individual galaxies in our sample at z ∼ 7 and 8, respectively.The red double squares show the FWHMs of CEERS_01019 and GN-z11 (Bunker et al. 2023b) in our sample, both of which have signatures of AGN (Bunker et al. 2023b; Larson et al. 2023; Isobe et al. 2023b).The blue and white triangles indicate the FWHMs of the galaxies at z ∼ 7 (Ouchi et al. 2010; Endsley et al. 2022) and 2-3 (Hashimoto et al. 2017; Kerutt et al. 2022), respectively.

Figure 8 .
Figure 8. Evolution of Lyα line profiles.The top and bottom axes represent the velocity offset and rest-frame wavelength.The y-axis corresponds to the fluxes normalized with the continuum fluxes.The blue and red solid lines represent the composite spectra of the 〈z〉 = 7 and 8 galaxies, respectively.The lightshaded regions show the 1σ errors of the composite spectra.The vertical blackdashed line denotes the rest-frame Lyα wavelength.
claim that the Lyα flux of one source measured from NIRSpec MSA is only 20% of those measured from MOSFIRE, Tang et al. (2023) find that the Lyα EWs of

Figure 9 .
Figure9.Rest-frame Lyα EW as a function of absolute UV magnitude.The blue, green, and orange circles (arrows) show the measurements (3σ upper limits) of the galaxies at z ∼ 7, 8, and 9-13 of our sample, respectively.The light blue, light green, and light orange symbols are the same as the blue, green, and orange symbols, but for the galaxies obtained from the literature(Ono et al. 2012;Schenker et al. 2012;Oesch et al. 2015;Song et al. 2016;Shibuya et al. 2018;Matthee et al. 2019;Tilvi et al. 2020;Endsley et al. 2022;Jung et al. 2022;Kerutt et al. 2022; see Appendix C ofJones et al. 2024).The top and bottom panels present the results of the linear and cenken function fittings, respectively.In the top panel, the solid (dashed) lines denote the bestfit linear functions obtained by adopting a normal (uniform) distribution for the EW 0,Lyα upper limits.See Section 3.4 for the details of the fittings.

Figure 10 .
Figure 10.Lyα fraction as a function of redshift.Top: the red star mark and black circles indicate the a X Ly 25 values obtained from the galaxies, all of which are spectroscopically confirmed in this work and Jones et al. (2024), respectively.The red arrow denotes the same as the red star mark, but for the upper limit of a X Ly 25 .The white circles represent the a X Ly 25values obtained from the galaxies, including those confirmed via the photometry in the literature(Stark et al. 2011;Ono et al. 2012;Schenker et al. 2012Schenker et al. , 2014;;Pentericci et al. 2018;Mason et al. 2019).Bottom: the red star marks show the EW i represents EW measurements of individual galaxies.Here, p(EW i |x H I ) is a likelihood function, and p(x H I ) is a uniform prior with 0 x H I 1.The likelihood functions are calculated, including the measured 1σ errors σ i by the equation, EW|x H I ) represents the EW distribution model as a function of x H I .For galaxies with no Lyα detections, we instead use the likelihood functions, including the 3σ upper limits of EW lim given by

Figure 11 .
Figure 11.Probability distribution functions of the rest-frame Lyα EW.The purple, blue, green, and orange solid lines indicate the EW distribution models (Dijkstra et al. 2011) with a neutral hydrogen column density N H I = 10 20 cm −2 , an outflow velocity v wind = 200 km s −1 , a scale factor EW c = 50 Å, and neutral hydrogen fractions x H I = 0.13, 0.60, 0.80, and 0.91, respectively.In the top panel, the dashed lines represent the same as the solid lines, but for the models with v wind = 25 km s −1 .In the bottom panel, dashed (dashed-dotted) lines denote the same as the solid lines, but for the models with EW c = 40 (30) Å.

Figure 12 .
Figure12.Absolute UV magnitude distribution of our sample galaxies as a function of redshift.The blue, green, and orange-filled (open) circles denote galaxies with (no) Lyα detections at z ∼ 7, 8, and 9-13, respectively.The gray solid lines represent the selection criteria for our sample (Section 4).

Figure 14 .
Figure 14.Comparison of the redshift evolution of x H I .The red star marks, black circles, and white squares represent x H I estimates obtained with the Dijkstra et al. (2011) models+corrected Lyα EW, the Dijkstra et al. (2011) models+uncorrected Lyα EW, and the Mason et al. (2018) models+corrected Lyα EW, respectively.The purple solid, blue-dotted, and orange-dashed lines represent the reionization scenarios suggested by Naidu et al. (2020), Ishigaki et al. (2018), and Finkelstein et al. (2019), respectively.The light-shaded regions show the uncertainties of the reionization scenarios.

Figure 15 .
Figure 15.Redshift evolution of x H I .The red star marks denote fiducial x H I estimates obtained in this work.The purple solid, blue-dotted, and orange-dashed lines represent the reionization scenarios suggested by Naidu et al. (2020), Ishigaki et al. (2018), and Finkelstein et al. (2019), respectively.The light-shaded regions show the uncertainties of the reionization scenarios.The black squares, circles, and diamonds present the x H I estimates derived from Lyα damping wing absorption of GRBs (Totani et al. 2006, 2014), QSOs (Schroeder et al. 2013; Davies et al. 2018; Greig et al. 2019; Wang et al. 2020), and LBGs (Curtis-Lake et al. 2023; Hsiao et al. 2023; Umeda et al. 2023), respectively.The gray triangles and circles indicate the x H I estimates using an LAE clustering analysis (Ouchi et al. 2018; Umeda 2023) and Lyα luminosity function (Ouchi et al. 2010; Konno et al. 2014; Inoue et al. 2018; Morales et al. 2021; Goto et al. 2021; Ning et al. 2022; Umeda 2023), respectively.The gray diamonds denote the x H I estimates obtained from the Lyα EW distribution of LBGs (Hoag et al. 2019; Mason et al. 2019; Jung et al. 2020; Whitler et al. 2020; Bruton et al. 2023; Morishita et al. 2023).The white triangles and circles represent the x H I estimates derived from the Gunn-Peterson through of QSOs (Fan et al. 2006) and Lyα and Lyβ forest dark fractions of QSOs (McGreer et al. 2015; Zhu et al. 2022; Jin et al. 2023), respectively.The white pentagon shows the x H I estimates obtained from the CMB observations under the assumption of instantaneous reionization (Planck Collaboration et al. 2020).
as the FWHM of the PSF of FORS2.In Figure B1, we present an example of the forward modeling for a source at z ∼ 5.For NIRSpec, we measure Lyα and UV fluxes by summing the fluxes within the

Figure A1 .
Figure A1.Comparison of Lyα EW measurements with our grating spectra and those with our prism spectra.The left panel presents a wide range of the figure, while the right panel zooms in to show the plots more clearly.The circles denote Lyα EWs for sources observed with both the grating and prism.The error bars and arrows represent the 1σ errors and 3σ upper limits, respectively.The black line indicates that Lyα EWs obtained from grating spectra are equivalent to those obtained from prism spectra.

Figure B1 .
Figure B1.Example of the forward modeling of a source at z ∼ 5. From left to right in the top row, each panel indicates the intrinsic Lyα spatial profile, PSF of NIRSpec, Lyα emission convolved with the PSF of NIRSpec, PSF of FORS2, and Lyα emission convolved with the PSF of FORS2.The bottom row shows the same as the top row but for UV emission.The radius of the black-dashed and red-dashed circles correspond to the scale lengths of the core and halo components, respectively.The white lines represent the NIRSpec MSA three-shutter slits (middle panel) and FORS2 slits (rightmost panel).

Figure B2 .
Figure B2.The fraction of Lyα EW as a function of the source position.Left: the red (blue) circles represent the fractions of Lyα EW measured within the NIRSpec MSA (FORS2) slit to Lyα EW measured with total Lyα and UV fluxes.The left, middle, and right gray-dashed lines indicate the corner, center, and corner of the slit, respectively.Right: the black circles denote the fractions of Lyα EW measurements with NIRSpec to those with FORS2.The gray-dashed lines show the same as the top panel.

Table 1
Sample in This Study