New z $>$ 7 Lyman-alpha Emitters in EGS: Evidence of an Extended Ionized Structure at z $\sim$ 7.7

Targets were selected from the photometrically-selected high-redshift galaxy catalog of Finkelstein et al. (2022), which is based on the updated HST CANDELS photometry. The photometric selection of high-redshift galaxies is done as described in Section 3.2 in Finkelstein et al. (2015), using the photometric redshift ($z_{p}$) probability distribution functions (PDFs) of $P(z)$ calculated by EAZY (Brammer et al., 2008). Then, we created a target list of $J_{\text{125}}\lesssim 27$ galaxies with $z_{p}>6$ in the CANDELS/EGS field, which was used for designing optimized slitmask configurations in MAGMA³³3https://www2.keck.hawaii.edu/inst/mosfire/magma.html for our Keck/MOSFIRE observations. In our slitmask design, we prioritized targets on slits based on the galaxy brightness ($J_{\text{125}}$) and the integrals of $P(z)$ in $7.0<z<8.2$, which corresponds to the MOSFIRE $Y$-band wavelength coverage for Ly$\alpha$ emission. This resulted in 61 Ly$\alpha$ targets across our four MOSFIRE pointings.

We use the photometric catalog of Finkelstein et al. (2022), which includes the HST ACS and WFC3 broadband photometry ($V_{606}$, $I_{814}$, $Y_{105}$, $J_{125}$, $JH_{140}$, and $H_{160}$) in addition to Spitzer/IRAC 3.6$\mu$m and 4.5$\mu$m band fluxes in the CANDELS/EGS field. We also use photometric redshift measurements that have been obtained with EAZY in Finkelstein et al. (2022) based on the updated CANDELS photometry.

To derive galaxy physical properties, we performed spectral energy distribution (SED) fitting with the photometric data to galaxy SED models. In the construction of galaxy model SEDs, we assume a Salpeter (1955) initial mass function with a stellar mass range of 0.1-100$M_{\odot}$. We allow a range of metallicity from 0.005$Z_{\odot}$ to 1.0$Z_{\odot}$, and exponential models of star formation histories are used with exponentially varying timescales, parameterized with $\tau=$ 10 Myr, 100 Myr, 1 Gyr, 10 Gyr, 100 Gyr, $-$300 Myr, $-$1 Gyr, $-$10 Gyr. We use the Calzetti (2001) dust attenuation description for a ranging from 0 to 0.8 mag in $E(B-V)$ values. Nebular emission lines are added, based on the Inoue (2011) emission-line ratio, through the same process as done in Salmon et al. (2015). The IGM attenuation was applied to model the galaxy SEDs, following Madau (1995).

Fiducial values of SED-derived physical properties, such as stellar masses, the absolute UV magnitudes ($M_{\text{UV}}$), and the UV continuum slope ($\beta$), were obtained from the best-fit models, which minimize $\chi^{2}$ to the observed photometry. We estimated the uncertainties of physical quantities from SED fitting with 1000 Monte Carlo (MC) realizations of the simulated photometric fluxes, which we perturbed the observed fluxes with their photometric errors. The 1$\sigma$ uncertainties denote the upper and lower limits of the central 68% range taken from the 1000 MC simulations. We repeated the process for all individual targets. We fixed galaxy redshifts with Ly$\alpha$-derived spectroscopic redshifts for emission-detected objects and with the best-fit photometric redshifts for non-detection objects. We derived $M_{\text{UV}}$ by averaging fluxes over a 100Å-bandpass (at the rest-frame 1450 – 1550Å) from SED models, which are not dust-corrected. The rest-frame UV continuum ($\beta$) was measured in the rest-frame UV bandpass of 1300 – 2600Å from the best-fit SED models as well, where $\beta$ is the spectral index in the form of $f_{\lambda}\propto\lambda^{\beta}$.

