The Role of Dust, UV Luminosity and Large-scale Environment on the Escape of Lyα Photons: A Case Study of a Protocluster Field at z = 3.1

Among the 93 LAEs, S19 reported that 24 sources reside in a significant LAE overdensity near a spectroscopically confirmed protocluster at z = 3.13 (Toshikawa et al. 2016). The observed overdensity parameter is comparable to that measured for several confirmed protoclusters (e.g., Kurk et al. 2000; Venemans et al. 2007; Lee et al. 2014b; Dey et al. 2016a). Moreover, a luminous and extended Lyα nebula (L_Lyα ≈ 2 × 10⁴³ erg s⁻¹, with a half-light radius of 40 kpc) was also discovered near the overdensity, as is the case for several known protoclusters (Steidel et al. 2000; Matsuda et al. 2004; Bădescu et al. 2017). These findings lend credence to the possibility that these LAEs trace a protocluster structure at z = 3.13 that will evolve into a galaxy cluster with a total mass of (1.0–1.5) × 10¹⁵ M_⊙ by z = 0.

We use the the protocluster LAE sample defined by S19. Briefly, S19 constructed the LAE surface density map by smoothing the LAE positions by a Gaussian kernel with the FWHM of 10 Mpc (51). The highest overdensity region outlined by an isodensity contour contains 21 LAEs within an area 72.8 arcmin² corresponding to a surface overdensity δ_Σ = 3.3 ± 0.9 (${\delta }_{{\rm{\Sigma }}}\equiv ({\rm{\Sigma }}-\bar{{\rm{\Sigma }}})/\bar{{\rm{\Sigma }}}$, where $\bar{{\rm{\Sigma }}}$ and Σ are mean and local surface overdensities). Additionally, we include three sources that reside in an LBG overdensity and were later spectroscopically confirmed by Toshikawa et al. (2016) to lie at z = 3.13; they are also recovered by our LAE selection. Thus, our protocluster LAE sample consists of 24 members; the remainder are referred to as "field LAEs." While these definitions are simplistic, they provide us with a rare opportunity to conduct a quantitative and fair comparison of the two environmental subsamples. For more details of the protocluster, its galaxy members, and their distribution, we refer interested readers to S19 (their Section 4.1 and Figure 6).

We make use of the VANDELS public spectroscopic survey DR2 data (McLure et al. 2018; Pentericci et al. 2018). The majority of VANDELS targets are drawn from three categories: bright star-forming galaxies (i_AB ≤ 25) at z = 2.4–5.5; massive quiescent galaxies at z = 1.0–2.5 with H_AB ≤ 22.5; and faint star-forming galaxies at a higher redshift range (25 ≤ H_AB ≤27, i_AB ≤ 27.5 at z = 3–7). To the combined galaxy sample, we apply further selection criteria on the VANDELS galaxies to construct an LAE-like sample, which include: (1) the redshift quality flag value is 3, 4, or 9 yielding ≳80% confidence in the redshift determination; (2) Lyα line in emission is visible in the 1D spectrum. Our selection results in 18 galaxies at z = 3.0–4.9.

Starting from the publicly available 1D spectra, we compute Lyα luminosities by integrating the flux density over a wavelength window of 50 Å, centered at the redshifted Lyα emission (λ_Lyα = 1215.67(1 + z) Å). The continuum flux density in the window is determined by the interpolation from the flux density blueward and redward to the window, which is then subtracted to obtain the line flux. ⁷

The power-law slope of the UV spectrum β (where f_λ ∝ λ^β ) and UV continuum luminosity at rest-frame wavelength 1700 Å (L₁₇₀₀) are determined from the VANDELS multiwavelength photometric catalog. For each galaxy, β is computed through the linear fitting of the flux densities of the bands that sample the UV portion of the spectral energy distribution (SED) but not Lyα: i.e., the izH bands at z > 3.2 and the iz bands at z ≤ 3.2. Once β is determined, L₁₇₀₀ is estimated by interpolating the i-band flux density. The rest-frame Lyα equivalent width (W₀) is derived from the ratio of Lyα flux to the continuum flux density at the same wavelength. The rest-frame W₀ ranges from 0.4 to 440 Å, with a median of 16 Å.

We measure the physical size of the Lyα and continuum emission from the stacked images. Because the r band only samples the continuum emission, the 1D surface brightness profile of the r band is approximated as an exponentially declining function:

$$\begin{eqnarray}&&{I}_{\mathrm{obs}}(r)=\mathrm{PSF}\,\ast \,{I}_{{\rm{c}}}\exp (-r/{r}_{{\rm{s}},{\rm{c}}})\end{eqnarray} \tag{ 2 }$$

where the r_s,c is the exponential scale length of the galaxy's continuum emission, I_c is the normalized surface brightness, and the symbol ∗ denotes the convolution. We use the Python function scipy.optimize.curve_fit and obtain the best-fit scale length 0.3 ± 0.3 kpc: this is consistent with our expectation that our LAEs are unresolved in the r band. The error accounts for the full variance of the parameter space.

The intrinsic Lyα surface brightness is modeled as a superposition of two exponential functions, representing the emission from the galaxy and the extended halo:

$$\begin{eqnarray}&&{S}_{\mathrm{Ly}\alpha ,\mathrm{intrinsic}}={S}_{{\rm{c}}}\exp (-r/{r}_{{\rm{s}},{\rm{c}}})+{S}_{{\rm{h}}}\exp (-r/{r}_{{\rm{s}},{\rm{h}}}).\end{eqnarray} \tag{ 3 }$$

