JWST spectroscopy of z ∼ 5–8 UV-selected galaxies: new constraints on the evolution of the Ly α escape fraction in the reionization era Open Access

ABSTRACT

We describe JWST/NIRSpec prism measurements of Ly α emission in z ≳ 5 galaxies. We identify Ly α detections in 10 out of 69 galaxies with robust rest-optical emission-line redshift measurements at 5 ≤ z < 7 in the Cosmic Evolution Early Release Science (CEERS) and DDT-2750 observations of the Extended Groth Strip field. Galaxies at z ≃ 6 with faint continuum (F150W=27–29 mag) are found with extremely large rest-frame Ly α equivalent widths (EWs; ranging up to 286 Å). Likely Ly α detections are also seen in two new z > 7 galaxies (z = 7.49 and 7.17) from the second epoch of CEERS observations, both showing large Ly α EWs that likely indicate significant transmission through the intergalactic medium (IGM). We measure high Ly α escape fractions in the 12 Ly α emitters in our sample (median 0.28), two of which show |$f_{\rm esc}^{ {\rm Ly}\alpha }$| near unity (>0.80). We find that |$50_{-11}^{+11}$| per cent of z ≃ 6 galaxies with [O iii] + H β EW>1000 Å have |$f_{\rm esc}^{ {\rm Ly}\alpha }$| >0.2, consistent with the fractions found in lower redshift samples with matched [O iii] + H β EWs. While uncertainties are still significant, we find that only |$10_{-5}^{+9}$| per cent of z > 7 galaxies with similarly strong rest optical emission lines show such large |$f_{\rm esc}^{ {\rm Ly}\alpha }$|⁠, as may be expected if IGM attenuation of Ly α increases towards higher redshifts. We identify photometric galaxy overdensities near the z ≳ 7 Ly α emitters, potentially providing the ionizing flux necessary to create large ionized sightlines that facilitate Ly α transmission. Finally, we investigate the absence of Ly α emission in a comparable (and spectroscopically confirmed) galaxy overdensity at z = 7.88 in the Abell 2744 field, discussing new prism spectra of the field obtained with the UNCOVER program.

1 INTRODUCTION

The reionization of intergalactic hydrogen marks the important stage where early structure formation has impacted the vast majority of baryons in the Universe. Through a variety of observational efforts over the past two decades, substantial progress has been made to constrain the timeline of this reionization process and the nature of the early ionizing sources (see Stark 2016; Robertson 2022; Fan, Bañados & Simcoe 2023 for reviews). Measurements of hydrogen absorption features presented in the quasar spectra indicate that the universe is partially neutral at z ≃ 7 (Greig et al. 2017; Davies et al. 2018; Wang et al. 2020; Yang et al. 2020a) and becomes mostly ionized by z ≃ 5.5–6 (e.g. Fan et al. 2006; McGreer, Mesinger & D’Odorico 2015; Yang et al. 2020b; Qin et al. 2021; Zhu et al. 2022; Jin et al. 2023). Measurements of the cosmic microwave background have suggested a similar picture, where constraints on the electron scattering optical depth imply a reionization mid-point at z = 7.8 ± 0.7 (Planck Collaboration XLVII 2016; Planck Collaboration VI 2020). With the rapidly declining quasar luminosity function at z > 3 (e.g. Matsuoka et al. 2018; Kulkarni, Worseck & Hennawi 2019; Jiang et al. 2022), young star-forming galaxies are thought to be the likely dominant contributors to the ionizing photons necessary for this process (Robertson et al. 2015; Stanway, Eldridge & Becker 2016; Dayal & Ferrara 2018; Finkelstein et al. 2019; Naidu et al. 2020).

The Ly α emission from early star-forming galaxies has been shown to be another useful probe in constraining the timeline of reionization (e.g. Dijkstra 2014; Ouchi, Ono & Shibuya 2020). While strong Ly α emission has been commonly observed in z ≃ 6 galaxies (e.g. Stark, Ellis & Ouchi 2011; Curtis-Lake et al. 2012; Cassata et al. 2015; De Barros et al. 2017), deep spectroscopic surveys show that it becomes increasingly rare at higher redshifts (z ≃ 7–8; e.g. Fontana et al. 2010; Stark et al. 2010; Treu et al. 2013; Caruana et al. 2014; Pentericci et al. 2014; Tilvi et al. 2014; Hoag et al. 2019; Mason et al. 2019; Jung et al. 2020). This decline of Ly α visibility is consistent with a significant neutral hydrogen fraction of the intergalactic medium (IGM) at z ≳ 7 (X_HI ≳ 0.6; e.g. Mesinger et al. 2015; Zheng et al. 2017; Mason et al. 2018a, 2019; Hoag et al. 2019; Jung et al. 2020; Whitler et al. 2020; Bolan et al. 2022).

JWST has recently ushered in a new era of high-redshift Ly α studies (e.g. Jones et al. 2023; Jung et al. 2023; Nakane et al. 2023; Napolitano et al.2024; Saxena et al. 2023a, b; Tang et al. 2023; Witten et al. 2023), with detections now being made out to z ≃ 11 (Bunker et al. 2023). Ly α investigations with JWST provide several advantages compared to earlier efforts. The absence of atmospheric OH lines greatly increases Ly α completeness and improves the reliability of upper limits on fluxes when lines are not detected. Access to strong rest-optical emission lines provides spectroscopic redshifts whether or not Ly α is detected, allowing the environment around Ly α emitters to be more efficiently characterized. Detection of H β (or H α) allows the intrinsic luminosity of Ly α to be calculated, which when combined with the observed Ly α flux yields estimates of the Ly α escape fraction (e.g. Ning et al. 2023; Roy et al. 2023; Saxena et al. 2023b; Tang et al. 2023,, 2024; Lin et al. 2024; Tang et al. 2024 in preparation). The rest-optical emission lines also constrain the systemic redshifts for measurement of Ly α velocity profiles, providing a key input to reionization calculations (e.g. Bunker et al. 2023; Prieto-Lyon et al. 2023; Saxena et al. 2023b; Tang et al. 2023).

The Cosmic Evolution Early Release Science (CEERS;¹ ERS-1345, PI: S. Finkelstein; Finkelstein et al. in preparation, also see Finkelstein et al. 2022, 2023) has recently obtained JWST imaging and spectroscopy of the Extended Groth Strip (EGS) field, providing one of first opportunities for a statistical investigation of Ly α in the reionization era. This region is known to have a large number of z ≳ 7 Ly α emitting galaxies, with previously confirmed structures at z = 7.48, z = 7.73, and z = 8.68 (Oesch et al. 2015; Zitrin et al. 2015; Roberts-Borsani et al. 2016; Stark et al. 2017; Tilvi et al. 2020; Jung et al. 2022; Larson et al. 2022; Cooper et al. 2023). The CEERS spectra have revealed additional galaxies at these redshifts (Arrabal Haro et al. 2023b; Fujimoto et al. 2023; Harikane et al. 2024; Tang et al. 2023; Sanders et al. 2024), potentially suggesting the Ly α emitters trace large-scale overdensities spanning several physical Mpc. Such overdensities of galaxies could generate significant amounts of ionizing photons, creating large ionized bubbles that facilitate the transmission of Ly α photons from the member galaxies (e.g. Barkana & Loeb 2004; Furlanetto, Zaldarriaga & Hernquist 2004; Wyithe & Loeb 2005; Iliev et al. 2006; Dayal & Ferrara 2018; Weinberger et al. 2018).

The CEERS spectra have also revealed insight into the nature of the z ≳ 7 Ly α emitters in the field. Tang et al. (2023) demonstrated that the systems with Ly α are atypical in their properties, with hard spectra and more efficient ionizing photon production than is common at this redshift. This may reflect very young stellar populations formed during a rapid burst (Endsley et al. 2021b; Saxena et al. 2023b; Tang et al. 2023) or the presence of an active galactic nucleus (AGN; Larson et al. 2023). Regardless of the origin, the ionizing properties will result in boosted intrinsic Ly α luminosities, enhancing the likelihood of detecting the line in the face of IGM attenuation even within small ionized bubbles. The Ly α detections at z ≳ 7 also appear to be redshifted significantly from systemic, with velocity offsets commonly in excess of 500 km s⁻¹ (Bunker et al. 2023; Tang et al. 2023). Given that the IGM attenuation of Ly α is strongest close to the line centre, the large velocity offsets will further contribute to the visibility of Ly α.

In spite of this progress, the size of the ionized bubbles around the z ≳ 7 Ly α emitting galaxies in the EGS field remains poorly constrained, making it unclear the relative role of overdensities and intrinsic galaxy properties play in explaining their strong observed Ly α (e.g. Jung et al. 2022; Tang et al. 2023; Whitler et al. 2023a). If overdensities have carved out large ionized regions (i.e. ≳1 physical Mpc), we expect the transmission of Ly α through the IGM to be enhanced (e.g. Mason et al. 2018b; Mason & Gronke 2020). If the overdensities and bubble sizes are smaller, we would expect to see significant attenuation in the Ly α line from the damping wing of the neutral IGM. The CEERS H β detections allow Ly α escape fractions (⁠|$f_{\rm esc}^{ {\rm Ly}\alpha }$|⁠) to be derived for the Ly α emitters in the field. Typical escape fractions are low (⁠|$f_{\rm esc}^{ {\rm Ly}\alpha }$| = 0.03–0.09), suggesting the majority of Ly α photons are not making their way through the NIRSpec slitlets (Tang et al. 2023). Galaxies with similar rest-optical spectral properties at lower redshifts (z ≃ 0.3–3) tend to have much stronger Ly α (e.g. Yang et al. 2017; Tang et al. 2021; Flury et al. 2022a), with Ly α equivalent widths (EWs) that are 6–12× greater than the z ≳ 7 Ly α emitters in CEERS. This may be indicating that the z ≳ 7 Ly α emitters trace relatively small ionized bubbles, with excess attenuation from the neutral IGM on the outside of the bubbles. Alternatively, the interstellar medium (ISM) and circumstellar medium (CGM) of the z ≳ 7 galaxies may scatter the Ly α more than the comparison samples at lower redshift.

In this paper, we seek to extend the Ly α investigations presented in Tang et al. (2023) to a larger sample of galaxies. We have two primary goals. First, we seek to characterize the Ly α properties of galaxies at slightly lower redshifts (5 ≤ z < 7) where the impact of the IGM on the line is reduced. At these redshifts, we expect to see Ly α with larger escape fractions than those studied to date at z ≳ 7. This comparison sample will be important in interpreting the emerging body of Ly α detections deep in the reionization era. Secondly, we aim to use new measurements of Ly α in faint z ≳ 7 galaxies to better characterize whether the previously known Ly α emitters in the EGS field are likely to trace large ionized regions. We will investigate whether there are significant overdensities of photometrically selected galaxies in the EGS field, as may be expected if large bubbles are present.

The connection between ionized bubbles and galaxy overdensities has been further tested in recent JWST observations of the z = 7.88 protocluster discovered in GLASS (ERS-1324, PI: Treu; Treu et al. 2022), JWST Director’s Discretionary Time (DDT-2756, PI: W. Chen; Roberts-Borsani et al. 2023), and GO-1840 (PI: Hashimoto; Hashimoto et al. 2023) observations of the Abell 2744 field (Morishita et al. 2023). The galaxies in this protocluster are significantly overdense (20×) over a relatively small radius of 60 physical kpc. While such a region may be expected to have carved out a large ionized bubble, JWST spectroscopy presented in Morishita et al. (2023) has revealed no Ly α emission. In this paper, we use new JWST spectra of galaxies in the z = 7.88 protocluster to understand why Ly α does not appear to be enhanced in this overdense volume.

The organization of the paper is as follows. In Section 2, we first describe JWST/NIRSpec observations (Section 2.1), our Ly α emitting galaxy sample (Section 2.2), and the photoionization modelling of their NIRCam spectral energy distributions (SEDs; Section 2.3). In Section 3, we then discuss the Ly α properties inferred for these Ly α emitting galaxies, comparing them to literature samples spanning a wide range of redshifts. Based on the photometric samples identified from the NIRCam imaging, we characterize the large-scale environments of the newly identified Ly α emitting galaxies at z > 7 in Section 4. In Section 5, we discuss the dependence on galaxy local environment for Ly α visibility, considering the detection of a large number of z ∼ 7–8 Ly α emitting galaxies in the EGS field and the new JWST observations of Ly α emission in the z = 7.88 protocluster presented in Morishita et al. (2023). Finally, we summarize our findings in Section 6. Throughout this paper, we adopt a flat Lambda cold dark matter cosmology with H₀ = 70 km s⁻¹ Mpc⁻¹, Ω_m = 0.3, and |$\Omega _\Lambda = 0.7$|⁠. All magnitudes are measured in the AB system (Oke & Gunn 1983), and the emission-line EWs are calculated in the rest frame.

