Conditions for Reionizing the Universe with a Low Galaxy Ionizing Photon Escape Fraction

To understand the contribution from galaxies to the ionizing emissivity, we adopt the "reference" luminosity functions of Finkelstein (2016, hereafter F16), which were the result of an MCMC Schechter function fit to all recently published data at z = 4–10. Rather than fitting luminosity functions separately at each redshift, they fit all data simultaneously, solving for the linear relations of the characteristic magnitude M*(z), faint-end slope α(z), and characteristic number density ϕ*(z). To incorporate the uncertainties in these fits into our analysis, we make use of the MCMC chains from F16, using a randomly chosen sample of 10⁴ chain steps (which were verified to be representative of the full chain). We note that while the F16 analysis did not include results from lensed galaxies in the Hubble Frontier Fields, the faint-end slopes used here are consistent with studies of those data, which reach to M_UV ≳ −15 at z = 7, finding α ≈ −2 (e.g., Bouwens et al. 2017; Livermore et al. 2017; Atek et al. 2018).

To fully explore the epoch of reionization, it is necessary to extrapolate these results to higher redshift. In our analysis, we consider redshifts from z = 4 to 20. As the data used to derive these luminosity functions were limited to z ≤ 10, it is unknown if this extrapolation is valid. Specifically, the faint-end slope α is observed to evolve somewhat steeply with redshift (dα/dz ∝ −0.11), implying α(z = 15) = −2.90. While this may be possible (indeed, the model results from Mason et al. 2015 find α[z = 16] = −3.51), it is unknown if this is actually the case. To cover this possibility, in our analysis, we consider two scenarios: our fiducial model is one in which the faint-end slope ceases to evolve at z > 10 and remains at the z = 10 value of −2.35, while in Section 4.2 we also explore the case where α continues to steepen at z > 10. The fiducial luminosity function parameters at the redshifts considered here are given in Table 1, and they are shown in the left panel of Figure 1.

Zoom In Zoom Out Reset image size

Table 1.  UV Luminosity Function Parameters

Redshift	M*	α	log ϕ*
	(mag)		(Mpc−3)
4	−21.05${}_{-0.06}^{+0.05}$	−1.69${}_{-0.04}^{+0.03}$	−2.99${}_{-0.04}^{+0.04}$
6	−20.79${}_{-0.04}^{+0.05}$	−1.91${}_{-0.03}^{+0.04}$	−3.37${}_{-0.04}^{+0.05}$
8	−20.52${}_{-0.04}^{+0.06}$	−2.13${}_{-0.03}^{+0.05}$	−3.75${}_{-0.04}^{+0.06}$
10	−20.25${}_{-0.06}^{+0.07}$	−2.35${}_{-0.04}^{+0.06}$	−4.13${}_{-0.06}^{+0.08}$
12	−19.98${}_{-0.08}^{+0.09}$	−2.57${}_{-0.06}^{+0.08}$	−4.50${}_{-0.07}^{+0.10}$
14	−19.71${}_{-0.10}^{+0.11}$	−2.79${}_{-0.07}^{+0.10}$	−4.88${}_{-0.09}^{+0.12}$
16	−19.44${}_{-0.13}^{+0.14}$	−3.01${}_{-0.09}^{+0.11}$	−5.25${}_{-0.11}^{+0.14}$
18	−19.17${}_{-0.15}^{+0.16}$	−3.23${}_{-0.11}^{+0.13}$	−5.63${}_{-0.12}^{+0.16}$
20	−18.90${}_{-0.18}^{+0.19}$	−3.45${}_{-0.12}^{+0.15}$	−6.00${}_{-0.14}^{+0.18}$

Note. The assumed rest-frame UV luminosity functions used in this work, following the evolutionary trend derived via observations as discussed in F16. Our fiducial model keeps the faint-end slope α fixed at z > 10 to the z = 10 value of −2.35.

Download table as:  ASCIITypeset image

As discussed in the following two subsections, we must project the observed UV luminosity function onto the underlying dark matter halo mass function. We use this mapping to obtain both the limiting magnitude and the escape fraction for a given UV luminosity. We follow the abundance-matching methods of Behroozi et al. (2013) to map the halo mass to each point on our luminosity functions. Following Finkelstein et al. (2015b), we assume a lognormal UV magnitude scatter at a fixed halo mass of 0.2 dex, though we note that a scatter as high as 0.4 dex does not affect the M_halo–M_UV relation at log(M_h/M_⊙) < 11. Our derived M_halo–M_UV relations are shown in the middle panel of Figure 1 for both luminosity function evolution cases we consider.

A complete accounting of the available photon budget requires us to include star formation in all galaxies, including those that are too faint to be observed directly. A recent analysis indicates that current observations using lensing at z = 6 probe galaxies hosted by log(M/M_⊙) = 9.5 halos (Finlator et al. 2016), and theoretical models generally predict that star formation in even lower-mass systems is expected (e.g., Paardekooper et al. 2013; Xu et al. 2016). We thus must extrapolate beyond what is observed, yet as the z ≥ 6 luminosity function has a very steep faint-end slope (e.g., Bouwens et al. 2015b; Finkelstein et al. 2015a; Livermore et al. 2017), small changes in the minimum luminosity can have a large impact on the total luminosity density. Star formation should occur in any halo that can both retain its gas and cool it to temperatures where it can condense and form stars. Efficient cooling via collisional excitation of H i can occur in galaxies with halo virial temperatures below ∼10⁴ K. This corresponds to log(M_h/M_⊙) ≈ 8 at z = 6 (Okamoto et al. 2008; Finlator et al. 2012), and this critical mass shifts to lower values at higher redshifts, as halos that collapse at earlier times have steeper dark matter potential wells and thus correspondingly higher virial velocities (e.g., Barkana & Loeb 2001). Lower-mass halos can efficiently cool if they have metals, as predicted by recent simulations that find that significant star formation is happening down to log(M_h/M_⊙) ≈ 7 at z ∼ 10–15 due to the availability of metal-line cooling in the immediate aftermath of the formation of Population III stars (e.g., Wise et al. 2014; Xu et al. 2016). Lacking metals, gas can cool inefficiently via molecular hydrogen cooling, which is believed to be the dominant cooling mechanism for the first generation of Population III stars forming at z = 15–30 (e.g., Yoshida et al. 2004; Maio et al. 2010; Johnson et al. 2013; Wise et al. 2014; Jaacks et al. 2018c). Due to the inefficiency of this method, Population III stars are not predicted to contribute significantly to the reionizing budget (e.g., Ricotti & Ostriker 2004; Greif & Bromm 2006; Ahn et al. 2012; Paardekooper et al. 2013); thus, we do not consider star formation in molecular-cooling halos in this work (though see Jaacks et al. 2018b).

Once the IGM begins to be photoheated, even atomic cooling halos will begin to have their star formation suppressed. For the lowest-mass halos, ionization fronts in reionized regions will suppress star formation in minihalos with log(M_h/M_⊙) ≲ 8 (e.g., Shapiro et al. 2004). While more massive halos may self-shield against this process, gas will not accrete onto dark matter halos with virial temperatures less than the IGM temperature through Jeans filtering. Simulations predict that the halo mass where this process begins to dominate is around log(M_h/M_⊙) = 9, though the predictions are quite uncertain (e.g., Gnedin 2000; Iliev et al. 2007; Mesinger & Dijkstra 2008; Okamoto et al. 2008; Alvarez et al. 2012; Ocvirk et al. 2016, 2018; Dawoodbhoy et al. 2018). Feedback likely also plays a strong role (see Somerville and Davé 2015, and references therein), as these small halos have relatively shallow potential wells, allowing gas to easily be lost. For example, Ceverino et al. (2017) found that stellar feedback causes a flattening in the UV luminosity function at M_UV > −14, or log(M_h/M_⊙) ≈ 9.

The physics here are complicated, but in this analysis, we wish only to capture the broad trend of an evolving halo mass where star formation is suppressed. We allow star formation to occur in halos above the redshift-dependent atomic cooling limit, M_h,atomic, which is given by Equation (26) in Barkana & Loeb (2001), assuming a critical virial temperature for atomic cooling of 10,000 K. After reionization begins, we implement photosuppression below a threshold halo mass due to the rising UV background. However, there are a range of plausible limiting halo masses for this photosuppression to take effect (e.g., Iliev et al. 2007; Mesinger & Dijkstra 2008; Okamoto et al. 2008; Alvarez et al. 2012). In addition, even once gas halts accreting, these galaxies may still form stars for a period of time until they use up all of their previously accreted gas (Sobacchi & Mesinger 2013). We approximate these uncertainties by adding this photosuppression mass as a free parameter, M_h,supp, with an adopted flat prior of log(M_h,supp/M_⊙) ∈ (8.5, 10.5) encompassing the range found in the literature. The lower bound was set so that this mass threshold was never lower than the atomic cooling limit at the redshifts considered here; this could have been avoided by allowing the photosuppression mass to be redshift-dependent, but we elected to choose a fixed value in the absence of evidence that this redshift dependence was needed and to avoid adding another free parameter to our model. The nonionizing specific UV luminosity density (ρ_UV) is calculated at each redshift as the integral of the UV luminosity function down to the magnitude corresponding to this limit, shown in the right panel of Figure 1. In Section 3.1 we describe how our model transitions from the atomic cooling limit to the photosuppression mass as reionization progresses.

We reiterate that while our model does not include star formation beyond the limits specified here, a number of recent simulations show star formation, especially in the pre-reionization universe, in very low-mass halos of 7 < log M_h/M_⊙ < 8. However, modern high-resolution simulations still predict a turnover in the UV luminosity function at magnitudes corresponding to approximately the atomic cooling limit (e.g., Wise et al. 2014; Jaacks et al. 2018b). While the flat luminosity function beyond this turnover implies that star formation activity is occurring in lower-mass halos, the shallowing of the luminosity function slope results in these small systems contributing little to the integrated UV luminosity density. As Jaacks et al. (2018b) showed in their Figure 17, although their UV luminosity function continues to M_UV > −8, the UV luminosity density asymptotes to a constant value when integrating to M_UV > −13. Future iterations of our model can better match these theoretical results by including a turnover in our luminosity function, and the results may not be inconsequential, as these extremely low-mass halos could have high ionizing photon escape fractions.

To convert the total nonionizing UV luminosity density to the ionizing emissivity, a value for the ionizing photon production efficiency (${\xi }_{\mathrm{ion}}$) needs to be assumed. This parameter encompasses all of the physics of the underlying stellar population, many of which likely evolve with redshift. For example, the mean metallicity of young stars in galaxies likely decreases from low to high redshift, observationally tracked by a decrease in the typical dust attenuation (e.g., Bouwens et al. 2012, 2014; Finkelstein et al. 2012b), leading to hotter stellar photospheres as the metal opacity (mostly due to iron) is lower and thus a higher ionizing-to-nonionizing UV photon ratio. Another factor to be considered is the effect of binary stars. Stellar population synthesis models that include binary stars (Eldridge & Stanway 2009) show that the ionizing flux is boosted by ∼60% (at low metallicities of Z < 0.3 Z_⊙) compared to models with isolated stars only (Stanway et al. 2016) due to both a harder ionizing spectrum from the primary star (which has its envelope stripped) and an increase in mass for the secondary star, allowing more massive stars to exist at later ages.

