A Comprehensive Study of Ly$\alpha$ Emission in the High-redshift Galaxy   Population

Grecco A. Oyarz\'un; Guillermo A. Blanc; Valentino Gonz\'alez; Mario; Mateo; John I. Bailey III

arXiv:1706.01886·astro-ph.GA·July 17, 2017

A Comprehensive Study of Ly$\alpha$ Emission in the High-redshift Galaxy Population

Grecco A. Oyarz\'un, Guillermo A. Blanc, Valentino Gonz\'alez, Mario, Mateo, John I. Bailey III

PDF

TL;DR

This study provides a detailed analysis of Ly$ extalpha$ emission in high-redshift galaxies, revealing how observational biases and survey limits influence the detection and interpretation of Ly$ extalpha$ properties and their evolution.

Contribution

It introduces a Bayesian model linking Ly$ extalpha$ equivalent widths to galaxy properties, accounting for survey selection effects and informing reionization studies.

Findings

01

Ly$ extalpha$ emission anti-correlates with galaxy mass and SFR.

02

Survey limits significantly bias Ly$ extalpha$ detection towards low-mass, low-SFR galaxies.

03

The observed drop in Ly$ extalpha$ fraction at high redshift can be explained by survey incompleteness.

Abstract

We present an exhaustive census of Lyman alpha (Ly $α$ ) emission in the general galaxy population at $3 < z < 4.6$ . We use the Michigan/Magellan Fiber System (M2FS) spectrograph to study a stellar mass (M $_{*}$ ) selected sample of 625 galaxies homogeneously distributed in the range $7.6<\log{\mbox{M$ _* $/ M$ _{\odot} $}}<10.6$ . Our sample is selected from the 3D-HST/CANDELS survey, which provides the complementary data to estimate Ly $α$ equivalent widths ( $W_{L y α}$ ) and escape fractions ( $f_{esc}$ ) for our galaxies. We find both quantities to anti-correlate with M $_{*}$ , star-formation rate (SFR), UV luminosity, and UV slope ( $β$ ). We then model the $W_{L y α}$ distribution as a function of M $_{U V}$ and $β$ using a Bayesian approach. Based on our model and matching the properties of typical Lyman break galaxy (LBG) selections, we conclude that the $W_{L y α}$ …

Figures17

Click any figure to enlarge with its caption.

Tables1

Table 1. Table 1 Ly α 𝛼 \alpha Equivalent Width distribution model comparison

Likelihood	Low Mass			Medium Mass			High Mass			Complete Sample
Output	Exp	Gaussian	Log n	Exp	Gaussian	Log n	Exp	Gaussian	Log n	Exp	Gaussian	Log n
Model odds^a^aObtained by integrating the likelihood over the whole parameter space.	0.93	0.06	0.01	0.97	0.01	0.02	1	$10^{- 3}$	$10^{- 5}$	0.98	$10^{- 6}$	0.02
Peak odds^b^bCalculated with the maximum of the likelihood.	0.25	0.02	0.73	0.23	$10^{- 3}$	0.77	0.97	$10^{- 3}$	0.03	0.19	$10^{- 7}$	0.81
$A$ ^c^cModel parameters for the corresponding likelihood maximum. Note that these values are different from the ones in Oyarzún et al. (2016), since here we are just working with the likelihood.	0.75	0.55	0.55	0.55	0.4	0.5	0.35	0.2	0.2	0.4	0.3	0.45
$W_{0}$ ^c^cModel parameters for the corresponding likelihood maximum. Note that these values are different from the ones in Oyarzún et al. (2016), since here we are just working with the likelihood.	46	74	0.7	26	46	0.85	14	26	0.75	38	64	1.05
$μ$ ^c^cModel parameters for the corresponding likelihood maximum. Note that these values are different from the ones in Oyarzún et al. (2016), since here we are just working with the likelihood.	…	…	3.85	…	…	3.05	…	…	3	…	…	3.05

Equations78

f_{r es t} (λ) = \frac{F}{2 π σ _{λ}^{2}} \mbox e x p [- \frac{( λ - λ _{0} ) ^{2}}{2 σ _{λ}^{2}}]

f_{r es t} (λ) = \frac{F}{2 π σ _{λ}^{2}} \mbox e x p [- \frac{( λ - λ _{0} ) ^{2}}{2 σ _{λ}^{2}}]

f_{r es t} (λ) = f_{b l u e} (λ) + f_{r e d} (λ)

f_{r es t} (λ) = f_{b l u e} (λ) + f_{r e d} (λ)

σ_{λ} \equiv σ_{b} \mbox i f λ < λ_{0}

σ_{λ} \equiv σ_{b} \mbox i f λ < λ_{0}

σ_{λ} \equiv σ_{r} \mbox i f λ > λ_{0}

W_{L y α} = \frac{F}{f _{λ}} \frac{1}{( 1 + z _{L y α} )}

W_{L y α} = \frac{F}{f _{λ}} \frac{1}{( 1 + z _{L y α} )}

f_{esc} (L y α) = \frac{L ( L y α ) _{o b s}}{L ( L y α ) _{in t r in s i c}} =

f_{esc} (L y α) = \frac{L ( L y α ) _{o b s}}{L ( L y α ) _{in t r in s i c}} =

\frac{\mbox S F R ( L y α )}{\mbox S F R ( \mbox U V ) _{cor r}} = \frac{4.4 \times 1 0 ^{- 42} L _{L y α} /8.7}{\mbox S F R ( \mbox U V ) _{cor r}}

p (W_{n} ∣ {F}) = \frac{p ({ F } ∣ W _{n} ) p ( W _{n} )}{p ({ F })}

p (W_{n} ∣ {F}) = \frac{p ({ F } ∣ W _{n} ) p ( W _{n} )}{p ({ F })}

p (F_{i} ∣ W_{n}) = \int_{0}^{\infty} p (F_{i} ∣ F) p (F ∣ W_{n}) d F

p (F_{i} ∣ W_{n}) = \int_{0}^{\infty} p (F_{i} ∣ F) p (F ∣ W_{n}) d F

p (F_{i} < F_{i}^{*} ∣ W_{n}) =

p (F_{i} < F_{i}^{*} ∣ W_{n}) =

\int_{0}^{\infty} (1 - p (F_{i} > F_{i}^{*} ∣ F)) p (F ∣ W_{n}) d F

p (W_{n} ∣ {F}) =

p (W_{n} ∣ {F}) =

\frac{p ( W _{n} )}{p ({ F })} D \prod p (F_{i} ∣ W_{n}) N D \prod p (F_{i} < F_{i}^{*} ∣ W_{n})

p (W_{L y α} ∣ W_{n}) = A \times H (W_{L y α}) \times p_{0} (W_{L y α} ∣ W_{n})

p (W_{L y α} ∣ W_{n}) = A \times H (W_{L y α}) \times p_{0} (W_{L y α} ∣ W_{n})

+ (1 - A) \times δ (W_{L y α})

p (W_{L y α} ∣ \mbox e x p) = \frac{A}{W _{0}} e^{- W_{L y α} / W_{0}} H (W_{L y α})

p (W_{L y α} ∣ \mbox e x p) = \frac{A}{W _{0}} e^{- W_{L y α} / W_{0}} H (W_{L y α})

+ (1 - A) δ (W_{L y α})

p (W_{L y α} ∣ \mbox G a u ss) = \frac{2 A}{2 π W _{0}^{2}} e^{- W_{L y α}^{2} /2 W_{0}^{2}} H (W_{L y α})

p (W_{L y α} ∣ \mbox G a u ss) = \frac{2 A}{2 π W _{0}^{2}} e^{- W_{L y α}^{2} /2 W_{0}^{2}} H (W_{L y α})

+ (1 - A) δ (W_{L y α})

p (W_{L y α} ∣ \mbox l o g n) = \frac{A}{2 π W _{0}^{2} W _{L y α}^{2}} e^{- (l n (W_{L y α}) - μ)^{2} /2 W_{0}^{2}}

p (W_{L y α} ∣ \mbox l o g n) = \frac{A}{2 π W _{0}^{2} W _{L y α}^{2}} e^{- (l n (W_{L y α}) - μ)^{2} /2 W_{0}^{2}}

+ (1 - A) δ (W_{L y α})

p (M_{i} ∣ {F}) = \frac{p ({ F } ∣ M _{i} ) p ( M _{i} )}{p ({ F })}

p (M_{i} ∣ {F}) = \frac{p ({ F } ∣ M _{i} ) p ( M _{i} )}{p ({ F })}

p (M_{i} ∣ {F}) \propto p ({F} ∣ M_{i}) = \int_{n} p (W_{n}) p ({F} ∣ W_{n}) d W_{n}

p (M_{i} ∣ {F}) \propto p ({F} ∣ M_{i}) = \int_{n} p (W_{n}) p ({F} ∣ W_{n}) d W_{n}

p (M_{i} ∣ {F}) = \frac{\int _{n} p ( W _{n} ) p ({ F } ∣ W _{n} ) d W _{n}}{\sum p ( M _{i} ∣ { F })}