2 DATA AND SAMPLE PROPERTIES

The goal of this paper is to characterize the Ly α emission strengths in galaxies observed with JWST/NIRSpec (Jakobsen et al. 2022; Böker et al. 2023). We primarily focus on the CEERS program in the EGS field (Finkelstein et al. in preparation, also see Finkelstein et al. 2022, 2023). We first describe the NIRSpec spectra and the spectroscopic sample in Section 2.1. We identify the z ≥ 5 Ly α emitting galaxies and constrain their Ly α properties in Section 2.2. Based on broad-band SEDs of the Ly α emitting galaxies, we then derive their galaxy physical properties through photoionization modelling in Section 2.3.

2.1 Spectroscopy of z ≳ 5 galaxies with JWST/NIRSpec

Our work builds on the previous efforts in characterizing the Ly α emission at z ≳ 7 from the CEERS NIRSpec spectra (e.g. Jung et al. 2023; Tang et al. 2023; Witten et al. 2023). To date, six Ly α emitting galaxies have been investigated at z ≳ 7 using the spectra taken in the CEERS epoch 2 (2022 December) observations. Two of these are newly identified Ly α detections, whereas the other four were previously reported in the literature based on ground-based spectroscopy. Here, we extend the Ly α emitting galaxy search to slightly lower redshifts (5 ≤ z < 7), where the IGM is expected to be significantly more ionized. This requires us to focus on the prism observations, where the wavelength coverage extends blue enough to recover Ly α at z ≥ 5. To obtain a complete sample of z ≥ 5 Ly α emitters in CEERS, we augment this with an updated search for z > 7 Ly α emitting galaxies from the new NIRSpec prism spectra taken since the first Ly α analysis in CEERS by Tang et al. (2023). This includes the CEERS epoch 3 observations (taken in 2023 February; see Arrabal Haro et al. in preparation) and the Director’s Discretionary Time (DDT) observations (DDT-2750, PI: P. Arrabal Haro; Arrabal Haro et al. 2023a, b), which focus on NIRCam-selected targets and allows us to detect fainter z ≥ 5 galaxies. In Section 5, we characterize the Ly α properties within a spectroscopically confirmed overdensity at z ∼ 7.88. Here, we utilize the NIRSpec prism observations for the overdensity recently taken as part of the Ultradeep NIRSpec and NIRCam ObserVations before the Epoch of Reionization (UNCOVER;² GO-2561, PIs: I. Labbé & R. Bezanson; Bezanson et al. 2022; Goulding et al. 2023; Wang et al. 2023; S. Price et al., in preparation) program. Next, we will briefly describe the spectroscopic observations and identification of z ≳ 5 galaxies in the CEERS field, but we note that the UNCOVER spectra are reduced and analysed in the same manner as the CEERS data.

The CEERS NIRSpec multi-object spectroscopy observations employ the Micro Shutter Array (MSA; Ferruit et al. 2022), with a detailed description of the configuration given in the CEERS phase-2 public PDF³ and Arrabal Haro et al. (in preparation). Briefly, the epoch 2 (2022 December) observations consist of six pointings with medium-resolution grating using the G140M/F100LP, G235M/F170LP, and G395M/F290LP grating/filter pairs and eight pointings with the low-resolution prism. The prism observations that are the focus of this paper deliver simultaneous spectral coverage spanning 0.6–5.4 μm, with a spectral resolution of R ∼ 100 (σ ≃ 1300 km s⁻¹). The observed targets on the prism pointings include 466 galaxies, spanning a wide (photometric or spectroscopic) redshift range from z = 0.1–12 (Arrabal Haro et al. 2023a, b). Each of the targets was observed using the 3-shutter MSA slitlet (1.5 arcsec × 0.2 arcsec), with a total exposure time of 3107 s split into three exposures of 14 groups per disperser and per pointing. We also take advantage of the two new prism pointings taken in CEERS epoch 3 (2023 February) observations and the additional prism pointing taken in the DDT observations (DD-2750; Arrabal Haro et al. 2023b). The two CEERS epoch 3 pointings target in total 283 NIRCam-selected galaxies also across a wide range of redshift, with a total integration time of 6127 and 2042 s, respectively. The DDT observation employs a single 5.1 h pointing, mainly targeting z ∼ 12–16 galaxy candidates selected with NIRCam but also with galaxies at lower redshifts in the MSA (145 galaxies in total). All the NIRSpec spectra are reduced in the same manner as presented in Tang et al. (2023) using the JWST data reduction pipeline.⁴ In this step, we have applied slit loss corrections assuming a point source given that the majority of sources analysed in this work are not significantly extended. More information on the data reduction is provided in Tang et al. (2023).

Each extracted 2D spectrum was visually inspected by three co-authors (M. Tang, M. Topping, and ZC) in search of robust emission-line detections indicating a redshift of z ≥ 5. Our redshift identification process is similar to that described in Tang et al. (2023), requiring multiple emission-line detections (often the [O iii] doublet and H β or H α). After investigating the full data base, we identify 69 emission-line galaxies at 5 ≤ z < 7.⁵ We further identify 10 new galaxies at 7 ≤ z ≤ 9 from the more recent CEERS observations conducted since Tang et al. (2023). In the following sections, we will add these systems to the 21 z ≳ 7 galaxies in Tang et al. (2023), resulting in a total sample of 100 galaxies at 5 ≤ z ≤ 9. A significant subset of these sources overlap with those presented in the literature, and we have verified that our redshifts are in excellent agreement with those published elsewhere (e.g. Arrabal Haro et al. 2023a; Fujimoto et al. 2023; Harikane et al. 2024; Nakajima et al. 2023; Sanders et al. 2023). The full prism spectroscopic catalogue will be presented in a future work (Chen et al., in preparation).

The 1D spectra are extracted with a boxcar aperture following the procedures in Tang et al. (2023), which is set to match the emission-line profile along the spatial direction for each target (typical aperture is 6 pixels). We determine the systemic redshifts by fitting Gaussian profiles (plus a constant factor for the underlying continuum) to the brightest rest-optical lines (i.e. [O iii] + H β or H α). The derived redshifts for the 69 prism-selected galaxies at 5 ≤ z < 7 range from z = 4.98 to 6.98, with a median of z = 5.73. For the 10 z > 7 galaxies identified with the new data sets taken in 2023 February–March, their redshifts span from 7.00 to 8.75, with a median of 7.48.

We measure optical emission-line fluxes through Gaussian profile fitting, following the method of Sanders et al. (2023). We allow the centroid of each line to vary from the determined redshift within half of the spectral velocity resolution (∼650 km s⁻¹ for the prism spectra). To account for the wavelength dependence of the spectral resolution, we also restrict the line width to be within 50 per cent to 150 per cent of the line widths determined from the strongest optical lines (i.e. [O iii] + H β or H α). The uncertainties are derived by perturbing the fluxes according to the error spectra and repeating the measurements 500 times. Dust corrections will be conducted via the Balmer decrement where available (and the SED where the Balmer decrement is not measured). We will discuss the corrections in Section 2.3 once SEDs have been characterized for our sample.

We measure HST/ACS + JWST/NIRCam photometry for our spectroscopic sample of z > 5 galaxies, following the photometric procedures detailed in Endsley et al. (2023b) (also see Tang et al. 2023; Whitler et al. 2023a). We utilize the NIRCam imaging in the EGS field observed as part of the CEERS program in two scheduling windows in 2022 June and December (see Bagley et al. 2023), and the HST/ACS (WFC3/IR as well for those without NIRCam coverage) imaging assembled and reduced with grizli (Brammer et al. 2022) as part of the Complete Hubble Archive for Galaxy Evolution (CHArGE) project (Kokorev et al. 2022; Kokorev et al. in preparation). Among the 79 new sources considered in this paper, 59 are covered by CEERS NIRCam observations. We measure their Kron photometry (Kron 1980) in ACS F435W, F606W, and F814W filters, and in NIRCam F115W, F150W, F200W, F275W, F356W, F410M, and F444W filters. For the remaining 20 sources, the flux density in the HST WFC3/IR filters (F125W, F140W, and F160W) are measured instead. Fig. 1 shows the NIRCam/F150W (or WFC3/F160W when not covered with NIRCam) magnitude distribution for the 79 galaxies we present in this paper, ranging from 24.8 to 28.9 with a median of 27.1 mag.

2.2 Ly α detections at z ≳ 5 in CEERS spectra

We now seek to identify and characterize Ly α emission from the CEERS spectroscopic sample established in the previous section, building on our earlier work characterizing z ≳ 7 Ly α emitting galaxies in CEERS (Tang et al. 2023). We first will focus on Ly α emitting galaxies at 5 ≤ z < 7 in the prism spectra before discussing new systems at z > 7 from data obtained since Tang et al. (2023).

We visually search for Ly α in each of the 69 galaxies with prism redshifts at 5 ≤ z ≤ 7. We detect emission features at the wavelength of Ly α in 10 galaxies (see Table 1), with redshifts spanning z = 4.999–6.749 (median 5.67) and SNR = 5.6–26.9 (median 11.9). We have verified that these emission features are not due to cosmic rays by visually inspecting the individual exposures, thus most likely corresponding to Ly α detection. For one source (CEERS-80925), the Ly α SNR is below 7 (5.6), and we consider it to be a tentative detection. We will include this in our analysis below, but our results are not significantly altered if we remove this source. We show the prism spectra in Fig. 2, highlighting the rest-optical lines and Ly α detection. The UV continuum is also detected for the 10 Ly α emitters, at SNR = 5.8–132 (median 13.6) per rest-frame 150 Å bin. The continuum (measured in consistent aperture as the Ly α flux) allows direct constraints to be placed on the Ly α EW. We first fit the UV continuum with a linear function over the window of [100, 400] Å redwards of Ly α in the rest frame. We then compute the total Ly α flux by directly integrating the continuum-subtracted line profile over the window of [−50, 80] Å (in the rest frame) around the line centre. This aperture is chosen to match the observed width of the line. The measured Ly α fluxes range from 1.59 × 10⁻¹⁸ to 3.69 × 10⁻¹⁷ erg s⁻¹ cm⁻², corresponding to luminosities of 5.56 × 10⁴¹ to 1.17 × 10⁴³ erg s⁻¹. The resulting Ly α EWs are presented in Table 2, ranging from 19 to 286 Å with a median EW = 134 Å.

For the 10 galaxies with Ly α at 5 ≤ z < 7, their NIRCam H-band (F150W) magnitudes (or ACS F160W if not covered by NIRCam) range from 26.1 to 28.7 (median 27.3, see Table 1), which are on average fainter than the full CEERS spectroscopic sample (Fig. 1). Following Tang et al. (2023), we estimate their absolute UV magnitudes from the NIRCam (or WFC3) photometry in the filter where the central wavelength is closest to the rest-frame 1500 Å. The resulting absolute UV magnitudes span M_UV =−18.3 to −20.2 (median −19.5). These values overlap with a subset of the CEERS z > 7 Ly α emitting galaxies recently identified in Tang et al. (2023), but are fainter than the z > 7 Ly α emitting galaxies identified with ground-based observations in this field (M_UV <−20.0; e.g. Oesch et al. 2015; Zitrin et al. 2015; Stark et al. 2017; Jung et al. 2022; Larson et al. 2022). We note that given the flux limit of the CEERS prism spectra, only large EW Ly α (> 100 Å) will be detected in galaxies with faint continuum (M_UV >−20.0; see Fig. 3). The sources with high Ly α EW identified in this work are similar to the strong Ly α emitting galaxies found at lower redshifts (e.g. Malhotra & Rhoads 2002; Cowie, Barger & Hu 2011; Nakajima et al. 2012; Yang et al. 2017; Izotov et al. 2018; Tang et al. 2021; Flury et al. 2022a; Naidu et al. 2022). For context, the CEERS Ly α emitting galaxies are up to ≈3 mag fainter in M_UV than GN-z11 (M_UV = −21.5 and Ly α EW = 18 Å; Oesch et al. 2015; Bunker et al. 2023), powering up to 10 times larger Ly α EW. Their large Ly α EWs and low UV luminosities approach what is measured for the ultra-faint Ly α emitting galaxy at z = 7.3, JADES-GS-z7-LA (M_UV = −16.7 and Ly α EW = 400 Å; Saxena et al. 2023b).

