# Ly$\alpha$ Line Properties at $z=3.78$ and Their Environmental Dependence: A Case Study around a Massive Protocluster

Nicola Malavasi¹,², Kyoung-Soo Lee², Arjun Dey³, Rui Xue⁴, Yun Huang², and Ke Shi⁵

¹Université Paris-Saclay, CNRS, Institut d'astrophysique spatiale, F-91405, Orsay, France; nicola.malavasi@universite-paris-saclay.fr

²Department of Physics and Astronomy, Purdue University, 525 Northwestern Avenue, West Lafayette, IN 47907, USA

³NSF's National Optical-Infrared Astronomy Research Laboratory, 950 N. Cherry Ave., Tucson, AZ 85719, USA

⁴National Radio Astronomy Observatory, 520 Edgemont Road, Charlottesville, VA 22903, USA

⁵Department of Astronomy and Institute for Theoretical Physics and Astrophysics, Xiamen University, Xiamen, Fujian 361005, People's Republic of China

Received 2021 March 23; revised 2021 July 1; accepted 2021 August 9; published 2021 November 5

## Abstract

Ly$\alpha$-emitting galaxies (LAEs) are easily detectable in the high-redshift universe and are potentially efficient tracers of large-scale structure at early epochs, as long as their observed properties do not depend strongly on environment. We investigate the luminosity and equivalent width functions of LAEs in the overdense field of a protocluster at redshift $z \simeq 3.78$. Using a large sample of LAEs (many spectroscopically confirmed), we find that the Ly$\alpha$ luminosity distribution is well represented by a Schechter function with $log(L*/erg s^{-1}) = 43.26_{-0.22}^{+0.20}$ and $log(\phi*/Mpc^{-3}) = -3.40_{-0.04}^{+0.03}$ with $\alpha = -1.5$. Fitting the equivalent width distribution as an exponential, we find a scale factor of $\omega = 79_{-15}^{+15} \AA$. We also measured the Ly$\alpha$ luminosity and equivalent width functions using the subset of LAEs lying within the densest cores of the protocluster, finding similar values for $L^*$ and $\omega$. Hence, despite having a mean overdensity more than $2 \times$ that of the general field, the shapes of the Ly$\alpha$ luminosity function and equivalent width distributions in the protocluster region are comparable to those measured in the field LAE population by other studies at similar redshift. While the observed Ly$\alpha$ luminosities and equivalent widths show correlations with the UV continuum luminosity in this LAE sample, we find that these are likely due to selection biases and are consistent with no intrinsic correlations within the sample. This protocluster sample supports the strong evolutionary trend observed in the Ly$\alpha$ escape fraction and suggests that LAEs at lower redshift can be on average significantly more dusty than their counterparts at higher redshift.

*Unified Astronomy Thesaurus concepts: High-redshift galaxies (734); Protoclusters (1297); Luminosity function (942); Narrow band photometry (1088); High-redshift galaxy clusters (2007)*

## 1. Introduction

Studies of galaxy formation, evolution, and cosmology necessarily require that we investigate the high-redshift universe. The class of systems known as Ly$\alpha$ emitters (LAEs, Partridge & Peebles 1967) is particularly helpful to build samples of galaxies in the distant universe. LAEs are objects whose spectrum presents a fairly strong Ly$\alpha$ emission line with respect to the ultraviolet (UV) continuum. As the Ly$\alpha$ emission line is redshifted to optical wavelengths at $z > 2$, it is relatively easy to detect these galaxies in ground-based narrowband images. Large samples of LAEs have been selected photo-metrically and confirmed spectroscopically (see, e.g., Cowie & Hu 1998; Kudritzki et al. 2000; Ajiki et al. 2006; Gronwall et al. 2007; Finkelstein et al. 2009; Nilsson et al. 2009b; Guaita et al. 2010; Hayes et al. 2010; Ono et al. 2010; Ciardullo et al. 2012; Ning et al. 2020) at redshift greater than two, with a few samples also in the redshift range $z \sim 6-9$ (Hu et al. 2002, 2019, 2021; Kodaira et al. 2003; Iye et al. 2006; Willis et al. 2006; Ouchi et al. 2008; Ono et al. 2010; Zheng et al. 2017). In recent times several large spectroscopic or photo-metric surveys have been developed with the explicit aim of specifically targeting LAEs, such as the Hobby-Eberly Telescope Dark Energy Experiment (HETDEX, Hill et al. 2008) and its Pilot Survey (Adams et al. 2011), the Large Area Lyman Alpha survey (LALA, Rhoads et al. 2000; Malhotra & Rhoads 2002), the Multiwavelength Survey by Yale-Chile (MUSYC, Gawiser et al. 2006b), the Subaru Deep Survey, and Subaru/XMM-Newton Deep Survey (SXDS, Ouchi et al.

2003, 2008). Large LAE samples resulting from these programs have enabled systematic investigations of their physical properties and their evolutionary connection to the present-day galaxies.

Studies have shown that LAEs are a population of low-mass ($10^8-10^9 M_\odot$), star-forming ($1-10 M_\odot \text{ yr}^{-1}$) galaxies, abundant in the high-redshift universe ($z > 2$) (Cowie & Hu 1998; Ouchi et al. 2003, 2008; Gawiser et al. 2006a, 2006b, 2007; Nilsson et al. 2007, 2009b; Pirzkal et al. 2007; Lai et al. 2008; Finkelstein et al. 2009). They are a numerous class that contribute significantly to the cosmic star formation rate density (SFRD) at high redshift (see, e.g., Madau & Dickinson 2014).

The complex process that leads to the emission of Ly$\alpha$ radiation in galaxies can be studied through the correlations that are present between the strength of the Ly$\alpha$ emission line (e.g., its equivalent width, EW) and other properties that trace the various characteristics of a galaxy's interstellar medium (ISM; e.g., its kinematics, the amount and the distribution of neutral hydrogen gas and dust). For example, Ly$\alpha$ luminosity is expected to be correlated with a galaxy's star formation rate (SFR), but this is modulated by the effects of radiative transfer through a complex gaseous and dusty medium. Ly$\alpha$ EW is expected to correlate with ISM kinematics, gas column density, and dust reddening (Charlot & Fall 1993; Dijkstra et al. 2006, 2007; Hansen & Oh 2006; Verhamme et al. 2006, 2008; Adams et al. 2009; Laursen 2010; Zheng et al. 2010; Hashimoto et al. 2013), in line with observations (e.g., Shapley et al. 2003).
The statistics of the LAE population are very helpful in understanding the physical properties of these objects. Since the amount of Lyα radiation leaking out of a galaxy depends strongly on the dust distribution and the ISM physics, studying how LAEs in a given volume of space are distributed in Lyα luminosity and EW and how these quantities change with redshift and/or environmental density can help in understanding the evolution and the changes in the physical properties of the ISM of LAEs. The Lyα luminosity function (LF) is usually described with a Schechter (1976) function, in which two regimes can be identified: a power-law region at low luminosities, and an exponential decline for the high-luminosity tail. The low-luminosity power-law region is dominated by faint, mostly low-mass galaxies, and is both sensitive to the physical conditions in the ISM (see, e.g., Santos et al. 2004; Rauch et al. 2008) and subject to incompleteness biases. The high-luminosity tail is more completely sampled, but is made up of both star-forming galaxies and active galactic nuclei (see, e.g., Sobral et al. 2017). The Lyα LF and the EW distribution of LAEs (the latter usually described with an exponential or a Gaussian function) have been explored by many studies (see, e.g., Gronwall et al. 2007; Ouchi et al. 2008; Guaita et al. 2010; Konno et al. 2016; Herenz et al. 2019; Fossati et al. 2021). Still, the redshift evolution of the parameters describing the two functions is not well constrained, nor their possible dependence on large-scale variations in the environmental density.

In this paper, we study the LAE population in the overdense field of a spectroscopically confirmed protocluster at $z \approx 3.78$ and present estimates of its Lyα luminosity function and EW distribution. By focusing on a high-density region of space where star formation in galaxies has likely been accelerated, we can study whether the properties of the LAE population depend on environment (see, e.g., Dey et al. 2016; Lemaux et al. 2018). Moreover, we study the correlation among the Lyα luminosity, EW, and UV continuum luminosity for these sources in order to be able to derive conclusions on the physical properties of their ISM, with particular regard to their dust distribution.

The rest of the paper is organized as follows: in Section 2 we describe the sample used for our work and in Section 3 we describe the procedure we used to measure the Lyα luminosity function and EW distribution. In Section 4 we present our results on the LF, while in Section 5 we describe the environment inhabited by the LAEs in our sample. In Section 6 we analyze the correlation between Lyα luminosity, EW, and UV continuum luminosity. Finally, a discussion of the results and our conclusions are presented in Sections 7 and 8, respectively. Throughout this work we use a concordance cosmology with $\Omega_m = 0.27$, $\Omega_\Lambda = 0.73$, $H_0 = 70 \text{ km s}^{-1} \text{ Mpc}^{-1}$.

## 2. Data

In this work, we use an LAE sample at $z \approx 3.78$. This sample has already been thoroughly described in other papers (see Lee et al. 2014; Dey et al. 2016; Xue et al. 2017, for a detailed description). Here we briefly summarize the information that is most important for our analysis.