The observed Lyα surface brightness is expressed as:

$$\begin{eqnarray}&&{S}_{\mathrm{Ly}\alpha ,\mathrm{obs}}(r)=\mathrm{PSF}\,\ast \,{S}_{\mathrm{Ly}\alpha ,\mathrm{intrinsic}}.\end{eqnarray} \tag{ 4 }$$

In Equation (3), r_s,h is the scale length of the halo component; S_c and S_h normalize the galaxy and halo component, respectively. We determine the best-fit scale length r_s,h while fixing the r_s,c value measured in the r band and obtain r_s,h = 4.9 ± 0.7 kpc.

In the left panel of Figure 2, we show the azimuthally averaged radial profiles of the Lyα (blue) and r band (red) together with the best-fit exponential models as solid lines. For ease of comparison, the profiles are normalized at the innermost angular bin. For the Lyα band, we also show both galaxy and LAH components as a dashed and a dotted line, respectively. The image PSF is indicated by a thick gray line. The uncertainties include both Poisson and bootstrap sampling errors.

Zoom In Zoom Out Reset image size

In order to explore how the LAH size changes with galaxies' large-scale environment, we repeat the same procedure for the two subsamples. In Figure 2 (right), these measurements are shown as pink circles and lines. The uncertainty in the protocluster subsample is ≈95% (≈28% in the field) larger than those of the full sample as the number of galaxies in each sample is smaller. When normalized, these measures are consistent with our result for the full sample within uncertainties. The LAH scale lengths are r_s,h = 5.1 ± 0.9 kpc and r_s,h = 4.7 ± 0.9 kpc for the protocluster and field subsample, respectively. Our analysis suggests that the large-scale environment is not a major determinant of LAH sizes, in agreement with the finding of X17.

Having separated the total Lyα emission into the halo and galaxy components, we estimate the fractional contribution of the LAH to the total line flux. The halo fraction is defined as:

$$\begin{eqnarray}&&{X}_{\mathrm{Ly}\alpha }\equiv {F}_{{\rm{h}}}/{F}_{\mathrm{tot}}\end{eqnarray} \tag{ 5 }$$

where the halo and total fluxes, F_h and F_tot, are computed as:

$$\begin{eqnarray}&&{F}_{{\rm{h}}}={\int }_{{r}_{\mathrm{halo}}}^{\infty }2\pi {{rS}}_{{\rm{h}}}\exp (-r/{r}_{{\rm{s}},{\rm{h}}}){dr};\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{F}_{\mathrm{tot}}={\int }_{0}^{\infty }2\pi r[{S}_{{\rm{c}}}\exp (-r/{r}_{{\rm{s}},{\rm{c}}})+{S}_{{\rm{h}}}\exp (-r/{r}_{{\rm{s}},{\rm{h}}})]{dr}.\end{eqnarray} \tag{ 7 }$$

We note that, unlike for F_tot, the integration interval for F_h in Equation (6) is set to [r_halo, ∞ ]; i.e., F_h encloses the observed flux outside the central region where r = r_halo serves as the inner "edge" of the LAH, or equivalently, the outer edge of the galaxy. This definition is motivated by the fact that, at a distance comparable to or smaller than the size of the host galaxy, the physical meaning of the LAH (as a component separate from the galaxy emission) becomes ambiguous in all but a mathematical sense. However, our definition of F_h would be identical to that used in Wisotzki et al. (2016) and Leclercq et al. (2017) if r_halo is set to zero.

In the middle panel of Figure 3, we show X_Lyα as a function of r_halo. The 1σ uncertainties (gray shade) are computed by bootstrapping the radial profile measurements. At r_halo = 0, our halo fraction is ${73}_{-12}^{+8} \%$, in good agreement with those estimated by Wisotzki et al. (2016) and Leclercq et al. (2017); (light gray and dark gray histograms, respectively, in the left panel of Figure 3), and also consistent with the expectation from simulations with a realistic treatment of the ISM (Verhamme et al. 2012). When evaluated at r_halo = 05 (3.9 kpc at z = 3.13), which lies safely outside typical galaxy sizes, the halo fraction decreases to ${61}_{-12}^{+9} \%$. Our analysis indicates that, for LAEs, the flux from the Lyα halo easily dominates over that originating from the galaxy itself, even in the most conservative estimate.

Zoom In Zoom Out Reset image size

From the observational viewpoint, it is often useful to know what fraction of the flux falls outside a given aperture. The calculation is similar to that shown in Equation (6) except that we integrate the observed surface brightness instead of the intrinsic one. In the right panel of Figure 3, we show the fractional loss of the Lyα (blue) and continuum r-band (red) flux as a function of aperture radius, r₀. We start by convolving the best-fit radial profiles with the Gaussian kernels with the FWHM of 09 (dotted), 12 (solid), and 15 (dashed–dotted), respectively, to simulate the realistic range of seeing in ground-based observations, which also affect the flux loss.

In a 3'' aperture in diameter (r₀ = 15), the flux loss is 32% ± 6% in the Lyα image, much larger than 19.0% ± 0.4% in the r band. The aperture radius that encloses 90% of the total flux, r₉₀, is 24 ± 03 for the Lyα band, considerably smaller than r₉₀ = 37 ± 02 measured for a sample of galaxies with both Hα and Lyα emission (Matthee et al. 2016). The discrepancy is too large to be explained by the difference in image quality.

