A Comprehensive Study of Lyα Emission in the High-redshift Galaxy Population

Our sample was initially composed of 629 galaxies in the COSMOS, GOODS-S, and UDS fields. Every object was observed under the 3D-HST/CANDELS program, providing HST and Spitzer photometry from 3800 Å to 7.9 μm (9 bands for COSMOS, 14 for GOODS-S, and 9 for UDS). We construct our sample using 3D-HST outputs (Skelton et al. 2014). According to these, our 629 photometric redshifts satisfy $3.25\lt {z}_{3{\rm{D}} \mbox{-} {\text{}}{\rm{HST}}}\lt 4.25$ and have a 95% probability of 2.9 < z < 4.25. Every galaxy also complies with a photometric redshift reliability parameter ${Q}_{z}\leqslant 3$ selection to remove catastrophic outliers (Brammer et al. 2008). In terms of M_*, our galaxies are homogeneously distributed in the range $8\lt \mathrm{log}{({M}_{* }[{M}_{\odot }])}_{3{\rm{D}} \mbox{-} {\text{}}{\rm{HST}}}\lt 10.4$. These 3D-HST output values are based on the Bruzual & Charlot (2003) stellar population synthesis model library with a Chabrier (2003) IMF and solar metallicity. Exponentially declining star formation histories (SFHs) with a minimum e-folding time of ${\mathrm{log}}_{10}(\tau /\mathrm{yr})=7$ and a Calzetti et al. (2000) dust attenuation law were also assumed for the calculation (Skelton et al. 2014). We further emphasize that sample selection was performed based on the values of ${z}_{3{\rm{D}} \mbox{-} {\text{}}{\rm{HST}}}$ and M_* detailed in this section, whereas all further analysis in this paper is based on our own estimates (Section 2.3).

Spectroscopy of the full sample was conducted at the Magellan Clay 6.5 m telescope during 2014 December and 2015 February. To this end, we used the M2FS, a multi-object fiber-fed spectrograph. The 12 fibers of this instrument allow the observation of 256 targets within 30 arcmin in a single exposure. We are then able to observe 210 targets in each field simultaneously, while using 40 fibers for sky apertures and five for calibration stars. We used this spectrograph in LoRes mode, which features an expected resolution of R = 2000 and a continuum sensitivity of V = 24 with S/N = 5 in 2 hr (Mateo et al. 2012). The final data set consists of six exposure hours on each of the three fields with an average seeing of 06.

For data reduction, we developed a custom M2FS pipeline. This routine features standard bias subtraction and dark correction. For wavelength calibration, we use HgArNe lamps observed on each night. The wavelength solution is obtained separately for each fiber, with a typical rms uncertainty of 0.03 Å. We further correct the solution for each fiber using the sky lines in the science spectra. For flat-fielding, we use sky-flats obtained during either twilight or dawn. The correction is calculated separately for each fiber, and features illumination, fiber profile, and a correction to account for the shifting of fiber spectra on the detector due to thermal effects. Despite the fact that we are not using dome flats, we estimate the uncertainties in our flat-fielding to be <5%.

Sky subtraction is performed using the 40 sky fibers homogeneously distributed over the field of view and across the detector, where fibers are grouped in blocks. We find the sky solution to be more dependent on fiber location in the CCD than in the sky, leading us to perform sky subtraction separately for each frame fiber block. The solution is computed using a non-parametric spline fit, yielding satisfactory sky subtraction for most sky lines. For the bright sky lines that we are not able to properly subtract, we build emission line masks. For consistency, we use the same mask for all spectra, except for particularly noisy fibers. In those rare cases (11), we mask broad sections of the spectrum. Due to fiber malfunction, we could not obtain spectra for 4 of the 629 targets, leaving our final sample in 625 objects.

Flux calibration is performed using five ${M}_{V}=19\mbox{--}22$ calibration stars on each exposure. Due to atmospheric differential refraction (ADR) comparable to the 12 fiber size, we need to recover the intrinsic spectrum of each star. To do so, we fit stellar templates from the Pickles library (Pickles 1998) to the continuum-normalized instrumental spectra of the calibration stars. We then use the five stars on each exposure to obtain an average sensitivity curve. We estimate the rms uncertainty of our method to be about ∼15%. After this calibration, we correct the fluxes in our spectra for Galactic extinction. Samples of sky-subtracted spectra are shown in Figure 1.

Zoom In Zoom Out Reset image size

Considering the three fields, the science area we survey is ∼550 arcmin². The resulting spectral FWHM line resolution is ∼2 Å, and we reach a 1σ continuum flux density limit of $\sim 4\times {10}^{-19}$ erg s⁻¹ cm⁻² Å⁻¹ per pixel in our 6 hr of exposure. We estimate a 5σ emission line-flux sensitivity of $\sim 8\times {10}^{-18}$ erg s⁻¹ cm⁻² in our final spectra (see Figure 1). We identify 120 Lyα emission lines with S/N ≥ 5.5 in our data (details in Section 3.1). Given that we have 625 observed targets, we are then recovering a Lyα-emitting galaxy with S/N ≥ 5.5 every ∼5 objects.

As stated by Skelton et al. (2014), 3D-HST outputs for these galaxies were obtained using FAST (Kriek et al. 2009). These calculations assume exponentially declining SFHs with a minimum e-folding time of ${\mathrm{log}}_{10}(\tau /\mathrm{yr})=7$ (Skelton et al. 2014). However, since recent studies suggest that it is more adequate to reproduce high-redshift observations with constant SFHs (cSFHs; González et al. 2014), or even rising SFHs (Maraston et al. 2010), we perform our own executions of FAST assuming cSFHs. For a detailed discussion on this topic, we refer the reader to Conroy (2013). Our FAST executions also adopt the Bruzual & Charlot (2003) stellar population synthesis model library, a Chabrier (2003) IMF, and solar metallicity, similarly to Skelton et al. (2014). We do not account for nebular emission lines in our SED fitting, which, in principle, can overestimate our reported M_* by a factor of up to 4 (Atek et al. 2011; Conroy 2013; Stark et al. 2013) or even ∼7 for strong emitters (de Barros et al. 2014). However, our galaxies are at z ∼ 4. At this redshift, the overestimate is within a factor of ≲1.1 (Stark et al. 2013; Salmon et al. 2015), since the 4.5 μm Spitzer/IRAC band is unaffected.⁷ Our FAST outputs also include extinction-corrected SFRs. The SFRs are derived from the SED, and the extinction correction assumes a Calzetti et al. (2000) attenuation law with ${R}_{V}=4.05$.

We are also required to run EAZY (Brammer et al. 2008) on the photometry from CANDELS/IRAC, since FAST requires redshifts as input. The executions of EAZY yield a most probable redshift z_EAZY. Our 625 objects satisfy $3\lt {z}_{\mathrm{EAZY}}\lt 4.25$, with a median uncertainty ${\sigma }_{\mathrm{EAZY}}=0.1$ (Figure 2). EAZY outputs also include 2σ constraints, which for our sample are limited to 2.95 < z < 4.5. Our FAST executions yield a mass coverage of $7.6\lt \mathrm{log}({M}_{* }[{M}_{\odot }])\lt 10.6$ (Figure 3), with a characteristic uncertainty of $\mathrm{log}({M}_{* }[{M}_{\odot }])\sim 0.3$. Most of our SFR values are in the range $1\mbox{--}100\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$ (further analysis and plots below). From now on, we use the FAST outputs based on the spectroscopic redshifts (${z}_{\mathrm{Ly}\alpha }$) for the 120 detections (S/N ≥ 5.5; Section 3.1) and the outputs based on z_EAZY for the 505 non-detections.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The flux calibration procedure performed in our spectra is based on using stars; therefore, it corrects for a ∼32% fiber flux loss, which corresponds to a point-source Lyα surface brightness distribution. In cases of extended Lyα emission halos (Steidel et al. 2011; Matsuda et al. 2012; Feldmeier et al. 2013; Hayes et al. 2013; Momose et al. 2014; Caminha et al. 2016; Patrício et al. 2016; Wisotzki et al. 2016), the fluxes we derive are mostly associated with the galaxies themselves and the inner parts of their Lyα halos. Furthermore, Lyα emission from a galaxy can show significant misalignment with the UV continuum (Rauch et al. 2011), although such cases are not the norm (Finkelstein et al. 2011; Jiang et al. 2013a). At the cosmic time of our sample, the fiber diameter corresponds to a scale of ∼8 kpc, roughly a factor of four larger than the typical effective diameter of galaxies (Bond et al. 2012; Law et al. 2012).

We find a median redshift offset of ${\rm{\Delta }}z={z}_{\mathrm{Ly}\alpha }\mbox{--}{z}_{\mathrm{EAZY}}\,=0.24$ in our detections (see Figure 2). We show in Oyarzún et al. (2016) that this offset does not have any noticeable effects on our sample dependence on M_*, which is the primary selection criterion for our galaxies. A very similar offset has also been found in the MUSE-Wide Survey (Herenz et al. 2017). We find the offset to correlate with ${W}_{\mathrm{Ly}\alpha }$ (see Figure 2), hinting at biases in photometric redshift fitting when a strong Lyα emission line is present. A plausible scenario is one in which the Lyman break from the template is fitted slightly blueshifted to account for the flux excess in the redder band. Thorough simulations beyond the scope of this work are required to explore the causes behind this bias.

We stress that the 3D-HST/CANDELS mass incompleteness is restricted to $\mathrm{log}({M}_{* }[{M}_{\odot }])\lt 8.5$ at z ∼ 4 (at least for GOODS-S; Duncan et al. 2014), which corresponds to about one-quarter of the sample. We want to stress that this is a homogeneously M_*-selected sample, designed to study Lyα emission statistics dependence on galaxy properties. As a consequence, it is by no means representative of the M_* or Lyα LFs at 3 < z < 4.6. This must be taken into account when comparing this sample to analogues directly drawn from the galaxy population (i.e., LBGs or narrowband samples).