We now investigate the CEERS spectra obtained after the first epoch of observations with the goal of identifying Ly α emitters at z ≳ 7. The data reveal two likely Ly α detections, CEERS-80374 at z = 7.17, and CEERS-80239 at z = 7.48. Both Ly α emitting galaxies are identified in the CEERS epoch 3 observations taken in 2023 February. The spectra show multiple rest-optical emission lines (i.e. H β and [O iii]) as well as an emission-line feature at the position of Ly α at SNR = 7.8 and 8.2 that supports our identification (see Fig. 2). Both Ly α emitting galaxies show the detection of UV continuum in the spectra. We compute Ly α EWs in the same manner as above. The derived values are |$205_{-27}^{+48}$| Å for CEERS-80374 and |$334_{-62}^{+109}$| Å for CEERS-80239, among the highest EWs detected at z > 7 (e.g. Jung et al. 2020; Bunker et al. 2023; Saxena et al. 2023b; Tang et al. 2023). We measure M_UV = −18.5 for CEERS-80374 and −18.2 for CEERS-80239 as above, which are among the faintest UV continuum magnitudes in our sample (see Fig. 3). The UV magnitudes are also fainter than the majority of the z > 7 Ly α emitting galaxies found previously, approaching the values measured for the strongest Ly α emitting galaxies found in other surveys (M_UV ∼ −15 to −17; e.g. Maseda et al. 2018, 2023; Saxena et al. 2023b). We will come back to discuss this population in later sections.

To assess the robustness of our Ly α EW measurements from the prism spectra, we simulate the observed Ly α spectra for galaxies spanning the redshift and M_UV range relevant to our Ly α emitting galaxy sample. We create mock 1D Ly α spectra for galaxies spanning z = 5–10 in steps of Δz = 0.2 and M_UV = −16 to −22 in steps of ΔM_UV = 0.5. Each Ly α spectrum consists of a skewed Gaussian profile as the intrinsic Ly α emission line before IGM attenuation. We also adopt a flat continuum that matches the assumed M_UV, where the flat continuum is consistent with the median UV slope (β = −2.3, see Section 2.3) measured in our sample of Ly α emitting galaxies. For the Ly α line, we assume a narrow [full width at half-maximum (FWHM) = 300 km s⁻¹] profile with the peak redshifted by 200 km s⁻¹ relative to the systemic redshift, consistent with observations of Green Pea galaxies in the local universe where the IGM effect is negligible (e.g. Orlitová et al. 2018). However, we note that the results do not strongly depend on the assumed FWHM and velocity offset given the much poorer prism spectral resolution (∼ 6000–10 000 km s⁻¹ at the observed Ly α wavelength). By adjusting the amplitude of the line, we obtain the intrinsic Ly α profile. We choose amplitudes that correspond to EW = 10–400 Å in steps of 5 Å. In this process, we have assumed that emission bluewards of 1216 Å is completely attenuated by IGM. Next, we convolve the intrinsic Ly α profiles to match the wavelength-dependent resolution of the prism observations and rebin the spectra to the observed wavelength grid.⁶ We also add wavelength-dependent noise that is estimated from the median noise of our spectroscopic sample at the same observed wavelength grid. An example of the resulting mock spectra is shown in Fig. 4. For each simulated observation of Ly α spectrum, we then compute the observed Ly α EW using the same method as we described above, adopting the same integration window of [−50, 80] Å relative to the line centre in the rest-frame to compute the line flux. We repeat this process 100 times to obtain the median observed Ly α EW and SNR at each redshift, M_UV, and input Ly α EW. Considering only the spectra where the Ly α is observed at SNR > 7 (similar to our adopted threshold for identifying Ly α in CEERS observations), we measure the observed Ly α EW from simulated spectra and compare them with the input values (Fig. 5).

At z = 6, for mock spectra with input Ly α EW >90 Å and M_UV <−18.5, we recover similarly high EW values, which in general are only slightly smaller (by ≲20 per cent) than the input EWs. The smaller ‘observed’ EWs are expected if a fraction of the Ly α photons fall on pixels dominated by the continuum break at the prism resolution, reducing the recovered line flux (Fig. 4; also see Jones et al. 2023). The percentage difference between the measured and input Ly α EW increases towards smaller input EW, with the measured value potentially underestimated by 30–50 per cent for more moderate EW (≲70 Å) lines. Among the eight strongest Ly α emitting galaxies (>80 Å) in our sample, we find that the EWs may be underestimated by a small fraction (≲25 per cent) due to the instrumental effects described above. For the remaining four sources with weaker Ly α (23–80 Å), the left panel of Fig. 5 suggests the measured values may be underestimated by up to 50 per cent. We note that the left panel of Fig. 4 indicates that the EW could also be slightly overestimated (by ≲15 per cent) at extremely faint magnitudes (M_UV=−17.0). This is due to fluctuations of EW measurements with the low SNR continuum at such faint magnitudes, and all galaxies in our sample are brighter than this magnitude and thus have their Ly α EWs underestimated. In what follows, we will report the measured values, but we will also consider how our results may change if the Ly α EWs are underestimated.

As the majority of galaxies in our spectroscopic sample do not show Ly α, we now consider what EW limits we can place on them given the prism non-detections. The right panel of Fig. 5 shows the limiting input Ly α EW that leads to a > 7σ emission line in the mock spectra. For bright (M_UV <−20) sources with intrinsic EW = 40 Å, we are able to detect Ly α at z ∼ 5–6. This translates to an observed EW ≈ 20 Å that corresponds to the smallest Ly α EW value measured in our CEERS Ly α emitting galaxies. We also see in Fig. 5 that the Ly α EW limit increases rapidly at lower luminosities. At z = 5–8, these tests suggest that a Ly α non-detection at M_UV = −18 would correspond to EW < 90–180 Å, whereas at M_UV = −20, the non-detection would correspond to EW < 30–50 Å. We note that these tests are consistent with Ly α being detected in GN-z11 at moderate resolution (R = 1000) but not with the prism (Bunker et al. 2023): given its absolute magnitude (M_UV = −21.5) and redshift (z = 10.6), we find that the R = 100 prism can only recover Ly α with EW>35 Å (assuming continuum per pixel SNR = 20 and no damping wing effect of IGM absorption), above the value recovered at higher resolution (EW = 18 ± 2 Å).

For individual galaxies with non-detections of Ly α in our sample, we have derived the Ly α EW upper limits by direct integration of the error spectra. We use the same integration window (i.e. [−50, 80] Å relative to Ly α in the rest frame) used for individual Ly α emitting galaxies in our sample. The resulting 7σ upper limits on Ly α fluxes range from 4.2 × 10⁻¹⁸ to 6.0 × 10⁻¹⁷ erg s⁻¹ cm⁻², translating into EW spanning from 15 to 560 Å. However, as already shown above with the mock Ly α observations, when the Ly α is detected, the observed EW could be underestimated (by up to 50 per cent) compared to the input Ly α EW. As a result, the Ly α upper limits derived from the observed spectra could still be underestimated given the poor resolution of the prism observations. In our following analysis, we will include these non-detections of Ly α and consider the impact of underestimated upper limits on our results. By combining these non-detections with the 12 Ly α emitting galaxies, we will be able to assess how common strong Ly α emission lines are at different redshift ranges, which we will discuss in the following sections.

We note again that the Ly α measurements reported in this section may be missing diffuse flux scattered from gas surrounding the UV-bright components. Studies with the VLT/MUSE Integral Field Spectrograph have characterized the surface brightness profile of Ly α haloes for large numbers of faint Ly α emitting galaxies up to z ∼ 6 (e.g. Wisotzki et al. 2016; Leclercq et al. 2017). To assess roughly how much of the total Ly α flux we may be missing in these NIRSpec microshutter observations, we compare the three compact z ∼ 6 Ly α emitting galaxies that have been observed by both NIRSpec prism and VLT/MUSE IFS (Saxena et al. 2023a). This analysis will be presented in Tang et al. (in preparation). The NIRSpec line flux is taken from Saxena et al. (2023a), and the MUSE Ly α flux is from Bacon et al. (2017, 2023) extracted with the ‘origin’ software (Mary et al. 2020), and more details will also be described in Tang et al. (in preparation). For the three systems, we find that the NIRSpec measured Ly α fluxes are 70–80 per cent of the fluxes measured from VLT/MUSE. If we assume the z ∼ 5 Ly α surface brightness profile from Leclercq et al. (2017), the MUSE ‘ORIGIN’ aperture in turn recovers ∼75 per cent of the total Ly α after accounting for the diffuse flux in the Ly α halo. This suggests that, for the z ≤ 5 galaxies in our sample, the total Ly α flux could be up to a factor of 1.67 larger than our NIRSpec measurements. In this paper, we are mostly focused on a comparison of Ly α properties in galaxies observed with NIRSpec MSA, somewhat mitigating the impact of the slit loss associated with Ly α haloes. Future studies with NIRSpec Integral Field Spectroscopy are required to assess slit loss fractions in larger samples of Ly α emitters in the reionization era.

2.3 SED fitting and photoionization modelling

The stellar population properties of our spectroscopic sample, including the new z ≥ 5 Ly α emitting galaxies, are derived from SEDs. The CEERS NIRCam footprint covers 10 of the 12 Ly α emitters (49 of the 67 galaxies not detected in Ly α) in our new prism sample. We will focus primarily on these sources in this section given the significant improvement NIRCam provides relative to HST and Spitzer. We show the NIRCam SEDs of the 10 newly identified Ly α emitting galaxies in Fig. 6. Their SEDs show characteristic flux excesses (0.3–1.6 mag) from [O iii] + H β and H α in the NIRCam long-wavelength filters that indicate the presence of young massive stars. Through power-law fitting to the broad-band SED sampling the rest-frame UV (see Topping et al. 2022b, 2023; Cullen et al. 2023), we estimate that their UV slopes (β) range from −1.2 to −2.7, with the median −2.3 (see Table 3). The bluest UV slopes are also consistent with substantial leakage of ionizing photons (e.g. Bouwens et al. 2010; Ono et al. 2012; Topping et al. 2022b, 2023; Furtak et al. 2023), which we will discuss next. We infer the physical properties of the galaxies in our sample by fitting the NIRCam photometry with the BayEsian Analysis of GaLaxy sEds (Beagle; Chevallard & Charlot 2016) code. We largely follow the modelling process described in Endsley et al. (2023b) and Tang et al. (2023), which we summarize next. For the 10 Ly α emitting galaxies, we also consider a set of models allowing the escape of ionizing photons from the host galaxies, which will also be described next. In both cases, we employ the Gutkin, Charlot & Bruzual (2016) photoionization models that self-consistently combine the stellar emission from the latest Bruzual & Charlot (2003) stellar population synthesis models with the nebular emission computed with the photoionization code cloudy (Ferland et al. 2013). We fix the redshift to the spectroscopic values determined in Section 2.1, and we only consider the NIRCam filters redwards of the Ly α emission line.

We assume a constant star formation history (CSFH), adopting the Chabrier (2003) initial mass function with the mass range of 0.1–300 M_⊙. We place a log-uniform prior on the galaxy age (defined as the duration of the constant star formation) from 1 Myr to the age of the universe at the given redshift, and on the stellar mass over the range of 5 ≤ log (M_*/M_⊙) ≤ 12. However, we note that in cases where the star formation history has undergone a recent upturn, the stellar masses inferred from CSFH could substantially underestimate the total stellar mass if there is a hidden older stellar population (e.g. Topping et al. 2022a; Endsley et al. 2023a; Tacchella et al. 2023; Whitler et al. 2023b, c). We assume log-uniform priors on both stellar metallicities and ionization parameters over a broad range of the parameter space allowed by the models (i.e. −2.2 ≤ log (Z/Z_⊙) ≤ −0.3; −4.0 ≤ log U ≤ −1.0). The interstellar (the gas-phase + dust) metallicity is kept the same as the stellar metallicities through a fixed depletion factor (ξ_d = 0.3) of metals into dust grains, and we assume the gas density n_H = 100 cm⁻³. We use the SMC curve (Pei 1992) with V-band optical depth (τ_V) adjusted from 0.001 ≤ τ_V ≤ 5. Finally, we assume the Inoue et al. (2014) model to account for the IGM attenuation.