The origin of this disagreement is likely intrinsic. Although the galaxies in both samples are both detected with significant Lyα line excess, the Hα-Lyα emitters tend to have larger stellar masses and higher star formation rates (SFRs) than our LAEs. Matthee et al. (2016) made the stellar mass estimates based on SED fitting and reported $\mathrm{log}({M}_{\mathrm{star}}/{M}_{\odot })=8.6\mbox{--}11.1$ with a median of 10.3 (see their Table 1). Although we do not have a direct estimate of stellar mass for our LAE sample, photometrically selected LAEs tend to be low-mass galaxies in the range (10⁸–10⁹) M_⊙ (Gawiser et al. 2006; Guaita et al. 2011; Nilsson et al. 2011; Shi et al. 2019b). As for the SFR, Matthee et al. (2016) reported the range (3–50) M_⊙ yr⁻¹ with a median of 20 M_⊙ yr⁻¹, once again, much larger than ≈5 M_⊙ yr⁻¹ for our sample by using dust-corrected SFRs generated by observed UV luminosity (Shi et al. 2019b, also see Section 4.3). In Section 3.3, we discuss how the LAH sizes change with physical properties.

In Figure 4, we show the compilation of Lyα size measurements in the literature as a function of line and continuum luminosity and W₀. The measurements of the full LAE sample and the protocluster and field subsamples are indicated as a blue star, open pentagon, and open diamond, respectively. Overlaid are similar measurements based on image stacking for the photometrically selected LAEs at z ∼ 2.7 and ∼3.8 (X17, orange and purple, respectively) and at z ∼ 2.2 (Momose et al. 2016, open triangles). We also show individual measurements for z ∼ 3–6 star-forming galaxies from deep MUSE observations from Wisotzki et al. (2016, light gray circles) and Leclercq et al. (2017, dark gray circles). The latter measurements are binned ⁹ for clarity. A green circle in each panel represents the stacked measurement of UV-selected star-forming galaxies (Steidel et al. 2011), with a median equivalent width of 0.9 Å (and the range −27 Å ≤ W₀ ≤ 89 Å), the majority of which would not be classified as LAEs. The remainder of the literature samples shown in Figure 4 employed the classical definition of LAEs, i.e., W₀ ≥ 20 Å with the exception of the z ∼ 2.7 sample of X17, which used W₀ ≥ 50 Å.

Zoom In Zoom Out Reset image size

While limited, our size measurements are in line with the existing data. In particular, the agreement is excellent when compared with the measurements that employed a similar decomposition method that simultaneously fits the galaxy and the LAH component (see Section 3.1), namely, those in Leclercq et al. (2017) and X17. Disagreement with the Momose et al. (2014) and Steidel et al. (2011) samples may be in part due to the difference in the fitting method in that they only considered a single component. For an in-depth discussion of how the size measurement can be affected by the fitting method and image point-spread function, we refer interested readers to Appendix C of X17.

Figure 4 showcases the overall trend that the LAH sizes increase with both Lyα luminosity and UV luminosity while showing little correlation with galaxy environment. The very large LAH size (≈25 kpc; Steidel et al. 2011) measured for UV continuum selected star-forming galaxies lies far above these trends, suggesting that a different scaling relation may apply to non-LAEs. The larger LAH sizes for more UV luminous galaxies imply that the flux loss we estimate in Section 3.2 is only applicable to typical LAEs (which are UV faint, M_UV ≳ −20) and cannot be generalized to all star-forming populations.

Discerning how LAH sizes and Lyα surface brightness profiles change with the properties of host galaxies and their large-scale structure can place strong constraints on the dominant physical mechanism that powers the Lyα halo (see, e.g., Matsuda et al. 2012; Momose et al. 2016; X17; Leclercq et al. 2017). Possible scenarios include resonant scattering of the Lyα photons originating from star formation in the CGM (Laursen & Sommer-Larsen 2007; Dijkstra & Loeb 2009; Verhamme et al. 2012), gravitational cooling radiation (Haiman et al. 2000; Fardal et al. 2001; Dijkstra & Loeb 2009; Lake et al. 2015), and star formation in ultra-faint satellite galaxies (Zheng et al. 2011; Lake et al. 2015). While the limited nature of our measurement (based on a single stack) prevents us from placing a new meaningful constraint on these scenarios, we stress that the mild disagreement between existing measurements, the large scatter observed in the individual measurements, and the apparent dichotomy between LAEs and non-LAEs in LAH sizes seen in Figure 4 highlight the incompleteness of the current observational picture. Larger samples spanning a wider range of parameter space (in luminosities and large-scale environments) will be crucial in establishing clearer trends and in discriminating different physical scenarios.

For our analyses in this and subsequent sections, we only consider 62 LAEs for which we have robust measurements of the UV slope β. To this end, we only use the galaxies with Δβ < 0.9. A majority of the sources that do not meet this criterion are simply too faint in the images and typically have inferred UV luminosities log(L₁₇₀₀) ≲ 27.3 erg s⁻¹ Hz⁻¹. Following the convention, we compute the Lyα escape fraction as:

$$\begin{eqnarray}&&{f}_{\mathrm{esc}}=\displaystyle \frac{\mathrm{SFR}(\mathrm{Ly}\alpha )}{\mathrm{SFR}(\mathrm{UV})\times {10}^{(0.4{k}_{1500}E(B-V))}}\end{eqnarray} \tag{ 8 }$$