For line detection, we use an automated maximum likelihood fitting routine after continuum subtraction. We assume intrinsic Gaussian profiles of the form

$\begin{eqnarray}&&{f}_{\mathrm{rest}}(\lambda )=\displaystyle \frac{F}{\sqrt{2\pi {\sigma }_{\lambda }^{2}}}{\rm{\exp }}\left[-\displaystyle \frac{{(\lambda -{\lambda }_{0})}^{2}}{2{\sigma }_{\lambda }^{2}}\right],\end{eqnarray} \tag{ 1 }$

where F, ${\lambda }_{0}$, and ${\sigma }_{\lambda }$ compose the parameter space explored by the maximum likelihood.

In order to account for false positives, we run the line detection routine on the 115 sky fibers. The results are shown in Figure 4. We detect four lines above 4σ and none above 5σ. Therefore, down to 5σ, we expect at most two false detections in our 120 lines, translating into ≲2% contamination using signal-to-noise S/N* = 5.5 as our threshold. We also characterize our detection completeness (see Figure 4). To obtain it, we use p(S/N${}_{i}\,\gt$ S/N ${}^{* }\,|$S/N), with S/N_i the measured signal-to-noise ratio and S/N* the imposed detection threshold. We define the simulated signal to noise as ${\rm{S}}/{{\rm{N}}}_{\mathrm{sim}}=F/\sqrt{\sum d{\lambda }^{2}{\sigma }_{k}^{2}}$, with $d\lambda$ the wavelength dispersion in the spectrum and ${\sigma }_{k}$ the flux uncertainty for pixel k. We find the most accurate representation by summing over an interval of 20 Å centered at 5500 Å. To recover the actual completeness shown in Figure 4, we simulate $\sim {10}^{3}$ lines on the 115 sky-spectra sampling fluxes of ${10}^{-19}\mbox{--}{10}^{-17}$ ergs s⁻¹ cm⁻², FWHMs between 5–13 Å, and wavelengths of 4800–6700 Å. The fraction for which we measure S/N_i > S/N* is our completeness. We use S/N* = 5.5 as our detection threshold, which corresponds to a line-flux sensitivity in the final spectra of $\sim 1\times {10}^{-17}$ erg s⁻¹ cm⁻² (see Figure 1).

Zoom In Zoom Out Reset image size

The radiative transfer and escape of Lyα radiation from galaxies can be highly complicated. As a matter of fact, the resonant nature of this line has led to thorough modeling of its radiative escape (e.g., Verhamme et al. 2006; Dijkstra & Kramer 2012). Such complications imply that the flux profile of a Lyα line is not always well reproduced by the usual Gaussian profile (Chonis et al. 2013; Trainor et al. 2015). Hence, for flux measurements, we adopt a more sophisticated model. Similarly to McLinden et al. (2011) and Chonis et al. (2013), we fit double-peaked Gaussian profiles of the form

$\begin{eqnarray}&&{f}_{\mathrm{rest}}(\lambda )={f}_{\mathrm{blue}}(\lambda )+{f}_{\mathrm{red}}(\lambda ),\end{eqnarray} \tag{ 2 }$

where ${f}_{\mathrm{blue}}(\lambda )$ represents the blue emission component and ${f}_{\mathrm{red}}(\lambda )$ represents the red component. We assume each of these components to be asymmetric, i.e., they follow Equation (1), with σ defined as

$\begin{eqnarray}{\sigma }_{\lambda } & \equiv & {\sigma }_{b}\ {\rm{if}}\ \lambda \lt {\lambda }_{0},\\ {\sigma }_{\lambda } & \equiv & {\sigma }_{r}\ {\rm{if}}\ \lambda \gt {\lambda }_{0}.\end{eqnarray} \tag{ 3 }$

Before fitting, we convolve the profiles given by Equations (2) and (3) with the spectral resolution. This allows us to properly characterize the errors in our measurements and methodology. In case there is some sky contamination (which happens for 28 of the 120 detections), we only fit single-peaked profiles. The resulting fits to six of our emission lines are shown in Figure 1.

There are typically two diagnostics used to characterize the prominence of Lyα emission in galaxies: the equivalent width (${W}_{\mathrm{Ly}\alpha }$) and the Lyα escape fraction (f_esc). The rest-frame ${W}_{\mathrm{Ly}\alpha }$ is defined as the fraction between line flux and UV continuum flux in the rest frame of the galaxy. Explicitly,

$\begin{eqnarray}&&{W}_{\mathrm{Ly}\alpha }=\displaystyle \frac{F}{{f}_{\lambda }}\displaystyle \frac{1}{(1+{z}_{\mathrm{Ly}\alpha })},\end{eqnarray} \tag{ 4 }$

with F the Lyα flux we measure in the spectra and ${f}_{\lambda }$ the observed flux at rest frame 1700 Å from 3D-HST rest-frame colors (Skelton et al. 2014). It must be noted that the 3D-HST rest-frame colors come from the best-fit template to the photometry, i.e., they are not direct measurements. Still, this allows us to use, in principle, the same rest-frame wavelength for every object, regardless of its actual redshift. In addition, we can also derive a value of ${f}_{\lambda }$ even if the rest-frame 1700 Å photometry is missing. For reference, we find that the templates match the photometry at ∼1700 Å of every object to a typical agreement of ∼90%. For the uncertainties on ${f}_{\lambda }$, however, we use the errors on the photometry. We do not use any ${f}_{\lambda }$ directly measured from our data, since every galaxy has a continuum fainter or comparable to the 1σ errors in the spectra.

On the other hand, f_esc is the fraction of the number of Lyα photons that escape the galaxy from the number produced. This diagnostic is typically indirectly recovered using SFRs derived from the Lyα line and intrinsic SFRs. The latter are typically calculated from UV continuum measurements or Hα fluxes, subject to extinction correction. In our case, we use the intrinsic (i.e., extinction-corrected) SFRs from FAST. The observed SFRs are derived from the SED, whereas the extinction corrections assume a Calzetti et al. (2000) attenuation law with ${R}_{V}=4.05$. The explicit definition is then (e.g., Blanc et al. 2011)

$\begin{eqnarray}{f}_{\mathrm{esc}}(\mathrm{Ly}\alpha ) & = & \displaystyle \frac{L{(\mathrm{Ly}\alpha )}_{\mathrm{obs}}}{L{(\mathrm{Ly}\alpha )}_{\mathrm{intrinsic}}}\\ & = & \displaystyle \frac{{\rm{SFR}}(\mathrm{Ly}\alpha )}{{\rm{SFR}}{({\rm{UV}})}_{\mathrm{corr}}}=\displaystyle \frac{4.4\times {10}^{-42}{L}_{\mathrm{Ly}\alpha }/8.7}{{\rm{SFR}}{({\rm{UV}})}_{\mathrm{corr}}},\end{eqnarray} \tag{ 5 }$

with ${L}_{\mathrm{Ly}\alpha }$ the luminosity associated with the Lyα line and SFR(UV)_corr the extinction-corrected SFRs. This assumes a Lyα to Hα ratio of 8.7 and the Hα SFR calibration for a Chabrier (2003) IMF of $4.4\times {10}^{-42}$ (Matthee et al. 2016). In case of non-detections, we use the corresponding S/N* = 5.5 line fluxes to derive upper limits on ${W}_{\mathrm{Ly}\alpha }$ and f_esc. Our calculations do not account for the UV excess associated with binaries, which induce uncertainties on f_esc that can be very difficult to account for (Stanway et al. 2016).

We do not measure any objects to have more than one emission line with ${\rm{S}}/{{\rm{N}}}_{i}\gt {\rm{S}}/{{\rm{N}}}^{* }$, dismissing the existence of evident active galactic nuclei (AGNs) in our sample. We do not find any of the 120 Lyα-emitting galaxies to have S/N > 1 emission potentially associated with the C iv λ1549 line. Moreover, we find no broad-line Lyα emission (${\rm{\Delta }}v\,\gt 1000$ km s⁻¹), consistent with the rarity of Type I AGNs at high redshift (Dawson et al. 2004), and in strong contrast with the AGN fractions at low redshift (∼0.4; Finkelstein et al. 2009). This is consistent with the fact that bright AGNs are found in ≲1% of z > 3 LAEs (Malhotra et al. 2003; Wang et al. 2004; Gawiser et al. 2007; Ouchi et al. 2008; Zheng et al. 2010). Low-luminosity AGN contamination in LAE samples is more difficult to constrain, but Wang et al. (2004) estimate the total AGN fraction to be <17%. We use our f_esc measurements to further dismiss the existence of evident AGNs. We perform cross-matching with the NASA/IPAC Extragalactic Database (NED⁸ ) for all objects with ${f}_{\mathrm{esc}}\gt 1$, and we find none of them to be a reported X-ray source.

Measurement of ${W}_{\mathrm{Ly}\alpha }$ for a galaxy sample yields the ${W}_{\mathrm{Ly}\alpha }$ distribution. Since ${W}_{\mathrm{Ly}\alpha }$ is directly measured from the data, characterization of the ${W}_{\mathrm{Ly}\alpha }$ distribution can be naively considered straightforward. However, careful consideration of uncertainties and completeness can yield important insights into the underlying information. Therefore, for proper characterization of uncertainties, significance, and trends in our results, we use Bayesian statistics. In this section, we explain how to recover the ${W}_{\mathrm{Ly}\alpha }$ distribution within this framework, complementary to the one introduced in Treu et al. (2012).