These effects certainly play a role in high-redshift galaxies, where we cannot directly probe the ionizing flux. However, the production rate of ionizing photons can be inferred via the detection of nebular emission lines. Bouwens et al. (2016b) inferred Hα emission line fluxes from Spitzer/IRAC photometry at 3.8 < z < 5.0, finding log ξ_ion = ${25.34}_{-0.08}^{+0.09}$ erg⁻¹ Hz, consistent with typically assumed values of log ξ_ion ∼ 25.2–25.3 erg⁻¹ Hz in previous reionization studies (e.g., Madau et al. 1999; Finkelstein et al. 2012a; Kuhlen & Faucher-Giguère 2012; Robertson et al. 2015). At 5.1 < z < 5.4, Bouwens et al. (2016b) found log ξ_ion = ${25.48}_{-0.23}^{+0.29}$ Hz erg⁻¹, hinting at evolution toward larger values at higher redshift, though not at a significant level given the observational uncertainties. Bouwens et al. (2016b) also found evidence that the bluest galaxies exhibit even higher values of ξ_ion, with log ξ_ion = 25.9${}_{-0.2}^{+0.4}$ erg⁻¹ Hz for galaxies in their 5.1 < z < 5.4 sample with β < −2.3 (similar results are found for the faintest galaxies in that sample, which, as shown by Bouwens et al. 2014, are likely also the bluest).

Stark et al. (2015b, 2016) measured ξ_ion via ionized carbon emission, finding ξ_ion = ${25.68}_{-0.19}^{+0.27}$ erg⁻¹ Hz in a lensed galaxy at z = 7.045 with an intrinsic M_UV = −19.3 and ξ_ion = 25.6 for three luminous (M_UV = −22) galaxies at z = 7.15, 7.48, and 7.73. Lastly, Wilkins et al. (2016b) investigated the range of ξ_ion expected from galaxies in the epoch of reionization based on the BlueTides simulation, finding that simulated galaxies spanned the range 25 < ξ_ion < 26, with the highest values obtained when assuming low-metallicity stellar population models that include binaries.

Taken together, this evidence implies that ξ_ion likely depends on redshift and luminosity, which we allow in our model via two free parameters, a redshift dependence dlog ξ_ion/dz and a magnitude dependence dlog ξ_ion/dM_UV. We assume that at our lowest redshift considered of z = 4, galaxies brighter than M_UV = −20 have log ξ_ion = 25.34 erg⁻¹ Hz, consistent with the results from Bouwens et al. (2016b) for this redshift and luminosity. Galaxies at higher redshifts and/or fainter luminosities have values of ξ_ion corresponding to

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{log}\,{\xi }_{\mathrm{ion}}(z,{M}_{\mathrm{UV}}) & = & 25.34+(z-4)\displaystyle \frac{d\mathrm{log}{\xi }_{\mathrm{ion}}}{{dz}}\\ & & +({M}_{\mathrm{UV}}-{M}_{\mathrm{UV},\mathrm{ref}})\displaystyle \frac{d\mathrm{log}{\xi }_{\mathrm{ion}}}{{{dM}}_{\mathrm{UV}}},\end{array}\end{eqnarray} \tag{ 1 }$$

where M_UV,ref is the reference magnitude of −20. We assume flat priors on both free parameters of dlog ξ_ion/dz ∈ (0.0, 0.4) and dlog ξ_ion/dM_UV ∈ (0.0, 0.2). Also, ξ_ion(z, M_UV) was constrained to have a maximum value of 26.0, which corresponds to the highest value seen observationally or from simulations (e.g., Bouwens et al. 2016b; Wilkins et al. 2016b; Izotov et al. 2017).

To derive the ionizing emissivity available to ionize the IGM, we must combine our intrinsic ionizing emissivity (${\rho }_{\mathrm{UV}}\times {\xi }_{\mathrm{ion}}$) with a model for f_esc, which we anchor in the results from the high-resolution First Billion Years (FiBY) simulations. Previous observational studies typically consider only single values, ranging from 10%–50% (e.g., Finkelstein et al. 2010, 2012a, 2015b; Robertson et al. 2013, 2015; Bouwens et al. 2015b). However, as discussed in Section 1, essentially all observations of escaping ionizing radiation from galaxies (albeit at lower redshifts) imply smaller escape fractions. Therefore, rather than assume a single arbitrary value, here we draw on information provided by simulations. While a number of simulations over the past several years have derived this quantity (e.g., Razoumov & Sommer-Larsen 2006; Gnedin et al. 2008; Yajima et al. 2011; Kim et al. 2013; Kimm & Cen 2014; Wise et al. 2014), here we use the high-resolution radiative transfer simulations of ionizing photon escape from Paardekooper et al. (2015). These simulations were post-processed on outputs from the FiBY simulation suite, which follows the formation of the first stars and galaxies from cosmological initial conditions and leads to a realistic galaxy population at z = 6 (S. Khochfar et al. 2019, in preparation). The escape fraction of ionizing photons was determined in more than 75,000 halos by post-processing the highest-resolution FiBY simulations with high-resolution radiative transfer simulations. The radiative transfer simulations were run at the same resolution as the hydrodynamics, allowing the average densities within giant molecular clouds in which stars are born to be resolved. This is essential for the determination of the escape fraction.

Comparing to a number of galaxy properties, this study found that the ionizing escape fraction is strongly anticorrelated with the halo mass. We use their results from all halos that are forming stars, which results in an escape fraction versus halo mass relation that is independent of redshift (we will address photoionization feedback below). At each halo mass (in steps of log[M_h/M_⊙] = 0.5), we compute the distribution of the escape fraction via the kernel density estimation (KDE), using cross-validation to compute the optimal KDE bandwidth. We used 20-fold cross-validation, optimizing how well the KDE fits the remaining data. In every step of the cross-validation, the KDE is constructed on 19/20 parts of the data, and the log-likelihood of the remaining 1/20 part of the data fitting this KDE is computed. That is done 20 times (every time changing which part of the data is left out), and the result is averaged. This procedure is repeated for different values of the bandwidth, and the bandwidth with the best score has been chosen. While a larger bandwidth would result in a smoother distribution, it would not fit the edges of the distribution as well as our adopted bandwidth. Our adopted escape fraction distributions are shown in Figure 2. These escape fractions are effectively time-averaged, as the distributions shown are the average of the instantaneous escape fractions of every halo in the simulation. This figure highlights that the escape fraction distributions are quite broad, but only halos with log(M_h/M_⊙) < 8.5 have more than half of their distribution at f_esc > 1%.

Zoom In Zoom Out Reset image size

At log(M_h/M_⊙) = 9, there is a small probability that ${f}_{\mathrm{esc}}$ ≫ 10%; these few halos are undergoing an extreme starburst, and the supernova feedback is able to evacuate almost all of the gas. As the simulation does not have a representative sample of halos with log(M_h/M_⊙) > 9.5, we assume that halos with log(M_h/M_⊙) > 11 have f_esc = 0, and that halos with log(M_h/M_⊙) = 9–11 have a similar distribution as those at halos with log(M_h/M_⊙) = 9 but without the small peak at high f_esc (due to the increased potential making it more difficult for supernovae to remove all the gas), as shown by the gray dashed line in Figure 2. We note that if we treat halos with log(M_h/M_⊙) > 11 the same as those with log(M_h/M_⊙) = 10, we find almost no differences in the resulting ionization history (completing at z = 5.7, compared to z = 5.5 for our fiducial model), although in the post-reionization universe, the galaxy emissivity drops off slightly more shallowly, with a corresponding slight decrease in the needed AGN emissivity (Section 2.2).

The normalization of the escape fraction may be inaccurate in the Paardekooper et al. (2015) simulations for several reasons. The resolution of the simulations is insufficient to resolve the birth cloud of the star particles in great detail, potentially missing physics on the scale of individual stars that can affect the escape fraction. Simulations have shown that better resolving the interstellar medium (ISM) around the stars results in a higher escape fraction because the porosity of the gas is better accounted for (Paardekooper et al. 2011). In addition, the stellar population model in their simulations does not include the effects of binary interaction, such as mass transfer between stars and mergers of binaries. These processes have been shown to affect the average escape fraction in a halo because massive stars in binaries live longer and thus emit many ionizing photons when the birth cloud of the stars has been dissolved by supernova explosions of the single massive stars in the population (Ma et al. 2016).

We thus adopt an escape fraction "scale factor," where in a given iteration of our model, the escape fractions at all halo masses are scaled by the same factor, preserving the halo mass dependence of the escape fraction. We do not allow this scale factor to vary with redshift, as the simulations find roughly constant escape fractions at fixed halo mass through the epoch 6 < z < 15, and the expected physical reasons for this scale factor do not depend on redshift. We denote this parameter below as f_esc,scale and adopt a flat prior on f_esc,scale over the range f_esc,scale ∈ (0, 10).

The total ionizing emissivity from galaxies is thus calculated as ${\dot{N}}_{\mathrm{ion},\mathrm{gal}}={\rho }_{\mathrm{UV}}\times {\xi }_{\mathrm{ion}}\times {f}_{\mathrm{esc}}$, where the latter term includes this scale factor.

The black lines in Figure 5 show the posterior distributions of our free parameters for our fiducial model. The median values and 68% confidence ranges are listed in Table 2, and we show the covariances between our model parameters as a triangle plot in Figure 6. Of these seven distributions, three have a clearly defined peak in the posterior distribution within our prior range. The parameter f_esc,scale prefers a value near 5, with a 68% confidence range of 3.3 < f_esc,scale < 7.5. This implies that our model best matches observations when the escape fractions from the simulations are scaled up by a factor of ∼3–8. As discussed in Section 2.1.5, this upscaling is not surprising, as the simulations were unable to resolve the birth cloud of the star particles, so the porosity of the gas may be underestimated (though see also Gnedin & Kaurov 2014, who invoked a subunity scale factor to account for unresolved systems providing excess absorption). In Section 5.2 we will discuss the implications of this upscaling for the average galaxy escape fraction and its evolution with redshift.

Zoom In Zoom Out Reset image size

The parameters related to ${\xi }_{\mathrm{ion}}$ also show a clear peak, with dlog ξ/dz ∼ 0.13 and dlog ξ/dM_UV ∼ 0.07, albeit with broad tails to higher values. These results indicate that, under our model, galaxies must have higher ionizing photon production efficiencies both at higher redshifts and at lower luminosities. These results are broadly consistent with both the scarce observations at high redshift and the observations from local analogs for high-redshift galaxies, which we discuss in Section 5.5.