Our default BEAGLE models are ionization-bounded, assuming that all ionizing photons emitted from stars have been reprocessed by the ISM into nebular emission lines or absorbed by dust grains. Previous work has shown that extremely strong Ly α emitting galaxies (i.e. > 50–100 Å) could leak a significant amount of Lyman continuum (LyC) photons (e.g. Dijkstra, Gronke & Venkatesan 2016; Verhamme et al. 2017; Izotov et al. 2020; Pahl et al. 2021; Flury et al. 2022b), diminishing the nebular continuum and emission-line flux. This can result in very blue UV slopes (β < −2.8) that are not seen in the standard ionization-bounded models (e.g. Raiter, Schaerer & Fosbury 2010; Zackrisson et al. 2017; Yamanaka et al. 2020; Chisholm et al. 2022; Topping et al. 2022b, 2023; Kim et al. 2023). Given the presence of very blue sources in our sample, we also consider a set of BEAGLE models that account for the escape of LyC photons (⁠|$f_{\rm esc}^{\rm LyC}$| models). The models we adopt use a ‘picket fence’ geometry where a fraction of LyC photons escape through low H i density sightlines from the host galaxy. (e.g. Heckman et al. 2001, 2011; Reddy et al. 2016; Gazagnes et al. 2020). We place a uniform prior on the LyC escape fraction (⁠|$f_{\rm esc}^{\rm LyC}$|⁠) over [0,1], with the remaining BEAGLE parameters kept the same as in our fiducial ionization-bounded models.

For the majority (nine-tenths) of the Ly α emitting galaxies, the |$f_{\rm esc}^{\rm LyC}$| models provide equally accurate fits to the NIRCam SEDs, with the implied |$f_{\rm esc}^{\rm LyC}$| in the range 0.19–0.47. In both sets of Beagle fits, the minimum χ² values are found to have small differences between the two models. Their fitted parameters from the two models are mostly consistent, except that the |$f_{\rm esc}^{\rm LyC}$| models prefer slightly younger CSFH ages (∼0.36 dex) and higher stellar masses (∼0.43 dex), similar to what is found in Tang et al. (2023). However, for CEERS-82069, where the observed SEDs cannot be fitted well with the fiducial BEAGLE models (best χ² = 36.0), the |$f_{\rm esc}^{\rm LyC}$| models improve the fitted SED significantly (best χ² = 1.71, see its SED in Fig. 6). Its potential high escape fraction of LyC photons is also consistent with the very blue UV continuum slope (⁠|$\beta = -2.47_{-0.34}^{+0.33}$|⁠) and the relatively weak [O iii] + H β (⁠|$682_{-315}^{+344}$| Å, the smallest value in our sample), both of which may be expected when a fraction of the ionizing photons escape rather than being used to produce nebular continuum and emission lines by the ISM (e.g. Topping et al. 2022b, 2023). We also obtain an O32 (= [O iii]λλ4959,5007/[O ii]λ3728) lower limit of >2.6 (3σ) from the weakly detected [O iii] and non-detection of [O ii] (Fig. 2), but deeper spectroscopy is required for a more robust measure of O32. Furthermore, as we will show in Section 3, we also infer a high Ly α escape fraction (⁠|$f_{\rm esc}^{ {\rm Ly}\alpha }$| = |$0.787_{-0.048}^{+0.053}$|⁠) that is also potentially consistent with a high escape fraction of LyC photons (e.g. Dijkstra, Gronke & Venkatesan 2016; Flury et al. 2022b). Therefore, the SED fitting results derived from the |$f_{\rm esc}^{\rm LyC}$| models are preferred over those from the default ionization-bounded models for CEERS-82069. In the following analysis, we will primarily focus on the results from the |$f_{\rm esc}^{\rm LyC}$| models for CEERS-82069, and the results from our default models for the remaining nine galaxies, which are reported in Table 3.

The Beagle models of the Ly α emitters show large flux excesses in the NIRCam SEDs at 3–5 μm, resulting in the models preferring large EW [O iii] + H β (682 Å up to 4300 Å in the rest frame). The median EW [O iii] + H β is 1982 Å, significantly higher than the typical values inferred for z ∼ 6–8 galaxies (median ≈ 780 Å; Endsley et al. 2021a, 2023a, b). The very high [O iii] + H β EWs also translate into extremely young stellar population ages (assuming CSFH) ranging from 2 to 54 Myr (median 6 Myr), younger than the more general populations at these redshifts (median 30–69 Myr; Endsley et al. 2023a; Whitler et al. 2023c) but comparable to the z > 7 Ly α emitting galaxies found in Tang et al. (2023). The stellar masses associated with these young stellar populations are estimated to be 1.3 × 10⁷–1.4 × 10⁸M_⊙ (median 2.3 × 10⁷M_⊙). The corresponding specific star formation rates range between 19 and 473 Gyr⁻¹ and occupy the high end of the sSFR distribution of UV-bright galaxies at this redshift (median |$18_{-5}^{+7}$| Gyr⁻¹; Topping et al. 2022a), consistent with rapidly rising SFHs occurring in these Ly α emitting galaxies (Endsley et al. 2023a, b). The rest-optical spectra also reveal O32 (= [O iii]λλ4959,5007/[O ii]λ3728) ratios up to >58.2 (median >7.3), suggesting gas under extreme ionization conditions. We also note that when the rest-optical continuum is also well detected (SNR>7 over a 1200 Å window in the rest-frame) in the prism spectra, we derive values of EW [O iii] + H β that are consistent with the Beagle-based values within the 1–2σ uncertainties.

With the constraints on the dust attenuation from SED modelling, we now reassess the dust correction to the optical emission-line fluxes measured from prism spectra. We perform dust correction considering both the Balmer decrement derived directly from the spectra and the effective optical depth estimated with Beagle. Among the 12 Ly α emitting galaxies, we detect both H α and H β at SNR(H β) = 16–106 for six systems, allowing us to estimate the dust attenuation using the H α/H β flux ratios. In one case (CEERS-1334), we use the Hγ/H β ratio (Hγ detected at SNR = 45) owing to an anomalously high H α flux. Assuming the case B recombination with the gas temperature T_e = 10⁴ K and density n_e = 250 cm⁻³, we expect an intrinsic H α/H β (Hγ/H β) = 2.86 (0.47; Draine 2011). Comparing this with the observed ratios (median H α/H β = 3.20, and Hγ/H β consistent with zero dust for CEERS-1334), we derive the corresponding optical depth in the V band (τ_V), assuming the SMC attenuation curve (Pei 1992). We find that inferred τ_V range from 0.0 (CEERS-1334 and CEERS-80916) to 0.56 (CEERS-1561), corresponding to small correction factors to the H α flux of 1.00–1.55 (median 1.24) to the six galaxies.

However, for the remaining six galaxies, the Balmer emission lines (H α and H β at z < 6.8, and H β and Hγ at z > 6.8) are poorly detected (SNR<7 for at least one line), leading to large uncertainties when deriving the dust attenuation based on the Balmer decrement. In order to estimate their dust attenuation, we rely on the optical depths in the V band derived with NIRCam SED modelling above, which span τ_V = 0.005–0.20 (median 0.009). To investigate whether these SED-based attenuation factors are applicable to emission lines, we compare the τ_V derived from the continuum SED and that from the Balmer decrement for four galaxies where both measurements are available. We find that the continuum SED suggests a τ_V that is 3.1–4.9σ smaller than the value inferred from Balmer emission lines. Since the SED-based attenuation is mostly driven by the reddening of the continuum, it is often dominated by the fraction of dust in the diffuse ISM, affecting stars that have freed themselves from their birth clouds. If the age of the stellar population is very young throughout the galaxy, we may expect the stars and nebular-emitting gas to face similar attenuation. It is unclear whether the offset of τ_V we find is due to H ii regions facing greater attenuation than stellar continuum (e.g. Calzetti et al. 2000; Förster Schreiber et al. 2009; Wuyts et al. 2011; Price et al. 2014; Reddy et al. 2015, 2020; Shivaei et al. 2020), or due to uncertainty in estimating attenuation from SEDs of very young and blue galaxies. Based on these results, when computing attenuation facing emission lines, we will adopt the Reddy et al. (2020) relation, where the nebular attenuation is 4.331× (assuming the SMC curve) of that faced by the stellar continuum. However, we note that our analysis in the following is not significantly impacted if we adopt the original τ_V inferred with Beagle fits to NIRCam SEDs. This requires small correction factors to the H α (H β) flux of 1.02–1.96 (median 1.03) for these six z < 7 galaxies without robust Balmer decrement measurements.

In the following section, we will quantify the escape fraction of Ly α for the Ly α emitting galaxies in our sample, using the H α or H β detections to compute the intrinsic Ly α luminosity. The Beagle models described above provide an independent estimate of the ionizing photon production rate, which also translates into an intrinsic Ly α luminosity and EW. In this process, Beagle takes the Ly α luminosity output by Cloudy (Ferland et al. 2013), which is supposed to account for the absorption of ionizing photons by dust before they ionize hydrogen (Charlot & Longhetti 2001; Gutkin, Charlot & Bruzual 2016). In Fig. 7, we show the intrinsic Ly α EW implied by Beagle as a function of both stellar population age (assuming CSFH) and stellar metallicity. Here, we have assumed the default ionization bounded H ii regions with an ionization parameter of log U = −2.5, and adopted the same Chabrier (2003) IMF and treatment to the gas-phase metallicity as described for our SED fitting. It is clear that the intrinsic Ly α EWs are largest at the youngest stellar population ages and lowest metallicity, reaching ≥200 Å at age ≲30 Myr and Z < 0.1Z_⊙. At the young stellar population ages that are common for our Ly α emitting galaxy sample (2–54 Myr), the Beagle models infer very high ionizing photon production efficiency log (ξ_ion/erg⁻¹ Hz) = 25.5–25.9, which corresponds to the hydrogen-ionizing photon production rate per unit intrinsic UV luminosity at rest-frame 1500 Å including nebular and stellar continuum (e.g. Chevallard et al. 2018; Tang et al. 2019). The efficient production of ionizing photons results in the production of very large Ly α EWs in the sample, with values ranging from 113 to 601 Å (median 256 Å). Compared with the observed Ly α EWs (20–334 Å), the larger intrinsic EW implies that only a fraction of the Ly α photons are able to escape from the host galaxies (on the scale of the NIRSpec microshutter) and through the IGM, which we will discuss in the following section.

3 REDSHIFT EVOLUTION OF Ly α ESCAPE FRACTION

JWST has recently opened a new window on Ly α emitting galaxies in the reionization era (e.g. Bunker et al. 2023; Jones et al. 2023; Jung et al. 2023; Saxena et al. 2023a, b; Tang et al. 2023). Initial investigations have demonstrated that many of the strongest known Ly α emitters at z ≳ 7 are not effective leakers of Ly α radiation, with typical Ly α escape fractions (0.073) that suggest significant attenuation relative to similar systems at z ≃ 0.3–3 (Bunker et al. 2023; Tang et al. 2023). It is not clear whether the line photons are mostly scattered by H i in the galaxy or in the surrounding IGM. More recent work has begun to uncover the first z ≳ 7 galaxies with very large Ly α escape fractions (≳0.3, Jung et al. 2023; Saxena et al. 2023a, b; Tang et al. 2023), leaving little room for attenuation from the IGM. These Ly α emitters potentially provide signposts of large ionized bubbles expected around overdensities of faint galaxies.

Future NIRSpec observations promise to build on the early work described above, delivering the distribution of Ly α escape fractions (as a function of galaxy properties) at z ≳ 7. If we are to robustly link these measurements to quantitative constraints on the IGM ionization state, we need knowledge of the distribution of Ly α escape fractions in similar galaxies just after reionization, where the impact of the damping wing from the neutral IGM is less important. Our primary goal in this section is to use our CEERS NIRSpec prism data base to begin investigating the Ly α escape fractions at the tail end of reionization. We seek to determine what fraction of galaxies have very high Ly α escape fractions at z ≃ 5–7 (similar to the galaxy reported in Saxena et al. 2023b) and to quantify how this fraction evolves with redshift into the reionization era.

