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

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

[    Anson D’Aloisio    Jan-Pieter Paardekooper    [    Peter Behroozi    [    Rachael Livermore    Phoebe R. Upton Sanderbeck    [    Sadegh Khochfar
Abstract

We explore scenarios for reionizing the intergalactic medium with low galaxy ionizing photon escape fractions. We combine simulation-based halo-mass dependent escape fractions with an extrapolation of the observed galaxy rest-ultraviolet luminosity functions to solve for the reionization history from 20 4. We explore the posterior distributions for key unknown quantities, including the limiting halo mass for star-formation, the ionizing photon production efficiency, and a potential contribution from active galactic nuclei (AGN). We marginalize over the allowable parameter space using a Markov Chain Monte Carlo method, finding a solution which satisfies the most model-independent constraints on reionization. Our fiducial model can match observational constraints with an average escape fraction of 5% throughout the bulk of the epoch of reionization if: i) galaxies form stars down to the atomic cooling limit before reionization and a photosuppression mass of log (/M) 9 during/after reionization (13 11); ii) galaxies become more efficient producers of ionizing photons at higher redshifts and fainter magnitudes, and iii) there is a significant, but sub-dominant, contribution by AGN at 7. In this model the faintest galaxies (15) dominate the ionizing emissivity, leading to an earlier start to reionization and a smoother evolution of the ionized volume filling fraction than models which assume a single escape fraction at all redshifts and luminosities. The ionizing emissivity from this model is consistent with observations at 4–5 (and below, when extrapolated), in contrast to some models which assume a single escape fraction. Our predicted ionized volume filling fraction at 7 of 78% ( 8%) is in modest (1–2) tension with observations of Ly emitters at 7 and the damping wing analyses of the two known 7 quasars, which prefer 40–50%.

early universe — galaxies: reionization — galaxies: formation — galaxies: evolution

0000-0001-8519-1130]Steven L. Finkelstein \move@AU\move@AF\@affiliationDepartment of Astronomy, The University of Texas at Austin, Austin, TX 78712, USA stevenf@astro.as.utexas.edu \move@AU\move@AF\@affiliationUniversity of California, Riverside, CA, USA \move@AU\move@AF\@affiliationUniversität Heidelberg, Zentrum für Astronomie, Institut für Theoretische Astrophysik, Albert-–Ueberle-–Str. 2, 69120 Heidelberg, Germany 0000-0003-0894-1588]Russell Ryan Jr. \move@AU\move@AF\@affiliationSpace Telescope Science Institute, Baltimore, MD, USA \move@AU\move@AF\@affiliationUniversity of Arizona, Tucson, AZ, USA 0000-0002-0496-1656]Kristian Finlator \move@AU\move@AF\@affiliationNew Mexico State University, Las Cruces, NM, USA \move@AU\move@AF\@affiliationDAWN Center for Reionization, Niels Bohr Institute, University of Copenhagen, Copenhagen, Denmark \move@AU\move@AF\@affiliationUniversity of Melbourne, Melbourne, Australia \move@AU\move@AF\@affiliationARC Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D) \move@AU\move@AF\@affiliationUniversity of Washington, Seattle, WA, USA 0000-0002-2620-7056]Claudio Dalla Vecchia \move@AU\move@AF\@affiliationInstituto de Astrofísica de Canarias, La Laguna, Tenerife, Spain \move@AU\move@AF\@affiliationDepartamento de Astrofísica, Universidad de La Laguna, La Laguna Tenerife, Spain \move@AU\move@AF\@affiliationInstitute for Astronomy, University of Edinburgh, Royal Observatory, Edinburgh, UK

1 Introduction

The reionization of the intergalactic medium (IGM) was the last major phase change in the universe, when high energy ultraviolet (UV) photons from the first luminous sources in the universe ionized hydrogen (and singly ionized helium) in the IGM. Observational constraints on this epoch come from a variety of complementary techniques, and are continuously improving in accuracy and growing in number. Present-day observations constrain the bulk of reionization to be completed by 6 (e.g., Fan et al., 2006; McGreer et al., 2015), though some lines of sight may remain somewhat neutral to 5.5 (e.g., McGreer et al., 2015; Kulkarni et al., 2018; Pentericci et al., 2018). The beginning of reionization is less well constrained, and depends sensitively on the nature of the ionizing sources. If rare objects such as quasars provide the bulk of the ionizing photons, reionization likely didn’t get well underway until 10 (e.g, Madau & Haardt, 2015). On the other hand, if young, massive stars dominated the ionizing photon budget, reionization may have started much sooner, although the constraints on the electron-scattering optical depth to the cosmic microwave background (CMB) measured by Planck Collaboration et al. (2016b) imply that the halfway point came at 8. The apparent dichotomy between the sharp decline in the number density of bright quasars at 2 (e.g., Richards et al., 2006; Hopkins et al., 2007) and the relatively shallower decline in the UV luminosity density from galaxies (e.g., Madau & Dickinson, 2014, and references therein) has led to the predominant theory that the bulk of the ionizing photon budget came from massive stars.

Better understanding both the temporal and spatial evolution of the process of reionization is key to understanding a variety of unknown physical processes in the early universe, including the time of the onset of the first stars and galaxies, the effects of reionization heating on galaxy formation and growth, and the escape of ionizing photons from galaxies. Present-day efforts to reconstruct the progress of hydrogen reionization involves several major uncertainties around the contribution of both massive stars in galaxies and quasars. Over the past decade advances in the capabilities of near-infrared imaging on the Hubble Space Telescope, necessary to measure rest-frame UV light in the epoch of reionization, have led to robust constraints on the observable non-ionizing UV (1500 Å) luminosity density from galaxies in this epoch.

To understand how these galaxies contribute to reionization, one needs convert this to an ionizing emissivity (; the number of ionizing photons produced per unit time per unit volume which escape the galaxy) as a function of redshift (e.g., Finkelstein et al., 2012a, 2015b; Robertson et al., 2013, 2015; Bouwens et al., 2016b, 2015b), which is dependent on three factors: the rest-UV non-ionizing specific luminosity density (), the ionizing photon production efficiency (), and the escape fraction of ionizing photons (). The product of the first two quantities produces the intrinsic ionizing emissivity produced within galaxies (), which when multiplied by produces the escaping ionizing emissivity . This quantity can be used to infer the evolution of the IGM ionized volume filling fraction (denoted as ) by solving a set of ordinary differential equations, which depend on this emissivity, the density of hydrogen, and the recombination time (dependent itself on the clumping factor of the gas and the temperature-dependent recombination coefficient; e.g., Madau et al. 1999; Robertson et al. 2013).

The value of is measured by integrating the observed rest-UV luminosity function to some observationally unknown limiting magnitude. This limiting magnitude is crucial as the steepening faint-end slope with increasing redshift means that the faintest galaxies dominate (Bunker et al., 2004; Yan & Windhorst, 2004; Bouwens et al., 2015b; Finkelstein et al., 2015b, e.g.,). Exactly how dominant these faint sources are depends on the shape of the extreme faint-end, which should reverse its steep rise due to stellar feedback, the ability of halos to atomically cool, and Jeans filtering due to the reionization-driven UV background, as shown by a variety of simulations (e.g. Bullock et al., 2000; Gnedin, 2000; Iliev et al., 2007; Okamoto et al., 2008; Mesinger & Dijkstra, 2008; Finlator et al., 2011, 2012; Alvarez et al., 2012; Oñorbe et al., 2017; Jaacks et al., 2018b). Informed by this theoretical work, observational studies have commonly used 13 as this integration limit (e.g., Robertson et al., 2015; Finkelstein et al., 2015b).