Different probability distribution models can be adopted to reproduce ${W}_{\mathrm{Ly}\alpha }$ distribution measurements. For instance, studies use Gaussian (e.g., Guaita et al. 2010), exponential (e.g., Jiang et al. 2013b; Zheng et al. 2014), and log-normal distributions. For our analysis, we define the probability distribution as $p({W}_{\mathrm{Ly}\alpha }| {W}_{n})$. Let W_n be the parameter space associated with the model. From now on, we describe the Bayesian approach to recover the posterior distribution of W_n. By means of this approach, we include the uncertainties in sample size, flux measurements, and photometry in the estimation of the posterior. The description we provide is not limited to this particular work, allowing for further application in similar data sets. We only present here the fundamental equations, as this procedure is already described in detail in Oyarzún et al. (2016).

Our Bayesian analysis is based on the Lyα line flux F instead of ${W}_{\mathrm{Ly}\alpha }$, as introduced in Equation (4). This approach simplifies the equations, since we can assume F and ${f}_{\lambda }$ to be normally distributed, which cannot be done for ${W}_{\mathrm{Ly}\alpha }$. According to Bayes' theorem, the posterior distribution $p({W}_{n}| \{F\})$, i.e., the parameter space probability distribution given our data set $\{F\}$, is

$\begin{eqnarray}&&p({W}_{n}| \{F\})=\displaystyle \frac{p(\{F\}| {W}_{n})p({W}_{n})}{p(\{F\})}.\end{eqnarray} \tag{ 6 }$

The likelihood is just the product of the individual likelihood for every galaxy, i.e., $p(\{F\}| {W}_{n})=\prod p({F}_{i}| {W}_{n})$. For a detection, it is given by

$\begin{eqnarray}&&p({F}_{i}| {W}_{n})={\int }_{0}^{\infty }p({F}_{i}| F)p(F| {W}_{n}){dF},\end{eqnarray} \tag{ 7 }$

where $p({F}_{i}| F)$ is the line-flux probability distribution for the corresponding galaxy valued at a flux F, which we consider to distribute normally. On the other hand, the term $p(F| {W}_{n})$ is just the probability distribution of F given W_n. Using the definition of ${W}_{\mathrm{Ly}\alpha }$ from Equation (4), the term translates to the product distribution $p(F| {W}_{n})=p({W}_{\mathrm{Ly}\alpha }| {W}_{n})p({f}_{\lambda })$. At this point, we include the probability distribution for the continuum, which we also assume to be Gaussian.

The limiting line flux ${F}_{i}^{* }$ for discerning detections from noise is given by our S/N threshold, i.e., ${F}_{i}^{* }={\rm{S}}/{{\rm{N}}}^{* }{\sigma }_{i}$. For galaxies with no detections above ${F}_{i}^{* }$, we adopt the following value for the likelihood:

$\begin{eqnarray}&&p({F}_{i}\lt {F}_{i}^{* }| {W}_{n})={\int }_{0}^{\infty }(1-p({F}_{i}\gt {F}_{i}^{* }| F))p(F| {W}_{n}){dF},\end{eqnarray} \tag{ 8 }$

with $p({F}_{i}\gt {F}_{i}^{* }| F)$ the detection completeness at a line flux F (see Section 3.1).

Using the expressions for detections and non-detections, the posterior distribution takes the final form:

$\begin{eqnarray}&&p({W}_{n}| \{F\})=\displaystyle \frac{p({W}_{n})}{p(\{F\})}\displaystyle \prod _{D}p({F}_{i}| {W}_{n})\displaystyle \prod _{{ND}}p({F}_{i}\lt {F}_{i}^{* }| {W}_{n}),\end{eqnarray} \tag{ 9 }$

with $p({W}_{n})$ the priors of the model parameters and $p(\{F\})$ a normalization constant reflecting the likelihood of the model. The form of Equation (9) is general and will be used as a starting point for multiple analysis throughout.

As we introduced in the previous section, multiple probability distributions can be adopted for the representation of the ${W}_{\mathrm{Ly}\alpha }$ distribution. In order to perform model selection, several elements are taken into account, such as model complexity and number of parameters. In this section, we describe our methodology to perform such a selection from a quantitative standpoint. By means of a Bayesian approach, we recover probability ratios for the different models, providing insight into how we perform the selection given our measurements. The analysis presented here is a quantitative implementation of Occam's Razor and is not unique to our data set, i.e., it can be applied to any data set modeling.

Every model we discuss here is composed of a scaled probability distribution and a Dirac delta, as proposed in Treu et al. (2012). If we define the standard probability distributions as ${p}_{0}({W}_{\mathrm{Ly}\alpha }| {W}_{n})$, the modified counterparts we consider are given by

$\begin{eqnarray}p({W}_{\mathrm{Ly}\alpha }| {W}_{n}) & = & A\times H({W}_{\mathrm{Ly}\alpha })\times {p}_{0}({W}_{\mathrm{Ly}\alpha }| {W}_{n})\\ & & +(1-A)\times \delta ({W}_{\mathrm{Ly}\alpha }),\end{eqnarray} \tag{ 10 }$

where the first term is the scaled probability distribution. It is multiplied by the Heaviside $H({W}_{\mathrm{Ly}\alpha })$ to ensure it only represents positive ${W}_{\mathrm{Ly}\alpha }$ values. Hence, this scaled term integrates A, i.e., A can only adopt values between 0 and 1. Note that this term is different from the fraction of detections in our sample, as we consider upper limits for our non-detections. The second term groups the fraction of galaxies that do not emit in Lyα (i.e., no line and/or absorption). As our data is restricted to emission lines, we represent this term using the Dirac delta.

We explore here exponential-, Gaussian-, and log-normal-based distributions. The first two are two-parameter models, while the log-normal is three-parameter dependent. Hence, we generalize our parameter space as ${W}_{n}=(A,{W}_{0},\mu )$. Then, the expressions for our exponential, Gaussian, and log-normal models are, respectively,

$\begin{eqnarray}p({W}_{\mathrm{Ly}\alpha }| {\rm{\exp }}) & = & \displaystyle \frac{A}{{W}_{0}}{e}^{-{W}_{\mathrm{Ly}\alpha }/{W}_{0}}H({W}_{\mathrm{Ly}\alpha })\\ & & +(1-A)\delta ({W}_{\mathrm{Ly}\alpha })\end{eqnarray} \tag{ 11 }$

$\begin{eqnarray}p({W}_{\mathrm{Ly}\alpha }| {\rm{Gauss}}) & = & \displaystyle \frac{2A}{\sqrt{2\pi {W}_{0}^{2}}}{e}^{-{W}_{\mathrm{Ly}\alpha }^{2}/2{W}_{0}^{2}}H({W}_{\mathrm{Ly}\alpha })\\ & & +(1-A)\delta ({W}_{\mathrm{Ly}\alpha })\end{eqnarray} \tag{ 12 }$

$\begin{eqnarray}p({W}_{\mathrm{Ly}\alpha }| \mathrm{log}\,n) & = & \displaystyle \frac{A}{\sqrt{2\pi {W}_{0}^{2}{W}_{\mathrm{Ly}\alpha }^{2}}}{e}^{-{(\mathrm{ln}({W}_{\mathrm{Ly}\alpha })-\mu )}^{2}/2{W}_{0}^{2}}\\ & & +(1-A)\delta ({W}_{\mathrm{Ly}\alpha }).\end{eqnarray} \tag{ 13 }$

We now describe our approach for comparing the three models. For a set of measurements, Bayes' theorem gives the probability of model M_i in the model space $\{M\}$:

$\begin{eqnarray}&&p({M}_{i}| \{F\})=\displaystyle \frac{p(\{F\}| {M}_{i})p({M}_{i})}{p(\{F\})},\end{eqnarray} \tag{ 14 }$

with $p({M}_{i})$ the prior for model M_i in the set $\{M\}$, which we assume to be equal for the three models. Once again, $p(\{F\})$ is a normalization constant. Therefore, the probability for model M_i is proportional to the likelihood of the model, i.e.,

$\begin{eqnarray}&&p({M}_{i}| \{F\})\propto p(\{F\}| {M}_{i})={\int }_{n}p({W}_{n})p(\{F\}| {W}_{n}){{dW}}_{n}.\end{eqnarray} \tag{ 15 }$

The absolute probability of each model M_i within the set $\{M\}$ is obtained by imposing that the models explored cover all possible choices, i.e., $\sum p({M}_{i}| \{F\})=1$. Hence, the probability of M_i given our data set is

$\begin{eqnarray}&&p({M}_{i}| \{F\})=\displaystyle \frac{{\int }_{n}p({W}_{n})p(\{F\}| {W}_{n}){{dW}}_{n}}{\sum p({M}_{i}| \{F\})}.\end{eqnarray} \tag{ 16 }$

For an analytically correct model comparison, analysis of Equation (16) is required. Still, the term in Equation (15) is strongly dependent on the priors assumed for the parameter space W_n of every model. Therefore, different prior selections can have significant effects on the odds for each model. As a workaround, we rewrite the model probabilities as

$\begin{eqnarray}&&p({M}_{i}| \{F\})=\displaystyle \frac{{\int }_{n}p(\{F\}| {W}_{n}){{dW}}_{n}}{\sum p({M}_{i}| \{F\})}.\end{eqnarray} \tag{ 17 }$

Effectively, this simplification conveniently limits our analysis to a pure likelihood comparison, i.e., we adopt constant, uninformative priors. This is equivalent to assuming ignorance in linear scales for every parameter. Since our distributions are smooth and single peaked, we are confident in this assumption. The dimensions of our parameter spaces do not go beyond three, and the uncertainties in our parameters are of the order of the most probable values, providing further assurance to our assumptions.