where the color excess E(B − V) is converted from the UV slope β by assuming the Calzetti et al. (2000) extinction law. For the UV- and Lyα-based SFR, we adopt the Kennicutt (1998) calibration. The term k₁₅₀₀ denotes the effective dust extinction at rest-frame 1500 Å.

The median (mean) value is 60 (${101}_{-101}^{+107}$)% for our LAE sample. Of the 62 LAEs, nearly one-third (21) have f_esc values greater than unity, nine of which lie within the protocluster (see Section 2.1). Nine LAEs (two in the protocluster) do so at a ≥1.5σ level. In comparison, Blanc et al. (2011) reported 24% for spectroscopically selected LAEs with f_esc greater than unity at z = 1.9–3.8. Both numbers are based on dust-corrected UV continuum to infer SFRs and thus are subject to similar systematics. Excluding the protocluster LAEs brings down the fraction of LAEs with f_esc ≥ 1 to 27%, in better agreement with Blanc et al. (2011). Comparison of field versus protocluster LAEs is given in Section 6.

Larger-than-unity f_esc values in these two samples may in part arise from photometric scatter, particularly in the bands used for estimating β. To quantify how robustly the estimates of β and f_esc can be made for individual galaxies, we create image simulations calibrated to closely reflect the brightness and the colors of the real LAEs. While the details of our simulations are given in the Appendix, our result suggests that the photometric recovery of β values is reasonably good at β ≈ −2 (≈0.2 dex). Photometric scatter tends to result in the recovered β values being more positive than the intrinsic ones, leading to the underestimation—and not overestimation—of f_esc values. Furthermore, many of our LAEs with f_esc > 1 have relatively small Δβ, placing their f_esc values >(2–3)σ outside the nominal f_esc = 1 line (see Figure 5). Based on these considerations, we argue that, while we cannot rule out the contribution of photometric scatter on f_esc > 1 sources, it is unlikely to be the sole driver of the f_esc > 1 sources.

Zoom In Zoom Out Reset image size

There are multiple physical factors that may result in the larger-than-unity f_esc. First, radiative transfer through an ISM with a complex geometry may cause Lyα emission to have a preferential direction (e.g., Verhamme et al. 2012) unlike continuum emission; in this case, the ratio of Lyα to UV emission would have little physical meaning. While possible, this is unlikely to be the dominant cause given the relatively tight correlation of f_esc and other galaxy properties (see Section 5). Second, low-level AGN may contaminate our sample. As discussed in Section 2, few of our LAEs have spectroscopic observations. It is also possible that the underlying assumptions we make about these galaxies are wrong: they have had relatively continuous and prolonged star formation histories and have solar or moderately subsolar metallicities. While Lyα luminosity traces instantaneous star formation (<10 Myr), far-UV continuum luminosity represents the SF activity averaged over the last ≈100 Myr (Kennicutt 1998). Similarly, substantially subsolar metallicities would result in a higher ionizing radiation at a fixed mass. Thus, highly episodic SF activities, extremely young ages, and low metallicities will all lead to a higher Lyα output relative to the UV and can result in seemingly unphysical f_esc values (Charlot & Fall 1993; Malhotra & Rhoads 2002; Kashikawa et al. 2012; Hashimoto et al. 2017).

Alternatively, the dust law we assume to convert measured UV flux density to SFR may not be appropriate. Existing studies (e.g., Reddy et al. 2006; Siana et al. 2008) found that UV-selected star-forming galaxies dominated by younger stellar populations may be better characterized by a Small Magellanic Cloud−like dust law (Gordon et al. 2003). Kusakabe et al. (2015) reached a similar conclusion for a large sample of LAEs. However, an SMC-like dust law would exacerbate the problem at hand. Assuming the SMC law, a given UV slope β would correspond to a smaller color excess E(B − V), leading to an even larger f_esc than previously.

After excluding the galaxies with f_esc ≥ 1, the median (mean) value of the Lyα escape fraction for our sample is 36%(40%) ± 26%. Our estimation is in good agreement with other LAE samples in the literature, ¹⁰ including ${29}_{-29}^{+40} \%$ (Blanc et al. 2011), 38% ± 11% (Oteo et al. 2015) and 37% ± 7% (Sobral et al. 2017). One exception is the Hα-emitting LAE sample (11% ± 11%: Matthee et al. 2016), which shows considerably lower f_esc values than the rest. As discussed in Section 3.2, these Hα-Lyα emitters tend to be more UV luminous and have larger stellar masses than those identified based on Lyα emission alone, likely signaling that f_esc correlates with these physical parameters (see Section 4.3).

So far, we have considered the f_esc values measured within the galaxy-sized apertures as has been done in the literature. If Lyα emission from the LAH originates from the same H ii regions as that from the galaxy, the f_esc estimates would need to be revised to account for the additional contribution from the LAH. Using the radial surface brightness profiles of the stacked Lyα and UV images (Section 3.2, Figure 3), we estimate the total Lyα escape fraction by statistically correcting for the expected flux loss. For the correction, we use the radius of the effective circular aperture, defined as ${r}_{0}=\sqrt{{A}_{\mathrm{iso}}/\pi }$ of each LAE. Radius r₀ ranges from 06 to 20 with a median of 10 (7.8 kpc). The Lyα flux loss ranges 6%–77% with a median of 54%, mainly due to the extended nature of the halo. While smaller, the UV flux loss due to PSF blurring (seeing 12) is also not negligible, ranging 10%–69% with a median of 40%.