The remaining four parameters had more one-sided distributions. The filtering mass has a clear preference for the lower end of our allowed range, with a one-sided 84% upper limit of log(M_h,supp/M_⊙) < 9.5. This indicates that the model requires star formation in halos as small as possible for as long as possible. We note that our prior of log(M_h,supp/M_⊙) > 8.5 was set so that the filtering mass did not drop below the atomic cooling mass throughout our redshift range. It is plausible that if we reduced this prior our model would prefer even lower values, though, as discussed below in Section 5.4, lower values would be in greater conflict with previous simulation work on the post-UV background filtering mass.

The final three parameters govern the AGN contribution. Combined, they prefer an AGN contribution that evolves relatively shallowly downward with increasing redshift (though not as shallowly as the Madau & Haardt 2015 model), with AGNs present to as high a redshift as allowed and with a high escape fraction for ionizing photons produced by AGNs. In Section 5.3.1, we compare the resultant emissivity from these combined parameters to previous observations. The posterior results from this fiducial model correspond to a DIC value of 5.3.

Figure 7 shows the distribution of both ${Q}_{{{\rm{H}}}_{\mathrm{II}}}$(z) and ${Q}_{{\mathrm{He}}_{\mathrm{III}}}(z)$ for our fiducial model. The blue shaded region shows the 68% confidence range on our evolution of ${Q}_{{{\rm{H}}}_{\mathrm{II}}}(z)$. This result is consistent with the observations from McGreer et al. (2015) used to constrain our model, as the lower 68% of our results fall near the lower 68% confidence on their lower limits at z = 5.6 and 5.9. Our model achieves ${Q}_{{{\rm{H}}}_{\mathrm{II}}}=1$ by z = 5.6 ± 0.5. Although previous studies have adopted a prior that reionization ended by z = 6, it is plausible that reionization ended later. Indeed, significant neutral patches may be necessary to explain the most opaque stretches of the Lyα forest at z = 5.5−6—in particular, the ∼110h⁻¹ comoving Mpc Lyα trough observed by Becker et al. (2015; Kulkarni et al. 2019).

Zoom In Zoom Out Reset image size

The shape of the evolution of ${Q}_{{{\rm{H}}}_{\mathrm{II}}}(z)$ from our fiducial model shows a roughly linear evolution with redshift from z ∼ 6 to 10, with a slight acceleration in the evolution at z > 10 as the faint-end slope of the galaxy luminosity function in our fiducial model stops evolving. This evolution is in sharp contrast to the results of R15, which imply that reionization starts very slowly at z > 10 and then undergoes a rapid acceleration at z < 10. This stark difference can be understood by the differences in these two models. The R15 model fits a variable SFR density to the observations of the galaxy UV luminosity function from a variety of sources and uses as an additional constraint the 2015 value of the Planck optical depth. They convert their derived SFR density to an ionizing emissivity by assuming a single value for the escape fraction (20%) and ionizing photon production efficiency (log[ξ_ion] = 53.14 Lyc photons s⁻¹ M_⊙⁻¹ yr⁻¹; 25.24 in the units used here of Hz erg⁻¹). These assumed values are comparable to those used for the purple shaded region in Figure 21 in the Appendix (f_esc = 13%; log[ξ_ion] = 25.34), which give near-identical results to those from Robertson et al. when using our assumed luminosity functions, highlighting that the minor differences in the luminosity functions assumed (or resultant SFR density) do not play a large role in the differences in ${Q}_{{{\rm{H}}}_{\mathrm{II}}}(z)$.

The differences in ξ_ion and f_esc must therefore be responsible. While both certainly play a role, the differences in ${Q}_{{{\rm{H}}}_{\mathrm{II}}}(z)$ are easy to understand when simply exploring f_esc. In the Robertson et al. model, galaxies at all redshifts and luminosities have the same escape fraction of 20%. This means that the ionizing emissivity from a given luminosity range is directly proportional to the nonionizing UV luminosity density in that same range. The result of this assumption is that the faintest galaxies (M_UV < −14; < 10⁻³L* at z = 6) do not play a major role. This is illustrated in Figure 8, where the gray lines show the cumulative ionizing emissivity using the assumptions from Robertson et al. (albeit with our luminosity functions used here; as discussed in the preceding paragraph, this results in minimal differences). At z = 6, at the end of reionization, ∼half of the ionizing emissivity comes from rather bright galaxies, with L > 0.1L* in the Robertson et al. model. While the evolving faint-end slope changes this with redshift, even by z = 10, a near majority of the ionizing emissivity comes from galaxies with L > 0.01L*. As massive/bright galaxies are building up with time, the relative insignificance of L ≪ 0.01L* galaxies to the ionizing photon budget means that reionization gets a late start. As the brighter/more massive galaxies build up, the ionizing emissivity increases rapidly, resulting in a reionization history that starts late but finishes quickly.

Zoom In Zoom Out Reset image size

The blue shaded lines in Figure 8 show the results of our fiducial model. While the same luminosity functions are assumed as the gray lines, in our model, the escape fraction is halo mass–dependent, so even modestly bright galaxies contribute very little, while the faintest galaxies dominate. This, combined with the steepening faint-end slope at higher redshift, allows reionization to begin earlier. This difference is enhanced by our evolution in ξ_ion to higher values at higher redshifts and lower luminosities, allowing those galaxies that have significant escape fractions to have larger ionizing photon budgets. However, as the faint-end slope shallows with decreasing redshift, the emissivity from galaxies from our model decreases, resulting in a very constrained ionizing photon budget toward the end of reionization, discussed in Section 4.1.4.

While both models can successfully complete reionization, the different assumptions on the escape fraction result in modest differences in the reionization history. These differences are greatest at z ≈ 9, where the Robertson et al. model, driven by modestly bright galaxies, predicts ${Q}_{{{\rm{H}}}_{\mathrm{II}}}$ ∼ 0.2, and our model, driven by the faintest galaxies, predicts ${Q}_{{{\rm{H}}}_{\mathrm{II}}}$ ∼ 0.5. Future observations may be able to distinguish between these scenarios and can thus potentially constrain the luminosity range of the galaxies driving reionization. In Section 5.1 and Figure 11, we further compare our results for ${Q}_{{{\rm{H}}}_{\mathrm{II}}}(z)$ to several observational and theoretical results in the literature.

While we did not include any constraints on the ionization of He ii in our model, here we comment on the resultant distribution of ${Q}_{{\mathrm{He}}_{\mathrm{III}}}(z)$, shown as the red shaded region in Figure 7. Our fiducial model results in an He ii reionization history that gets started at a low level at z ∼ 6, as the AGN helps to complete hydrogen reionization. The volume-ionized fraction of He iii hits 50% at z = 3.4 ± 0.6, and He ii reionization completes at z = 2.7 ± 0.4. This is consistent with observations of the H i and He ii Lyα forests. Current observations of the latter show a strong evolution in its mean opacity and dispersion at z ≳ 2.8 (e.g., Worseck et al. 2011, 2016a). Additionally, the analogs of Gunn–Peterson troughs observed in the z ≳ 2.7 He ii Lyα forest have been used to place limits of Q_{He ii} ≳ 10% at the cosmic mean density, which implies that He ii reionization was likely still in progress at z ∼ 2.7 (McQuinn 2009; Shull et al. 2010; Syphers & Shull 2014). Finally, recent H i Lyα forest measurements have found evidence for a bump in the thermal history of the IGM at z ∼ 2.8, which has been interpreted to coincide with the end of the He ii reionization process (e.g., Schaye et al. 2000; Becker et al. 2011; Puchwein et al. 2015; Upton Sanderbeck et al. 2016; Hiss et al. 2018).

Together, these measurements provide evidence that He ii reionization may be ending at z ∼ 2.5, consistent with our findings (see, however, Davies & Furlanetto 2014; McQuinn & Worseck 2014; Davies et al. 2017, for alternative interpretations of the He ii Lyα forest data). We note that our results also suggest a more extended He ii reionization process than is found in existing simulations—with a 1σ upper bound of as much as 18% complete by z ≈ 6 (McQuinn et al. 2009; Compostella et al. 2013; La Plante et al. 2017). This is qualitatively consistent with the conclusion of Worseck et al. (2016a), who found evidence for extended reionization in the statistics of the He ii Lyα forest opacity.¹⁴

In the left panel of Figure 9, we show the evolution of the ionizing emissivity. The blue shaded region denotes the total H i ionizing emissivity, which is consistent with the observational constraints used for our model. Notably, the emissivity in our model exhibits a slight rise from z = 4 to 10, consistent with the observations from Becker & Bolton (2013) of an increasing emissivity from z ∼ 2 to 5. The purple and red shaded regions denote the components of this emissivity contributed by galaxies and AGNs, respectively. Although the nonionizing emissivity from galaxies increases with decreasing redshift from z = 10 to 4 (Figure 1), the ionizing emissivity does the opposite. This is due to the faint-end slope of the UV luminosity function becoming shallower, resulting in less luminosity coming from faint galaxies, combined with our model where it is those very faint galaxies that have the highest escape fractions. This can be better understood by examining Figure 8. While the cumulative ionizing emissivity at z = 4 is not highly dissimilar to that at z = 10, Figure 1 highlights that the z = 10 universe harbors a greater abundance of extremely faint (M_UV > −15) galaxies than the z = 4 universe. These factors combine to create a galaxy ionizing emissivity that counterintuitively increases with increasing redshift at z < 10. This emissivity peaks at z ∼ 10 and then proceeds to decline to higher redshift. This transition point is set by the assumption in our fiducial model that the faint-end slope does not get any steeper at z > 10, while M* and ϕ* continue to decline to higher redshifts, lowering the overall emissivity.

Zoom In Zoom Out Reset image size

The gray curve shows the emissivity if one assumed a fixed limiting magnitude of M_UV = −13 and escape fraction of 13%. The difference in the ionization history from our model and that from previous results (e.g., Finkelstein et al. 2015a; R15) can be understood by comparing this to the blue curve. Our model has a greater emissivity at z > 9, allowing an earlier start to reionization. The emissivity at z ≤ 7 flattens out just enough to complete reionization by z ∼ 5.6 but not enough to exceed the emissivity observations at lower redshifts. We note that this gray curve, when extrapolated to z < 4, will exceed the emissivity measurements from Becker & Bolton (2013). This indicates that some previous models with fixed large escape fractions may not have matched all observational data when considering z < 6 (see also Stanway et al. 2016).

The AGN emissivity, shown by the red shaded region, increases with decreasing redshift. This is by construction, as our method assumes an AGN emissivity that rises with decreasing redshift, although the exact slope of this increase and the normalization are set by the posterior distributions of these two free parameters. The AGN emissivity, and specifically, how it compares to that from galaxies, is discussed further in Section 5.3.1 below.