In Oyarzún et al. (2016), we show that a galaxy sample with a broad M_* range yields a composite ${W}_{\mathrm{Ly}\alpha }$ distribution. For the rest of this section, we divide our sample into three M_* bins, as shown in Figure 3. We use the complete sample and these three subsamples to contrast the models. The outcomes are presented in Table 1 and Figure 5. Table 1 gives evidence that the best likelihood is obtained with the log-normal model for three of the four distributions. This can be verified in our distribution simulations for the complete sample in Figure 5, especially toward the high ${W}_{\mathrm{Ly}\alpha }$ tail. Still, when integrating the likelihoods, the log-normal distribution is the least probable. This is a consequence of the extra parameter needed by the model, which penalizes the likelihood when integrating over the parameter space. Then, according to our analysis, the preferred model is the exponential. While it models the distribution better than the Gaussian, it also reproduces our ${W}_{\mathrm{Ly}\alpha }$ measurements fairly well, despite depending on only two parameters. In addition, the uncertainties in Figure 5, especially for low ${W}_{\mathrm{Ly}\alpha }$, confirm this model is the most adequate to reproduce our measurements. We remark that the procedure for model selection described here considers the lower ${W}_{\mathrm{Ly}\alpha }$ end of the distribution, which includes our completeness and non-detections. Nonetheless, further ${W}_{\mathrm{Ly}\alpha }$ distribution analysis in this paper is mostly focused on the higher ${W}_{\mathrm{Ly}\alpha }$ end and is not strongly dependent on model preference.

Zoom In Zoom Out Reset image size

From now on, we perform our ${W}_{\mathrm{Ly}\alpha }$ distribution analysis using the exponential model of Equation (11). We stress that this expression is dependent on the parameters A and W₀, with the first being the fraction of galaxies showing emission and the second the e-folding scale of the distribution. We advise caution when using A as a proxy for the fraction of line emitters in the parent population, since it is tied to the adequacy of an exponential profile.

Evidence suggests Lyα emission is strongly dependent on the M_* of galaxies. Galaxies with higher M_* have been forming stars for longer, leading to greater ISM dust that presumably forms in supernovae and AGB stars (Silva et al. 1998). A greater dust content leads to more Lyα photon absorption, decreasing ${W}_{\mathrm{Ly}\alpha }$. This effect has already been observed, at least for high ${W}_{\mathrm{Ly}\alpha }$, in Blanc et al. (2011) and Hagen et al. (2014). Similarly, the bulk of M_* is dominated by older stars, which do not contribute significantly to the Lyα photon budget of galaxies. As a matter of fact, Lyα emission decreases steadily with the age of stellar populations, as seen in Charlot & Fall (1993) and Schaerer (2003). Lyα radiative transfer is also severely affected by the neutral gas structure and kinematics of the ISM and circumgalactic medium (Verhamme et al. 2006). Since more massive star-forming galaxies are bound to have higher gas mass (e.g., Kereš et al. 2005; Finlator et al. 2007), Lyα photons should be subject to more resonant scattering, therefore decreasing their ${W}_{\mathrm{Ly}\alpha }$. The trends we find in Oyarzún et al. (2016) confirm this qualitative scheme at $z\sim 4$. In this section, we perform a more detailed and robust characterization of Lyα emission dependence on M_*.

Our sample is especially designed to study the dependence of Lyα emission in M_*. As shown in Figure 3, our objects are selected in redshift and M_*, homogeneously covering the range $7.6\lt \mathrm{log}({M}_{* }[{M}_{\odot }])\lt 10.6$. For further clarity, we plot in Figure 6 the M_* and SFRs of our complete sample and detections. Comparison of the M_* histograms slightly hints at an anti-correlation between the LAE fraction and M_*, at least down to our detection limit. However, since our completeness depends on M_*, a more thorough analysis is required (see below). The existence of any Lyα emission dependence on M_* becomes clearer in Figure 7, where we plot ${W}_{\mathrm{Ly}\alpha }$ and f_esc as a function of M_*. We also plot in Figure 7 our upper limits for non-detections. The regions sampled by our non-detections reveal how our completeness is not independent of M_*. Our detections are flux limited, so we achieve lower ${W}_{\mathrm{Ly}\alpha }$ for galaxies with brighter UV continuum. Therefore, even though we observe lower M_* galaxies to have higher ${W}_{\mathrm{Ly}\alpha }$ and f_esc, any qualitative conclusions we can draw involving the fraction of LAEs as a function of M_* are affected by our completeness. Still, there is a clear upper envelope to the distribution of galaxies in this plot, where we are not affected by incompleteness. For ${M}_{* }\lt {10}^{8.5}$, there is a clear anti-correlation between ${W}_{\mathrm{Ly}\alpha }$ and M_*. For more massive galaxies, however, the trend is mostly flat, except for the presence of a few interlopers. Still, our qualitative result is that both Lyα emission diagnostics show an anti-correlation with M_*. For the rest of this section, we focus on how to thoroughly quantify the effect of M_* on ${W}_{\mathrm{Ly}\alpha }$.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The overall dependence of Lyα emission on M_* implies that ${W}_{\mathrm{Ly}\alpha }$ distributions in the literature (e.g., Gronwall et al. 2007; Zheng et al. 2014) are influenced by the M_* distribution of the sample. In comparison to deeper M_UV surveys, shallower samples are bound to observe lower ${W}_{\mathrm{Ly}\alpha }$. This can lead to incorrect contrast of surveys and misinterpretation of trends. To verify these claims, we divide our sample into three M_* bins (see Figure 3) and plot the resulting ${W}_{\mathrm{Ly}\alpha }$ distributions in Figure 8. As expected, there is an apparent anti-correlation between M_* and both the tail of the ${W}_{\mathrm{Ly}\alpha }$ distribution and the normalization. We perform our first quantitative characterization of the ${W}_{\mathrm{Ly}\alpha }$ distribution dependence on M_* in Oyarzún et al. (2016). In that study, we divide our sample into three M_* bins and obtain the posterior distribution for the exponential parameters separately. This procedure allowed us to fit a linear relation to the final parameters, recovering A(M${}_{* }$) and W₀(M${}_{* }$) using expressions of the form

$\begin{eqnarray}&&A({M}_{* })={A}_{{M}_{* }}\mathrm{log}({M}_{* }[{M}_{\odot }])+{A}_{C},\end{eqnarray} \tag{ 18 }$

$\begin{eqnarray}&&{W}_{0}({M}_{* })={W}_{{M}_{* }}\mathrm{log}({M}_{* }[{M}_{\odot }])+{W}_{C}.\end{eqnarray} \tag{ 19 }$

In this section, we recover these linear relations directly from the complete sample, i.e., we recover the posterior distribution for the four-parameter space composed of the linear coefficients in Equations (18) and (19). This more robust methodology does not rely on binning, while also allowing us to constrain the errors on the coefficients directly from the model and measurements. We once again start from Equation (9). As mentioned, our parameter space is now ${W}_{n}=({A}_{{M}_{* }},{A}_{C},{W}_{{M}_{* }},{W}_{C})$. As these linear coefficients represent the exponential parameters of Equation (11), deriving non-informative priors is highly complicated. Therefore, in order to determine our priors, we only consider linear scales ignorance. Therefore, the posterior translates to

$\begin{eqnarray}&&p({A}_{{M}_{* }},{A}_{C},{W}_{{M}_{* }},{W}_{C}| \{F\})=C\times \displaystyle \prod _{D}p({F}_{i}| A({M}_{* }),{W}_{0}({M}_{* }))\\ &&\quad \times \displaystyle \prod _{{ND}}p({F}_{i}\lt {F}_{i}^{* }| A({M}_{* }),{W}_{0}({M}_{* })),\end{eqnarray} \tag{ 20 }$

with C a normalization constant.

Zoom In Zoom Out Reset image size

We use MCMC simulations to characterize this four-parameter posterior. Its maximum gives the best solution for $A({M}_{* })$ and ${W}_{0}({M}_{* })$, while the collapsed posteriors yield the uncertainties on the parameters. We can then write Equations (18) and (19) as

$\begin{eqnarray}&&A({M}_{* })=-{0.28}_{-.02}^{+.12}\mathrm{log}({M}_{* }[{M}_{\odot }])+{3.1}_{-1.1}^{+0.15},\end{eqnarray} \tag{ 21 }$

$\begin{eqnarray}&&{W}_{0}({M}_{* })=-{17.7}_{-2.0}^{+2.6}\mathrm{log}({M}_{* }[{M}_{\odot }])+{190}_{-25}^{+20}.\end{eqnarray} \tag{ 22 }$

In the framework of an exponential profile, these relations we recover yield the ${W}_{\mathrm{Ly}\alpha }$ probability distribution for an object with known M_*. Hence, we can simulate the expected ${W}_{\mathrm{Ly}\alpha }$ distribution for each of the three M_* subsamples and compare with our direct measurements. The results are presented in Figure 8. Our constraints are consistent with the observed ${W}_{\mathrm{Ly}\alpha }$ distributions.

Our results regarding the dependence of Lyα on M_* are conclusive. In the range ${10}^{8}\mbox{--}{10}^{10.5}\,{M}_{\odot }$, the ${W}_{\mathrm{Ly}\alpha }$ probability distribution extends to higher ${W}_{\mathrm{Ly}\alpha }$ for lower M_* galaxies. In other words, more massive galaxies tend to have lower ${W}_{\mathrm{Ly}\alpha }$. A similar trend is observed for f_esc, further highlighting the role of dust and gas mass in the escape of Lyα photons. At z ∼ 2.2, Matthee et al. (2016) also observe the f_esc anti-correlation from Figure 7, although only when stacking galaxies. Their massive objects showing high f_esc, which they associate with dusty gas outflows, seem to lie below the M_*–SFR sequence at z ∼ 2. These inferences, combined with the significantly higher f_esc they measure for larger apertures, make sense in a Lyα diffuse halo scheme, which we do not observe due to our aperture size. Studies on the dependence of f_esc on M_* also explore 1.9 < z < 3.6 (Hagen et al. 2014). Their ${W}_{\mathrm{Ly}\alpha }$-selected LAEs follow a trend similar to the one we find at z ∼ 4. In summary, the evidence for an anti-correlation between ${W}_{\mathrm{Ly}\alpha }$ (or f_esc) and M_* is significant, but the scatter seems to depend on measurement methodology and sample selection.