We present the $M_{\text{UV}}$ distribution of our targets in Figure 1. Our targets are somewhat biased toward UV-brighter galaxies, and a majority of the targets have $M_{\text{UV}}\lesssim-20$, comparable to the CANDELS/HST $J_{125}$ imaging depth in the EGS field. Figure 2 presents our targets in the $M_{\text{star}}$–$M_{\text{UV}}$ plane (top) and their rest-UV continuum slope ($\beta$) versus $M_{\text{UV}}$ (bottom). Although our sample contains limited coverage of UV-faint ($M_{\text{UV}}>-20$) sources, our spectroscopic targets are broadly consistent with the $z\sim 7$ $M_{\text{star}}$–$M_{\text{UV}}$ relation (Song et al., 2016b), which is representative of the typical high-redshift galaxy population. Also, we find the median value of the rest-UV continuum slopes at $\beta=-1.96$ from $M_{\text{UV}}<-20$ galaxies. This is comparable to the typical range of the UV slope measurements at this redshift (e.g., Finkelstein et al. 2012 find $\beta=-2.15^{+0.25}_{-0.16}$ for $M_{\text{UV}}<-20$ at $z=7$).

Spectroscopic observations of our sample were obtained over two nights in April 2021 using the Keck MOSFIRE spectrograph (McLean et al., 2012). This observing time was awarded through the NASA/Keck allocation (PI: I. Jung). We created four slitmask configurations that accommodate 61 high-redshift candidate galaxies for Ly$\alpha$ emission within a total effective sky area of $\sim 10^{\prime}\times 10^{\prime}$ (Figure 3). We observed two pointings each night, resulting in $\sim$3.5 hr of total exposure time per mask. We used the $Y$-band filter to cover Ly$\alpha$ at $7.0<z<8.2$. The spectral resolution of the $Y$-band filter is $\sim 3$Å ($R=3500$), and the slit width was set to be 0$\farcs$7, which corresponds to the typical seeing level at Mauna Kea. During the observations, individual science frames were taken with 180-sec exposures, and we used a standard ABAB dither pattern (+1$\farcs$25, -1$\farcs$25, +1$\farcs$25, -1$\farcs$25). The seeing level varies through the nights from 0$\farcs$7 to 1$\farcs$2. The observational details are listed in Table 1.

We used the recent version of the public MOSFIRE data reduction pipeline (DRP)⁴⁴4https://keck-datareductionpipelines.github.io/MosfireDRP/ to reduce the raw data. The public DRP provides a sky-subtracted, flat-fielded, and rectified two-dimensional (2D) slit spectrum per slit object. The reduced spectra are wavelength-calibrated using telluric sky emission lines. Reduced 2D spectra have the spectral resolution of 1.09Å pixel^-1 and the spatial resolution of 0$\farcs$18 pixel^-1.

It has been reported that there is significant slit drift in the spatial direction (up to $\sim$pixel hr^-1) in previous MOSFIRE observations (e.g., Kriek et al., 2015; Song et al., 2016a; Jung et al., 2019; Hutchison et al., 2020; Larson et al., 2022), which needs to be handled separately if observations last longer than a couple of hours of exposure time. However, the general use of the public DRP is not aimed to correct the known slit drift in the spatial direction for long-exposure science. We corrected the slit drifts found in our observations, following Jung et al. (2020). Briefly, we reduced each adjacent pair of science frames with the public DRP separately, generating the reduced 2D spectra of 360 sec exposure time. In our observations, we placed slits on two faint stars per slitmask for flux calibration and used them to trace the slit drifts in individual MOSFIRE pointings as well. The amount of slit drift is estimated by tracing the spatial positions of slit continuum sources on the DRP-reduced 2D spectra of 360-sec exposure. We corrected the measured slit drifts when combining 360-sec DRP-reduced 2D spectra to generate a single science frame for each slit target. Cosmic ray rejection and/or bad pixel cleaning are not feasible with the DRP runs on a pair of science frames. Thus, we cleaned them by taking sigma-clipped means in the 2D combination step. Also, to maximize a resulting signal-to-noise ratio (SNR), we weight the DRP-reduced frames with the Gaussian peak fluxes of the slit stars, which reflect observing conditions.

