EoR Inference from \lya

The Universe is Reionizing at : Bayesian Inference of
the IGM Neutral Fraction Using \lya Emission from Galaxies

Charlotte A. Mason, Tommaso Treu, Mark Dijkstra, Andrei Mesinger, Michele Trenti,
Laura Pentericci, Stephane de Barros, and Eros Vanzella,
Department of Physics and Astronomy, UCLA, Los Angeles, CA, 90095-1547, USA Institute of Theoretical Astrophysics, University of Oslo, P.O. Box 1029, N-0315 Oslo, Norway Scuola Normale Superiore, Piazza dei Cavalieri 7, I-56126 Pisa, Italy School of Physics, University of Melbourne, Parkville, Victoria, Australia ARC Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D) INAF Osservatorio Astronomico di Roma, Via Frascati 33, I-00040 Monteporzio (RM), Italy Observatoire de Genève, Universitè de Genève, 51 Ch. des Maillettes, 1290 Versoix, Switzerland INAF Osservatorio Astronomico di Bologna, via Ranzani 1, 40127 Bologna, Italy cmason@astro.ucla.edu
Abstract

We present a new flexible Bayesian framework for directly inferring the fraction of neutral hydrogen in the intergalactic medium (IGM) during the Epoch of Reionization (EoR, ) from detections and non-detections of Lyman Alpha (\lya) emission from Lyman Break galaxies (LBGs). Our framework combines sophisticated reionization simulations with empirical models of the interstellar medium (ISM) radiative transfer effects on \lya. We assert that the \lya line profile emerging from the ISM has an important impact on the resulting transmission of photons through the IGM, and that these line profiles depend on galaxy properties. We model this effect by considering the peak velocity offset of \lya lines from host galaxies’ systemic redshifts, which are empirically correlated with UV luminosity and redshift (or halo mass at fixed redshift). We use our framework on the sample of LBGs presented in Pentericci et al. (2014) and infer a global neutral fraction at of , consistent with other robust probes of the EoR and confirming reionization is on-going Myr after the Big Bang. We show that using the full distribution of \lya equivalent width detections and upper limits from LBGs places tighter constraints on the evolving IGM than the standard \lya emitter fraction, and that larger samples are within reach of deep spectroscopic surveys of gravitationally lensed fields and JWST NIRSpec.

Subject headings:
dark ages, reionization, first stars – galaxies: high-redshift – galaxies: evolution – intergalactic medium

1. Introduction

In the first billion years of the universe’s history, intergalactic hydrogen atoms, formed at Recombination, were ionized (e.g., Robertson et al., 2015; Mesinger, 2016; Planck Collaboration et al., 2016). This reionization of the intergalactic medium (IGM) was driven by the first sources of light: stars, and accretion disks around black holes, in galaxies. By understanding the process and timeline of reionization we can learn about the nature of these nascent populations of galaxies.

Ground-breaking observations within the last decade have provided significant information about this Epoch of Reionization (EoR, ). With the largest near-IR instruments in space and on the ground we have now discovered large populations of galaxies at (e.g., McLure et al., 2010; Trenti et al., 2011; Bradley et al., 2012; Illingworth et al., 2013; Schenker et al., 2013b; Schmidt et al., 2014; Yue et al., 2014; Bouwens et al., 2015b; Finkelstein et al., 2015b; Calvi et al., 2016). Young stars in these galaxies are likely the primary sources of reionizing photons (e.g., Lehnert & Bremer, 2003; Bouwens et al., 2003; Yan & Windhorst, 2004; Bunker et al., 2004; Finkelstein et al., 2012a; Robertson et al., 2013; Schmidt et al., 2014), though a contribution from AGN cannot be excluded (Giallongo et al., 2015; Madau & Haardt, 2015; Onoue et al., 2017): we do not know if sufficient hard ionizing photons escape from galaxies as we do not fully understand the interactions between these early galaxies and their surrounding media.

Absorption features in quasar spectra suggest reionization was largely complete by ( Gyr after the Big Bang, e.g., Fan et al., 2006; Schroeder et al., 2013; McGreer et al., 2014; Venemans et al., 2015), whilst the electron scattering optical depth to the CMB (Planck Collaboration et al., 2015, 2016; Greig & Mesinger, 2017b) indicates significant reionization was occurring at . A robust constraint, albeit from a single sightline, on on-going reionization comes from the absorption spectrum of the quasar ULAS J1120+0641, where Greig et al. (2017) recently inferred a neutral fraction of .

To produce a timeline of reionization consistent with the evolution suggested by observations generally requires optimistic assumptions about the numbers of as-yet undetected ultra-faint galaxies (Shull et al., 2012; Robertson et al., 2013; Mason et al., 2015) – which are likely the hosts of high redshift gamma-ray bursts (Kistler et al., 2009; Trenti et al., 2012), and/or the production efficiency and escape fraction of hard ionizing photons (Bouwens et al., 2015a; Ma et al., 2015; Grazian et al., 2017; Vanzella et al., 2017). However, the timeline of reionization is not well-constrained, especially beyond where quasars become extremely rare (Fan et al., 2001; Manti et al., 2016).

Into the EoR, a powerful probe of the IGM is the Lyman alpha (\lya, 1216Å) emission line from galaxies, which is strongly attenuated by neutral hydrogen (Haiman & Spaans, 1999; Malhotra & Rhoads, 2004; Santos, 2004; Verhamme et al., 2006; McQuinn et al., 2007a; Dijkstra, 2014). Observing \lya at high redshift gives us key insights into both the IGM ionization state and galaxy properties, and, whilst quasars probably live in the densest regions of the early universe (Mesinger, 2010), observing galaxies enables us to trace reionization in cosmic volumes in a less biased way.

Dedicated spectroscopic follow-up of young star-forming galaxies at high redshift, identified as photometric dropouts (Lyman break galaxies, hereafter LBGs) combined with low redshift comparison samples (Hayes et al., 2013; Yang et al., 2016) show that the fraction of LBGs emitting \lya increases with redshift (Stark et al., 2010; Curtis-Lake et al., 2012; Hayes et al., 2011; Cassata et al., 2015), likely because the dust fraction in galaxies decreases (e.g., Finkelstein et al., 2012b; Bouwens et al., 2014) which reduces the absorption of \lya (Hayes et al., 2011). However, there is a potential smoking gun signature of reionization at : recent observations show a declining fraction of \lya emitters in the LBG population with redshift (e.g. Fontana et al., 2010; Stark et al., 2010; Caruana et al., 2012; Treu et al., 2013; Caruana et al., 2014; Faisst et al., 2014; Tilvi et al., 2014; Schenker et al., 2014; Pentericci et al., 2014), as well as an evolving \lya luminosity function (e.g., Ouchi et al., 2010; Zheng et al., 2017; Ota et al., 2017; Konno et al., 2017), suggesting an increasingly neutral, but inhomogeneous, IGM (Dijkstra et al., 2013; Mesinger et al., 2015).

Robust conversions from observations to the IGM state are challenging, however, and current constraints from \lya emission measurements show some tension. The sudden drop in \lya emission from LBGs suggests a high neutral fraction at , (Dijkstra et al., 2013; Choudhury et al., 2015; Mesinger et al., 2015), whereas measurements from clustering of \lya emitters at imply a lower neutral fraction (, Ouchi et al., 2010; Sobacchi & Mesinger, 2015). These constraints are consistent within 1 but the qualitative tension motivates a more thorough treatment of the properties of \lya emitters during reionization. Given this, and that tight constraints on the reionization history can constrain properties of the sources of reionization (e.g., the minimum mass/luminosity of galaxies, and the escape fraction of ionizing photons, Bouwens et al., 2016a; Mitra et al., 2015; Greig & Mesinger, 2017a, b), we aim to develop a robust framework for inferring the ionization state of the IGM from observations \lya from galaxies.

The conversion from the evolving transmission of \lya emission from galaxies to a constraint on the IGM ionization state is non-trivial and involves physics from pc to Gpc scales. Multiple observations (e.g., Treu et al., 2013; Pentericci et al., 2014; Becker et al., 2015) and simulations (Furlanetto et al., 2006; McQuinn et al., 2007b) suggest reionization of the IGM is likely a ‘patchy’ process, with large ionized bubbles growing faster in overdense regions filled with star-forming galaxies. An accurate model of reionization must include realistic large-scale IGM structure (Trac et al., 2008; Iliev et al., 2014; Sobacchi & Mesinger, 2014).