This limit is 100 fainter than that achievable in even the Hubble Ultra Deep Field (Beckwith et al., 2006) at these redshifts. However, recent observations of much fainter galaxies rendered detectable via gravitational lensing in the Hubble Frontier Fields (Lotz et al., 2017) have begun to provide empirical justification, with evidence that the observed luminosity functions maintain their steep slopes down to 16 at 6 (Atek et al., 2015), and possibly to 15 (Livermore et al., 2017; Bouwens et al., 2017; Atek et al., 2018). We note the concept of a limiting magnitude is an approximation, as the luminosity function should gradually roll over rather then exhibit a steep cut-off (e.g. Jaacks et al., 2013; Weisz et al., 2014; Boylan-Kolchin et al., 2015; Jaacks et al., 2018b), and any cut-off or turnover point will evolve with redshift as the halo masses evolve, the UV background ramps up, and feedback effects manifest.

The ionizing photon production efficiency converts the (dust-corrected) rest-UV non-ionizing specific luminosity density [erg s Hz Mpc] to the intrinsic ionizing emissivity [s Mpc]. This efficiency depends on the surface temperatures of the massive stars, which in turn depends on the stellar metallicity, age, and binarity, as well as the initial mass function (e.g., Eldridge & Stanway, 2009; Stanway et al., 2016; Stanway & Eldridge, 2018). If these quantities were known, one could then measure directly. There are however large uncertainties, thus until recently most studies assumed a value of 25.2, expected from modestly metal-poor, but otherwise-normal, single-star population models (e.g. Finkelstein et al., 2012a; Robertson et al., 2015). This is consistent with the observations that stellar populations in faint 7 galaxies are non-primordial (e.g., Finkelstein et al., 2010, 2012b; Wilkins et al., 2011; Bouwens et al., 2012, 2014; Dunlop et al., 2013), though it is possible that fainter galaxies have much lower metallicities (Dunlop et al., 2013; Jaacks et al., 2018a).

Recent work has shown that the typically assumed conversion from observed non-ionizing to ionizing UV is consistent with the inferred strength of H emission deduced from IRAC photometric colors for bright galaxies at 4 (e.g., Bouwens et al., 2016a). However, this same work shows evidence that fainter/bluer galaxies, and galaxies at 5 have higher values of . This implies that may vary both with galaxy luminosity (and/or perhaps halo mass) as well as redshift, consistent with the high values of inferred from the few 7 galaxies with detectable C iii] emission (Stark et al., 2016).

Finally, one needs to assume an escape fraction () for ionizing photons, which is the dominant source of uncertainty. A variety of analyses have shown that when assuming a limiting magnitude of 13 and 25.2, an escape fraction of 10–20% produces the requisite number of ionizing photons to complete reionization by 6 with no contribution from other sources (e.g., Finkelstein et al., 2012a; Robertson et al., 2013; Finkelstein et al., 2015b; Robertson et al., 2015; Bouwens et al., 2015a). The assumption of a relatively high escape fraction at 6 is impossible to directly verify, as even a predominantly ionized IGM produces an ionizing optical depth sufficient to absorb all ionizing UV radiation at 4 (e.g., Vanzella et al., 2018).

We must observe galaxies at to directly measure , where there is unambiguous observational evidence that most studied galaxies have low escape fractions (e.g., Siana et al., 2010; Sandberg et al., 2015; Rutkowski et al., 2017; Grazian et al., 2017, though see Steidel et al. 2018). Recent observational programs have improved at identifying galaxies likely to exhibit higher escape fractions, specifically those which exhibit intense ionizing environments as traced by ratios of nebular emission lines, resulting a few dozen direct detections of escaping ionizing photons (e.g., Shapley et al., 2016; de Barros et al., 2016; Bian et al., 2017; Vanzella et al., 2018; Izotov et al., 2018). However, the lack of significant ionizing photon escape from the bulk of galaxies strongly implies that the escape fraction from all galaxies at all redshifts cannot be as high as 10–20%.

One way to reconcile this, suggested by a variety of simulations, is if the escape fraction is dependent on the halo mass, where lower-mass halos have higher escape fractions due to lower gas covering fractions and an increased susceptibility of starburst-driven escape routes (e.g., Paardekooper et al., 2013; Wise et al., 2014; Paardekooper et al., 2015; Anderson et al., 2017; Xu et al., 2016), while massive halos occasionally exhibit high escape fractions for short periods due to extreme starburst-driven winds clearing channels in the ISM (e.g. Paardekooper et al., 2015).

\H@refstepcounter

figure \hyper@makecurrentfigure

Figure 0. \Hy@raisedlink\hyper@@anchor\@currentHrefLeft) The reference luminosity functions from Finkelstein (2016) used for this work, with the light-to-dark blue shading denoting 4 to 14. The dashed lines show the case when we fix the faint-end slope at 10 to be equal to the value at 10. Center) The relation between halo mass and UV magnitude obtained by abundance matching the luminosity functions from the left panel following Behroozi et al. (2013). The dark gray region denotes the halo mass range where H i cooling likely does not take place. The lighter gray region denotes the regime when the post-reionization UV background likely suppresses star formation. Right) The non-ionizing UV luminosity density, highlighting the much shallower evolution when the faint-end slope is allowed to evolve to extremely steep values at 10.

In this work, we make use of halo-mass dependent escape fractions predicted from simulations to explore scenarios for completing reionization with low escape fractions for most observable galaxies. In §2, we focus on a critical examination of all assumptions needed to solve for the reionization history, while in §3 we discuss our MCMC framework which we use to probe the full parameter space for all such assumptions, using predominantly model-independant reionization observations to constrain our analysis. These results are given in §4, and discussed in §5. In §6 we explore the implications on the cosmic star-formation rate density to faint luminosities and extremely high redshifts. Throughout this paper, we assume AB magnitudes (Oke & Gunn, 1983) and a Planck 2015 cosmology: H=67.74 km s Mpc, 0.309, 0.691, 0.0486 and Y 0.2453 (Planck Collaboration et al., 2016a). We use the variable to denote both halo mass and absolute magnitude, thus we distinguish these by and , respectively.

2 Defining the Total Ionizing Emissivity

To model reionization we must calculate the evolution of the ionizing emissivity with redshift, which depends on a number of variables discussed above. Here we attempt to gain new insight into reionization by pairing galaxy observations with simulated, halo-mass dependent escape fractions. As shown in the Appendix, simply replacing a flat 10–20% escape fraction with results from simulations is destined to fail due to the low escape fraction values for all but the smallest halos. We must thus re-examine the assumptions for all critical variables. We allow this by flexibly exploring the dependence of reionization on assumptions about the ionizing photon production efficiency, limiting halo mass for star formation, and a potential contribution from AGN to the ionizing emissivity.

As described in this section, our model includes seven free parameters which, when combined with the observed UV luminosity function, define the emissivity as a function of redshift. In §3 we describe how we constrain the posterior distribution of these parameters within a MCMC framework constrained by several robust observations. We restrict our model to 4, as this more than encompasses the full epoch of hydrogen reionization, and at lower redshifts dusty star-forming galaxies (some of which could be absent from our UV luminosity functions) may contribute to the ionizing photon budget (e.g., Gruppioni et al. 2013; Cowie et al. 2017; Koprowski et al. 2017, though see Casey et al. 2018).

2.1 The Galaxy Ionizing Emissivity