The one-dimensional (1D) spectra of our slit objects were extracted from the combined 2D spectra using an optimal extraction scheme (Horne, 1986) with a 1$\farcs$4 spatial window twice the typical seeing level of Mauna Kea. We model a spatial weight profile that follows the spatial profile of the slit stars, thus the pixels near the peak of the stellar spatial profile are maximally weighted. This enables us to correct the offsets of the actual spatial locations of slit objects from the expected positions, which are found up to a couple of pixels in a spatial direction.

The reduced 1D spectra were used to search for emission-line candidates. For the detection candidates, we repeated 1D extraction by shifting the centers of the optimal 1D extraction within $\pm$3 pixels from the corrected spatial locations of our slit objects. This accounts for the uncertainties in the centering of the objects’ spatial locations, allowing us to obtain maximum SNRs of the emission-line candidates.

For absolute flux calibration and telluric absorption correction, we used long-slit observations of a spectro-photometric standard star (HIP56147, the spectral type A0V) and Kurucz (1993) model stellar spectra. We estimated a wavelength-dependent response curve each night by dividing the model stellar spectrum with the reduced long-slit stellar spectra. The response curves were scaled to match the known photometric magnitude of the standard star. However, our science observations were obtained in different observing conditions, such as seeing and airmass, to the standard star observations. Thus, we estimated additional scaling factors using slit stars in the science slitmasks to refine the absolute flux calibration by matching their $Y$-band magnitudes measured from our spectra to the known $Y_{\text{105}}$ magnitudes from the existing HST photometry (Finkelstein et al., 2022). The additional scaling factors were at the level of $<10\%$. The slit losses due to the narrow slit width of 0$\farcs$7 in our observations are corrected in this step, considering the seeing conditions. We assume our high-redshift target galaxies are point sources, as they are unresolved in our observations. Overall, we obtained 3$\sigma$ detection limits of emission lines at $\sim$5$\times 10^{-18}$ erg s^-1 cm^-2 between sky-emission lines (with $\sim$3.5hr integration), and this is comparable to typical detection limits from previous MOSFIRE $Y$-band observations (e.g., Finkelstein et al., 2013; Song et al., 2016a; Jung et al., 2020).

We implemented an automated search scheme to capture plausible emission-line candidates consistently, similar to the method in Jung et al. (2020). We first collected emission-line candidates that were selected via our automated search on both 1D and 2D spectra by performing Gaussian line fitting on the 1D and Source Extractor (Bertin & Arnouts, 1996) runs on the 2D spectra. We required a 3$\sigma$ detection threshold in both the 1D and 2D searches. Then, we manually inspected individual emission-line candidates to rule out (i) sky-emission residuals, (ii) spurious sources, and (iii) contaminants from nearby sources. We conservatively removed emission-line features that are found close to the edge of sky-emission lines. To rule out spurious sources, we inspected the 2D spectra to ensure that there are clear negative peaks shown at the expected locations, $\pm$ 2$\farcs$5 apart from source positions, caused by the dither pattern of MOSFIRE. Additionally, we inspected the HST images to see if there are potential nearby contaminants whose emission lines could be captured at the same spatial locations of the MOSFIRE slits. Lastly, we performed tailored asymmetric (for the extended emission) or Gaussian (for the sharp/unresolved) emission-line fitting in reduced 1D spectra to calculate the line fluxes of emission-line candidates.

For emission-line-detected objects, we further checked their possibility of being low-redshift interlopers to ensure their nature as Ly$\alpha$. First, we checked if multiple emission lines are found in the same object. These would originate from a combination of emission lines from low-redshift galaxies (e.g., [O iii] $\lambda\lambda$4959, 5007; H$\beta$; [N ii] $\lambda\lambda$6548, 6584; H$\alpha$). We manually inspected the wavelengths of the possible companion lines, and we find no evidence of multiple emission lines in our sample. Second, we checked the possible low-redshift solution of being an [O ii] $\lambda\lambda$3727, 3729 emitter which can mimic the Lyman-break feature with the Balmer break of low-redshift galaxies. If that is the case, the [O ii] doublet should be resolved with the spectral resolution of Keck/MOSFIRE ($R=3500$ or $\sim$3Å). However, none of our emission-line candidates display the doublet emission lines with a $\sim$7–8Å separation (an expected peak separation of the [O ii] doublet at $z\sim$ 1.7–1.8).