The right panel of Figure 9 shows the posterior distribution on τ_es from our fiducial model. The median value from our model is τ_es = 0.071 ± 0.005 (integrated to z = 20).

This is higher than the recent 2016 value published by the Planck Collaboration of ${\tau }_{\mathrm{es}}$ = 0.055 ± 0.009, which we used as our constraint. However, the tension is not high, with the difference being significant at only the 1.55σ level (we note that our result is 0.4σ higher than the 2015 Planck value of τ_es,Planck15 = 0.066 ± 0.012 and 1.2σ lower than the final WMAP9 value of τ_es,WMAP9 = 0.088 ± 0.013). Comparing to R15 (who used the 2015 Planck value), while our model prefers an earlier start to reionization, the large uncertainties on τ_es result in both our model and the results from R15 maintaining consistency with the observational constraints.

Lowering ${\tau }_{\mathrm{es}}$ would require removing electrons somewhere along the line of sight. One could do this by, for example, slightly lowering the galaxy emissivities, resulting in a slightly later end to reionization. However, the emissivities from the model are already at the lower end of the observations, and the ionization history is also pushing the limits of the observations. Thus, while this would decrease the tension between the model τ_es and the observations, it would increase the tension on our other two constraints. Clearly, reduced observational uncertainties on τ_es will help, though the outlook for significant future improvement is uncertain. We present the tabulated values of the results from this fiducial model in Table 3. As listed in the table, the ratio of the galaxy-to-AGN ionizing emissivity becomes infinite at redshifts when there is no AGN contribution.

Table 3.  Results of Fiducial Model

Redshift	QH ii			τe			log ${\dot{N}}_{\mathrm{total}}$ (s−1 Mpc−3)			log ${\dot{N}}_{\mathrm{gal}}$/${\dot{N}}_{\mathrm{AGN}}$
	16%	Median	84%	16%	Median	84%	16%	Median	84%	16%	Median	84%
4.0	1.000	1.000	1.000	0.0222	0.0225	0.0226	50.42	50.64	50.76	−0.81	−0.24	0.32
4.2	1.000	1.000	1.000	0.0237	0.0240	0.0242	50.40	50.62	50.76	−0.78	−0.21	0.44
4.4	1.000	1.000	1.000	0.0253	0.0255	0.0258	50.38	50.62	50.76	−0.72	−0.12	0.53
4.6	1.000	1.000	1.000	0.0268	0.0271	0.0274	50.36	50.59	50.74	−0.78	−0.07	0.63
4.8	1.000	1.000	1.000	0.0284	0.0287	0.0290	50.35	50.60	50.74	−0.72	0.01	0.78
5.0	1.000	1.000	1.000	0.0300	0.0304	0.0307	50.38	50.60	50.74	−0.66	0.06	0.92
5.2	0.993	1.000	1.000	0.0317	0.0320	0.0323	50.34	50.59	50.74	−0.58	0.13	1.01
5.4	0.963	1.000	1.000	0.0333	0.0337	0.0340	50.37	50.60	50.75	−0.60	0.18	1.16
5.6	0.938	1.000	1.000	0.0350	0.0354	0.0358	50.35	50.58	50.76	−0.57	0.23	1.29
5.8	0.905	0.967	1.000	0.0366	0.0370	0.0375	50.39	50.60	50.75	−0.53	0.33	1.41
6.0	0.870	0.930	1.000	0.0382	0.0387	0.0392	50.41	50.62	50.79	−0.44	0.42	1.51
6.2	0.832	0.899	0.996	0.0398	0.0404	0.0409	50.42	50.62	50.78	−0.33	0.44	1.69
6.4	0.795	0.866	0.959	0.0413	0.0420	0.0425	50.44	50.64	50.80	−0.29	0.56	1.93
6.6	0.756	0.836	0.925	0.0428	0.0435	0.0442	50.38	50.61	50.78	−0.27	0.58	2.16
6.8	0.722	0.809	0.892	0.0442	0.0450	0.0458	50.39	50.63	50.79	−0.20	0.66	2.44
7.0	0.694	0.781	0.862	0.0456	0.0465	0.0475	50.43	50.64	50.82	−0.14	0.81	$\infty$
7.2	0.659	0.753	0.836	0.0469	0.0479	0.0490	50.42	50.65	50.82	−0.04	0.93	$\infty$
7.4	0.628	0.724	0.808	0.0482	0.0493	0.0506	50.41	50.65	50.82	0.03	1.11	$\infty$
7.6	0.598	0.699	0.780	0.0494	0.0507	0.0521	50.43	50.66	50.83	0.07	1.21	$\infty$
7.8	0.570	0.672	0.753	0.0506	0.0520	0.0535	50.44	50.67	50.84	0.16	1.42	$\infty$
8.0	0.540	0.647	0.730	0.0517	0.0534	0.0550	50.47	50.69	50.86	0.21	1.74	$\infty$
8.2	0.513	0.617	0.703	0.0529	0.0546	0.0564	50.45	50.68	50.87	0.24	1.92	$\infty$
8.4	0.485	0.590	0.678	0.0539	0.0559	0.0577	50.49	50.70	50.86	0.30	2.25	$\infty$
8.6	0.458	0.563	0.654	0.0549	0.0571	0.0590	50.46	50.68	50.85	0.38	2.46	$\infty$
8.8	0.432	0.533	0.629	0.0558	0.0582	0.0603	50.49	50.72	50.90	0.45	2.83	$\infty$
9.0	0.403	0.502	0.603	0.0568	0.0593	0.0615	50.51	50.72	50.90	0.51	3.24	$\infty$
9.2	0.373	0.476	0.574	0.0577	0.0603	0.0628	50.51	50.74	50.92	0.58	4.42	$\infty$
9.4	0.344	0.443	0.543	0.0584	0.0613	0.0639	50.52	50.74	50.92	0.63	$\infty$	$\infty$
9.6	0.314	0.413	0.507	0.0591	0.0623	0.0650	50.55	50.77	50.95	0.80	$\infty$	$\infty$
9.8	0.285	0.377	0.468	0.0598	0.0632	0.0661	50.57	50.79	50.98	0.90	$\infty$	$\infty$
10.0	0.254	0.342	0.428	0.0605	0.0639	0.0670	50.55	50.79	50.98	1.04	$\infty$	$\infty$
10.2	0.227	0.309	0.389	0.0610	0.0646	0.0679	50.53	50.79	50.96	1.20	$\infty$	$\infty$
10.4	0.202	0.278	0.352	0.0615	0.0653	0.0687	50.51	50.75	50.94	1.32	$\infty$	$\infty$
10.6	0.180	0.247	0.318	0.0620	0.0659	0.0694	50.48	50.73	50.92	1.58	$\infty$	$\infty$
10.8	0.159	0.220	0.287	0.0624	0.0664	0.0701	50.44	50.68	50.90	2.15	$\infty$	$\infty$
11.0	0.140	0.196	0.256	0.0628	0.0669	0.0707	50.42	50.68	50.87	2.78	$\infty$	$\infty$
11.2	0.125	0.174	0.230	0.0632	0.0674	0.0712	50.38	50.62	50.83	4.34	$\infty$	$\infty$
11.4	0.112	0.156	0.207	0.0634	0.0677	0.0717	50.36	50.60	50.81	$\infty$	$\infty$	$\infty$
11.6	0.099	0.139	0.185	0.0637	0.0681	0.0721	50.31	50.57	50.77	$\infty$	$\infty$	$\infty$
11.8	0.087	0.124	0.166	0.0639	0.0684	0.0725	50.29	50.55	50.75	$\infty$	$\infty$	$\infty$
12.0	0.077	0.110	0.147	0.0641	0.0687	0.0729	50.24	50.50	50.71	$\infty$	$\infty$	$\infty$
12.2	0.068	0.098	0.131	0.0643	0.0689	0.0732	50.17	50.44	50.66	$\infty$	$\infty$	$\infty$
12.4	0.060	0.088	0.117	0.0645	0.0691	0.0735	50.18	50.44	50.63	$\infty$	$\infty$	$\infty$
12.6	0.054	0.078	0.104	0.0646	0.0693	0.0738	50.11	50.39	50.60	$\infty$	$\infty$	$\infty$
12.8	0.048	0.070	0.094	0.0647	0.0695	0.0740	50.08	50.33	50.54	$\infty$	$\infty$	$\infty$
13.0	0.043	0.062	0.084	0.0648	0.0697	0.0742	50.05	50.30	50.52	$\infty$	$\infty$	$\infty$
13.2	0.038	0.056	0.075	0.0649	0.0698	0.0744	49.99	50.25	50.47	$\infty$	$\infty$	$\infty$
13.4	0.034	0.051	0.067	0.0650	0.0699	0.0745	49.97	50.22	50.45	$\infty$	$\infty$	$\infty$
13.6	0.031	0.045	0.060	0.0651	0.0701	0.0747	49.98	50.23	50.44	$\infty$	$\infty$	$\infty$
13.8	0.027	0.040	0.053	0.0652	0.0702	0.0748	49.94	50.19	50.40	$\infty$	$\infty$	$\infty$
14.0	0.024	0.036	0.048	0.0653	0.0703	0.0749	49.90	50.14	50.36	$\infty$	$\infty$	$\infty$
14.2	0.021	0.032	0.043	0.0654	0.0704	0.0751	49.88	50.14	50.34	$\infty$	$\infty$	$\infty$
14.4	0.019	0.028	0.038	0.0655	0.0704	0.0752	49.87	50.10	50.29	$\infty$	$\infty$	$\infty$
14.6	0.017	0.025	0.033	0.0655	0.0705	0.0752	49.81	50.05	50.23	$\infty$	$\infty$	$\infty$
14.8	0.015	0.022	0.030	0.0655	0.0706	0.0753	49.76	50.00	50.20	$\infty$	$\infty$	$\infty$
15.0	0.013	0.020	0.026	0.0656	0.0706	0.0754	49.75	50.01	50.19	$\infty$	$\infty$	$\infty$
⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯
20.0	0.000	0.000	0.000	0.0658	0.0710	0.0759	48.70	48.93	49.10	$\infty$	$\infty$	$\infty$

Download table as:  ASCIITypeset image

First, we consider constraints inferred via observations of Lyα emission. It can be an excellent tracer of a neutral IGM, as neutral H i gas resonantly scatters Lyα photons, attenuating their observable surface brightness (e.g., Miralda-Escudé and Rees 1998; Malhotra & Rhoads 2004; Santos 2004; Verhamme et al. 2006; McQuinn et al. 2007; Dijkstra 2014). Constraints on the evolution of the Lyα luminosity function from z = 5.7 to 6.5 have been debated for over a decade, as Malhotra & Rhoads (2004) found no evidence for significant evolution, while Ouchi et al. (2010) found evidence for a significant, albeit mild (∼30% in line luminosity), decrease to z = 6.6.