Our CEERS prism data base consists of 69 galaxies with confident rest-optical emission-line redshifts at 5 ≤ z < 7. As reported in Section 2, 10 of these systems have Ly α detections in the prism spectra (Fig. 2). The implied EWs are often quite large, with seven in excess of 70 Å and a further three in the range 100–300 Å, comparable to the most extreme systems uncovered in ground-based surveys (e.g. Trainor et al. 2015; Hashimoto et al. 2017a, b; Maseda et al. 2018; Ning et al. 2020; Vanzella et al. 2020; Matthee et al. 2021; Kerutt et al. 2022; Torralba-Torregrosa et al. 2023). These values approach the maximum intrinsic Ly α EWs expected from stellar population synthesis models for very young and metal-poor stellar populations (200 to >600 Å; Fig. 7), suggesting a significant fraction of the line is likely to have been transmitted through both the galaxy and IGM. The majority of strong Ly α emitters have intense [O iii] + H β emission (Fig. 8), with values largely following the trend between the two quantities seen in lower redshift samples (Du et al. 2020; Tang et al. 2021). Here, we see that the typical Ly α EWs start to increase above [O iii] + H β EW = 1000 Å, with the majority having Ly α EWs between 50 and 300 Å at all redshifts considered.

To infer the Ly α escape fractions in the z ≃ 5–7 galaxies, we first quantify the intrinsic Ly α luminosity passing through the NIRSpec microshutter. We calculate the intrinsic Ly α luminosity from the dust-corrected H α flux and basic recombination assumptions. In CEERS-1334, we use the dust-corrected H β flux instead owing to its anomalously high H α flux. We also use the dust-corrected H β flux for the two z > 7 galaxies in our sample, because the H α emission line falls outside the prism wavelength coverage. We note that this definition of the Ly α escape fraction will focus on the transmission through the microshutter, as is common when Ly α escape fractions are calculated with slit-based spectrographs (e.g. Begley et al. 2024; Bunker et al. 2023; Roy et al. 2023; Saxena et al. 2023a; Tang et al. 2023). In this paper, we will primarily focus on the evolution of Ly α escape fractions derived self-consistently from NIRSpec MSA observations at z ≳ 5. While potential corrections for Ly α scattered into a diffuse halo may be somewhat more important when comparing to ground-based observations at lower redshifts (given the slightly larger slits), such comparisons are not central to this paper.

We first consider case B recombination where the intrinsic line ratio Ly α/H α (Ly α/H β) = 8.6 (24.6) at gas temperature T_e = 10⁴ K and electron density n_e = 250 cm⁻³ using the Pyneb (Luridiana, Morisset & Shaw 2015) package with the Storey & Hummer (1995) atomic data. Alternative assumptions on gas properties (T_e = 5000–20 000 K and n_e = 100–1000 cm⁻³) only change the Ly α/H α (Ly α/H β) ratios by ∼10 per cent, with larger ratios at higher T_e and larger n_e. However, we note that the assumption of case B recombination with optically thick H ii regions is not always valid, in particular for galaxies with substantial leakage of Ly α or LyC photons (e.g. Gazagnes et al. 2020). Given the large Ly α EWs quoted above, this may be relevant for a subset of our sample. Here, the escape fraction may be better estimated by case A recombination (Osterbrock & Ferland 2006), where the intrinsic Ly α/H α (Ly α/H β) ratio is 12.0 (33.7) for T_e = 10⁴ K and n_e = 250 cm⁻³ gas. This is a factor of 1.4 larger than that of case B, resulting in smaller Ly α escape fractions for a given Ly α and H α flux when case A is adopted. Unless otherwise stated, we will focus on the case B escape fractions in the following. For galaxies lacking Ly α detections, we place 7σ upper limits on |$f_{\rm esc}^{ {\rm Ly}\alpha }$|⁠, combining the Ly α flux upper limits derived in Section 2.2 and the dust-corrected H α or H β flux.

The derived Ly α escape fractions are reported in Table 2 for the sources with Ly α detections. As expected, the results confirm the presence of large Ly α escape fractions in galaxies with large EWs. The median value in our z ≃ 5–7 Ly α emitting galaxies is |$f_{\rm esc}^{ {\rm Ly}\alpha }$| = 0.28, well above typical values of UV-selected galaxies at lower redshifts (Hayes et al. 2011). If we instead adopt the case A value, we find |$f_{\rm esc}^{ {\rm Ly}\alpha }$| decreases by about 30 per cent (median |$f_{\rm esc}^{ {\rm Ly}\alpha }$| decreases to 0.20), but the values still indicate that significant transmission of Ly α photons is fairly common at z ≃ 5–7. One of our goals was identifying galaxies similar to the strong Ly α emitter with near unity |$f_{\rm esc}^{ {\rm Ly}\alpha }$| identified by Saxena et al. (2023b) at z ≃ 7.3. Our sample contains four Ly α emitting galaxies at 5 ≤ z < 7 with ≳ 50 per cent Ly α escape fractions. Two of these systems (CEERS-82069 and 476) have escape fractions (and EWs) that are similar to the z = 7.3 galaxy in Saxena et al. (2023b). Each of these has prominent Ly α coupled with weak H α lines in the rest-optical (Fig. 2), implying Ly α/H α ratios approaching the theoretical limit. The Ly α escape fractions quoted above do not account for the impact of the low-resolution prism on the recovered line flux (see Section 2.2). To estimate the impact of this effect on our results, we apply the redshift and M_UV-dependent flux corrections (Fig. 4) to our observed |$f_{\rm esc}^{ {\rm Ly}\alpha }$| values. We find that this typically leads to small increases (median 1.14×) to the Ly α escape fractions described above. We provide these corrected values in Table 2.

To establish a useful baseline for comparison against z ≳ 7 samples, we need to consider the distribution of escape fractions in the galaxy population near the end of reionization. Here, we primarily focus on the Ly α escape fractions derived at slightly lower redshifts (5.0 ≤ z < 6.8). The Ly α escape fractions derived from NIRSpec MSA observations depend on both galaxy properties (i.e. dust, H i covering fraction) and the neutral content of the IGM. To reliably link evolution in the escape fractions to the IGM, we must compare samples across cosmic time with matched galaxy properties. Here, we attempt to do this by considering only galaxies with intense rest-optical nebular emission ([O iii] + H β EW > 1000 Å based on Beagle SED modelling or direct measurements from the prism spectra) and lacking significant dust attenuation (τ_V < 0.7). This has the effect of isolating sources that are moderately metal-poor and dominated by very young stellar populations (e.g. Tang et al. 2019; Curti et al. 2023; Endsley et al. 2023a; Nakajima et al. 2023; Boyett et al. 2024; Sanders et al. 2024). There is one source in our sample (CEERS-82069) that appears to have extremely young stellar populations (CSFH age = |$6_{-4}^{+48}$| Myr) coupled with a much lower [O iii] + H β EW. As we discussed in Section 2.3, this source has a very blue UV slope that may suggest the emission lines may be weakened by the leakage of ionizing photons. To ensure we include all the young and metal-poor sources, we include this source in our sample for calculating the distribution of Ly α escape fractions. This leaves us 24 5.0 ≤ z < 6.8 galaxies satisfying our selection criteria, including nine Ly α emitting galaxies (median |$f_{\rm esc}^{ {\rm Ly}\alpha }$| = 0.25) and an additional 15 systems with Ly α non-detections (see Table 4).

To infer the |$f_{\rm esc}^{ {\rm Ly}\alpha }$| distribution, we adopt a Bayesian approach following the method outlined in Schenker et al. (2014) and Boyett et al. (2022). Briefly, we assume a lognormal distribution for |$f_{\rm esc}^{ {\rm Ly}\alpha }$| with parameters (μ, σ) in the following functional form:

where |$A= 2/\Bigl (1 + {\rm erf} \bigl (\frac{- \mu }{\sigma \sqrt{2}}\bigr) \Bigr)$| is the normalization factor, and x = |$f_{\rm esc}^{ {\rm Ly}\alpha }$|⁠. Here, we have adopted a |$f_{\rm esc}^{ {\rm Ly}\alpha }$| upper limit of 1. For systems with derived |$f_{\rm esc}^{ {\rm Ly}\alpha }$| (or the upper limits) >1, we will assume |$f_{\rm esc}^{ {\rm Ly}\alpha }$| (or the upper limits) = 1 in deriving the distribution. We use Bayes’ Theorem to derive the posterior probability distribution of model parameters (μ, σ) given the observed data (D)

where P(μ, σ|D) is the likelihood of observing the data given the model parameters, and P(μ, σ) are the priors. We adopt a uniform prior on the lognormal location parameter μ from −9.2–0 (corresponding to median x = 0.0001–1 in the linear space for the probability density function in equation 1), and Gaussian prior on the scale parameter σ (centred at 0.6 with standard deviation 0.3).

For a given set of model parameters (μ, σ), we compute the likelihood of observing the data considering both |$f_{\rm esc}^{ {\rm Ly}\alpha }$| measurements and those with upper/lower limits. For each system with |$f_{\rm esc}^{ {\rm Ly}\alpha }$| measurements (detected both in Ly α and Balmer emission lines), the Gaussian measurement uncertainty of |$f_{\rm esc}^{ {\rm Ly}\alpha }$| is given by

where |$\mu _{{\scriptscriptstyle {D}}_{i}}$| and |$\sigma _{{\scriptscriptstyle {D}}_{i}}$| are the measured |$f_{\rm esc}^{ {\rm Ly}\alpha }$| and uncertainty for the i-th source. The individual likelihood for each source is thus

For sources not detected in Ly α and with |$f_{\rm esc}^{ {\rm Ly}\alpha }$| upper limits, the individual likelihood is given by

where the x_7σ corresponds to the observed |$f_{\rm esc}^{ {\rm Ly}\alpha }$| 7σ upper limit. Similarly, in the case of |$f_{\rm esc}^{ {\rm Ly}\alpha }$| lower limits due to non-detection of Balmer emission lines (e.g. CEERS-80239), the individual likelihood is

where the x_7σ corresponds to the observed |$f_{\rm esc}^{ {\rm Ly}\alpha }$| 7σ lower limit. The likelihood over the entire data set is thus taken as the product of the individual likelihood of each source, including both |$f_{\rm esc}^{ {\rm Ly}\alpha }$| measurements and those with |$f_{\rm esc}^{ {\rm Ly}\alpha }$| upper/lower limits.

To derive the Ly α escape fraction distribution at 5.0 ≤ z < 6.8, we adopt the Markov chain Monte Carlo approach to sample the posteriors of the model parameters (μ, σ) with the Emcee package (Foreman-Mackey et al. 2013). We first consider the raw Ly α escape fractions derived from the prism spectra, and then we will consider the distribution changes if we make flux corrections accounting for the impact of the low resolution prism on the line recovery (Fig. 5). Considering the raw escape fractions, we compute the marginalized posterior distributions of lognormal parameters, finding |$\mu = -1.92_{-0.36}^{+0.37}$| and |$\sigma = 1.10_{-0.16}^{+0.21}$|⁠. In what follows, we define the percentage of the galaxies in this sample with moderately high escape fractions of Ly α photons (⁠|$f_{\rm esc}^{ {\rm Ly}\alpha }$| > 0.2) as X(⁠|$f_{\rm esc}^{ {\rm Ly}\alpha }$|>0.2). This analysis demonstrates that |$36_{-10}^{+11}\,\,\mathrm{per\,cent}$| of the high-EW ([O iii] + H β > 1000 Å) z ≃ 6 systems have |$f_{\rm esc}^{ {\rm Ly}\alpha }$|>0.2. If we apply the flux corrections to the prism measurements described above (see Fig. 5), we find a slightly larger fraction of the sample (⁠|$50_{-11}^{+11}\,\,\mathrm{per\,cent}$|⁠) has Ly α escape fractions in excess of 0.2 (lognormal parameters |$\mu = -1.53_{-0.30}^{+0.34}$|⁠, |$\sigma = 1.02_{-0.17}^{+0.20}$|⁠).