Irrespective of reionization, as a highly resonant line, \lya photons experience significant scattering within the interstellar medium (ISM) of their host galaxies, and absorption within the circumgalactic medium (CGM) which affects the visibility of emission (Verhamme et al., 2006, 2008; Dijkstra et al., 2007; Laursen et al., 2011). ISM effects on \lya are likely to correlate with galaxy mass and star formation rate (SFR) via dust absorption, neutral hydrogen column density and covering fraction, and outflows (Erb et al., 2014; Erb, 2015; Oyarzún et al., 2016; Hayward & Hopkins, 2017; Yang et al., 2017).

UV faint galaxies () tend to be the strongest \lya emitters at all redshifts due to lower dust masses and neutral hydrogen column densities (Yang et al., 2016, 2017). However, a small sample of UV bright galaxies at with strong Spitzer/IRAC excesses have recently been observed with \lya (Finkelstein et al., 2013; Roberts-Borsani et al., 2016; Oesch et al., 2015; Zitrin et al., 2015; Stark et al., 2017), at a redshift when the IGM is expected to be significantly neutral (Planck Collaboration et al., 2016; Greig & Mesinger, 2017b). Are these objects a new class of highly ionizing galaxies (Stark et al., 2017), emitting \lya with very high EW so some flux is still observable even after attenuation in the IGM? Do they inhabit large ionized bubbles in the IGM at high redshift? How do different halo environments and ISM properties affect the impact on reionization on galaxies?

Dijkstra et al. (2011) first considered the effects of the ISM on \lya photons during reionization, using shell models (e.g. Verhamme et al., 2006; Gronke et al., 2015a) to mimic the ISM radiative transfer, and showed ISM effects had a large impact on the transmission of \lya photons through the reionizing IGM. As described above, the \lya photons’ journey through the ISM depends on galaxy properties. However, previous constraints on the evolving transmission of \lya emission at have limited treatment of this effect: Dijkstra et al. (2011) and Mesinger et al. (2015) parametrically accounted for the ISM but assumed the LBG galaxy population is homogeneous; Jensen et al. (2013) obtained similar results combining cosmological hydro-simulations of reionization with a different sub-grid prescription for \lya radiative transfer in the ISM; simpler models do not treat the ISM but consider two bins of UV bright and faint galaxies (e.g., Treu et al., 2012).

In this paper we introduce a flexible modeling framework to enable Bayesian inference of the IGM neutral fraction from detections and non-detections of \lya from LBGs. Our framework includes realistic cosmological IGM simulations which contain the large-scale structure of the reionizing IGM. We generate 1000s of sightlines through these simulations to halos born from the same density field as the IGM, and populate these halos with simple, but realistic, ISM properties drawn from empirical models, which, for the first time in a reionization model, are linked to observable galaxy properties.

Our model asserts the impact of the ISM on the \lya line profile is the most important galaxy property to consider when trying to make accurate inferences about reionization. In our model we include this effect via the peak velocity offset of the \lya line profile from systemic (), which correlates with galaxy mass (or UV magnitude at fixed redshift), for which there are a handful of measurements at (Pentericci et al., 2016; Bradač et al., 2017; Mainali et al., 2017; Stark et al., 2017). Galaxies with high \lya velocity offsets have higher probabilities of transmitting \lya photons through the IGM. This effect is robustly accounted for in our model as a nuisance parameter in our inference.

The paper is structured as follows: in Section 2 we explain the ISM, CGM, and IGM radiative transfer modeling components of our model; in Section 3 we describe our flexible Bayesian framework for inferring the neutral fraction ; in Section 4 we give our results including key insights from the model, the inferred value of from current observations and forecasts for spectroscopic surveys with the James Webb Space Telescope (JWST); we discuss our results in Section 5 and present a summary and conclusions in Section 6.

We use the Planck Collaboration et al. (2015) cosmology where (0.69, 0.31, 0.048, 0.97, 0.81, 68 km s Mpc), and all magnitudes are given in the AB system.

2. ISM, CGM, and IGM Radiative Transfer Modeling

\lya

photons are significantly affected by the neutral hydrogen they encounter within the ISM of their source galaxies, their local CGM, and the IGM through which they travel to our telescopes. To make constraints in the Epoch of Reionization we must model \lya radiative transfer in all three media. Here we describe the combination of empirical formalisms and numerical simulations to model the effect of the ISM (Section 2.1) and the CGM and IGM (Section 2.2) on \lya transmission.

2.1. Ism \lya radiative transfer

\lya

photons are produced predominantly via recombination in H 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 \lya 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 \lya 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 \lya lineshape after transmission through the ISM as a Gaussian, centered at a velocity offset 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 below in Section 2.2 even after reionization residual neutral gas in the IGM and CGM will absorb all blue flux at .

Figure 1.— The effect of the IGM on simulated line profiles. We show two example intrinsic line profiles (black dotted lines), with peak velocity offsets of 75 and 300 km s, with flux densities normalized to that of the line at 75 km s. This is the line after transmission through the ISM. The solid black shows the lineshape in an ionized universe at where all flux bluer than the halo’s circular velocity is resonantly absorbed by neutral hydrogen in the local CGM/IGM (i.e. they experience only ). The colored lines show the emission lines after transmission through a reionizing IGM with damping wing optical depths , where the median IGM attenuation is also plotted (lighter line, corresponds to right axis). Lines emitted with high velocity offsets are less attenuated by the IGM: for the line with km s of the emitted flux is observed for (green), for the line with km s this fraction rises to . For (purple) of the total flux is transmitted from the line with km s whilst is emitted for the line at km s.

2.1.1 Modeling \lya velocity offsets

Figure 2.— \lya velocity offset as a function of UV absolute magnitude (left), halo mass (right, derived from the Mason et al. (2015) UV magnitude - halo mass relation) for a collection of data from the literature (Erb et al., 2014; Steidel et al., 2014; Willott et al., 2015; Stark et al., 2015, 2017; Inoue et al., 2016; Mainali et al., 2017; Pentericci et al., 2016; Bradač et al., 2017). The gray squares show data from a sample, whilst the colored points are at . We take the distribution as complete and intrinsic and fit a log-normal distribution to the points as shown in Equation 1. The median fit is shown as a black solid line, and the gray shaded region shows the scatter. We add a 0.2 mag uncertainty to the UV magnitude measurements to account for scatter in the UV magnitude – halo mass relation and use the propagated uncertainties in halo mass in the . The hashed region in the right panel indicates the galaxies with which are discarded from fitting due to large uncertainties in assigning their halo masses. We plot the circular velocities, of halos at (dashed orange) and (dashed blue) for comparison. The relation closely traces the circular velocities, suggesting galaxy mass is a key mediator of \lya radiative transfer.