Our LAE sample is composed of galaxies in a region encompassing two protocluster structures (PC 217.96+32.3, hereafter simply referred to as “PCF”) at a redshift of $z \approx 3.786$ (see Lee et al. 2014; Dey et al. 2016) and those in the general field. This sample covers an area of roughly $1.2 \times 0.6 \text{ deg}^2$ (i.e., 8428 Mpc}^2$, see Table 1) at the southern end of the Boötes field of the NOAO Deep Wide-Field Survey (NDWFS, Jannuzi & Dey 1999). The optical data used to define the LAE sample were taken with the Mosaic 1.1 wide-field imaging camera at the Mayall 4 m telescope at the Kitt Peak National Observatory (KPNO).

The three broadband filters ($B_W R_I$) are centered at wavelengths of $\lambda_c = 4222, 6652$, and $8118$ Å while the narrowband filter (WRC4, KPN0 filter number k1024) has a central wavelength of $\lambda_c = 5819$ Å and FWHM of $\Delta\lambda = 42$ Å. The narrowband transmission samples the Lyα emission at redshifts of $z = 3.775-3.810$, corresponding to a line-of-sight distance of $l_{\text{field}} \approx 27$ Mpc. The volume sampled by these observations is derived as $V = A_{\text{eff}} \times l_{\text{field}}$ and is reported in Table 1.

We isolate LAEs as those galaxies with blue narrow-to-broadband colors (WRC4-R), which correspond to an excess of Lyα emission in the WRC4 band. The selection criteria are

$$ (WRC4 - R) < -0.8 \cap S/N(WRC4) \ge 7 \cap ((B_W - R) > 1.8 \cup S/N(B_W) < 2) \quad (1) $$

where S/N indicates the signal-to-noise ratio (measured within the isophotal area) in a given band. When these criteria are applied to our sources they produce a sample of 171 objects. Of these, 100 galaxies were observed with the Deep Imaging Multi-Object Spectrograph at the 10 m Keck II telescope (DEIMOS, see Faber et al. 2003), resulting in the confirmation of 62 galaxies and an estimated spectroscopic success rate of 89% (as measured by Dey et al. 2016).

### 2.1. Measuring Lyα Luminosity and Equivalent Width

We compute the Lyα luminosity and the rest-frame Lyα equivalent width of each LAE candidate based on the photometry following the procedure described in Xue et al. (2017). Briefly, we define the $Q$ parameter for each passband:

$$ Q(z, \beta) = \frac{\int e^{-\tau}(\lambda/\lambda_0)^{2+\beta}\lambda^{-1}R(\lambda)d\lambda}{\int \lambda^{-1}R(\lambda)d\lambda} \quad (2) $$

where $\lambda_0 = 1215.67(1+z)$ Å, $R(\lambda)$ is the total filter throughput, $\tau$ is the effective optical depth of the intergalactic medium (IGM) as given in Inoue et al. (2014), and $\beta$ is the slope of the UV continuum ($f_\lambda \propto \lambda^\beta$). Then, the integrated line flux $F_{\text{Ly}\alpha}$ and continuum flux density $f_{\text{cont}}$ are expressed as

$$ \begin{aligned} F_{\text{Ly}\alpha} &= \frac{B_{\text{BB}} B_{\text{NB}} (Q_{\text{BB}} f_{\text{AB,NB}} - Q_{\text{NB}} f_{\text{AB,BB}})}{Q_{\text{BB}} B_{\text{BB}} - Q_{\text{NB}} B_{\text{NB}}} \\ f_{\text{cont}} &= \frac{B_{\text{BB}} f_{\text{AB,BB}} - B_{\text{NB}} f_{\text{AB,NB}}}{Q_{\text{BB}} B_{\text{BB}} - Q_{\text{NB}} B_{\text{NB}}}. \end{aligned} \quad (3) $$

The $Q$ parameters and bandpass $B$ for each filter are given in Table 5 of Xue et al. (2017). In this equation, $f_{\text{AB}}$ denotes the flux density in units of erg s$^{-1}$ cm$^{-2}$ Hz$^{-1}$. The labels BB and NB refer to the broadband and narrowband filters that sample the Lyα wavelength, respectively. The monochromatic flux density and the Lyα rest-frame equivalent width are expressed in terms of the above quantities as

$$ f_{\text{AB}} = 10^{-0.4(m_{\text{AB}}+48.6)} = \frac{F_{\text{Ly}\alpha}}{B} + f_{\text{cont}} Q \quad (4) $$
<table><caption>**Table 1**
Summary of Relevant Quantities for the LAE Sample</caption><thead><tr><th>Quantity</th><th>Value</th></tr></thead><tbody><tr><td>Field center (J2000)</td><td>14:31:35.7 +32:20:05</td></tr><tr><td>Redshift</td><td>3.775–3.810</td></tr><tr><td>A<sub>eff</sub> [deg<sup>2</sup>]</td><td>0.526</td></tr><tr><td>N<sub>LAE</sub></td><td>171</td></tr><tr><td>l<sub>field</sub> [Mpc]</td><td>26.75</td></tr><tr><td>V [Mpc<sup>3</sup>]</td><td>2.255 × 10<sup>5</sup></td></tr></tbody></table>

**Note.** A<sub>eff</sub> is the effective area of the field, l<sub>field</sub> is the comoving line-of-sight distance corresponding to the narrowband filter transmission FWHM, and V is the effective comoving volume of the field defined as A<sub>eff</sub> × l<sub>field</sub>.

$$W_0 = \frac{W_{\text{obs}}}{1+z} = \frac{\lambda_0^2}{c} \times \frac{F_{\text{Ly}\alpha}}{(1+z)f_{\text{cont}}}. \quad (5)$$

In computing the line luminosity and EW for each LAE in our sample, we assume that all galaxies lie at the same redshift (taken to be z = 3.790) for the PCF LAEs. We assume all sources to be at the redshift corresponding to the center of the narrowband filter. Hence, the dominant source of uncertainty on the Lyα luminosity function will be the fact that some sources may be at the filter edge and not at the center of the filter itself (which we take into account by means of our simulations, see Section 2.2). The filter has FWHM = 42 Å, and therefore the redshift error resulting from our assumption is at most ±0.0173, which results in a luminosity error of <1%.

## 2.2. Completeness Estimate

Our ability to identify LAEs and measure their line properties depends on the imaging depth, the color selection method we adopt, as well as the redshift of the source. For example, brighter and compact LAEs are easier to identify than fainter and diffuse ones. For galaxies lying near the edge of the filter transmission curve, the EWs and line luminosities would be underestimated. Here, we describe the extensive image simulations we carry out in order to fully characterize the measurement systematics and their impact on the observed number counts and color distribution.

We generate a list in which, for each entry, a redshift, UV continuum slope β, Lyα luminosity, and EW are randomly chosen from the ranges W₀ = 0–350 Å, log L<sub>Lyα</sub> = 41.25–44.25, and β = -(1.5–2.0); for the redshift, we assume the range z = 3.730–3.853. The assumed parameter ranges are much wider than those expected in real data, to ensure that robust statistics can be obtained for the galaxies near the selection cutoffs (in source detection, color selection, and redshift). In particular, the chosen redshift range is much larger than the width of the WRC4 filter.

For each entry, we create a spectrum; the base continuum spectrum is generated using the Bruzual & Charlot (2003) stellar population synthesis model assuming a constant star formation history and a population age of 100 Myr. A Salpeter (1955) initial mass function and solar metallicity are assumed. The continuum spectrum is normalized to the input flux density at the rest-frame 1700 Å; a Lyα emission line normalized to match the line luminosity is added as a Gaussian profile centered at 1215.67 Å and with an intrinsic line width of 3 Å. Our results are insensitive to the assumed line width as long as they reasonably reproduce the galaxy colors and line widths as

observed. As mentioned in Lee et al. (2014), the line flux is computed independently from the continuum spectrum, to take into account the lower resolution of the Bruzual & Charlot (2003) galaxy templates. The spectrum is then attenuated using the Inoue et al. (2014) prescription to account for IGM absorption by neutral hydrogen. Finally, synthetic photometry is computed by convolving each mock spectrum with the throughputs of the relevant passbands.

Based on the synthetic photometric catalog we generated above, we carry out image simulations as follows. In each run, we randomly draw 100 entries from the catalog and insert them into the real images. We assume a log-normal distribution for their angular sizes as measured from high-redshift observations (Bouwens et al. 2004; Ferguson et al. 2004; Law et al. 2012; Shibuya et al. 2015; Ribeiro et al. 2016). However, our results are insensitive to the morphologies and sizes because all are essentially point sources in our ground-based data. The images are convolved with the appropriate point-spread function and added to the real data. A total of 1.5 × 10⁶ sources are simulated in the PCF data set. Source detection, photometry, LAE selection, and derivation of line luminosities and EWs are performed in the identical manner as in the real data.

In Figure 1, we illustrate how both Lyα EWs and luminosities are recovered in our simulations of the PCF data set when the results are split into the “UV-bright” and “UV-faint” categories (m<sub>1700</sub> = 23.40 AB mag is used for the cut). This figure shows the recovered (output) Lyα luminosity and EW as a function of the input ones. It is apparent that UV-bright sources are more likely to be recovered with the correct Lyα luminosity and EW values. In comparison, the same values for UV-faint sources scatter over a wide range.

