The Universe Is Reionizing at z ∼ 7: Bayesian Inference of the IGM Neutral Fraction Using Lyα Emission from Galaxies

Lyα photons are produced predominantly via recombination in ${\rm{H}}\,{\rm{II}}$ regions around young stars and have a high cross-section for resonant scattering (for a detailed review, see Dijkstra 2014). As the ISM of individual galaxies contains a large amount of neutral hydrogen gas to escape the ISM, Lyα photons must diffuse both spatially and spectrally (e.g., Shapley et al. 2003; McLinden et al. 2011; Chonis et al. 2013; Mostardi et al. 2013; Song et al. 2014). This produces the fiducial double-peaked Lyα lineshape, for which the red (blue) peak is enhanced for outflows (inflows; Zheng & Miralda-Escude 2002; Verhamme et al. 2006).

In this work, we model the Lyα lineshape after transmission through the ISM as a Gaussian, centered at a velocity offset Δv from the systemic redshift of the galaxy (due to scattering through the ISM, described in Section 2.1.1) with a velocity dispersion σ_α (due to scattering and thermal broadening in the ISM, described in Section 2.1.2). We refer to this lineshape as "intrinsic"; examples are shown as dotted black lines in Figure 1. As described in Section 2.2, even after reionization residual neutral gas in the IGM and CGM will absorb all blue flux at z ≳ 6.

Zoom In Zoom Out Reset image size

2.1.1. Modeling Lyα Velocity Offsets

Numerous studies of star-forming galaxies at z ≲ 4 have identified the column density of neutral hydrogen (N_hi) within the ISM as a key mediator of Lyα radiative transfer. Lyα photons traveling through highly dense neutral ISM scatter more frequently and emerge with larger velocity offsets than galaxies with lower N_hi (Shibuya et al. 2014; Hashimoto et al. 2015; Yang et al. 2016, 2017; Guaita et al. 2017).

Low mass galaxies, especially at high redshifts, are less likely to contain significant fractions of neutral gas due to enhanced photoionization feedback. Additionally, strong star formation feedback may drive outflows and/or reduce the covering fraction of neutral gas in the ISM, which can facilitate Lyα escape (Jones et al. 2013; Trainor et al. 2015; Leethochawalit et al. 2016).

Recently, a correlation has been suggested between UV magnitude and Lyα velocity offset (Schenker et al. 2013a; Erb et al. 2014; Song et al. 2014; Stark et al. 2015, 2017; Mainali et al. 2017), again indicating that galaxy mass and/or SFR strongly affects Lyα escape. However, galaxies with the same UV magnitudes at different redshifts likely have very different masses because of increasing SFR at high redshift (e.g., Behroozi et al. 2013b; Barone-Nugent et al. 2014; Finkelstein et al. 2015b; Mason et al. 2015; Harikane et al. 2016, 2017), so one should be cautious of comparing galaxies with the same UV magnitudes at different redshifts. We plot a compilation of ${M}_{\mathrm{UV}}$–Δv measurements from the literature (Erb et al. 2014; Steidel et al. 2014; Stark et al. 2015, 2017; Willott et al. 2015; Inoue et al. 2016; Pentericci et al. 2016; Bradač et al. 2017; Mainali et al. 2017) in Figure 2 (left), where it is clear the high redshift galaxies have lower Δv at given ${M}_{\mathrm{UV}}$ compared to the low redshift galaxies, probably because they have lower mass.

Zoom In Zoom Out Reset image size

To model the effect of the ISM on Lyα escape, we assume the column density of neutral hydrogen within the ISM is determined by halo mass and is the most important quantity for understanding the emerging Lyα line profile. This is likely an over-simplification—for example, "shell" models take ∼6 parameters to model Lyα lines (Verhamme et al. 2006, 2008; Gronke & Dijkstra 2016)—but is an efficient first-order approach. With this in mind, we assume a correlation between Δv and halo mass of the form ${\rm{\Delta }}{\text{}}v\sim {({M}_{h})}^{m}$, where we determine m empirically from observations, as described below.

We take a sample of 158 z ∼ 2–3 galaxies with both UV magnitudes and Lyα velocity offsets from Erb et al. (2014) and Steidel et al. (2014). The Steidel et al. (2014) sample (from the KBSS-MOSFIRE survey) is effectively complete at ${M}_{\mathrm{UV}}$ ≲ −20, where ∼90% of their photometrically selected LBGs have rest-frame optical emission lines detected in deep near-IR spectroscopy with Keck/MOSFIRE (McLean et al. 2012). The Erb et al. (2014) sample comprises 36 galaxies selected as Lyα emitters in narrow-band photometry with −21 ≲ ${M}_{\mathrm{UV}}$ ≲ −18; all these objects had rest-frame optical lines detected in MOSFIRE observations. We note the Erb et al. (2014) sample does not include faint Lyα emitters (W ≲ 25 Å), which may have higher velocity offsets, given observed anti-correlations between Lyα EW and Δv (Hashimoto et al. 2013; Shibuya et al. 2014; Erb et al. 2014). While these samples are the largest available to measure an ${M}_{\mathrm{UV}}$–Δv correlation, future rest-frame optical follow-up of large samples of galaxies with detected Lyα emission (e.g., from the HETDEX and MUSE-Wide spectroscopic surveys; Song et al. 2014; Herenz et al. 2017) will provide more complete information about the relationships between Lyα radiative transfer in the ISM and galaxy properties.

As described above, it is difficult to directly compare galaxies at fixed UV magnitude across cosmic time, so we map UV magnitudes to halo mass. To first order, the depth of a halo's gravitational potential well is the dominant influence on galaxy properties independent of redshift (Behroozi et al. 2013a; Mason et al. 2015; Moster et al. 2017). We assume no redshift evolution between halo mass and velocity offset. We convert UV magnitude to halo mass using the successful model derived by Mason et al. (2015), which assumes the SFR is proportional to the halo mass assembly rate at a given halo mass and redshift, and is consistent with ${M}_{\mathrm{UV}}$–M_h measurements from clustering at z ∼ 7 (Barone-Nugent et al. 2014; Harikane et al. 2016, 2017). We add a 0.2 mag uncertainty to the UV magnitude measurements to account for scatter in the UV magnitude–halo mass relation (e.g., Finkelstein et al. 2015b) and use the propagated uncertainties in halo mass in Δv − M_h.

In the right panel of Figure 2 we plot the literature Δv measurements as a function of the estimated halo masses. Due to the uncertainties in mapping from UV magnitude to halo mass for very bright galaxies at z ≲ 4, which may be significantly more starbursty than average, we discard the z ∼ 2 galaxies with ${M}_{\mathrm{UV}}$ < −21 from further analysis. Likewise, we exclude from this inference the galaxies at z ∼ 7 with ${M}_{\mathrm{UV}}$ < −22, deferring their analysis to a later paper (Mason et al. 2018).

When we transform to halo mass, the high redshift literature points clearly lie within the low redshift data space. This suggests halo mass is a useful approximately redshift independent indicator of Lyα escape routes. We note gravitationally lensed objects at intermediate redshifts suggest these trends hold at low mass/luminosity (e.g., a lensed ${M}_{\mathrm{UV}}$ = −17 galaxy at z ∼ 3 was recently observed with a Lyα velocity offset of 51 km s⁻¹; Vanzella et al. 2016). Further studies, using NIRSpec on JWST, will be able to investigate these trends at high redshifts.