Inferences on the M_* distribution of LAEs are not as evident. Hagen et al. (2014) do not find their 1.9 < z < 3.6 Lyα luminosity-selected LAE number distribution to depend on M_*. Their results agree with the z ∼ 3.1 narrowband-selected survey of McLinden et al. (2014). However, we show in Figure 7 that spectroscopic completeness is not independent of M_*. Therefore, most LAEs survey follow-ups could have higher incompleteness toward lower M_*. Since our Bayesian analysis takes into account our completeness for every object, we can test the significance of this claim. The coefficient ${A}_{{M}_{* }}$ in Equation (18) represents the exponential fraction of LAE dependence on M_*. As evidenced by Equation (21), our measurements are more than $2\sigma$ consistent with a decrease in the fraction going to higher M_*. In the scheme of an exponential model, this translates to an LAE distribution dominated by lower M_* galaxies. This result complements the much more significant anti-correlation between W₀ and M_* we find (see Equation (22)).

High-redshift galaxies have been observed to follow a correlation between SFR and M_*, known as the star-forming main sequence (e.g., Kereš et al. 2005; Finlator et al. 2007; Stark et al. 2009; González et al. 2011; Whitaker et al. 2012; see our Figure 6). In terms of the underlying physics, more massive objects dominate gas accretion in their neighborhood, feeding and triggering star formation. Such gas infall seems to dominate over galaxy growth at high redshift (Kereš et al. 2005; Finlator et al. 2007). This scheme implies that more massive objects form stars at higher rates, at least down to our observational limitations and modeling of high-redshift ISM. Given our results on M_* from the previous section, we expect similar trends between ${W}_{\mathrm{Ly}\alpha }$ (f_esc) and SFR (Figure 6). Even more, star-forming galaxies have a higher neutral gas mass, which can hamper the escape of Lyα photons from galaxies (Verhamme et al. 2006). In fact, it has also been suggested that photoelectric absorption rules Lyα depletion, even over dust attenuation (Reddy et al. 2016). In this section, we explore any Lyα dependence on SFR within our data set. We remark that our SFRs come from SED fitting of 3D-HST photometry using FAST (see Section 2.3), i.e., they have typical associated timescales of 100 Myr (Kennicutt 1998). We stress that our derived SFRs differ from 3D-HST SFRs, since our calculation assumes cSFHs instead of exponentially declining SFHs (Skelton et al. 2014).

We show the ${W}_{\mathrm{Ly}\alpha }$ and f_esc dependence on SFR in Figures 6 and 9. In the latter, we include upper limits for our non-detections to give an insight into how our incompleteness depends on SFR. A clear anti-correlation between ${W}_{\mathrm{Ly}\alpha }$ (f_esc) and SFR is observed. These results come as no surprise, as they have been previously reported. Most studies of ${W}_{\mathrm{Ly}\alpha }$ dependence on SFR involve uncorrected SFRs (Pettini et al. 2002; Shapley et al. 2003; Yamada et al. 2005; Gronwall et al. 2007; Tapken et al. 2007; Ouchi et al. 2008; all compiled in Verhamme et al. 2008). Even without dust correction, the anti-correlation is still present in these studies (refer to Figure 19 in Verhamme et al. 2008 and Section 5.4 of this work). Based on a z ∼ 2 Hα emitters sample, Matthee et al. (2016) also observe a clear anti-correlation between Lyα f_esc and SFR. Interestingly, they do not only observe such a trend in their individual objects, but likewise on their stacks when using different apertures (galaxy diameters of 12 and 24 kpc). As their data set includes Hα fluxes, they can recover SFRs and f_esc using Hα luminosities. The fact that they observe similar trends with such a different sample suggests that the anti-correlation between ${W}_{\mathrm{Ly}\alpha }$ (f_esc) and SFR is not only independent of redshift, but also observational constraints like aperture and methodology for recovering SFRs. Their comparison, however, is restricted to SFRs higher than $\sim 5\,{M}_{\odot }$ yr⁻¹. Most of our low-mass objects have SFRs lower than $\sim 5\,{M}_{\odot }$ yr⁻¹, but they seem to follow the same regime as the rest of our sample. Even though f_esc uncertainties and incompleteness increase toward lower-SFR, UV-fainter galaxies (Figure 9), our results suggest that f_esc reaches values of 100% toward $\mathrm{SFR}\sim 1\mbox{--}3\,{M}_{\odot }$ yr⁻¹. These numbers are consistent with the analysis by Atek et al. (2014). They compare their z < 0.5 SFR(Lyα)/SFR(UV) measurements with the literature at z > 2 (Taniguchi et al. 2005; Gronwall et al. 2007; Guaita et al. 2010; Curtis-Lake et al. 2012; Jiang et al. 2013b). Our finding that f_esc reaches values of 100% toward $\mathrm{SFR}\sim 1\mbox{--}3\,{M}_{\odot }$ yr⁻¹ overestimates SFR(Lyα)/SFR(UV) for z < 0.5, but is consistent with higher redshift. As pointed by Atek et al. (2014), it seems that f_esc at fixed SFR(UV) increases with redshift. Then again, it must be kept in mind that these literature results consider uncorrected SFRs.

Zoom In Zoom Out Reset image size

In this section, we analyze the M_UV distribution of our sample, while also exploring any correlations between ${W}_{\mathrm{Ly}\alpha }$ and UV luminosities. It must be noted, though, that our sample is not representative of the galaxy population at 3 < z < 4. First, our galaxies are homogeneously distributed in M_*, i.e., they are not a random sample from 3 < z < 4 CANDELS objects. Second, we are affected by CANDELS completeness, which decreases toward lower M_* galaxies. Since more massive galaxies tend to have higher UV luminosities (Stark et al. 2009; González et al. 2014), our sample has a higher contribution of bright M_UV galaxies than a population-representative subsample.

In order to determine M_UV for our objects, we use the CANDELS i₇₇₅ band. We show in Figure 10 the corresponding distribution of the complete sample and detections. We also present the dependence of ${W}_{\mathrm{Ly}\alpha }$ on M_UV in this figure. As expected from the $\mathrm{SFR}\mbox{--}{W}_{\mathrm{Ly}\alpha }$ anti-correlations we recover in Section 5.2, a similar trend is observed for UV luminosities. This anti-correlation comes as no surprise, since brighter UV galaxies tend to have higher M_* at the cosmic time of our sample (Stark et al. 2009; González et al. 2014). Galaxies brighter in the UV have been subject to more intensive star formation events in 100 Myr timescales. Typically, higher neutral gas, turbulence, and bulk gas motions are associated with higher SFRs, boosting the scatter of Lyα photons. In combination, higher dust extinction, older age of stellar populations, and greater neutral gas mass in UV brighter galaxies seem to rule Lyα statistics dependence on SFR.

Zoom In Zoom Out Reset image size

Most analyses of Lyα emission dependence on M_UV have been performed using uncorrected SFRs. The literature compilation shown in Verhamme et al. (2008) reveals how observed UV SFRs anti-correlate with ${W}_{\mathrm{Ly}\alpha }$ from z ∼ 2.7 to 5.7 (Pettini et al. 2002; Shapley et al. 2003; Yamada et al. 2005; Gronwall et al. 2007; Tapken et al. 2007; Ouchi et al. 2008). Analyses explicitly using M_UV have also been performed in the surveys of Shimasaku et al. (2006), Ouchi et al. (2008), Vanzella et al. (2009), Balestra et al. (2010), Stark et al. (2010), Schaerer et al. (2011), and Cassata et al. (2011), yielding similar trends up to z ∼ 6. More recent studies confirm these trends (Jiang et al. 2013b, 2016; Zheng et al. 2014). Simulations likewise predict such correlations (see Shimizu et al. 2011). However, regardless of the methodology, the scatter in these correlations is non-negligible. Moreover, Atek et al. (2014) question the existence of the correlation in their sample. We argue that the scatter observed in the ${W}_{\mathrm{Ly}\alpha }$ dependence on M_UV/SFR_obs can be a consequence of the role played by dust. Galaxies brighter in the UV naturally have a greater Lyα photon production, but the anti-correlation with M_* affects the escape fraction (probably through increased dust extinction). In this scenario, Lyα photon escape is a complex process simultaneously ruled by different properties of high-redshift galaxies. We explore such a property space in our analysis of the ${W}_{\mathrm{Ly}\alpha }$ distribution dependence on the ${M}_{\mathrm{UV}}\mbox{--}\beta$ relationship in Section 5.5.

Implications involving the fraction of detections as a function of M_UV are not straightforward. In principle, this fraction seems to correlate with the UV luminosities of our galaxies, as opposed to what is observed for characteristic ${W}_{\mathrm{Ly}\alpha }$. Nevertheless, our Lyα measurements are flux limited. Therefore, our detection completeness in ${W}_{\mathrm{Ly}\alpha }$ is higher for brighter objects, leading to biases difficult to account for, as noted in Nilsson et al. (2009). Under these circumstances, the ideal approach is to consider both detections and non-detections, while taking into account the uncertainties for line and continuum fluxes (F and f_λ, respectively). Hence, we encourage further interpretations of these results to focus on the analysis performed on the ${M}_{\mathrm{UV}}\mbox{--}\beta$ plane (Section 5.5).