Using the image simulations we account for the effects of the LAE color selection, detection completeness, and photometric scatter on the observed galaxy statistics. To this end, we define the likelihood that a galaxy with intrinsic luminosity and equivalent width [L<sub>i</sub>, W<sub>j</sub>] has of being identified as an LAE and of being observed with the observed quantities [L<sub>i'</sub>, W<sub>j'</sub>] as

$$P(L_{i'}, W_{j'}|L_i, W_j) = \frac{N_{ij}^{\text{OUT}}(L_{i'}, W_{j'})}{N_{ij}^{\text{IN}}} \quad (6)$$

where N<sub>ij</sub><sup>IN</sup> is the number of sources simulated in the range [L<sub>i</sub> – ΔL/2, L<sub>i</sub> + ΔL/2] and [W<sub>j</sub> – ΔW/2, W<sub>j</sub> + ΔW/2]. The quantity N<sub>ij</sub><sup>OUT</sup>(L<sub>i'</sub>, W<sub>j'</sub>) denotes the number of sources identified as LAEs in the same bin measured with the quantities [L<sub>i'</sub> – ΔL/2, L<sub>i'</sub> + ΔL/2] and [W<sub>j'</sub> – ΔW/2, W<sub>j'</sub> + ΔW/2]. Throughout this work, all intrinsic quantities are denoted as unprimed running indices, while all observed quantities are shown with primed indices.

Figure 2 reports an estimate of the ratio of observed to expected LAE counts in the Lyα–EW space. The distributions shown in Figure 2 are obtained when Equation (6) is summed over all intrinsic Lyα luminosity and EW values (i, j indices) and shown as a function of the observed Lyα luminosity and EW values (i', j' indices) for continuum-bright and continuum-faint sources. The sum of the probability from Equation (6) over the intrinsic values describes how the detected fraction at a given Lyα luminosity depends on the photometric uncertainties. For the majority of continuum-bright sources at least 30% of the sources in a given intrinsic Lyα luminosity and EW bin are observed as belonging to the same bin. Moreover, the distribution in the observed space of the bins for which the
<figure><img src="image_1.png" /><figcaption>Figure 1. Relation between input and output Lyα luminosity and EW. The color coding reflects the fraction of all simulated galaxies inside each bin and is consistent between the left and right columns. The left column is for UV-luminous simulated galaxies and right column is for UV-faint ones, where the distinction between the two samples is made considering those below and above the median UV apparent magnitude of the full simulated sample, $m_{1700,median} = 23.40$ AB mag.</figcaption></figure>

<figure><img src="image_2.png" /></figure>

probability of observing an LAE with the correct Lyα luminosity and EW is high is rather compact, restricted to a smaller range of observed Lyα luminosities and EWs than in the continuum-faint case.

For continuum-faint sources, several bins are also present where we observe more galaxies than expected based on the intrinsic number (i.e., an observed fraction greater than one). This is due to galaxies being scattered into a given Lyα luminosity and EW bin from neighboring ones due to photometric uncertainties. Moreover for continuum-faint sources, the number of bins with a very low fraction of observed simulated galaxies is larger and scattered also to higher EW values and Lyα luminosities.

### 3. Methodology

In this work, we aim to determine the intrinsic distributions of Lyα luminosity and equivalent width. Following the convention in the literature (e.g., Gronwall et al. 2007; Guaita et al. 2010; Blanc et al. 2011; Ciardullo et al. 2012), we model the luminosity function and EW function as a Schechter (1976) function and an exponentially declining function, respectively. They are expressed as

$$ \Phi(L_i) = \ln(10)\log(\phi^*)e^{-L_i/L^*}\left(\frac{L_i}{L^*}\right)^{\alpha+1} \quad (7) $$

$$ \Psi(W_j) = \psi^*e^{-W_j/\omega}. \quad (8) $$

In the above equations, both functions represent comoving number densities in a unit bin of log-luminosity $[\log(L_i), \log(L_i) + \Delta \log(L)]$ and EW $[W_j, W_j + \Delta W]$. $L^*$ is the characteristic luminosity while $\alpha$ and $\phi^*$ denote the faint-end slope and the normalization parameter, respectively. The

equivalent width function (EWF) is a single-parameter function characterized by the width $\omega$ as the normalization itself $\psi^*$ is tied to the LF.

Given the true LF and EWF, we can predict the expected luminosity and EW distributions “as observed” by essentially convolving the observational selection and biases with these quantities as

$$ N_{\text{exp}}(L_{i'}, W_{j'}) = \sum_{i,j} P(L_{i'}, W_{j'}|L_i, W_j) \Phi(L_i) \Delta \log(L) \Psi(W_j) V \quad (9) $$

where $V$ is the comoving volume of the survey as listed in Table 1. By comparing this expectation with the observed distribution, one can determine the best-fit parameters $L^*, \alpha, \phi^*, \text{and } \omega$.

To this end, we employ a maximum likelihood method similar to that presented in Bouwens et al. (2007). The log-likelihood is defined as

$$ p_{i'j'} \equiv N_{\text{obs}}(L_{i'}, W_{j'}) \ln \left( \frac{N_{\text{exp}}(L_{i'}, W_{j'})}{\sum_{m,n} N_{\text{exp}}(L_m, W_n)} \right) \quad (10) $$

where the sum in the denominator is performed for all $m, n$ while the $p_{i'j'}$ values are computed only on the set of indices $A = \{i', j': N_{\text{obs}}(L_{i'}, W_{j'}) \neq 0\}$. The log-likelihood is then computed as $\mathcal{L} = \sum_{i', j' \in A} p_{i'j'}$. The best-fit solution is the one that maximizes the log-likelihood $\mathcal{L}$.

Our fitting method departs from the convention in the literature in that we make full use of the 2D distribution rather than fitting the Lyα luminosity and EW functions separately. If the two are related in a known fashion, the methodology can be easily generalized to fit such a model. In the case of no intrinsic
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Two-dimensional distribution of the integrated probability of observing an LAE in a given output Lyα luminosity-EW bin. The color scale parameterizes the fraction of simulated LAEs that will be observed in a given output Lyα luminosity-EW bin with respect to the number of LAEs originally in the same input bin. Note that this fraction can be higher than one if more LAEs are scattered into a bin than out of it. The left panel is for UV-luminous simulated galaxies and the right panel is for UV-faint ones, where the distinction between the two samples is made considering those below and above the median UV apparent magnitude of the full simulated sample, $m_{1700,median} = 23.40$ AB mag.</figcaption></figure>

relation between the two quantities, our method should yield a result equivalent to the conventional method. While we proceed assuming that there is no correlation, in Section 6, we return to the possibility that the two quantities are linked.

A wide range of parameters are considered in the fitting procedure: $log(L^*/(erg s^{-1})) = 40-45$, $\alpha = -(1.50-2.25)$, and $\omega = (1-700)$ Å, respectively. Using the best-fit parameters, we compute $N_{best}(L_{i'}, W_{j'})$ and define the expected number densities per log-luminosity and per EW as

$$ \Phi_{\text{exp}}(L_i') = \frac{\sum_{j'} N_{\text{best}}(L_{i'}, W_{j'})}{\Delta \log(L) \times V}; \\ \Psi_{\text{exp}}(W_{j'}) = \frac{\sum_{i'} N_{\text{best}}(L_{i'}, W_{j'})}{\Delta W \times V} \quad (11) $$

and similarly for the observed number densities, replacing $N_{best}(L_{i'}, W_{j'})$ with $N_{obs}(L_{i'}, W_{j'})$.

Finally, the normalization parameters for both functions ($\phi^*$, $\psi^*$) can be obtained as

$$ \phi^* = \frac{\sum_{i'} \Phi_{\text{obs}}(L_{i'})}{\sum_{i'} \Phi_{\text{exp}}(L_{i'})}; \quad \psi^* = \frac{\sum_{j'} \Psi_{\text{obs}}(W_{j'})}{\sum_{j'} \Psi_{\text{exp}}(W_{j'})}. \quad (12) $$

The confidence levels of the parameters are estimated assuming that the likelihood ratio $\Lambda = -2(\mathcal{L} - \mathcal{L}_{max})$ obeys a $\chi^2$ distribution with a number of degrees of freedom equal to the number of free parameters in the fit.

## 4. Lyα Luminosity and Equivalent Width Functions

Following the procedure described in Section 3, we measure the LF and EWF using the full LAE sample at $z \approx 3.8$. In Figure 3, we show our measurements together with the best-fit models. In the top panel, both data and models show the comoving number densities “as observed,” while in the bottom panel the same quantities are shown as “intrinsic.” The two are related through Equation (10).

We are unable to place robust constraints on the faint-end slope, $\alpha$. Formally, the best-fit $\alpha$ values is $-1.5$. This value is in a reasonable agreement with other measures within uncertainties (e.g., Sobral et al. 2017). However, all $\alpha$ values listed in Table 2 provide similarly good fits to the data as

<figure><img src="image_5.png" /><figcaption>**Figure 3.** The comoving number density of LAEs as observed (top) and as corrected for measurement and selection biases (bottom). Top: both data (red) and best-fit model (blue) include measurement and selection biases as described and characterized in Section 2.2. Bottom: observed densities corrected for incompleteness are shown in red; the green line represents the best-fit Schechter function.</figcaption></figure>

illustrated by the $\mathcal{L}_{max}$ values. We instead opt to fix $\alpha$ to five different values, linearly spaced in the range $\alpha = -(1.50-2.25)$, and determine $L^*$, $\Phi^*$, and $\omega$ for each case. The results are tabulated in Table 2 listed under “PCF.” For ease of comparison with other results in the literature, we fix the faint-end slope to $\alpha = -1.50$ for the $z \approx 3.8$ sample.