The distribution is well-described by a log-normal distribution with a peak that increases with increasing luminosity, and approximately constant variance:

$$\begin{eqnarray}&&p({\rm{\Delta }}{\text{}}v\,| \,{M}_{h})=\displaystyle \frac{\exp \left[-\tfrac{{({\mathrm{log}}_{10}{\rm{\Delta }}{\text{}}v-V({M}_{h}))}^{2}}{2{\sigma }_{v}^{2}}\right]}{{\rm{\Delta }}{\text{}}v\mathrm{ln}10\sqrt{2\pi }{\sigma }_{v}},\end{eqnarray} \tag{ 1 }$$

where V is a linear relation corresponding to the most likely ${\mathrm{log}}_{10}({\rm{\Delta }}{\text{}}v)$ at a given halo mass,

$$\begin{eqnarray}&&V({M}_{h})=m{\mathrm{log}}_{10}\left(\displaystyle \frac{{M}_{h}}{1.55\times {10}^{12}{M}_{\odot }}\right)+c.\end{eqnarray} \tag{ 2 }$$

To find the parameters m, c, and σ_v, we take Equation (1) as the likelihood function and perform a Bayesian inference on the z ∼ 2 galaxies with ${M}_{\mathrm{UV}}$ > −21, with uniform priors on the parameters. The inferred parameters are m = 0.32 ± 0.07, c = 2.48 ± 0.03, and σ_v = 0.24 ± 0.02. We show this relation in Figure 2.

We can obtain an approximate relation between velocity offset, UV magnitude, and redshift by approximating the Mason et al. (2015) UV magnitude–halo mass relation as broken linear: ${\mathrm{log}}_{10}{M}_{h}[{M}_{\odot }]\approx \gamma ({M}_{\mathrm{UV}}+20.0+0.26z)+11.75$, where γ = −0.3 for ${M}_{\mathrm{UV}}\geqslant -20.0-0.26z$, and γ = −0.7 otherwise. The mean velocity offset in km s⁻¹ can then be approximated as

$$\begin{eqnarray}&&{\mathrm{log}}_{10}{\rm{\Delta }}{\text{}}v({M}_{\mathrm{UV}},z)\approx 0.32\gamma ({M}_{\mathrm{UV}}+20.0+0.26z)+2.34.\end{eqnarray} \tag{ 3 }$$

In this work we sample directly from the distribution in Equation (2) to calculate velocity offsets directly for simulated halos, including scatter.

In Figure 2 we also plot the circular velocity (${v}_{c}={[10{{GM}}_{h}H(z)]}^{1/3}$) at z = 2 and z = 7 for comparison with the observed data. The circular velocities are higher at low redshifts, as halos are less dense and more extended, but there is a clear similarity in our derived trend ${\rm{\Delta }}{\text{}}v\sim {M}_{h}^{0.32}$ and the circular velocity ${v}_{c}\sim {M}_{h}^{1/3}$. Investigating these trends with larger samples at low redshifts with dynamical mass measurements (e.g., Erb et al. 2014; Trainor et al. 2015) could determine to what extent Lyα radiative transfer depends on the gravitational potential of the halo.

2.1.2. Modeling Lyα Line Widths

2.1.3. EW Distribution in an Ionized Universe

The key observable of Lyα emission lines at high redshift is their equivalent width (EW or W), a measure of the brightness of the emission line relative to the UV continuum. As Lyα photons from high redshift galaxies are attenuated by neutral gas in the intervening CGM and IGM, we observed only a fraction, ${{ \mathcal T }}_{\mathrm{IGM}}$ (the Lyα transmission fraction), of the emitted EW (i.e., ${W}_{\mathrm{obs}}={W}_{\mathrm{em}}\times {{ \mathcal T }}_{\mathrm{IGM}}$), where W_em is the emitted equivalent width without any damping due to reionization.

In this work we consider the differential evolution of Lyα equivalent widths between z ∼ 6 and z ∼ 7, and assume the distribution of equivalent widths changes only because of the increasing neutrality of the IGM due to reionization. This is likely a simplification, as trends at lower redshifts show increasing EW with redshift as dust decreases in galaxies (Hayes et al. 2011), but the time between z ∼ 6 and z ∼ 7 is short (<200 Myr). If the underlying EW distribution does evolve significantly during that time, it will likely be to increase the emitted EW (due to decreasing dust; Hayes et al. 2011); thus the reduction due to reionization would need to be greater to match the observed EW distribution at z ∼ 7 (Dijkstra et al. 2011). Papovich et al. (2011) suggests there may be an increase in gas reservoirs with increasing redshifts due to rapid accretion rates, which could also reduce the emission of strong Lyα.

Thus observed equivalent widths at z ∼ 7 are ${W}_{7}={W}_{6}\,\times {{{ \mathcal T }}_{\mathrm{IGM}}}_{,7}/{{{ \mathcal T }}_{\mathrm{IGM}}}_{,6}$, where ${{{ \mathcal T }}_{\mathrm{IGM}}}_{,z}$ is the transmission fraction of Lyα emission for a single object at redshift z. In Section 2.2 we describe the calculation of transmission fractions along thousands of lines of sight using state-of-the-art cosmological reionization simulations.

A key input to the model, then, is the z ∼ 6 distribution of EW as a function of galaxy properties. Lyα EWs for UV continuum-selected galaxies have an observed equivalent width distribution with a peak at zero and some tail to high EW—which is usually parameterized as an exponential function (Dijkstra & Wyithe 2012), log-normal (Schenker et al. 2014), or truncated normal distribution plus a delta function (Treu et al. 2012). At z ≲ 2, where large samples exist, including the local "Green Peas," Lyα EW is observed to anti-correlate strongly with UV luminosity (Shapley et al. 2003; Stark et al. 2011; Hashimoto et al. 2013) SFR (Verhamme et al. 2008), H i covering fraction (Shibuya et al. 2014), and Lyα escape fraction (Yang et al. 2017), all indicating Lyα photons are significantly absorbed by neutral hydrogen gas and dust inside the ISM of massive, highly star-forming galaxies (e.g., Verhamme et al. 2008; Erb et al. 2014; Yang et al. 2017). At high redshift, the Lyα EW distribution is usually parameterized as a conditional probability of $p(W\,| \,{M}_{\mathrm{UV}})$ (Dijkstra & Wyithe 2012; Treu et al. 2012), though dependence on UV spectral slope β has also been considered (Schenker et al. 2014).

We take the z ∼ 6 EW distribution from De Barros et al. (2017) and L. Pentericci et al. (2018, in preparation) from a Large Program with VLT/FORS2. This sample contains 127 objects, with UV magnitudes between −22.5 ≲ ${M}_{\mathrm{UV}}$ ≲ −17.5, of which 63% have Lyα detections. We parameterize it as an exponential distribution plus a delta function:

$$\begin{eqnarray}{p}_{6}(W\,| \,{M}_{\mathrm{UV}}) & = & \displaystyle \frac{A({M}_{\mathrm{UV}})}{{W}_{c}({M}_{\mathrm{UV}})}{e}^{-\tfrac{W}{{W}_{c}({M}_{\mathrm{UV}})}}H(W)\\ & & +[1-A({M}_{\mathrm{UV}})]\delta (W),\end{eqnarray} \tag{ 4 }$$

where A and W_c account for the fraction of non-emitters, and for the anti-correlation of EW with ${M}_{\mathrm{UV}}$. H(W) is the Heaviside step function, and δ(W) is a Dirac delta function. A implicitly includes contamination by low redshift interlopers in the photometric selection (the interloper fraction is ≤29% for this sample, assuming all non-detections were low redshift contaminants; De Barros et al. 2017)—that is, we do not distinguish between non-emitters at z ∼ 6 and low redshift contaminants when accounting for non-detections in fitting the parameters (see below). Within our framework, this means we assume a similarly small total interloper and non-emitter fraction at z ∼ 7. Recent work by Vulcani et al. (2017) supports this assumption: they found comparably low contamination fractions at z ∼ 6 and z ∼ 7 in an evaluation of photometric selections.

To find these parameters, we divided the sample into three bins: ${M}_{\mathrm{UV}}$ ≤ −21; −21 < ${M}_{\mathrm{UV}}$ < −20; and ${M}_{\mathrm{UV}}$ ≥ −20. We used Equation (4) as a likelihood (${p}_{6}(W| \,A,{W}_{c})$) and performed a Bayesian inference to infer A and W_c for each bin, similar to the methods of Oyarzún et al. (2017), using uniform priors with 0 < A < 1 and 0 < W_c < 100. In the inference, we fully account for the the non-detections of Lyα (using ${p}_{6}(W\lt {W}_{\mathrm{lim}}\,| \,A,{W}_{c})$ as the likelihood given an EW limit W_lim, in the same way as described in more detail in Section 3). The uncertainties and EW limits calculated by De Barros et al. (2017) are obtained using simulations that fully account for incompleteness and wavelength sensitivity. We note that in this framework the EW limits are a conservative minimum that could be measured over the entire wavelength range (see also Section 4.2); future work could incorporate the full wavelength-dependent line flux sensitivity. To allow these parameters to smoothly vary with magnitude in the range of −21 < ${M}_{\mathrm{UV}}$ < −20, we use a hyperbolic tangent function to connect our inferred parameters, without extrapolating beyond the range of the data. We find $A=0.65\,+0.1\tanh [3({M}_{\mathrm{UV}}+20.75)]$ and ${W}_{c}=31+12\tanh [4({M}_{\mathrm{UV}}\,+20.25)]$Å from fitting to the data. A and W_c vary smoothly with magnitude.

We choose this exponential parameterization of the data because it gives a good description of the data and is easy to treat analytically, and has previously been shown to be an excellent fit to Lyα EW PDFs (e.g., Oyarzún et al. 2017). We do not include uncertainties in these parameters, and we note the parameterization of ${p}_{6}(W| \,{M}_{\mathrm{UV}})$ is fairly arbitrary but does not qualitatively affect Lyα modeling during the EoR (Treu et al. 2012; Gronke et al. 2015b). Indeed we get the same results, within the uncertainties, if we use the p₆(W) truncated Gaussian distribution from Treu et al. (2012) based on the sample presented in Stark et al. (2011).

Example PDFs given by Equation (4) are plotted in Figure 3 for two values of ${M}_{\mathrm{UV}}$. We show both the intrinsic PDF and the distribution convolved with a 5 Å typical measurement error which introduces at "bump" around W = 0, where the underlying distribution is a delta function. We also show histograms of the EW observations of De Barros et al. (2017) and L. Pentericci et al. (2018, in preparation) in two bins corresponding to UV bright and faint LBGs. As shown by, for example, Verhamme et al. (2008), Stark et al. (2010), and Oyarzún et al. (2017), Lyα EW strongly depends on UV magnitude.

Zoom In Zoom Out Reset image size

A Lyα emission line is significantly attenuated by the CGM and IGM as its photons redshift into resonance with abundant neutral hydrogen along the line of sight. Effectively, for a Lyα line at z ≳ 6, all photons emitted blueward of the Lyα resonance (1216 Å) are absorbed by the IGM, as even after reionization there is still is a fraction of neutral hydrogen within ${\rm{H}}\,{\rm{II}}$ regions (Gunn & Peterson 1965). Infalling overdense gas around halos can also increase the opacity of the IGM near the Lyα resonance and onto the red side of the Lyα line (Santos 2004; Dijkstra et al. 2007; Laursen et al. 2011).

For simplicity we assume all Lyα photons emitted below the circular velocity of a halo are absorbed in the CGM, and all redder photons are transmitted (Dijkstra et al. 2011; Laursen et al. 2011). This treatment of the CGM may be crude, but it enables us to investigate the relative difference between observations at z ∼ 6 and z ∼ 7, assuming any evidence of a difference is driven by reionization. Future work could incorporate more complex modeling of CGM absorption (e.g., Kakiichi & Dijkstra 2017). Figure 1 shows example model Lyα emission lines, where the dotted black lines correspond to the intrinsic line profile after transmission through the ISM and the black solid lines correspond to the lineshape after resonant absorption in the CGM/IGM, which absorbs the flux blueward of v_circ.

During reionization, there is an additional opacity to Lyα caused by the presence of cosmic diffuse neutral hydrogen patches, which attenuate the damping wing of the Lyα line cross-section (Miralda-Escude 1998). The transmission of Lyα photons through the reionizing IGM is driven by the global fraction of neutral hydrogen, ${\overline{x}}_{{\rm{H}}{\rm{I}}}$(z).

Thus the total opacity to Lyα due to neutral hydrogen within the IGM is given by

$$\begin{eqnarray}&&{\tau }_{{\rm{IGM}}}(z,v)={\tau }_{{\rm{d}}}(z,v)+{\tau }_{{\rm{H}}{\rm{II}}}(z,v),\end{eqnarray} \tag{ 5 }$$

where ${\tau }_{{\rm{d}}}(z,v)$ is the damping wing optical depth which is only present during the EoR, and ${\tau }_{{\rm{H}}{\rm{II}}}$ (z, v) is the optical depth due to resonant absorption within the CGM of galaxies (infalling gas) and any neutral hydrogen within the local ${\rm{H}}\,{\rm{II}}$ region of a galaxy. For simplicity, we assume ${e}^{-{\tau }_{{\rm{H}}{\rm{II}}}}=H(v-{v}_{\mathrm{circ}})$ at both z ∼ 6 and z ∼ 7.

In this model, we assume the universe is fully ionized at z ∼ 6; thus the damping wing opacity only becomes important at z > 6. This may not be exactly the case, but current constraints on ${\overline{x}}_{{\rm{H}}{\rm{I}}}$ at z ∼ 6 suggest the neutral fraction is low (${\overline{x}}_{{\rm{H}}{\rm{I}}}$ < 0.1; McGreer et al. 2014), so the reionization effect on Lyα emission will be small.