p (M_{i} ∣ {F}) = \frac{\int _{n} p ( W _{n} ) p ({ F } ∣ W _{n} ) d W _{n}}{\sum p ( M _{i} ∣ { F })}

p (M_{i} ∣ {F}) = \frac{\int _{n} p ({ F } ∣ W _{n} ) d W _{n}}{\sum p ( M _{i} ∣ { F })}

p (M_{i} ∣ {F}) = \frac{\int _{n} p ({ F } ∣ W _{n} ) d W _{n}}{\sum p ( M _{i} ∣ { F })}

A (\mbox M_{*}) = A_{\mbox M_{*}} lo g (\mbox M_{*} [\mbox M_{⊙}]) + A_{C}

A (\mbox M_{*}) = A_{\mbox M_{*}} lo g (\mbox M_{*} [\mbox M_{⊙}]) + A_{C}

W_{0} (\mbox M_{*}) = W_{\mbox M_{*}} lo g (\mbox M_{*} [\mbox M_{⊙}]) + W_{C}

p (A_{\mbox M_{*}}, A_{C}, W_{\mbox M_{*}}, W_{C} ∣ {F}) = C \times

p (A_{\mbox M_{*}}, A_{C}, W_{\mbox M_{*}}, W_{C} ∣ {F}) = C \times

D \prod p (F_{i} ∣ A (\mbox M_{*}), W_{0} (\mbox M_{*})) N D \prod p (F_{i} < F_{i}^{*} ∣ A (\mbox M_{*}), W_{0} (\mbox M_{*}))

A (\mbox M_{*}) = - 0.2 8_{- .02}^{+ .12} lo g (\mbox M_{*} [\mbox M_{⊙}]) + 3. 1_{- 1.1}^{+ 0.15}

A (\mbox M_{*}) = - 0.2 8_{- .02}^{+ .12} lo g (\mbox M_{*} [\mbox M_{⊙}]) + 3. 1_{- 1.1}^{+ 0.15}

W_{0} (\mbox M_{*}) = - 17. 7_{- 2.0}^{+ 2.6} lo g (\mbox M_{*} [\mbox M_{⊙}]) + 19 0_{- 25}^{+ 20}

f_{λ} = f_{0} λ^{β} \propto λ^{- 2.23} 1 0^{0.4 k_{λ} E (B - V)}

f_{λ} = f_{0} λ^{β} \propto λ^{- 2.23} 1 0^{0.4 k_{λ} E (B - V)}

A (M_{U V}, β) = A_{M_{U V}} M_{U V} + A_{β} β + A_{C}

A (M_{U V}, β) = A_{M_{U V}} M_{U V} + A_{β} β + A_{C}

W_{0} (M_{U V}, β) = W_{M_{U V}} M_{U V} + W_{β} β + W_{C}

p (W_{n} ∣ {F}) \propto D \prod p (F_{i} ∣ W_{n}) N D \prod p (F_{i} < F_{i}^{*} ∣ W_{n})

p (W_{n} ∣ {F}) \propto D \prod p (F_{i} ∣ W_{n}) N D \prod p (F_{i} < F_{i}^{*} ∣ W_{n})

A (M_{U V}, β) =

A (M_{U V}, β) =

Peer Reviews

No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.

Videos

No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.

Full text

A Comprehensive Study of Ly $\alpha$ Emission in the High-redshift Galaxy Population

Grecco A. Oyarzún11affiliationmark: 22affiliationmark: , Guillermo A. Blanc11affiliationmark: 33affiliationmark: 44affiliationmark: , Valentino González11affiliationmark: 44affiliationmark: 55affiliationmark: , Mario Mateo66affiliationmark: , and John I. Bailey III66affiliationmark:

1Departamento de Astronomía, Universidad de Chile, Casilla 36-D, Santiago, Chile; [email protected]

2Astronomy Department, University of California, Santa Cruz, CA 95064, USA

3The Observatories of the Carnegie Institution for Science, 813 Santa Barbara Street, Pasadena, CA 91101, USA

4Centro de Astrofísica y Tecnologías Afines (CATA), Camino del Observatorio 1515, Las Condes, Santiago, Chile

5Chinese Academy of Sciences South America Center for Astronomy, China-Chile Joint Center for Astronomy,

Camino del Observatorio 1515, Las Condes, Santiago, Chile

6Department of Astronomy, University of Michigan, Ann Arbor, MI 48109, USA

Abstract

We present an exhaustive census of Lyman alpha (Ly $\alpha$ ) emission in the general galaxy population at $3<z<4.6$ . We use the Michigan/Magellan Fiber System (M2FS) spectrograph to study a stellar mass (M∗) selected sample of 625 galaxies homogeneously distributed in the range $7.6<\log{\mbox{M$ {}{*} $/M$ {}{\odot} $}}<10.6$ . Our sample is selected from the 3D-HST/CANDELS survey, which provides the complementary data to estimate Ly $\alpha$ equivalent widths ( $W_{Ly\alpha}$ ) and escape fractions ( $f_{esc}$ ) for our galaxies. We find both quantities to anti-correlate with M∗, star-formation rate (SFR), UV luminosity, and UV slope ( $\beta$ ). We then model the $W_{Ly\alpha}$ distribution as a function of MUV and $\beta$ using a Bayesian approach. Based on our model and matching the properties of typical Lyman break galaxy (LBG) selections, we conclude that the $W_{Ly\alpha}$ distribution in such samples is heavily dependent on the limiting MUV of the survey. Regarding narrowband surveys, we find their $W_{Ly\alpha}$ selections to bias samples toward low M∗, while their line-flux limitations preferentially leave out low-SFR galaxies. We can also use our model to predict the fraction of Ly $\alpha$ -emitting LBGs at $4\leqslant z\leqslant 7$ . We show that reported drops in the Ly $\alpha$ fraction at $z\geqslant 6$ , usually attributed to the rapidly increasing neutral gas fraction of the universe, can also be explained by survey MUV incompleteness. This result does not dismiss reionization occurring at $z\sim 7$ , but highlights that current data is not inconsistent with this process taking place at $z>7$ .

Subject headings:

dark ages, reionization, first stars - galaxies: evolution - galaxies: formation - galaxies: high-redshift - galaxies: statistics - ISM: lines and bands

††slugcomment: Accepted for publication in ApJ

1. Introduction

In hand with observational progress, the understanding of high redshift galaxies has progressed immensely in the last few decades. Almost 20 years ago, the only efficient way of observing these galaxies was by detecting breaks in their spectra with broadband photometry. For example, the two main classes of rest-frame UV-selected galaxies are Lyman Break Galaxies (LBGs; Steidel et al. 1996, 1999, 2003; Shapley et al. 2003; Stark et al. 2009; Bouwens et al. 2015) and Lyman-alpha Emitters (LAEs; Cowie et al. 1996; Cowie & Hu 1998; Hu et al. 1998; Kashikawa et al. 2006; Shimasaku et al. 2006; Gronwall et al. 2007; Ouchi et al. 2008; Blanc et al. 2011; Sobral et al. 2015). While the former are observed by their Lyman break at 912 Å, the latter are detected by their Lyman alpha (Ly $\alpha$ ) emission at 1216 Å, which is produced in star forming regions and is subject to resonant scattering in the neutral hydrogen of the ISM. After being, at first, used as a galaxy detection method, Ly $\alpha$ emission from galaxies is now widely used to derive a wide variety of information about the physical properties of galaxies and the universe as a whole.

A lot of effort has been dedicated to study the process of Ly $\alpha$ emission within galaxies. Radiative transfer analysis focuses on the dependence of Ly $\alpha$ escape on ISM kinematics, clumpiness, and outflows (Dijkstra & Kramer 2012; Duval et al. 2014; Rivera-Thorsen et al. 2015; Gronke & Dijkstra 2016; Gronke et al. 2016). These studies yield elementary insights that, when combined with observations of the Ly $\alpha$ escape fraction ( $f_{esc}$ ; Hayes et al. 2011; Ciardullo et al. 2014), provide a coherent picture for Ly $\alpha$ emission. Still, despite the complex modeling required to explain this particular type of radiation at small scales, Ly $\alpha$ is actively used as a technique for solving cosmological scale paradigms. Among other things, Ly $\alpha$ emission is being used to study the properties of neutral hydrogen in the interstellar, circumgalactic and intergalactic mediums. For instance, the fraction of LAEs as a function of redshift is used as a proxy for the fraction of neutral gas in the IGM (Stark et al. 2011; Ono et al. 2012; Tilvi et al. 2014). If successful, this insight can be used to trace the epoch of reionization, which we have been able to constrain using QSO sightlines (Fan et al. 2006) and cosmic microwave background (CMB) measurements (Bennett et al. 2013; Planck Collaboration et al. 2016a, b). Moreover, with the hope of tracing the sources responsible for reionization, Ly $\alpha$ emission is also conveniently used for characterizing luminosity functions (LFs) at high redshift (Ouchi et al. 2008; Sobral et al. 2009; Dressler et al. 2015; Santos et al. 2016). Nevertheless, the cosmological questions that we are trying to answer by means of this complex emission are not restricted to the reionization of the universe. For instance, Ly $\alpha$ emission can be used to trace the cosmic web (Gould & Weinberg 1996; Cantalupo et al. 2005, 2012, 2014; Kollmeier et al. 2010) and young, metal-poor stellar populations in the early universe (Sobral et al. 2015). Even more, Ly $\alpha$ emission will also be used to study dark energy through baryonic acoustic oscillations in the clustering of LAEs (Adams et al. 2011; Blanc et al. 2011).