Lower redshift samples of similarly selected galaxies also show large escape fractions of Ly α emission (e.g. Henry et al. 2015; Izotov et al. 2016; Yang et al. 2017; Tang et al. 2021; Flury et al. 2022b). Applying a similar selection ([O iii] + H β EW > 1000 Å) to the z ∼ 0.3 Green Pea sample of Yang et al. (2017), we derive the |$f_{\rm esc}^{ {\rm Ly}\alpha }$| distribution for the resulting sample following our method above. We obtain a lognormal location parameter |$\mu = -1.50_{-0.36}^{+0.44}$| and scale parameter |$\sigma = 1.30_{-0.15}^{+0.18}$|⁠, which suggests significant leakage of Ly α photons (⁠|$f_{\rm esc}^{ {\rm Ly}\alpha }$| > 0.2) in |$46_{-9}^{+10}\,\,\mathrm{per\,cent}$| of the z ∼ 0.3 sample. In addition, at z ∼ 2, we find |$f_{\rm esc}^{ {\rm Ly}\alpha }$| lognormal distribution with |$\mu = -1.89_{-0.41}^{+0.42}$| and |$\sigma = 0.78_{-0.20}^{+0.22}$| with the Tang et al. (2021) extreme emission-line galaxies (EELGs) satisfying the same [O iii] + H β EW cut. The modelled |$f_{\rm esc}^{ {\rm Ly}\alpha }$| distribution translates into a fraction of |$34_{-16}^{+20}\,\,\mathrm{per\,cent}$| of the EELGs that show |$f_{\rm esc}^{ {\rm Ly}\alpha }$| > 0.2. This suggests that the fractions of large escape fraction systems in z ≃ 0.3–2 samples are broadly comparable to what we have found using NIRSpec at z = 5.0–6.8 (see Fig. 9).

We now seek to determine if existing constraints on the Ly α escape fraction at z ≳ 7 reveal evolution with respect to our benchmark measurement at z ≃ 5.0–6.8. Here, we build on the z > 7 spectroscopic catalogue assembled in Tang et al. (2023) (see also Jung et al. 2023; Larson et al. 2023), adding the 10 new galaxies we confirmed at z > 7 in this paper (see Section 2.2). The z > 7 data base contains 10 galaxies with [O iii] + H β EW > 1000 Å. We also include two galaxies that lack NIRCam SEDs but have very large O32 ratios (>5) that are often linked to large [O iii] + H β EWs (see Tang et al. 2023). (We note that our main z > 7 results do not change by including these sources.) This results in 12 z > 7 galaxies that satisfy our selection. We note that the |$f_{\rm esc}^{ {\rm Ly}\alpha }$| derived in Tang et al. (2023) assumed zero dust attenuation based on Balmer decrements from composite spectra. For consistency with our z ≃ 5–6 measurements, we apply small dust corrections to the z ≳ 7 values following the same method described in Section 2.3. This results in a median 14 per cent decrease of the original Ly α escape fractions, as the predicted intrinsic Ly α luminosity is modestly larger following the dust correction. The derived lognormal |$f_{\rm esc}^{ {\rm Ly}\alpha }$| distribution (⁠|$\mu = -2.96_{-0.43}^{+0.40}$| and |$\sigma = 1.06_{-0.19}^{+0.21}$|⁠) is markedly different from that at z ≃ 5–6.

The results demonstrate that large escape fractions (⁠|$f_{\rm esc}^{ {\rm Ly}\alpha }$| > 0.2) are only found in |$9_{-5}^{+8}$| per cent of the z > 7 systems in CEERS. Considering the corrections to the prism-based |$f_{\rm esc}^{ {\rm Ly}\alpha }$| values, we find this high-|$f_{\rm esc}^{ {\rm Ly}\alpha }$| fraction is |$10_{-5}^{+9}$| per cent (lognormal parameters |$\mu = -2.94_{-0.46}^{+0.40}$|⁠, |$\sigma = 1.07_{-0.18}^{+0.21}$|⁠). This fraction is lower than that at z ≃ 0.3, z ≃ 2, and z ≃ 5–6 (see Fig. 9), as would be expected given the inferred neutral fraction of the IGM at z > 7 (Mesinger et al. 2015; Zheng et al. 2017; Mason et al. 2018a, 2019; Hoag et al. 2019; Whitler et al. 2020; Bolan et al. 2022). Comparison between the z ∼ 6 and z > 7 samples admittedly faces limited statistics and uncertainties in the prism-based EWs, but the data are none the less consistent with a significant decline in |$f_{\rm esc}^{ {\rm Ly}\alpha }$| towards earlier times (Fig. 9). Furthermore, in Section 5, we will show that the CEERS EGS field may have large ionized sightlines that amplify Ly α relative to the global average at z ≃ 7. In this sense, the downturn in the Ly α escape fractions is likely to be greater in other deep fields. As sample sizes increase at z ≳ 7, it will be possible to build up the statistics while also considering sources matched in a wider range of galaxy properties.

4 LY α ASSOCIATIONS AT z > 7 IN THE EGS FIELD

We have demonstrated in Section 3 that large Ly α escape fractions and Ly α EWs are relatively common in faint z ≃ 6 galaxies, with |$f_{\rm esc}^{ {\rm Ly}\alpha }$|>0.2 presented in 50 per cent of those with [O iii] + H β EW>1000 Å. At z ≳ 7, we expect high |$f_{\rm esc}^{ {\rm Ly}\alpha }$| systems to become less common as the IGM grows more neutral, with detections mostly limited to galaxies in large ionized regions of the IGM. The discovery of Ly α emission in numerous z ≃ 7.5–9 galaxies in the EGS (Oesch et al. 2015; Zitrin et al. 2015; Roberts-Borsani et al. 2016; Stark et al. 2017; Tilvi et al. 2020; Jung et al. 2022; Larson et al. 2022; Cooper et al. 2023) with ground-based telescopes has long hinted at the possible presence of several ionized bubbles in the field. The first epoch of CEERS NIRSpec spectroscopy targeted a number of these previously known Ly α emitters, generally revealing very low escape fractions (Jung et al. 2023; Tang et al. 2023) that indicate Ly α photons face significant attenuation. This may suggest that these systems are situated in relatively small ionized bubbles (with significant IGM attenuation) or it may reflect large column densities of H i on galaxy scales. As these first galaxies targeted tend to be very luminous (and fairly massive) and often with large Ly α velocity offset (e.g. Bunker et al. 2023; Tang et al. 2023) implying significant scattering of Ly α photons, it is perhaps likely that the H i within the galaxy is largely responsible for reducing the measured escape fraction through the NIRSpec microshutter. In this case, the Ly α emitters could still trace large ionized structures where the IGM attenuation is minimal.

The more recent CEERS NIRSpec observations have targeted fainter z ≳ 7 galaxies identified in NIRCam images of the EGS field (e.g. Arrabal Haro et al. 2023b). In this paper, we present likely detections of Ly α in two of these galaxies, CEERS-80239 (z = 7.49) and CEERS-80374 (z = 7.10). Because of the faint continuum of both galaxies (M_UV = −18.2 and −18.5, respectively), both detections require extremely large EW Ly α emission (⁠|$334_{-62}^{+109}$| Å and |$205_{-27}^{+48}$| Å). These large Ly α EWs suggest these galaxies are situated in ionized regions that facilitate significant transmission of Ly α through the IGM (also see Witstok et al. 2024). In one of these systems (CEERS-80374), we are able to measure the Ly α escape fraction (⁠|$f_{\rm esc}^{ {\rm Ly}\alpha }$| = |$0.472_{-0.048}^{+0.066}$|⁠) via detection of H β, indicating that at least ∼50 per cent of the Ly α photons are transmitted through the ISM and IGM (with a larger IGM transmission implied if some of the attenuation occurs in the ISM and CGM of the galaxy).⁷ To achieve this large |$f_{\rm esc}^{ {\rm Ly}\alpha }$| may require that the nearest patch of mostly neutral IGM be located at least 1–2 physical Mpc (pMpc) away. Here, we follow the prescription detailed in Mason & Gronke (2020) to compute the Ly α transmission assuming a range of ionized bubble sizes and fully neutral IGM outside the bubble. Specifically, for a source observed at redshift z_s with photons observed at wavelength λ_obs = λ_em(1 + z_obs), the corresponding Ly α optical depth is

where n_H is the total hydrogen number density, x_HI is the hydrogen neutral fraction, and σ_α is the Ly α scattering cross-section through an ensemble of hydrogen atoms. To estimate the bubble size required for a given Ly α escape fraction, we have assumed the Ly α velocity offset is 100–200 km s⁻¹, a typical range for faint galaxies with strong Ly α (Saxena et al. 2023b).⁸Adopting slightly different assumptions will alter the values quoted above, but it is unlikely to change the fact that these two galaxies likely must sit in large ionized regions to facilitate observation of Ly α with EW in excess of 200–300 Å.⁹

The next step is identifying the galaxy population responsible for carving out the ionized regions that are facilitating such efficient Ly α escape. Given that both CEERS-80374 and CEERS-80239 are very faint continuum sources (and unlikely to create significant bubbles on their own), it is probable that they are part of overdense structures that are responsible for creating the ionized region that is enabling efficient Ly α escape. We first search for spectroscopically confirmed galaxies at the redshifts of both galaxies. A map is shown in Fig. 10. We note that the new Ly α emitter CEERS-80239 (z = 7.49) is at a similar redshift as CEERS-698 (z = 7.469), an extremely UV-luminous (M_UV = −21.86) galaxy shown previously to have Ly α emission with Keck (Roberts-Borsani et al. 2016; Stark et al. 2017). Taking the systemic redshifts of both, we find the two are separated by only 3.1 pMpc. CEERS NIRSpec observations also confirm an additional six z = 7.45–7.54 galaxies between CEERS-698 and CEERS-80239, hinting at an extended structure that could be contributing to the ionized bubble. The other new Ly α emitter we identify (CEERS-80374) is at a similar redshift as CEERS-44 (separated by 4.7 pMpc), another strong Ly α emitter (EW = 77Å) in the same field (Tang et al. 2023). CEERS-44 also shows a large Ly α escape fraction (⁠|$f_{\rm esc}^{ {\rm Ly}\alpha }$| = 0.339 ± 0.044), perhaps suggesting reduced IGM attenuation to Ly α photons at z ≃ 7.1–7.2 across the field. Another four galaxies with very similar spectroscopic redshifts and close separations (0.7–3.9 pMpc) from CEERS-80374 and CEERS-44 have been identified in the CEERS program (see Tang et al. 2023 and Fig. 10). We present these spectroscopically confirmed galaxies at each Ly α emitting galaxy redshift in Table 5. Unfortunately, their Ly α line either falls outside the detector or has upper limits on Ly α EW (<185 Å), making them also potential Ly α emitting galaxies within the same ionized regions.

The NIRCam imaging obtained with CEERS provides a more complete census of the galaxies likely to be at the same redshift as the two z ≳ 7 Ly α emitters identified in this study. In the remainder of this section, we characterize colour-selected galaxies likely to lie at z ≃ 7.0–7.6 and assess implications for overdense structures in the vicinity of CEERS-80374 and CEERS-80239. We use an updated catalogue of F814W-dropout galaxies spanning the entire footprint of CEERS NIRCam observations, expanding from the catalogue published in Endsley et al. (2023b). These dropout galaxies are selected with the same colour selection as in Endsley et al. (2023b) but include the full data base of CEERS imaging (whereas the earlier paper only included the first epoch). Photometric redshifts are inferred from Beagle SED fitting to the HST/ACS + JWST/NIRCam multiband photometry. The reader is directed Endsley et al. (2023b) for more details. In total, our selection yields 269 galaxies with F200W<28 mag at z ∼ 6.5–8 across the 92.1 arcmin² CEERS footprint with overlapping ACS imaging.

Our goal is to identify galaxies in this sample that are likely to be at similar redshifts as CEERS-80374 (z ≃ 7.17) and CEERS-80239 (z ≃ 7.49). We adopt an empirically motivated approach utilizing the F410M medium-band photometry to isolate galaxies at z ∼ 7.0–7.6 (see Endsley et al. 2023a for full details). In this redshift range, the [O iii] + H β emission lines fall in the F410M filter, leading to red colours in F356W − F410M. By selecting dropouts in our sample that are red in F356W − F410M (>0.6), we are likely to pick out those galaxies with large EW [O iii] + H β (> 400 Å). Given the [O iii] + H β EW distribution in CEERS (Endsley et al. 2023b), this corresponds to 82 per cent of the population at these redshifts. This colour cut results in a sample of 54 galaxies likely to lie at z ∼ 7.0–7.6 based on the F410M excesses (see Table 6). The Beagle SED fits described above indicate that the photometric redshifts of all 54 sources lie in the desired redshift range.