To obtain the damping wing optical depths, ${\tau }_{{\rm{d}}}(z=7,v)$ requires a model of the IGM topology during reionization. While observation papers of Lyα emission with reionization inferences have used simple "patchy" or "smooth" IGM topologies (Treu et al. 2012, 2013; Pentericci et al. 2014; Tilvi et al. 2014), for this work, we consider realistic reionization topologies from state-of-the-art theoretical modeling. We obtain Lyα damping optical depths from the public Evolution of 21 cm Structure (EoS) suite of reionization simulations described by Mesinger et al. (2015, 2016).⁹

Due to the strong clustering of the first galaxies, spatial fluctuations in the IGM neutral fraction during reionization existed on scales of tens of Mpcs. Accurately modeling these fluctuations and the growth of ionized ${\rm{H}}\,{\rm{II}}$ bubbles in the IGM requires cosmological simulations at least 100 Mpc in size (Trac & Gnedin 2011; Mesinger et al. 2015). The EoS reionization simulations use 21cmfastv2 (Sobacchi & Mesinger 2014), where inhomogeneous recombinations and ionizations in the IGM are treated at a sub-grid level on a density field in a box with sides 1.6 Gpc with a resolution 1024³. The simulations produce ${\overline{x}}_{{\rm{H}}{\rm{I}}}$ maps at different redshifts and superimpose them on the z ∼ 7 halo field to produce cubes of the z ∼ 7 IGM for a range of neutral fractions. For the bulk of reionization, this is analogous to changing the ionization efficiency at a fixed redshift (e.g., McQuinn et al. 2007b; Mesinger & Furlanetto 2008b).

The timeline and topology of reionization is determined by the mass of galaxies that dominate reionization and the ionization efficiency, $\zeta \propto {f}_{\mathrm{esc}}\times {f}_{\star }$, where f_esc is the fraction of ionizing photons which escape galaxies into the IGM, and f_⋆ is the stellar mass fraction in galaxies. As both of these parameters are expected to scale with halo mass in complementary ways (i.e., low mass halos host galaxies with a low stellar mass fraction and high escape fraction; e.g., Kimm et al. 2017; Trebitsch et al. 2017), over the relevant range of halo masses that host galaxies that dominate reionization (M_h ≲ 10¹¹ M_⊙; e.g., Kimm et al. 2017), ζ is assumed to be constant in the EoS simulations. The simulations use a free parameter that sets the minimum mass of halos capable of hosting star formation, and then adjust ζ to produce a Thompson scattering optical depth to the CMB consistent with the Planck Collaboration (2015).

We use the fiducial "Faint Galaxies" run, which corresponds the primary drivers of reionization (being low mass star-forming galaxies) with an atomic cooling threshold of T_vir ≳ 10⁴ K, with ζ = 20, producing IGM morphologies characterized by small HII regions. While the EoS simulations have another run, "Bright Galaxies," where reionization is driven by more massive galaxies, producing larger HII regions, it has been shown that information about the reionization morphology is mostly smeared out when using galaxies spread in redshift (Δz ≳ 0.1 bin; e.g., Sobacchi & Mesinger 2015, though with large spectroscopic samples, Δz ∼ 0.01, the sensitivity does increase), as is the case for our sample (see Section 4.2), so we do not expect a significant change in our results if we were to use an alternative run. Indeed, Mason et al. (2018) uses both simulation runs but shows that the transmission of Lyα from galaxies M_h ∼ 10¹⁰–10¹² is relatively independent of the reionization morphology. Similarly, Greig et al. (2017) show QSO damping wing effects are not particularly sensitive to the reionization morphology.

Halos are located in the same density field as the IGM simulation. We ignore absorption from Damped Lyα Absorbers (DLAs) inside the cosmic ${\rm{H}}\,{\rm{II}}$ regions (Bolton & Haehnelt 2013), which has been shown to have a minor impact on the Lyα fraction when self-shielding is calculated more accurately (Mesinger et al. 2015). We refer the reader to Mesinger et al. (2016) for more details of the simulation. For this work we focus on z ∼ 7, where large samples of LBGs have spectroscopic follow-up (Pentericci et al. 2014; Schmidt et al. 2016), but it is easy to extend the work to any other redshift.

We take thousands of sightlines emanating from halos with masses ∼10^10–12 M_⊙ (comparable to typical z ≳ 5 halo masses for −22 ≲ ${M}_{\mathrm{UV}}$ ≲ −18 galaxies; Barone-Nugent et al. 2014; Harikane et al. 2016, 2017) and compute the damping wing optical depth, ${\tau }_{{\rm{d}}}$, for Lyα emission as a function of velocity offset from the systemic redshift of the source halos in boxes with a range of global neutral fractions. Median values of $\exp [-{\tau }_{{\rm{d}}}]$ along ∼50 (to the rarest high mass halos) to ≳4000 (to typical 10^10.5 M_⊙ halos) sightlines are plotted in Figure 4 for a range of halo masses and ${\overline{x}}_{{\rm{H}}{\rm{I}}}$. The optical depths are smooth functions of velocity and clearly damp Lyα more effectively for higher ${\overline{x}}_{{\rm{H}}{\rm{I}}}$. In general, higher mass halos have lower optical depths to Lyα, as their large bias means they are more likely to live in the centers of large ${\rm{H}}\,{\rm{II}}$ regions, relatively more distant from the cosmic H i patches that produce the damping wing absorption during the EoR.

Zoom In Zoom Out Reset image size

For a given sightline, the final fraction of Lyα photons emitted by a galaxy in a halo with mass M_h, which are transmitted through the IGM, ${{ \mathcal T }}_{\mathrm{IGM}}$, is given by

$$\begin{eqnarray}&&{{ \mathcal T }}_{\mathrm{IGM}}({\overline{x}}_{{\rm{H}}{\rm{I}}},{M}_{h},{\rm{\Delta }}{\text{}}v)={\int }_{0}^{\infty }{dv}\,{J}_{\alpha }({\rm{\Delta }}{\text{}}v,{M}_{h},v){e}^{-{\tau }_{{\rm{i}}{gm}}({\overline{x}}_{{\rm{H}}{\rm{I}}},{M}_{h},v)},\end{eqnarray} \tag{ 6 }$$