We find the median (mean) value of the total Lyα escape fraction, 〈f_esc,tot〉, is 46% (51%) ± 32%, in good agreement with 42% ± 24% reported by Kusakabe et al. (2018), which includes the LAH flux in their estimation. Similar to ours, Kusakabe et al. (2018) also defined their LAEs to be galaxies with W₀ ≥ 20 Å. While our value is higher than those in the literature, the relative increase is modest at ≈29% when the flux loss correction is made consistently for both UV and Lyα fluxes. Additionally, the LAH size for typical LAEs is small enough (4.5 kpc corresponds to 057 at z = 3.13) that galaxy-sized apertures can contain much of the flux.

We explore how f_esc values correlate with the galaxies' UV and Lyα properties by combining our own measurements with those in the literature. Since most of these measures do not include the LAH component, we opt to use our f_esc values before the LAH correction. Given the modest change that it brings, our conclusion should not change substantially.

In Figure 5, we show the correlations of f_esc with W₀, UV slope (β), and line and continuum luminosities. In all cases, we show protocluster LAEs with larger symbols than those in the field. Other measurements include 89 galaxies from the HETDEX Pilot Survey (HPS hereafter; Blanc et al. 2011), 17 LAEs with Hα detection (Matthee et al. 2016), and 18 spectroscopic Lyα detections from VANDELS. The average measurements via a stacking analysis of z ∼ 2.2 LAEs (Sobral et al. 2017) are also shown.

The HPS LAEs span a comparable range of UV slope and W₀ to our sample, but tend to have higher line and UV luminosities; the Hα-Lyα emitters studied by Matthee et al. (2016) are much more UV luminous than the rest. The VANDELS sources are comparable to our LAEs in both luminosities but have considerably bluer UV slopes and smaller W₀.

In order to test how these physical properties affect f_esc, we use two nonparametric ranked correlations, namely, the Kendall τ and Spearman ρ rank correlation coefficients (Spearman 1904; Kendall 1938). Both tests are run for each parameter shown in Figure 5 using two Python scripts scipy.stats.kendalltau and scipy.stats.spearmanr. To check the consistency, we run these tests on our data set and the HPS data set separately, and then on the combined data set (our LAEs+HPS and LAEs+HPS+VANDELS). We do not include the Kusakabe et al. (2018) measures on our tests because they represent the stacked averages. In Table 1, we list the correlation coefficients (τ and ρ) and the probabilities of null hypotheses (p-values). We consider the case for which either τ or ρ is larger than 0.6 as a robust correlation, while ≈0 values in these parameters indicate no correlation.

We find a clear anticorrelation between f_esc and the UV continuum slope β, in agreement with existing studies (e.g., Blanc et al. 2011; Song et al. 2014). This is not surprising considering that we use the β value as a proxy for dust reddening to correct the SFR_UV in our estimation of f_esc. A weaker correlation with observed UV luminosity is notable: the trend is clear when the combined data set is considered, which spans two orders of magnitude in luminosities. In this parameter space, the Hα-Lyα emitters are no longer significant outliers but are part of the correlation at the highest luminosity end. The trend is likely related to the fact that lower-luminosity star-forming galaxies tend to have bluer β values at high redshift (e.g., Bouwens et al. 2012, 2014). We further quantify the luminosity dependence in Section 5.2.

No correlation is found with line luminosity or W₀. This is in contrast to the anticorrelation between L_Lyα and f_esc reported by Sobral et al. (2017), which is based on stacking analyses. While their measurements are not inconsistent with the mean values of the combined data set, the intrinsic scatter is too large to discern any meaningful correlation.

All in all, we conclude that the galaxy's dust reddening and UV luminosity are two primary determinants of the escape fraction of Lyα photons. In the next section, we further quantify these dependencies further to shed light on the physical driver of their escape.

We consider our result in the context of several ISM models in the literature. First, in a medium composed of uniform mixtures of dust and H i gas, the effective attenuation of Lyα photons should be much higher than for continuum photons (e.g., Charlot & Fall 1991), as the former traverse much longer path lengths than the latter, leading to f_esc values much lower than the continuum expectation. Our results suggest the opposite.

Alternatively, in a "multi-phase" medium where dust is mostly confined within H i clumps embedded in an otherwise ionized dust-free medium (Neufeld 1991), Lyα photons spend most of their time traveling in the intercloud medium while occasionally scattering at cloud surfaces. Thus, f_esc primarily depends on the cloud albedo (the likelihood of surface absorption by dust) and the number of Lyα scatterings before their eventual escape. The consequence is that the dependence of f_esc on the total dust optical depth is mild at best (Hansen & Oh 2006; Duval et al. 2014). By contrast, continuum photons travel through clumps and undergo more severe attenuation. The selective attenuation can enhance the Lyα W₀ to, in some cases, above the case B expectation. The degree of W₀ enhancement increases with increasing extinction.

One possible exception would be if star formation mainly occurs deep inside cold H i clouds where dust resides: Lyα photons are decimated by dust absorption before they reach the cloud surface, leaving the role of preferential resonant scattering less consequential (Verhamme et al. 2012). All in all, the strong anticorrelation between f_esc and dust reddening, the lack of LAEs with excessively high W₀, and the larger-than-expected f_esc are at odds with the multi-phase ISM.