The advent of wider-area narrowband searches has recently led to improved statistics, notably the HyperSuprimeCam (HSC) Subaru Strategic Program (SSP). The larger areas covered now allow constraints on reionization from not only the luminosity function but also the angular clustering. The first HSC results at z > 6 were published by Ouchi et al. (2018), who studied the evolution of ∼2000 Lyα-emitting galaxies (LAEs) at z = 5.7 and 6.6 over 14–21 deg². They specifically compared the evolution in the angular clustering across this redshift range, finding that LAEs at z = 6.6 appeared slightly more biased than at z = 5.7, inferring constraints on the ionized IGM gas fraction of ${Q}_{{{\rm{H}}}_{\mathrm{II}}}$ = 0.85 ± 0.15 at z = 6.6, consistent with model inferences from previous clustering results (Ouchi et al. 2010; Sobacchi & Mesinger 2015). This same sample was used by Konno et al. (2018) to study the evolution of the Lyα luminosity function. They found a slight evolution in the luminosity function (primarily a factor of ∼2 decrease in ϕ*) and inferred ${Q}_{{{\rm{H}}}_{\mathrm{II}}}$ = 0.7 ± 0.2 at z = 6.6.

Studies of the Lyα luminosity function at higher redshift are being pursued with a variety of telescopes, though they are still nascent. Krug et al. (2012) published results from a narrowband search for LAEs at z = 7.7, finding four candidate galaxies. They concluded that if at least one candidate was a true z = 7.7 LAE, there was no evidence for significant evolution of the LAE luminosity function. Similar results were found at this redshift using data with similar depths by Tilvi et al. (2010) and Hibon et al. (2010). However, Konno et al. (2014) used a Subaru SuprimeCam survey for LAEs at z = 7.3, a factor of ∼2 deeper and wider that the z = 7.7 surveys, to find seven LAEs over 0.5 deg², while 65 would have been expected in the case of no evolution from z = 6.6. They concluded that this was evidence for significant evolution, inferring ${Q}_{{{\rm{H}}}_{\mathrm{II}}}$ = 0.55 ± 0.25 at z = 7.3. Similar results are found when modeling the evolution of the observed Lyα luminosity functions and correlation functions across 5.7 < z < 7.3 from Inoue et al. (2018).

While tenuous, these surveys combine to suggest that the IGM is not completely ionized by z ∼ 7, consistent with our model predictions. However, while the Lyα luminosity function appears to decline, this may not be uniform at all luminosities, as there has been evidence for a bright-end "bump" seen in multiple surveys, where the abundance of LAEs at log L ≳ 43 erg s⁻¹, most spectroscopically confirmed, is higher than expected from a Schechter function fit to lower luminosities (e.g., Matthee et al. 2015; Bagley et al. 2017; Hu et al. 2017; Zheng et al. 2017). One proposed physical explanation is that the neutral fraction inferred by the decline in the overall luminosity function is real, while these bright LAEs live in ionized bubbles, and thus the photons suffer reduced attenuation. Another proposed scenario is that these extremely bright emitters are powered by AGNs. At a much lower redshift of z = 0.3 and 2.2, Wold et al. (2017) and Konno et al. (2016), respectively, found that with a wide enough dynamic range, one could simultaneously measure the Lyα luminosity function from star-forming galaxies as a Schechter function and those from AGN as a power law. Whether this extends to such distant redshifts is unknown, though interestingly, three highly luminous (log L > 43 erg s⁻¹) z ≳ 7 LAEs have tenuous detections of N v (Hu et al. 2017; Laporte et al. 2017; Mainali et al. 2018).

Constraints on the IGM ionization state have also been made using spectroscopic follow-up of Lyα emission from continuum-selected distant star-forming galaxies. For brevity, we focus our discussion on more recent results that have an ever-increasing statistical confidence but acknowledge that a significant amount of work has been done in this area (e.g., Fontana et al. 2010; Stark et al. 2010, Pentericci et al. 2011; 2011; Ono et al. 2012; Rhoads et al. 2012, 2013; Schenker et al. 2012; Caruana et al. 2014). Pentericci et al. (2014) observed 12 z ∼ 7 galaxy candidates in the CANDELS EGS field with the FORS2 optical spectrograph on the Very Large Telescope (VLT). They included data on similarly selected sources from previous observations (Castellano et al. 2010; Fontana et al. 2010; Pentericci et al. 2011; Vanzella et al. 2011; Bradač et al. 2012; Ono et al. 2012; Schenker et al. 2012) to amass a total sample of 68 candidate z ∼ 7 galaxies spectroscopically observed with 8–10 m telescopes. Over this combined sample, they found that Lyα is significantly detected with a rest-frame equivalent width (EW) >25 Å in seven of 39 galaxies at −21.25 < M_UV < −20.25 and five of 25 galaxies at −20.25 < M_UV < −18.75. By comparing to reionization models, they found that this detection fraction is consistent with ${Q}_{{{\rm{H}}}_{\mathrm{II}}}$ < 0.49 (1σ) at z = 7, under the assumption that the IGM is fully ionized at z = 6.

More recently, Pentericci et al. (2018) used spectroscopic observations of 167 z ∼ 6–7 candidates with a large ESO FORS2 program to revisit the Lyα detectability at z ∼ 7 but also to improve constraints at z ∼ 6. Interestingly, this study finds a similar detection fraction in bright galaxies (M_UV < −20.25) at z ∼ 6 and 7 of 10%, while the fraction drops from 35% at z ∼ 6% to 20% at z ∼ 7 in fainter galaxies. Taken together with previous results at z ∼ 5 from Stark et al. (2010), this implies a flat evolution in the Lyα detectability from z = 5 to 6, followed by a milder decline to z = 7 (see also De Barros et al. 2017). This is in contrast to previous results, which found a sharper decrease at z > 6 due to the previous rise in the detectability from z = 5 to 6 (e.g., Stark et al. 2010). This is consistent with their observation of only a slight decrease in the residual flux on the blue side of the line at z ∼ 6 compared to z ∼ 7, as well as only a modest evolution in the EW distribution from z = 6 to 7. Taken together with previous measurements of a rise in the Lyα detectability from z = 4 to 5, this result indicates that the IGM may not be fully reionized by z ∼ 6, necessitating a smaller change in the neutral fraction to z = 7 than when assuming the IGM is fully ionized at z = 6.

Mason et al. (2018b), building on a Bayesian framework introduced in Treu et al. (2012) and Treu et al. (2013), used this new sample of z ∼ 7 galaxies from Pentericci et al. (2018) to infer the neutral fraction of the IGM. In an advance over previous studies, this paper measures the magnitude-dependent evolution of the Lyα EW distribution, rather than just the Lyα detection fraction (see also Jung et al. 2018). They also included the effects of Lyα velocity offsets, which may evolve to lower values at higher redshifts or for lower-mass galaxies (e.g., Song et al. 2014; Stark et al. 2015a; Pentericci et al. 2016; Bradač et al. 2017), allowing a smaller amount of neutral IGM gas to obscure a given Lyα line. They used the sample of z ∼ 6 galaxies from De Barros et al. (2017) and Pentericci et al. (2018) to set a baseline z ∼  6 EW distribution measurement and then explore the evolution to z ∼ 7 using the sample from Pentericci et al. (2014), finding Q_{H ii} = 0.41${}_{-0.11}^{+0.15}$ at z = 7. This result is model-dependent, as it draws sight lines from cosmological simulations to realistically model the impact on Lyα on all scales. Additionally, Mason et al. (2018b) assumed that the IGM was fully ionized at z ∼ 6; thus, if the implication from Pentericci et al. (2018) is correct, the difference in the neutral fraction to z ∼ 7 needed to explain the observations may be lower. Additionally, this paper relied on the relatively few observations of the Lyα velocity offset at high redshift to calibrate their model; larger samples of nonresonant line redshifts are needed to increase the robustness of these calibrations. Finally, improved statistical power on the IGM at z ∼ 7 can be obtained by fitting the evolution of EWs down to lower redshift while also making use of larger samples of z ∼ 7 galaxies (e.g., Pentericci et al. 2018).

As the number of Lyα detections at z > 7.5 is small, less work has been done on the implied neutral fraction at such epochs. Using deep Keck/MOSFIRE integrations of 48 candidate z ∼ 7–8 galaxies in the CANDELS GOODS-N field, Finkelstein et al. (2013) explored how the single detection at z = 7.51 compared to the expectation from the predicted Lyα EW distribution from Stark et al. (2011), finding that this EW distribution could be ruled out at 2.6σ significance. Song et al. (2016) used Keck/MOSFIRE observations of 12 z ∼ 7–8 candidates in the CANDELS GOODS-S field to perform a similar analysis, finding that their single detection at z = 7.66 also ruled out no evolution in the EW distribution from z = 6 (at 1.3σ significance). A more detailed Bayesian analysis with the data from Finkelstein et al. (2013) was performed by Tilvi et al. (2014), conservatively finding that Q_{H ii} < 0.7 at z ∼ 7.5.

Together, these observational results on Lyα emission (luminosity functions, clustering, and spectroscopic follow-up) tell a coherent story: the IGM must be significantly neutral by z ∼ 7. Taken at face value, they imply tension with our fiducial model, which gives Q_{H ii} = 0.78 ± 0.08 at z = 7, though we note that our result at z = 6 of Q_{H ii} > 0.87 is consistent with the slightly later end to reionization discussed by Pentericci et al. (2018). This comparison is highlighted in Figure 11, which shows our fiducial model compared to a number of the Lyα-based constraints discussed here.

However, there are a few reasons to be cautious about these interpretations, especially those that infer very large neutral fractions (Q_{H ii} <  0.5 at z = 7) from an apparent deficit of observable Lyα emission. First, not all galaxies have been difficult to detect. Several publications have noted that bright (M_UV < −21.5) galaxies with a nonzero IRAC color indicative of strong [O iii] emission have a very high Lyα success rate (e.g., Finkelstein et al. 2013; Oesch et al. 2015; Zitrin et al. 2015). Notably, Roberts-Borsani et al. (2016) identified four such galaxies, which, together with spectroscopy from Zitrin et al. (2015) and Stark et al. (2016), all have detected Lyα emission, with Stark et al. (2016) arguing that these galaxies have boosted Lyα transmission due to inhabiting ionized bubbles, consistent with similar inferences from the high-luminosity "bump" seen in z > 6 Lyα luminosity functions. Similar results are also seen for bright galaxies at z = 6–6.5 (Curtis-Lake et al. 2012). Such a physical scenario may also explain the Lyα line discovered by Larson et al. (2018) at z = 7.45, which has an EW of 140 Å, more than double that of any other detected Lyα line at z > 7. Mason et al. (2018a) recently considered this from a modeling perspective, finding that the data are consistent with bright galaxies having a greater Lyα transmission fraction through the IGM (by ∼2×). They found that if the samples are large enough, these bright galaxies can accurately constrain the neutral fraction, which presents a compelling argument for targeted follow-up of bright galaxies at z > 7. However, while the modeling results by Weinberger et al. (2018) agree that the IGM transmission fraction is higher for brighter galaxies, they also found that self-shielded neutral gas in the circumgalactic medium (CGM) can attenuate Lyα and that this occurs preferentially in more massive halos. They thus argued that fainter LAEs are a more reliable probe of the IGM state, and that current observations are consistent with ${Q}_{{{\rm{H}}}_{\mathrm{II}}}$ =  45%–75% at z = 7 (for their preferred "very late" and "late" reionization models, respectively).