We show the distribution of photometric targets in Fig. 10. We identify three regions (region A, B, and C) in the vicinity of CEERS-80374 and CEERS-80239 with high surface densities of z ∼ 7.0–7.6 galaxies. Each of the two regions contains 9–12 photometric galaxies within a radius of 2 arcmin (0.6 projected pMpc), corresponding to a surface area of 7.5–10.9 arcmin² with overlapping NIRCam imaging. The derived surface density (1.1–1.3 galaxies / arcmin²) is ≈2× greater than the average over the full footprint. In particular, region A overlaps with CEERS-80239 and is found only 0.5 arcmin away from CEERS-80374, suggesting that one (or both) of the Ly α emitters may be associated with this overdensity. The other regions (B and C) are both found ≈ 5.5 arcmin away from CEERS-80239, which could be expected if the structures extend to larger distances over the field.

In order to quantify the strength of these photometric overdensities, we compare the observed galaxy number with that predicted from the UV luminosity function and our selection efficiency. Our method to compute the expected number of galaxies (〈N〉) broadly follows Endsley et al. (2023a). We summarize the approach below (also see Endsley et al. 2021a; Whitler et al. 2023a). As a first step, we consider the selection efficiency imposed by both the dropout selection and the F356W − F410M colour cut before estimating the predicted number of observable galaxies. We start by modelling the probability of a galaxy being selected by both criteria as a function of both redshift and magnitude using a model SED with a flat rest-UV continuum in f_ν (β = −2.0). At each redshift and magnitude, we add [O iii] + H β emission lines with the EWs randomly drawn from the Endsley et al. (2023b) distribution at z ∼ 7, assuming a fixed [O iii]λ5007/H β ratio = 6 (e.g. Tang et al. 2019). We compute the mock photometry in each ACS + NIRCam filter, taking into account the IGM transmission model of Inoue et al. (2014) and perturbing the photometry with the typical noise estimated from the observed values for our sample. The selection efficiency at each redshift (z = 6.5–8.0 in Δz = 0.02 unit spacing) and magnitude (F150W = 25.0–28.0 in 0.1 mag unit spacing) is then estimated as the fraction of 1000 realizations passing both our dropout criteria for z ∼ 6.5–8 galaxies and the F356W − F410M colour cut.

We combine the selection efficiency with the redshift-dependent UV luminosity function of Bouwens et al. (2022) to derive the strength of the photometric overdensities. To compute the cosmic mean galaxy number surface density, we convolve the UV luminosity function with the selection efficiency and integrate it over z = 7.0–7.6 and magnitudes brighter than M_UV = −19, the M_UV limit of our sample (also see Endsley et al. 2023b). Comparing the resulting mean surface density (0.4 galaxies arcmin⁻²) with what we observe in Regions A, B, and C, we found the three regions correspond to overdensities at |$3.1_{-0.9}^{+0.9}$|σ, |$3.1_{-1.0}^{+1.0}$|σ, and |$3.7_{-1.2}^{+1.1}$|σ, respectively. We note that our selection includes candidate galaxies over a relatively wide redshift range (z = 7.0–7.6). It is plausible each of the regions could be considerably more overdense if the redshifts of the candidate galaxies are more narrowly distributed around a given redshift.

The photometric data suggest that overdensities (in region A) may surround the two intense Ly α emitters identified in this study. The potential overdensities in regions B and C suggest that there may be extended structures across the field at the redshifts of one or both of the Ly α emitters. A single overdensity of this amplitude (∼3σ) over a radius of R ∼ 0.6 pMpc alone can create moderately sized bubbles (⁠|$R_{{\rm H\, {\small II}}} \simeq$| 0.8 pMpc) with nominal assumptions on source properties (see Endsley & Stark 2022; Whitler et al. 2023a). Given the presence of three such overdensities in this field, it is conceivable that collectively they create a much larger ionized region (⁠|$R_{{\rm H\, {\small II}}}\gtrsim$| 1–2 pMpc) as is likely required for explaining the observed escape fractions of the Ly α emitting galaxies. Given that the galaxies in this sample are selected on F356W − F410M colours that indicate strong [O iii] + H β emission, it should be trivial in the future to map out rest-optical emission-line redshifts (with NIRSpec or NIRCam grism) and better characterize whether the galaxies in the overdense regions are associated with the z ≃ 7.2 or z ≃ 7.5 extreme Ly α emitters.

5 DISCUSSION

The spectroscopic capabilities of JWST are ushering in a new era of investigating reionization. With rest-optical emission-line redshifts now easily accessible at z ≃ 6–9, it is becoming possible to fully map out the spatial distribution of galaxies in 3D across known deep fields, identifying regions that are overdense and those that are underdense. Deep Ly α follow-up will in turn offer insight into which of these regions have been able to carve out large (≳1 pMpc) ionized bubbles. In so doing, these observations provide a local view of reionization, taking a new step beyond the global reionization studies of the past two decades. In this section, we contrast observations of two z ≳ 7 overdensities with Ly α measurements.

The EGS field presents one of the most interesting cases of a likely ionized region, with 16 robust Ly α detections at 7 < z < 8 (Oesch et al. 2015; Roberts-Borsani et al. 2016; Stark et al. 2017; Tilvi et al. 2020; Jung et al. 2022, 2023; Tang et al. 2023; Witten et al. 2023). No other deep field has proven as rich in z ≳ 7 Ly α emitters, in spite of considerable observational investment. We described new Ly α detections in the EGS field in Section 4, building the case for three associations of Ly α emitters along the line of sight: at z ≃ 7.2, z ≃ 7.5, and z ≃ 7.7 (see Figure 11). The radial distances between these associations are 7 pMpc (z ≃ 7.2 and z ≃ 7.5) and 12 pMpc (z ≃ 7.5 and z ≃ 7.7), comparable to the diameters of large ionized bubbles expected at z ≃ 7 (Lu et al. 2024). The new discovery of Ly α with EW>200 Å in these associations confirms that there is likely to be large ionized sightlines (>1–2 pMpc) at z ≃ 7.1–7.7 in the EGS. Both photometric and spectroscopic samples are now revealing significant overdensities of galaxies in this redshift range (see Section 5).

These observations raise the possibility that the EGS contains an extended structure oriented radially along the line of sight over part of the region spanning z ≃ 7.2 to z ≃ 7.7. Such filament-like structures are expected to be common topological features in the reionization era (e.g. Elbers & van de Weygaert 2019, 2023). If viewed along the line of sight, ionized filaments may create the optimal conditions for detecting Ly α at z ≳ 7, providing large ionized sightlines for galaxies throughout the structure. In particular, galaxies at the far end of the structure can have their Ly α transmission significantly boosted. For example, we find that galaxies viewed at the far side of large (>10 pMpc) filaments will transmit 2.4–2.9× more Ly α flux than if they were situated in relatively small (R = 0.5 pMpc) ionized bubbles. Here, we follow the procedures described in Mason & Gronke (2020) to estimate the Ly α transmission at different sizes of the ionized region using equation (7), with the range of factors dictated by our assumed distribution of Ly α velocity offsets (50–200 km s⁻¹). Clearly, radial-extended ionized filaments will amplify many Ly α lines that otherwise would not be visible.

The future prospects for characterizing the ionized regions in the EGS are good. To date, the majority of existing Ly α detections have come from brighter galaxies identified in HST imaging. The NIRCam photometric samples extend significantly fainter, providing much-improved statistics of the galaxy distribution. Thus far, only a very small fraction of the NIRCam samples have been followed up with spectroscopy (fewer than 6 per cent of the galaxies reported in Section 4 have NIRSpec observations, and none of them have R ≥ 1000 spectra required for bubble size inferences). Once a larger subset of them have rest-optical spectra, it will be possible to map the 3D spatial distribution of galaxies through the field. Such a data set will reveal whether overdensities are peaked near the known Ly α emitters (as may be expected if we are observing a series of bubbles along the line of sight) or whether the entire redshift range spanning z ≃ 7.1–7.7 is uniformly overdense (as may be expected if there is a single extended ionized sightline). Deep R ≥ 1000 Ly α spectroscopy of the NIRCam-identified galaxies will enable direct tests of these scenarios. Here, the fainter magnitudes of the NIRCam data base are critical. As shown in Section 3, Ly α EWs often reach extremely large values (>100 Å) in UV-faint galaxies (M_UV > −18). When such objects are identified at z ≳ 7 with small Ly α velocity offsets (<100 km s⁻¹), it immediately points to large radial ionized sightlines (e.g. Mason & Gronke 2020; Saxena et al. 2023b). Comprehensive high-resolution Ly α spectroscopy of the NIRCam-selected galaxies in Fig. 10 will allow these ionized regions to be mapped across the EGS.

The connection between z ≳ 7 Ly α visibility and galaxy overdensities is now starting to be seen in several other fields (e.g. Hu et al. 2021; Endsley et al. 2023a; Witten et al. 2023; Witstok et al. 2024). However, one recently confirmed overdensity at z = 7.88 in the Abell 2744 field presents a conflicting picture (Hashimoto et al. 2023; Morishita et al. 2023). NIRSpec observations of this field from the GLASS (ERS-1324, PI: Treu; Treu et al. 2022), DDT-2756 (PI: W. Chen; Roberts-Borsani et al. 2023), and GO-1840 (PI: Hashimoto; Hashimoto et al. 2023) programs have confirmed redshifts for eight galaxies clustered within a small region with a projected radius of 60 physical kpc (pkpc). The density of galaxies in this small area implies a strong galaxy overdensity (⁠|$24_{-8}^{+12} \sigma$|⁠; Roberts-Borsani et al. 2023), albeit over much smaller physical scales than the CEERS overdensities discussed above. Five of the eight galaxies have Ly α constraints from NIRSpec. None show Ly α emission, with small reported 2σ Ly α EW upper limits ranging from 15.9 to 27.6 Å. The Ly α non-detections raise the question of whether this dense region has been able to ionize its surroundings.

The UNCOVER program (GO-2561, PIs: I. Labbé & R. Bezanson; Bezanson et al. 2022) has recently obtained NIRSpec prism spectra of two of the seven spectroscopically confirmed galaxies (GLASSz8-2 and YD4) and one of eight photometric candidates (ZD4) associated with the z = 7.88 overdensity (Morishita et al. 2023). Following the procedures described in Section 2.1, we obtained reduced spectra for the three galaxies, which are shown in Fig. 12. The photometric candidate ZD4 is confirmed to lie at z = 7.88, placing it as part of the same overdensity. In addition to strong rest-optical emission lines ([O iii] + H β EW = |$794_{-220}^{+310}$| Å), we identify an emission-line feature (SNR = 5.9) at a wavelength consistent with Ly α in the galaxy. We measure a relatively large Ly α EW (⁠|$68_{-9}^{+10}$| Å), suggesting the overdensity may indeed have carved out an ionized region. For the other two galaxies where we do not detect Ly α, we obtain 7σ upper limits on the EWs (EW < 85 Å for GLASSz8-2 and EW < 91 Å for YD4) that are consistent with the previous constraints (Morishita et al. 2023).