We also compared the $\chi^{2}$ values from best-fit SEDs between high-redshift (with Ly$\alpha$) and low-redshift (with [O ii]) solutions for Ly$\alpha$-emission candidates as supplementary check, removing Ly$\alpha$-emission candidates disfavored in SED fitting analysis. The best-fit model SEDs of the host galaxies are shown in Figure 4. Our SED fitting analysis presents that the high-redshift solutions with Ly$\alpha$ are preferred over the low-redshift solutions for our Ly$\alpha$-detected galaxies in agreement with the photometric redshift PDFs.

From our emission-line search, we find eight galaxies with Ly$\alpha$ emission detected in the spectra (SNR $\geq$ 4). The 1D and 2D spectra of the detected emission lines are displayed in Figure 5. We caution that we are unable to completely rule out the chance of being low-redshift interlopers. The secondary lines that might have served to reject a high-redshift interpretation could be below the detection limit, particularly when coincident with a strong sky emission line. The strongly varying signal-to-noise can be seen in Figure 5. Nevertheless, based on our robust identification of emission lines and thorough diagnosis of low-redshift interlopers, we conclude that the detected emission lines are consistent with Ly$\alpha$ emission at $z>7$.

The process of reionization is expected to be inhomogeneous as overdensity regions with clustered galaxies were to be ionized earlier than field areas (e.g., Finlator et al., 2009; Mesinger et al., 2015; Katz et al., 2019). The currently most accessible tool to probe ionized structures in the middle of reionization is to search for Ly$\alpha$ emission from the reionization-era galaxies. A recent effort to spectroscopic search for Ly$\alpha$ resulted in the discoveries of clustered LAEs in the middle/late phase of reionization at $z\sim$ 7 – 8 (e.g., Zheng et al., 2017; Castellano et al., 2018; Tilvi et al., 2020; Jung et al., 2020; Endsley et al., 2021a). Particularly, the EGS field has a couple of known clustered structures with LAEs at $z\gtrsim 7.5$ (Zitrin et al., 2015; Oesch et al., 2015; Roberts-Borsani et al., 2016; Tilvi et al., 2020; Larson et al., 2022), and their additional membership candidates are discussed based on photometric selection (Leonova et al., 2022). Our MOSFIRE $Y$-band program detected Ly$\alpha$ emission from eight sources, and five of them are potentially associated with the known $z\sim 7.7$ LAE cluster (Tilvi et al., 2020)⁵⁵5Ly$\alpha$ from the brightest galaxy (z8_5) was first detected in Oesch et al. (2015) and Roberts-Borsani et al. (2016), and its C iii] emission was detected in Stark et al. (2017).. This demonstrates a possible extension of the LAE structure with up to eight LAEs at $z\sim 7.7$. This is currently the largest measured LAE cluster system in this early Universe at $z>7$. We summarize the $z\sim 7.7$ LAEs in Table 3.

With the discovery of the clustered LAEs at $z\sim 7.7$, we calculate 3-dimensional (3D) separations of individual $z\sim 7.7$ LAEs from the brightest (and potentially central) galaxy (z8_5). The left panel of Figure 7 shows the distribution of our $z\sim 7.7$ LAEs (green) in addition to the three $z\sim 7.7$ LAEs that were previously discovered (white; Oesch et al., 2015; Roberts-Borsani et al., 2016; Tilvi et al., 2020). The estimated 3D physical distances from z8_5 are listed in parentheses, ranging from 0.7pMpc at the nearest to 2.5pMpc at the farthest.