Second, there is also evidence from recent observations that the previous evolution of Lyα detectability from Stark et al. (2010), which has been used to invoke large neutral fraction evolution, may need to be revised. In addition to the result from Pentericci et al. (2018), Caruana et al. (2018) revisited the Lyα fraction evolution from z = 3–5 with deep VLT/MUSE integral field unit (IFU) spectroscopy, finding no significant evolution from z ∼ 3 to 5 (X_Lyα ∼ 30% across this range). Several spectroscopic surveys are underway to improve these measurements; thus, the near future should yield more well-characterized Lyα fractions and EW distributions at z = 3–6.

Finally, it is worth discussing the several factors working against detecting Lyα in this epoch, even if the galaxies are copiously producing these photons and the IGM is ionized. First, the galaxies targeted spectroscopically can have broad photometric redshift probability distributions, which at z = 7 results in a significant probability that Lyα would be observed both in the optical and in the near-infrared. This would necessitate two observations with different instruments to fully sample the redshift PDF, which is not commonly done. Even within one instrument, an individual grating/filter pair may not fully sample the PDF. While this incompleteness effect can be accounted for, it is reliant on the shape and width of the photometric redshift PDFs, which have not been calibrated at the redshifts of interest. Jung et al. (2018) explored this effect and found that if the redshift uncertainties are increased by 50%, the number of expected Lyα detections for a fixed EW distribution is reduced by more than a factor of 2. This does not include the possibility that some galaxies are contaminants from lower redshift, though this is explored by Pentericci et al. (2011, 2014), who found that significant contamination is unlikely (though this depends strongly on the quality and amount of photometric data available). Second, at z > 5, Lyα is observed at wavelengths significantly contaminated by night-sky emission, rendering much (>50% at R ∼ 2000) of the wavelength space unavailable (except for extremely bright lines). Finkelstein et al. (2013) found that the presence of night-sky emission lines alone reduces the expected number of detections by a factor of ∼3 (see their Figure S5). Last, and easiest to model, is the effect of the depth of the observations, though this could be impacted by inaccurate flux calibration or slit-loss correction.

Nonetheless, as the current observations do indicate some tension with our model, in Section 5.6 we explore how our results change if we include recent Lyα-based constraints during our fitting process.

In addition to the statistical measures of the Lyα and Lyβ dark pixel fraction over multiple combined quasars used by McGreer et al. (2015) to constrain the neutral fraction, individual quasars can also be constraining. Spectroscopy of the region near the Lyα transition can allow one to measure excess absorption redward of the Lyα transition indicative of the presence of neutral gas and thus a strong Lyα absorption damping wing. Here we consider the two presently available damping-wing measurements at z > 7.

The first quasar discovered at z > 7 was published by Mortlock et al. (2011) at z = 7.085. The spectrum of this object, specifically the Lyα damping wing, was further analyzed by Bolton et al. (2011), who found a lower limit on the neutral fraction of >10%, or Q_{H ii} < 90%. These data were further analyzed by Greig et al. (2017), who found 0.39 < Q_{H ii} < 0.79 at 1σ (and 0.19 < Q_{H ii} < 0.92 at 2σ). These constraints are consistent with the 1σ confidence range from our fiducial model at this redshift of 0.68 < Q_{H ii} < 0.85, although the allowable range of ionized fractions from the quasar analysis extends to lower values, as shown in Figure 11. More recently, Bañados et al. (2018) published the discovery of a quasar at z = 7.54, finding a spectral damping wing consistent with 0.23 < Q_{H ii} < 0.62 at 1σ (and 0.06 < Q_{H ii} < 0.83 at 2σ). While the z = 7.085 quasar implies a neutral fraction consistent with our fiducial model, this higher-redshift quasar has a damping wing that implies a higher neutral fraction than our model, although the discrepancy is only at 1.3σ.

However, we note that these measurements are for single lines of sight to these cosmologically biased sources, and with only two lines of sight, it is difficult to make robust conclusions about the global neutral fraction. Additionally, these measurements rely on the ability to model the intrinsic quasar spectrum near the Lyα line, which cannot be directly observed. Finally, the observed signature of neutral gas could also be created if there was a nearby damped Lyα absorber along the line of sight. As discussed by Bolton et al. (2011), if the z = 7.085 quasar had a nearby absorber, the spectrum would be consistent with Q_{H ii} > 0.999 (though see Finlator et al. 2013). Combined with the low significance of the difference between our fiducial model and the quasar results, we conclude that the quasar observations do not rule out our fiducial model, though future observations may necessitate a revision.

Further confidence will be gained as the sample of z > 7 quasars increases, in particular, if they all show a significant damping wing. Similar studies can be pursued with fast follow-up of gamma-ray bursts (GRBs; e.g., Totani et al. 2014), potentially to even higher redshifts. A variety of wide-field surveys are underway, and the future for quasar discovery is bright with the launch in the next decade of the wide-field near-infrared survey telescopes Euclid and WFIRST.

The past few decades have seen tremendous advances in our ability to theoretically model the process of reionization using a variety of methods (e.g., Shapiro et al. 1994; Trac & Cen 2007; Haardt & Madau 2012; Kuhlen & Faucher-Giguère 2012; Gnedin et al. 2017). In Figure 11, we compare our fiducial results to a few recent simulation results using different techniques. The gray line shows the result from R15, who showed that one could match both the Planck Collaboration et al. (2016a) optical depth values and the observed SFR density from high-redshift galaxies with a model where galaxies all have a uniformly high escape fraction of 20% and modest ionizing photon production efficiency of log ξ_ion = 25.2. As discussed in Section 4.1.2, this results in L > 0.1L* galaxies dominating the photon budget, a late beginning to reionization, followed by a rapid ramp-up in the ionized volume-filling fraction at z < 8. While this model has the benefit of exhibiting larger consistency with the Lyα and quasar constraints, it relies on a uniformly high escape fraction, which is unlikely (see Section 5.2). Also shown are the results from the model of Greig & Mesinger (2017), using only Planck τ_es and the McGreer et al. (2015) dark pixel fraction as constraints. While the Greig & Mesinger (2017) model is less empirical than our own, using a set of analytical formulae with three free parameters, they found similar reionization histories as our model when they used similar constraints.

Bouwens et al. (2015a) used a two-parameter model to model the emissivity from galaxies, constraining the redshift evolution of this quantity using a variety of observational constraints. They found that the emissivity must decrease by 0.15 dex per unit of increasing redshift over 6 < z < 10, and that this emissivity is consistent with being fully produced by galaxies if log(f_escξ_ion) = 24.53. For the commonly used fiducial value of log ξ_ion = 25.2, this corresponds to an average escape fraction of ∼20%, similar to that used by R15. For a higher log ξ_ion = 25.5, this can be accomplished with f_esc = 10%. One key difference between the results of Bouwens et al. (2015a) and our own is that our inferred emissivity increases with increasing redshift over 6 < z < 10, which is a consequence of our lower average escape fraction over this epoch (see Section 5.2).

Recently, Rosdahl et al. (2018) explored reionization with the SPHINX simulations, which simultaneously resolve the small-scale physics regulating the ionizing emissivity and begin to approach the larger scales needed to solve for global reionization. This paper focuses specifically on the impact of binary stars, finding that by including binary stars, their full volume ionizes by z ≈ 7; this is due to both the increase in ionizing photon production by massive binaries and the increase in escape fraction (see Section 5.2 for further discussion). Similar to our fiducial model, their simulation begins reionization relatively early. However, their reionization history completes by z = 7. Examining their results, this is likely because their typical escape fractions are ∼10%, whereas the globally averaged escape fractions in our model are rapidly evolving from >10% at z ∼ 15 to <5% at z < 10 (Section 5.2). Therefore, while the early phases of reionization are similar between this model and our own, the end of reionization is more extended in our fiducial model due to the declining galaxy emissivity. Finally, the Technicolor Dawn simulations (Finlator et al. 2018) combine an updated model for galaxy formation and feedback with a multifrequency moment-based radiation transport solver that models reionization and photoionization heating in detail. The redshift-dependent ionizing escape fraction is calibrated to match observations of the optical depth to Thomson scattering, as well as ionizing emissivity at z = 5. As shown in Figure 11, this calibration leads to a reionization history that predicts a similar ionized fraction as our model at z ∼ 7, though it predicts a somewhat lower value of ${Q}_{{{\rm{H}}}_{\mathrm{II}}}$ ∼ 35% at z = 9 (compared to ∼50% for our model).

One of the primary differences of our model from previous analyses is the escape fraction distribution, which is skewed highly toward low halo masses. As shown in Figure 8, this results in the bulk of the ionizing photons coming from extremely faint galaxies. Here we examine what effect this has on the redshift evolution of the global escape fraction, which is shown in Figure 12. The blue shading in this figure shows the posterior distribution at each redshift of the total number of escaping ionizing photons divided by the total number created over all galaxies (i.e., over the full luminosity function) and highlights that this globally averaged escape fraction evolves significantly, from <1% at z = 4 to >10% at z = 15. In this figure, we also separately show the global escape fraction from faint galaxies with ${M}_{\mathrm{UV}}$ > −15 and brighter galaxies with −20 < M_UV < −16.

Zoom In Zoom Out Reset image size

This evolution is entirely driven by the coevolution of the luminosity and halo mass functions. At the highest redshifts, the luminosity function is steep, and these faint galaxies live in very low-mass halos. For example, at z = 10, the bulk of the escaping ionizing photons are produced in halos with M_UV > −15 (shown in purple in Figure 12). These faint galaxies, which, according to our abundance-matching results shown in Figure 1, have log(M_h/M_⊙) ≲ 9, have globally averaged escape fractions of 6%–10%. At this same redshift, brighter galaxies (−20 < M_UV < −16; log M_h/M_⊙ > 9.5) have escape fractions of only 1%–3%, shown by the red shading in Figure 12. While many of these halos should have negligible escape fractions, this nonzero global value is driven by high individual escape fractions (f_esc > 50%) in a small number of starbursting galaxies (see red line in Figure 2 and discussion in Section 2.1.5).