In Figure 4, we show our measurements together with those at $z \approx 3.7$ (Ouchi et al. 2008; Cassata et al. 2011). The bright
<table><caption>**Table 2**
Best-fit Parameters for the Rest-frame Ly$\alpha$ Luminosity and Equivalent Width Functions</caption><thead><tr><th>Field</th><th>α</th><th>log(L*/erg s<sup>-1</sup>)</th><th>log(φ*/Mpc<sup>-3</sup>)</th><th>ω (Å)</th><th>L<sub>max</sub></th></tr></thead><tbody><tr><td rowspan="6">PCF<br>171 LAEs</td><td>-1.50</td><td>43.26<sup>+0.20</sup><sub>-0.22</sub></td><td>-3.40<sup>+0.03</sup><sub>-0.04</sub></td><td>79<sup>+15</sup><sub>-15</sub></td><td>-620.4</td></tr><tr><td>-1.65</td><td>43.37<sup>+0.34</sup><sub>-0.23</sub></td><td>-3.57<sup>+0.03</sup><sub>-0.04</sub></td><td>79<sup>+16</sup><sub>-14</sub></td><td>-620.9</td></tr><tr><td>-1.70</td><td>43.37<sup>+0.47</sup><sub>-0.20</sub></td><td>-3.60<sup>+0.03</sup><sub>-0.04</sub></td><td>79<sup>+16</sup><sub>-14</sub></td><td>-621.1</td></tr><tr><td>-1.75</td><td>43.47<sup>+0.53</sup><sub>-0.26</sub></td><td>-3.74<sup>+0.03</sup><sub>-0.04</sub></td><td>79<sup>+17</sup><sub>-14</sub></td><td>-621.2</td></tr><tr><td>-2.00</td><td>44.08<sup>+0.92</sup><sub>-0.62</sub></td><td>-4.61<sup>+0.03</sup><sub>-0.04</sub></td><td>79<sup>+19</sup><sub>-13</sub></td><td>-622.3</td></tr><tr><td>-2.25</td><td>...</td><td>-6.11<sup>+0.03</sup><sub>-0.03</sub></td><td>79<sup>+22</sup><sub>-11</sub></td><td>-624.3</td></tr><tr><td rowspan="6">PCF, HD<br>43 LAEs</td><td>-1.50</td><td>43.26<sup>+0.97</sup><sub>-0.31</sub></td><td>...</td><td>72<sup>+40</sup><sub>-21</sub></td><td>-149.2</td></tr><tr><td>-1.65</td><td>43.47<sup>+1.53</sup><sub>-0.43</sub></td><td>...</td><td>72<sup>+42</sup><sub>-20</sub></td><td>-149.2</td></tr><tr><td>-1.70</td><td>43.47<sup>+1.53</sup><sub>-0.41</sub></td><td>...</td><td>72<sup>+43</sup><sub>-20</sub></td><td>-149.2</td></tr><tr><td>-1.75</td><td>43.57<sup>+1.43</sup><sub>-0.47</sub></td><td>...</td><td>72<sup>+43</sup><sub>-20</sub></td><td>-149.2</td></tr><tr><td>-2.00</td><td>...</td><td>...</td><td>72<sup>+47</sup><sub>-19</sub></td><td>-149.3</td></tr><tr><td>-2.25</td><td>...</td><td>...</td><td>79<sup>+43</sup><sub>-25</sub></td><td>-149.8</td></tr><tr><td rowspan="6">PCF, Structure<br>72 LAEs</td><td>-1.50</td><td>43.26<sup>+0.59</sup><sub>-0.24</sub></td><td>...</td><td>72<sup>+24</sup><sub>-17</sub></td><td>-258.8</td></tr><tr><td>-1.65</td><td>43.47<sup>+1.53</sup><sub>-0.35</sub></td><td>...</td><td>72<sup>+25</sup><sub>-16</sub></td><td>-259.1</td></tr><tr><td>-1.70</td><td>43.57<sup>+1.43</sup><sub>-0.43</sub></td><td>...</td><td>72<sup>+26</sup><sub>-16</sub></td><td>-259.2</td></tr><tr><td>-1.75</td><td>43.57<sup>+1.43</sup><sub>-0.39</sub></td><td>...</td><td>72<sup>+26</sup><sub>-16</sub></td><td>-259.3</td></tr><tr><td>-2.00</td><td>...</td><td>...</td><td>72<sup>+29</sup><sub>-16</sub></td><td>-260.0</td></tr><tr><td>-2.25</td><td>...</td><td>...</td><td>72<sup>+31</sup><sub>-15</sub></td><td>-261.3</td></tr></tbody></table>

**Notes.** Columns are: the sample for which the LF and EWF have been computed (Field, with PCF referring to the full LAE sample, PCF HD referring to the LAEs in high densities defined through the Delaunay tessellation, and PCF Structure referring to LAEs within spectroscopically identified structures), the α, L*, and φ* parameters of the LF, the ω parameter of the EWF, and the maximum likelihood parameter of the fit (L<sub>max</sub>). Below the name of the sample is the number of LAEs used to constrain the fit.

a In this case the fitting procedure is not able to constrain the L* parameter.

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Comparison of measured Ly$\alpha$ LFs at z ≈ 3.8 with those reported by Ouchi et al. (2008) and Cassata et al. (2011). Dashed lines represent the best-fit Schechter function as tabulated in Table 2. Data from the literature have been corrected to the cosmology adopted in this paper. We find that the overall normalization of our LAE sample is a factor of ~2.0–2.5 higher than that in the average field at the bright end.</figcaption></figure>

end of our LF is not well sampled, due to a small sample size, leading to large uncertainties in characteristic luminosity L*. Nevertheless, there is a clear enhancement of the overall amplitude of the LF (of ~0.3 dex) compared to the Ouchi et al. (2008) determination based on field LAEs at the same redshift, in excellent agreement with our earlier estimate that the observed number of LAEs in the entire survey field is nearly twice the average expected in the field environments (Lee et al. 2014). Although our high-luminosity end is not well sampled,

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Distribution of Ly$\alpha$ equivalent widths are shown for the data in red for our LAE sample. The best-fit EWF is shown in green, together with the expected EW distribution in blue given the selection effects and photometric scatter (see Equation (6)).</figcaption></figure>

there is also an indication for the shape to be different from that measured by Ouchi et al. (2008), with an enhancement in our LF with respect to that of field LAEs. If confirmed, this could point toward an enhancement and an acceleration of galaxy formation in protocluster cores, as suggested, e.g., by Shimakawa et al. (2018) and Ito et al. (2020), albeit through analyses performed with different kinds of data sets.

In Figure 5, we show the observed EW distribution for our sample along with the best-fit intrinsic EWF (green). The EW distribution is well fitted by an exponential function, which shows how the majority of sources have EW < 250 Å with a tail extending up to 350 Å and beyond.

As illustrated in Figure 6, our measurements are in good agreement with those in the literature at similar redshift (van
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Best-fit LF and EWF parameters for our LAE sample are shown in red. The contours represent $1\sigma$, $2\sigma$, and $3\sigma$ confidence levels. The faint-end slope $\alpha$ is fixed to the value assumed by Ouchi et al. (2008), which is -1.50 for our $z \approx 3.8$ sample. Other colors represent the measurements from the existing literature at $z = 3-5$: Ouchi et al. (2008, O08), Gronwall et al. (2007, Gr07), Ciardullo et al. (2012, C12), Blanc et al. (2011, B11), van Breukelen et al. (2015, vB05), Sobral et al. (2018, So18), and Guo et al. (2020, Guo20). Letters in the legend refer to they way in which $\omega$ has been computed, whether it is from a Gaussian or exponential fit to the distribution (G or E, respectively), the average of the distribution (A), or a lower limit (L). The limits are displayed as arrows whenever appropriate. These literature data are corrected for different cosmology whenever applicable.</figcaption></figure>

Breukelen et al. 2005; Gronwall et al. 2007; Ouchi et al. 2008; Blanc et al. 2011; Ciardullo et al. 2012; Sobral et al. 2018; Guo et al. 2020) and are consistent with the limits set by Sobral et al. (2018). In the case of the EWF parameter $\omega$, Ouchi et al. (2008) parameterized the EW distribution with a Gaussian function, whereas we model it with an exponential function. The value of $\omega$ we report in Figure 6 for Ouchi et al. (2008) is therefore the standard deviation of a Gaussian (and not the e-folding of an exponential), therefore it is in principle not directly comparable to our measurement of $\omega$. However, we checked that values extracted from a range of Gaussian functions centered on zero and with standard deviations corresponding to the range covered by our parameter space are also well fitted by exponential functions with e-folding parameters ranging over the same range. Therefore, our measurement for $\omega$ and the one from Ouchi et al. (2008) can be directly compared. The contours in the $L^*-\omega$ space represent the $1\sigma$, $2\sigma$, and $3\sigma$ confidence levels; it is evident that the small size of the PCF $z \approx 3.8$ sample translates to large uncertainties, particularly for the determination of the characteristic luminosity $L^*$.

## 5. The Impact of the Environment on Ly$\alpha$ Emitters

