The Evolution of the Lyman-alpha Luminosity Function during Reionization

We model the evolution of the Lyα LF as a function of redshift and x _Hɪ by convolving models for the Lyα emission from LBGs and the UV LF during reionization. This enables us to disentangle the effects of redshift evolution from the evolution due to IGM absorption (Mason et al. 2015b; Mesinger 2016; Mason et al. 2018).

The LF of galaxies shows the number density of galaxies in a certain luminosity interval and is typically described using the Schechter (1976) function:

$$\begin{eqnarray}&&\phi (L){dL}=\phi * {\left(\displaystyle \frac{L}{L* }\right)}^{\alpha }\exp \left(-\displaystyle \frac{L}{L* }\right)\,d\left(\displaystyle \frac{L}{L* }\right),\end{eqnarray} \tag{ 1 }$$

where ϕ* is the normalization constant, α is the power law for faint-end slope for L<L^*, and there is an exponential cutoff at L<L^*. These parameters are known to be conditional on the observed wavelength and cosmic time, as well as the type of galaxy (e.g., Dahlen et al. 2005). The rest-frame UV LF has been measured in detail out to z ∼ 10 and is one of our best tools for studying the evolution of galaxy populations (e.g., Bouwens et al. 2016, 2015; Finkelstein et al. 2015; Oesch et al.2015).

Following Dijkstra & Wyithe (2012) and Gronke et al. (2015) we can predict the number density of LAEs by using the UV LF and Lyα luminosity probability distribution for LBGs to model the Lyα LF:

$$\begin{eqnarray}&&\begin{array}{l}{\phi }_{\mathrm{LAE}}({L}_{\alpha },{x}_{\mathrm{H\unicode{x0026A}}},z){{dL}}_{\alpha }={{dL}}_{\alpha }F\\ \quad \times {\int }_{{M}_{\mathrm{UV},\min }}^{{M}_{\mathrm{UV},\max }}{{dM}}_{\mathrm{UV}}\,\phi ({M}_{\mathrm{UV}},z)P({L}_{\alpha }\,| \,{M}_{\mathrm{UV}},{x}_{\mathrm{H\unicode{x0026A}}},z).\end{array}\end{eqnarray} \tag{ 2 }$$

Here, ϕ(M _UV , z)dM _UV is the UV LF in the range M _UV ± dM _UV /2 and is described in Section 2.2. P(L_α ∣ M _UV , x _Hɪ ), is the conditional probability of galaxies that have a Lyα luminosity L_α in L_α ± dL_α /2, given a M _UV value and neutral fraction and is described in Section 2.1.1. We integrate our Lyα LF over the M _UV range − 24 < M _UV < − 12 covering the observed range of the UV LF, − 23 < M _UV < − 17.

The factor F in Equation (16) is a normalization constant to fit the LF model to observations and can be thought of as the ratio of predicted LAEs versus the total number of LAEs recorded. If the Lyα luminosity distribution, P(L_α ∣ M _UV , x _Hɪ , z), accurately describes the luminosities of the same LBGs measured in the UV LF this factor should be F = 1 (See Section 4.3 for further discussion).

2.1.1. Lyα Luminosity Probability Distribution

The probability distribution for Lyα luminosity is derived from the Lyα rest-frame equivalent width (EW) probability density function P(EW ∣ M _UV ), where $P({L}_{\alpha }\,| \,{M}_{{\rm\small{UV}}},{x}_{{\rm\small{H\unicode{x0026A}}}})\,=P({EW}\,| \,{M}_{{\rm\small{UV}}},{x}_{{\rm\small{H\unicode{x0026A}}}})\tfrac{\partial {EW}}{\partial {L}_{\alpha }}$. We use the rest-frame EW probability distribution models by Mason et al. (2018) (based on observations by De Barros et al. 2017) who forward-model observed EW after transmission through 1.6 Gpc³ inhomogeneous reionization simulations (Mesinger 2016) at fixed average neutral fraction x _{H ɪ} = 0.01 − 0.95, with a spacing of Δx _Hɪ ∼ 0.01–0.03. We use the following relationship between Lyα luminosity in erg s⁻¹ and EW to obtain ∂EW/∂L_α :

$$\begin{eqnarray}&&{L}_{\alpha }=\mathrm{EW}\times {L}_{\mathrm{UV},\nu }\times \frac{c}{{\lambda }_{\mathrm{Ly}\alpha }^{2}}{\left(\frac{{\lambda }_{\mathrm{Ly}\alpha }}{{\lambda }_{\mathrm{UV}}}\right)}^{\beta +2}\end{eqnarray} \tag{ 3 }$$

Here, λ_α ∼ 1216 Å is the wavelength of the Lyα resonance, the rest-frame wavelength of the UV continuum is typically measured at λ _UV ∼ 1500 Å, β is the UV continuum slope, where we assume β = −2 typical for high redshifts (e.g., Bouwens et al. 2014), and c is the speed of light. L _UV,ν is the UV LD:

$$\begin{eqnarray}&&{L}_{{\rm\small{UV}},\nu }=4\pi {\left(10\right)}^{2}\times {10}^{-0.4({M}_{{\rm\small{UV}}}+48.6)}\,\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{Hz}}^{-1}.\end{eqnarray} \tag{ 4 }$$

We normalize P(L_α ∣ M _UV , x _{H ɪ} ) over the luminosity range L_α = 0 − 10^44.5 erg s⁻¹ to encompass a large Lyα luminosity interval, and within our defined M _UV range between − 24 <M _UV < − 12 (see Section 2.1).