The UNCOVER spectra also hint at why Ly α is not seen in many galaxies of the Morishita et al. (2023) overdensity. In each of the three spectra, we see evidence for damped Ly α absorption (see Fig. 12) that will significantly attenuate any Ly α emission that is present. Strong H i absorption is perhaps not surprising given the close separations (within a projected radius of 60 pkpc) of many of the galaxies in the protocluster structure (see Figure 13). In contrast, the EGS field presents a galaxy overdensity on much larger scales (>1 pMpc), where the obscuring effect of DLAs is likely to be less common. To estimate the viability of the damped Ly α interpretation for the continuua of the UNCOVER spectra, we perform a simple analysis fitting a DLA transmission profile to the spectra over the window of 1100–2000 Å in the rest frame. We compute the Ly α optical depth using the cross-section approximation by Tasitsiomi (2006) and assume a continuum spectrum (prior to absorption by the damped Ly α system) obtained by extrapolating the power-law fit to the observed continuum over 1400–3000 Å in the rest frame. We perform the fit in the observed frame, convolving the model spectrum with the instrumental resolution of the prism. We fit using the H i column density, covering fraction, and peculiar velocity of the DLA, and the Ly α EW (fixing a Ly α FWHM of 200 km s⁻¹) prior to attenuation by the DLA as free parameters. As this region is a dense candidate protocluster, redshifted DLA systems may arise from foreground protocluster members with line-of-sight motion towards the observed galaxies, motivating our inclusion of peculiar velocity as a free parameter. We find degeneracies between the peculiar velocity of the absorber, the H i covering fraction, and the Ly α EW (prior to attenuation by the DLA). The parameters also are sensitive to different assumptions about the local ionization state of the IGM in the vicinity of the galaxies. A detailed characterization of these properties is outside the scope of this study, given these degeneracies and the resolution of the spectrum (∼4000 km s⁻¹ per pixel). Higher resolution spectroscopy is required for a more robust view of the DLA properties.

Nevertheless, in all cases, we find large H i column densities (N_HI ≈ 10^21.5–10^23.5 cm⁻²)¹⁰ are required to reproduce the observed spectra (see Fig. 12). The inferred column densities are high but similar to that inferred for other z > 8 galaxies with recent NIRSpec prism observations (e.g. Heintz et al. 2023). We suggest such column densities may be more commonly seen in protocluster environments at z ≳ 7 where a large number of galaxies are closely separated. Such dense neutral H i in the local surroundings could heavily attenuate the Ly α photons, contributing to the apparent lack of Ly α detection in the majority (five-sixths) of the member galaxies associated with the Morishita et al. (2023) overdensity. We finally note that the UNCOVER spectrum of ZD4 has recently been presented in Cameron et al. (2023). They suggest nebular continuum may provide an alternative explanation to damped Ly α absorption. Higher resolution data will ultimately be required to fully disentangle the relative contribution of these effects. However, the weak rest-optical emission lines in ZD4 (H β EW = |$217_{-85}^{+122}$| Å) suggest that such strong nebular emission is not likely present in this system. The ionizing photon production efficiency implied by the H β detection in ZD4 (log (ξ_ion/erg⁻¹ Hz) = |$25.5_{-0.2}^{+0.1}$|⁠) is well below that required for the nebular continuum to dominate the far-UV continuum (e.g. Topping et al. 2022b). Also, given the presence of such apparent absorption features redwards of Ly α wavelength in all three galaxies observed with UNCOVER in this z = 7.88 overdensity, we suggest that the damped Ly α absorption remains a very likely explanation for ZD4 and the other two systems.

Strong galaxy overdensities that span small physical size scales (within a projected radius of 60 pkpc) are more likely to have small ionized bubbles. In particular, for the Abell 2744 overdensity at z = 7.88, Lu et al. (2024) demonstrated that the range of bubble sizes was likely to be ∼0.5 pMpc (with a range between 0.2 and 1.8 Mpc), lower than that expected in the CEERS field (1.0 pMpc, with a range between 0.2–2.7 pMpc). The small bubble sizes potentially explain the weak Ly α emission reported in Morishita et al. (2023). The observations here suggest an additional factor regulating Ly α visibility in this region may be damped Ly α absorption. This is likely to be most impactful in attempts to characterize Ly α in large overdensities spanning small scales, where the close galaxy separations (<50–390 pkpc) lead to the common presence of foreground absorbers. Efforts to characterize ionized bubbles may be best focused on regions that are overdense over larger scales (>1 pMpc). Systematic surveys of Lyman alpha visibility as a function of galaxy overdensity will enable us to directly connect ionized regions to both their local environment and galaxy properties, paving the way for a detailed understanding of how galaxies achieved reionization.

6 SUMMARY

We investigate Ly α properties of galaxies at 5 ≤ z < 7 in JWST/NIRSpec prism spectroscopy taken as part of CEERS (Finkelstein et al. 2022, 2023) and DDT-2750 (Arrabal Haro et al. 2023a, b) observations of the EGS field. We summarize our main findings below.

• We have assembled a spectroscopic catalogue of 69 galaxies at 5.00 < z < 6.98 confirmed via rest-optical emission-line detections in JWST prism spectroscopy. Many of the galaxies are faint continuum sources (F150W = 27–29 mag) selected in NIRCam imaging, allowing characterization of Ly α output in very low-mass star-forming systems. We identify likely Ly α emission (SNR = 5.6–26.9) in 10 of the 69 galaxies. Given the faint continuum flux densities, the implied Ly α EWs extend to very large values (ranging up to 286 Å with a median of 134 Å). These are comparable to the Ly α EWs found in faint galaxies with MUSE (e.g. Hashimoto et al. 2017b; Maseda et al. 2018; Vanzella et al. 2020; Kerutt et al. 2022) and other narrow-band surveys (e.g. Trainor et al. 2015; Hashimoto et al. 2017a; Harikane et al. 2018; Ning et al. 2020; Matthee et al. 2021).

• We present likely detections of Ly α (SNR = 7.8 and 8.2) in two z > 7 galaxies, CEERS-80374 (z = 7.49) and CEERS-80239 (z = 7.10). These prism spectra were obtained in the second epoch of CEERS observations, after the first papers describing z ≳ 7 Ly α spectroscopy were released (Jung et al. 2023; Tang et al. 2023). Both galaxies are very faint in the continuum (F150W = 28.7 and 28.9 mag), with the Ly α detections implying extremely large Ly α EWs (⁠|$334_{-62}^{+109}$| and |$205_{-27}^{+48}$| Å). Such intense Ly α likely requires near-unity transmission through the IGM, possibly indicating the two galaxies are situated in large ionized regions, as recently found in the faint Ly α emitter JADES-GS-z7-LA (Saxena et al. 2023b).

• We explore the physical properties of our sample using NIRCam SEDs and NIRSpec rest-optical prism spectroscopy. The galaxies in our sample with Ly α emission tend to be dominated by extremely young stellar populations (median age = 6 Myr for constant star formation), as is common for galaxies that have undergone a recent burst of star formation. This results in very efficient ionizing photon production (median ξ_ion = 10^25.8 erg⁻¹ Hz) that translates into intense intrinsic Ly α emission (median EW = 256 Å). As expected for such young systems, we see intense [O iii] + H β emission (median EW = 1982 Å), well above the average for z ≃ 5–7 galaxies (Endsley et al. 2023a).

• We compute Ly α escape fractions for our sample using flux ratios of the Ly α and the Balmer emission lines (H α or H β) in the NIRSpec data. We find that large escape fractions are still common in the galaxy population at z ≃ 5–7, with a Case B median value of |$f_{\rm esc}^{ {\rm Ly}\alpha }$| = 0.28. This is well above what is common in the first z ≳ 7 spectroscopic samples (median |$f_{\rm esc}^{ {\rm Ly}\alpha }$| = 0.073; Tang et al. 2023). We find two galaxies in our z ≃ 5–7 sample that show Ly α/H α ratios suggesting close-to-unity (>80 per cent) escape fractions of Ly α. These efficient Ly α leakers appear similar to the extreme Ly α emitters now being found at z ≳ 7 (e.g. Saxena et al. 2023b).

• The distribution of Ly α escape fractions at z ≃ 5–7 in our sample provides a baseline to compare against emerging z > 7 samples (e.g. Jung et al. 2023; Tang et al. 2023), allowing insight into the impact of the neutral IGM on Ly α. We find that |$50_{-11}^{+11}$| per cent of the z ∼ 6 galaxies with large EW [O iii] + H β (> 1000 Å) have moderately large |$f_{\rm esc}^{ {\rm Ly}\alpha }$| (>0.2). We show that lower redshift (z ≃ 0–2) galaxies with similar properties show comparable Ly α escape fractions. If the IGM significantly attenuates z ≳ 7 Ly α emission, we expect this percentage to go down. While current statistics are limited, we find only |$10_{-5}^{+9}$| per cent of the CEERS z > 7 galaxies (from both Tang et al. 2023 and this work) satisfying the EW [O iii] + H β selection have such large escape fractions. As JWST builds larger z ≳ 7 spectroscopic samples, it will become possible to reduce the uncertainties on such measurements, improving efforts to link measurements of the Ly α escape fraction to the IGM.

• The two z ≳ 7 galaxies with likely detections of strong Ly α emission are found to be located at the same redshifts as other Ly α emitters in the CEERS EGS field, potentially hinting at extended ionized structures at both redshifts. We characterize the distribution of faint galaxies around both Ly α emitters using NIRCam imaging. We find evidence for photometric overdensities of z ≃ 7.0–7.6 galaxies (at the 3.1–3.7σ level), as would be expected if there are large ionized bubbles in the field. The detection of a large number of high-EW Ly α emitting galaxies combined with the presence of significant photometric overdensities also suggest the presence of long ionized sightlines that boost the Ly α emission in this field. Future spectroscopic observations are required to better characterize the overdensity structures and ionized region sizes at the redshifts of these Ly α emitters.

• We further investigate the dependence of z > 7 Ly α transmission on the local environment, considering the z = 7.88 galaxy overdensity with 7 previously confirmed member galaxies, z7p9OD, in the Abell 2744 field. Utilizing the new UNCOVER deep NIRSpec prism observations covering this overdensity, we confirm one more photometrically selected member to be associated with the overdensity, which is also the only source detected in Ly α (with EW = |$69_{-9}^{+10}$| Å) in this region. We also identify damped Ly α absorption in all three member galaxies with deep UNCOVER spectra, suggesting the presence of high column density H i in the local surroundings of the overdensity. Such H i gas may obscure much of the Ly α photons from the member galaxy, leading to the low fraction (one-sixth) of Ly α detection in this overdensity. JWST NIRSpec observations mapping a larger sample of z > 7 overdensities will be required to statistically investigate how Ly α visibility may be affected by overdense environments spanning different physical scales.

ACKNOWLEDGEMENTS

The authors thank the referee, Dr Takuya Hashimoto, for helpful comments that improved the quality of this paper. We thank Jacopo Chevallard for the use of the Beagle tool used for much of our SED fitting analysis, Gabe Brammer for providing the optical imaging of the EGS as part of CHArGE program, John Chisholm for providing data of LzLCS sample, and Ting-Yi Lu for helpful discussion on ionized bubbles. This work is based in part on observations made with the NASA/ESA/CSA JWST. The data were obtained from the Mikulski Archive for Space Telescopes at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5–03127 for JWST. These observations are associated with programs #ERS-1345 (CEERS), DDT-2750, and GO-2561 (UNCOVER). The authors acknowledge the CEERS, DDT-2750, and UNCOVER teams led by Steven L. Finkelstein, Pablo Arrabal Haro, and I. Labbé & R. Bezanson for developing their observing program with a zero-exclusive-access period.

DPS acknowledges support from the National Science Foundation through the grant AST-2109066. CM acknowledges support by the VILLUM FONDEN under grant 37459 and the Carlsberg Foundation under grant CF22-1322. The Cosmic Dawn Center (DAWN) is funded by the Danish National Research Foundation under grant DNRF140. LW acknowledges support from the National Science Foundation Graduate Research Fellowship under Grant No. DGE-2137419. MT acknowledges funding from the JWST Arizona/Steward Postdoc in Early galaxies and Reionization (JASPER) Scholar contract at the University of Arizona. This material is based in part upon High Performance Computing (HPC) resources supported by the University of Arizona TRIF, UITS, and Research, Innovation, and Impact (RII) and maintained by the UArizona Research Technologies department.

This work used astropy,¹¹ a community-developed core Python package and an ecosystem of tools and resources for astronomy (Astropy Collaboration 2013, 2018, 2022); beagle (Chevallard & Charlot 2016); emcee (Foreman-Mackey et al. 2013); jupyter (Kluyver et al. 2016); matplotlib (Hunter 2007); numpy (Harris et al. 2020); photutils, an astropy package for detection and photometry of astronomical sources (Bradley et al. 2022); scikit-image (van der Walt et al. 2014); scipy (Virtanen et al. 2020); and sedpy (Johnson 2021).

DATA AVAILABILITY

The JWST and HST imaging data used in this work are available through the MAST Portal (https://mast.stsci.edu/portal/Mashup/Clients/Mast/Portal.html). The data product and analysis code used in this work will be shared upon reasonable request to the corresponding author.

Footnotes

References