Measurement of the UV slopes of high-redshift galaxies is a direct way of tracing the amount of dust inside galaxies, given the assumption of an extinction law and an intrinsic spectral shape. This is of particular interest for Lyα surveys, since simulations (e.g., Verhamme et al. 2008) and observations at low redshift (Hayes et al. 2011; Atek et al. 2014) suggest that dust plays an important role in the escape of Lyα photons. Since we have rest-frame UV photometry from CANDELS for all our objects, we can determine their UV slopes and study their effect on Lyα emission at z ∼ 4. In this section, we detail our method to estimate the UV slopes for our objects and show our results on the Lyα dependence on this galaxy property.

To determine UV slopes, we fit a power law ${f}_{\lambda }={f}_{0}{\lambda }^{\beta }$ (Calzetti et al. 1994) to the photometry of each object. For fitting, we just use a standard least-squares routine on the photometry between rest frame 1400 and 3500 Å. These calculations correspond, in principle, to 8–15 bands between 1500 and 3600 Å. However, since we only use fluxes with S/N > 3, our median number of bands is eight. Naturally, we further require at least two bands to associate a slope to our targets, which means we can measure the UV slope for 611 of the 625 observed objects composing our sample.

Our FAST outputs include dust extinction in the V-band, A_V, for every object. However, if we want to associate a reddening $E(B-V)$ with every galaxy, we need to use our derived UV slopes. We use the relation

$\begin{eqnarray}&&{f}_{\lambda }={f}_{0}{\lambda }^{\beta }\propto {\lambda }^{-2.23}{10}^{0.4{k}_{\lambda }E(B-V)},\end{eqnarray} \tag{ 23 }$

where we assume a pristine slope of β = −2.23 (Meurer et al. 1999). We also require the adoption of an attenuation law ${k}_{\lambda }$ (e.g., Calzetti et al. 2000). Reddy et al. (2015) recover a more appropriate attenuation law than Calzetti et al. (2000) for high-redshift galaxies, so we use the ${k}_{\lambda }$ of the former. It is worth noting, nevertheless, that both attenuation laws yield almost identical results.

Since our targets are M_* selected, we have a higher contribution of massive objects in comparison to the M_* distribution of the galaxy population. As more massive objects tend to have higher $E(B-V)$ and redder UV slopes, we expect z ∼ 4 samples representative of the galaxy population to have a lower contribution from such galaxies. In any case, we show in Figure 11 the UV slope histogram of our sample and our Lyα emitters. Just by comparing both distributions, it is clear that the fraction of emitters increases toward bluer galaxies. We present a quantitative analysis of this claim in Section 5.5.

Zoom In Zoom Out Reset image size

We also show in Figure 11 our results on ${W}_{\mathrm{Ly}\alpha }$ and f_esc as a function of β and $E(B-V)$. Since lower mass galaxies have bluer UV slopes than more massive ones, the trends we find are complementary to our previous results. There is a correlation between the steepness of the UV spectrum and ${W}_{\mathrm{Ly}\alpha }$/f_esc, although with significant scatter. As extinction seems to play a major role in Lyα photon escape from galaxies, mainly through scattering and absorption (Blanc et al. 2011; Hagen et al. 2014), these correlations come as no surprise. Qualitatively, our results agree with the measurements from Shapley et al. (2003), Pentericci et al. (2009), Blanc et al. (2011), and Atek et al. (2014). Regarding the scatter we find at fixed β (see also Blanc et al. 2011), it is consistent with a scenario where the observed Lyα flux is mostly affected by the dependence of the dust-covering fraction on the line of sight. We discuss this picture in detail in Section 5.5.

When comparing the more dusty Lyα emitters in our sample with results from the literature, however, some differences show up. The 1.9 < z < 3.8 survey from Blanc et al. (2011) and the z < 0.5 study from Atek et al. (2014) find Lyα emitters up to $E(B-V)\sim 1$. Similarly, Matthee et al. (2016) find a population of dusty LAEs at z ∼ 2.23, and they speculate on how dusty gas outflows might be the feature driving the escape of Lyα radiation. Lyα sources up to z ∼ 3 are Herschel (Oteo et al. 2011, 2012; Casey et al. 2012; Sandberg et al. 2015) and SCUBA (Geach et al. 2005; Hine et al. 2016) detected. Moreover, Lyα emission has also been measured in submillimeter galaxies (e.g., Chapman et al. 2005). On the other hand, our 625 object sample features ∼5 objects with $E(B-V)\gt 0.4$, but none of them qualifies as a detection. We cannot state whether the absence of such LAEs in our sample is representative of z ∼ 4. True enough, the fraction of dusty galaxies at z ∼ 4 is expected to be lower than at z < 4. Along such lines, cosmic evolution in the fraction of dusty LAEs has already been discussed in the literature. Blanc et al. (2011) study any dust evolution in their 1.9 < z < 3.8 LAEs sample and find no significant trend. Hagen et al. (2014) study the same sample, and show that there is little anti-correlation, if any, between $E(B-V)$ and redshift. Shapley et al. (2003) do observe such evolution in the range 2 < z < 3.5, but question the validity given the selection effects associated with LBG surveys. We conclude that even though Lyα emitting galaxies are mostly low-dust objects (e.g., Song et al. 2014), there is also a population of dustier, low-${W}_{\mathrm{Ly}\alpha }$ LAEs. Still, the significance of their numbers at z ≳ 4 is still an open question.

Our results also confirm at z ∼ 4 a trend in f_esc that has already been observed at lower redshift. Atek et al. (2014) perform a f_esc study at z < 0.5 and recover a similar anti-correlation to the one we show in Figure 11. Hayes et al. (2011) at z ∼ 2.2, Blanc et al. (2011) at 1.9 < z < 3.8, Song et al. (2014) at z ∼ 2.1–2.5, and Matthee et al. (2016) at z ∼ 2.23 also observe the same trends. Verhamme et al. (2008) replicate such trends using radiative simulations of galaxies in the range 2.8 < z < 5, confirming that qualitative explanations for these observational relations are well supported by theory. We show several best-fit relations from the literature in our lower plot of Figure 11. Considering that our results are dominated by upper limits, the two higher redshift relations are roughly consistent with our measurements. Still, of particular interest might be the potential redshift evolution suggested by these relations. There is a clear decrease in f_esc at high dust contents when going from low to high redshift. If real, this trend could back our previous analysis on the fraction of more dusty LAEs and their evolution as a function of cosmic time. At higher redshift (Verhamme et al. 2008; this work), low-${W}_{\mathrm{Ly}\alpha }$, very dusty LAEs do not seem to be common, driving the $E(B-V)-{f}_{\mathrm{esc}}$ relation down significantly for $E(B-V)\gt 0.2$. However, at lower redshift, such objects are actually observed, driving the relation up for $E(B-V)\gt 0.2$.

We have characterized in this paper the dependence of Lyα emission on M_*, SFR, UV luminosity, and UV slope. We find ${W}_{\mathrm{Ly}\alpha }$ and f_esc to anti-correlate with these four properties. However, these results might not be independent, since these properties are correlated (see Figure 6). For instance, M_* and SFR are known to follow a relation at high redshift referred to as the main sequence (Kereš et al. 2005; Finlator et al. 2007; Noeske et al. 2007; Daddi et al. 2007). Similarly, a relation between M_UV and β has been studied in LBGs (Bouwens et al. 2009, 2014). As we discuss throughout this work, Lyα escape from galaxies is likely to be a process ruled by many parameters, such as age of the population, gas column density, extinction, and SFR. Therefore, it is interesting to explore Lyα emission dependence on multi-dimensional spaces. In this section, we study ${W}_{\mathrm{Ly}\alpha }$ dependence on M_UV and β. There are two reasons to justify exploring this parameter space in particular. First, both properties are observables, i.e., the amount of assumptions involved in their calculation is kept to a minimum (unlike, for example, M_* and SFR). Second, these two properties are directly related to elements that rule Lyα escape. The UV luminosity of a galaxy traces its ionizing output and gas content, whereas the UV slope is a proxy for dust. Hence, in this section, we focus on characterizing Lyα emission in the ${M}_{\mathrm{UV}}\mbox{--}\beta$ plane. Similarly to previous analyses in this paper, we take advantage of a Bayesian approach to obtain our results.

We first present in Figure 12 the location of our detections and non-detections in this plane. In qualitative terms, our detections sample most of the sequence initially covered by our targets. For further insight into how our observations depend on this sequence, we construct four subsamples for visualization. The anti-correlation we observe between ${W}_{\mathrm{Ly}\alpha }$ and M_UV is still clearly observed when comparing the low M_UV with the high M_UV subsamples, independent of β. For the UV slope, however, the trend between ${W}_{\mathrm{Ly}\alpha }$ and β does not seem that clear anymore. The normalization of the distribution, however, looks to be highly dependent on the UV slope. To characterize the significance of these insights, we simultaneously model the exponential profile parameters of Equation (11) as a function of M_UV and β. We start again with linear expressions of the form

$\begin{eqnarray}&&A({M}_{\mathrm{UV}},\beta )={A}_{{M}_{\mathrm{UV}}}{M}_{\mathrm{UV}}+{A}_{\beta }\beta +{A}_{C},\end{eqnarray} \tag{ 24 }$

$\begin{eqnarray}&&{W}_{0}({M}_{\mathrm{UV}},\beta )={W}_{{M}_{\mathrm{UV}}}{M}_{\mathrm{UV}}+{W}_{\beta }\beta +{W}_{C}.\end{eqnarray} \tag{ 25 }$