The crowd of ionizing sources could create the extended ionized structure beyond a $\sim$1pMpc scale of individual ionized bubbles and eventually enhance Ly$\alpha$ transmission in the IGM. To examine whether these clustered objects are situated in a connected ionized structure, we estimated the ionized bubble sizes which could be created by individual LAEs. Following Tilvi et al. (2020) and Jung et al. (2020), we used the relation between Ly$\alpha$ luminosities and ionized bubble sizes, predicted in the theoretical models from Yajima et al. (2018). Briefly, Yajima et al. (2018) model LAEs and the ionized bubble sizes based on individual halo merger trees using star formation history which is modeled to provide a reionization history consistent with the Planck observations (Planck Collaboration et al., 2016). Based on the Yajima et al. (2018) models, we used the measured Ly$\alpha$ luminosities to derive the predicted sizes of individual ionized bubbles around LAEs. The estimated ionized bubble sizes are ranging from $\sim$0.7 to 1.0 pMpc, as listed in Table 3. As the models in Yajima et al. (2018) predict the growth of isolated ionized bubbles around LAEs, it does not consider the additional expansion due to the overlapping ionized bubbles. Thus, the derived ionized bubble sizes may indicate the lower limits of ionized bubble sizes around LAEs; a much larger ionized structure could be created by the overlaps of multiple ionized bubbles around these LAEs. Such an extended ionized structure may promote Ly$\alpha$ escape from galaxies (e.g., Mason & Gronke, 2020; Park et al., 2021; Qin et al., 2021; Smith et al., 2021), resulting in enhanced Ly$\alpha$ detection rate in our observations at this redshift.

The detectability of Ly$\alpha$ emission in targeted spectroscopic observations is affected by target selection functions (which considers photometric redshift measurement PDFs and galaxy $M_{UV}$ distribution) and detection limits (depending on e.g., observing conditions and the presence of sky-emission lines) in addition to the Ly$\alpha$ IGM transmission particularly during the epoch of reionization. Thus, this makes it complicated to interpret a Ly$\alpha$ detection rate at its face value.

Instead, we performed Ly$\alpha$ EW distribution modeling to estimate the expected number of Ly$\alpha$ detections in our observations, which also consider target selection as well as observational conditions. Following Jung et al. (2020), in our Ly$\alpha$ EW modeling, we assume the Ly$\alpha$ EW distribution in its exponential functional form, $dN/dEW\propto\text{exp(-EW)}/W_{0}$, characterized with an $e$-folding scale ($W_{0}$). We populate mock Ly$\alpha$ emission lines for spectroscopic targets with (i) EW values that are randomly taken from the assumed EW distributions and (ii) wavelength locations, also randomly chosen based on galaxy photometric redshift PDFs. We calculate the expected detection rates above the detection limits of these simulated Ly$\alpha$ emission lines. To sum up, we quantify the Ly$\alpha$ detectability into the expected number of detections above detection limits as a function of $W_{0}$ via our EW modeling.

In the right panel of Figure 7, we compare the actual Ly$\alpha$ detection to what is estimated in our EW modeling. The figure shows the number of Ly$\alpha$ detection above detection limits per unit volume, a 1cMpc-think slice in the line-of-sight (LOS) direction in the sky area covered in our observations. We present the 1$\sigma$ range of the expected Ly$\alpha$ detections for the high and low Ly$\alpha$ EW cases with $W_{0}=200$Å and 20Å, respectively. The choice of $W_{0}=200$Å represents the distribution of extremely large EW Ly$\alpha$ emission lines whereas the low Ly$\alpha$ EW case of $W_{0}=20$Å is comparable to the statistical measurement of $W_{0}$ from $M_{\text{UV}}<-20$ galaxies in this redshift (Jung et al., 2022). In the figure, the redshift distribution of our actual Ly$\alpha$ detections are shown as the blue star symbols, and the black histogram shows the estimated detection number density per unit volume. The spike of our actual Ly$\alpha$ detection at $z\sim 7.7$ exceeds the expectation of the extreme case of the high Ly$\alpha$ EW distribution ($W_{0}=200$Å  red shades) whereas non/rare detections of Ly$\alpha$ at other redshift ranges are more consistent with the low Ly$\alpha$ EW case ($W_{0}=20$Å). Even without the three known LAEs of Tilvi et al. (2020), our Ly$\alpha$-detection-rate analysis demonstrates that we observe significantly stronger Ly$\alpha$ from the clustered galaxies compared to that from the rest of the galaxies.