Finally, Scarlata et al. (2009) considered a "clumpy dust screen" through which Lyα photons take the paths of the least resistance (i.e., escaping through the holes with the lowest opacity between clumps). The process allows f_esc to increase while keeping W₀ unchanged without invoking preferential resonant scattering (also see Atek et al. 2014). One defining characteristic of this model is that the relative Lyα enhancement scales with dust extinction, resulting in k_Lyα < k₁₂₁₆. In Figure 6, the effect would manifest itself as a shallower slope. Resonant scattering would increase the Lyα optical depth, bringing the scaling relation closer to the uniform dust screen expectation by steepening the slope k_Lyα .

The clumpy dust screen scenario also does not align well with the existing measurements. In all three LAE samples in Table 2, we consistently find k_Lyα ≈ k₁₂₁₆ within uncertainties. There is only one exception: if we only consider sources with f_esc ≤ 1, the formal fit favors a considerably lower k_Lyα value (≈6; column 4 in Table 2). A similar reanalysis of the HPS sources, however, results in a marginal change in the slope, from 10.0 to 8.2, suggesting that the dramatic change we observe in our result is caused by a small sample size. Regardless, the removal of those with f_esc ≥ 1 without a firm physical basis is arbitrary at best. If the underlying assumptions we made in deriving the escape fraction are incorrect (as discussed in Section 4.1), the correction needed to recover the true quantities would have to be uniformly applied to all LAEs. Larger sample sizes would improve the measurements, particularly, by getting a better handle on the f_esc values of highly reddened LAEs and their intrinsic scatter.

In summary, our results show that Lyα photons undergo interstellar dust attenuation in a similar manner to continuum radiation (i.e., k_Lyα ≈ k₁₂₁₆). However, the overall Lyα transmission is roughly twice higher than that expected for continuum photons. The constraints derived from one of the largest compilations of LAEs to date appear to favor some form of multi-phase media conducive to the preferential transfer of Lyα photons. Yet, none of the models we have considered is entirely in line with our findings, hinting that the distribution of interstellar gas and dust in distant galaxies is complex. One possibility may be a hybrid of the two aforementioned models in which H i clouds are rendered partially porous due to strong stellar feedback. This would allow for a substantial fraction of Lyα photons to escape from their birth cloud relatively easily and proceed to undergo selective attenuation via resonant scattering. Detailed calculations will require full radiative transfer calculations based on realistic simulations at subparsec scale resolution.

The range of C_Lyα changes substantially from ≈0.45 for the Hayes et al. (2011) sample to ≈1.4 for the HPS sources to ≈2 for the other LAE samples, while the slope k_Lyα remains unchanged: i.e., at a fixed reddening, f_esc varies up to ≈4 in these samples. The high C_Lyα is real, as the majority of our LAEs lie well above the fiducial Calzetti expectation (Figure 6, top right).

We evaluate the goodness of fit by computing the total chi-square value as ${\chi }^{2}={{\rm{\Sigma }}}_{i}{({f}_{\mathrm{esc},\mathrm{model}}({\epsilon }_{i})-{f}_{\mathrm{esc},i})}^{2}/{\sigma }_{\mathrm{esc},{\rm{i}}}^{2}$, where f_esc,i and σ_esc are the estimates of the Lyα escape fraction and its error, _i is the estimated E(B − V) value for the galaxy i, and f_esc,model is computed using Equation (9). Although the measurement errors are unlikely to obey a Gaussian distribution and therefore a χ² value does not carry the usual statistical significance, it should still provide us with a way to evaluate whether a given model gives a better description of the data over another. When we refit the data while forcing C_Lyα to be unity, we obtain unrealistically shallow slopes k_Lyα with a much poorer agreement with the data, with the total χ² increasing from ∼30 to 60. The results for our single-parameter fits are listed in Table 2.

The differences in C_Lyα appear to be linked to galaxies' UV luminosities. Figure 5 shows that the HPS sources tend to be more UV luminous than our LAEs by nearly an order of magnitude, while the W₀ and β distributions of the two samples largely overlap. As for the Hα-Lyα emitters from Matthee et al. (2016), which have the lowest f_esc values, all but four lie at the bright end (L_UV > 10^28.3 erg s⁻¹ Hz⁻¹; M_UV < −20.7), but they too have a W₀ and L_Lyα range similar to that of our LAEs. Thus, we conclude that f_esc depends not only on the galaxy's dust content but also on its UV luminosity.

Although we are limited to analyzing the Lyα-emitting population only, our conclusions are in broad agreement with several studies of the general star-forming galaxy population: Kim et al. (2020) found an anticorrelation between UV luminosity and Lyα escape fraction in a sample of local Lyman Break Galaxy Analogs. Stark et al. (2010) found that the fraction of Lyα-emitting galaxies, X, increases with decreasing UV luminosity among the UV-selected galaxy samples. ¹² Oyarzún et al. (2017) studied a stellar mass selected sample of galaxies at z = 3.0–4.5 and found that both X and f_esc strongly anticorrelate with stellar mass (see also Oyarzún et al. 2016). Given that UV luminosity and stellar mass broadly track each other through the star formation main sequence (González et al. 2011; Lee et al. 2012; Stark et al. 2013; Song et al. 2016), the latter correlation is consistent with the L_UV–f_esc scaling we find. Similarly, Weiss et al. (2021) found a similar correlation between stellar mass and f_esc of [O iii] emitters at z = 1.9–2.4.