As outlined in Section 2, our sample contains 54 spectroscopically confirmed members of two $z \approx 3.78$ structures (PC217.96+32.3-C and PC217.96+32.3-NE, hereafter PCF; see Lee et al. 2014; Dey et al. 2016). The estimated total masses of these structures are $\approx 1 \times 10^{15} M_{\odot}$ and $6 \times 10^{14} M_{\odot}$, respectively.

We quantify the environment traced by the LAEs at $z \approx 3.78$ by performing a two-dimensional Delaunay tessellation (Zaninetti 1990; van de Weygaert 1994; Marinoni et al. 2002), a geometrical dual operation of the Voronoi tessellation (see Marinoni et al. 2002). The line-of-sight distance sampled by the narrowband filter WRC4 is $\approx 20$ Mpc, roughly

comparable to that of the angular extent of the LAE overdensity. We use a 2D approach because the precise positions of individual LAEs are unknown in the line-of-sight direction, due to the fact that Ly$\alpha$ line centroids are typically redshifted with respect to systemic redshift, with the offset sensitive to the geometry, motion, and physical properties of interstellar gas.

The 2D plane is divided into a series of triangles having LAEs as vertices with a property that, for every triangle, the enclosing circle uniquely defined by its vertices does not contain any other galaxy. Unlike in the case of Voronoi tessellation, the polygons that tessellate the space are all triangles. Each galaxy is connected with a number of triangles $N_T$ each of area $A_t$ with $t = 0, 1, \dots, N_T$. The surface density of each LAE is determined as an average of the inverse of the areas of the triangles:

$$ \Sigma_D = \frac{1}{N_T} \sum_t \frac{1}{A_t}. \qquad (13) $$

The surface overdensity is then defined as

$$ \delta = \frac{\Sigma_D - \langle \Sigma \rangle}{\langle \Sigma \rangle} \qquad (14) $$

where $\langle \Sigma \rangle$ is the mean LAE density computed by dividing the total number of LAEs by the effective field area, which is listed in Table 1.

The use of the Delaunay tessellation is a slightly different approach than the one adopted by Dey et al. (2016), which relied on Voronoi tessellation. In this case, we wanted to take advantage of the fact that the measurement of the density field through the Delaunay tessellation allows for a slight smoothing of the density at the position of each galaxy, by means of the average of the areas of the triangles that connect it to the
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Delaunay tessellation of the $z \approx 3.8$ sample. Black lines show the tessellation polygons while circles represent the positions of candidate LAEs in the field. Galaxies are color-coded according to the percentile of the estimated LAE surface density. Two gray circles mark the positions of PC217.96+32.3C and PC217.96+32.3-NE (Lee et al. 2014), and have radii of 25 Mpc and 15 Mpc (comoving), respectively.</figcaption></figure>

neighboring galaxies. In this way, the density at the position of each galaxy is measured taking into account the density at the position of the neighboring ones, because the same triangle will contribute to the density measurement of multiple LAEs.

In Figure 7, we show the Delaunay tessellation of the PCF field, where LAEs are color-coded according to the percentile of the surface density. The positions of the two spectroscopically confirmed protoclusters are marked by dashed gray circles (with radii of 25 Mpc and 15 Mpc, comoving, for the central and NE structure, respectively); the majority of the LAEs therein are at the 80th percentile in surface overdensity. The radii of the overdensities are very approximate and are based on visual estimates (see Lee et al. 2014). We define the high-density LAE sample in two different ways: one according to the tessellation map and the other to the spectroscopic overdensity. These samples will be referred to as “HD” and “Structure” samples, respectively. The “HD” sample is defined as those LAEs with a density greater than the 75th percentile of the density distribution (43 LAEs belong to this sample), the

“Structure” sample as those LAEs within a radius of 25 and 15 comoving Mpc from the positions of the central and NE overdensities, respectively (72 LAEs belong to this sample).

We determine the best-fit LF and EWF parameters for the high-density LAE samples. The faint-end slope is fixed to the same values as the full sample. The results are listed in Table 2 under “Structure” for the spectroscopy-defined overdensity sample and “HD” for the tessellation-based sample. At $z=3.78$, the characteristic luminosity and EWF width (assuming $\alpha = -1.5$) are $\log L^* = 43.26_{-0.22}^{+0.20}$, $\omega = 79_{-15}^{+15} \AA$ for the full sample, $\log L^* = 43.26_{-0.24}^{+0.59}$, $\omega = 72_{-17}^{+24} \AA$ for the “HD” sample, and $\log L^* = 43.26_{-0.31}^{+0.97}$, $\omega = 72_{-21}^{+40} \AA$ for the “Structure” sample, respectively. All parameters obtained for the high-density LAEs are consistent with those of the full sample within errors. It appears that there is no significant difference in the formation process of LAEs in the average and overdensity field.

We chose to compare the LF and EWF in overdensities with the total LF and EWF because we wanted to check whether the LAE population in very dense environments showed any difference in the Ly$\alpha$ luminosity distribution with respect to the general LAE population. As the sky field where observations of our protocluster have been performed is rather small and essentially dominated by the presence of the PCF structures, we do not consider the option of constructing a sample of low-density LAEs representative of the field population as viable.

## 6. The Correlation between Ly$\alpha$ Luminosity, EW, and UV Continuum Luminosity

In our analyses thus far, we have assumed that there is no statistical correlation between Ly$\alpha$ equivalent width and line luminosity. By doing so, we have obtained results that are consistent with existing studies in the literature. However, if an intrinsic correlation exists between a galaxy’s UV luminosity ($M_{UV}$), line luminosity ($L_{Ly\alpha}$), and EW ($W_0$), it can substantially bias the best-fit LF and EWF parameters.

Indeed, Ando et al. (2006) noted a significant deficiency of high-EW sources in UV-luminous Lyman break galaxies (LBGs) at $z \sim 4-5$. Subsequent studies found a similar trend not only for LBGs (Vanzella et al. 2009; Oyarzún et al. 2017) but also for LAEs (e.g., Shioya et al. 2009); possible causes include changing dust geometry, metallicity, and population age with the galaxy’s UV luminosity (and with SFR), which alter the production rate of Ly$\alpha$ photons or modulate the likelihood of their escape through the interstellar medium (e.g., Kobayashi et al. 2010).

Given that the “Ando effect” is statistical in nature, a rigorous examination of such a trend, or the lack of one, requires a large sample of galaxies covering a sufficiently wide dynamic range in Ly$\alpha$ and UV continuum luminosity. Nilsson et al. (2009a) conducted a series of Monte Carlo realizations in which different scaling laws between the two were used to populate the galaxies according to their luminosity function. By comparing their result with the distribution of 232 LAEs and 128 LBGs at $z \sim 3$ in the EW-$L_{Ly\alpha}$ space, they concluded that no correlation is necessary to explain the observed distribution of galaxy luminosities and EWs. Furthermore, they argued that the lack of sources with high EW and high $L_{UV}$ is simply due to them being rare objects at the extreme end of both line and continuum luminosities. Based on their study of 130 LAEs at $z \sim 3$, Ciardullo et al. (2012) stated that photometric scatter plays a crucial role in creating an artificial correlation between
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Distribution of real and simulated LAEs on the planes defined by observable quantities. Top: UV continuum luminosity-EW plane. Bottom: UV continuum luminosity-Lyα luminosity plane. In both panels, cyan points represent real LAEs from our PCF ($z \sim 3.790$) sample, orange points represent random extractions from our simulated LAE sample. Black contours encompass the given percentage of the simulated LAE distribution. The agreement between the simulations and real data suggests that the observed LAE sample is consistent with the Lyα luminosity and EW distributions being independent of UV continuum magnitude.</figcaption></figure>

$M_{\text{UV}}$ and EW because they are anticorrelated, with EW decreasing with increasing $M_{\text{UV}}$ (when other galaxy properties are fixed).

Having measured the relevant quantities (EW, $L_{\text{Ly}\alpha}$, and $M_{\text{UV}}$) of our 171 LAEs and having characterized how photometric scatter affects these measurements through our simulations (Sections 2.1 and 2.2), we revisit the possibility of any intrinsic correlation between $M_{\text{UV}}$ and $W_0$.

We begin by evaluating the agreement between the observed and simulated data in the case of no correlation in order to establish the baseline. From our simulated LAE catalog (Section 2.2), we randomly extract sources matching the number of real LAEs. For each simulated galaxy (with values $L_{\text{Ly}\alpha,j}$, $W_{0,j}$), the product $\Phi_{\text{exp}}(L_{\text{Ly}\alpha,j}) \cdot \Psi_{\text{exp}}(W_{0,j})$ determines the likelihood of extraction, where we assume the best-fit $\Phi(L_{\text{Ly}\alpha})$ and $\Psi(W_0)$ for the full sample as given in Table 2; the procedure ensures that a mock galaxy sample is characterized overall by identical number counts and color distributions to the real LAE sample within Poisson uncertainties. The extraction is repeated 1000 times. In Figure 8, we show the positions of real LAEs (blue) and the distribution of simulated sources (orange dots; black contour lines enclose given

percentages of all simulated galaxies). It is evident that the observed galaxy distribution is similar to the simulated one.

We bin both real and simulated data using a bin size of $\Delta M_{\text{UV}} = 0.47$ mag, $\Delta W_0 = 0.2$ Å, and $\Delta \log L_{\text{Ly}\alpha} = 0.15$ for PCF. The simulated counts are normalized to match the real data, and thus can be considered as the expected value. The reduced chi-square $\chi_r^2$ (indicated in Figure 8) is close to unity in all cases, confirming our visual impression of a good agreement.