As discussed above, Ly$\alpha$ transmission can be enhanced in an extended ionization structure formed by clustered galaxies because the Ly$\alpha$ transmissivity increases with the distance between the source and the boundary of the ionized region (e.g., Dijkstra, 2014). It is remarkable that the brightest member (z8_5) of the clustered LAEs has the lowest redshift in the cluster indicating that it is likely located in front of all the other LAEs toward the observer (see the right panel of Figure 8). That is, we may be witnessing a bright galaxy boosting the Ly$\alpha$ visibility preferentially on its rear side.

The left panel of Figure 8 illustrates how Ly$\alpha$ transmissivity can be boosted in the rear side of a bright galaxy. The bright galaxy may be located at the center of the local ionized bubble. The ionized bubble opens a wavelength window on the red side of Ly$\alpha$ where photons can avoid being resonantly scattered by the neutral IGM. Ly$\alpha$ photons from galaxies on the rear side experience a larger cosmological redshift before reaching the neutral region and have a better chance of transmitting through the window. Indeed, a simulation study by Park et al. (2021) reported that fainter galaxies tend to have a larger variation in Ly$\alpha$ transmission as their transmissivity sensitively depends on where they are located with respect to their brighter neighbors.

A peculiar motion of galaxies can contribute to the line shift as well. However, the redshift range of these clustered LAEs corresponds to $\gtrsim 1000$ km s^-1 in velocity, which is too large to be explained by the gravitational dynamics of galaxies alone. In theoretical studies, galaxies that are similar to z8_5 in UV magnitude have their total masses around $\sim 10^{12}M_{\odot}$ (e.g., Ocvirk et al., 2020). The gravitational infall velocity of such a galaxy peaks at $\sim$250–300 km s^-1 near the virial radius ($\sim$0.04 pMpc at $z=7.7$) and decreases as $r^{-0.5}$ with increasing distance to the galaxy ($r$). As the clustered LAEs spread out over a $\gtrsim$1pMpc distance, we conclude that the LAEs observed here are far from forming a dynamically relaxed system in this early Universe.

Recent surveys have produced a higher yield of Ly$\alpha$ detections from UV-luminous ($M_{\text{UV}}\lesssim-20$) galaxies than earlier attempts (e.g., Castellano et al., 2018; Jung et al., 2019, 2020, and in this work). In particular, spectroscopic observations of luminous sources with IRAC color excess, which reflects intense [O iii]+H$\beta$ emission, delivered a much higher Ly$\alpha$ detection rate (e.g., Oesch et al., 2015; Zitrin et al., 2015; Roberts-Borsani et al., 2016; Stark et al., 2017; Endsley et al., 2021b; Laporte et al., 2021). In contrast, Ly$\alpha$ has rarely been detected in follow-up spectroscopic observations of UV-faint ($M_{\text{UV}}\gtrsim-20$) galaxies (Hoag et al., 2019; Mason et al., 2019). Moreover, the unusually high detection rate of Ly$\alpha$ from galaxies with the IRAC excess is again somewhat diminished when it comes to targeting faint galaxies (Roberts-Borsani et al., 2022b). Additionally, Ly$\alpha$ emission has not been detected in the recent JWST NIRSpec observations which targeted UV-faint galaxies (Roberts-Borsani et al., 2022a; Williams et al., 2022; Morishita et al., 2022).

This notable difference in a Ly$\alpha$ detection rate between UV-luminous and UV-faint populations could be understandable in the sense that an inhomogeneous nature of reionization suggests an earlier process of reionization around overdense regions where bright galaxies are preferentially located (e.g., Finlator et al., 2009; Katz et al., 2019; Ocvirk et al., 2021; Kannan et al., 2022). Although Morishita et al. (2022) report non-detection of Ly$\alpha$ even from clustered galaxies, those galaxies are less luminous (with $M_{\text{UV}}\gtrsim-20$) than those studied here. Thus they may not create ionized regions large enough to allow the escape of Ly$\alpha$.