At progressively lower redshifts, the continued evolution of the luminosity function, specifically the flattening of the faint-end slope, results in fewer luminous halos with such low masses, reducing the overall escape fraction. This effect is also observed when exploring the ionizing emissivity in the left panel of Figure 9, as the emissivity from galaxies decreases by nearly 1 dex from z = 10 to 4, even though the nonionizing ρ_UV increases by nearly 1 dex over this same redshift range.

There have been a large number of observational attempts to directly measure the ionizing photon escape fraction. These are difficult endeavors for a variety of reasons. First, corrections need to be made for both the galaxy SED shape and dust attenuation to convert the observed ratio of Lyman continuum–to–1500 Å flux into the total fraction of escaping ionizing photons. Second, the observed ionizing photons have traveled through the IGM, which, while ionized, still absorbs some fraction of the photons, requiring a statistical correction with a large variance, which propagates through into significant uncertainties on the escape fraction. Finally, the IGM optical depth becomes so great at z > 4 as to render direct measures of the escape fraction in the epoch of reionization impossible, with the most distant detection presently at z = 4.0 (Vanzella et al. 2018).

Nonetheless, the importance of observationally constraining the escape fraction has led to a large number of ambitious observational programs. A trend in recent years has been an increasing success rate of robust direct detections of Lyman continuum radiation with f_esc ≳ 10%, both at z ∼ 2–4 (e.g., Steidel et al. 2001; de Barros & Vanzella 2016; Shapley et al. 2016; Bian et al. 2017; Fletcher et al. 2018; Vanzella et al. 2018), and locally (Izotov et al. 2016a, 2016b, 2018). These galaxies were typically selected for follow-up based on very high observed or inferred ionizing environments, typically constrained via the ratio of [O iii]/[O ii] line emission, which has been proposed as a predictor of strong ionizing photon escape (e.g., Jaskot & Oey 2013; Nakajima & Ouchi 2014; Stasińska et al. 2015). However, other studies have explored the characteristic global escape fraction of galaxies by stacking larger numbers of individually undetected galaxies, typically finding no detection even in the stack, setting strong upper limits on the typical escape fraction to be as low as <2% (e.g., Siana et al. 2010; Sandberg et al. 2015; Grazian et al. 2017; Japelj et al. 2017; Rutkowski et al. 2017; Hernandez et al. 2018).