In both $M_{\text{UV}}-W_0$ and $M_{\text{UV}}-L_{\text{Ly}\alpha}$ planes, the two parameters are clearly correlated in a similar manner to that shown in Figure 6 of Ciardullo et al. (2012); the lack of UV-luminous high-EW sources is also apparent. However, the same trend seen in mock galaxies suggests that these correlations are not physical but are rather driven by the LF, EWF, LAE selection, and photometric scatter.

Our base model provides an excellent fit ($\chi_r^2 \approx 1$) to the data, suggesting that any correlation must be weak at best. Nevertheless, we investigate whether a better or an equally good fit can be obtained when we require that galaxies avoid a certain zone in the parameter space. Introducing a real correlation between Lyα luminosity and UV luminosity would require us to drastically change our simulation setup, so we instead opted to obtain it artificially through our random extractions. Specifically, we assume that, at a fixed $m_{1700}$, the maximum or minimum $L_{\text{Ly}\alpha}$ allowed for a galaxy is given by $L_{\text{Ly}\alpha} = p \times m_{1700} + q$. For a given set of $(p, q)$, two models can exist in which the relevant properties of galaxies can lie above or below the line, which we will refer to as $(p, q, 0)$ and $(p, q, 1)$, respectively. We vary the parameter range in $p = [-5, 5]$ and we derive $q$ by forcing the lines to pass through the same point, making this essentially a one-parameter model. The Ando effect can be emulated by a family of models with $(p, q, 1)$.

We repeat the same extraction procedure but this time choose only sources that satisfy the new restriction, and create the expected galaxy distribution. In Figure 9, we show the best-fit model, which is $p = -0.24_{-0.14}^{+0.07}$ and $q = 48.67_{-3.24}^{+1.62}$ Å for our $z \sim 3.8$ sample,<fnref n="6" /> yielding a $\chi_r^2$ of 1.57, considerably worse than the best fit with no restriction. This best-fit model is in the configuration $(p, q, 0)$. In the same figure, we also show the chi-square distribution of all tested models (top panel). This distribution shows how there is a large number of solutions that yield a poor fit, extending to high chi-square values.

Using the goodness-of-fit values obtained from the series of restrictive models, we identify the parameter space with the most constraining power. In practice, we define it as the region that most often yields high reduced $\chi^2$ values, i.e., greater than $3\sigma$. We identify this region through the fraction of solutions ($f(\chi^2)$) that yield a reduced $\chi^2$ value greater than $3\sigma = 9$. Figure 10 shows the UV continuum-EW plane and the UV continuum-Lyα luminosity plane color-coded according to $f(\chi^2 > 9)$. It can be seen how the regions of these planes corresponding to bright LAEs are systematically disfavoured. Gathering new data that systematically better sample the high-luminosity end of the UV luminosity function and the faint end of the Lyα luminosity function should provide better insight

<footnote n="6">⁶ The errors on the fit parameter are computed as the minimum and maximum values that *p* takes within a $1\sigma$ range of $\chi_r^2$ values centered on the best-fit *p* value. The errors on *q* are computed by taking the *q* value corresponding to the maximum and minimum values of *p* because this is essentially a one-parameter fit. Notice that an increase in *p* corresponds to a decrease in *q* and vice versa.</footnote>
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Best-fit solution to the distribution of real and simulated LAEs on the planes defined by observable quantities. Top: the distribution of reduced $\chi_r^2$ values. Bottom: UV continuum luminosity-EW plane and UV continuum luminosity-Ly$\alpha$ luminosity plane. Cyan points represent real LAEs from our PCF ($z \sim 3.790$) sample; orange points represent random extractions from our simulated LAE sample. Black contours encompass the given percentages of the simulated LAE distribution. The introduction of any intrinsic correlation yields a poorer fit to the observations.</figcaption></figure>

into the possible presence of a correlation between Ly$\alpha$ luminosity and UV continuum luminosity.

# 7. Discussion

## 7.1. Implications for Galaxy Formation in a Protocluster Environment

In this work, we have examined how the Ly$\alpha$ properties vary as a function of large-scale environment, finding no significant difference in the shape of the Ly$\alpha$ LF and EWF between overdense regions and the general protocluster environment (aside from an obvious increase in the LF normalization in the denser environments). Though limited by small-number

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Regions of the UV continuum-EW plane (bottom) and of the UV continuum-Ly$\alpha$ luminosity plane (top) color-coded according to the fraction of solutions that have a reduced $\chi_r^2 > 9$.</figcaption></figure>

statistics, doing so allows us to contemplate the future prospect of Ly$\alpha$ mapping as a way to study galaxies in a dense protocluster environment, or alternatively, to explore how to use Ly$\alpha$-emitting galaxies to search for massive protoclusters. Such insights will serve well in the era of the Vera C. Rubin Observatory Legacy Survey of Space and Time and HETDEX (Hill et al. 2008) to develop realistic expectations.

The present study strongly suggests that the line properties (both LF and EWF) measured for the protocluster LAEs are broadly consistent with that determined in the average field, in good agreement with our previous studies based on multiple protoclusters (Dey et al. 2016; Shi et al. 2019). This is in mild disagreement, however, with other works (see, e.g., Lemaux et al. 2018, and references therein), which found a reduced Ly$\alpha$ emission and smaller EWs for galaxies in the protocluster environment (although, in the case of Lemaux et al. 2018 this could be due to their galaxies being more luminous and evolved sources than those investigated here). Xue et al. (2017) showed that the similarity in the line properties between protocluster and field LAEs extends to the circumgalactic scale (but see Matsuda et al. 2012).

To date, only a handful of high-redshift protoclusters have received deep wide-field imaging with a narrowband filter sampling their redshifted Ly$\alpha$ emission (Matsuda et al. 2005; Lee et al. 2014; Dey et al. 2016; Bädescu et al. 2017; Xue et al. 2017; Ouchi et al. 2018; Harikane et al. 2019; Higuchi et al. 2019). Despite their small number, these case studies have demonstrated that LAEs can effectively map the large-scale structure in and around these protoclusters, which show spatial morphologies reminiscent of dark matter structures (e.g., Boylan-Kolchin et al. 2009).
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Redshift evolution of the best-fit parameters for the Lyα LF and EW distributions. From top to bottom the panels show the redshift evolution of the parameters (Φ, L*) of the Schechter (1976) function and the exponential characteristic EW (ω). In each panel, black circles are values from the literature, sorted by redshift, and the red square corresponds to our high-redshift LAE sample. Values from the literature come from Konno et al. (2016, z = 2.2), Ouchi et al. (2008, z = 3.1, 3.7, 5.7), Gronwall et al. (2007, z = 3.095), Ciardullo et al. (2012, z = 3.113), Guaita et al. (2010, z = 2.063), Blanc et al. (2011, z = 2.85, 2.35, 3.3), Sobral et al. (2017, z = 2.23), van Breukelen et al. (2005, z = 3.45), Shimasaku et al. (2006, z = 5.7), Kashikawa et al. (2011, z = 5.7, 6.5), Sobral et al. (2018, z = 2.2, 2.5, 3.1, 3.9, 4.7, 5.4), and Guo et al. (2018, z = 3.1). These literature data are corrected for different cosmology whenever applicable. In the middle panel, the orange line shows our best fit to the log L*–z relation. In the middle and bottom panels, the green diamond and blue star show the best-fit parameters obtained for our “HD” and “Structure” samples of high-density LAEs, respectively.</figcaption></figure>

Taken together, these recent results appear to reinforce the notion that a deep and wide-field LAE survey will be effective in finding galaxies above a given threshold of SFR with little to no bias in regard to the large-scale environment. In the era of wide-field imaging surveys such as the Vera C. Rubin Observatory Legacy Survey of Space and Time, a sensitive wide-field (≥50 deg²) LAE survey of suitable depth will provide a competitive and relatively economical avenue to discover a large sample of protoclusters.

## 7.2. Redshift Evolution of the Lyα Emitters

The results we obtained at z = 3.78 are consistent with recent studies in the literature performed at several different redshifts. Sobral et al. (2018) compiled large samples of LAEs at z ~ 2–6 and reported a continuous rise of L* with redshift by a factor of ~5 (see also Konno et al. 2016). During this period, the φ* parameter decreases with redshift by a factor of 7, most of which occurs at z ≳ 3.

Figure 11 shows the LF and EWF best-fit parameter values from the literature sorted by increasing redshift and also includes the measurements obtained in this work. This figure shows that our estimate is within the general trend with redshift as described by the literature and that it is consistent with measurements of the LF and EWF performed at similar redshifts to our LAE sample. In particular, the L* parameter slowly decreases with increasing redshift. Our measurement is in agreement (within the errors) with values at similar redshift. We tested the possibility of fitting a linear evolution of L* with

redshift and obtained a best-fit solution of

$$ \log L^* = -(0.17 \pm 0.02)z + (43.6 \pm 0.1). \quad (15) $$

The fact that the characteristic luminosity L* has had a shallow evolution over a broad redshift range may be an indication that the properties of the LAE population have not changed dramatically across this period of cosmic time.

The LF normalization and the characteristic EW seem to have a less clear trend (although this could be due to cosmic variance or, in the case of characteristic EW, to the different parameterizations by different works). In particular, the normalization of the LF may be a more difficult parameter to constrain, which presents more variance among the samples (e.g., several of the works mentioned in Figure 11 rely on samples drawn from fields where structures may be present).