As evidenced by small-scale studies, Ly $\alpha$ emission is highly stochastic and complex. The observed equivalent width ( $W_{Ly\alpha}$ ) of Ly $\alpha$ lines is highly dependent on neutral gas opacity and kinematics, clumpiness of the ISM, dust distribution, and, therefore, line of sight toward the observer. Furthermore, typical Ly $\alpha$ emitters are extreme objects in terms of their luminosity, metallicity, and emission line diagnostics (Trainor et al. 2016). As a consequence, if we are to use Ly $\alpha$ emission as a tracer, a careful and thorough approach is required. However, limitations in survey design make such an approach difficult. For instance, spectroscopic studies of UV continuum detected galaxies, typically LBGs (e.g. Shapley et al. 2003; Stark et al. 2010; Ono et al. 2012), reject galaxies with no significant Lyman break and/or red rest-frame UV spectral energy distributions (SEDs), i.e., these studies do not account for passive and heavily extincted galaxies. Observationally, this technique is limited by the MUV sensitivity of the survey, which translates to incompleteness at low-SFR. On the other hand, samples where galaxies are detected directly in Ly $\alpha$ using narrowband imaging (e.g. Gronwall et al. 2007; Ouchi et al. 2008) or blind spectroscopy (e.g. Blanc et al. 2011; Dressler et al. 2011) select on Ly $\alpha$ equivalent width (WLyα) and emission line flux. Such methodology does not require continuum selection, which allows for line detection in faint objects. However, follow-up rarely allows for more than spectroscopic confirmation of the line, limiting the use of this technique to the measurement of line fluxes, $W_{Ly\alpha}$ , and emission profiles. Even more, such surveys can have non-negligible amounts of line interlopers (Sobral et al. 2017).

Comprehensive studies of Ly $\alpha$ emission in the overall galaxy population are, therefore, of crucial importance. Correlations between galaxy properties and emission are the statistical manifestation of the radiative escape of Ly $\alpha$ photons. Properties such as stellar mass (M∗), star formation rate (SFR), UV slope ( $\beta$ ), and merger rate reveal information of how the escape of Ly $\alpha$ radiation is affected by stellar population ages, gas fraction, dust content, and ISM turbulence, respectively. Acknowledgment of such relations is required to draw conclusions from high-redshift Ly $\alpha$ surveys. Oyarzún et al. (2016) show how the normalization and e-folding scale of the WLyα distribution of $3<z<4.6$ galaxies anti-correlate with M∗. In other words, higher M∗ objects typically have lower WLyα. More massive galaxies typically have lower gas fractions, but higher gas mass (Tacconi et al. 2013; Song et al. 2014; Genzel et al. 2015). More neutral gas contributes to increase the scatter of Ly $\alpha$ photons, spreading this radiation throughout the galaxy and toward the circumgalactic medium. Higher M∗ galaxies also have more dust extinction (Franx et al. 2003; Daddi et al. 2004, 2007; Förster Schreiber et al. 2004; Pérez-González et al. 2008), which leads to more Ly $\alpha$ photon absorption. The escape fraction of this radiation is, therefore, severely affected by the stellar mass of galaxies. These results even explain some discrepancies between LBG and narrowband studies of Ly $\alpha$ emission that originated in sample selection effects. This example highlights the need for proper assessment.

In this work, we present a comprehensive analysis of Ly $\alpha$ emission in the $3<z<4.6$ galaxy population. To this end, we further analyze the data introduced by Oyarzún et al. (2016). This sample is composed of 625 galaxies in the M∗ range $7.6<\log{(M_{*}[M_{\odot}])}<10.6$ from the 3D-HST/CANDELS survey (Grogin et al. 2011; Koekemoer et al. 2011). Our work is based on spectroscopic observations of the sample using the Michigan/Magellan Fiber System (M2FS; Mateo et al. 2012), allowing us to measure Ly $\alpha$ fluxes and use 3D-HST/CANDELS ancillary data. In particular, we study Ly $\alpha$ emission dependence on the M∗, SFR, UV luminosity, and dust extinction of galaxies. To do so, we introduce a Bayesian approach to properly compare $W_{Ly\alpha}$ distribution models, account for incompleteness, and quantify the significance of observed correlations. We also use our results to simulate high-redshift Ly $\alpha$ emission surveys, allowing us to infer their biases and imprints on the resulting samples. We finally use the correlations we recover to predict the fraction of LAEs as a function of redshift. This semi-analytic approach provides a baseline for comparing observational drops in the fraction of LAEs at high redshift. This work is structured as follows. In Section 2, we describe our sample and data set. In Section 3, we describe our line detection and measurement methodologies. We explain our Bayesian $W_{Ly\alpha}$ distribution analysis in Section 4. We show our results on the $W_{Ly\alpha}$ distribution dependence on different properties and selection techniques in Sections 5 and 6. We state in Section 7 our inferences on higher redshift Ly $\alpha$ emission. We present our conclusions in Section 8. Throughout this work, all magnitudes are in the AB system (Oke & Gunn 1983). A CDM cosmology with $H_{0}$ = 70 km s*-1* Mpc*-1*, $\Omega_{m}$ = 0.3, and $\Omega_{\Lambda}$ = 0.7 was assumed whenever needed.

2. Data Set

2.1. Sample Selection

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 $\mu 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<z_{3D-HST}<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<\log({\mbox{M}_{*}[\mbox{M}_{\odot}]})_{3D-HST}<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 $log_{10}(\tau/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_{3D-HST}$ and M∗ detailed in this section, whereas all further analysis in this paper is based on our own estimates (Section 2.3).

2.2. Data

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 $1\farcs 2$ 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 $0\farcs 6$ .

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-22$ calibration stars on each exposure. Due to atmospheric differential refraction (ADR) comparable to the $1\farcs 2$ 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 $\sim 15\%$ . After this calibration, we correct the fluxes in our spectra for Galactic extinction. Samples of sky-subtracted spectra are shown in Figure 1.

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

2.3. Sample Properties

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 $log_{10}(\tau/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 $\sim 7$ for strong emitters (de Barros et al. 2014). However, our galaxies are at $z\sim 4$ . At this redshift, the overestimate is within a factor of $\lesssim 1.1$ (Stark et al. 2013; Salmon et al. 2015), since the 4.5 $\mu m$ Spitzer/IRAC band is unaffected111For the redshift range of our galaxies, the H $\beta$ , O II , and O III emission lines fall between the $H_{160}$ and 3.6 $\mu m$ bands. H $\alpha$ only contributes to the 3.6 $\mu m$ Spitzer/IRAC band for $z>3.75$ .. 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<z_{EAZY}<4.25$ , with a median uncertainty $\sigma_{EAZY}=0.1$ (Figure 2). EAZY outputs also include $2\sigma$ constraints, which for our sample are limited to $2.95<z<4.5$ . Our FAST executions yield a mass coverage of $7.6<\log{(\mbox{M}_{*}[\mbox{M}_{\odot}])}<10.6$ (Figure 3), with a characteristic uncertainty of $\log{(\mbox{M}_{*}[\mbox{M}_{\odot}])}\sim 0.3$ . Most of our SFR values are in the range $1-100\mbox{M}_{\odot}/yr$ (further analysis and plots below). From now on, we use the FAST outputs based on the spectroscopic redshifts ( $z_{Ly\alpha}$ ) for the 120 detections (S/N $\geqslant 5.5$ ; Section 3.1) and the outputs based on $z_{EAZY}$ for the 505 non-detections.

2.4. Considerations

The flux calibration procedure performed in our spectra is based on using stars; therefore, it corrects for a $\sim 32$ % fiber flux loss, which corresponds to a point-source Ly $\alpha$ surface brightness distribution. In cases of extended Ly $\alpha$ 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 $\alpha$ halos. Furthermore, Ly $\alpha$ 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. 2013b). At the cosmic time of our sample, the fiber diameter corresponds to a scale of $\sim 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 $\Delta z=z_{Ly\alpha}-z_{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 (see Herenz et al. 2017). We find the offset to correlate with $W_{Ly\alpha}$ (see Figure 2), hinting at biases in photometric redshift fitting when a strong Ly $\alpha$ 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 $\log{(\mbox{M}_{*}[\mbox{M}_{\odot}])}<8.5$ at $z\sim 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 $\alpha$ emission statistics dependence on galaxy properties. As a consequence, it is by no means representative of the M∗ or Ly $\alpha$ 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).

3. Ly $\alpha$ Measurements

3.1. Line Detection

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

[TABLE]

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 $\sigma$ and none above 5 $\sigma$ . Therefore, down to 5 $\sigma$ , we expect at most two false detections in our $120$ lines, translating into $\lesssim 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}>$ S/N ${}^{*}|$ S/N), with S/Ni the measured signal-to-noise ratio and S/N*∗* the imposed detection threshold. We define the simulated signal to noise as S/Nsim=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}-10^{-17}$ ergs s*-1* cm*-2*, 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*-1* cm*-2* (see Figure 1).