where Δv is the velocity offset of the Lyα line center from the systemic redshift of the source galaxy (which depends on the galaxy's ISM, as described in Section 2.1) and J_α(Δv, M_h, v) is the line profile of Lyα photons escaping from the galaxy as function of velocity v.

As any photons emitted bluer than the halo circular velocity will be resonantly absorbed by intervening neutral hydrogen (Gunn & Peterson 1965; Dijkstra et al. 2007; Zheng et al. 2010; Laursen et al. 2011; Schroeder et al. 2013), J_α is described as

$$\begin{eqnarray}{J}_{\alpha }({\rm{\Delta }}{\text{}}v,{M}_{h},v)\propto \left\{\begin{array}{ll}\displaystyle \frac{1}{\sqrt{2\pi }{\sigma }_{\alpha }}{e}^{-\tfrac{{(v-{\rm{\Delta }}v)}^{2}}{2{\sigma }_{\alpha }^{2}}} & \mathrm{if}\ v\geqslant {v}_{\mathrm{circ}}({M}_{h})\\ 0 & \mathrm{otherwise}\end{array}\right..\end{eqnarray} \tag{ 7 }$$

If J_α is normalized, then ${{{ \mathcal T }}_{\mathrm{IGM}}}_{,6}=1$, as we assume ${\tau }_{{\rm{d}}}=0$ at z ∼ 6. Compared to the intrinsic emitted line, ${{{ \mathcal T }}_{\mathrm{IGM}}}_{,6}$ can be very low (Dijkstra et al. 2007; Zheng et al. 2010; Laursen et al. 2011). For ease of notation we refer to the differential transmission at z ∼ 7, ${{{ \mathcal T }}_{\mathrm{IGM}}}_{,7}/{{{ \mathcal T }}_{\mathrm{IGM}}}_{,6}$, as ${{ \mathcal T }}_{\mathrm{IGM}}$.

Example intrinsic and transmitted emission lines are plotted in Figure 1. Sightline median values of ${{ \mathcal T }}_{\mathrm{IGM}}({\overline{x}}_{{\rm{H}}{\rm{I}}},{\rm{\Delta }}{\text{}}v)$ at fixed halo mass are plotted in Figure 5. As expected, as the neutral fraction increases, the transmission fraction of Lyα decreases smoothly. While at low neutral fractions the velocity offset of Lyα has little impact, in a predominantly neutral universe, (${\overline{x}}_{{\rm{H}}{\rm{I}}}$ ≳ 0.6) lines are more easily transmitted if they were emitted at high velocity offset.

Zoom In Zoom Out Reset image size

In Figure 6 we plot probability distribution functions for ${{ \mathcal T }}_{\mathrm{IGM}}$ for three different values of ${M}_{\mathrm{UV}}$, where we have transformed from halo mass to ${M}_{\mathrm{UV}}$ using the Mason et al. (2015) LF model as above and drawn Δv values for halos using the distribution in Equation (2). The transmission distributions evolve smoothly with neutral fraction and UV magnitude. Transmission of Lyα evolves more slowly for the brightest galaxies, due to a combination of their increased velocity offsets and their locations in the most overdense regions, far from the cosmic H i regions, which cause the damping wing absorption.

Zoom In Zoom Out Reset image size

Galaxies in high mass halos (${M}_{h}\gt {10}^{12}{M}_{\odot }$, corresponding to approximately ${M}_{\mathrm{UV}}$ < −22) require special attention. First, they are rare, and lines of sight to such objects in the simulations are not well-sampled, leading to large statistical errors. Second, the correlation between ${M}_{\mathrm{UV}}$ and M_h is particularly uncertain in this regime. Third, such bright galaxies have been observed to buck the trend in the declining Lyα emission fraction at z ≳ 7 at z > 7.5 (Curtis-Lake et al. 2012; Stark et al. 2017). For these reasons, they require special attention, especially because they are prime targets for detailed spectroscopic follow-up. Since they are intrinsically rare, they would contribute negligibly to the analysis presented in this paper. Thus we leave their analysis for future work (Mason et al. 2018) and exclude them from the sample considered here.

Our observed data are a set of Lyα equivalent widths (and limits) and absolute magnitudes from galaxies at a given redshift: {W, ${M}_{\mathrm{UV}}$}. Due to the complexity of the IGM topology, there is no simple analytic model to express the likelihood of obtaining these data, given a neutral fraction ${\overline{x}}_{{\rm{H}}{\rm{I}}}$. Thus we use our model to generate large samples of mock observations that provide a non-analytic likelihood.

We take IGM simulations with global neutral fractions 0.01 ≤ ${\overline{x}}_{{\rm{H}}{\rm{I}}}$ ≤ 0.95 (Δ ${\overline{x}}_{{\rm{H}}{\rm{I}}}$ ∼ 0.02) and a population of halos with masses 10¹⁰ ≲ M_h[M_⊙] ≲ 10¹² with ${\rm{\Delta }}\mathrm{log}{M}_{h}\sim 0.1$. This mass range corresponds to UV magnitudes of −16 ≳ ${M}_{\mathrm{UV}}$ ≳ −22 at z ∼ 7 (Mason et al. 2015). The likelihood is computed in the following way:

• 1.  
  Obtain the Lyα damping wing optical depths (see Section 2.2) along thousands of different sightlines to individual halos in each simulation, to account for the inhomogeneous nature of reionization.
• 2.  
  For a grid of UV magnitudes $-22\leqslant {M}_{\mathrm{UV}}\leqslant -17$, we nearest-neighbor match the simulation halo masses with UV magnitudes at z ∼ 7, given by the relation in Mason et al. (2015), which is consistent with ${M}_{\mathrm{UV}}$–M_h measurements from clustering at z ∼ 7 (Barone-Nugent et al. 2014; Harikane et al. 2016, 2017). We do not add scatter to this matching, but note the halo mass step in the simulations (∼0.13 dex) is not dissimilar to the scatter in the inferred ${M}_{\mathrm{UV}}$–M_h relation for galaxies around ${M}_{\mathrm{UV}}^{\star }$ (e.g., 0.3 dex; Finkelstein et al. 2015b), so some ${M}_{\mathrm{UV}}$–M_h scatter is included. Furthermore we note the optical depth scatter between sightlines for a given halo mass is far greater than the scatter among sightlines between halos of different masses (compare lines and shaded region in Figure 4); thus the ${M}_{\mathrm{UV}}$–M_h scatter is sub-dominant.
• 3.  
  Populate these model galaxies with Lyα line velocity offsets from the distribution $p({\rm{\Delta }}{\text{}}v\,| \,{M}_{h})$, as described by Equation (2), including the scatter σ_v, and the Lyα equivalent widths for an ionized universe (we which assumed to be the same as at z ∼ 6), ${p}_{6}({W}_{\mathrm{em}}\,| \,{M}_{\mathrm{UV}})$ described in Section 2.1, creating 10⁶ realizations of a galaxy in each halo.
• 4.  
  We compute the differential Lyα transmission fraction, ${{ \mathcal T }}_{\mathrm{IGM}}$, with Equation (6) along sightlines through the IGM to every model galaxy and the observed equivalent width, where ${W}_{\mathrm{obs}}={{ \mathcal T }}_{\mathrm{IGM}}\times {W}_{\mathrm{em}}$.
• 5.  
  The distributions of model observed W_obs at fixed (${\overline{x}}_{{\rm{H}}{\rm{I}}}$, ${M}_{\mathrm{UV}}$) are described by the form

  $$\begin{eqnarray}p(W\,| \,{\overline{x}}_{{\rm{H}}{\rm{I}}},{M}_{\mathrm{UV}}) & = & A({M}_{\mathrm{UV}})f(W,{\overline{x}}_{{\rm{H}}{\rm{I}}})H(W)\\ & & +[1-A({M}_{\mathrm{UV}})]\delta (W),\end{eqnarray} \tag{ 13 }$$

  where f(W, ${\overline{x}}_{{\rm{H}}{\rm{I}}}$) describes the evolution of the equivalent width distribution as the neutral fraction evolves and is fitted with a Gaussian Kernel Density Estimator (Rosenblatt 1956; Parzen 1962), and A(${M}_{\mathrm{UV}}$) denotes the fraction of non-emitters and contaminants, as described in Equation (4), that does not change as the neutral fraction increases (${{ \mathcal T }}_{\mathrm{IGM}}\ne 0$ exactly).

These distributions $p(W\,| \,{\overline{x}}_{{\rm{H}}{\rm{I}}},{M}_{\mathrm{UV}})$ are the likelihoods for the observed data. Some examples are plotted in Figure 7. For increasing neutral fraction, the EW distribution becomes steeper, as more Lyα is damped by cosmic neutral patches. The evolution of $p(W\,| \,{\overline{x}}_{{\rm{H}}{\rm{I}}},{M}_{\mathrm{UV}})$ is slower for more UV bright (more massive) galaxies because the transmission functions evolve more slowly with increasing neutral fraction (see Section 2.2 and Figure 6).

We chose to marginalize out Δv at this stage to ease computation by reducing a degree of freedom, but it is possible to produce the likelihood conditional on Δv: $p(W\,| \,{\overline{x}}_{{\rm{H}}{\rm{I}}},{M}_{\mathrm{UV}},{\rm{\Delta }}{\text{}}v)$. It is then possible to infer Δv for an individual observed galaxy, or if Δv is already known, recover a narrower posterior on the neutral fraction.

To test our inference framework we perform simulated surveys of LBG follow-up. We draw a realistic sample of LBGs at z ∼ 7 from the Mason et al. (2015) UV luminosity function model (which is consistent with all observations, including new deep data from the Hubble Frontier Fields at z ≳ 7; e.g., Atek et al. 2015b; McLeod et al. 2016; Bouwens et al. 2017a; Livermore et al. 2017). We populate these galaxies with an EW given by our simulated $p(W\,| \,{\overline{x}}_{{\rm{H}}{\rm{I}}},{M}_{\mathrm{UV}})$ (see Section 3.1) for several test values of the neutral fraction.

We assume an apparent magnitude limit of ${m}_{{\rm{ab}}}=28.5$, corresponding to ${M}_{\mathrm{UV}}$ ∼ −18.5 and a 5σ flux limit of 10⁻¹⁸ erg s⁻¹ cm⁻². We draw samples of 100 and 1000 total galaxies, and perform the inference on the full samples, including upper limits.

In Figure 8 we plot the resulting posterior distributions for ${\overline{x}}_{{\rm{H}}{\rm{I}}}$. With large samples we can clearly recover the input neutral fraction well. With small samples the posterior distribution is broader as we sample less of the likelihood, but the posteriors still include the input value within 1σ.

Zoom In Zoom Out Reset image size

We use the inference framework described above to infer the neutral fraction from current observations. We take the largest published sample of LBGs at z ∼ 7 with spectroscopic follow-up to date, presented in Pentericci et al. (2014). These data comprise 68 galaxies spanning UV magnitudes −22.75 ≲ ${M}_{\mathrm{UV}}$ ≲ −17.8 and include 10 intrinsically faint objects gravitationally lensed behind the Bullet Cluster (Bradač et al. 2012) as well as observations in deep HST legacy fields (Fontana et al. 2010; Vanzella et al. 2011; Ono et al. 2012; Schenker et al. 2012).

In total, the sample comprises eight independent lines of sight with field areas ∼50–100 arcmin² each. The detections are spread over these fields. Pentericci et al. (2014) quantified the cosmic variance in this sample is very low (∼6% uncertainty in the optical depth to Lyα; see also Trenti & Stiavelli 2008). Of the 68 LBGs, 12 Lyα lines were spectroscopically confirmed. In Figure 9 we plot the UV magnitude, and in Figure 10 the EW distributions for this sample. As for the De Barros et al. (2017) z ∼ 6 sample, the EW limits are obtained by inserting simulated lines of varying flux, FWHM, and redshift into raw data, and then trying to recover them. A conservative minimum flux that could be measured over the entire wavelength range is used for the EW limit. Our framework utilizes the fact that the non-detections must be fainter than this conservative limit: fainter lines could be observed, for example, in regions free of sky lines. Future work could include the wavelength-dependent line flux sensitivity in the likelihood for non-detections (Equation (11)).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The majority of targets were z₈₅₀-dropouts selected primarily using a color criteria technique, described in detail in Grazian et al. (2012), to find targets with a high probability of having redshifts 6.5 ≲ z ≲ 7.5. The median redshift for this selection function was z = 6.9 (see Figure 1 in Grazian et al. 2012). For literature targets not directly observed by the Pentericci et al. (2014) group, but included in the sample, they only included z-dropouts with colors consistent with the color selection criteria. As noted by Pentericci et al. (2014), the probability of galaxies being outside of the observable range for their setup (z ∼ 7.3) is negligible, except for the 10 objects in the Bullet Cluster (Bradač et al. 2012), where ∼48% of objects could be above this redshift due to the broad J filter used for selection of that sample (Hall et al. 2012). To test the impact of these few objects potentially being at higher redshifts, we ran the inference excluding the Bullet Cluster and found it does not significantly impact the results.

For each galaxy in this sample, we compute the likelihoods for obtaining the observed equivalent width or upper limit using Equations (10) and (11) for every value of the neutral fraction in our simulations. We exclude the brightest objects (${M}_{\mathrm{UV}}$ < −22, 1 object) due to the insufficient sampling of very massive halos in the simulations (see Section 2.2) and the uncertainty in their intrinsic EW evolution (Stark et al. 2017), but note this does not affect the inferred neutral fraction for our sample because the UV bright objects are so rare. We use an MCMC sampler (Foreman-Mackey et al. 2013) to infer the posterior distribution of ${\overline{x}}_{{\rm{H}}{\rm{I}}}$ from these data, which is shown in Figure 11. We infer a neutral fraction of ${\overline{x}}_{{\rm{H}}{\rm{I}}}={0.59}_{-0.15}^{+0.11}$ (16–84%).

Zoom In Zoom Out Reset image size

This constraint is much tighter than previous measurements of the neutral fraction from Lyα observations (e.g., Pentericci et al. 2014; Mesinger et al. 2015) because we use the full distribution of equivalent widths, $p(W| {M}_{\mathrm{UV}})$, in our inference. Previous analyses used only the fraction of galaxies emitting Lyα with W > 25 Å, f_Lyα, to constrain the neutral fraction. In Figure 11 we also plot the posterior distribution obtained if we had used only f_Lyα, that is, the posterior is $p({\overline{x}}_{{\rm{H}}{\rm{I}}}| {f}_{\mathrm{Ly}\alpha })$, where we compare the simulation f_Lyα(${\overline{x}}_{{\rm{H}}{\rm{I}}}$) derived from Equation (13) with the fraction obtained in Pentericci et al. (2014): ${f}_{\mathrm{Ly}\alpha }={0.29}_{-0.15}^{+0.20}$ (for their faint sample, 20.25 < ${M}_{\mathrm{UV}}$ < 18.75). With just the Lyα fraction, we infer a neutral fraction of ${\overline{x}}_{{\rm{H}}{\rm{I}}}$ = 0.46 ± 0.29. Clearly, using the full distribution of EW enables us to constrain the neutral fraction much more accurately, and now that large samples of LBGs with spectroscopic follow-up are available, it should become the statistic of choice for Lyα reionization inferences.

Where does this constraint sit in our consensus picture of reionization? In Figure 12 we plot constraints derived from observations of Lyα emission from galaxies (Mesinger et al. 2015); the clustering of Lyα emitting galaxies (Ouchi et al. 2010; Sobacchi & Mesinger 2015); Lyα and Lyβ forest dark fraction (McGreer et al. 2014); and QSO ULASJ1120+0641 damping wings (Greig & Mesinger 2017b). We also plot the neutral hydrogen fraction as a function of redshift, using the Mason et al. (2015) UV luminosity function model, assuming galaxies are the source of ionizing photons and using two limiting magnitudes for the galaxy population: ${M}_{\mathrm{UV}}$ < −17 (currently detectable galaxies) and ${M}_{\mathrm{UV}}$ < −12 (ultra-faint undetected galaxies). The uncertainties in the Mason et al. (2015) reionization histories come from the range of possible reionization parameters (e.g., ionizing photon escape fraction, IGM clumping factor, number of ionizing photons per UV photon).

Zoom In Zoom Out Reset image size

Our constraint is consistent within 1σ with the other constraints at z ∼ 7, providing more strong evidence that reionization is ongoing at z ∼ 7. Our constraint lies Δ ${\overline{x}}_{{\rm{H}}{\rm{I}}}$ ∼ 0.2 higher than the constraint from the z = 7.1 QSO ULASJ1120+0641 damping wings (Greig et al. 2017), but is still consistent within the uncertainties.

JWST will be uniquely equipped to observe Lyα and rest-frame optical emission lines into Cosmic Dawn, with extremely sensitive spectrometers NIRSpec and NIRISS covering 1–5 μm in a large field of view (Gardner et al. 2006). This will enable direct measurement of the Lyα Δv and detailed studies of the ISM properties of galaxies during reionization.

Using our inferred value of ${\overline{x}}_{{\rm{H}}{\rm{I}}}={0.59}_{-0.15}^{+0.11}$ for the neutral fraction at z ∼ 7, we predict the number of Lyα emitters detectable in one NIRSpec pointing (∼9 sq arcmins) by drawing galaxies from the Mason et al. (2015) UV luminosity function and populating them with EW given by our simulated $p(W\,| \,{\overline{x}}_{{\rm{H}}{\rm{I}}},{M}_{\mathrm{UV}})$. We transform Lyα equivalent width W to flux using the relation

$$\begin{eqnarray}&&f(W,m,z)={{Wf}}_{0}{10}^{-0.4{m}_{\mathrm{UV}}}\displaystyle \frac{c}{{\lambda }_{\alpha }^{2}(1+z)}{\left(\displaystyle \frac{{\lambda }_{\mathrm{UV}}}{{\lambda }_{\alpha }}\right)}^{-\beta -2},\end{eqnarray} \tag{ 14 }$$

where ${f}_{0}=3.631\times {10}^{-20}$ erg s⁻¹ Hz⁻¹ cm⁻², ${m}_{\mathrm{UV}}$ is the apparent magnitude of the UV continuum, c is the speed of light, λ_α is the rest-frame wavelength of Lyα, ${\lambda }_{\mathrm{UV}}$ is the rest-frame wavelength of the UV continuum (usually 1500 Å), and β is the UV slope. For simplicity we assume β = −2, consistent with observations of z ∼ 7 objects (e.g., Bouwens et al. 2012), though very UV faint galaxies likely have steeper slopes due to extremely low metallicities (Vanzella et al. 2016).

We plot the predicted number counts in Figure 13, where we assume a 5σ UV continuum flux limit of ${m}_{{\rm{ab}}}\gt 29$ (${M}_{\mathrm{UV}}$ ∼ −18, corresponding to ∼1 hr integration in JWST NIRCam). We predict a 3 hr exposure in one pointing (∼9 sq arcmins) with JWST NIRSpec will detect ∼6 ± 3 z ∼ 7 Lyα lines with a 5σ flux limit of ∼3 × 10⁻¹⁸ erg s⁻¹ cm⁻² (calculated using the JWST ETC), from a total of ∼80 LBG dropout detections. We also show the forecast for a cluster lensing survey (e.g., GLASS, Treu et al. 2015; Schmidt et al. 2016), assuming a simple uniform magnification factor of μ = 2 due to gravitational lensing (i.e., p(μ) = δ(μ −2)). In this case, all fluxes are magnified by μ while the area decreases by 1/μ, and assuming the same flux limit as above, we predict ∼10 ± 2 Lyα lines from a total of ∼90 LBG detections. The NIRSpec field of view is still small compared to large-scale structure at z ∼ 7, so wide area random pointing surveys will be essential to estimate the global ${\overline{x}}_{{\rm{H}}{\rm{I}}}$.

Zoom In Zoom Out Reset image size

We simulate a 10 pointing NIRSpec survey with F070LP/G140M (R = 1000), with 3 hr exposures in each field, by again sampling the Mason et al. (2015) luminosity function in a larger area. We perform the inference on these mock JWST observations at z ∼ 7, assuming ${\overline{x}}_{{\rm{H}}{\rm{I}}}$ = 0.59. This yields ∼60 detections from ∼800 LBGs. Again, we assume a 5σ flux limit of >3 × 10⁻¹⁸ erg s⁻¹ cm⁻². The posterior distribution obtained from the JWST mock observations is shown in Figure 14, with the posterior from the current observations (Section 4.2) shown for comparison. We obtain ${\overline{x}}_{{\rm{H}}{\rm{I}}}={0.60}_{-0.06}^{+0.02}$, a ∼70% reduction in uncertainty compared to the current sample. We note this is an average forecast, and a more realistic survey forecast would require sampling the simulation directly (e.g., Mesinger & Furlanetto 2008a). We also caution our mock survey assumes 100% completeness, and maximized filling of NIRSpec slits, but, nevertheless, observations with NIRSpec will constrain the neutral fraction much more tightly than current observations.

Zoom In Zoom Out Reset image size

Robust constraints on the reionization history are challenging. While quasars provide high S/N information about individual (but rare) lines of sight, they are likely to be biased to overdense and more ionized regions (Barkana & Loeb 2004; Mesinger 2010; Decarli et al. 2017), and the number densities of bright quasars drop dramatically at z > 6 (Fan et al. 2001; Manti et al. 2016; Parsa et al. 2018). Constraining reionization with large samples of galaxies clearly avoids these problems; with the help of gravitational lensing in clusters (e.g., the Frontier Fields; Lotz et al. 2017), we know there are large populations of faint galaxies at z > 6 (Yue et al. 2014; Atek et al. 2015a; Livermore et al. 2017; Vanzella et al. 2017), and GRB host galaxy searches indicate far fainter galaxies must also exist (Kistler et al. 2009; Trenti et al. 2012).

Lyα emission from galaxies has long been recognized as a probe of reionization (Haiman & Spaans 1999; Malhotra & Rhoads 2004; Santos 2004; Verhamme et al. 2006; McQuinn et al. 2007a; Dijkstra 2014), and the framework presented in this paper provides a direct constraint on the IGM neutral fraction from observations of Lyα emission from galaxies, incorporating both realistic galaxy properties and realistic IGM topologies for the first time.

Our constraint on the neutral fraction, ${\overline{x}}_{{\rm{H}}{\rm{I}}}={0.59}_{-0.15}^{+0.11}$, is consistent with other robust probes of IGM neutrality at z ∼ 7 (Mesinger et al. 2015; Greig et al. 2017), demonstrating the power of Lyα follow-up of LBGs to constrain the neutral fraction, and providing more strong evidence that the IGM is undergoing significant reionization at z ∼ 7. Using the full distribution of observed {W, ${M}_{\mathrm{UV}}$} as inputs to our inference provides much tighter constraints than using the standard "Lyα fraction," as we demonstrated in Figure 11.

Our median value lies Δ ${\overline{x}}_{{\rm{H}}{\rm{I}}}$ ∼ 0.2 higher than that inferred by Greig et al. (2017) from the QSO ULASJ1120+0641 damping wings at z = 7.1, which was obtained using the same IGM simulations, though our posterior distribution is marginally skewed to lower values (see Figure 11). This offset is not significant given the uncertainties, and does not require us to invoke any additional evolution in galaxy properties. Within the next few years larger samples, as demonstrated in our mock survey with JWST described in Section 4.3, will greatly reduce the uncertainties in our constraints from Lyα detections and non-detections.

With large samples, it will be possible to measure the variations in ${\overline{x}}_{{\rm{H}}{\rm{I}}}$ over the sky, and cross-correlate with other constraints from quasars and eventually 21 cm observations (Lidz et al. 2009; Mesinger et al. 2016; Mirocha et al. 2016; Sobacchi et al. 2016; Vrbanec et al. 2016; Greig & Mesinger 2017a) to directly observe the inhomogeneous process of reionization. Furthermore, with tighter constraints on the timeline of reionization, it will be possible to better constrain the sources of ionizing photons, as the ionizing photon budget from galaxies depends on, for example, the minimum mass/luminosity of galaxies and the rate of ionizing photons per unit UV luminosity.

In our model, we include empirically calibrated relations for both the intrinsic dependence of Lyα EW on UV magnitude and ISM radiative transfer in galaxies of a given halo mass (UV magnitude at fixed redshift), which builds in a simple redshift evolution assuming galaxies of the same UV magnitude live in less massive halos at higher redshifts. In this framework, UV faint galaxies have intrinsically higher EW than UV bright galaxies and lower Lyα velocity offsets.

These correlations are motivated by numerous observations of Lyα emission from galaxies at a range of redshifts, including very low redshift samples where detailed spatial and spectral observations are possible (Hayes et al. 2013; Yang et al. 2016). It is likely the density and distribution of neutral gas in the ISM plays a key role in the mediation of Lyα propagation through galaxies: an ISM with high column densities of neutral hydrogen, N_hi, scatters Lyα photons more significantly, spectrally, and spatially (Verhamme et al. 2006; Zheng et al. 2010). Observations of z < 4 galaxies confirm high N_hi correlates with high Lyα velocity offset (Hashimoto et al. 2015; Henry et al. 2015; Yang et al. 2016), and more Lyα extended halos (Hayes et al. 2013; Guaita et al. 2017).

With increasing redshift, where galaxies are less massive (Lacey & Cole 1993), Lyα should escape more easily with high EW: these galaxies will contain less dust (as ALMA and Plateau de Bure Interferometer [PdBI] results are suggesting; Walter et al. 2012; Ouchi et al. 2013; Ota et al. 2014; Capak et al. 2015; Maiolino et al. 2015; Schaerer et al. 2015; Bouwens et al. 2016a; Pentericci et al. 2016) and neutral gas than at low redshifts. Additionally, the covering fraction of neutral hydrogen may evolve with galaxy mass, star formation rate, stellar populations, and/or redshift. Hard ionizing spectra from low metallicity stars (which may be significant at high redshifts; Stark et al. 2015, 2017; Mainali et al. 2017; Schmidt et al. 2017) can create more ionized holes through the ISM, reducing the covering fraction, an effect that is enhanced for low mass galaxies (Trebitsch et al. 2017). A low covering fraction would facilitate Lyα escape closer to the galaxy systemic velocity, and some observations have indicated a decreasing covering fraction with redshift (Jones et al. 2013; Leethochawalit et al. 2016).

All these factors, and the correlation of Δv with halo mass, as shown in Figure 2, suggest velocity offsets should decrease with increasing redshift for galaxies at fixed UV magnitude. These low velocity offsets are correlated with reduced scattering within the ISM and thus a higher EW of Lyα. This should increase the visibility of Lyα until the IGM starts to become neutral and these low Δv lines are easily attenuated by nearby neutral hydrogen. As was noted by Mesinger et al. (2015) and Choudhury et al. (2015), this offers a simple explanation for evolving galaxy properties which may accelerate the decline in Lyα in UV faint galaxies.

A high fraction of Lyα observed in some UV bright (${M}_{\mathrm{UV}}$ < −21.5) galaxies at z > 6 (Curtis-Lake et al. 2012; Stark et al. 2017, though cf. Treu et al. 2013 for non-detections of Lyα in slightly fainter galaxies) is surprising for several reasons: the electron scattering optical depth from the Planck Collaboration (2016) favors a significant IGM neutral fraction at these redshifts, with instantaneous reionization occurring at z = 7.8–8.8; and the observed fraction of UV faint galaxies appears to steadily decrease at the same redshifts (Pentericci et al. 2014; Schenker et al. 2014). Why can we more easily see Lyα in some UV bright galaxies into the Epoch of Reionization?

The most UV bright galaxies at high redshift probably reside in halos with mass ≳10¹¹ M_⊙, which may already have stable gaseous disks, as suggested by recent ALMA observations of two UV bright galaxies at z ∼ 6 (Smit et al. 2018) and observations of stable rotation in low mass galaxies at z ∼ 1–2 (Mason et al. 2016; Stott et al. 2016). Thus it is likely Lyα photons traveling from these galaxies will experience significant radiative transfer effects with the ISM.

The enhanced Doppler shift of the emerging Lyα photons in UV bright galaxies provides some explanation for the high fraction of Lyα observations for these populations compared to UV faint galaxies at z ∼ 7. As shown in Figures 6 and 7, we predict the transmission of UV bright galaxies evolves more slowly with the evolving IGM compared to fainter objects, making them visible far into the epoch of reionization and thus prime targets for spectroscopic confirmation. With this in mind, note their underlying EW distribution is likely much steeper and has a higher peak of non-emitters than for UV faint galaxies. When Lyα is emitted from UV bright objects, it is likely to have low EW, as the photons are dispersed spatially and spectrally.

However, this effect is also highly correlated with the large-scale environment in which these galaxies reside; assessing the relative contributions of evolving galaxy properties and environment to this apparent increase in the Lyα fraction is explored by Mason et al. (2018). The high Lyα transmission of UV bright galaxies make them ideal targets for spectroscopic follow-up to understand the star formation processes occurring in the early universe.