The Lyα luminosity probability distribution, P(L_α ∣ M _UV , x _Hɪ ), is shown in Figure 1 for a few x _{H ɪ} and M _UV values. We note that as the Mason et al. (2018) models assume no evolution of the intrinsic Lyα EW distribution with redshift—the only redshift evolution is due to the increasing neutral fraction, our Lyα luminosity distribution models also depend only on x _{H ɪ} with no additional redshift evolution. Recent works by Hashimoto et al. (2017), Jung et al. (2018), and Shibuya et al. (2018) confirm a suitable approach to the Lyα EW distribution models we incorporate into our work. Each paper ultimately notes no significant evolution in the Lyα EW distribution with respect to redshift for z ∼ 5–7.

Zoom In Zoom Out Reset image size

Figure 1 demonstrates large differences in the probability distribution for UV-bright and UV-faint galaxies. For UV-faint galaxies, we expect a higher probability of Lyα luminosity emission at all x _Hɪ values but with a Lyα luminosity cutoff at L_α ≲ 10⁴² erg s⁻¹. For UV-bright galaxies, at all x _Hɪ we expect galaxies to have higher Lyα luminosity values up to L_α ≲10⁴⁴ erg s⁻¹, but this is rarer. For both UV-bright and UV-faint galaxies, as the neutral fraction increases, the probability of galaxies emitting strong Lyα overall decreases.

The M _UV dependence of P(L_α ∣ M _UV ) is a direct consequence of the EW distribution model by Mason et al. (2018). This model can be described as an exponential distribution plus a delta function:

$$\begin{eqnarray}&&\begin{array}{l}P(\mathrm{EW}| {M}_{\mathrm{UV}})=\frac{A({M}_{\mathrm{UV}})}{{\mathrm{EW}}_{c}({M}_{\mathrm{UV}})}{e}^{-\tfrac{\mathrm{EW}}{{\mathrm{EW}}_{c}({M}_{\mathrm{UV}})}}H(\mathrm{EW})\\ \quad +[1-A({M}_{\mathrm{UV}})]\delta (\mathrm{EW}),\end{array}\end{eqnarray} \tag{ 5 }$$

where

$$\begin{eqnarray}&&A({M}_{{\rm\small{UV}}})=0.65+0.1\,\tanh [3({M}_{{\rm\small{UV}}}+20.75)],\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{\mathrm{EW}}_{c}({M}_{\mathrm{UV}})=31+12\,\tanh [4({M}_{\mathrm{UV}}+20.25)]\,\mathring{\rm A} \end{eqnarray} \tag{ 7 }$$

account for the fraction of emitters and the anticorrelation of EW with M _UV , respectively. H(EW) is the Heaviside step function and δ(EW) is a Dirac delta function (see Section 2.1.3 of Mason et al. 2018 for further details). The intrinsic, emitted, distribution (i.e., x _Hɪ = 0) is an empirical model fit to observations by De Barros et al. (2017) at z ∼ 6, where it was found that UV-bright galaxies had a lower probability of being emitters, and had lower average EWs (consistent with previous findings; e.g., Ando et al. 2006; Stark et al. 2010). The model EW distribution is then painted onto galaxies in inhomogeneous reionization simulations, with different average neutral fractions, and the observed EW distribution in each of those simulations is recovered by sampling the transmission along thousands of sightlines.

Although bright galaxies have low EW compared to faint galaxies, they are, on average, less affected by neutral gas in the IGM: UV-bright galaxies are typically more massive and reside in dense regions of the universe, in large IGM bubbles that have already reionized. In the simulations we use, reionization occurs first in over-dense regions due to the excursion set formalism (Mesinger & Furlanetto 2007; Mesinger 2016). Observationally, clustering analyses show that UV-bright galaxies typically live in massive halos in dense regions (e.g., Figure 15 of Harikane et al. 2018). So, Lyα photons can escape more easily and EWs decrease at a slower rate (transmissions are already high). UV-faint galaxies have a higher intrinsic EW, on average, that decreases more rapidly than for UV-bright galaxies as the universe becomes more neutral. This is because they can be more typically found in neutral patches of IGM.

As described above in Section 2.1, we generate our Lyα LF by integrating the UV LF over the M _UV range − 24 < M _UV < − 12, covering the observed range of the UV LF, − 23 ≲ M _UV ≲ − 17. As the Mason et al. (2018) EW models were defined for − 23 ≤ M _UV ≤ − 17 to include galaxies outside of this range we set galaxies brighter than M _UV = − 23 to have the same P(L_α ∣ M _UV = − 23), and all galaxies fainter than M _UV = − 17 have the same P(L_α ∣ M _UV = − 17) (for a given x _Hɪ ).

In this paper, we use the Mason et al. (2015b) UV LF model. In this model galaxy evolution is dependent on star formation that is associated with the construction of dark matter halos, with the assumption that these halos have a star formation efficiency that is mass dependent but redshift independent, which successfully reproduces observations over 13 Gyr (see also, e.g., Trenti et al. 2010; Tacchella et al. 2013, 2018; Mirocha et al. 2020).

The UV LF is plotted in Figure 2. It is well described by a Schechter (1976) function (Equation (1)). The steep drop-off of UV-bright galaxies can be explained in terms of rare high mass halos and their star formation efficiency: high mass halos are not efficient at forming stars, likely due to strong negative active galactic nuclei (AGN) feedback. The drop in number density with increasing redshift indicates a shift in star formation toward fainter, less massive galaxies.

Zoom In Zoom Out Reset image size

To obtain an accurate model of the Lyα LF to compare with observations, we must calibrate the Lyα LF by finding the normalization constant, F. This factor accounts for any over-prediction in the number density of LAEs caused by the Lyα luminosity distribution, P(L_α ∣ M _UV , x _Hɪ , z) (Dijkstra & Wyithe 2012; Gronke et al. 2015). If the Lyα luminosity distribution accurately describes the luminosities of the LBGs measured in the UV LF we should obtain F = 1.