3.2. Line Profiles

The radiative transfer and escape of Ly $\alpha$ 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 $\alpha$ 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

[TABLE]

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

[TABLE]

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 that 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.

3.3. Ly $\alpha$ Equivalent Width and Escape Fraction

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

[TABLE]

with $F$ the Ly $\alpha$ 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 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 to match the photometry at $\sim 1700$ Å of every object to a typical agreement of $\sim 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 $\sigma$ errors in the spectra.

On the other hand, $f_{esc}$ is the fraction of the number of Ly $\alpha$ photons that escape the galaxy from the number produced. This diagnostic is typically indirectly recovered using SFRs derived from the Ly $\alpha$ line and intrinsic SFRs. The latter are typically calculated from UV continuum measurements or H $\alpha$ 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)

[TABLE]

with $L_{Ly\alpha}$ the luminosity associated with the Ly $\alpha$ line and SFR(UV)corr the extinction-corrected SFRs. This assumes a Ly $\alpha$ to H $\alpha$ ratio of 8.7 and the H $\alpha$ 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_{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).

3.4. AGN contamination

We do not measure any objects to have more than one emission line with S/N ${}_{i}>$ S/N*∗, dismissing the existence of evident active galactic nuclei (AGNs) in our sample. We do not find any of the 120 Ly $\alpha$ -emitting galaxies to have S/N $>1$ emission potentially associated with the C IV $\lambda 1549$ line. Moreover, we find no broad-line Ly $\alpha$ emission ( $\Delta v>1000$ km s-1*), 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 ( $\sim 0.4$ ; Finkelstein et al. 2009). This is consistent with the fact that bright AGNs are found in $\lesssim$ 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 (NED222The NASA/IPAC Extragalactic Database (NED) is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration.) for all objects with $f_{esc}>1$ , and we find none of them to be a reported X-ray source.

4. The Ly $\alpha$ Equivalent Width Distribution

4.1. Bayesian Inference