Motivated by the dependence of f_esc on UV luminosity as shown in Figures 5–6, we parameterize f_esc as:

$$\begin{eqnarray}&&{f}_{\mathrm{esc}}={C}_{\mathrm{UV}}\cdot {10}^{-0.4{k}_{\mathrm{Ly}\alpha }E(B-V)}.\end{eqnarray} \tag{ 10 }$$

The equation is identical to Equation (9) except for the difference that the normalization, C_UV, is expressed as ${C}_{\mathrm{UV}}\equiv \mathrm{log}{\left(\tfrac{{L}_{\mathrm{UV}}}{{L}_{0}}\right)}^{\alpha }$ where L_UV is given in units of erg s⁻¹ Hz⁻¹ and α and L₀ are constants. We repeat the fitting procedure using the combined sample of our LAEs and the HPS sources. The best-fit k_Lyα = 7.37 ± 0.98, α = −1.62 ± 0.20, and L₀ = (2.90 ± 1.15) × 10²⁹ erg s⁻¹ Hz⁻¹. The α value is firmly in the negative, affirming the fact that the normalization C indeed decreases with increasing luminosity.

In Figure 7, we illustrate our result. The left panel shows our measurements as observed, while, in the right panel, we "correct" for the luminosity dependence by showing ${f}_{\mathrm{esc},\mathrm{renorm}}\equiv {C}_{\mathrm{UV}}^{-1}$ f_esc. Relative to the power law that best describes each set of data points (dashed lines), the scatter decreases from 0.26 to 0.18 dex. In short, Equation (10) allows us to predict a galaxy's f_esc with a ≈50% accuracy provided that both the UV slope β and L_UV are known.

Zoom In Zoom Out Reset image size

The efficacy of the luminosity-dependent f_esc calibration is better illustrated in Figure 8, where we combine L_UV and E(B − V) into a single parameter, i.e., dust-corrected UV SFR. In the top left panel, we once again show the measurements as observed. A clear separation between our LAEs and the HPS sources is visible. The gray swath marks the correlation reported by Oyarzún et al. (2017) for Lyα-emitting galaxies, which shows an excellent agreement with our own measurements. Although their stellar mass range is $\mathrm{log}[{M}_{\mathrm{star}}/{M}_{\odot }]=7.6\mbox{--}10.6$, the majority lies below 10^9.5 M_⊙ and thus should be well matched to the stellar mass range for LAEs (Gawiser et al. 2006; Hathi et al. 2016; Arrabal Haro et al. 2020; Santos et al. 2020).

Zoom In Zoom Out Reset image size

In the two right panels, we renormalize the SFR as SFR_UV,renorm ≡ C_UV · SFR_UV. Both panels show identical data that are color coded by L_UV (top) and E(B − V); (bottom), respectively. The renormalization moves our LAEs to the right, in better alignment with the HPS sources and the Hα-Lyα emitters. Only the former are used in the fitting, but the latter are rescaled consistently in the figure.

A small fraction (9%) of galaxies lie significantly outside an otherwise tight sequence formed by the rest, which we demarcate with a dotted–dashed line. The outlier group largely consists of high UV luminosity galaxies (≳10²⁹ erg s⁻¹ Hz⁻¹, or M_UV ≲ −22.4) with uncharacteristically blue colors (E(B − V) ≲ 0.2) and EWs near the EW cutoff (≈20 Å). Eight of this group belong to the HPS sample, four are Hα-Lyα emitters, and the additional three are in our sample (one in the protocluster region and the other two in the field). Given that both HPS sources and our LAE selection target single line emission, we cannot rule out that some are low-redshift interlopers such as [O ii] emitters at z = 0.34. In Section 6, we discuss how the exclusion or inclusion of these three LAE candidates affects our conclusion on the environmental dependence on f_esc.

After excluding these outliers, we fit a power law to the data and find that the scatter is reduced from 0.29 dex, for the unaltered data, to 0.21 dex. When we account for the expected effect of photometric scatter using image simulations, the intrinsic scatter in the f_esc–C_UVSFR_UV scaling relation is less than 0.19 dex.

The physical origin of the luminosity-dependent C_Lyα is unclear. However, it is worth noting that our result does not mean that C_UVSFR_UV is the true SFR; in such a case, f_esc (≡SFR_Lyα /SFR_UV) in the abscissa should also scale accordingly, which would undo the alignment.

One possibility is that dust laws change with UV luminosity such that more UV luminous galaxies obey the Calzetti et al. (2000) law, while less luminous ones gradually transition to a law more similar to the SMC law. Kusakabe et al. (2015) reported that the average IR luminosity of LAEs at z = 2.2 measured from the Spitzer/MIPS and Herschel/PACS image stacks lie well below that expected under the Calzetti et al. (2000) extinction, evidence that the LAE population may be better described by the SMC law. Similarly, Reddy et al. (2012) found that UV-selected star-forming galaxies at z ∼ 2 with young stellar population ages (≲100 Myr) appear to obey an SMC-like dust law, while older galaxies follow the Calzetti et al. (2000) law.