We estimate F using a maximum-likelihood approach to fit our model at z = 5.7, x _Hɪ ∼ 0 to the Konno et al. (2018) and Ouchi et al. (2008) observations at z = 5.7. We set the calibration at this redshift and neutral fraction as it is likely to be after the end of reionization (e.g., McGreer et al. 2015). Note, due to our reionization simulation grid (see Section 2.1.1), we use x _Hɪ = 0.01 for the calibration, rather than x _Hɪ = 0.0, but note that the difference in P(L_α ∣M _UV , x _Hɪ ) should be negligible for such a small change in neutral fraction (as shown by Mason et al. 2018).

We maximize the likelihood for the observed Lyα LFs in each luminosity bin L_i : P(ϕ_i ∣ θ = F, L_i , σ_i ) given our model ${\phi }_{\mathrm{mod}}(\theta =F,L)$. This estimation using binned LFs may not be the most optimal: a more accurate likelihood would be obtained using individual source information and the survey selection function (see, e.g., Trenti & Stiavelli 2008; Kelly et al. 2008; Schmidt et al. 2014; Mason et al. 2015a); however, collating these data is not feasible within the scope of this project, and we leave this for future works. We note that Trenti & Stiavelli (2008) demonstrated that LFs estimated from binned data are generally in good agreement with those measured from unbinned data, but can bias the faint-end slope toward steeper values. However, as the observed high-redshift Lyα LFs are mostly ≳L*, we do not expect this to have a large impact on our results as the faint end will already have large uncertainties.

Following Gillet et al. (2020) we use a split-norm likelihood (Equation (9)), to take into account asymmetric error bars. Assuming each observation and luminosity bin are independent, the total likelihood is

$$\begin{eqnarray}&&P({\phi }_{\mathrm{obs}}\,| \,\theta =F,L)=\prod _{i}S({\phi }_{\mathrm{mod}}(F,{L}_{i}),{\mu }_{i},{\sigma }_{1,i},{\sigma }_{2,i}).\end{eqnarray} \tag{ 8 }$$

Here, ${\phi }_{\mathrm{mod}}(F,{L}_{i})$ is the model LF at luminosity L_i for parameter θ = F, μ_i is the observed number density value at L_i , and σ_1,i and σ_2,i are the respective lower and upper errors of the observed number density. In a single luminosity bin