Once the parameterization and priors are set, we can obtain the posterior distribution of $\{{W}_{n}\}=({A}_{{M}_{\mathrm{UV}}},{A}_{\beta },{A}_{C},{W}_{{M}_{\mathrm{UV}}},{W}_{\beta },{W}_{C})$ using Equation (9). We assume the priors to be independent, which translates to $p(\{{W}_{n}\})\propto \prod p({W}_{i})$. Once again, we impose ignorance on the parameters, i.e., we adopt non-informative priors. We study the posterior distribution assuming parameter ignorance in both linear and logarithmic scales. We decide for the first, since logarithmic priors diverge for parameters than can adopt values close to zero (${A}_{{M}_{\mathrm{UV}}}$, A_C and ${W}_{\beta }$). Therefore, our prior is simply a constant $p(\{{W}_{n}\})$. Considering that $p(\{F\})$ is just a normalization factor in Equation (9), the posterior distribution of the six-parameter space can be obtained:

$\begin{eqnarray}&&p({W}_{n}| \{F\})\propto \displaystyle \prod _{D}p({F}_{i}| {W}_{n})\displaystyle \prod _{{ND}}p({F}_{i}\lt {F}_{i}^{* }| {W}_{n}).\end{eqnarray} \tag{ 26 }$

Zoom In Zoom Out Reset image size

The value of the parameters A and W₀ is, by definition, constrained to the intervals $A=[0,1]$ and ${W}_{0}=(0,\infty )$. Still, in the case of extreme luminosities and/or slopes, our linear parameterizations can yield values outside these intervals. As a solution, we just impose A and W₀ to saturate outside their corresponding ranges. This condition can be considered simply as an indirect prior on the parameters. For the particular case of ${W}_{0}\lt 0$, we impose A = 0, i.e., the ${W}_{\mathrm{Ly}\alpha }$ distribution is a Dirac delta (i.e., no Lyα emission and/or absorption).

We recover our best solution from Monte Carlo simulations and obtain the uncertainties on our parameters using the collapsed distributions. The results are as follows:

$\begin{eqnarray}&&A({M}_{\mathrm{UV}},\beta )={0.08}_{-0.06}^{+0.1}{M}_{\mathrm{UV}}-{0.6}_{-1.2}^{+0.1}\beta +{1.1}_{-1.4}^{+1.5},\end{eqnarray} \tag{ 27 }$

$\begin{eqnarray}&&{W}_{0}({M}_{\mathrm{UV}},\beta )={11}_{-1.8}^{+2.0}{M}_{\mathrm{UV}}-{7.8}_{-5.5}^{+5.2}\beta +{235}_{-42}^{+45}.\end{eqnarray} \tag{ 28 }$

We show with more detail the results for these coefficients in Figure 13. The early analysis based on the sample binning from Figure 12 is verified in Figure 13. The normalization factor A anti-correlates with β with a significance ≳2σ. There also seems to be an anti-correlation between A and M_UV, but to a level of ∼1σ. For W₀, the e-folding scale of the distribution, the behavior is the opposite. The extent of the distribution correlates with M_UV and shows a weak anti-correlation with β.

We now focus on the interpretation of these results. The fact that beta correlates mostly with A suggests that the UV slope is somehow related to a stochastic probability of being able to observe or not a particular object in Lyα. On the other hand, the stronger correlation of M_UV with ${W}_{\mathrm{Ly}\alpha }$ suggests that the UV luminosity is associated with physical processes that determine the magnitude of the resulting ${W}_{\mathrm{Ly}\alpha }$. A possible scenario is one in which M_UV traces SFR and, therefore, the total cold gas mass of galaxies. Hence, higher UV luminosities translate into increased levels of Lyα photon scattering into an extended halo, smoothly reducing the ${W}_{\mathrm{Ly}\alpha }$ of the central source. At the same time, in a significantly clumpy ISM in which dust is well mixed with the gas, a large contrast in gas density/column will imply that dust can effectively block the bulk of both UV and Lyα radiation in some parts of the disk. Such impact will be significant toward certain lines of sight, with little to no effect on others. Such a "covering factor" scenario could explain that β correlates more significantly with the normalization (A) than with the shape of the distribution (W₀). This is consistent with the finding that Lyα and UV suffer from similar levels of extinction by dust in LAEs (e.g., Blanc et al. 2011). Moreover, the scatter in the f_esc of Lyα photons at fixed reddening (see Figure 11 and Blanc et al. 2011) is also consistent with this scenario where wide ranges of dust absorption lines of sight and photon scattering halos rule the escape of Lyα radiation.

Naturally, this interpretation is a simplified description of Lyα escape from galaxies. Going no further, Figure 13 reveals how both distribution parameters do not solely depend on one property. Ideally, studies of ${W}_{\mathrm{Ly}\alpha }$ and f_esc in three-parameter spaces (e.g., M_*, M_UV, β) can yield predictive and more accurate parameterization of Lyα emission. However, such analysis must be performed in much larger data sets than ours, at least if significant enough results are to be obtained. Furthermore, hydrodynamical simulations can also give insight into how Lyα escape depends on the line of sight and the distribution of gas, dust, and star-forming regions (e.g., Verhamme et al. 2012), especially if performed in a statistically significant sample.

Zoom In Zoom Out Reset image size

The Lyman break selection technique has proven to be a very efficient method for detecting high-redshift galaxies (e.g., Steidel et al. 1996; Shapley et al. 2003; Stark et al. 2010; Ono et al. 2012). The fact that the Lyman break is in the optical region of the observed spectrum for galaxies at redshift z ≥ 3 allows for efficient detection from ground telescopes. By only requiring the use of a few broadband filters, several galaxies can be detected in a deep single exposure. Still, to avoid aliasing with the Balmer break, there are unavoidable biases associated with this technique. First, galaxies with no prominent Lyman break are, by construction, not detected. As a consequence, either extremely young or passive, heavily extincted galaxies are underrepresented in LBG surveys. Second, these surveys also impose color restrictions on the slope of galaxies, further increasing the selection toward, in principle, bluer UV objects. Third, this technique is limited by the M_UV sensitivity of the survey, creating incompleteness at low SFR. In the case of Lyα emission, this observational limit can have a significant downside. As shown in Section 5.2, there is a clear anti-correlation between UV luminosity and ${W}_{\mathrm{Ly}\alpha }$, leading to the possibility that Lyα studies in LBG samples are missing the highest ${W}_{\mathrm{Ly}\alpha }$ galaxies in the universe. Indeed, the deepest surveys nowadays reach ∼50% completeness at observed ${m}_{\mathrm{UV}}\sim 29.5$ (e.g., Bouwens et al. 2015; Bowler et al. 2015, 2017; Ishigaki et al. 2017). These limits translate to higher redshift ∼50% completeness at ${M}_{\mathrm{UV}}\sim -18.3$, −19, −19.5, and −19.7 (z ∼ 4, 5, 6, and 7, respectively). These limitations become problematic when comparing the results between LBG and narrowband samples. Presumably due to not requiring a continuum detection, the ${W}_{\mathrm{Ly}\alpha }$ observed in narrowband surveys are systematically higher (e.g., Zheng et al. 2014).

Our data set and results can be useful for characterizing the effects that selection techniques can have on Lyα emission. Even though it is possible to identify these selections directly in our sample, any comprehensive analysis must consider our galaxy selection procedure. Since our sample follows neither the M_* or M_UV functions of the galaxy population, correcting requires assumption of M_* and/or M_UV distributions. Nevertheless, we can still simulate Lyα emission samples on CANDELS galaxies using our results of ${W}_{\mathrm{Ly}\alpha }$ dependence on the ${M}_{\mathrm{UV}}\mbox{--}\beta$ plane given by Equations (27) and (28). In this section, we describe our simulation of a z ∼ 4 LBG sample from the 3D-HST catalogs. We then compare the properties of LBGs with the parent distribution, focusing on M_*, SFR, M_UV, and β. We finally simulate Lyα fractions for each sample and conclude on the effects of LBG selection at z ∼ 4.

LBGs at z ∼ 4 are typically selected using B-dropouts and imposing color selections on redder filters. As we are simulating this survey in CANDELS, our detection limits are given by the depth of their H₁₆₀ images. For the z ∼ 4 LBG selection, we adopt the same methodology applied in Bouwens et al. (2012):

$\begin{eqnarray}&&{\rm{B}} \mbox{-} \mathrm{dropouts}:\ ({B}_{435}-{V}_{606}\gt 1.1)\wedge ({B}_{435}-{V}_{606}\\ &&\quad \gt ({V}_{606}-{z}_{850})+1.1)\wedge ({V}_{606}-{z}_{850}\lt 1.6).\end{eqnarray} \tag{ 29 }$

After color selection, we impose 5σ detections in the V₆₀₆ filter, to which we now refer to as the detection band. As commonly applied in LBG surveys, we also require all candidates to have <1.5σ measurements in the U band, since the Lyman break must be redshifted out of this filter. Before performing any analysis on the CANDELS LBG sample, we run EAZY and FAST according to our prescriptions (see Section 2.3). Of the ∼3200 CANDELS galaxies that comply with our selection, ∼800 classify as low-redshift interlopers (${z}_{\mathrm{EAZY}}\lt 3$), according to our outputs. We check EAZY 3σ constraints for these interlopers and find most to have ${z}_{99}\lt 4$ . Hence, we find a low-redshift contamination in z ∼ 4 CANDELS LBGs that is higher than the typically reported number of ∼10% (Bouwens et al. 2007). From now on, we remove these presumed ${z}_{\mathrm{EAZY}}\lt 3$ contaminants and work with a z ∼ 4 LBG sample of nearly 2300 objects.