The same figure also reports our measurements of L* and characteristic EW for the LAEs identified in different environments, showing no detectable difference from the general case. Based on this analysis it can be concluded that from the viewpoint of Lyα properties there seems to be no difference between LAEs in the most overdense regions and those in the general protocluster environment (at least at a level that can be detected with our current sample size).

## 7.3. Lyα Escape Fraction

The current description of the physical properties of LAEs is that of dusty, low-mass, star-forming galaxies. Star-forming regions in these galaxies are responsible for the production of both the Lyα radiation and the UV continuum emission. Of these two quantities, Lyα emission is the most affected by dust
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Redshift evolution of the ratio of the SFRD derived from Lyα radiation and that from the UV. The red square corresponds to our high-redshift LAE sample. Black points are derived from Table 6 of Sobral et al. (2018) divided by the dust-corrected values from Tables 7 and 5 of Bouwens et al. (2015) and Reddy & Steidel (2009), respectively.</figcaption></figure>

extinction, while UV continuum emission can be considered as a better tracer of the total star formation activity, both obscured and unobscured. On the other hand, any measurement of the SFR derived from Lyα radiation will trace only the unobscured star formation activity. For this reason, the ratio of the Lyα luminosity density to the UV luminosity density for a population of LAEs is a proxy for the fraction of unobscured star formation characterizing a galaxy population and can yield information on its typical dust properties.

For this reason, we have derived the ratio of the Lyα and UV luminosity densities (converted to SFRD) using the Lyα LF derived in our work and UV LFs from the literature. Using the best-fit parameters for the Lyα LF derived in this work, we have integrated the Schechter (1976) function down to a luminosity of $1.75 \times 10^{41} \text{ erg s}^{-1}$ (i.e., down to $0.04L_{z=3}^*$, adopting the value from Gronwall et al. 2007) to be consistent with what was done by Sobral et al. (2018).

We use the best-fit parameters for the UV LF from Table 6 of Bouwens et al. (2015, at $z \approx 3.8$ consistent with our $z \approx 3.790$ sample) together with the mean dust extinction correction factor of 2.4 quoted in the text for the relevant redshift. We integrate the UV luminosity functions down to $0.04L_{z=3}^*$ as derived by Steidel et al. (1999).

Following Equation (11) of Sobral et al. (2018), we derive the ratio of the Lyα-to-UV SFRD as

$$ \frac{\text{SFRD}_{\text{Ly}\alpha}}{\text{SFRD}_{\text{UV}}} = \xi_{\text{ion,N}} \times f_{\text{esc}} \qquad (16) $$

where $\xi_{ion,N} = \xi_{ion}/(1.3 \times 10^{25} \text{ Hz erg}^{-1})$. Following Equation (9) of Sobral et al. (2018),

$$ \xi_{\text{ion}} \times f_{\text{esc}} = \frac{\rho_{\text{Ly}\alpha}}{8.7 c_{\text{H}\alpha} \rho_{\text{UV}}} \qquad (17) $$

where $\xi_{ion}$ is the production efficiency of hydrogen-ionizing photons and $c_{H\alpha} = 1.37 \times 10^{12} \text{ erg}$ is the recombination coefficient.

Figure 12 shows the evolution with redshift of SFRD$_{Lyα}$/SFRD$_{UV} \propto \xi_{ion} \times f_{esc}$ for our work, together with values from Sobral et al. (2018). The overall trend is for the SFRD ratio to slowly increase with redshift from below 5% to

above ~10% in the redshift range $z \sim 2-4$, reaching 35% at redshift $z \sim 5.5$. Our estimated ratio of SFRD$_{Lyα}$/SFRD$_{UV}$ is consistent with other measurements from the literature at a similar redshift, with a value of ~10%. We fit the redshift evolution of SFRD$_{Lyα}$/SFRD$_{UV}$ with an exponential function of $1+z$ and obtain a best-fit result in the form of

$$ \frac{\text{SFRD}_{\text{Ly}\alpha}}{\text{SFRD}_{\text{UV}}} = 8.3 \times 10^{-4} \times (1 + z)^{3.134} \qquad (18) $$

An increase in the SFRD with redshift corresponds to both a change in the amount of ionizing radiation available and a change in the ISM properties. In particular, Sobral et al. (2018) argue that modeling the change in $\xi_{ion}$ as a linear evolution with $1+z$ results in substantial evolution of $f_{esc}$ with redshift. The general trend seems to suggest that with cosmic time a larger fraction of SFR activity happens in dust-obscured environments. At $z \approx 3.8$, the majority of stars are formed in a dust-rich environment, allowing for a very small amount of Lyα radiation to escape from these galaxies. Several high-redshift protoclusters show evidence of massive dusty galaxies; e.g., dusty submillimeter galaxies are detected in a $z=4.3$ protocluster (Miller et al. 2018; Hill et al. 2020; Rotermund et al. 2021). The detection of an increase in the escape fraction of LAEs with redshift suggests the scenario in which the amount of dust in these lower-mass galaxies increases with cosmic time.

## 8. Conclusions and Summary

In this work, we present the Lyα luminosity function and equivalent width distribution of a sample of LAEs at $z \approx 3.8$ in the field of a dense protocluster. The field contains one of the largest and most massive protoclusters known to date with more than 60 spectroscopically confirmed members, providing a rare opportunity to examine how these properties change with the large-scale environment of galaxies. Based on the analysis of this sample, we find:

1. In terms of normalization, the LF of our entire LAE sample is a factor of ~2.0-2.5 higher than that in the average field (e.g., Ouchi et al. 2008; Cassata et al. 2011) while the overall shape remains similar, reflecting the significant galaxy overdensity in the protocluster region.

2. There is no significant difference in the Lyα luminosity function and equivalent width distributions of the LAE population inhabiting the most overdense regions with respect to those of the LAE population inhabiting the general protocluster environment in the PCF field (other than the overall normalization). From the Lyα viewpoint, the implication is that LAEs in the dense protocluster environment form in a manner indistinguishable from that in the average field. However, our conclusion may be tempered by the reduced size of our sample.

3. The measured Lyα luminosity function at $z \approx 3.8$ is broadly in agreement with the existing measures in the literature. The characteristic luminosity $L^*$ shows a mild evolution, decreasing with increasing redshift. We find that the EW distribution does not evolve with redshift, in agreement with the existing measures.

4. By measuring the ratio of the integrated Lyα LF to the integrated UV LF from the literature we are able to estimate the ratio of the unobscured-to-total star
formation rate density. This is a proxy for the escape fraction of our LAE sample. Our measurement of $\xi_{\text{ion}} \times f_{\text{esc}} \sim 10\%$ is consistent with values from the literature. We find the escape fraction to increase with redshift, with a very clearly defined trend.

5. Using a large suite of image simulations, we examine the possibility that an intrinsic correlation exists between a galaxy's UV luminosity, line luminosity, and EW. Such correlations would introduce a bias in the measurement of the Ly$\alpha$ luminosity function and EWF parameters. We find that the apparent correlation in our measurements can be explained by photometric selection effects. Introduction of any intrinsic correlation yields a poorer fit to the observations.

This research has been supported by the funding for the ByoPiC project from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program grant agreement ERC-2015-AdG 695561. Based in part on observations at Kitt Peak National Observatory at NSF's NOIRLab (NOIRLab Prop. ID 2012A-0454, 2014A-0164; PI: Kyoung-Soo Lee), which is managed by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation. The authors are honored to be permitted to conduct astronomical research on Iolkam Du'ag (Kitt Peak), a mountain with particular significance to the Tohono O'odham. This work also made use of images and data products provided by the NOAO Deep Wide-Field Survey (Jannuzi and Dey 1999), which was supported by the National Optical Astronomy Observatory (NOAO, now NOIRLab). The work of A. Dey is supported by NOIRLab. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.

## ORCID iDs

Nicola Malavasi https://orcid.org/0000-0001-9033-7958

Kyoung-Soo Lee https://orcid.org/0000-0003-3004-9596

Arjun Dey https://orcid.org/0000-0002-4928-4003

Rui Xue https://orcid.org/0000-0001-7689-9305

Yun Huang https://orcid.org/0000-0002-2073-5325

Ke Shi https://orcid.org/0000-0003-4747-9656

## References

Adams, J. J., Blanc, G. A., Hill, G. J., et al. 2011, ApJS, 192, 5

Adams, J. J., Hill, G. J., & MacQueen, P. J. 2009, ApJ, 694, 314

Ajiki, M., Mobasher, B., Taniguchi, Y., et al. 2006, ApJ, 638, 596

Ando, M., Ohta, K., Iwata, I., et al. 2006, ApJL, 645, L9

Blanc, G. A., Adams, J. J., Gebhardt, K., et al. 2011, ApJ, 736, 31

Bouwens, R. J., Illingworth, G. D., Blakeslee, J. P., Broadhurst, T. J., & Franx, M. 2004, ApJL, 611, L1

Bouwens, R. J., Illingworth, G. D., Franx, M., & Ford, H. 2007, ApJ, 670, 928

Bouwens, R. J., Illingworth, G. D., Oesch, P. A., et al. 2015, ApJ, 803, 34

Boylan-Kolchin, M., Springel, V., White, S. D. M., Jenkins, A., & Lemson, G. 2009, MNRAS, 398, 1150