Measurement of $W_{Ly\alpha}$ for a galaxy sample yields the $W_{Ly\alpha}$ distribution. Since $W_{Ly\alpha}$ is directly measured from the data, characterization of the $W_{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_{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_{Ly\alpha}$ distribution measurements. For instance, studies use Gaussian (e.g. Guaita et al. 2010), exponential (e.g. Jiang et al. 2013a; Zheng et al. 2014), and log-normal distributions. For our analysis, we define the probability distribution as $p(W_{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 datasets. 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 $\alpha$ line flux $F$ instead of $W_{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_{Ly\alpha}$ . According to Bayes’ theorem, the posterior distribution $p(W_{n}|\{{\mbox{$ F $}}\})$ , i.e., the parameter space probability distribution given our data set $\{F\}$ , is

[TABLE]

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

[TABLE]

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_{Ly\alpha}$ from Equation (4), the term translates to the product distribution $p(F|W_{n})=p(W_{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}=$ S/N ${}^{*}\sigma_{i}$ . For galaxies with no detections above $F^{*}_{i}$ , we adopt the following value for the likelihood:

[TABLE]

with $p(F_{i}>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:

[TABLE]

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.

4.2. Model Comparison

As we introduced in the previous section, multiple probability distributions can be adopted for the representation of the $W_{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 dataset, i.e., it can be applied to any dataset 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_{Ly\alpha}|W_{n})$ , the modified counterparts we consider are given by

[TABLE]

where the first term is the scaled probability distribution. It is multiplied by the Heaviside $H(W_{Ly\alpha})$ to ensure it only represents positive $W_{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 $\alpha$ (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,

[TABLE]

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\}$

[TABLE]

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.,

[TABLE]

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 dataset is

[TABLE]

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

[TABLE]

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_{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_{Ly\alpha}$ tail. Still, when integrating the likelihoods, the lognormal 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_{Ly\alpha}$ measurements fairly well, despite depending on only two parameters. In addition, the uncertainties in Figure 5, especially for low $W_{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_{Ly\alpha}$ end of the distribution, which includes our completeness and non-detections. Nonetheless, further $W_{Ly\alpha}$ distribution analysis in this paper is mostly focused on the higher $W_{Ly\alpha}$ end and is not strongly dependent on model preference.

From now on, we perform our $W_{Ly\alpha}$ distribution analysis using the exponential model of Equation (11). We stress that this expression is dependent on the parameters $A$ and $W_{0}$ , 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.

5. Ly $\alpha$ Emission Dependence on Galaxy Properties

5.1. Stellar Mass

Evidence suggests Ly $\alpha$ 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 $\alpha$ photon absorption, decreasing $W_{Ly\alpha}$ . This effect has already been observed, at least for high $W_{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 $\alpha$ photon budget of galaxies. As a matter of fact, Ly $\alpha$ emission decreases steadily with the age of stellar populations, as seen in Charlot & Fall (1993) and Schaerer (2003). Ly $\alpha$ 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 $\alpha$ photons should be subject to more resonant scattering, therefore decreasing their $W_{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 $\alpha$ emission dependence on M∗.

Our sample is especially designed to study the dependence of Ly $\alpha$ emission in M∗. As shown in Figure 3, our objects are selected in redshift and M∗, homogeneously covering the range $7.6<log(\mbox{M}_{*}[\mbox{M}_{\odot}])<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 $\alpha$ emission dependence on M∗ becomes clearer in Figure 7, where we plot $W_{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_{Ly\alpha}$ for galaxies with brighter UV continuum. Therefore, even though we observe lower M∗ galaxies to have higher $W_{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 ${}_{*}<10^{8.5}$ , there is a clear anti-correlation between $W_{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 $\alpha$ 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_{Ly\alpha}$ .

The overall dependence of Ly $\alpha$ emission on M∗ implies that $W_{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 MUV surveys, shallower samples are bound to observe lower $W_{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_{Ly\alpha}$ distributions in Figure 8. As expected, there is an apparent anti-correlation between M∗ and both, the tail of the $W_{Ly\alpha}$ distribution and the normalization. We perform our first quantitative characterization of the $W_{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_{0}$ (M∗) using expressions of the form

[TABLE]

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_{\mbox{\scriptsize M}_{*}},A_{C},W_{\mbox{\scriptsize 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

[TABLE]

with $C$ a normalization constant.

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

[TABLE]

In the framework of an exponential profile, these relations we recover yield the $W_{Ly\alpha}$ probability distribution for an object with known M∗. Hence, we can simulate the expected $W_{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_{Ly\alpha}$ distributions.

Our results regarding the dependence of Ly $\alpha$ on M∗ are conclusive. In the range $10^{8}-10^{10.5}\mbox{M}_{\odot}$ , the $W_{Ly\alpha}$ probability distribution extends to higher $W_{Ly\alpha}$ for lower M∗ galaxies. In other words, more massive galaxies tend to have lower $W_{Ly\alpha}$ . A similar trend is observed for $f_{esc}$ , further highlighting the role of dust and gas mass in the escape of Ly $\alpha$ photons. At $z\sim 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\sim 2$ . These inferences, combined with the significantly higher $f_{esc}$ they measure for larger apertures, make sense in a Ly $\alpha$ 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_{Ly\alpha}$ -selected LAEs follow a trend similar to the one we find at $z\sim 4$ . In summary, the evidence for an anti-correlation between $W_{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 $\alpha$ luminosity-selected LAE number distribution to depend on M∗. Their results agree with the $z\sim 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_{\mbox{\scriptsize 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_{0}$ and M∗ we find (see Equation 22).

5.2. SFR

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_{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 $\alpha$ photons from galaxies (Verhamme et al. 2006). In fact, it has also been suggested that photoelectric absorption rules Ly $\alpha$ depletion, even over dust attenuation (Reddy et al. 2016). In this section, we explore any Ly $\alpha$ dependence on SFR within our dataset. 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_{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_{Ly\alpha}$ ( $f_{esc}$ ) and SFR is observed. These results come as no surprise, as they have been previously reported. Most studies of $W_{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\sim 2$ H $\alpha$ emitters sample, Matthee et al. (2016) also observe a clear anti-correlation between Ly $\alpha$ $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 dataset includes H $\alpha$ fluxes, they can recover SFRs and $f_{esc}$ using H $\alpha$ luminosities. The fact that they observe similar trends with such a different sample suggests that the anti-correlation between $W_{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\mbox{M}_{\odot}$ yr*-1*. Most of our low-mass objects have SFRs lower than $\sim 5\mbox{M}_{\odot}$ yr*-1*, 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 SFR $\sim 1-3\mbox{M}_{\odot}$ yr*-1*. These numbers are consistent with the analysis by Atek et al. (2014). They compare their $z<0.5$ SFR(Ly $\alpha$ )/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. 2013a). Our finding that $f_{esc}$ reaches values of 100% toward SFR $\sim 1-3\mbox{M}_{\odot}$ /yr overestimates SFR(Ly $\alpha$ )/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.

5.3. UV Luminosity

In this section, we analyze the MUV distribution of our sample, while also exploring any correlations between $W_{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 MUV galaxies than a population-representative subsample.

In order to determine MUV for our objects, we use the CANDELS $i_{775}$ band. We show in Figure 10 the corresponding distribution of the complete sample and detections. We also present the dependence of $W_{Ly\alpha}$ on MUV in this figure. As expected from the SFR- $W_{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 $\alpha$ photons. In combination, higher dust extinction, older age of stellar populations, and greater neutral gas mass in UV brighter galaxies seem to rule Ly $\alpha$ statistics dependence on SFR.

Most analyses of Ly $\alpha$ emission dependence on MUV have been performed using uncorrected SFRs. The literature compilation shown in Verhamme et al. (2008) reveals how observed UV SFRs anti-correlate with $W_{Ly\alpha}$ from z $\sim 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 MUV 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\sim 6$ . More recent studies confirm these trends (Jiang et al. 2013a, 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_{Ly\alpha}$ dependence on MUV/SFRobs can be a consequence of the role played by dust. Galaxies brighter in the UV naturally have a greater Ly $\alpha$ photon production, but the anti-correlation with M∗ affects the escape fraction (probably through increased dust extinction). In this scenario, Ly $\alpha$ photon escape is a complex process simultaneously ruled by different properties of high-redshift galaxies. We explore such property space in our analysis of the $W_{Ly\alpha}$ distribution dependence on the M ${}_{UV}-\beta$ relationship in Section 5.5.

Implications involving the fraction of detections as a function of MUV 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_{Ly\alpha}$ . Nevertheless, our Ly $\alpha$ measurements are flux limited. Therefore, our detection completeness in $W_{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_{\lambda}$ , respectively). Hence, we encourage further interpretations of these results to focus on the analysis performed on the M ${}_{UV}-\beta$ plane (Section 5.5).

5.4. UV slope

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 $\alpha$ 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 $\alpha$ 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 $\alpha$ emission at $z\sim 4$ . In this section, we detail our method to estimate the UV slopes for our objects and show our results on the Ly $\alpha$ 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, AV, 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

[TABLE]

where we assume a pristine slope of $\beta=-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\sim 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 $\alpha$ 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.

We also show in Figure 11 our results on $W_{Ly\alpha}$ and $f_{esc}$ as a function of $\beta$ 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_{Ly\alpha}$ / $f_{esc}$ , although with significant scatter. As extinction seems to play a major role in Ly $\alpha$ 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 $\beta$ (see also Blanc et al. 2011), it is consistent with a scenario where the observed Ly $\alpha$ 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 $\alpha$ 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 $\alpha$ emitters up to $E(B-V)\sim 1$ . Similarly, Matthee et al. (2016) find a population of dusty LAEs at $z\sim 2.23$ , and they speculate on how dusty gas outflows might be the feature driving the escape of Ly $\alpha$ radiation. Ly $\alpha$ sources up to $z\sim 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 $\alpha$ emission has also been measured in sub-millimeter galaxies (e.g. Chapman et al. 2005). On the other hand, our 625 object sample features $\sim 5$ objects with $E(B-V)>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\sim 4$ . True enough, the fraction of dusty galaxies at $z\sim 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 $\alpha$ emitting galaxies are mostly low-dust objects (e.g. Song et al. 2014), there is also a population of dustier, low- $W_{Ly\alpha}$ LAEs. Still, the significance of their numbers at $z\gtrsim 4$ is still an open question.

Our results also confirm at $z\sim 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\sim 2.2$ , Blanc et al. (2011) at $1.9<z<3.8$ , Song et al. (2014) at $z\sim 2.1-2.5$ , and Matthee et al. (2016) at $z\sim 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_{Ly\alpha}$ , very dusty LAEs do not seem to be common, driving the $E(B-V)-f_{esc}$ relation down significantly for $E(B-V)>0.2$ . However, at lower redshift, such objects are actually observed, driving the relation up for $E(B-V)>0.2$ .

5.5. MUV- $\beta$ sequence

We have characterized in this paper the dependence of Ly $\alpha$ emission on M∗, SFR, UV luminosity, and UV slope. We find $W_{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 MUV and $\beta$ has been studied in LBGs (Bouwens et al. 2009, 2014). As we discuss throughout this work, Ly $\alpha$ 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 $\alpha$ emission dependence on multi-dimensional spaces. In this section, we study $W_{Ly\alpha}$ dependence on MUV and $\beta$ . 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 $\alpha$ 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 $\alpha$ emission in the M ${}_{UV}-\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 location 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_{Ly\alpha}$ and MUV is still clearly observed when comparing the low MUV with the high MUV subsamples, independent of $\beta$ . For the UV slope, however, the trend between $W_{Ly\alpha}$ and $\beta$ 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 MUV and $\beta$ . We start again with linear expressions of the form

[TABLE]

Once the parameterization and priors are set, we can obtain the posterior distribution of $\{W_{n}\}=(A_{M_{UV}},A_{\beta},A_{C},W_{M_{UV}},W_{\beta},W_{C})$ using Equation (9). We assume the priors to be independent, which translates to $p(\{W_{n}\})\propto\displaystyle\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_{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:

[TABLE]

The value of the parameters $A$ and $W_{0}$ 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_{0}$ 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}<0$ , we impose $A=0$ , i.e., the $W_{Ly\alpha}$ distribution is a Dirac delta (i.e., no Ly $\alpha$ 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:

[TABLE]

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 $\beta$ with a significance $\gtrsim 2\sigma$ . There also seems to be an anti-correlation between $A$ and MUV, but to a level of $\sim 1\sigma$ . For $W_{0}$ , the e-folding scale of the distribution, the behavior is the opposite. The extent of the distribution correlates with MUV and shows a weak anti-correlation with $\beta$ .

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 $\alpha$ . On the other hand, the stronger correlation of MUV with $W_{Ly\alpha}$ suggests that the UV luminosity is associated with physical processes that determine the magnitude of the resulting $W_{Ly\alpha}$ . A possible scenario is one in which MUV traces SFR and, therefore, the total cold gas mass of galaxies. Hence, higher UV luminosities translate into increased levels of Ly $\alpha$ photon scattering into an extended halo, smoothly reducing the $W_{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 $\alpha$ 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 “covering factor” scenario could explain that $\beta$ correlates more significantly with the normalization ( $A$ ) than with the shape of the distribution ( $W_{0}$ ). This is consistent with the finding that Ly $\alpha$ 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 $\alpha$ 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 $\alpha$ radiation.

Naturally, this interpretation is a simplified description of Ly $\alpha$ escape from galaxies. Going no further, Figure 13 reveals how both distribution parameters do not solely depend on one property. Ideally, studies of $W_{Ly\alpha}$ and $f_{esc}$ in three-parameter spaces (e.g. M∗, MUV, $\beta$ ) can yield predictive and more accurate parameterization of Ly $\alpha$ emission. However, such analysis must be performed in much larger datasets than ours, at least if significant enough results are to be obtained. Furthermore, hydrodynamical simulations can also give insight into how Ly $\alpha$ 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.

6. Ly $\alpha$ Dependence on Sample Selection

6.1. LBG Samples

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\geqslant 3$ allows for efficient detection from ground telescopes. By only requiring the use of a few broadband filters, several galaxies that 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 MUV sensitivity of the survey, creating incompleteness at low-SFR. In the case of Ly $\alpha$ 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_{Ly\alpha}$ , leading to the possibility that Ly $\alpha$ studies in LBG samples are missing the highest $W_{Ly\alpha}$ galaxies in the universe. Indeed, the deepest surveys nowadays reach $\sim 50\%$ completeness at observed m ${}_{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 $\sim 50\%$ completeness at M ${}_{UV}\sim$ -18.3, -19, -19.5, and -19.7 ( $z\sim$ 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_{Ly\alpha}$ observed in narrowband surveys are systematically higher (e.g. Zheng et al. 2014).

Our dataset and results can be useful for characterizing the effects that selection techniques can have on Ly $\alpha$ 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 MUV functions of the galaxy population, correcting requires assumption of M∗ and/or MUV distributions. Nevertheless, we can still simulate Ly $\alpha$ emission samples on CANDELS galaxies using our results of $W_{Ly\alpha}$ dependence on the M ${}_{UV}-\beta$ plane given by Equations (27) and (28). In this section, we describe our simulation of a $z\sim 4$ LBG sample from the 3D-HST catalogs. We then compare the properties of LBGs with the parent distribution, focusing on M∗, SFR, MUV, and $\beta$ . We finally simulate Ly $\alpha$ fractions for each sample and conclude on the effects of LBG selection at $z\sim 4$ .

LBGs at $z\sim 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_{160}$ images. For the $z\sim 4$ LBG selection, we adopt the same methodology applied in Bouwens et al. (2012):

[TABLE]

After color selection, we impose $5\sigma$ detections in the $V_{606}$ 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\sigma$ 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 $\sim 3200$ CANDELS galaxies that comply with our selection, $\sim 800$ classify as low-redshift interlopers ( $z_{EAZY}<3$ ), according to our outputs. We check EAZY 3 $\sigma$ constraints for these interlopers and find most to have $z_{99}<4$ . Hence, we find a low-redshift contamination in $z\sim 4$ CANDELS LBGs that is higher than the typically reported number of $\sim 10\%$ (Bouwens et al. 2007). From now on, we remove these presumed $z_{EAZY}<3$ contaminants and work with a $z\sim 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 ${}_{*}>10^{8.5}\mbox{M}_{\odot}$ . We construct two $3<z<4.2$ CANDELS samples for comparison: all galaxies with M ${}_{*}>10^{8.5}\mbox{M}_{\odot}$ and all $V$ -band detected M ${}_{*}>10^{8.5}\mbox{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 to color requirements are minor. It is still worth noting, however, that these color criteria seem to neglect blue rather than red galaxies at $z\sim 4$ (center of Figure 14). We find the $(B-V>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, b, 2016). In analyses that explore Ly $\alpha$ emission strength, UV continuum properties, morphologies, and sizes, they find no significant differences between LBGs and LAEs at $z\geqslant 6$ .

The insights we present can be further tested by simulating the fraction of LAEs above an $W_{Ly\alpha}$ threshold. We use our results from Section 5.5 to estimate the Ly $\alpha$ 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_{Ly\alpha}$ . To deal with this selection, it has been proposed to compare the Ly $\alpha$ 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 $\alpha$ fraction with cosmic time. Regarding color selections, we find the Ly $\alpha$ 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.

6.2. Narrowband Samples

Ly $\alpha$ -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_{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_{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_{Ly\alpha}$ cuts. For instance, Vargas et al. (2014) select sources with $W_{Ly\alpha}>20$ Å at $z\sim 2.1$ , Zheng et al. (2014) select $W_{Ly\alpha}>9$ Å at $z\sim 4.5$ , and Sobral et al. (2017) go as low as $W_{Ly\alpha}>5$ Å at $z\sim 2.23$ . If $W_{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 ${}_{UV}-\beta$ model. We show in Figure 15 the outcome of $W_{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 $>10^{8.5}\mbox{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.

To explore the effect of $W_{Ly\alpha}$ selections (top panel of Figure 15), we base our analysis on the $z\sim 2.23$ narrowband survey of Sobral et al. (2017), whose galaxies are limited to $W_{Ly\alpha}>5$ Å and fluxes $>2\times 10^{-17}$ erg cm*-2* s*-1* (5 $\sigma$ ). We simulate a flux-limited-only Ly $\alpha$ sample and use it as baseline (black). We then show the region where $W_{Ly\alpha}>5$ Å (1 $\sigma$ ; blue) and $W_{Ly\alpha}>20$ Å (1 $\sigma$ ; green) objects lie. The $W_{Ly\alpha}>5$ Å contours are intended to reproduce Sobral et al. (2017) selections, whereas the more restrictive selection $W_{Ly\alpha}>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_{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_{Ly\alpha}\leqslant 20$ Å. Still, selections based solely on $W_{Ly\alpha}$ systematically fail to remove AGNs and neglect bright Ly $\alpha$ emitters (Sobral et al. 2017). All these insights remark on the importance of using low $W_{Ly\alpha}$ cuts and thorough interloper controls.

After narrowband outcomes are used for sample selection, spectroscopic observations follow. However, these follow-ups of Ly $\alpha$ 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_{Ly\alpha}>L_{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\sim 2-3$ LAEs alongside the $z\sim 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 $>7-10\times 10^{-17}$ erg cm*-2* s*-1* (5 $\sigma$ ) and $W_{Ly\alpha}>20$ Å. In contrast, the Vargas et al. (2014) flux depth is higher than that of HETDEX, reaching $2\times 10^{-17}$ erg cm*-2* s*-1* (5 $\sigma$ ), while also selecting sources with $W_{Ly\alpha}>20$ Å. Therefore, Hagen et al. (2014) are comparing two $W_{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*⊙* 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.2Z_{\odot}$ .

Vargas et al. (2014) and Hagen et al. (2014) find M ${}_{*}<10^{9}$ M*⊙* LAEs to lie above the M∗-SFR relation. A similar trend is reported in the work of Karman et al. (2017), who study Ly $\alpha$ emitters down to $10^{6}$ 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 ${}_{*}<10^{9}$ M⊙* Ly $\alpha$ 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 $\alpha$ emitters at $z\sim 4.5$ . To answer whether low M∗ LAEs lie above starforming galaxies or the slope of the sequence changes toward M ${}_{*}<10^{9}$ M*⊙, different approaches are required. For instance, no discrepancies are found between $z\sim 2$ LAEs and H $\alpha$ sources (see Matthee et al. 2016), although most of their galaxies have M ${}_{*}>10^{9}$ M⊙. In summary, evidence suggests that $z>2$ , M ${}_{*}<10^{9}$ M⊙* 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.

7. Inferences on the $4\leqslant z\leqslant 7$ Ly $\alpha$ Fraction

Probably the most important use for Ly $\alpha$ emission is tracing the neutral hydrogen fraction in the IGM. Several studies over the last decade have constrained the fraction of Ly $\alpha$ -emitting galaxies as a function of redshift (Stark et al. 2010, 2011; Ono et al. 2012; Schenker et al. 2012; Tilvi et al. 2014; Cassata et al. 2015), bringing us closer to the goal of constraining the epoch of reionization. The Ly $\alpha$ fraction is fairly well understood at $z<5$ (Stark et al. 2011; Cassata et al. 2015), with most efforts nowadays focusing on $z\sim 7,8$ (Ono et al. 2012; Treu et al. 2013; Tilvi et al. 2014; Furusawa et al. 2016). Along these lines, our characterization of $W_{Ly\alpha}$ on the M ${}_{UV}-\beta$ plane can be used to simulate $W_{Ly\alpha}$ distributions at higher redshift. In this section, we apply the observed M ${}_{UV}-\beta$ relations from Bouwens et al. (2014) and LFs from Bouwens et al. (2015) to simulate high-redshift $W_{Ly\alpha}$ distributions. By means of these simulations, we can predict the Ly $\alpha$ fraction in galaxies up to $z\sim 7$ , providing the first semi-analytical constraint to this tracer toward the reionization epoch. Furthermore, we also simulate dropouts from CANDELS photometry to explore the effects of observational limitations on the inferred Ly $\alpha$ fractions at high redshift.

Our results from this section are based on assuming the same $z\sim 4$ $W_{Ly\alpha}$ dependence on M ${}_{UV}-\beta$ at every redshift. We remark that testing this assumption with currently available datasets is challenging, since different sample selections and incompleteness levels can affect the observed $W_{Ly\alpha}$ dependence on the M ${}_{UV}-\beta$ plane. It is also worth noting that the analysis described in this section does not account for any effects related to changes in the merger fraction, IGM opacity (e.g. Gunn & Peterson 1965; Becker et al. 2001; Fan et al. 2006), and/or conditions of the ISM with cosmic time (e.g. Carilli & Walter 2013). As the number of complete, unbiased Ly $\alpha$ surveys grows (e.g. Cassata et al. 2015; Hathi et al. 2016; this work), we will be able to test these assumptions and explore their dependence on redshift.

We now detail our simulation of high-redshift LBG samples. For the LFs, we use Bouwens et al.’s 2015 best Schechter parameters (excluding CANDELS-EGS). We start by drawing objects following the LF and then associate a UV slope according to the best-fit relations from Table 3 of Bouwens et al. (2014). We finish by adding an intrinsic scatter of 0.35 to the slopes at every redshift (Bouwens et al. 2012; Castellano et al. 2012; Bouwens et al. 2014). Similarly, we perform an analogue procedure following the MUV distribution of CANDELS LBGs, which allows us to assess the effect of magnitude incompleteness. In order to do so, we make use of the optical and IR photometry publicly available from the 3D-HST catalogs (Skelton et al. 2014). For every redshift, we use the selections from Bouwens et al. (2015), since they are based on CANDELS photometry:

[TABLE]

Similarly to our Section 6.1 selection, we impose S/N cuts in the detection bands. For $z\sim 4$ and $5$ dropouts, we require $3\sigma$ detections in the $V_{606}$ , $z_{850}$ bands. For $z\sim 6$ and 7 dropouts, we impose at least $3\sigma$ detections in the $Y_{102}$ and $J_{125}$ bands, respectively. To associate a dropout redshift with every galaxy, however, we use the photometric redshifts from 3D-HST instead of the actual dropout from where it was selected. By doing so, we are in agreement with the methodology of Bouwens et al. (2015). We finally complete these CANDELS LBG samples by associating UV slopes to every galaxy. Just as for the complete sample, we associate the slopes following the M ${}_{UV}-\beta$ relations from Bouwens et al. (2014).

We present in Figure 16 the assumed M ${}_{UV}-\beta$ sequences and the UV luminosity distribution for both sets under consideration. Since we are starting from the Bouwens et al. (2014) M ${}_{UV}-\beta$ characterization, the steepness of the linear relation increases with redshift. The contours we show in Figure 16 (left panel) are obtained after incorporating the intrinsic scatter, which we assume to be of 0.35. Note that the addition of such a scatter does not render the redshift evolution of the M ${}_{UV}-\beta$ distribution negligible. According to our Equations (27) and (28), this implies that the probability of observing $W_{Ly\alpha}$ in emission at fixed MUV should increase with redshift. The right panel of Figure 16 shows the absolute magnitude distribution of $4\leqslant z\leqslant 7$ LBGs as found by Bouwens et al. (2015, solid lines). Note the growth of the LF as a function of redshift. Since the LF becomes steeper as the redshift increases, the probability of observing a high- $W_{Ly\alpha}$ galaxy should also increase. One more important point regarding Figure 16 must be noted. As initially suggested, CANDELS incompleteness affects the observed MUV distribution of LBGs (shown as crosses). The fact that there is more contribution from brighter galaxies could, in principle, bias samples toward lower $W_{Ly\alpha}$ when compared to the LFs. We now focus on quantifying all these effects within the framework of our model.

We simulate the $W_{Ly\alpha}$ distributions for the two sample sets we constructed. To do so, we use the MUV and $\beta$ of each object to draw $W_{Ly\alpha}$ distributions following our sampling of Equations (27) and (28). This methodology allows us to simulate the $W_{Ly\alpha}$ distribution for the number of galaxies in every redshift subsample, properly taking into account the uncertainties on each M ${}_{UV}-\beta$ domain and the errors in our modeling. We present the corresponding $4\leqslant z\leqslant 7$ Ly $\alpha$ fractions for such galaxies in Figure 17. To this end, we select using the thresholds $W_{Ly\alpha}^{*}$ = 25 Å and $W_{Ly\alpha}^{*}$ =55 Å (Stark et al. 2011). Still, as we note in Section 5.3, LBGs with MUV distributions dominated by fainter objects can yield higher fractions of Ly $\alpha$ . In recent works, characterization has been performed separately in galaxies with M ${}_{UV}<-20.25$ and M ${}_{UV}>-20.25$ (Stark et al. 2011; Ono et al. 2012; Schenker et al. 2012; Tilvi et al. 2014; Cassata et al. 2015). Following the same approach, we now focus our analysis on samples in the ranges $-21.75<$ M ${}_{UV}<-20.25$ and $-20.25<$ M ${}_{UV}<-18.75$ (Figure 16), which are well covered by our dataset (see Figure 12).

We plot in Figure 17 the resulting Ly $\alpha$ fractions. We find the complete fractions (i.e., derived from the LFs; black circles) to be consistent with a steady increase up to $z\sim 7$ , which has been previously assumed as baseline for measuring drops in the Ly $\alpha$ fraction (Stark et al. 2010, 2011; Ono et al. 2012). Apart from backing this assumption, our result provides a more robust baseline for comparison. We reiterate that a closer look at Figure 16 reveals the explanation for such trend, at least within our approach. First, the LF of LBGs becomes steeper as redshift increases (Bouwens et al. 2015). Second, the M ${}_{UV}-\beta$ relation favors bluer slopes as redshift increases. The combination of both relationships explain an increase in the Ly $\alpha$ fraction as a function of redshift for LBGs. It has already been discussed (Blanc et al. 2011; Cassata et al. 2015) how the growth in the Ly $\alpha$ fraction is likely tied to changes in the $f_{esc}$ of Ly $\alpha$ photons, at least for $3<z<5$ . This change in the $f_{esc}$ is likely tied to evolution in the dust content and ISM of galaxies, which we are empirically accounting for in this model. Under this picture, Lyman continuum (LyC) photon escape should follow a similar trend. Indeed, Faisst (2016) empirically predicts that LyC photon escape fraction increases with redshift at fixed M∗. In summary, our results are in strong agreement with the picture that UV photon escape fraction from galaxies tends to decrease with cosmic time.

From the literature, we include in Figure 17 the measurements from Fontana et al. (2010), Pentericci et al. (2011), Stark et al. (2011), Vanzella et al. (2011), Ono et al. (2012), Schenker et al. (2012), and Treu et al. (2013). The overall literature trend is an increase in the Ly $\alpha$ fraction up to $z\sim 5$ , which flattens at $z\sim 6$ and drops toward $z\sim 7$ . Departures from this increasing regime are typically associated with changes in the transparency of the IGM, with drops toward $z\sim 7$ associated with a more neutral universe (Stark et al. 2010; Ono et al. 2012). Even though such inferences might be correct, the contrast between the completeness-corrected (black circles) and -uncorrected (black crosses) fractions in Figure 17 suggests that the driving factor behind these $z\sim 7$ drops might be sample MUV incompleteness. Our simulated $z\geqslant 6$ CANDELS LBG samples do not follow the corresponding LF as a result of survey depth (see Figure 16), leading them to include more luminous objects than a population-representative sample. Therefore, as seen in Figure 17, drops in the Ly $\alpha$ fraction are induced. As we point out in Section 6.1, the deepest LBG surveys nowadays (e.g Bouwens et al. 2015; Bowler et al. 2015, 2017; Ishigaki et al. 2017) reach $\sim 50\%$ completeness at M ${}_{UV}\sim$ -18.3, -19, -19.5, and -19.7 for $z\sim$ 4, 5, 6, and 7, respectively. Therefore, Ly $\alpha$ studies of even the deepest LBG surveys are sufficiently incomplete to bias the reported Ly $\alpha$ fractions (see Figures 16 and 17). Taking into account this effect is essential if we are to use Ly $\alpha$ to further constrain the reionization epoch during the next decade.

It is important to clarify the limitations of our modeling. First, our characterization is based on the whole $W_{Ly\alpha}$ distribution, which may misrepresent the actual fraction of high- $W_{Ly\alpha}$ galaxies. In fact, the analysis we perform in Section 4 reveals that a log-normal profile might be better for reproducing the tail of the distribution (Figure 5). For example, Zheng et al. (2014) find an exponential profile to only adequately represent the low- $W_{Ly\alpha}$ end of the distribution, leading them to model the high- $W_{Ly\alpha}$ tail separately. As a consequence, we encourage the reader to focus most of their analysis on our $W_{Ly\alpha}^{*}$ = 25 Å results (top panel in Figure 17). Along those lines, we expect our modeling to underestimate the fraction for populations with significantly higher- $W_{Ly\alpha}$ tails. Indeed, our fractions for galaxies with $-20.25<$ M ${}_{UV}<-18.75$ are systematically lower than Stark et al. (2011) measurements (see Figure 17).

In summary, our results suggest that the Ly $\alpha$ fraction steadily increases up to $z\sim 7$ in an ionized universe. However, drops in this fraction are not necessarily related to changes in the opacity of the IGM, since we find that sample MUV incompleteness can alternatively explain them. In fact, the uncertainties in both literature measurements and our predictions are still large enough to set assertive constraints on the Ly $\alpha$ trends toward $z\sim 7$ (see Robertson et al. 2015). Our results, therefore, further highlight that Ly $\alpha$ constraints on $z_{re}$ (the redshift at which the fraction of ionized hydrogen is 0.5) are still inconclusive.

It is worth discussing whether $z_{re}>7$ is in agreement with observational reionization constraints besides Ly $\alpha$ . Recently, the Planck Collaboration et al. (2016b) analysis of CMB anisotropies has yielded $z_{re}\sim 8.2^{+1.0}_{-1.2}$ . Systematic detection of Ly $\alpha$ -emitting galaxies in the range $z\sim 7-9$ (Zitrin et al. 2015; Furusawa et al. 2016; Stark et al. 2017) still leads to the question whether reionization is already in place at $z\sim 7$ . Several studies have also explored Ly $\alpha$ LFs toward $z\sim 9$ , although observational limitations require deeper/wider surveys to obtain more meaningful conclusions (Sobral et al. 2009). Therefore, and since reionization is predicted to happen so rapidly (Finlator et al. 2009), claims of $z_{re}>7$ are not in tension with measurements of the Ly $\alpha$ fraction at $z\sim 7$ . Interestingly enough, it is not even clear whether reionization happens first around fainter, lower density environments (“outside-in,” Kashikawa et al. 2006), or brighter, higher density regions (“inside-out,” Santos et al. 2016). This overall picture suggests that the role of MUV in reionization studies will no longer restrict itself to sample control, but it will also extend to parameterizing property. Therefore, the natural approach for future studies would be to constrain reionization as a smooth function of both redshift and UV luminosity/environment (e.g. Robertson et al. 2013).

8. Conclusions

In this work, we present an exhaustive analysis of Ly $\alpha$ emission at $3<z<4.6$ . To this end, we M∗-select 625 galaxies from the CANDELS survey, allowing us to study Ly $\alpha$ emission over a diverse and heterogeneous galaxy sample. We conduct spectroscopic observations of our targets with the M2FS, a multi-fiber spectrograph at the Clay 6.5m telescope. We then use a Bayesian approach for proper statistical handling of our results. By means of this framework, we are capable of characterizing Ly $\alpha$ emission in the high-redshift galaxy population. In summary, our conclusions are the following.

We present a Bayesian methodology to measure the $W_{Ly\alpha}$ distribution considering the completeness in fluxes measured and the uncertainties in both spectroscopy and photometry. We combine this approach with Monte Carlo simulations for robust modeling of observed trends. Combining all of these features allows us to properly state the significance and uncertainties in our results.
We take advantage of the Bayesian framework to give an insight into how to compare multiple $W_{Ly\alpha}$ distribution models from a quantitative standpoint. In this paper, we explore the extent to which exponential, Gaussian, and log-normal probability distributions can reproduce $W_{Ly\alpha}$ measurements. We conclude that, in our case, an exponential profile is the most adequate. However, as we find in Section 5, this profile can struggle to reproduce the tail of high- $W_{Ly\alpha}$ populations (e.g. samples with low M∗, SFR, UV luminosity, or $\beta$ ).
Our measured $W_{Ly\alpha}$ and $f_{esc}$ strongly anti-correlate with M∗. We associate these trends with the higher dust fraction and gas mass in more massive galaxies, boosting the scattering and absorption of Ly $\alpha$ photons. We model both exponential $W_{Ly\alpha}$ distribution parameters ( $A,W_{0}$ ) using linear relations dependent on M∗, and find them both to anti-correlate with M∗. Our modeling is also capable of reproducing our observed $W_{Ly\alpha}$ distributions when binning the sample (Figure 8).
We also explore the dependence of Ly $\alpha$ emission on MUV. We find $W_{Ly\alpha}$ and $f_{esc}$ to be typically higher for UV faint objects, which have been previously observed in the literature (Stark et al. 2010; Schaerer et al. 2011). Since $z\sim 4$ galaxies seem to follow an M ${}_{*}-$ MUV sequence (González et al. 2014), this result is consistent with our M∗ trends. This confirms that the role played by higher dust fraction and gas mass rules $W_{Ly\alpha}$ dependence on UV luminosity. We also observe $W_{Ly\alpha}$ and $f_{esc}$ to anti-correlate with the SFR and UV slope, confirming the qualitative Ly $\alpha$ escape scheme presented.
The relatively uniform location of our targets in the M ${}_{UV}-\beta$ plane allows us to characterize the $W_{Ly\alpha}$ distribution in this space. Once more, we use a Bayesian framework and Monte Carlo simulations to obtain linear representations of exponential $W_{Ly\alpha}$ distribution parameters. Our results suggest that the probability of observing Ly $\alpha$ emission from a galaxy is determined by its dust content, whereas the actual magnitude of $W_{Ly\alpha}$ is mostly ruled by the UV luminosity of the object. This characterization also allows us to simulate $W_{Ly\alpha}$ distributions from M ${}_{UV}-\beta$ galaxy samples (see below).
Using CANDELS, we mimic a $z\sim 4$ LBG survey based on B-dropouts. We simulate $W_{Ly\alpha}$ distributions from the LBG sample and verify that the lower probability of measuring high $W_{Ly\alpha}$ when compared to the population is mostly a consequence of detection band limitations. In other words, shallower LBG surveys are bound to include fewer low-M∗ galaxies in their samples, decreasing the probability of including high- $W_{Ly\alpha}$ objects. We find color constraints imposed by the LBG technique to also lower the Ly $\alpha$ fraction, although the effect is minor. These results highlight the importance of comparing surveys with similar sensitivities and MUV distributions.
Our results on the M ${}_{UV}-\beta$ plane also allow us to explore the effects $W_{Ly\alpha}$ an line-flux limitations induce on Ly $\alpha$ narrowband surveys. We find $W_{Ly\alpha}$ cuts to bias samples toward low-M∗ objects, whereas flux limitations seem to preferentially leave out low-SFR galaxies. These insights contribute to explain the location of LAEs within the main sequence (e.g. some of the differences between Vargas et al. 2014 and Hagen et al. 2014).
We generate $4\leqslant z\leqslant 7$ LBG samples following reported M ${}_{UV}-\beta$ relations and LFs from the literature (Bouwens et al. 2014, 2015). Assuming the same $W_{Ly\alpha}$ dependence on M ${}_{UV}-\beta$ that we find at $z\sim 4$ , we estimate the fraction of Ly $\alpha$ emitting galaxies in LBG samples above an $W_{Ly\alpha}$ threshold. Our findings are consistent with observational measurements in the literature that suggest an increase in the fraction up to $z\sim 6$ (Stark et al. 2010, 2011; Cassata et al. 2015). This result constitutes the first semi-analytical constraint to the Ly $\alpha$ fraction at $z\sim 7$ , replacing extrapolations of lower redshift regimes.
We simulate the $4\leqslant z\leqslant 7$ Ly $\alpha$ fraction in LBG samples drawing from CANDELS MUV distribution. We conclude that MUV incompleteness can lower the Ly $\alpha$ fraction, reproducing drops at $z\sim 7$ . Claims of drops at this redshift need to thoroughly account for this effect if they are to be used to constrain the reionization epoch. This result highlights that reported drops in the Ly $\alpha$ fraction are not inconsistent with $z_{re}>7$ . Given the uncertainties in both predictions and measurements, more work is needed to conclude on the actual Ly $\alpha$ trends toward $z\sim 7$ and beyond.

We thank the referee for useful feedback that improved the quality of this paper. G.O. was supported by CONICYT, Beca Magíster Nacional 2014, Folio 22140924. G.B. was supported by CONICYT/FONDECYT, Programa de Iniciación, Folio 11150220. V.G. was supported by CONICYT/FONDECYT, Programa de Iniciación, Folio 11160832. We also acknowledge the NSF/MRI grant for the M2FS, number AST 0923160. This work is based on observations taken by the 3D-HST Treasury Program (GO 12177 and 12328) with the NASA/ESA HST, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555. This research has made use of the NASA/IPAC Extragalactic Database (NED), which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. This paper includes data gathered with the 6.5 meter Magellan Telescopes located at Las Campanas Observatory, Chile.

Bibliography169

The reference list from the paper itself. Each links out to its DOI / PubMed record.

1Adams et al. (2011) Adams, J. J., Blanc, G. A., Hill, G. J., et al. 2011, Ap JS, 192, 5 · doi ↗
2Atek et al. (2014) Atek, H., Kunth, D., Schaerer, D., et al. 2014, A&A, 561, A 89 · doi ↗
3Atek et al. (2011) Atek, H., Siana, B., Scarlata, C., et al. 2011, Ap J, 743, 121 · doi ↗
4Balestra et al. (2010) Balestra, I., Mainieri, V., Popesso, P., et al. 2010, A&A, 512, A 12 · doi ↗
5Becker et al. (2001) Becker, R. H., Fan, X., White, R. L., et al. 2001, AJ, 122, 2850 · doi ↗
6Bennett et al. (2013) Bennett, C. L., Larson, D., Weiland, J. L., et al. 2013, Ap JS, 208, 20 · doi ↗
7Blanc et al. (2011) Blanc, G. A., Adams, J. J., Gebhardt, K., et al. 2011, Ap J, 736, 31 · doi ↗
8Bond et al. (2012) Bond, N. A., Gawiser, E., Guaita, L., et al. 2012, Ap J, 753, 95 · doi ↗

TL;DR

Contribution

Findings

Abstract

Peer Reviews

Videos

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

Abstract

Subject headings:

1. Introduction

2. Data Set

2.1. Sample Selection

2.2. Data

2.3. Sample Properties

2.4. Considerations

3. Lyα\alphaα Measurements

3.1. Line Detection

3.2. Line Profiles

3.3. Lyα\alphaα Equivalent Width and Escape Fraction

3.4. AGN contamination

4. The Lyα\alphaα Equivalent Width Distribution

4.1. Bayesian Inference

4.2. Model Comparison

5. Lyα\alphaα Emission Dependence on Galaxy Properties

5.1. Stellar Mass

5.2. SFR

5.3. UV Luminosity

5.4. UV slope

5.5. MUV-β\betaβ sequence

6. Lyα\alphaα Dependence on Sample Selection

6.1. LBG Samples

6.2. Narrowband Samples

7. Inferences on the 4⩽z⩽74\leqslant z\leqslant 74⩽z⩽7 Lyα\alphaα Fraction

8. Conclusions

A Comprehensive Study of Ly $\alpha$ Emission in the High-redshift Galaxy Population

3. Ly $\alpha$ Measurements

3.3. Ly $\alpha$ Equivalent Width and Escape Fraction

4. The Ly $\alpha$ Equivalent Width Distribution

5. Ly $\alpha$ Emission Dependence on Galaxy Properties

5.5. MUV- $\beta$ sequence

6. Ly $\alpha$ Dependence on Sample Selection

7. Inferences on the $4\leqslant z\leqslant 7$ Ly $\alpha$ Fraction