Using our outputs, we plot relevant properties of CANDELS LBGs in Figure 14. To ensure that we do not venture way beyond CANDELS completeness, we restrict the upcoming analysis to galaxies with ${M}_{* }\gt {10}^{8.5}\,{M}_{\odot }$. We construct two $3\lt z\lt 4.2$ CANDELS samples for comparison: all galaxies with ${M}_{* }\gt {10}^{8.5}\,{M}_{\odot }$ and all V-band detected ${M}_{* }\gt {10}^{8.5}\,{M}_{\odot }$ objects. As seen in the central panel, the major LBG selection effect is associated with the detection band threshold. These left-out objects are either red, heavily extincted galaxies, or blue, intrinsically faint objects (Quadri et al. 2007). Indeed, as revealed by the comparison between LBGs and the detection-band-selected sample, any selection imprints associated with color requirements are minor. It is still worth noting, however, that these color criteria seem to neglect blue rather than red galaxies at z ∼ 4 (center of Figure 14). We find the $(B-V\gt 1.1)$ Lyman break cut to be the driving criterion behind this selection effect (see Equation (29)). We remark that these insights might not be valid at different redshifts or lower M_*, however. Some studies have indeed explored differences at higher redshift (e.g., Jiang et al. 2013a, 2013b, 2016). In analyses that explore Lyα emission strength, UV continuum properties, morphologies, and sizes, they find no significant differences between LBGs and LAEs at z ≥ 6.

Zoom In Zoom Out Reset image size

The insights we present can be further tested by simulating the fraction of LAEs above an ${W}_{\mathrm{Ly}\alpha }$ threshold. We use our results from Section 5.5 to estimate the Lyα fraction for each sample and plot our results in the right panel of Figure 14. Indeed, depending on the depth of the detection band, the fraction of low UV luminosity galaxies composing the sample can vary significantly. As a consequence, deeper surveys can potentially recover a higher fraction of line emitters at high ${W}_{\mathrm{Ly}\alpha }$. To deal with this selection, it has been proposed to compare the Lyα fractions over narrow UV luminosity samples (Stark et al. 2011; Mallery et al. 2012; Ono et al. 2012; Schenker et al. 2012). Our results in this paper further emphasize the need to compare galaxies of the same luminosity and properly characterize completeness when studying the evolution of the Lyα fraction with cosmic time. Regarding color selections, we find the Lyα fraction to be slightly lower for LBGs than for the V-band-detected sample. As we suggested, this is a consequence of LBG color cuts leaving out some of the bluest galaxies.

Lyα-emitting galaxies can also be selected using narrowband imaging or blind spectroscopy (e.g., Gronwall et al. 2007; Ouchi et al. 2008; Adams et al. 2011). By establishing a ${W}_{\mathrm{Ly}\alpha }$ detection threshold between the narrow- and broadband flux measurement, this technique allows for efficient line-emitter selection. Since only line detection is required, such surveys can trace fainter objects than the LBG technique, possibly leading to selection of the youngest and faintest galaxies at high redshift. However, even though most narrowband measurements of LAEs are followed up by spectroscopy, any sample selection effects induced by the narrowband technique are already present in the sample. The ${W}_{\mathrm{Ly}\alpha }$ threshold used for detection, which is determined by the ability to separate low-redshift interlopers from high-redshift LAEs, can adopt a wide range of values. Depending on redshift, observer frame thresholds translate into different rest-frame ${W}_{\mathrm{Ly}\alpha }$ cuts. For instance, Vargas et al. (2014) select sources with ${W}_{\mathrm{Ly}\alpha }\gt 20$ Å at $z\sim 2.1$, Zheng et al. (2014) select ${W}_{\mathrm{Ly}\alpha }\gt 9$ Å at z ∼ 4.5, and Sobral et al. (2017) go as low as ${W}_{\mathrm{Ly}\alpha }\gt 5$ Å at z ∼ 2.23. If ${W}_{\mathrm{Ly}\alpha }$ selections induce important biases on galaxy samples, the comparison of different surveys is not straightforward.

In this section, we explore the effects that narrowband selections have on the population of LAEs, focusing on the M_*–SFR plane. The insights we present here are based on Equations (27) and (28), i.e., our ${M}_{\mathrm{UV}}\mbox{--}\beta$ model. We show in Figure 15 the outcome of ${W}_{\mathrm{Ly}\alpha }$ (top) and line-flux selections (bottom) on the M_*–SFR sequence. The red contours show 3.5 < z < 4.5 CANDELS galaxies with $M\gt {10}^{8.5}\,{M}_{\odot }$. For clarity, however, we concentrate our analysis on the main sequence, removing objects with young and old ages (red dots). This approach allows our results to be dominated by objects optimally fitted by our FAST executions.

Zoom In Zoom Out Reset image size

To explore the effect of ${W}_{\mathrm{Ly}\alpha }$ selections (top panel of Figure 15), we base our analysis on the z ∼ 2.23 narrowband survey of Sobral et al. (2017), whose galaxies are limited to ${W}_{\mathrm{Ly}\alpha }\gt 5$ Å and fluxes $\gt 2\times {10}^{-17}$ erg cm⁻² s⁻¹ (5σ). We simulate a flux-limited-only Lyα sample and use it as baseline (black). We then show the region where ${W}_{\mathrm{Ly}\alpha }\,\gt 5$ Å (1σ; blue) and ${W}_{\mathrm{Ly}\alpha }\gt 20$ Å (1σ; green) objects lie. The ${W}_{\mathrm{Ly}\alpha }\gt 5$ Å contours are intended to reproduce Sobral et al. (2017) selections, whereas the more restrictive selection ${W}_{\mathrm{Ly}\alpha }\gt 20$ Å is representative of multiple surveys (Gronwall et al. 2007; Hagen et al. 2014; Vargas et al. 2014). As confirmed by this plot, galaxies with higher characteristic ${W}_{\mathrm{Ly}\alpha }$ (i.e., low M_*, low SFR) are more likely to be selected by narrowband samples, even though the effect is minor for the typical narrowband cuts of ${W}_{\mathrm{Ly}\alpha }\leqslant 20$ Å. Still, selections based solely on ${W}_{\mathrm{Ly}\alpha }$ systematically fail to remove AGNs and neglect bright Lyα emitters (Sobral et al. 2017). All these insights remark on the importance of using low ${W}_{\mathrm{Ly}\alpha }$ cuts and thorough interloper controls.

After narrowband outcomes are used for sample selection, spectroscopic observations follow. However, these follow-ups of Lyα emitting galaxies are flux limited, just like the ones presented in this work. Still, as we consider completeness in our modeling and simulations of Section 5.5, we can still make inferences on flux limited studies. We use our Monte Carlo simulation outputs to also assess the effects of line luminosity selections. Using the corresponding M_*, SFR, and f_UV for every object, we obtain the probability of ${L}_{\mathrm{Ly}\alpha }\gt {L}_{\mathrm{Ly}\alpha }^{* }$ for every galaxy. The outcome (bottom panel of Figure 15) reveals that flux selections bias samples toward high-SFR LAEs.

Discrepancies in the location of LAEs in the M_*–SFR plane have already been observed in the literature. In Figure 10 of their paper, Hagen et al. (2014) plot the M_*–SFR relation of their z ∼ 2–3 LAEs alongside the z ∼ 2.1 counterparts of Vargas et al. (2014). The Hagen et al. (2014) objects are part of the HETDEX Pilot Survey (Adams et al. 2011), for which detections are constrained to line fluxes $\gt 7\mbox{--}10\,\times {10}^{-17}$ erg cm⁻² s⁻¹ (5σ) and ${W}_{\mathrm{Ly}\alpha }\gt 20$ Å. In contrast, the Vargas et al. (2014) flux depth is higher than that of HETDEX, reaching $2\times {10}^{-17}$ erg cm⁻² s⁻¹ (5σ), while also selecting sources with ${W}_{\mathrm{Ly}\alpha }\gt 20\,\mathring{\rm A}$. Therefore, Hagen et al. (2014) are comparing two ${W}_{\mathrm{Ly}\alpha }$-selected surveys with different line-flux depths. Our results in this section would, then, predict Vargas et al. (2014) LAEs to sample lower M_* and SFRs because they go deeper (see $\sim {10}^{9}\,{M}_{\odot }$ galaxies in Figure 15). This is in fact the pattern observed in Figure 10 of Hagen et al. (2014), i.e., we can qualitatively reproduce their comparison with our simulations from the results of Section 5.5. This explanation is backed by the fact that these discrepancies are not caused by inconsistent M_* or SFR derivations. Both studies assume a Salpeter (1955) IMF, cSFHs, and a fixed metallicity of $Z=0.2\,{Z}_{\odot }$.

Vargas et al. (2014) and Hagen et al. (2014) find ${M}_{* }\lt {10}^{9}\,{M}_{\odot }$ LAEs to lie above the M_*–SFR relation. A similar trend is reported in the work of Karman et al. (2017), who study Lyα emitters down to 10⁶ M_⊙. Karman et al. (2017) point out that this offset can be a consequence of how uncertain SFRs are for low-M_* starburst galaxies. This explanation is consistent with the fact that most of our ${M}_{* }\lt {10}^{9}\,{M}_{\odot }$ Lyα emitters are constrained to the lowest age locus of the M_*–SFR plane (see Figure 6). Moreover, Finkelstein et al. (2015) observe the same lowest age, low -M_* Lyα emitters at $z\sim 4.5$. To answer whether low -M_* LAEs lie above star-forming galaxies or the slope of the sequence changes toward ${M}_{* }\lt {10}^{9}\,{M}_{\odot }$, different approaches are required. For instance, no discrepancies are found between z ∼ 2 LAEs and Hα sources (see Matthee et al. 2016), although most of their galaxies have ${M}_{* }\gt {10}^{9}\,{M}_{\odot }$. In summary, evidence suggests that z > 2, ${M}_{* }\lt {10}^{9}\,{M}_{\odot }$ LAEs lie above extrapolations of the M_*–SFR relation. It is unclear whether this offset is real or the position of the relation at low M_* is far from certain.