Bruzual, G., & Charlot, S. 2003, MNRAS, 344, 1000

Bädescu, T., Yang, Y., Bertoldi, F., et al. 2017, ApJ, 845, 172

Cassata, P., Le Fèvre, O., Garilli, B., et al. 2011, A&A, 525, A143

Charlot, S., & Fall, S. M. 1993, ApJ, 415, 580

Ciardullo, R., Gronwall, C., Wolf, C., et al. 2012, ApJ, 744, 110

Cowie, L. L., & Hu, E. M. 1998, AJ, 115, 1319

Dey, A., Lee, K.-S., Reddy, N., et al. 2016, ApJ, 823, 11

Dijkstra, M., Haiman, Z., & Spaans, M. 2006, ApJ, 649, 37

Dijkstra, M., Lidz, A., & Wyithe, J. S. B. 2007, MNRAS, 377, 1175

Faber, S. M., Phillips, A. C., Kibrick, R. I., et al. 2003, Proc. SPIE, 4841, 1657

Ferguson, H. C., Dickinson, M., Giavalisco, M., et al. 2004, ApJL, 600, L107

Finkelstein, S. L., Rhoads, J. E., Malhotra, S., & Grogin, N. 2009, ApJ, 691, 465

Fossati, M., Fumagalli, M., Lofthouse, E. K., et al. 2021, MNRAS, 503, 3044

Gawiser, E., Francke, H., Lai, K., et al. 2007, ApJ, 671, 278

Gawiser, E., van Dokkum, P. G., Gronwall, C., et al. 2006a, ApJL, 642, L13

Gawiser, E., van Dokkum, P. G., Herrera, D., et al. 2006b, ApJS, 162, 1

Gronwall, C., Ciardullo, R., Hickey, T., et al. 2007, ApJ, 667, 79

Guaita, L., Gawiser, E., Padilla, N., et al. 2010, ApJ, 714, 255

Guo, Y., Jiang, L., Egami, E., et al. 2020, ApJ, 902, 137

Hansen, M., & Oh, S. P. 2006, MNRAS, 367, 979

Harikane, Y., Ouchi, M., Ono, Y., et al. 2019, ApJ, 883, 142

Hashimoto, T., Ouchi, M., Shimasaku, K., et al. 2013, ApJ, 765, 70

Hayes, M., Östlin, G., Schaerer, D., et al. 2010, Natur, 464, 562

Herenz, E. C., Wisotzki, L., Saust, R., et al. 2019, A&A, 621, A107

Higuchi, R., Ouchi, M., Ono, Y., et al. 2019, ApJ, 879, 28

Hill, G. J., Gebhardt, K., Komatsu, E., et al. 2008, in ASP Conf. Ser. 399, Panoramic Views of Galaxy Formation and Evolution, ed. T. Kodama, T. Yamada, & K. Aoki (San Francisco, CA: ASP), 115

Hill, R., Chapman, S., Scott, D., et al. 2020, MNRAS, 495, 3124

Hu, E. M., Cowie, L. L., McMahon, R. G., et al. 2002, ApJL, 568, L75

Hu, W., Wang, J., Infante, L., et al. 2021, NatAs, 5, 485

Hu, W., Wang, J., Zheng, Z.-Y., et al. 2019, ApJ, 886, 90

Inoue, A. K., Shimizu, I., Iwata, I., & Tanaka, M. 2014, MNRAS, 442, 1805

Ito, K., Kashikawa, N., Toshikawa, J., et al. 2020, ApJ, 899, 5

Iye, M., Ota, K., Kashikawa, N., et al. 2006, Natur, 443, 186

Jannuzi, B. T., & Dey, A. 1999, in ASP Conf. Ser. 191, Photometric Redshifts and the Detection of High Redshift Galaxies, ed. R. Weymann et al. (San Francisco, CA: ASP), 111

Kashikawa, N., Shimasaku, K., Matsuda, Y., et al. 2011, ApJ, 734, 119

Kobayashi, M. A. R., Totani, T., & Nagashima, M. 2010, ApJ, 708, 1119

Kodaira, K., Taniguchi, Y., Kashikawa, N., et al. 2003, PASJ, 55, L17

Konno, A., Ouchi, M., Nakajima, K., et al. 2016, ApJ, 823, 20

Kudritzki, R.-P., Méndez, R. H., Feldmeier, J. J., et al. 2000, ApJ, 536, 19

Lai, K., Huang, J.-S., Fazio, G., et al. 2008, ApJ, 674, 70

Laursen, P. 2010, PhD Thesis, Niels Bohr Institute, Univ. Copenhagen

Law, D. R., Steidel, C. C., Shapley, A. E., et al. 2012, ApJ, 759, 29

Lee, K.-S., Dey, A., Hong, S., et al. 2014, ApJ, 796, 126

Lemaux, B. C., Le Fèvre, O., Cucciati, O., et al. 2018, A&A, 615, A77

Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415

Malhotra, S., & Rhoads, J. E. 2002, ApJL, 565, L71

Marinoni, C., Davis, M., Newman, J. A., & Coil, A. L. 2002, ApJ, 580, 122

Matsuda, Y., Yamada, T., Hayashino, T., et al. 2005, ApJL, 634, L125

Matsuda, Y., Yamada, T., Hayashino, T., et al. 2012, MNRAS, 425, 878

Miller, T. B., Chapman, S. C., Aravena, M., et al. 2018, Natur, 556, 469

Nilsson, K. K., Møller, P., Møller, O., et al. 2007, A&A, 471, 71

Nilsson, K. K., Møller-Nilsson, O., Møller, P., Fynbo, J. P. U., & Shapley, A. E. 2009a, MNRAS, 400, 232

Nilsson, K. K., Tapken, C., Møller, P., et al. 2009b, A&A, 498, 13

Ning, Y., Jiang, L., Zheng, Z.-Y., et al. 2020, ApJ, 903, 4

Ono, Y., Ouchi, M., Shimasaku, K., et al. 2010, ApJ, 724, 1524

Ouchi, M., Harikane, Y., Shibuya, T., et al. 2018, PASJ, 70, S13

Ouchi, M., Shimasaku, K., Akiyama, M., et al. 2008, ApJS, 176, 301

Ouchi, M., Shimasaku, K., Furusawa, H., et al. 2003, ApJ, 582, 60

Oyarzún, G. A., Blanc, G. A., González, V., Mateo, M., & Bailey, J. I., III 2017, ApJ, 843, 133

Partridge, R. B., & Peebles, P. J. E. 1967, ApJ, 147, 868

Pirzkal, N., Malhotra, S., Rhoads, J. E., & Xu, C. 2007, ApJ, 667, 49

Rauch, M., Haehnelt, M., Bunker, A., et al. 2008, ApJ, 681, 856

Reddy, N. A., & Steidel, C. C. 2009, ApJ, 692, 778

Rhoads, J. E., Malhotra, S., Dey, A., et al. 2000, ApJL, 545, L85

Ribeiro, B., Le Fèvre, O., Tasca, L. A. M., et al. 2016, A&A, 593, A22

Rotermund, K. M., Chapman, S. C., Phadke, K. A., et al. 2021, MNRAS, 502, 1797

Salpeter, E. E. 1955, ApJ, 121, 161

Santos, M. R., Ellis, R. S., Kneib, J.-P., Richard, J., & Kuijken, K. 2004, ApJ, 606, 683

Schechter, P. 1976, ApJ, 203, 297

Shapley, A. E., Steidel, C. C., Pettini, M., & Adelberger, K. L. 2003, ApJ, 588, 65

Shi, K., Huang, Y., Lee, K.-S., et al. 2019, ApJ, 879, 9

Shibuya, T., Ouchi, M., & Harikane, Y. 2015, ApJS, 219, 15

Shimakawa, R., Kodama, T., Hayashi, M., et al. 2018, MNRAS, 473, 1977

Shimasaku, K., Kashikawa, N., Doi, M., et al. 2006, PASJ, 58, 313

Shioya, Y., Taniguchi, Y., Sasaki, S. S., et al. 2009, ApJ, 696, 546
Sobral, D., Matthee, J., Best, P., et al. 2017, MNRAS, 466, 1242

Sobral, D., Santos, S., Matthee, J., et al. 2018, MNRAS, 476, 4725

Steidel, C. C., Adelberger, K. L., Giavalisco, M., Dickinson, M., & Pettini, M. 1999, ApJ, 519, 1

van Breukelen, C., Jarvis, M. J., & Venemans, B. P. 2005, MNRAS, 359, 895

van de Weygaert, R. 1994, A&A, 283, 361

Vanzella, E., Giavalisco, M., Dickinson, M., et al. 2009, ApJ, 695, 1163

Verhamme, A., Schaerer, D., Atek, H., & Tapken, C. 2008, A&A, 491, 89

Verhamme, A., Schaerer, D., & Maselli, A. 2006, A&A, 460, 397

Willis, J., Courbin, F., Kneib, J.-P., & Minniti, D. 2006, NewAR, 50, 70

Xue, R., Lee, K.-S., Dey, A., et al. 2017, ApJ, 837, 172

Zaninetti, L. 1990, A&A, 233, 293

Zheng, Z., Cen, R., Trac, H., & Miralda- Escudé, J. 2010, ApJ, 716, 574

Zheng, Z.-Y., Wang, J., Rhoads, J., et al. 2017, ApJL, 842, L22