Numerous studies of star-forming galaxies at have identified the column density of neutral hydrogen () within the ISM as a key mediator of \lya radiative transfer. \lya photons traveling through highly dense neutral ISM scatter more frequently and emerge with larger velocity offsets than galaxies with lower (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 \lya escape (Jones et al., 2013; Trainor et al., 2015; Leethochawalit et al., 2016).

Recently, a correlation has been suggested between UV magnitude and \lya velocity offset (Schenker et al., 2013a; Erb et al., 2014; Song et al., 2014; Stark et al., 2015; Mainali et al., 2017; Stark et al., 2017), again indicating galaxy mass and/or SFR strongly affects \lya 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; Mason et al., 2015; Finkelstein et al., 2015a; 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 measurements from the literature (Erb et al., 2014; Steidel et al., 2014; Bradač et al., 2017; Inoue et al., 2016; Mainali et al., 2017; Pentericci et al., 2016; Stark et al., 2015, 2017; Willott et al., 2015) in Figure 2 (left) where it is clear the high redshift galaxies have lower at given compared to the low redshift galaxies, probably because they have lower mass.

To model the effect of the ISM on \lya 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 \lya line profile. This is likely an over-simplification, e.g. ‘shell’ models take parameters to model \lya 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 and halo mass of the form , where we determine empirically from observations, as described below.

We take a sample of 158 galaxies with both UV magnitudes and \lya 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 where 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 \lya emitters in narrow-band photometry with , all these objects had rest-frame optical lines detected in MOSFIRE observations. We note the Erb et al. (2014) sample does not include faint \lya emitters (Å), which may have higher velocity offsets given observed anti-correlations between \lya EW and (Hashimoto et al., 2013; Shibuya et al., 2014; Erb et al., 2014). Whilst these samples are the largest available to measure a correlation, future rest-frame optical follow-up of large samples of galaxies with detected \lya 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 \lya 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 measurements from clustering at (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., 2015a) and use the propagated uncertainties in halo mass in the .

In the right panel of Figure 2 we plot the literature 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 , which may be significantly more starbursty than average, we discard the galaxies with from further analysis. Likewise, we exclude from this inference the galaxies at with , 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 \lya escape routes. We note gravitationally lensed objects at intermediate redshifts suggest these trends hold at low mass/luminosity (e.g., a lensed galaxy at was recently observed with a \lya 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 which increases with increasing luminosity, and approximately constant variance:

(1)

where is a linear relation corresponding to the most likely at a given halo mass:

(2)

To find the parameters and we take Equation 1 as the likelihood function and perform a Bayesian inference on the galaxies with , with uniform priors on the parameters. The inferred parameters are: , and . We show this relation on 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: , where for , and otherwise. The mean velocity offset in km s can then be approximated as:

(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 () at and 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 and the circular velocity . Investigating these trends with larger samples at low redshifts with dynamical mass measurements (e.g., Trainor et al., 2015; Erb et al., 2014) could determine to what extent \lya radiative transfer depends on the gravitational potential of the halo.

2.1.2 Modeling \lya line widths

The widths of \lya lines are also likely dominated by radiative transfer effects which both shifts and broadens the line (Verhamme et al., 2006, 2008; Gronke et al., 2016). \lya velocity dispersions are also observed to be systematically higher than in nebular emission lines which are not resonantly scattered (Trainor et al., 2015).

For simplicity we model the FWHM of the \lya lines as equal to the velocity offset of the line, which accounts for the broadening of the lines through scattering and is a good approximation for the observed correlation between \lya FWHM and velocity offset (Yang et al., 2016, 2017, Verhamme et al. submitted).

2.1.3 EW distribution in an ionized universe

The key observable of \lya emission lines at high redshift is their equivalent width (EW or ), is a measure of the brightness of the emission line relative to the UV continuum. As \lya photons from high redshift galaxies are attenuated by neutral gas in the intervening CGM and IGM we observed only a fraction, (the \lya transmission fraction) of the emitted EW, i.e. , where is the emitted equivalent width without any damping due to reionization.

In this work we consider the differential evolution of \lya equivalent widths between and , 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 and is short ( 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 (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 \lya.

Thus observed equivalent widths at are , where is the transmission fraction of \lya emission for a single object at redshift . In Section 2.2 below 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 distribution of EW as a function of galaxy properties. \lya 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 , where large samples exist, including the local ‘Green Peas’, \lya 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 \lya escape fraction (Yang et al., 2017), all indicating \lya 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 \lya EW distribution is usually parameterized as a conditional probability of (Treu et al., 2012; Dijkstra & Wyithe, 2012), though dependence on UV spectral slope has also been considered (Schenker et al., 2014).

We take the EW distribution from De Barros et al. (2017) and Pentericci et al. (2018, in preparation) from a Large Program with VLT/FORS2. This sample contains 127 objects, with UV magnitudes between , of which 63% have \lya detections. We parameterize it as an exponential distribution plus a delta function:

(4)

and account for the fraction of non-emitters, and for the anti-correlation of EW with . is the Heaviside step function and is a Dirac delta function. implicitly includes contamination by low redshift interlopers in the photometric selection (the interloper fraction is for this sample assuming all non-detections were low redshift contaminants, De Barros et al., 2017), i.e. we do not distinguish between non-emitters at 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 . Recent work by Vulcani et al. (2017) supports this assumption: they found comparably low contamination fractions at and in an evaluation of photometric selections.

To find these parameters we divided the sample into three bins: ; ; and . We used Equation 4 as a likelihood () and performed a Bayesian inference to infer and for each bin, similar to the methods of Oyarzún et al. (2017), using uniform priors with and . In the inference we fully account for the the non-detections of \lya (using as the likelihood given an EW limit , in the same way as described in more detail in Section 3 below). The uncertainties and EW limits calculated by De Barros et al. (2017) are obtained using simulations which fully account for incompleteness and wavelength sensitivity. We note that in this framework the EW limits are a conservative minimum which 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 between we use a hyperbolic tangent function to connect our inferred parameters, without extrapolating beyond the range of the data. We find and Å from fitting to the data. and vary smoothly with magnitude.

We choose this exponential parameterization of the data because is gives a good description of the data and is easy to treat analytically, and has previously been shown to be an excellent fit to \lya EW PDFs (e.g., Oyarzún et al., 2017). We do not include uncertainties in these parameters and we note the parameterization of is fairly arbitrary but does not qualitatively affect \lya 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 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 . We show both the intrinsic PDF and the distribution convolved with a 5Å typical measurement error which introduces at ‘bump’ around where the underlying distribution is a delta function. We also show histograms of the EW observations of De Barros et al. (2017) and Pentericci et al. (2018, in preparation) in two bins corresponding to UV bright and faint LBGs. As shown by e.g., Verhamme et al. (2008); Stark et al. (2010) and Oyarzún et al. (2017), \lya EW strongly depends on UV magnitude.

Figure 3.— \lya equivalent width distributions for Lyman Break galaxies given by Equation 4. The dotted lines show the true distribution. For better comparison with the data, we show the PDFs convolved with a 5Å typical measurement error on as solid lines. We plot the PDFs for two values of UV magnitude: (blue, orange) . UV faint objects tend to have higher EW and a higher duty cycle of \lya emission, whereas UV bright galaxies are less likely to emit \lya and have lower EWs. We also plot the observed EW from De Barros et al. (2017) and Pentericci et al. (2018, in preparation) in UV bright (orange) and UV faint (blue) bins. In these histograms we plot all upper limits at , though note we fully account for upper limits in fitting the EW distribution and the reionization inferences (see Equation 3).

2.2. IGM and CGM \lya Radiative Transfer

A \lya 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 \lya line at , all photons emitted blue-ward of the \lya resonance (1216Å) are absorbed by the IGM as even after reionization there is still is a fraction of neutral hydrogen within H ii regions (Gunn & Peterson, 1965). Infalling overdense gas around halos can also increase the opacity of the IGM near the \lya resonance and onto the red side of the \lya line (Santos, 2004; Dijkstra et al., 2007; Laursen et al., 2011).

For simplicity we assume all \lya 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 and 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 show example model \lya 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 .

During reionization, there is an additional opacity to \lya caused by the presence of cosmic diffuse neutral hydrogen patches which attenuate the damping wing of the \lya line cross-section (Miralda-Escude, 1998). The transmission of \lya photons through the reionizing IGM is driven by the global fraction of neutral hydrogen, .

Thus the total opacity to \lya due to neutral hydrogen within the IGM is given by:

(5)

where is the damping wing optical depth which is only present during the EoR, and is the optical depth due to resonant absorption within the CGM of galaxies (infalling gas) and any neutral hydrogen within the local H ii region of a galaxy. For simplicity, we assume at both and .

In this model, we assume the universe is fully ionized at , thus the damping wing opacity only becomes important at . This may not be exactly the case, but current constraints on at suggest the neutral fraction is low (, McGreer et al., 2014) so the reionization effect on \lya emission will be small.

To obtain the damping wing optical depths requires a model of the IGM topology during reionization. Whilst observation papers of \lya 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 \lya damping optical depths from the public Evolution of 21cm Structure (EoS) suite of reionization simulations described by Mesinger et al. (2015, 2016)111http://homepage.sns.it/mesinger/EOS.html.

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 H 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 . The simulations produce maps at different redshifts and superimpose them on the halo field to produce cubes of the 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 which dominate reionization and the ionization efficiency, , where is the fraction of ionizing photons which escape galaxies into the IGM, and 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 which host galaxies which dominate reionization (, e.g., Kimm et al., 2017), is assumed to be constant in the EoS simulations. The simulations use a free parameter which 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 Planck Collaboration et al. (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 K, with , producing IGM morphologies characterized by small HII regions. Whilst 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 ( bin, e.g., Sobacchi & Mesinger, 2015, though with large spectroscopic samples, , 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 \lya from galaxies 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 \lya Absorbers (DLAs) inside the cosmic H ii regions (Bolton & Haehnelt, 2013) which has been shown to have a minor impact on the \lya 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 , 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.

Figure 4.— Median \lya IGM damping wing optical depths due to cosmic H i patches during reionization as a function of velocity offset from the center of the source halos. We plot optical depths for 5 different mass halos (indicated by tone of the line - where darkest lines are the highest mass halos) and for 4 volume-averaged neutral fractions (indicated by color). We plot the median optical depth for each halo from the 1000s of simulated sightlines. For we plot the 1 confidence region for the optical depths from all the sightlines to the halos as a shaded area, showing the large variation across sightlines.

We take 1000s of sightlines emanating from halos with masses (comparable to typical halo masses for galaxies, Barone-Nugent et al., 2014; Harikane et al., 2016, 2017) and compute the damping wing optical depth, , for \lya 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 along (to the rarest high mass halos) to (to typical halos) sightlines are plotted in Figure 4 for a range of halo masses and . The optical depths are smooth functions of velocity and clearly damp \lya more effectively for higher . In general, higher mass halos have lower optical depths to \lya as their large bias means they are more likely to live in the centers of large H ii regions, relatively more distant from the cosmic H i patches which produce the damping wing absorption during the EoR.

For a given sightline, the final fraction of \lya photons emitted by a galaxy in a halo with mass which are transmitted through the IGM, , is given by:

(6)

where is the velocity offset of the \lya line center from the systemic redshift of the source galaxy (which depends on the galaxy’s ISM, as described in Section 2.1) and is the line profile of \lya photons escaping from the galaxy as function of velocity .

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), is described as:

(7)

If is normalized , as we assume at . Compared to the intrinsic emitted line 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 , , as .

Figure 5.— Median fraction of \lya photons transmitted through the IGM, , as a function of and computed with Equation 6 from 5000 sightlines to halos with mass , assuming . Contours show transmission fractions of 25%, 50% and 75%. In a predominantly neutral universe \lya photons have higher probability of escape through predominately ionized IGM and if emitted at high velocity offsets from their originating galaxies.

Example intrinsic and transmitted emission lines are plotted in Figure 1. Sightline median values of at fixed halo mass are plotted in Figure 5. As expected, as the neutral fraction increases the transmission fraction of \lya decreases smoothly. Whilst at low neutral fractions the velocity offset of \lya has little impact, in a predominantly neutral universe, lines are more easily transmitted if they were emitted at high velocity offset.

In Figure 6 we plot probability distribution functions for for three different values of , where we have transformed from halo mass to using the Mason et al. (2015) LF model as above and drawn values for halos using the distribution in Equation 2. The transmission distributions evolve smoothly with neutral fraction and UV magnitude. Transmission of \lya 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.

Galaxies in high mass halos (, corresponding to approximately ) require special attention. First, they are rare and lines of sights to such objects in the simulations are not well-sampled leading to large statistical errors. Second, the correlation between and is particularly uncertain in this regime. Third, such bright galaxies have been observed to buck the trend in the declining \lya emission fraction at at (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.

Figure 6.— Distributions of differential \lya transmission fractions at for simulated galaxies of different UV luminosities (UV bright = darkest lines), for a range of IGM neutral fractions . As described in Section 2.2 this is the ratio of \lya transmission at to that at where there is already significant absorption within the ionized IGM (Dijkstra et al., 2007; Zheng et al., 2010; Laursen et al., 2011). The transmission fractions evolve smoothly with the neutral fraction, though the evolution is more gradual for UV bright galaxies.

3. Bayesian Inference

Figure 7.— Simulated observed distribution of \lya equivalent widths (the likelihoods for our model) for a range of neutral fractions (colors), for faint (solid line) and bright (dashed line) UV magnitudes. The intrinsic distributions (Equation 4) are shown as black lines. The EW distribution evolves significantly for the UV faint galaxies with increasing , whilst the distribution for UV bright galaxies evolves more slowly.

Bayes’ Theorem enables us to infer the posterior distribution of model parameters of interest, given our observed data from the likelihood of obtaining the data given our model and our prior information of the model parameters. The posterior probability of is written as:

(8)

where is the likelihood function, is the prior, and is the Bayesian Evidence which normalizes the posterior.

We want to obtain the posterior distribution of the volume averaged fraction of neutral hydrogen, , a global IGM property, given the observed data: measurements of \lya equivalent widths and galaxy rest-frame UV absolute magnitudes . As described in Section 2 we model both IGM and ISM effects on \lya transmission and produce forward models of the observed \lya equivalent widths for galaxies of a given UV magnitude.

Using Bayes’ Theorem we can write the posterior probability for inferred from one observation in the absence of noise as:

(9)

where is the likelihood of observing a \lya equivalent width given our forward model of the ISM and IGM, and is the prior on the neutral fraction which we assume is uniform between 0 and 1.

Usually, the likelihood function is obtained from a model with an analytic form - e.g. a normal distribution, however, due to including simulated IGM cubes, our model is complex and does not have an analytic parameterization. We therefore generate the likelihood by sampling realizations of galaxies in our model at a given and then perform a Kernel Density Estimation (Rosenblatt, 1956; Parzen, 1962) to fit a smooth probability density function to the sampled distribution. Examples of the likelihood function are shown in Figure 7. Generation of the likelihoods is described in more detail below in Section 3.1.

In reality, our observations will always have measurement uncertainties, and some observations can only place an upper limit on a measurement, given a noise level. When we include noise, our likelihood for measuring an equivalent width with Gaussian noise level becomes:

(10)

and the likelihood for upper limits, is given by:

where is the complementary error function for .

In this work we consider samples with good redshift completeness (i.e. the probability of a \lya line falling within the observable range is close to one, see Section 4.2). Thus, this inference framework does not include information about redshift in the likelihood, this is left for future work.

We can combine the inference from a set of independent observations (i.e. individual galaxies) by simply multiplying the posteriors:

(12)

3.1. Generating the likelihood

Our observed data are a set of \lya equivalent widths (and limits) and absolute magnitudes from galaxies at a given redshift: . 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 . Thus we use our model to generate large samples of mock observations which provide a non-analytic likelihood.

We take IGM simulations with global neutral fractions () and a population of halos with masses with . This mass range corresponds to UV magnitudes of at (Mason et al., 2015). The likelihood is computed in the following way:

  1. Obtain the \lya 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 we nearest-neighbor match the simulation halo masses with UV magnitudes at given by the relation in Mason et al. (2015) which is consistent with measurements from clustering at (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 ( dex) is not dissimilar to the scatter in the inferred relation for galaxies around (e.g., dex, Finkelstein et al., 2015a), so some scatter is included. Furthermore we note the optical depth scatter between sightlines for a given halo mass is far greater than the scatter between sightlines between halos of different masses (compare lines and shaded region in Figure 4), thus the scatter is sub-dominant.

  3. Populate these model galaxies with \lya line velocity offsets from the distribution as described by Equation 2, including the scatter , and the \lya equivalent widths for an ionized universe (we which assumed to be the same as at , described in Section 2.1, creating realizations of a galaxy in each halo.

  4. We compute the differential \lya transmission fraction, with Equation 6 along sightlines through the IGM to every model galaxy and the observed equivalent width, where .

  5. The distributions of model observed at fixed are described by the form:

    (13)

    where 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 denotes the fraction of non-emitters and contaminants as described in Equation 4 which does not change as the neutral fraction increases ( exactly).

These distributions 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 \lya is damped by cosmic neutral patches. The evolution of 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 at this stage to ease computation by reducing a degree of freedom, but it is possible to produce the likelihood conditional on : . It is then possible to infer for an individual observed galaxy, or, if is already known, recover a narrower posterior on the neutral fraction.

4. Results

In this section we describe the key results and predictions from our model. In Section 4.1 we show our method can accurately recover the neutral fraction for simulated datasets. We perform inference on current data from Pentericci et al. (2014) in Section 4.2. In Section 4.3 we make predictions for future surveys with JWST.

4.1. Large samples of galaxies can accurately constrain the neutral fraction

Figure 8.— Posterior distributions for from simulated samples of \lya detections from 1000 (solid lines) and 100 (dashed lines) galaxies, for a simulation input value of (blue, orange, green - the input value is shown by the vertical dotted line). With large samples input neutral fraction is recovered well. With smaller samples, the posterior is wider, but includes the true value within 1 uncertainty.

To test our inference framework we perform simulated surveys of LBG follow-up. We draw a realistic sample of LBGs at 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 e.g., Atek et al., 2015b; McLeod et al., 2016; Livermore et al., 2017; Bouwens et al., 2017a). We populate these galaxies with an EW given by our simulated (see Section 3.1) for several test values of the neutral fraction.

We assume an apparent magnitude limit of , corresponding to and a 5 flux limit of 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 . 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.

4.2. Inference from current data

We use the inference framework described above to infer the neutral fraction from current observations. We take the largest published sample of LBGs at with spectroscopic follow-up to-date, presented in Pentericci et al. (2014). These data comprise 68 galaxies spanning UV magnitudes 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).

Figure 9.— UV magnitude distributions for the sample used for the inference. We plot the median UV magnitude of the sample as a dashed line ()
Figure 10.— EW distribution for the sample used for the inference. We show both the \lya EW measurements (filled blue) and upper limits for the non-detections (orange line).

In total, the sample comprises 8 independent lines-of-sight with field areas arcmin each. The detections are spread over these fields. Pentericci et al. (2014) quantified the cosmic variance in this sample is very low ( uncertainty in the optical depth to \lya, see also Trenti & Stiavelli, 2008). Of the 68 LBGs 12 \lya 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) 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, e.g., in regions free of sky lines. Future work could include the wavelength-dependent line flux sensitivity in the likelihood for non-detections (Equation 3).

The majority of targets were -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 . The median redshift for this selection function was (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 -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 () is negligibly, except for the 10 objects in the Bullet Cluster (Bradač et al., 2012) where 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 3 for every value of the neutral fraction in our simulations. We exclude the brightest objects (, 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 from these data, which is shown in Figure 11. We infer a neutral fraction of ().

Figure 11.— Posterior distribution for from the dataset of 68 galaxies at (including 12 with detected \lya emission) from Pentericci et al. (2014). In red we plot the posterior distribution obtained from the full sample of measurements as described in Section 3, and infer a neutral fraction of (). The dashed line shows the median value, and shaded region shows the ( and ) percentile bounds. For comparison, in blue we plot the posterior for obtained if we used only the fraction of galaxies emitting \lya with Å, . In this case we infer . Using the full distribution of EW provides much more information about the evolving IGM compared to and allows for tighter constraints on the neutral fraction.
Figure 12.— The fraction of neutral hydrogen as a function of redshift. Our new constraint is plotted as a red open star. We plot constraints derived from observations of: previous estimates from the fraction of LBGs emitting \lya (open black star, Mesinger et al., 2015); the clustering of \lya emitting galaxies (square, Ouchi et al., 2010; Sobacchi & Mesinger, 2015); \lya and Ly forest dark fraction (circle, McGreer et al., 2014); and QSO damping wings (diamond, Greig & Mesinger, 2017b; Bañados et al., 2017). We offset the constraints at (QSO ULASJ1120+0641 damping wing, Greig et al., 2017, \lya fraction and our new constraint) by for clarity. We also plot the Planck Collaboration et al. (2016) redshift range of instantaneous reionization (black hatched region). We show as shaded regions the reionization history from integrating the Mason et al. (2015) UV luminosity function to two magnitude limits of (green) and (purple) and drawing from uniform distributions for the ionizing photon escape fraction % and clumping factor , and log-normal distribution for the ionizing efficiency with mean and standard deviation dex. There are many uncertainties in obtaining the reionization history from luminosity functions so these should not be taken as real constraints on the neutral fraction, but given that galaxies fainter than likely exist (e.g., Kistler et al., 2009; Weisz & Boylan-Kolchin, 2017; Livermore et al., 2017; Bouwens et al., 2017b) our result suggests high escape fractions may not be necessary for reionization.

This constraint is much tighter than previous measurements of the neutral fraction from \lya observations (e.g., Pentericci et al., 2014; Mesinger et al., 2015) because we use the full distribution of equivalent widths, in our inference. Previous analyses used only the fraction of galaxies emitting \lya with Å, , to constrain the neutral fraction. In Figure 11 we also plot the posterior distribution obtained if we had used only , i.e. the posterior is , where we compare the simulation derived from Equation 13 with the fraction obtained in Pentericci et al. (2014): (for their faint sample, ). With just the \lya fraction we infer a neutral fraction of . Clearly, using the full distribution of EW enables us to constrain the neutral fraction much more accurately and, now large samples of LBGs with spectroscopic follow-up are available, should become the statistic of choice for \lya reionization inferences.

Where does this constraint sit in our consensus picture of reionization? In Figure 12 we plot constraints derived from observations of: \lya emission from galaxies (Mesinger et al., 2015); the clustering of \lya emitting galaxies (Ouchi et al., 2010; Sobacchi & Mesinger, 2015); \lya and Ly forest dark fraction (McGreer et al., 2014); 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: (currently detectable galaxies) and (ultra-faint undetected galaxies). The uncertainties in the Mason et al. (2015) reionization histories comes from the range of possible reionization parameters (e.g., ionizing photon escape fraction, IGM clumping factor, number of ionizing photons per UV photon).

Our constraint is consistent within 1 with the other constraints at , providing more strong evidence that reionization is on-going at . Our constraint lies higher than the constraint from the QSO ULASJ1120+0641 damping wings (Greig et al., 2017), but is still consistent within the uncertainties.

4.3. Predictions for JWST

Figure 13.— Predicted cumulative number counts of LAEs with JWST NIRSpec at (gray), and (orange) using our recovered neutral fraction . Galaxies are drawn from the Mason et al. (2015) UV luminosity function model and populated with equivalent widths via - the likelihood described in Section 3.1. The number counts obtained within the () regions on are shown as dotted orange lines. We also show the cumulative number counts for a gravitationally lensed field where we assume a uniform magnification factor of (pink line), which would reveal more emission lines. We obtain the \lya fluxes using Equation 14. The dashed black line shows the flux limit for a hour exposure at with JWST NIRSpec F070LP/G140M at m calculated with the JWST ETC (https://jwst.etc.stsci.edu)
Figure 14.— Posterior distribution of for a simulated 10 pointing JWST NIRSpec survey (orange) which is able to tightly constrain the IGM neutral fraction compared to our inference on current observations (red - same as Figure 11). Dashed lines show the median of the distributions, the shaded regions show the regions. We take a 10 pointing JWST/NIRSpec mock survey at which assumes as described in Section 4.3, and perform Bayesian inference, assuming a flux limit of erg s cm. We show the posterior distribution for inferred from current data (as described in Section 4.2) for comparison. In this example JWST could reduce the uncertainty on the neutral fraction by .

JWST will be uniquely equipped to observe \lya and rest-frame optical emission lines into Cosmic Dawn, with extremely sensitive spectrometers NIRSpec and NIRISS covering m in a large field of view (Gardner et al., 2006; Stiavelli et al., 2007). This will enable direct measurement of the \lya  and detailed studies of the ISM properties of galaxies during Reionization.

Using our inferred value of for the neutral fraction at we predict the number of \lya emitters detectable in one NIRSpec pointing ( sq arcmins) by drawing galaxies from the Mason et al. (2015) UV luminosity function and populating them with EW given by our simulated . We transform \lya equivalent width to flux using the relation:

(14)

where erg s Hz cm, is the apparent magnitude of the UV continuum, is the speed of light, is the rest-frame wavelength of \lya, is the rest-frame wavelength of the UV continuum (usually 1500Å), and is the UV slope. For simplicity we assume , consistent with observations of 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 (, corresponding to hour integration in JWST NIRCam). We predict a 3 hour exposure in one pointing ( sq arcmins) with JWST NIRSpec will detect \lya lines with a flux limit of erg s cm (calculated using the JWST ETC), from a total of 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 due to gravitational lensing (i.e. ). In this case, all fluxes are magnified by whilst the area decreases by , and assuming the same flux limit as above we predict \lya lines from a total of LBG detections. The NIRSpec field-of-view is still small compared to large scale structure at so wide area random pointing surveys will be essential to estimate the global .

We simulate a 10 pointing NIRSpec survey with F070LP/G140M (), with 3 hour 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 , assuming . This yields detections from LBGs. Again, we assume a flux limit of 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 , a 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.

5. Discussion

In this section we discuss our result in the context of other probes of reionization (Section 5.1), and we discuss the implications of the mass-dependent \lya velocity offset on the evolving \lya fraction for average (Section 5.2) and UV bright (Section 5.3) galaxies.

5.1. The global reionization history

Robust constraints on the reionization history are challenging. Whilst 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 (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 (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).

\lya

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 provide a direct constraint on the IGM neutral fraction from observations of \lya emission from galaxies, incorporating both realistic galaxy properties and realistic IGM topologies for the first time.

Our constraint on the neutral fraction, , is consistent with other robust probes of IGM neutrality at (Mesinger et al., 2015; Greig et al., 2017) demonstrating the power of \lya follow-up of LBGs to constrain the neutral fraction, and providing more strong evidence that the IGM is undergoing significant reionization at . Using the full distribution of observed as inputs to our inference provides much tighter constraints than using the standard ‘\lya fraction’, as we demonstrated in Figure 11.

Our median value lies higher than that inferred by Greig et al. (2017) from the QSO ULASJ1120+0641 damping wings at , 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 \lya detections and non-detections.

With large samples, it will be possible to measure the variations in over the sky, and cross-correlate with other constraints from quasars and eventually 21cm observations (Lidz et al., 2009; Vrbanec et al., 2016; Sobacchi et al., 2016; Mirocha et al., 2016; Mesinger 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 e.g., the minimum mass/luminosity of galaxies and the rate of ionizing photons per unit UV luminosity.

5.2. A sudden drop in \lya emission
– redshift evolution of ?

In our model, we include empirically calibrated relations for both the intrinsic dependence of \lya 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 high EW than UV bright galaxies and lower \lya velocity offsets.

These correlations are motivated by numerous observations of \lya 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 \lya propagation through galaxies: an ISM with high column densities of neutral hydrogen, , scatters \lya photons more significantly, spectrally and spatially (Verhamme et al., 2006; Zheng et al., 2010). Observations of galaxies confirm high correlates with high \lya velocity offset (Yang et al., 2016; Hashimoto et al., 2015; Henry et al., 2015), and more \lya extended halos (Hayes et al., 2013; Guaita et al., 2017).

With increasing redshift, when galaxies were less massive (Lacey & Cole, 1993), \lya 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; Schaerer et al., 2015; Maiolino et al., 2015; Capak et al., 2015; Bouwens et al., 2016b; 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, Mainali et al., 2017; Stark et al., 2015; Schmidt et al., 2017; Stark et al., 2017) can create more ionized holes through the ISM, reducing the covering fraction, an effect which is enhanced for low mass galaxies (Trebitsch et al., 2017). A low covering fraction would facilitate \lya escape closer to the galaxy systemic velocity, and some observations have indicated a decreasing covering fraction with redshift (Leethochawalit et al., 2016; Jones et al., 2013).

All these factors, and the correlation of 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 \lya. This should increase the visibility of \lya until the IGM starts to become neutral and these low 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 \lya in UV faint galaxies.

5.3. \lya from UV bright galaxies
– redshifted away from resonance?

A high fraction of \lya observed in some UV bright () galaxies at (Curtis-Lake et al., 2012; Stark et al., 2017, though c.f. Treu et al. (2013) for non-detections of \lya in slightly fainter galaxies) is surprising for several reasons: the electron scattering optical depth from the Planck Collaboration et al. (2016) favors a significant IGM neutral fraction at these redshifts, with instantaneous reionization occurring at ; 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 \lya in some UV bright galaxies into the Epoch of Reionization?

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

The enhanced Doppler shift of the emerging \lya photons in UV bright galaxies provides some explanation for the high fraction of \lya observations for these populations compared to UV faint galaxies at . As shown in Figures 6 and 7, we predict 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. Though note their underlying EW distribution is likely much steeper and has a higher peak of non-emitters than for UV faint galaxies. When \lya is emitted from UV bright objects it is likely to have low EW as the photons are so dispersed spatially and spectrally.

However, this effect is also highly correlated to 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 \lya fraction is explored by Mason et al. (2018). The high \lya transmission of UV bright galaxies make them ideal targets for spectroscopic follow-up to understand the star formation processes occurring in the early universe.

6. Summary and Conclusions

We have developed a flexible Bayesian inference framework to infer the IGM neutral fraction during reionization by forward-modeling the observed equivalent width distribution of \lya emission from LBGs. Our model incorporates sightlines through realistic IGM simulations to model galaxies with realistic ISM properties.

Our main conclusions are as follows:

  1. The \lya line profile emerging from the ISM has a huge impact on the probability of transmission through the IGM (Dijkstra et al., 2011), and is related to the properties of the source galaxy. This must be systematically accounted for in reionization inference.

  2. We introduce a simple empirical relation between the halo mass of a galaxy (or UV luminosity at fixed redshift) and its \lya line peak velocity offset, where the most massive galaxies have the largest velocity offsets likely due to increased in the ISM, higher halo circular velocities and/or the presence of star-formation induced outflows.

  3. This relation predicts that with increasing redshift, \lya velocity offsets will decrease for galaxies at fixed UV luminosity, making \lya lines more susceptible to absorption in the IGM. This effect would accelerate the decline in \lya emission compared to other reionization probes and be a factor in explaining the sudden drop of \lya emission observed at .

  4. We conduct a Bayesian inference from current observations at from Pentericci et al. (2014) and infer the first direct constraint on the neutral fraction from \lya transmission of , which is consistent with other robust probes of the neutral fraction and confirms reionization is on-going at .

  5. Using the full distribution of \lya equivalent width measurements enables us to provide much tighter constraints on the neutral fraction compared to the standard ‘\lya fraction’, , used in previous analyses.

  6. We make predictions for spectroscopic surveys with JWST and find a hour LBG follow-up survey with JWST/NIRSpec could reduce the uncertainty in by .

Future near-IR spectrographs in space, such as JWST NIRSpec and NIRISS, will be able to observe both \lya and rest-frame optical lines for galaxies to and to measure SFRs and \lya velocity offsets for these objects, enabling us to further understand the interactions between star-forming regions, the ISM, and the reionizing IGM. It will soon be possible to apply our framework to large samples, free of cosmic variance, to get accurate universal constraints on the evolution of the neutral fraction.

The authors thank Dawn Erb and Dan Stark for providing their observational data. We thank Simon Birrer, Fred Davies, Max Gronke, Joe Hennawi and Crystal Martin for useful discussions. C.M. acknowledges support by NASA Headquarters through the NASA Earth and Space Science Fellowship Program Grant NNX16AO85H. A.M. acknowledges support from the European Research Council (ERC) under the European Unions Horizon 2020 research and innovation program (grant agreement No 638809 AIDA). This research was partially supported by the Australian Research Council through awards FT130101593 and CE170100013. This work was supported by the HST BoRG grants GO-12572, 12905, and 13767, and the HST GLASS grant GO-13459 This work made use of the following open source software: IPython (Pérez & Granger, 2007), matplotlib (Hunter, 2007), NumPy (Van Der Walt et al., 2011), SciPy (Oliphant, 2007), Astropy (Robitaille et al., 2013) and EMCEE (Foreman-Mackey et al., 2013).

References

  • Atek et al. (2015a) Atek, H., Richard, J., Jauzac, M., et al. 2015a, ApJ, 814, 69
  • Atek et al. (2015b) Atek, H., Richard, J., Kneib, J.-P., et al. 2015b, ApJ, 800, 18
  • Bañados et al. (2017) Bañados, E., Venemans, B. P., Mazzucchelli, C., et al. 2017, Nature
  • Barkana & Loeb (2004) Barkana, R., & Loeb, A. 2004, ApJ, 601, 64
  • Barone-Nugent et al. (2014) Barone-Nugent, R. L., Trenti, M., Wyithe, J. S. B., et al. 2014, ApJ, 793, 17
  • Becker et al. (2015) Becker, G. D., Bolton, J. S., Madau, P., et al. 2015, MNRAS, 447, 3402
  • Behroozi et al. (2013a) Behroozi, P. S., Marchesini, D., Wechsler, R. H., et al. 2013a, ApJ, 777, L10
  • Behroozi et al. (2013b) Behroozi, P. S., Wechsler, R. H., & Conroy, C. 2013b, ApJ, 770, 57
  • Bolton & Haehnelt (2013) Bolton, J. S., & Haehnelt, M. G. 2013, MNRAS, 429, 1695
  • Bouwens et al. (2017a) Bouwens, R. J., Illingworth, G. D., Oesch, P. A., et al. 2017a, ApJ, 843, 41
  • Bouwens et al. (2015a) Bouwens, R. J., Illingworth, G. D., Oesch, P. a., et al. 2015a, ApJ, 811, 140
  • Bouwens et al. (2017b) Bouwens, R. J., Oesch, P. A., Illingworth, G. D., Ellis, R. S., & Stefanon, M. 2017b, ApJ, 843, 129
  • Bouwens et al. (2016a) Bouwens, R. J., Smit, R., Labbé, I., et al. 2016a, ApJ, 831, 176
  • Bouwens et al. (2003) Bouwens, R. J., Illingworth, G. D., Rosati, P., et al. 2003, ApJ, 595, 589
  • Bouwens et al. (2012) Bouwens, R. J., Illingworth, G. D., Oesch, P., et al. 2012, ApJ, 754, 83
  • Bouwens et al. (2014) Bouwens, R. J., Illingworth, G. D., Oesch, P. A., et al. 2014, ApJ, 793, 115
  • Bouwens et al. (2015b) Bouwens, R. J., Illingworth, G. D., Oesch, P. a., et al. 2015b, ApJ, 803, 34
  • Bouwens et al. (2016b) Bouwens, R. J., Aravena, M., Decarli, R., et al. 2016b, ApJ, 833, 72
  • Bradač et al. (2012) Bradač, M., Vanzella, E., Hall, N., et al. 2012, ApJ, 755, L7
  • Bradač et al. (2017) Bradač, M., Garcia-Appadoo, D., Huang, K.-H., et al. 2017, ApJ, 836, L2
  • Bradley et al. (2012) Bradley, L., Trenti, M., Oesch, P. A., et al. 2012, ApJ, 760, 108
  • Bunker et al. (2004) Bunker, A. J., Stanway, E. R., Ellis, R. S., & McMahon, R. G. 2004, MNRAS, 355, 374
  • Calvi et al. (2016) Calvi, V., Trenti, M., Stiavelli, M., et al. 2016, ApJ, 817, 120
  • Capak et al. (2015) Capak, P. L., Carilli, C., Jones, G., et al. 2015, Nature, 522, 455
  • Caruana et al. (2012) Caruana, J., Bunker, A. J., Wilkins, S. M., et al. 2012, MNRAS, 427, 3055
  • Caruana et al. (2014) —. 2014, MNRAS, 443, 2831
  • Cassata et al. (2015) Cassata, P., Tasca, L. A. M., Le Fèvre, O., et al. 2015, A&A, 573, A24
  • Chonis et al. (2013) Chonis, T. S., Blanc, G. A., Hill, G. J., et al. 2013, ApJ, 775, 99
  • Choudhury et al. (2015) Choudhury, T. R., Puchwein, E., Haehnelt, M. G., & Bolton, J. S. 2015, MNRAS, 452, 261
  • Curtis-Lake et al. (2012) Curtis-Lake, E., McLure, R. J., Pearce, H. J., et al. 2012, MNRAS, 422, 1425
  • De Barros et al. (2017) De Barros, S., Pentericci, L., Vanzella, E., et al. 2017, A&A, 608, A123
  • Decarli et al. (2017) Decarli, R., Walter, F., Venemans, B. P., et al. 2017, Nature, 545, 457
  • Dijkstra (2014) Dijkstra, M. 2014, PASA, 31, e040
  • Dijkstra et al. (2007) Dijkstra, M., Lidz, A., & Wyithe, J. S. B. 2007, MNRAS, 377, 1175
  • Dijkstra et al. (2011) Dijkstra, M., Mesinger, A., & Wyithe, J. S. B. 2011, MNRAS, 414, 2139
  • Dijkstra & Wyithe (2012) Dijkstra, M., & Wyithe, J. S. B. 2012, MNRAS, 419, 3181
  • Dijkstra et al. (2013) Dijkstra, M., Wyithe, S., Haiman, Z., Mesinger, A., & Pentericci, L. 2013, MNRAS, 440, 3309
  • Erb (2015) Erb, D. K. 2015, Nature, 523, 169
  • Erb et al. (2014) Erb, D. K., Steidel, C. C., Trainor, R. F., et al. 2014, ApJ, 795, 33
  • Faisst et al. (2014) Faisst, A. L., Capak, P., Carollo, C. M., Scarlata, C., & Scoville, N. 2014, ApJ, 788, 87
  • Fan et al. (2001) Fan, X., Narayanan, V. K., Lupton, R. H., et al. 2001, Astron. J., 122, 2833
  • Fan et al. (2006) Fan, X., Strauss, M. A., Becker, R. H., et al. 2006, Astron. J., 132, 117
  • Finkelstein et al. (2012a) Finkelstein, S. L., Papovich, C., Ryan, R. E., et al. 2012a, ApJ, 758, 93
  • Finkelstein et al. (2012b) Finkelstein, S. L., Papovich, C., Salmon, B., et al. 2012b, ApJ, 756, 164
  • Finkelstein et al. (2013) Finkelstein, S. L., Papovich, C., Dickinson, M., et al. 2013, Nature, 502, 524
  • Finkelstein et al. (2015a) Finkelstein, S. L., Song, M., Behroozi, P., et al. 2015a, ApJ, 814, 95
  • Finkelstein et al. (2015b) Finkelstein, S. L., Ryan, R. E., Papovich, C., et al. 2015b, ApJ, 810, 71
  • Fontana et al. (2010) Fontana, A., Vanzella, E., Pentericci, L., et al. 2010, ApJ, 725, L205
  • Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306
  • Furlanetto et al. (2006) Furlanetto, S. R., Zaldarriaga, M., & Hernquist, L. 2006, MNRAS, 365, 1012
  • Gardner et al. (2006) Gardner, J. P., Mather, J. C., Clampin, M., et al. 2006, Space Sci. Rev., 123, 485
  • Giallongo et al. (2015) Giallongo, E., Grazian, A., Fiore, F., et al. 2015, A&A, 578, A83
  • Grazian et al. (2012) Grazian, A., Castellano, M., Fontana, A., et al. 2012, A&A, 547, A51
  • Grazian et al. (2017) Grazian, A., Giallongo, E., Paris, D., et al. 2017, A&A, 602, A18
  • Greig & Mesinger (2017a) Greig, B., & Mesinger, A. 2017a, MNRAS, 472, 2651
  • Greig & Mesinger (2017b) —. 2017b, MNRAS, 465, 4838
  • Greig et al. (2017) Greig, B., Mesinger, A., Haiman, Z., & Simcoe, R. A. 2017, MNRAS, 466, 4239
  • Gronke et al. (2015a) Gronke, M., Bull, P., & Dijkstra, M. 2015a, ApJ, 812, 123
  • Gronke & Dijkstra (2016) Gronke, M., & Dijkstra, M. 2016, ApJ, 826, 14
  • Gronke et al. (2016) Gronke, M., Dijkstra, M., McCourt, M., & Peng Oh, S. 2016, ApJ, 833, L26
  • Gronke et al. (2015b) Gronke, M., Dijkstra, M., Trenti, M., & Wyithe, S. 2015b, MNRAS, 449, 1284
  • Guaita et al. (2017) Guaita, L., Talia, M., Pentericci, L., et al. 2017, A&A, 606, A19
  • Gunn & Peterson (1965) Gunn, J. E., & Peterson, B. A. 1965, ApJ, 142, 1633
  • Haiman & Spaans (1999) Haiman, Z., & Spaans, M. 1999, ApJ, 518, 138
  • Hall et al. (2012) Hall, N., Bradač, M., Gonzalez, A. H., et al. 2012, ApJ, 745, 155
  • Harikane et al. (2016) Harikane, Y., Ouchi, M., Ono, Y., et al. 2016, ApJ, 821, 123
  • Harikane et al. (2017) —. 2017, PASJ, 00, 1
  • Hashimoto et al. (2013) Hashimoto, T., Ouchi, M., Shimasaku, K., et al. 2013, ApJ, 765, 70
  • Hashimoto et al. (2015) Hashimoto, T., Verhamme, A., Ouchi, M., et al. 2015, ApJ, 812, 157
  • Hayes et al. (2011) Hayes, M., Schaerer, D., Östlin, G., et al. 2011, ApJ, 730, 8
  • Hayes et al. (2013) Hayes, M., Östlin, G., Schaerer, D., et al. 2013, ApJ, 765, L27
  • Hayward & Hopkins (2017) Hayward, C. C., & Hopkins, P. F. 2017, MNRAS, 465, 1682
  • Henry et al. (2015) Henry, A., Scarlata, C., Martin, C. L., & Erb, D. 2015, ApJ, 809, 19
  • Herenz et al. (2017) Herenz, E. C., Urrutia, T., Wisotzki, L., et al. 2017, A&A, 606, A12
  • Hunter (2007) Hunter, J. D. 2007, Comput. Sci. Eng., 9, 99
  • Iliev et al. (2014) Iliev, I. T., Mellema, G., Ahn, K., et al. 2014, MNRAS, 439, 725
  • Illingworth et al. (2013) Illingworth, G. D., Magee, D., Oesch, P. A., et al. 2013, ApJS, 209, 6
  • Inoue et al. (2016) Inoue, A. K., Tamura, Y., Matsuo, H., et al. 2016, Science, 352, 1559
  • Jensen et al. (2013) Jensen, H., Laursen, P., Mellema, G., et al. 2013, MNRAS, 428, 1366
  • Jones et al. (2013) Jones, T. A., Ellis, R. S., Schenker, M. A., & Stark, D. P. 2013, ApJ, 779, 52
  • Kakiichi & Dijkstra (2017) Kakiichi, K., & Dijkstra, M. 2017, ArXiv e-prints, arXiv:1710.10053
  • Kimm et al. (2017) Kimm, T., Katz, H., Haehnelt, M., et al. 2017, MNRAS, 466, 4826
  • Kistler et al. (2009) Kistler, M. D., Yüksel, H., Beacom, J. F., Hopkins, A. M., & Wyithe, J. S. B. 2009, ApJ, 705, L104
  • Konno et al. (2017) Konno, A., Ouchi, M., Shibuya, T., et al. 2017, PASJ, 00, 1
  • Lacey & Cole (1993) Lacey, C., & Cole, S. 1993, MNRAS, 262, 627
  • Laursen et al. (2011) Laursen, P., Sommer-Larsen, J., & Razoumov, A. O. 2011, ApJ, 728, 52
  • Leethochawalit et al. (2016) Leethochawalit, N., Jones, T. A., Ellis, R. S., Stark, D. P., & Zitrin, A. 2016, ApJ, 831, 13
  • Lehnert & Bremer (2003) Lehnert, M. D., & Bremer, M. 2003, ApJ, 593, 630
  • Lidz et al. (2009) Lidz, A., Zahn, O., Furlanetto, S. R., et al. 2009, ApJ, 690, 252
  • Livermore et al. (2017) Livermore, R. C., Finkelstein, S. L., & Lotz, J. M. 2017, ApJ, 835, 113
  • Lotz et al. (2017) Lotz, J. M., Koekemoer, A., Coe, D., et al. 2017, ApJ, 837, 97
  • Ma et al. (2015) Ma, X., Kasen, D., Hopkins, P. F., et al. 2015, MNRAS, 453, 960
  • Madau & Haardt (2015) Madau, P., & Haardt, F. 2015, ApJ, 813, L8
  • Mainali et al. (2017) Mainali, R., Kollmeier, J. A., Stark, D. P., et al. 2017, ApJ, 836, L14
  • Maiolino et al. (2015) Maiolino, R., Carniani, S., Fontana, A., et al. 2015, MNRAS, 452, 54
  • Malhotra & Rhoads (2004) Malhotra, S., & Rhoads, J. E. 2004, ApJ, 617, L5
  • Manti et al. (2016) Manti, S., Gallerani, S., Ferrara, A., Greig, B., & Feruglio, C. 2016, MNRAS, 466, 1160
  • Mason et al. (2015) Mason, C. A., Trenti, M., & Treu, T. 2015, ApJ, 813, 21
  • Mason et al. (2016) Mason, C. A., Treu, T., Fontana, A., et al. 2016, ApJ, 838, 14
  • Mason et al. (2018) Mason, C. A., Treu, T., de Barros, S., et al. 2018, Submitted to ApJL, arXiv:1801.01891
  • McGreer et al. (2014) McGreer, I. D., Mesinger, A., & D’Odorico, V. 2014, MNRAS, 447, 499
  • McLean et al. (2012) McLean, I. S., Steidel, C. C., Epps, H. W., et al. 2012, in Proc. SPIE, ed. I. S. McLean, S. K. Ramsay, & H. Takami, Vol. 8446, 84460J
  • McLeod et al. (2016) McLeod, D. J., McLure, R. J., & Dunlop, J. S. 2016, MNRAS, 459, 3812
  • McLinden et al. (2011) McLinden, E. M., Finkelstein, S. L., Rhoads, J. E., et al. 2011, ApJ, 730, 136
  • McLure et al. (2010) McLure, R. J., Dunlop, J. S., Cirasuolo, M., et al. 2010, MNRAS, 403, 960
  • McQuinn et al. (2007a) McQuinn, M., Hernquist, L., Zaldarriaga, M., & Dutta, S. 2007a, MNRAS, 381, 75
  • McQuinn et al. (2007b) McQuinn, M., Lidz, A., Zahn, O., et al. 2007b, MNRAS, 377, 1043
  • Mesinger (2010) Mesinger, A. 2010, MNRAS, 407, 1328
  • Mesinger (2016) Mesinger, A., ed. 2016, Astrophysics and Space Science Library, Vol. 423, Understanding the Epoch of Cosmic Reionization (Cham: Springer International Publishing)
  • Mesinger et al. (2015) Mesinger, A., Aykutalp, A., Vanzella, E., et al. 2015, MNRAS, 446, 566
  • Mesinger & Furlanetto (2008a) Mesinger, A., & Furlanetto, S. R. 2008a, MNRAS, 385, 1348
  • Mesinger & Furlanetto (2008b) —. 2008b, MNRAS, 386, 1990
  • Mesinger et al. (2016) Mesinger, A., Greig, B., & Sobacchi, E. 2016, MNRAS, 459, 2342
  • Miralda-Escude (1998) Miralda-Escude, J. 1998, ApJ, 501, 15
  • Mirocha et al. (2016) Mirocha, J., Furlanetto, S. R., & Sun, G. 2016, MNRAS, 6, 1365
  • Mitra et al. (2015) Mitra, S., Roy Choudhury, T., & Ferrara, A. 2015, MNRAS, 454, L76
  • Mostardi et al. (2013) Mostardi, R. E., Shapley, A. E., Nestor, D. B., et al. 2013, ApJ, 779, 65
  • Moster et al. (2017) Moster, B. P., Naab, T., & White, S. D. M. 2017, ArXiv e-prints, arXiv:1705.05373
  • Oesch et al. (2015) Oesch, P. A., Dokkum, P. G. V., Illingworth, G. D., et al. 2015, ApJ, 804, L30
  • Oliphant (2007) Oliphant, T. E. 2007, Comput. Sci. Eng., 9, 10
  • Ono et al. (2012)