$$\begin{eqnarray}S({\phi }_{\mathrm{mod}})=\left\{\begin{array}{ll}B\exp [-\displaystyle \frac{1}{2}\displaystyle \frac{{\left({\phi }_{\mathrm{mod}}-\mu \right)}^{2}}{{\sigma }_{1}^{2}}] & \mathrm{if}{\phi }_{\mathrm{mod}}\leqslant \mu \\ B\exp [-\displaystyle \frac{1}{2}\displaystyle \frac{{\left({\phi }_{\mathrm{mod}}-\mu \right)}^{2}}{{\sigma }_{2}^{2}}] & \mathrm{if}{\phi }_{\mathrm{mod}}\geqslant \mu \end{array}\right.,\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&B={\left[\sqrt{2\pi }\left(\displaystyle \frac{{\sigma }_{1}+{\sigma }_{2}}{2}\right)\right]}^{-1}.\end{eqnarray} \tag{ 10 }$$

We minimize the logarithm of the likelihood Equation (8) to find F using the Python package SciPy minimize. The obtained minimum value is F = 0.974, which we then use for Lyα LF at all redshifts and x _Hɪ values.

Our recovered value of F ≈ 1 indicates our luminosity distribution is a good model for LBGs. Other works, such as Dijkstra & Wyithe (2012) and Gronke et al. (2015), found F ∼ 0.5 at z = 5.7 and F < 1 at all lower redshifts, which is due to the different EW distribution they employed. We discuss this further in Section 4.3.

In Section 3.5, we use our model to infer the IGM neutral fraction from observations.

Bayes' theorem allows us to establish a posterior distribution for x _Hɪ given observations. Bayes' theorem is defined as

$$\begin{eqnarray}&&\begin{array}{l}P({x}_{{\rm\small{H\unicode{x0026A}}}}\,| \,\{{\phi }_{\mathrm{obs},i}({L}_{i})\},z)\\ \quad =\displaystyle \frac{P(\{{\phi }_{\mathrm{obs},i}({L}_{i})\}\,| \,{x}_{{\rm\small{H\unicode{x0026A}}}},z)P({x}_{{\rm\small{H\unicode{x0026A}}}}| z)}{P(\{{\phi }_{\mathrm{obs},i}({L}_{i})\})}.\end{array}\end{eqnarray} \tag{ 11 }$$

Here, {ϕ_obs,i (L_i )} is the set of observed data in luminosity bins L_i (where μ_i is the observed number density value at L_i , and σ_1,i and σ_2,i are the respective lower and upper errors of the observed number density). P({ϕ_obs,i (L_i )} ∣ x _Hɪ , z) is the likelihood of obtaining our observed data given the model. P(x _Hɪ ∣z) = P(x _Hɪ ) is the prior for the model parameter, x _Hɪ , where we assume the neutral fraction is independent of redshift. More physically, this prior could be dependent on redshift, but we leave the prior independent of redshift to allow more flexibility when estimating the neutral fraction. Regardless, we still see the inferred neutral fraction increase with redshift. We use a uniform prior from 0–1. Although this is technically not needed, making our approach essentially a maximum-likelihood estimate, we keep the Bayesian formalism to allow more physical priors to be used in future works. P({ϕ_obs,i (L_i )}) is the Bayesian information that normalizes the posterior distribution.

We obtain the posterior distribution of the neutral fraction of hydrogen using our Lyα LF model given the observed Lyα luminosity values, L_α , number density, ϕ(L_α ), and the uncertainties in the number density. We use the same split-norm likelihood defined in Equation (8) with ${\phi }_{\mathrm{mod}}({x}_{\mathrm{H\unicode{x0026A}}},{L}_{i})$. Further explanation of the inference of the neutral fraction can be seen in Appendix B.

We include uncertainties in our model Lyα LF due to uncertainties in the UV LF via a Monte Carlo approach. We generate 100 UV LFs with a 0.2 dex uncertainty in number density (estimated from the Mason et al. 2015b UV LF model). We then calculate the standard deviation of the resulting Lyα LF, ${\sigma }_{\mathrm{mod}}$, as a function of Lyα luminosity. We find the standard deviation is well described by ${\sigma }_{\mathrm{mod}}({L}_{\alpha })\approx 0.1\times {\phi }_{\mathrm{mod}}({L}_{\alpha })$. We use this uncertainty in calculating the likelihood (Equation (9)), where

$$\begin{eqnarray}&&{\sigma }_{1}\to \sqrt{{\sigma }_{1}^{2}+{\sigma }_{\mathrm{mod}}^{2}},\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&{\sigma }_{2}\to \sqrt{{\sigma }_{2}^{2}+{\sigma }_{\mathrm{mod}}^{2}}.\end{eqnarray} \tag{ 13 }$$

In comparing our model to observations, we wanted to ensure we used data sets where the selection strategies were similar to each other and similar to the data sets used to calibrate our model (Section 2.3), as it is known that different survey selection techniques can produce different estimates of the Lyα LF (for further discussion see Taylor et al. 2020). This led to the inclusion or exclusion of certain surveys from the estimation of the neutral fraction. In general, we aimed to use surveys that covered the widest areas (to minimize cosmic variance) and deepest Lyα luminosity limits.

For the neutral fraction inference (Section 3.5) it was important to use observed LFs that were calculated consistently with each other and our calibration LF at z = 5.7 (Section 2.3). As Konno et al. (2018) covers the largest area, we used their LF for our calibration. Therefore, for the neutral fraction inference, we included additional data sets that covered the largest redshift range with similar flux measurements and corrections for their systematic uncertainties.

The observational data sets we used to infer the neutral fraction, and their survey areas are as follows: Lyα LFs at z = 5.7 and 6.6 by Konno et al. (2018), who surveyed ∼13.8 and ∼21.2 deg² areas of the sky using Subaru/Hyper Suprime-Cam (HSC) Subaru Strategic Program Survey for redshifts z = 5.7 and 6.6, respectively, and by Ouchi et al. (2008, 2010), who surveyed a 1 deg² area of the sky using Subaru/XMM-Newton Deep Survey (SXDS) fields for both redshifts z = 5.7 and 6.6. At z = 7.0, we used Lyα LFs observed by Ota et al. (2017), who measured the total effective area of the Subaru Deep Field (SDF) and SXDS survey images for LF candidates to be ∼0.5 deg². Itoh et al. (2018) conducted an ultra-deep and large-area HSC imaging survey under the Cosmic HydrOgen Reionization Unveiled with Subaru Program in a total of 3.1 deg² using two independent blank fields. Finally, Hu et al. (2019) implemented a large-area survey using the Lyman Alpha Galaxies in the Epoch of Reionization project's deep-fields Cosmic Evolution Survey (COSMOS) and Chandra Deep Field South (CDFS) covering an effective area of 2.14 deg². Lyα LFs z = 7.3 are identified by Konno et al. (2014), who surveyed a ∼0.5 deg² area in the SXDS and COSMOS fields and Shibuya et al. (2012), who surveyed a total area of $1719\,{\mathrm{arcmin}}^{2}$ (∼0.5 deg²) in the SDF and the Subaru/XMM-Newton Deep Survey Field (SXDF), using the Suprime-Cam.

We ultimately excluded LFs measured by Santos et al. (2016) at z = 5.7 and 6.6 when estimating the neutral fraction because these LFs were significantly higher than those by Konno et al. (2018) and Ouchi et al. (2008, 2010). This is most likely due to differences in incompleteness corrections and the methodology for taking Lyα flux measurements from narrow-band images (see Santos et al. 2016; Konno et al. 2018, for further discussion). We also excluded the LFs measured by Taylor et al. (2020) at z = 6.6 from the estimation of the neutral fraction, where unlike other works, they corrected for a selection incompleteness; however, our model only incorporates observations uncorrected for selection incompleteness (this decision is discussed further in Section 3.2).

Although there are Lyα LF measurements at higher redshift values (e.g., Hibon et al. 2010; Tilvi et al. 2010; Krug et al. 2012; Clément et al. 2012; Matthee et al. 2014, at z = 7.7), we decide not to include these works in comparison to our model. Ultimately, the areas of surveys greater than z = 7.3 are much smaller than surveys completed at lower redshifts. Therefore, surveys at z > 7.3 more likely to be biased because reionization is inhomogeneous (see, e.g., Figure 11 from Jensen et al. 2013, where they compare LFs for different survey areas). The highest redshift LAE candidates are also prone to higher rates of contamination (Matthee et al. 2014), so with only the inclusion of lower redshifts, we can obtain more robust estimations of the neutral fraction.

To understand the impact of environment and galaxy properties on the evolution of the Lyα LF during reionzation, we investigate the typical Lyα luminosity of LBGs. In Figure 3, we plot the expectation value of Lyα luminosity, 〈L_α 〉, as a function of UV magnitude at z = 6.0. The expectation value is defined as

$$\begin{eqnarray}&&\langle {L}_{\alpha }\rangle ={\int }_{{L}_{\min }}^{{L}_{\max }}{L}_{\alpha }\cdot {P}_{\mathrm{norm}}({L}_{\alpha }\,| \,{M}_{{\rm\small{UV}}})\,{{dL}}_{\alpha },\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&{P}_{\mathrm{norm}}({L}_{\alpha }\,| \,{M}_{{\rm\small{UV}}})=\displaystyle \frac{P({L}_{\alpha }\,| \,{M}_{{\rm\small{UV}}})}{{\int }_{{L}_{\min }}^{{L}_{\max }}P({L}_{\alpha }\,| \,{M}_{{\rm\small{UV}}})},\end{eqnarray} \tag{ 15 }$$

where we calculate the integrals over the range 10³⁶ < L_α <10^44.5 erg s⁻¹. This range (i.e., a range greater than zero) is chosen because we want to observe the typical 〈L_α 〉–M _UV relation for Lyα emitters. We also want to ensure coverage of the Lyα luminosity values over the M _UV range − 24 ≤M _UV ≤ − 12.

Zoom In Zoom Out Reset image size

Figure 3 demonstrates that for UV-bright galaxies, M _UV ≲ − 20, we expect an average Lyα luminosity of L_α ∼ 10⁴²–10^43.6 erg s⁻¹. For UV-faint galaxies we expect a lower typical Lyα luminosity of L_α ∼ 10³⁹–10⁴¹ erg s⁻¹. Here, we also show how 〈L_α 〉 compares for galaxies brighter or fainter than M _UV * (where we use M _UV * = − 20.9 at z = 6.0 from Mason et al. 2015b).

Figure 3 shows a decrease in Lyα luminosity for a given M _UV as x _Hɪ increases, as expected due to the reduced transmission in an increasingly neutral IGM (Mason et al. 2018). This effect is strongest for UV-faint galaxies, where the average Lyα luminosity decreases by a factor of ∼10 as x _Hɪ increases to 1. UV-bright galaxies do not show much decrease in Lyα luminosity at different x _Hɪ . The more sizeable impact of reionization on UV-faint galaxies is because they typically exist in the outskirts of dense IGM environments. Thus, a more neutral IGM shifts their Lyα luminosity toward even lower values.

The bump in the plot, around M _UV ∼ − 20, is due to the EW probability distribution threshold between UV-bright and UV-faint galaxies (Mason et al. 2018). We tested the importance of this bump by fixing the P(EW) distribution (in our case we tested at P(EW∣M _UV ) = P(EW∣M _UV = − 17)), which removes the bump. However, this drastically affected the Lyα LF model where it does not fit observations well on the bright end. This means that P(EW) must be shifted to lower EW values for galaxies brighter than M _UV < − 20.

We compare our model Lyα LF to observations at z = 5.7, 6.6, 7.0, and 7.3. In Figure 4, we plot our Lyα LF models for a range of x _Hɪ from a fully neutral to fully ionized IGM at a given redshift. We also plot observations by Ouchi et al. (2008), Ouchi et al. (2010), Shibuya et al. (2012), Konno et al. (2014), Santos et al. (2016), Ota et al. (2017), Konno et al. (2018), Itoh et al. (2018), Hu et al. (2019), and Taylor et al. (2020) for comparison. Note that we use the selection incompleteness uncorrected LFs by Hu et al. (2019) for the best comparison with other observations and our model, which is calibrated using data that do not account for this incompleteness. Our simple model reproduces the shape of the Lyα LF remarkably well. Note that our Lyα LF is slightly lower than the one observed by Santos et al. (2016). This mismatch between their Lyα LF and other works found in literature is known and discussed, e.g., in Taylor et al. (2020) and Hu et al. (2019).

Zoom In Zoom Out Reset image size

We see that at z = 5.7, 6.6, and 7.0 the observations are fairly consistent with x _Hɪ ∼ 0.15–0.36, whereas at z = 7.3 the data are more consistent with x _Hɪ ∼ 0.66–0.87, suggesting an increasingly neutral IGM environment as redshift increases, consistent with other observations at z ∼ 7 (e.g., Mason et al. 2018; Whitler et al. 2020; Mason et al. 2019; Hoag et al. 2019).

Our model predicts that there is not much decrease in number density at low neutral fractions, x _Hɪ ≲ 0.4, but the LF decreases more rapidly at higher neutral fractions. Based on Figure 3 the lack of evolution at x _Hɪ ≲ 40% can be explained by the fact that the bright end of the Lyα LF is dominated by UV-bright galaxies, which exist in over-dense regions of IGM that tend to reionize early (Mesinger & Furlanetto 2007; Mesinger 2016; Harikane et al.2018). Thus, only in reionization's earliest stages do these galaxies experience a significant reduction in transmission.

We fit Schechter parameters, α, L*, ϕ*, for our models using emcee (Foreman-Mackey et al. 2013) to predict how the shape of the Lyα LF evolves with redshift and neutral fraction. To compare with observations, we fit the LF over the luminosity range $42.5\lt {\mathrm{log}}_{10}{L}_{\alpha }/\mathrm{erg}\,{{\rm{s}}}^{-1}\lt 44$. We fit the Schechter function to our Lyα LF models with all possible combinations of redshift and x _Hɪ but point out that the resulting Lyα LF are not exact Schechter functions for realistic Lyα EW distributions (even if the input UV LF is)—our Lyα LFs are generally less steep at the bright end than a Schechter function. Further details about the fitting are provided in Appendix A.

In Figure 5, we plot the evolution of each parameter (where we show the median value from the fits) with respect to redshift and x _Hɪ . We see that the parameters decrease overall as the universe becomes more neutral due to Lyα photon attenuation and decrease with redshift because galaxies become fainter and rarer as redshift increases (as seen in the UV LF (Figure 2)).

Zoom In Zoom Out Reset image size

The left panel of Figure 5 shows redshift versus α, which is the power-law slope for very low luminosities. At fixed x _Hɪ , α decreases as redshift increases, within the slope range of approximately −2.5 <α<−1.8, as expected due the hierarchical build-up of galaxies producing an increasingly steep faint-end slope of the UV LF with redshift (Mason et al. 2015b). The points plotted at each redshift show the impact of neutral hydrogen. α decreases significantly more due to the neutral gas than it does with redshift because Lyα attenuation affects faint galaxies more and thus makes them fainter, forcing them further back into the Lyα LF.

The center panel of Figure 5 reveals that L* decreases in the range z = 5–7, but increases sharply toward z = 8, and declines toward higher redshifts at fixed neutral fractions. This upturn is a consequence of the evolving shape of the UV LF due to dust attenuation at these redshifts in our model. In Mason et al. (2015b), there is an overlapping between z = 6 and 8 for the UV LF model around M _UV = − 23, which corresponds to L_α = 10^43.6 erg s⁻¹, also seen in observations (e.g., Bouwens et al. 2016, 2015). This overlapping is consistent with a reduction in dust obscuration, such that younger, brighter galaxies at higher redshifts contain less dust, and so, there is a possibility of observing more of them and shifting the LF models toward higher luminosities. As the neutral fraction increases at each redshift, the characteristic Lyα luminosity, L*, decreases. This trend can be attributed to an increasing attenuation of Lyα photons from UV-bright galaxies as the neutral fraction increases.

The right panel of Figure 5 shows a decreasing number density of Lyα emitting galaxies as we look back to higher redshifts. For each redshift, as shown, assuming the IGM is ionized, more Lyα emitting galaxies are expected to be visible to us and thus the number density increases with decreasing redshift, compared to a neutral IGM. At higher redshifts and at a fixed neutral fraction, we see an overall decrease in number density of Lyα emitters—due to the overall reduction in the number of galaxies at high redshifts (as seen in the evolution of the UV LF, e.g., Bouwens et al. 2015, 2016; Mason et al. 2015b).

We did not compare our Schechter function parameters directly with observations as previous works used a fixed α to determine their best-fit Schechter function parameters (e.g., Konno et al. 2014, 2018; Ouchi et al. 2008, 2010; Itoh et al. 2018; Ota et al. 2017; Hu et al. 2019). As the Schechter function parameters are degenerate (see, e.g., the discussion in Herenz et al. 2019), it is difficult to compare with our model directly. Also note that our model for the Lyα LF is not well described by a Schechter fit—our model LFs are typically less steep at the bright end than a Schechter function's exponential drop-off (see Figure 9). In our approach (Equation (16)) we do not expect the Lyα LF to be a Schechter form as the integral of the Schechter UV LF does not have a Schechter form, as discussed in Section 3.1 of Gronke et al. (2015). Physical reasons for an observed bright-end excess in the Lyα LF are discussed in Section 4.1 of Konno et al. (2018), including the contribution of AGN, large ionized bubbles around bright LAEs, and gravitational lensing.

The Lyα LD is the total energy emitted in Lyα by all galaxies obtained by integrating the LF:

$$\begin{eqnarray}&&{\rho }_{\alpha }(z,{x}_{\mathrm{H\unicode{x0026A}}})=\int {L}_{\alpha }\phi ({L}_{\alpha },z,{x}_{\mathrm{H\unicode{x0026A}}}){{dL}}_{\alpha },\end{eqnarray} \tag{ 16 }$$

where ϕ(L_α , z, x _Hɪ ) is our Lyα LF model number density. We generate ρ_α (z, x _Hɪ ) from our model by integrating over our luminosity grid 10^42.4 ≤ L_α ≤ 10^44.5 erg s⁻¹ to compare with observations that use similar luminosity limits, as the LD is highly sensitive to the minimum luminosity, due to the power-law slope of the LF faint end.

Figure 6 shows the evolution of the LD along a range of redshifts, z = 5–10, and the neutral fraction, x _Hɪ , predicted by our model. The modeled LD decreases by a factor of ∼10 as the universe becomes more neutral and declines overall as redshift increases. Observations from Konno et al. (2014), Konno et al. (2018), Hu et al. (2019), Itoh et al. (2018), and Ota et al. (2017) at z = 5.7, 6.6, 7.0, and 7.3 show a decrease in LD to higher redshifts. We compare our model with observations and see that the LD observations are consistent with a mostly ionized IGM at z = 5.7, 6.6, and 7.0, and increase toward a more neutral IGM at z = 7.3. These observations, similar to the Lyα LF observations, are chosen based on their selection incompleteness being uncorrected further explained in Hu et al. (2019). We also chose to plot Lyα LD observations that were considered fiducial (some works also tested separate LD points at different α to compare with others).

Zoom In Zoom Out Reset image size

We perform a Bayesian inference for the neutral fraction based on our model, as explained in Section 2.4. We infer x _Hɪ at z = 6.6, 7.0, and 7.3 by fitting our model to LF observations by Shibuya et al. (2012), Konno et al. (2014), Konno et al. (2018), Itoh et al. (2018), Ota et al. (2017), and Hu et al. (2019). We infer ${x}_{{\rm\small{H\unicode{x0026A}}}}(z=6.6)={0.08}_{-0.05}^{+0.08},{x}_{{\rm\small{H\unicode{x0026A}}}}(z=7.0)=0.28\pm 0.05$ and ${x}_{{\rm\small{H\unicode{x0026A}}}}(z=7.3)={0.83}_{-0.07}^{+0.06}$ (all errors are 1σ credible intervals). Appendix B shows the posterior distributions for x _Hɪ at each redshift (Figure 10).

We also show the comparisons between the posterior distributions for the neutral fraction obtained using the Lyα LD observations (Konno et al. 2014; Ota et al. 2017; Konno et al. 2018; Itoh et al. 2018; Hu et al. 2019, see Section 3.4). As expected, we see larger uncertainties in the estimations of the neutral fraction from the Lyα LD data due to including fewer data points compared to the LFs. We find ${x}_{{\rm\small{H\unicode{x0026A}}}}(z=6.6)={0.22}_{-0.11}^{+0.12}$, x _Hɪ (z = 7.0) = 0.25 ± 0.08, and ${x}_{{\rm\small{H\unicode{x0026A}}}}(z=7.3)={0.69}_{-0.11}^{+0.12}$. Using the full LF data thus not only enables us to infer neutral fractions that are more robust to nonuniform Lyα transmission, but adds a statistical advantage over previous LD methods by reducing the uncertainty on the neutral fraction (for further discussion, we refer the reader to Section 4.2).

Figure 7 shows our new constraints on reionzation history along with other approaches to estimating the neutral fraction. Our results show clear upward trend in neutral fraction at higher redshifts, consistent with an IGM that reionizes fairly rapidly. Our measurement at z = 6.6 is consistent with previous upper limits on the neutral fraction at z ≤ 6.6 (McGreer et al. 2015; Ouchi et al. 2010; Sobacchi & Mesinger 2015). Our measurement at z = 7.0 is consistent with inferences from the Lyα damping wing in the quasar ULAS J1120+0641 (Davies et al. 2018), but is lower than than constraints from the Lyα EW distribution in LBGs by Mason et al. (2018) and Whitler et al. (2020), though is still consistent within 2σ. Our measurement at z = 7.3 is consistent with other constraints at z > 7 (Hoag et al. 2019; Mason et al. 2019; Davies et al. 2018) though, like the other z > 7 constraints, is higher than the QSO damping wing measurements at z = 7.5 for the quasar ULAS J1342+0928 by Greig et al. (2019).

Zoom In Zoom Out Reset image size

The NGRST High Latitude Survey (HLS) and the Euclid Deep Field Survey (DFS) will both be particularly important surveys that will probe into higher redshifts and detect Lyα emission lines further into the Epoch of Reionization with wide-area slit-less spectroscopy. JWST will also be able to detect high-redshift Lyα with high sensitivity, albeit in smaller areas. Here we make predictions for potential Lyα LFs with these telescopes.

The Euclid DFS will cover a 40 deg² area at a 5σ flux limit of ∼8.6 × 10⁻¹⁷ erg s⁻¹ cm⁻². It will cover 0.9 ≤ λ_obs ≤1.3 μm corresponding to a redshift range for Lyα of 6 ≤ z ≤ 10 (e.g., Laureijs et al. 2012; Bagley et al. 2017). The NGRST HLS will cover 1.00 ≤ λ_obs ≤ 1.95 μm corresponding to a redshift range for Lyα of 8 ≤ z ≤ 15, and survey a 2200 deg² area at a 5σ flux limit of ∼7.1 × 10⁻¹⁷ erg s⁻¹ cm⁻² (e.g., Ryan et al. 2019; Spergel et al. 2013). JWST's Near Infrared Imager and Slitless Spectrograph (NIRISS) will cover 0.7 ≤ λ_obs ≤ 5.0 μm, capable of detecting Lyα at z ≳ 5. While there is no dedicated wide-area survey with NIRISS, we investigate a mock pure-parallel survey of 50 pointings (∼240 arcmin²) with a 5σ flux limit of ∼5.0 × 10⁻¹⁸ erg s⁻¹ cm⁻², assuming 2 hr exposures with the F115W filter.

We make predictions for these surveys in Figure 8. We plot our model Lyα LF and the approximate median neutral fraction value based the reionization history allowed by the cosmic microwave background (CMB) optical depth and dark pixel fraction at each redshift (Mason et al. 2019) between 6 < z < 10 (shown as the gray shaded region in Figure 7). We see that these surveys will detect luminous Lyα emitters at higher redshifts. We predict the NGRST HLS will be able to discover galaxies L_α > 10^43.76 erg s⁻¹ at redshifts up to z = 10. The Euclid DFS will be able to detect bright galaxies at L_α > 10^43.84 erg s⁻¹ but only up to z ∼ 8. Using the JWST mock pure-parallel survey, we estimate it be able to detect bright galaxies at L_α > 10^42.61 erg s⁻¹ up to z ∼ 9–10.

Zoom In Zoom Out Reset image size

We note that the predicted number counts will likely be higher than shown in Figure 8 due to gravitational lensing magnification bias, which can increase the observed number of galaxies at the bright end of the LF in flux-limited surveys (Wyithe et al. 2011; Mason et al. 2015a; Marchetti et al. 2017).

In building our model, we make several assumptions that can affect the results, which are summarized here. These assumptions are also described by Mason et al. (2018) and Whitler et al. (2020), and we refer the reader there for additional details. First, we assume that the intrinsic emitted Lyα luminosity distribution does not evolve with redshift (only M _UV ) and is the same as the observed Lyα luminosity distribution at z ∼ 6 (as modeled from the EW distribution by Mason et al. 2018), and that only evolution in the observed luminosity distribution is due to reionization alone. However, we should expect some evolution in the Lyα luminosity distribution with redshift, as galaxy properties evolve. Physically, we could expect galaxies to have higher luminosities as redshift increases, due to, e.g., lower dust attenuation (Hayes et al. 2011), which leads to a steepening of the Lyα LF with redshift (Gronke et al. 2015; Dressler et al. 2015). In that case, more significant absorption by the IGM would be required to explain the observed Lyα LFs and we would thus infer a higher neutral fraction. Alternatively, a decrease in outflow velocities possibly associated with a decreasing specific star formation rate could lower Lyα escape from galaxies, decreasing the emitted luminosity (Hassan & Gronke 2021), in which case a lower neutral fraction would be inferred.

We also assume that Lyα visibility evolution between z = 6 and 7 is due only to the evolution of the Lyα damping wing optical depth (e.g., Miralda-Escudé 1998), due to the diffuse neutral IGM. We do not model redshift evolution of the Lyα transmission in the ionized IGM or circumgalactic medium (CGM) at fixed halo mass (e.g., Laursen et al. 2011; Weinberger et al. 2018). The amount of transmission due to these components is determined by the Lyα line shape. If the Lyα line velocity offset from systemic decreases with redshift at fixed halo mass, this would decrease the Lyα transmission throughout the ionized IGM and CGM (Dijkstra et al. 2011; Choudhury et al. 2015) and reduce the need for a highly neutral IGM. A full exploration of the degeneracies and systematic uncertainties due to the Lyα emission model is left to future work.

In this work we directly compare measurements of the Lyα LF to our model. Previous works, e.g., Ouchi et al. (2010), Zheng et al. (2017), Konno et al. (2014), Konno et al. (2018), Inoue et al. (2018), and Hu et al. (2019), estimate the neutral fraction by evaluating the Lyα luminosity density. While this provides a reasonable first-order estimate, this method can be difficult to interpret as it relies on models and observations using the same luminosity limit in the luminosity density integral, and it collapses any valuable information that is obtained in the evolving shape of the Lyα LF (including increasing the statistical uncertainty on x _Hɪ by reducing the number of data points—see Appendix B). As demonstrated in Section 3.3, the Lyα shape does evolve.

Furthermore, these works estimate the neutral fraction based on the assumption that the transmission fraction, T_IGM, is constant for all galaxies. However, as shown by, e.g., Mason et al. (2018) and Whitler et al. (2020) it is a broad distribution and depends on galaxy properties, through their large scale structure environment and the internal kinematics that sets the Lyα line shape. If the transmission fraction varies with Lyα or UV luminosity, the common method of calculating Lyα transmission by taking the ratio of the Lyα luminosity density at different redshifts is invalid, thus our work enables a more robust estimate of x _Hɪ from the Lyα LFs. Many of these works compare with simulations by McQuinn et al. (2007), which do model the impact of inhomogeneous reionization on the Lyα LF but take a more simplistic approach to modeling Lyα luminosity: each galaxy has a Lyα luminosity proportional to its mass. In our work, we use the EW probability distribution from Mason et al. (2015b) to take into account that galaxies have a range of Lyα EW at fixed UV magnitude. Our approach is thus more similar to that of Jensen et al. (2013), who model the Lyα luminosity probability distribution as a function of halo mass. However, Jensen et al. (2013) model the Lyα luminosity and EW distributions independently, whereas we have shown the Lyα LF can be self-consistently described by the same Lyα EW distribution that describes LBGs.

Finally, our semi-analytic model provides flexibility over approaches that model radiative transfer in N-body simulations (e.g., Dayal et al. 2011; Jensen et al. 2013; Hutter et al. 2014; Inoue et al. 2018) or sophisticated hydrodynamical simulations (e.g., Dayal et al. 2011; Weinberger et al. 2019) by keeping the IGM neutral fraction and redshift as free parameters, rather than assuming a fixed reionization history. Using this new model for the Lyα LF, we can separate redshift and the neutral fraction, fixing either parameter if needed, and see how observations compare to our model.

As described in Section 2.3, a normalization factor, F, is introduced to the Lyα LF to account for any mismatch in the number density of LAEs (Dijkstra & Wyithe 2012). If the Lyα luminosity distribution model accurately describes the LBG population observed in the UV LF, we expect F = 1. We find F = 0.974 in our model, which is considerably higher than previous work focusing on the Lyα LF at lower redshifts, which found F ∼ 0.5 (Dijkstra & Wyithe 2012; Gronke et al. 2015).

The key difference compared to this previous work that leads to our model successfully reproducing the Lyα LF without the need for additional normalization is due to the EW distribution we employed. While we, as well as the previous work, include non-emitters in the model, the EW distribution used by Dijkstra & Wyithe (2012), Gronke et al. (2015) (calibrated to the measurements of Shapley et al. 2003; Stark et al. 2010, 2011, at z ∼ 3–6) shows an increasing probability of Lyα emission for lower UV brightness galaxies, leading to a Lyα emitter (EW > 0 Å) fraction of unity for M _UV ≳ − 19. The EW distribution we used caps this emitter fraction at 65% (our function A(M _UV ) in Section 2.1.1).

Which parameterization of the EW distribution is most appropriate is still up for debate (see, e.g., the detailed discussion by Oyarzún et al. 2017, who study more complex distributions also dependent on the stellar mass and the UV slope). The recent study of Caruana et al. (2018) supports our non-emitter fraction, as they find a fraction of ∼0.5 ± 0.15 galaxies with EW > 0 Å at 3 < z < 6 in Hubble Space Telescope (HST) continuum selected galaxies for within the Multi Unit Spectroscopic Explorer (MUSE) wide field (Herenz et al. 2017; Urrutia et al. 2019) is present. However, Caruana et al. (2018) also find a non-evolution of this fraction with UV magnitude (in contrast with previous models) as well as typically lower Lyα fractions for larger EW cuts (e.g., the value for W > 50 Å seems to be in slight tension with measurements by Stark et al. (2010), who find ≈45% ± 10% of galaxies have Lyα EW > 55 Å and a strong anticorrelation with UV brightness).

Another factor to consider when comparing z ≲ 3 and z ≳ 5 EW distributions is the impact of the IGM on the Lyα line even at z ∼ 5–6. In fact, it is has been suggested by Weinberger et al. (2019) that F could stem from this effect but note that Dijkstra & Wyithe (2012) and Gronke et al. (2015) compare their modeled Lyα LFs at z ∼ 3 to data by Ouchi et al. (2008); Gronwall et al. (2007) and Rauch et al. (2008), respectively, and still require F ∼ 0.5 at this low z.

Future observational studies will constrain the Lyα EW distribution further both as a function of redshift and UV magnitude, and thus can quantify the fraction of non-emitters for UV-faint galaxies. We have shown that a constant non-emitter fraction of ∼35% for M_UV ≲ − 19 makes the fudge factor F ∼ 0.5 obsolete, which indicates that such a cutoff exists in reality.