These observations only probe presently observable galaxies and thus are probing relatively massive halos. However, these results qualitatively agree with our assumed escape fraction for such halos. These massive galaxies typically have very low escape fractions, consistent with the nondetections in the stacks of full galaxy samples. However, they do have a small (∼5%) probability of being observed with a very high escape fraction, due in the simulation to starbursts. This could be the origin of the small fraction of galaxies with observed high escape fractions. While this is a tidy explanation, it could also be coincidental. What is truly needed are much more stringent escape fraction measurements for fainter/lower-mass galaxies. This would require a leap in our space-based UV imaging or spectroscopic capabilities, which is the focus of several concepts for future space missions (McCandliss & O'Meara 2017; Scowen et al. 2017). However, fainter galaxies have begun to be probed. By analyzing the H i column densities in GRB host galaxies, Tanvir et al. (2019) found a typical escape fraction for GRB hosts of 0.005. As GRB hosts can be extremely faint, this result implies that the escape fractions from even faint galaxies in our model may be too high.

Finally, we comment on the recent result by Steidel et al. (2018), who derived the escape fraction from bright (${ \mathcal R }$ < 25.5) z ∼ 3 galaxies via very deep spectroscopic follow-up, combined with an intensive modeling effort. In contrast to some previous imaging results, this spectroscopic campaign detected escaping Lyman continuum radiation and derived an average escape fraction of f_esc =  0.09 ± 0.01 (similar results have been found via other methods by Kakiichi et al. 2018 and Fletcher et al. 2018). This average comes from a direct detection from 15 galaxies with f_esc = 0.60 ± 0.06 and a stacked detection from 109 individually undetected galaxies with f_esc = 0.04 ± 0.01. This trend is qualitatively similar to that shown for massive galaxies in Figure 2, where ∼10% of the probability density lies at high escape fractions, though in our fiducial model, the remaining probability density has f_esc ≪ 0.04. They also found a luminosity dependence, where the brightest 50% of their galaxies have escape fractions consistent with zero, and the faintest 50% have f_esc ∼ 0.3 (similar trends are also seen with Lyα EW). They used these observations to conclude that bright galaxies are dominating the emissivity at z ∼ 3 (∼3× that of quasars). While we do not track galaxies to z < 4, extrapolating our results shows that in our fiducial model, quasars should dominate the emissivity at z ∼ 3.

While more work is required to see whether contamination from line-of-sight interlopers is found to be minimal for this sample (available soon from an in-progress HST program [PI: Shapley]) and whether this result holds across different fields and redshifts, we explored how our model would change if we assumed a fixed escape fraction of 9% for all galaxies. Unsurprisingly, this model completes reionization earlier at z = 7.0 ± 0.8, with a total emissivity dominated by galaxies at all epochs that rises continuously with decreasing redshift, reaching log $\dot{N}$ = 51.25 ± 0.1 s⁻¹ Mpc⁻³ by z = 4. This is higher than but consistent within 1σ of the observational constraint at that redshift. Due to the high emissivity from galaxies, the AGN contribution is negligible, with the AGN emissivity matching the lower bound allowed by our model at z < 6. As the emissivity from galaxies should continue to rise to z ∼ 3 in this model, it may begin to be in tension with the measured value of the emissivity at that redshift of log $\dot{N}$ = 50.8 s⁻¹ Mpc⁻³, though our model would need to be extended to lower redshifts to explore this further.

While we have assumed results from a particular simulation for our escape fraction parameterization, here we explore how this differs from other recent simulations, which find escape fractions over nearly the full possible range, often with a mass and/or redshift dependence (e.g., Razoumov & Sommer-Larsen 2006; Gnedin et al. 2008; Wise & Cen 2009; Paardekooper et al. 2011, 2015; Yajima et al. 2011; Kim et al. 2013; Wise et al. 2014). Anderson et al. (2017) studied the escape fraction from several dozen galaxies in their 25 × 12 × 10 Mpc simulation box, finding that their faintest galaxies at −16 < M_UV < −14 had f_esc ∼ 35%, while galaxies at M_UV = −18 had f_esc ∼ 1%. While the trend is qualitatively similar to what we assume, their normalization is higher, such that utilizing the escape fractions from this simulation would put less of an emphasis on the extreme faintest galaxies. However, these results are based on very few (<10) galaxies fainter than M_UV = −18. Additionally, the resolution of this simulation is 350 pc, and unlike the simulations of Paardekooper et al. (2015), the halos of interest were not resimulated at higher resolution, which could imply that the important physical scales for ionizing photon escape were not resolved (perhaps not fully resolving important physical processes such as turbulence; Safarzadeh & Scannapieco 2016). Nonetheless, their conclusion that the faintest galaxies in their simulation dominate the ionizing emissivity is qualitatively consistent with our results.

Xu et al. (2016) measured the ionizing escape fraction using the Renaissance Simulation suite (O'Shea et al. 2015), finding very high escape fractions at z = 8–15 of 40%–60% at log(M_h/M_⊙) = 7, decreasing to ∼5% at log(M_h/M_⊙) = 8–9. Interestingly, they found that this rises to 10%–20% in their few halos log(M_h/M_⊙) ∼ 9.5, presumably due to starburst activity. These results are in very good agreement with our assumptions, highlighting again the dependence of extremely faint/low-mass galaxies to account for the ionizing budget. Alternatively, Sharma et al. (2016) found that if all ionizing photons escape when the SFR surface density exceeds a critical threshold, then it is the brighter galaxies that dominate the ionizing photon budget, as they exceed this critical threshold more frequently.

Ma et al. (2015) explored ionizing photon escape using the FIRE simulations, which include advanced treatment for feedback and utilize zoom-ins to achieve high (<1 pc) resolution on halos of interest. In this work, they found no dependence of f_esc on halo mass, but they did find that the time-averaged escape fraction is <5%. In a follow-up paper, Ma et al. (2016) used these same simulations but explored the impact of assuming binary stellar population synthesis models. Mass transfer between binaries can result in massive stars having significantly longer lifetimes, with the most massive ionizing photon–producing stars extended to >3 Myr, long enough for the "birth cloud," responsible for absorbing most of these photons, to disperse (e.g., Eldridge & Stanway 2009). They found that this effect results in an increase of the time-averaged escape fraction by a factor of ∼5–10, with some massive simulated halos exhibiting f_esc > 10%. This effect is amplified by the increased ionizing photon output of massive stars in binaries, which exhibit hotter stellar photospheres due to the mass transfer. Although these results are only available for a small number of simulated halos, it heavily suggests that future simulations should account for the impact of binary stars.

As described in Section 2.2, we allow a contribution to the ionizing emissivity from AGNs, with this emissivity allowed have a redshift evolution anywhere from the steep Hopkins et al. (2007) evolution, which implies a minimal contribution during reionization, to the shallower Madau & Haardt (2015) evolution, which would result in AGNs dominating reionization. The shaded region in Figure 3 shows the 68% confidence range of our fiducial model compared to both of these previously proposed trends, as well as data from the literature. Our fiducial model prefers an AGN ionizing emissivity in between these previous trends at z ∼ 4, with a redshift evolution slope similar to that of Madau & Haardt (2015), albeit at a lower normalization.

The posterior distribution of the three parameters that govern this emissivity evolution is shown by the black lines in Figure 5. The median values of these parameters are a maximum redshift for AGNs of z_AGN,max = 9.2, a slope of the emissivity with a redshift of AGN_slope = −0.39, and a normalization factor of AGN_scale = 0.77. However, none of these posterior distributions have a well-constrained central value, so we can more appropriately place one-sided 84% lower limits of z_AGN,max > 6.9, AGN_slope > −0.93, and AGN_scale > 0.47. While these are broad lower limits, they do constrain the emissivity from AGNs to have a slope shallower than that from Hopkins et al. (2007) kicking in at some point during the epoch of reionization.

The AGN scale factor could be nonunity for a variety of reasons, but this is most analogous to an ionizing photon escape fraction for AGNs. This is often assumed to be unity (e.g., Madau & Haardt 2015), which is reasonable for very bright quasars, as the energetics near the accreting supermassive black holes likely create channels for ionizing photon escape. However, if the AGN emissivity evolves as shallowly as suggested by our model, it is not these rare systems that are dominating the emissivity; thus, more pressing is the escape fractions from fainter AGNs. Most recently, Grazian et al. (2018) spectroscopically observed 16 faint AGNs blueward of the Lyman continuum break and significantly detected ionizing flux from every object, implying AGN ionizing photon escape fractions ranging from 44% to 100%, with a median of 74%, consistent with the posterior distribution from our model (see also Smith et al. 2018). It is worth noting that the escape fraction for ionizing photons from AGNs and massive stars need not be the same in a halo of fixed mass, since the central power source of the AGN may clear a channel for those centrally created photons, while massive stars farther from the center still may be subject to a high H i column density.

In Figure 13, we compare the emissivity from galaxies to that from AGNs in our fiducial model, with the shading denoting the 68% confidence range, which is somewhat broad, as it includes the uncertainty in both distributions. However, this plot makes clear that while the emissivity from AGNs is significant, galaxies are still the dominant driver behind reionization. The most conservative statement we can make is that at 68% confidence, galaxies have a higher ionizing emissivity than AGNs at z > 7.0. However, even at z = 6.0, our model prefers galaxies to dominate, with log(${\dot{N}}_{\mathrm{gal}}$/${\dot{N}}_{\mathrm{AGN}}$) = 0.5${}_{-0.4}^{+1.2}$. As low as z = 4.6, our model still implies that half of the ionizing photons are coming from galaxies, with log(${\dot{N}}_{\mathrm{gal}}$/${\dot{N}}_{\mathrm{AGN}}$) = 0.0${}_{-0.3}^{+0.7}$, and it is not until lower redshifts that the median of the posterior of our model crosses into the AGN-dominated regime. These results highlight that our fiducial model of reionization is "AGN-assisted" rather than "AGN-dominated."

Zoom In Zoom Out Reset image size

Qualitatively similar models of the AGN ionizing emissivity were considered in previous studies (D'Aloisio et al. 2017; Mitra et al. 2018; Puchwein et al. 2019). All of these authors found that AGN-dominated scenarios were disfavored by contemporaneous Lyα forest temperature measurements (we shall expand upon this topic below). This is consistent with the results of Finlator et al. (2016), who found that an AGN-dominated UV background would fail to reproduce the ratios of observed CGM metal absorber column density distributions at z ∼ 6, while a small AGN contribution could not be ruled out. Mitra et al. (2018) also examined the contribution of quasars to reionization and similarly ruled out an AGN-dominated scenario for hydrogen reionization. They utilized observations of the He ii Lyα forest to constrain their models, finding that the contribution of AGNs to the ionizing background must be negligible for z > 6, preferring a lower value of the AGN ionizing emissivity at z = 4–6 than that of our model.

A significant contribution of AGNs to the z > 4 ionizing background, while contrary to previous expectations, can potentially resolve some interesting recent observations. First, large-scale opacity fluctuations have been observed in the Lyα forest at z = 5–6 (Fan et al. 2006; Becker et al. 2015; Bosman et al. 2018). Chardin et al. (2017) proposed that this could be evidence of AGNs dominating the emissivity and found that a ≥50% contribution of AGNs at these redshifts could explain these results, which are similar to, though in excess of, our model results. However, this is not a unique explanation. For example, D'Aloisio et al. (2015) suggested that some of the opacity fluctuations could be due to residual temperature fluctuations imprinted by the patchy reionization process (see, however, Becker et al. 2018). Davies & Furlanetto (2016) showed that galaxies alone could generate large fluctuations in the ionizing background, if there are strong variations in the mean free path from location to location (see also D'Aloisio et al. 2018). Second, Worseck et al. (2016b) observed the He ii Lyα forest at 2.3 < z < 3.5 and found evidence that He ii must be significantly ionized by z ∼ 3.4, earlier than previous results, which suggested z < 3. As shown in Figure 7, our model predicts ${{\rm{Q}}}_{{\mathrm{He}}_{\mathrm{III}}}$ ∼ 0.60 at z = 3.4, consistent with this result (see discussion also in Section 4.1.3). Finally, several recent spectroscopic studies of high-redshift galaxies have detected potential N v emission, an energetic transition that cannot be produced via starlight and is thus indicative of AGN activity (Hu et al. 2017; Laporte et al. 2017; Mainali et al. 2018). While these detections are tenuous, the impending spectroscopy of a large sample of reionization-era galaxies with the James Webb Space Telescope (JWST) should further probe this line of evidence.

One way to constrain these AGN-assisted models is through the thermal history of the IGM. In general, models with a larger contribution from AGNs at earlier times produce an earlier onset of He ii reionization. The additional energy injected into the IGM by this process can, in principle, be measured by its impact on the small-scale statistics of the H i Lyα forest. Unfortunately, modeling thermal histories for a large grid of models is currently too computationally intensive to be incorporated into our MCMC analysis. In addition, Lyα forest temperature measurements are still subject to potentially large systematic uncertainties, which may be reflected in the level of discord among existing measurements (see below). Rather than incorporate these constraints into our pipeline, we provide some illustrative calculations comparing our models to some of the most recent temperature measurements.

To compute the IGM thermal history for a given model, we adopt the approach of D'Aloisio et al. (2017), which is based on the multizone model of Upton Sanderbeck et al. (2016). We refer the reader to those papers for technical details. In summary, we track the temperatures of an ensemble of gas parcels at various initial densities. These densities evolve according to the Zel'dovich approximation, and each gas parcel is impulsively heated at a different time to mimic the patchy reionization process. To model the effects of the ionizing background generated by AGNs, heating from He ii reionization is separated into two regimes. In addition to impulsive heating, which mimics the effects of He iii ionization fronts sweeping through the IGM, gas parcels are also slowly heated by a uniform EUV/X-ray background that is built up over time with the rise of the AGN population. This multizone model was designed as an approximation to the heating effects that are observed in radiative transfer simulations of the reionization process. Indeed, Garaldi et al. (2019) found good agreement between our modeling and their radiative transfer simulations of an AGN-driven reionization.

Figure 14 shows the results of this modeling.¹⁵ The top panel shows the evolution of the IGM temperature at the cosmic mean gas density, T₀. The blue curve and shaded region in the top panel correspond to our median and 68% C.L. values of the ionizing emissivity, respectively. The right and left blue regions in the bottom panel show the corresponding volume-filling factors of H ii and He iii. In the top panel, we compare our calculations to several recent Lyα forest temperature measurements in the literature. We note that the curvature method employed by Becker et al. (2011) and Boera et al. (2014) probes the temperature at a characteristic gas density that is different from the mean. To compare them against our T₀ calculations, we have extrapolated these measurements to the mean density using the temperature–density relation of the median model.

Zoom In Zoom Out Reset image size

The thermal history contains two peaks associated with the completion of the H i and He ii reionization processes. Figure 14 shows that a larger contribution of AGNs to the z > 4 ionizing background leads to an earlier onset of He ii reionization, which, in turn, leads to an earlier heating of the IGM. This earlier heating is discrepant with the lower z ∼ 5 temperatures measured by Becker et al. (2011), Boera et al. (2019), and Walther et al. (2018), which imply a later and more rapid rise of ionizing emissions from AGNs compared to our models—consistent with the conclusions of previous studies (D'Aloisio et al. 2017; Mitra et al. 2018; Puchwein et al. 2019). For reference, the red dashed curves show a minimal AGN model based on the QSO luminosity function of Hopkins et al. (2007; see next subsection for more details). Previous studies have noted that the Hopkins et al. (2007) luminosity function yields temperatures that are generally consistent with the Becker et al. (2011) and Boera et al. (2014) measurements (Puchwein et al. 2015; Upton Sanderbeck et al. 2016; D'Aloisio et al. 2017).

The higher temperatures implied by the fiducial AGN-assisted model are in better agreement with the recent z < 3.5 measurements of Hiss et al. (2018). In principle, all Lyα forest temperature measurements attempt to extract the effects of thermal broadening on the forest absorption features. However, in practice, the discrepant sets of measurements shown in Figure 14 utilize different techniques for isolating those effects. The curvature method of Becker et al. (2011) probes the shape of the absorption features using the second derivative of the flux, while the measurements of Hiss et al. (2018) are based on fitting Lyα absorption lines to extract Doppler parameters. On the other hand, the most recent measurements of Boera et al. (2019) and Walther et al. (2018) are obtained from the shape of the flux power spectrum. The disagreement between the measurements is likely due (at least in part) to different systematics between the techniques. Clearly, it will be important for future studies to determine the cause(s) of these discrepancies. In addition to constraining He ii reionization, pinning down the thermal history at redshifts z = 2–5 will have important implications for the epoch of hydrogen reionization.

With the thermal history results of the previous section in mind, it is interesting to explore how our model changes if we fix the AGN emissivity to follow the redshift evolution of Hopkins et al. (2007). These results are shown in Figure 15. The top panels show the distribution of the four non-AGN-related free parameters compared to our fiducial model. Understandably, without a significant contribution from AGNs, these parameters are skewed (albeit only slightly) toward values that promote higher ionizing emissivity from galaxies—higher escape fractions, higher ξ_ion values, and a lower photosuppression mass. The bottom row compares the results of this model to our constraints, which can be directly compared to our fiducial model in Figures 7 and 9, with the gray transparent shading showing those fiducial results. In the left panel, one can see that with a minimal AGN contribution, the reionization history is consistent with our fiducial result, though it is shifted to the upper end of the fiducial posterior. This model requires more ionizing photons from galaxies than the fiducial model, which results in a shift of the free parameters to create those higher emissivities, creating more ionizing photons at higher redshifts where galaxies have (on average) the highest escape fractions, contributing to the nudge to higher ionized fractions. However, with minimal AGNs, the slope of the reionization history exhibits a slower evolution at z < 8, resulting in Q_{H ii} = 1 at the lower redshift of z = 5.3${}_{-1.0}^{+0.7}$. This model has a tail extending to reionization completion down to the low redshift of z = 4.3.

Zoom In Zoom Out Reset image size

This very slightly earlier onset of reionization increases the electron scattering optical depth by Δ τ =  0.001. This model increases tension with the ionizing emissivity measurements at lower redshift, shown in the middle panel. Though the free parameters conspire to increase the ionizing emissivity at higher redshifts, the flattening faint-end slope toward lower redshifts results in very low escape fractions at z < 6, so galaxies cannot account for the loss of ionizing photons from AGNs in our fiducial model at 4 < z < 6, resulting in a modest tension with the observations at z = 4 and 4.75.

Given this slightly increased tension, it is instructive to compare the goodness of fit from this model to our fiducial model. We find DIC = 6.4 for the minimal AGN model, compared to 5.3 for our fiducial model. We conclude that altering our model to have a minimal level of AGN contribution to reionization is not ruled out, though it does increase the tension with observations slightly along multiple axes.

Specifically, the portion of the posterior that has reionization finishing at $z$ < 5 is robustly ruled out via observations (e.g., Becker et al. 2018). Since these observations were not used as constraints on our model, this tension does not have high statistical significance (with ΔDIC only 1.1), predominantly due to the large observational uncertainties.