In Figure 9, we show the same data but this time we apply the SMC law for our LAEs while the HPS sources and the Hα-Lyα emitters still obey the Calzetti law. The tightness of the correlation is comparable to that in Figure 8 but with fewer outliers. Modeling a luminosity-dependent change in dust laws could conceivably further tighten the scaling relation. The transformation shown in the figure effectively moves our LAEs (in the top left panel of Figure 8) to lower SFRs (to the left along the x-axis) and to higher f_esc values by the same factor.

Zoom In Zoom Out Reset image size

One caveat of this interpretation is that it drives already high f_esc values even higher. Nearly all of our LAEs would have Lyα escape fractions ≫100%. The problem may be mitigated if low UV luminosity (or SMC laws) go hand in hand with low metallicity, young ages, or level of burstiness in star formation histories: these traits either boost the production of ionizing radiation (and thus Lyα production) or enhance Lyα transmission. A comprehensive study of these properties in LAEs is needed to test this hypothesis.

Higher Lyα escape fraction for the protocluster LAEs is intriguing. Recent absorption-line studies found that the regions of galaxy overdensities are also richer in H i than average fields (Lee et al. 2014a, 2018; Cai et al. 2016; Liang et al. 2021). Liang et al. (2021) reported a cross-correlation length 4 ± 1 Mpc (comoving) between the LAE overdensity and H i line-of-sight optical depth. These results are broadly consistent with the expectation that both galaxy and gas overdensities track those of the underlying dark matter. Higher H i column densities at the galaxy scales would result in more frequent random walks of Lyα photons, thereby enhancing the chance of their destruction by dust grains along their paths and thus lowering f_esc values. Our result runs counter to this expectation.

One possible explanation is that different ISM conditions for protocluster galaxies facilitate Lyα photons to escape more easily relative to field galaxies. Gazagnes et al. (2020) studied local star-forming galaxies and found that the escape fraction (for both Lyα and Lyman continuum photons) negatively correlates with H i covering fraction and dust reddening, favoring the scenario in which low column density channels serve as privileged routes through which Lyα and LyC photons can escape the galaxy. Such photoionized "tunnels" can be carved out by supernova feedback (e.g., Kimm & Cen 2014) or by the turbulent early phase of the H ii regions and the surrounding molecular clouds (Kakiichi & Gronke 2019). In protocluster galaxies embedded in the H i richer large-scale environment, faster growth brought on by a higher star formation efficiency (e.g., Chiang et al. 2017) may work through these processes to effectively create a more porous ISM relative to that of their field counterparts. Such a scenario would have a strong implication for the role of protocluster galaxies in cosmic reionization.

Another possibility is that we may be witnessing simultaneous births of young primeval galaxies occurring in high-density regions. Extremely young stellar ages may lead not only to a higher level of turbulence in the ISM creating low-density channels (Kakiichi & Gronke 2019) but also to a higher ionizing photon production efficiency (thus increasing Lyα production), higher specific SFRs (sSFRs; Clarke & Oey 2002; Endsley et al. 2021), and low dust content. Enhanced sSFRs in galaxies residing in high-density environment have been reported recently (Lemaux et al. 2020; Shi et al. 2020). While the positive SFR−density relation is relatively weak for the general population, there is also tentative evidence that LAEs—as young, low-mass galaxies—lie above the M_star–SFR main sequence (Guaita et al. 2011; Hagen et al. 2016; but see Kusakabe et al. 2018 for a contradictory result) in which such effects may manifest more clearly.

Finally, protocluster LAEs may be more metal-poor than the field LAEs, producing hotter stellar photospheres and thus higher ionizing radiation efficiency (Sternberg et al. 2003; Leitherer et al. 2010). While lower metallicity in high-density regions is counterintuitive, it is plausible if these galaxies are hosted by low-mass (≲10¹¹ M_⊙) halos undergoing the first major starburst fueled by the H i rich environment. At z ∼ 3, the volume of a single protocluster is immense (Chiang et al. 2013) and the chemical enrichment therein is expected to be heterogeneous.

Distinguishing these different but possibly related scenarios would require considerably larger samples of LAEs spanning a range of large-scale environments and deeper imaging data to enable photometric measurements with improved precision. Optical spectroscopy can enable comparative analyses of the Lyα properties (i.e., the shape and velocity offset) and dust content to test if the distribution of H i gas in the galaxy indeed changes with their large-scale structure. Robust measurements of their SFR, stellar masses, and the overall shape of the SED will also help discern a significant difference in metallicity and ages. Deep James Webb Space Telescope NIRSPEC observations can place direct constraints on the metallicity of galaxies in a diverse range of environments.

The One-hundred square-degree DECam Imaging in Narrowbands (ODIN) survey will provide the largest samples of LAEs in protocluster and field environments alike, which can serve as the basis of these investigations. As a new NSF NOIRLab program (Program ID: 2020B-0201), the survey began in early 2021 to image an area of 91 deg² with three narrowband filters sampling redshifted Lyα emission at z = 4.5, 3.1, and 2.4 (cosmic ages of 1.3, 2.0, and 2.7 Gyr, respectively), straddling the epoch in which the stellar mass assembly rate of both cluster and field galaxies reached its peak (Madau & Dickinson 2014). Within the survey volume of ≈0.24 Gpc³ (comoving), ≈100,000 LAEs, ≈45 Coma progenitors, and ≈600 Virgo progenitor systems are expected. Combined with the existing and upcoming facilities, ¹⁴ these new data sets will considerably enrich our understanding of the galaxy growth occurring in the largest cosmic structures in direct comparison with that of average galaxies at the same epoch.