2.1.1 Luminosity Functions

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 a Markov Chain Monte Carlo (MCMC) Schechter function fit to all recently published data at 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(), faint-end slope () and characteristic number density (). 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 at 7, finding 2 (e.g., Livermore et al., 2017; Bouwens 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 4 to 20. As the data used to derive these luminosity functions were limited to 10, it is unknown if this extrapolation is valid. Specifically, the faint-end slope is observed to evolve somewhat steeply with redshift (0.11), implying 2.90. While this may be possible (indeed the model results from Mason et al. 2015 find 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 10, and remains at the 10 value of 2.35, while in § 4.2 we also explore the case where continues to steepen at 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.

\H@refstepcounter table \hyper@makecurrenttable Table 0. \Hy@raisedlink\hyper@@anchor\@currentHrefUV Luminosity Function Parameters

Redshift M log (mag) (Mpc) 4 21.05 1.69 2.99 6 20.79 1.91 3.37 8 20.52 2.13 3.75 10 20.25 2.35 4.13 12 19.98 2.57 4.50 14 19.71 2.79 4.88 16 19.44 3.01 5.25 18 19.17 3.23 5.63 20 18.90 3.45 6.00

Note—The assumed rest-frame UV luminosity functions used in this work, following the evolutionary trend derived via observations as discussed in Finkelstein (2016). Our fiducial model keeps the faint-end slope fixed at 10 to the 10 value of 2.35.

2.1.2 Abundance Matching

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. (2015a), we assume a log-normal UV magnitude scatter at fixed halo mass of 0.2 dex, though we note that a scatter as high as 0.4 dex does not affect the relation at log (/M) 11. Our derived relations are shown in the middle panel of Figure 1 for both luminosity function evolution cases we consider.

2.1.3 Minimum Halo Mass for Star Formation

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 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 6 luminosity function has a very steep faint-end slope (e.g., Bouwens et al., 2015b; Finkelstein et al., 2015b; 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 which can both retain its gas, and cool it\go 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) 8 at 6 (Okamoto et al., 2008; Finlator et al., 2012), and this critical mass shifts to lower values at higher redshifts, as halos which 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 which find that significant star formation is happening down to log (/M) 7 at 10–15 due to the availability of metal-line cooling in the immediate aftermath of the formation of Pop 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 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 photo-heated, 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 mini-halos with log (/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) 9, though the predictions are quite uncertain\go (e.g., Gnedin, 2000; Iliev et al., 2007; Mesinger & Dijkstra, 2008; Okamoto et al., 2008; Alvarez et al., 2012). Feedback likely also plays a strong role (see Somerville & 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) find that stellar feedback causes a flattening in the UV luminosity function at 14, or log (/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, , 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 photo-suppression below a threshold halo mass due to the rising UV background. However, there are a range of plausible limiting halo masses for this photo-suppression 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 photo-suppression mass as a free parameter, , with an adopted flat prior of log () (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 dependance was needed, and to avoid adding another free parameter to our model.\go The non-ionizing specific UV luminosity density () 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-hand panel of Figure 1. In § 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 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 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) show in their Figure 17, although their UV luminosity function continues to 8, the UV luminosity density asymptotes to a constant value when integrating to 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 extreme low-mass halos could have high ionizing photon escape fractions.

2.1.4 Ionizing Photon Production Efficiency

To convert the total non-ionizing ultraviolet luminosity density to the ionizing emissivity, a value for the ionizing photon production efficiency () needs to be assumed. This parameter encompasses all of the physics of the underlying stellar population, many 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; Finkelstein et al., 2012b; Bouwens et al., 2014), leading to hotter stellar photospheres as the metal opacity (mostly due to iron) is lower, and thus a higher ionizing to non-ionizing UV photon ratio. Another factor to be considered is the effect of binary stars.\go Stellar population synthesis models which include binary stars (Eldridge & Stanway, 2009) show that the ionizing flux is boosted by 60% (at low metallicities of 0.3) 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 their via the detection of nebular emission lines. Bouwens et al. (2016a) inferred H emission line fluxes from Spitzer/IRAC photometry at 3.8 5.0, finding log 25.34 erg Hz, consistent with typically assumed values of log 25.2–25.3 erg Hz in previous reionization studies (e.g., Madau et al., 1999; Kuhlen & Faucher-Giguère, 2012; Finkelstein et al., 2012a; Robertson et al., 2015). At 5.1 5.4, Bouwens et al. (2016a) found log 25.48 Hz erg, hinting at evolution towards larger values at higher redshift, though not at a significant level given the observational uncertainties. Bouwens et al. (2016a) also find evidence that the bluest galaxies exhibit even higher values of , with log 25.9 erg Hz for galaxies in their 5.1 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. (2015a) and Stark et al. (2016) measured via ionized carbon emission, finding 25.68 erg Hz in a lensed galaxy at 7.045 with an intrinsic 19.3, and 25.6 for three luminous (22) galaxies at 7.15, 7.48 and 7.73.\go Lastly, Wilkins et al. (2016b) investigated the range of expected from galaxies in the epoch of reionization based on the BlueTides simulation, finding that simulated galaxies spanned the range 25 26, with the highest values obtained when assuming low-metallicity stellar population models which include binaries.

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

(1)

where is the reference magnitude of 20. We assume flat priors on both free parameters of log (0.0,0.4) and log (0.0,0.2). was also 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., 2016a; Wilkins et al., 2016b; Izotov et al., 2017).

2.1.5 Escape Fractions

To derive the ionizing emissivity available to ionize the IGM, we must combine our intrinsic ionizing emissivity () with a model for , 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; Robertson et al., 2013, 2015; Bouwens et al., 2015b; Finkelstein et al., 2015a). However, as discussed in §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 6 (Khochfar et al in prep). The escape fraction of ionizing photons was determined in more than 75,000 haloes 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\go 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 anti-correlated 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 which is independant of redshift (we will address photoionization feedback below). At each halo mass (in steps of log[/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/20th part of the data, and the log-likelihood of the remaining 1/20th 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.\go 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.1.5. 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\go. This figure highlights that the escape fraction distributions are quite broad, but only halos with log(/M) 8.5 have more than half of their distribution at 1%.

\H@refstepcounter

figure \hyper@makecurrentfigure

Figure 0. \Hy@raisedlink\hyper@@anchor\@currentHrefThe probability distribution functions of the ionizing photon escape fraction for different halo masses. These distributions come from the simulations of Paardekooper et al. (2015), based on high-resolution radiative transfer modeling of 75,000 halos extracted from the First Billion Years (FiBY) simulation (Khochfar et al., in prep). While the escape fraction does not appear to be heavily redshift dependent in their explored epoch of 6 15, as shown here it varies quite strongly with halo mass, with only halo masses with log (/M) 8.0 having a majority of their probability distribution at 1%. Not shown in this figure is a small peak in the distribution at log 10 for log (/M) 8.5, comprising 10% of the probability density. The thin vertical lines denote the median value of each distribution, ranging from 16% for log (/M) 7, to 0.1% for log (/M) 8.5.

At log(/M) 9 there is a small probability that 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) 9.5, we assume that halos with log(/M) 11 have 0, and that halos with log(/M) 9-11 have a similar distribution as those at halos with log(/M) 9 but without the small peak at high (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.1.5. We note that if we had treated halos with log(/M) 11 the same as those with log(/M) 10, we find almost no differences in the resulting ionization history (completing at 5.7, compared to 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 (§2.2).\go

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\go. Simulations have shown that better resolving the 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 15, and the expected physical reasons for this scale factor do not depend on redshift. We denote this parameter below as , and adopt a flat prior on over the range (0,10).

The total ionizing emissivity from galaxies is thus calculated as , where the latter term includes this scale factor.

2.2 Inclusion of an Active Galactic Nuclei Contribution

While quasars have been disfavored as the dominant source of the reionization ionizing photon budget (e.g., Shapiro & Giroux, 1987)\go, the low observed galaxy escape fractions leave room for some contribution from active galactic nuclei. This is not in violation of previous results, as most observations at 4 probe the bright end of the AGN luminosity function (e.g., quasars only), thus, similar to galaxies, it may be that the AGN luminosity function has a steepening faint end slope, and that faint AGNs, and not the rare quasars, are significant contributors.

There is observational evidence in favor of this possibility, as Giallongo et al. (2015) discovered faint AGNs at 4–6 by searching the positions of known galaxies at those epochs in deep Chandra X-ray data in the GOODS-S field. At 4, Giallongo et al. found ionizing emissivities nearly an order of magnitude greater than those implied by the bolometric quasar luminosity function work of Hopkins et al. (2007), and a factor of a few higher than Glikman et al. (2011), with the evolution to 5–6 shallower than that of Hopkins et al. (2007). Taken at face value, ionizing photons generated from AGN could account for the entire reionization photon budget, with no contribution from galaxies at all (Madau & Haardt, 2015). The Giallongo et al. results have been met with some skepticism over the photometric redshifts of the sources (e.g., Parsa et al., 2018), and also the apparent lower emissivity at similar redshifts (e.g., McGreer et al., 2018). Additionally, Giallongo et al. note that they cannot rule out a significant stellar contribution to the X-ray luminosity. However, given the difficulties in isolating faint AGN at high redshift (e.g., Stevans et al., 2018), and the fact that current observations span a large range at 4, we allow a contribution from AGN to the ionizing budget in our fiducial model. Figure 2.2 shows the inferred evolution of the AGN comoving ionizing emissivity both from the previous work by Hopkins et al. (2007), and the “quasars can do it all” recent work by Madau & Haardt (2015), along with a number of observational results from the literature.

\H@refstepcounter

figure \hyper@makecurrentfigure

Figure 0. \Hy@raisedlink\hyper@@anchor\@currentHrefThe evolution of the AGN comoving monochromatic 912 Å emissivity with redshift. The dashed green line shows the results from the Hopkins et al. (2007) bolometric quasar luminosity function, while the dashed purple line shows the form proposed by Madau & Haardt (2015), which allows quasars to complete reionization with no contribution from star-forming galaxies. The circles denote results from the literature, using the compilation provided by Madau & Haardt (2015), with orange symbols denoting the recent results from McGreer et al. (2018) and Akiyama et al. (2018), and the orange bar denoting the range from Stevans et al. (2018). The 68% confidence range on our fiducial result is shown as the blue shaded region (with the shading density denoting the shape of the probability distribution function), which is consistent with the observed data, and roughly between the two previous evolutionary trends, at 4.

Our initial emissivity matches Madau & Haardt (2015) at 2.5, and at higher redshifts is a simple exponential with a slope constrained to be between those of Hopkins et al. (2007) at the low end, and Madau & Haardt (2015) at the high end, spanning the full range of observational results. Our emissivity is governed by three free parameters: a scale factor (0,1) applied to the emissivity allowing it to be lower than initially assumed (due to a range of physical effects, including a non-unity AGN ionizing photon escape fraction, which is likely the case for less-luminous AGNs; Trebitsch et al. 2018); a redshift-evolution exponential slope (1.05,0.34), which approximately reproduces the Hopkins et al. 2007 and Madau & Haardt 2015 respective AGN ionizing emissivity evolution in this formalism; and a maximum redshift , above which the AGN ionizing emissivity is assumed to be zero. The functional form for our monochromatic 912 Å emissivity is given by

(2)

where the first exponential term is the initial emissivity, and the term in parentheses is a normalization factor, normalizing our emissivity (prior to the application of a scale factor) to be equal to that of Madau & Haardt (2015) at 2.5 (whose functional form is given in the numerator). We note that while this emissivity is included in our analysis, it is allowed within our formalism to be negligibly low in the epoch of reionization, thus we are not “forcing” AGNs to contribute significantly. We discuss our fiducial results in §4, but they are shown in Figure 2.2, falling roughly in the middle of the allowed range.

The total emissivity from our model at a given redshift is the sum of that from galaxies (§ 2.1) and AGN: .

\H@refstepcounter table \hyper@makecurrenttable\hb@xt@ Table 0. \Hy@raisedlink\hyper@@anchor\@currentHrefMCMC Model Parameters

Parameter Name Initialization Initialization Flat prior Posterior Central Value constraints Median (68% Confidence) f 5.0 1.0  0, 10 5.2 (3.3 to 7.5) log (M/M) 9.0 0.5  8.5, 10.5 8.90 (9.5) log 0.10 0.05  0, 0.4 0.13 (0.05 to 0.25) log 0.05 0.03  0, 0.2 0.07 (0.03 to 0.13) AGN 0.8 0.2  0, 1 0.77 (0.47) 10.0 2.0  4, 12 9.20 (6.9) AGN 0.5 0.3  1.2, 0.1 0.39 (0.93)

Note—The free parameters for our fiducial model. The initialization central value and initialization define a normal distribution from which each walker draws an initial value. The scale factor applied to the halo-mass-dependent escape fractions from the Paardekooper et al. (2015) simulations. Post-reionization photosuppression halo mass. Evolution of ionizing photon production efficiency with redshift and absolute magnitude. Scale factor applied to the AGN emissivity (mimicking an AGN ionizing photon escape fraction). Exponential slope of the AGN emissivity with redshift, constrained to be zero at some maximum redshift. The last column gives the median of the posterior distribution, and the central 68% confidence range (or upper/lower 84% confidence limits when the distribution is one-sided).

\twocolumngrid

2.3 Calculating

We calculate the IGM volume ionized fraction by solving the differential equation

(3)

where is the comoving ionizing emissivity derived above, is the comoving hydrogen density, and is the IGM hydrogen recombination time. The comoving hydrogen density is calculated as the product of the hydrogen mass fraction (defined as , where is the helium mass fraction), the dimensionless cosmic baryon density , and the critical density (defined as 3/8). The IGM recombination time is given by

(4)

where is the temperature-dependent case B recombination coefficient for hydrogen using the functional form given by Hui & Gnedin (1997). Following Robertson et al. (2015, hereafter R15), we evaluate this at 20,000 K (had we assumed 15,000 K, would be higher by a factor of 1.29). We assume a redshift-dependent hydrogen clumping factor from the simulations of Pawlik et al. (2015), which evolves from 4.8 at 6, to =1.5 at 14. We solve for by integrating the ordinary differential equation in Equation 1 using the IDL routine from 20 to 4.

3 Exploring the Full Reionization Parameter Space with MCMC

Using the set of seven free parameters defined in §2 our model can describe the escaping ionizing emissivity from both star-forming galaxies and super-massive black hole accretion (AGN) activity. In this section we describe how we use a MCMC framework to derive the posteriors on these free parameters using a set of robust observational constraints. We used an IDL implementation of the affine-invariant sampler (Goodman & Weare, 2010) to sample the posterior, which is similar in production to the Python emcee package (Foreman-Mackey et al., 2013). We used the recommended stretch parameter of with 1000 walkers. Each walker was initialized by choosing a starting position for each of the free parameters, randomly drawn from a normal distribution with a central value and width given in Table 2. We assumed a flat prior on each of our seven free parameters, with the prior bounds also listed in Table 2. If the log likelihood of a given set of parameters was not finite (i.e., it violated the parameter flat priors), a new set of parameters was drawn, until a set which gave a finite probability was drawn to initialize each of the 1000 walkers. The exact initialization values are not crucial as the burn-in process ensures that the starting positions do not affect the results.

3.1 Method

In this sub-section we describe in detail our MCMC analysis. A flowchart of this procedure is shown in Figure 1. In each step of the chain, our routine used the chosen set of seven free parameters to complete the following calculations:

  1. For each redshift interval of =0.1 from 4 to 20, we use a randomly drawn set of Schechter function parameters from the available F16 MCMC chains, where the drawn parameters are the redshift-evolution terms , , and . We use these parameterizations to calculate the non-ionizing specific UV luminosity density in absolute magnitude bin intervals of =0.1 from 6 to 24. These were then corrected for dust attenuation () using the method described in Finkelstein (2016), which uses the relation between and from Bouwens et al. (2014), the relation between and A from Meurer et al. (1999), and the dust attenuation curve from Calzetti et al. (2000). A scatter in at fixed of 0.35 was assumed, and zero dust attenuation was assumed at 9 (Bouwens et al., 2014; Wilkins et al., 2016a). We note that given the halo mass dependency of the escape fractions the bulk of ionizing photons come from faint galaxies with minimal dust, thus our final results are not sensitive to this correction. To validate this, we performed a model run with no dust correction, and found no significant change in the evolution of the ionization history\go.

    \H@refstepcounter

    figure \hyper@makecurrentfigure

    Figure 0. \Hy@raisedlink\hyper@@anchor\@currentHrefA visual description of our Markov Chain Monte Carlo procedure for constraining the posteriors on our free parameters, described in full in § 3.1. All figures appear full-size elsewhere in the paper.

  2. The intrinsic ionizing emissivity was calculated by multiplying by the appropriate value of for the values of log and log in a given step. The escaping ionizing emissivity was then calculated as multiplied by the escape fraction, where the escape fraction is randomly drawn in each step of the chain for each absolute magnitude interval, from the probability distribution function (PDF) corresponding to the halo mass for the given absolute magnitude (from the relations described above). One feature of our process is that by randomly sampling these PDFs over many MCMC chain steps, we marginalize over the distribution of possible escape fractions, such that this scatter is encompassed in our final results.

  3. The IGM volume ionized fraction of hydrogen was calculated following § 2.3. While solving the differential equation, we emulated the effects of photosuppression by calculating the emissivity down to the limiting UV magnitude corresponding to both the atomic cooling limit at a given redshift (; applicable for neutral regions), and also to for the given chain step (; for ionized regions). The total value of for each redshift bin was then calculated as

    (5)

    where the first term accounts for ionizing photons from all galaxies above the photosuppression limit, while the second term accounts for those photons from halos between the photosuppression limit and the atomic cooling limit, but only in the fraction of the volume which is still neutral at a given redshift. We note that this is an approximation as we can only track the globally-averaged ionized fraction, and thus it does not account for the effects of spatial clustering of halos on the ionized fraction in their proximity (e.g., the topology of reionization). The total ionizing emissivity was then calculated as that from galaxies combined with that from AGN, , where the latter was calculated from as described in § 2.2, assuming an AGN H I ionizing spectral index of 1.7 (Lusso et al., 2015).

  4. While helium becomes singly ionized at a similar energy as hydrogen, high-energy photons from AGN can doubly-ionize helium. We thus calculate the emissivity of He ii ionizing photons (energies 4 Ryd) again assuming a spectral index of 1.7, and solve for the IGM volume ionized fraction of He iii using

    (6)

    and

    (7)

    The volume ionized fraction of He ii was then assumed to be . We assume (see § 5.3.2 for discussion).

  5. Using the calculated volume ionized fractions for H ii, He ii, and He iii, we calculated the electron scattering optical depth as measured from the CMB as

    (8)

    with

    (9)

    integrated from 0 to 20.

3.2 Observational Constraints

The outcomes of these calculations are the volume ionized fractions of H i, He i and He ii, the galaxy and AGN ionizing emissivities, and the electron scattering optical depth, all as a function of redshift. To calculate the likelihood for our model and constrain our free parameters, we used the following observations:

  1. The integrated hydrogen ionizing emissivity at 4.0 and 4.75 from Becker & Bolton (2013), synthesized from a variety of measurements of the IGM based on spectroscopy of high-redshift quasars. They find that this quantity rises from 3.2 to 4.75, which is consistent with the idea that a steepening galaxy UV luminosity function faint-end slope results in more ionizing photons at higher redshifts. We use their measurements of log 0.139 photons s Mpc at 4.0 and log 0.014 photons s Mpc at 4.75. These uncertainties include both the statistical errors, and the much larger systematic errors. We did not use measurements at 4, as our luminosity functions were calculated only at 4. We also did not include the upper limit on the integrated emissivity at 6 from Bolton & Haehnelt (2007), as the value of log (0.585 was derived assuming that the ionizing background is uniform, while recent measurements imply that there are substantial spatial variations (e.g., Fan et al., 2006; Becker et al., 2015; Bosman et al., 2018; Becker et al., 2018). Future results taking advantage of updated measurements of the spatial inhomogeneities in the ionizing background, the IGM temperature, and the mean free path of ionizing photons will both decrease these systematic uncertainties, and allow more robust results at higher redshifts, further constraining models such as the one we present here. For each step in our chain, we calculated the goodness-of-fit statistic between both the 4 and 4.75 observations and the summed galaxy and AGN ionizing emissivity from our model.

    \H@refstepcounter

    figure \hyper@makecurrentfigure

    Figure 0. \Hy@raisedlink\hyper@@anchor\@currentHrefEach panel shows the posterior distribution for each of our modeled free parameters. The different colored lines denote which steps out of our 100,000-step chain were used (denoted by the legend in the lower-right). As these distributions do not appear to significantly change after Step #70,000, we adopt as our fiducial posteriors that from the last 5,000 steps of the chain (corresponding to the last 5%), denoted by the black line. We performed the Gelman-Rubin test, which showed that this model is highly converged. The majority of these parameters have a relatively broad posterior distribution, which is propagated forward into the uncertainties in our reionization model. More robust observational constraints are needed if we wish to further constrain these parameters.

  2. The electron scattering optical depth to the CMB (). This quantity measures the integrated optical depth of Thomson scattering of CMB photons with free electrons as they travel from the surface of last scattering to the present day, and thus it possesses constraining power on models of reionization. The first measurements of from the Wilkinson Microwave Anisotropy Probe (WMAP) Year 1 results showed 0.170.06, which suggested an instantaneous “reionization redshift” of 175 (Spergel et al., 2003). Additional data from WMAP revised these estimates immediately downward, from the WMAP Year 3 result of 0.0880.03 (Spergel et al., 2007) to the final WMAP Year 9 result of 0.0880.013 (Hinshaw et al., 2013), and 10.51.1. The advent of the Planck satellite has revised these estimates again downward to 0.0660.012 and 8.81.1 (Planck Collaboration et al., 2016a). Given the relatively large uncertainties of both the WMAP 9 and Planck measures, the discrepancy is only significant at the 1.3 level. However, even more recent 2016 results have been published, highlighting improved removal of systematics from the Planck high-frequency data, showing 0.0550.009 (Planck Collaboration et al., 2016b), discrepant from the WMAP9 data at 2.1 significance.

    We elect use this newer 2016 Planck value of as our fiducial constraint. When comparing to R15, it is important to remember that they used the 2015 value, though we note that the 2015 and 2016 Planck values differ only at the 0.7 level. For each step in our chain, we computed between the observational value of and that calculated from our model given in Equation 8, which includes the contribution to both from ionizing hydrogen, and singly and doubly ionized helium.

  3. The model-independent lower limits on the IGM ionized fraction of Q(1) at 5.9, and Q (1) at 5.6 from McGreer et al. (2015). This study measured the fraction of “dark” pixels in the Ly and Ly forests of a sample of 22 bright quasars at 5.7 – 6.4. Regions of the IGM containing any pre-reionization neutral hydrogen should result in completely saturated absorption in both of these transitions. This method only provides a lower limit, as some absorption may also be due to collapsed systems, or residual H i in ionized gas. However, it is model independent as it does not depend on the intrinsic quasar spectral shape (see discussion in McGreer et al., 2015). Finally, as this method combines several objects (and, several locations along the line-of-sight to each object), it is far more robust against cosmic variance, and uncertainties due to inhomogeneous reionization than neutral fractions derived via effective optical depths to single quasars (e.g., Fan et al., 2006; Bolton et al., 2011; Greig et al., 2017; Bañados et al., 2018).

    For each step in our chain, we calculated the goodness-of-fit statistic between the ionized fraction in our model at each of these two redshifts, and these measurements. As these are one-sided lower limits on the ionized fraction, if the model value was above these measurements ( 0.94, 0.96 at 5.9, 5.6), was set to zero; if the model was below these values, was calculated in the usual way using the published uncertainties in each redshift bin.

3.3 Deriving Posteriors

Rather than choose a pre-defined number of steps for the burn-in period, we elected to run our chain for 10 steps, and then examine the results to explore where the chain has converged, and select a final set of samples to derive the posteriors. In Figure 1 we show the distributions of each of our seven free parameters for different groupings of 10 steps. One can see that over the first few iterations of 10 steps, the parameter distributions change, but that after 7 10 steps the changes begin to stabilize, such that the distributions only exhibit minor changes towards the end of the chain. For this reason, we defined the last 5000 steps of the chain as those used to sample the posterior distribution of our free parameters.

\H@refstepcounter

figure \hyper@makecurrentfigure

Figure 0. \Hy@raisedlink\hyper@@anchor\@currentHrefThe covariances between our model parameters from the posterior distribution of our fiducial model. The purple, yellow and blue contours denote the 68%, 95% and 99.5% confidence range between the listed model parameters, while the histograms show the distribution of a single parameter (similar to Figure 1).

The acceptance fraction for this fiducial model was 14.5%. We also ran this model ten independant times to check for convergence, finding that the acceptances fractions spanned 14.1% – 14.6%. We tested for convergence by comparing the posterior distributions of our free parameters using the Gelman-Rubin test. Specifically, we used the IDL routine, which computes the Gelman-Rubin statistic. We found 1.00 for all seven of our free parameters, showing that this model is highly converged.\go

When we compare versions of our model, we use the Bayesian Information Criteria (rather than the full Bayes factors, which are prohibitive to compute for our high-dimensional parameter space). This is defined as , where is the median likelihood from a given posterior, N is the number of data points, and k is the number of free parameters (Liddle, 2004). This is preferable over the simpler statistic, as it takes into account both the number of data points and the number of free parameters. For a model to be preferred over a competing model, it must have a BIC lower by at least 2.\go

4 Results

4.1 Fiducial Model

4.1.1 Posterior Distributions of Free Parameters

The black lines in Figure 1 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 3.3. Of these seven distributions, three have a clearly defined peak in the posterior distribution within our prior range. The parameter f prefers a value near 5, with a 68% confidence range of 3.3 f 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 § 2.1.5, this up-scaling is not surprising, as the simulations were unable to resolve the birth cloud of the star particles, and so the porosity of the gas may be underestimated (though see also Gnedin & Kaurov 2014, who invoke a sub-unity scale factor to account for unresolved systems providing excess absorption)\go. In § 5.2 we will discuss the implications of this up-scaling for the average galaxy escape fraction, and its evolution with redshift.

The parameters related to also show a clear peak, with log 0.13, and log 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, but also observations from local analogs for high-redshift galaxies, which we discuss in § 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 1-sided 84% upper limit of log (/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) 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 §5.4, lower values would be in greater conflict with previous simulation work on the post-UV background\go filtering mass.

The final three parameters govern the AGN contribution. Combined, they prefer an AGN contribution which 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 §5.3.1, we compare the resultant emissivity from these combined parameters to previous observations. The likelihood from the median value of these parameters corresponds to a BIC value of 16.9\go.

4.1.2 Comparison with Constraints

Figure 4.1.2 shows the distribution of both and for our fiducial model. The blue-shaded region shows the 68% confidence range on our evolution of . 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 5.6 and 5.9. Our model achieves 1 by 5.6 0.5. Although previous studies have adopted a prior that reionization ended by , 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 – in particular, the comoving Mpc Ly trough observed by Becker et al. 2015 (Kulkarni et al., 2018).

\H@refstepcounter

figure \hyper@makecurrentfigure

Figure 0. \Hy@raisedlink\hyper@@anchor\@currentHrefThe comparison of the ionized volume filling fraction of hydrogen (in blue) as a function of redshift to observations. The blue shading denotes the shape of the probability distribution function within the 68% central confidence range, with the darkest shading denoting the median (similar shading is used in many of the remaining figures). The observations from McGreer et al. (2015) used to constrain our model are shown in green, while in gray we show the results from the model of Robertson et al. (2015). Our model completes reionization by 5.6 0.5. Although it was not used to constrain our model, we also show the ionized volume filling fraction of He ii (in red), noting that although AGNs are included in our model, He ii does not reionize too early (see § 4.1.3).

\H@refstepcounter

figure \hyper@makecurrentfigure

Figure 0. \Hy@raisedlink\hyper@@anchor\@currentHrefA comparison of the cumulative ionizing emissivity from this work to that using the assumptions of a fixed escape fraction and ionizing photon production efficiency from Robertson et al. (2015). We indicate with red lines the values of [1, 0.1, 0.01. 0.001] L at 8 (values at other redshifts are comparable; this work assumes dM/dz=0.13). In our model, where the lowest-mass halos have the highest escape fractions, the extreme faint-end of the luminosity function dominates the ionizing emissivity.

The shape of the evolution of from our fiducial model shows a roughly linear evolution with redshift from 6 to 10, with a slight acceleration in the evolution at 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 Robertson et al. (2015), which imply that reionization starts very slowly at 10, and then undergoes a rapid acceleration at 10. This stark difference can be understood by the differences in these two models. The R15 model fits a variable star-formation rate (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[] 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 Conditions for Reionizing the Universe with A Low Galaxy Ionizing Photon Escape Fraction in the Appendix (=13%; log[] 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 .

The differences in and must therefore be responsible. While both certainly play a role, the differences in are easy to understand when simply exploring . 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 non-ionizing UV luminosity density in that same range. The result of this assumption is that the faintest galaxies (14; at 6) do not play a major role. This is illustrated in Figure 4.1.2, where the gray shaded 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 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 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.

The blue shaded lines in Figure 4.1.2 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 to higher values at higher redshifts and lower luminosities, allowing those galaxies which 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 towards the end of reionization, discussed in § 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 9 where the Robertson et al. model, driven by modestly bright galaxies, predicts 0.2, and our model, driven by the faintest galaxies, predicts 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 § 5.1 and Figure 5.1.2 we further compare our results for to several observational and theoretical results in the literature.

4.1.3 He II Reionization

While we did not include any constraints on the ionization of He ii in our model, here we comment on the resultant distribution of , shown as the red-shaded region in Figure 4.1.2. Our fiducial model results in a He ii reionization history which gets started at a low level at 6, as the AGN help to complete hydrogen reionization. The volume ionized fraction of He iii hits 50% at 3.4 0.6, and He ii reionization completes at 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 (e.g. Worseck et al., 2011, 2016a). Additionally, the analogue of Gunn-Peterson troughs observed in the He ii Ly forest have been used to place limits of at the cosmic mean density, which implies that He ii reionization was likely still in progress at (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 , 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; Hiss et al., 2017; Puchwein et al., 2015; Upton Sanderbeck et al., 2016).

Together, these measurements provide evidence that He ii reionization may be ending at , consistent with our findings (see however McQuinn & Worseck, 2014; Davies & Furlanetto, 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\go as much as complete by (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 find evidence for extended reionization in the statistics of the He ii Ly forest opacity222See however D’Aloisio et al. (2017) for a discussion of potential caveats in simulating the impact of He ii reionization on the transmission statistics of the He ii Ly forest..

4.1.4 Comparison with Constraints

In the left panel of Figure 4.1.2 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 4 to 10, consistent with the observations from Becker & Bolton (2013) of an increasing emissivity from 2 to 5. The purple and red shaded regions denote the components of this emissivity contributed by galaxies and AGN, respectively. Although the non-ionizing emissivity from galaxies increases with decreasing redshift from 10 to 4 (Figure 1), the ionizing emissivity does the opposite. This is due to a combination of the faint-end slope of the UV luminosity function becoming shallower, resulting in less luminosity coming from faint galaxies, while it is those very faint galaxies which have the highest escape fractions. This can be better understood by examining Figure 4.1.2. While the cumulative ionizing emissivity at 4 is not highly dissimilar to that at 10, Figure 1 highlights that the 10 universe harbors a greater abundance of extremely faint (15) galaxies than the 4 universe. These factors combine to create a galaxy ionizing emissivity which counterintuitively increases with increasing redshift at 10. This emissivity peaks at 10, and then procedes 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 10, while and continue to decline to higher redshifts, lowering the overall emissivity.

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

4.1.5 Comparison with Constraints

The right panel of Figure 4.1.2 shows the posterior distribution on from our fiducial model. The median value from our model is 0.071 0.005 (integrated to 20).

This is higher than the recent 2016 value published by the Planck collaboration of 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 0.066 0.012, and 1.2 lower than the final WMAP9 value of 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 result in both our model and the results from R15 maintaining consistency with the observational constraints.

\H@refstepcounter

figure \hyper@makecurrentfigure

Figure 0. \Hy@raisedlink\hyper@@anchor\@currentHrefResults from our model when we allow the faint-end slope of the luminosity function to continue evolving to steeper slopes at 10. The transparent gray shading denotes our fiducial model, where the faint-end slope stays fixed at 10 to 2.35. The evolving faint-end slope results in a higher emissivity at 10, which begins reionization slightly sooner. This earlier ionization results in a higher electron scattering optical depth of = 0.077 0.008 (0.7 higher than that from our fiducial model, which was already higher than observations, though not statistically discrepant). This model is a significantly worse fit to the observational constraints than our fiducial model, with the BIC comparison showing positive evidence that the fiducial model is preferred.

Lowering 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 and the observations, it would increase the tension on our other two constraints. Clearly reduced observational uncertainties on will help, though the outlook for significant future improvement is uncertain. We present 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.

4.2 Evolving Versus Flat Faint-End Slope

Our model uses a redshift-dependent parameterization of the UV luminosity function, which stipulates that the characteristic luminosity and number density decrease to increasing redshift, while the faint-end slope becomes steeper. While these parameterizations were developed by fitting to the observational data at 410, in our model we extrapolate to higher redshift. As discussed in § 2.1.1, our fiducial model assumes the faint-end slope ceases to become steeper at 10, assuming 2.35 (the 10 value) at higher redshifts. This choice was made as the redshift evolution of would result in extremely steep slopes of 3.5 by 20. As shown in Figure 1, this results in a higher abundance of 10 galaxies at 15 compared to lower redshifts. Another reason to potentially disfavor such steep values is that the halo mass function is expected to asymptote to a slope of 2 at very low masses. However, at 15 the halo mass function exponential decline begins at log (/M) 8, near to the atomic cooling limit, thus any potential luminous emission relevant here likely originates from a mass regime where the slope is steeper (though of course the luminosity function is further affected by galaxy physics).

Therefore as there is nothing preventing these steep slopes, we consider the results of our model if we allowed to continue to evolve at 10. In this model run, all free parameters were fit with the same priors as our fiducial run, and are free to develop different posteriors to satisfy (if possible) the observational constraints in light of the additional photons available at 10. The right panel of Figure 1 shows how the evolution of the non-ionizing 1500 Å luminosity density evolves for these two scenarios, highlighting that when the faint-end slope continues to steepen, the non-ionizing luminosity density evolves much more shallowly to higher redshift.

We show the results of this analysis in Figure 4.1.5. This model is still consistent with the observational constraints on the ionized volume filling fraction and ionizing emissivity, though the increased number of galaxies at higher redshifts results an overall emissivity which is essentially flat at all redshifts, resulting in a slightly earlier start to reionization, though with a similar midpoint of = 8.9 0.9. This earlier ionization results in a higher electron scattering optical depth of = 0.077 0.008 (0.7 higher than that from our fiducial model). The right-hand panel of Figure 4.1.5 highlights that this value is in tension with the Planck 2016 constraints at the 1.9 level. Unsurprisingly, this model with an evolving faint-end slope has a worse BIC value of 21.0 than our fiducial model, which has 16.9. This difference is large enough to statistically differentiate between these models, with our fiducial model with the fixed faint-end slope at 10 showing “positive” evidence of being preferred (BIC 2). We conclude that our fiducial model is a more plausible scenario, which we discuss in more detail below.\go

\H@refstepcounter table \hyper@makecurrenttable\hb@xt@ Table 0. \Hy@raisedlink\hyper@@anchor\@currentHrefResults of Fiducial Model

Redshift Q log (s Mpc) log / 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 7.2 0.659 0.753 0.836 0.0469 0.0479 0.0490 50.42 50.65 50.82 -0.04 0.93 7.4 0.628 0.724 0.808 0.0482 0.0493 0.0506 50.41 50.65 50.82 0.03 1.11 7.6 0.598 0.699 0.780 0.0494 0.0507 0.0521 50.43 50.66 50.83 0.07 1.21 7.8 0.570 0.672 0.753 0.0506 0.0520 0.0535 50.44 50.67 50.84 0.16 1.42 8.0 0.540 0.647 0.730 0.0517 0.0534 0.0550 50.47 50.69 50.86 0.21 1.74 8.2 0.513 0.617 0.703 0.0529 0.0546 0.0564 50.45 50.68 50.87 0.24 1.92 8.4 0.485 0.590 0.678 0.0539 0.0559 0.0577 50.49 50.70 50.86 0.30 2.25 8.6 0.458 0.563 0.654 0.0549 0.0571 0.0590 50.46 50.68 50.85 0.38 2.46 8.8 0.432 0.533 0.629 0.0558 0.0582 0.0603 50.49 50.72 50.90 0.45 2.83 9.0 0.403 0.502 0.603 0.0568 0.0593 0.0615 50.51 50.72 50.90 0.51 3.24 9.2 0.373 0.476 0.574 0.0577 0.0603 0.0628 50.51 50.74 50.92 0.58 4.42 9.4 0.344 0.443 0.543 0.0584 0.0613 0.0639 50.52 50.74 50.92 0.63 9.6 0.314 0.413 0.507 0.0591 0.0623 0.0650 50.55 50.77 50.95 0.80 9.8 0.285 0.377 0.468 0.0598 0.0632 0.0661 50.57 50.79 50.98 0.90 10.0 0.254 0.342 0.428 0.0605 0.0639 0.0670 50.55 50.79 50.98 1.04 10.2 0.227 0.309 0.389 0.0610 0.0646 0.0679 50.53 50.79 50.96 1.20 10.4 0.202 0.278 0.352 0.0615 0.0653 0.0687 50.51 50.75 50.94 1.32 10.6 0.180 0.247 0.318 0.0620 0.0659 0.0694 50.48 50.73 50.92 1.58 10.8 0.159 0.220 0.287 0.0624 0.0664 0.0701 50.44 50.68 50.90 2.15 11.0 0.140 0.196 0.256 0.0628 0.0669 0.0707 50.42 50.68 50.87 2.78 11.2 0.125 0.174 0.230 0.0632 0.0674 0.0712 50.38 50.62 50.83 4.34 11.4 0.112 0.156 0.207 0.0634 0.0677 0.0717 50.36 50.60 50.81 11.6 0.099 0.139 0.185 0.0637 0.0681 0.0721 50.31 50.57 50.77 11.8 0.087 0.124 0.166 0.0639 0.0684 0.0725 50.29 50.55 50.75 12.0 0.077 0.110 0.147 0.0641 0.0687 0.0729 50.24 50.50 50.71 12.2 0.068 0.098 0.131 0.0643 0.0689 0.0732 50.17 50.44 50.66 12.4 0.060 0.088 0.117 0.0645 0.0691 0.0735 50.18 50.44 50.63 12.6 0.054 0.078 0.104 0.0646 0.0693 0.0738 50.11 50.39 50.60 12.8 0.048 0.070 0.094 0.0647 0.0695 0.0740 50.08 50.33 50.54 13.0 0.043 0.062 0.084 0.0648 0.0697 0.0742 50.05 50.30 50.52 13.2 0.038 0.056 0.075 0.0649 0.0698 0.0744 49.99 50.25 50.47 13.4 0.034 0.051 0.067 0.0650 0.0699 0.0745 49.97 50.22 50.45 13.6 0.031 0.045 0.060 0.0651 0.0701 0.0747 49.98 50.23 50.44 13.8 0.027 0.040 0.053 0.0652 0.0702 0.0748 49.94 50.19 50.40 14.0 0.024 0.036 0.048 0.0653 0.0703 0.0749 49.90 50.14 50.36 14.2 0.021 0.032 0.043 0.0654 0.0704 0.0751 49.88 50.14 50.34 14.4 0.019 0.028 0.038 0.0655 0.0704 0.0752 49.87 50.10 50.29 14.6 0.017 0.025 0.033 0.0655 0.0705 0.0752 49.81 50.05 50.23 14.8 0.015 0.022 0.030 0.0655 0.0706 0.0753 49.76 50.00 50.20 15.0 0.013 0.020 0.026 0.0656 0.0706 0.0754 49.75 50.01 50.19 20.0 0.000 0.000 0.000 0.0658 0.0710 0.0759 48.70 48.93 49.10

\twocolumngrid

5 Discussion

5.1 Comparison of to Model-Dependent Analyses

Observationally constraining the ionization fraction of the IGM towards the end of reionization is a very active area of astrophysics. Here we compare our results for to those from several recent studies, focusing on complementary observational methods involving Ly emission and quasars at 6, as well as theoretical studies, which were not used as constraints on our analysis, summarized in Figure 5.1.2.

5.1.1 Contraints from Ly Emission – Clustering and Luminosity Functions

First we consider constraints inferred via observations of Ly emission. Ly 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é & Rees, 1998; Santos, 2004; Malhotra & Rhoads, 2004; Verhamme et al., 2006; McQuinn et al., 2007; Dijkstra, 2014). Constraints on the evolution of the Ly luminosity function from 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 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 6 were published by Ouchi et al. (2017), who studied the evolution of 2000 Ly emitting galaxies (LAEs) at 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 6.6 appeared slightly more biased than at 5.7, inferring constraints on the ionized IGM gas fraction of 0.85 0.15 at 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. (2017) 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 0.7 0.2 at 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 7.7, finding four candidate galaxies. They concluded that if at least one candidate was a true 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 7.3 a factor of 2 deeper and wider that the 7.7 surveys to find seven LAEs over 0.5 deg, while 65 would have been expected in the case of no evolution from 6.6. They concluded this was evidence for significant evolution, inferring 0.55 0.25 at 7.3. Similar results are found when modeling the evolution of the observed Ly luminosity functions and correlation functions across 5.7 7.3 from Inoue et al. (2018).

While tenuous, these surveys combine to suggest that the IGM is not completely ionized by 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 43 erg s, most spectroscopically confirmed\go, is higher than expected from a Schechter function fit to lower luminosities (e.g., Matthee et al., 2015; Bagley et al., 2017; Zheng et al., 2017; Hu 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 AGN. At a much lower redshift of 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 43 erg s) 7 LAEs have tenuous detections of N v (Tilvi et al., 2016; Hu et al., 2017; Laporte et al., 2017; Mainali et al., 2018).

5.1.2 Contraints from Ly Emission – Spectroscopy

Constraints on the IGM ionization state have also been made using spectroscopic followup of Ly emission from continuum-selected distant star-forming galaxies. For brevity we focus our discussion on more recent results which have ever-increasing statistical confidence, but acknowledge that a significant amount of work has been done in this area (e.g., Stark et al., 2010; Fontana et al., 2010; Stark et al., 2011; Pentericci et al., 2011; Schenker et al., 2012; Ono et al., 2012; Rhoads et al., 2012, 2013; Caruana et al., 2014). Pentericci et al. (2014) observed 12 7 galaxy candidates in the CANDELS EGS field with the FORS2 optical spectrograph on the 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; Ono et al., 2012; Schenker et al., 2012; Bradač et al., 2012) to amass a total sample of 68 candidate 7 galaxies spectroscopically observed with 8–10m telescopes. Over this combined sample, they find that Ly is significantly detected with rest-frame EW 25 Å in 7 of 39 galaxies at 21.25 20.25, and 5 of 25 galaxies at 20.25 18.75. By comparing to reionization models, they find that this detection fraction is consistent with 0.49 (1) at 7, under the assumption that the IGM is fully ionized at 6.

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

Mason et al. (2018a), building on a Bayesian framework introduced in Treu et al. (2012) and Treu et al. (2013), use this new sample of 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 include 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., 2015b; Pentericci et al., 2016; Bradač et al., 2017), allowing a smaller amount of neutral IGM gas to obscure a given Ly line. They use the sample of 6 galaxies from De Barros et al. (2017) and Pentericci et al. (2018) to set a baseline 6 EW distribution measurement, and then explore the evolution to 7 using the sample from Pentericci et al. (2014), finding Q0.41 at 7. This result is model dependent as it draws sightlines from cosmological simulations to realistically model the impact on Ly on all scales. Additionally, Mason et al. (2018a) assumed that the IGM was fully ionized at 6, thus if the implication from Pentericci et al. (2018) is correct, the neutral fraction at 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 non-resonant line redshifts are needed to increase the robustness of these calibrations. Finally, improved statistical power on the IGM at 7 can be obtained by fitting the evolution of EWs down to lower redshift, while also making use of larger samples of 7 galaxies (e.g. Pentericci et al., 2018).

\H@refstepcounter

figure \hyper@makecurrentfigure

Figure 0. \Hy@raisedlink\hyper@@anchor\@currentHrefA comparison of our fiducial model to recent observational and theoretical results. The dark and light shading for our model denote the 1 and 2 confidence ranges, respectively, and a similar shading is used in orange for the results of the Greig & Mesinger (2017) model when using constraints similar to our model. Purple symbols denote constraints from Ly