# [

## Abstract

The reionization of cosmic hydrogen marks a critical juncture in the history of structure formation in the universe. Here we present a new formulation of the standard reionization equation for the time evolution of the volume-averaged H fraction that is more consistent with the accepted conceptual model of inhomogeneous intergalactic absorption. The revised equation retains the basic terminology and simplicity of the classic calculation but explicitly accounts for the presence of the optically thick “Lyman-limit systems” that are known to determine the mean free path of ionizing radiation after overlap. Integration of the modified equation provides a better characterization of the timing of reionization by smoothly linking the pre-overlap with the post-overlap phases of such process. We confirm the validity of the quasi-instantaneous approximation as predictor of reionization completion/maintenance, and discuss new insights on the sources of cosmic reionization using the improved formalism. A constant emission rate into the intergalactic medium (IGM) of about 3 Lyman continuum (LyC) photons per atom per Gyr leads to a reionization history that is consistent with a number of observational constraints on the ionization state of the 5–9 universe and with the reduced Thomson scattering optical depth recently reported by the Planck Collaboration. While star-forming galaxies can dominate the reionization process if the luminosity-weighted fraction of LyC photons that escape into the IGM, , exceeds 15% (for a faint magnitude cut-off of the galaxy UV luminosity function of and a LyC photon yield per unit 1500 Å luminosity of ), simple models where the product of the two unknowns is not evolving with redshift fail to reproduce the changing neutrality of the IGM observed at these epochs.

###### Subject headings:

cosmology: theory — dark ages, reionization, first stars — diffuse radiation — intergalactic medium — galaxiesCosmic Reionization After Planck and Before JWST]Cosmic Reionization After Planck and Before JWST: An Analytic Approach

## 1. Introduction

The transformation of cold neutral intergalactic hydrogen into a highly ionized warm plasma marks the end of the cosmic dark ages and the beginning of the age of galaxies. Studies of resonant absorption in the spectra of distant quasars have shown that hydrogen reionization was still ongoing at and fully completed by redshift 5.5 (e.g., Becker et al., 2015; McGreer et al., 2015; Fan et al., 2006). An early onset of reionization appears to be strongly disfavoured by the latest Planck analysis: the combination of CMB temperature anisotropies with low-multipole polarization data yields an integrated Thomson scattering optical depth, (Planck Collaboration et al., 2016), which is lower and more precise than previous measurements. The ionization state of the 6–8 universe is being constrained by other tracers of reionization history, from the damping wing absorption profiles in the spectra of quasars (e.g., Schroeder et al., 2013) to the luminosity function (LF) and clustering properties of Ly emitting galaxies (LAEs) (e.g., Schenker et al., 2014; Ouchi et al., 2010). Such studies indicate that the intergalactic medium (IGM) was significantly neutral at redshifts between 6 and 7, in agreement with the Planck results.

While a broad consensus on the epoch and duration of the reionization process may be emerging from a variety of astrophysical probes, many key aspects, such as the very nature of the first sources of UV radiation, how they interacted with their environment, and the early thermodynamics of primordial baryonic gas, remain uncertain despite a considerable community effort over the last two decades (for recent reviews on this topic see Haiman, 2016; Lidz, 2016). A detailed modeling of reionization is very challenging, as it requires cosmological numerical simulations that self-consistently couple all the relevant physical processes – dark matter dynamics, gas dynamics, self-gravity, star formation/feedback, radiative transfer, nonequilibrium ionization/recombination, chemical enrichment, heating and cooling – over a large range of scales (e.g., Gnedin, 2014; So et al., 2014). To zeroth order, however, understanding reionization is mainly about a proper accounting of the production and absorption of ionizing LyC radiation as a function of epoch in a clumpy, expanding medium, and this is commonly achieved by solving a “reionization equation” of the form (Madau et al., 1999)

(1) |

for the time evolution of the volume-averaged hydrogen ionized fraction . Here, is the emission rate into the IGM of ionizing photons per unit proper volume, is the cosmological mean proper density of hydrogen, and is an effective recombination timescale. It is this simple ODE that statistically describes the transition from a neutral universe to a fully ionized one, turning reionization into a photon-counting exercise in which the growth rate of H regions is equal to the rate at which ionizing photons are produced minus the rate at which they are consumed by radiative recombinations in the IGM (Madau et al., 1999). Extensively used in the literature (e.g., Haardt & Madau, 2012; Kuhlen & Faucher-Giguère, 2012; Bouwens et al., 2015; Robertson et al., 2015; Madau & Haardt, 2015; Khaire et al., 2016; Ishigaki et al., 2017; Sharma et al., 2017) as it allows an estimate of the photon budget required to achieve reionization with a fast exploration of parameter space, Equation (1) has been shown to provide a good description of the results of radiative transfer simulations (Gnedin, 2000, 2016). With one major limitation: it assumes that all LyC photons escaping from individual galaxies are absorbed by the diffuse IGM, mathematically permitting values of that are above unity when reionization is completed and H regions overlap, which is physically impossible.

At , however, the diffuse IGM is observed to be highly ionized, with only a small residual amount of neutral hydrogen set by the balance between radiative recombinations and photoionizations from a nearly uniform UV radiation background, and provides negligibly small H photoelectric absorption. The continuum opacity is instead dominated by the optically thick (to ionizing radiation) “Lyman-limit systems” (LLSs), high density regions that trace non-linear and collapsed structures, occupy a small portion of the volume, and are able to keep a significant fraction of their hydrogen in neutral form (e.g. Miralda-Escudé et al., 2000; Gnedin & Fan, 2006; Furlanetto & Mesinger, 2009; Haardt & Madau, 2012). In this paper we present a new formulation of the standard reionization equation that explicitly accounts for the presence of these LLSs. We shall see how the LLS opacity causes the mean free path of LyC radiation to remain relatively small even after overlap, and never exceeds unity. The revised ODE retains the basic terminology and simplicity of the classic calculation but is more consistent with the accepted conceptual model of inhomogeneous intergalactic absorption, and provides a better characterization of the timing of reionization by connecting the pre-overlap with the post-overlap epochs. We use the improved formalism for a fresh reassessment of a scenario in which star-forming galaxies dominate the production of LyC radiation in the pre-overlap, overlap, and immediate post-overlap stages of cosmic hydrogen reionization.

Below, we shall adopt a flat cosmology with . The hydrogen and helium mass fractions are and , respectively.

## 2. Basic Formalism

We start by summarizing the basic theory describing the impact of ionizing radiation on an inhomogeneous, primordial IGM. The general idea of a gradual reionization process driven by a steadily increasing UV photon production rate can be cast into a quantitative framework by integrating the rate equations for the fractional abundances of the three species H , He , and He ,

(2) | ||||

(3) | ||||

(4) |

supplemented with three closure conditions for the conservation of charge and of the total abundances of hydrogen and helium. Here, is the photoionization rate of ion , is the radiative recombination coefficient (in units of volume per unit time) of species into , is the proper electron density, and the angle brackets denote an average over all space. The evolution of, e.g., the quantity with redshift fully specifies the global reionization history of cosmic hydrogen. The volume-averaged hydrogen photoionization rate can be written as

(5) |

where is the photoionization cross-section and is the mean monochromatic UV radiation intensity, averaged over all space and directions,

(6) |

Here, is the specific intensity of the radiation field at time and position , measured along the direction . The integral solution of the equation of radiative transfer is given by

(7) |

where , , is the specific photon emission rate into the IGM per unit proper volume, is the cosmic scale factor, and

(8) |

is an “effective optical depth” for photons traveling from to . A packet of photons will travel a proper mean free path before suffering an attenuation.

The mean opacity in Equation (8) is defined as , where the average is taken over all space and all directions (e.g., Gnedin & Ostriker, 1997). This is not the space average of the absorption coefficient, since it is weighted in the radiative transfer equation by the local value of the specific intensity . The relation holds only in the limit of a uniform and isotropic ionizing background, or when and are independent random variables. This is a reasonable approximation after overlap (the point at which H regions merge and the ionizing background rises by a large factor) or during the late stages of reionization, when ionized bubbles expand into low-density regions under the collective influence of many UV sources (e.g., Iliev et al., 2006), the photon mean free path is insensitive to the strength of the local radiation field, and variations in the LyC opacity are relatively modest. The assumption of a nearly uniform UV background, however, breaks down during the early stages of reionization, when rare peaks in the density field that contain more absorbers as well as more sources ionize first, and the photon mean free path is comparable to or smaller than the average source separation. A detailed modeling of the spatial correlations between sources and sinks of ionizing radiation can only be achieved using radiative transfer simulations (e.g., Gnedin, 2014; So et al., 2014) or semi-numerical schemes based on the excursion set formalism (e.g. Furlanetto et al., 2004; Mesinger et al., 2011; Zahn et al., 2011; Kaurov & Gnedin, 2013). It is possible to capture some of these complexities into a single parameter, the “photoionization clumping factor” , that accounts for the cross-correlation between the fractional abundance of H and the radiation field, and write

(9) |

In reionization studies it is conventional to include in the ionization balance equation only a fraction of the hydrogen photoionizations and radiative recombinations, i.e. those that take place in the low-density IGM. The rest are absorbed in-situ, in high-density regions within virialized halos, close to the ionizing sources. As shown by Kohler et al. (2007), when all local ionizations and absorptions are removed from cosmological simulations, the photoionization clumping factor is of order unity close to overlap, and decreases at earlier times as the neutral fraction is anticorrelated with the radiation field.

When the photon mean free path is much smaller than the horizon size, , cosmological effects such as source evolution and frequency shifts can ne neglected. In this local, approximation, the mean specific intensity relaxes to the source function ,

(10) |

and only emitters within the volume defined by an absorption length contribute to the background intensity (Zuo & Phinney, 1993; Madau et al., 1999). While, for a given ionizing photon emissivity, the local source approximation can overestimate the background intensity close to the hydrogen Lyman edge () at (this is because a significant fraction of the emitted LyC photons at these epochs gets redshifted beyond the Lyman limit and does not contribute to the ionizing flux), it quickly approaches the exact value of at the high redshifts, , of interest here (Becker & Bolton, 2013). In order to simplify our calculations, we shall adopt the source function approximation throughout the rest of this paper.

## 3. Continuum Absorption

The absorption opacity in the post-reionization universe is a crucial boundary condition for reionization models, and it is the highly ionized, post-reionization IGM that provides some of the most stringent empirical tests on the nature of cosmological ionizing sources. The technique of stacking quasar spectra provides a direct measurements of the photon proper mean free path at the hydrogen Lyman edge caused by the LLSs,

(11) |

in the interval (Worseck et al., 2014; O’Meara et al., 2013; Prochaska et al., 2009), where the average is taken over all possible quasar sightlines. The rapid evolution of this opacity exceeds that expected from cosmological expansion, indicating an increase in the number density and/or physical size of absorbing structures with redshift. In the limit of optically thick, Poisson-distributed clouds, and neglecting the cumulative photoelectric opacity of lower column density systems, the mean free path is equal to the average spacing between LLSs (“picket-fence” absorption),

(12) |

where is the mean number of absorbers per unit redshift with hydrogen columns , is the Hubble parameter, and is the proper length interval in a Friedmann cosmology. Recent estimates of the mean free path based on the incidence of LLSs along the line of sight are in reasonable agreement with Equation (11) (Faucher-Giguère et al., 2008; Songaila & Cowie, 2010; Haardt & Madau, 2012; Rudie et al., 2013; Prochaska et al., 2014).

To proceed further, we need a model of intergalactic absorption before the epoch of complete reionization.
Our formalism assumes that the early IGM is organized in three main phases: a uniform ionized medium, a uniform neutral medium, and the LLSs, which contain
only a small fraction of the cosmic baryons. UV photons, with their short mean free path, carve out cosmological H regions in the uniform component
that expand in an otherwise largely neutral phase.^{3}

(13) |

This quantity decreases rapidly during reionization as and ionized bubbles merge and overlap. Throughout the volume, however, Lyman-limit absorbers in the outskirt of galaxies will still be consuming UV photons. The presence of these dense neutral clumps will not affect the topology of H regions as long as their mean spacing exceeds the typical ionized bubble radius. Once H regions grow beyond the mean separation between these systems, however, it is the LLSs rather than diffuse gas that determine the mean free path. Under the premises that Equation (11) correctly represents the volume-averaged opacity of LLSs and can be extrapolated to higher redshifts without significant loss in accuracy (since LLSs only affect the late stages of reionization), we can now estimate the volume-averaged, total absorption probability per unit length as

(14) |

Additional intuition on this result can be obtained by writing the effective optical depth per unit redshift interval of a clumpy IGM as

(15) |

where is the joint frequency distribution in redshift and column density of Poisson-distributed absorbers along the line of sight, and is the LyC optical depth through an individual absorber. Expanding the exponential in Equation (15) to first order when and neglecting it compared to unity when gives

(16) |

In the limit in which the LLSs are embedded in a uniform diffuse IGM, the first term on the right-hand side can be written as . Using Equations (12) and (13), and writing the total average opacity as , one then recovers Equation (14).

An expression for the frequency-dependent opacity of LLSs can be derived by noting that the column density distribution of intervening absorbers can be fit with a set of power-laws, , with breaks at pivot points. Around neutral hydrogen columns of cm, uncertainties are still large, with recently estimated slopes at ranging from (Rudie et al., 2013) to (Prochaska et al., 2014). At frequencies , the absorption coefficient can be shown to decrease as (e.g., Madau et al., 1999)

(17) |

where we have approximated the photoionization cross section as . When , this scaling matches that expected in a uniform absorbing medium, . While we shall adopt in all our calculations below, in practice our results are rather insensitive to this choice because of the steep frequency dependence of the photoelectric cross section and galaxy UV spectrum, which also enter in the hydrogen photoionization rate.

We shall use this simple, semi-empirical approach to clumpy IGM absorption to model photon sink processes on large scales during reionization. A different treatment was presented by Miralda-Escudé et al. (2000), who used a volume-weighted density distribution of IGM gas, (where ), to describe the ionization state of an inhomogeneous universe. They assumed that reonization proceeds “outside-in”, first into the underdense voids and then more gradually into overdense regions, and that at any given epoch all the gas with density above a critical overdensity remains neutral, self-shielded from the ionizing background, while all lower density material is completely ionized in a fraction of the volume. Additionally, their model postulates that the absorption mean free path scales as , where is the volume filling factor of gas with . The main advantage of this technique is that it allows one to calculate the mean free path and the “recombination clumping factor” based on a realistic density distribution. The main disadvantage is that numerical simulations appear to support the opposing, “inside-out” view, where reionization proceeds from high- to low-density regions, and the gas density and fractional ionization are correlated on large scales (e.g., Ciardi & Madau, 2003; Iliev et al., 2006; Finlator et al., 2009; Friedrich et al., 2011; So et al., 2014; Bauer et al., 2015). In our model, high-density regions (the LLSs) control the photon mean free path only in the very late stages of reonization and after overlap, while the low-density neutral IGM provides the bulk of the opacity in the early stages. To improve further on our treatment, we will account for inhomogeneities in the ionized diffuse phase via an effective recombination timescale derived from cosmological hydrodynamics simulations.

## 4. Reionization Equation

### 4.1. Hydrogen Photoionization Rate

To make headway towards a reionization equation that can be solved analytically, we start by writing the volume-averaged hydrogen photoionization rate in the local approximation as

(18) |

The photon spectrum of star-forming galaxies between 1 and 4 Ryd can be approximated by a power-law, (e.g., Kewley et al., 2001). Ignoring spatial and temporal variations in the spectral shape of the ionizing photon emissivity, evaluating the integral in Equation (18) analytically for , and using the relation (e.g., Kohler et al., 2007; Oñorbe et al., 2017)

(19) |

we finally obtain an expression for the volume-averaged photoionization rate per hydrogen atom

(20) |

where . Note that the volume-averaged photoionization rate is independent of the photoionization clumping factor , and can therefore be computed without any knowledge of the cross-correlation between the fractional abundance of H and the radiation field.

### 4.2. Radiative Recombination Rate

Let us now recast the hydrogen recombination term in Equation (2) as

(21) |

where is the effective recombination timescale of a clumpy IGM.^{4}

(22) |

where accounts for the presence of photoelectrons from He , is the recombination coefficient at temperature ,
and is the recombination clumping factor that is evaluated with the help of numerical simulations of reionization (which may or may not satisfy all existing
observational constraints). In the above expression, corrects for the enhanced
recombination rate induced by structure formation relative to a highly ionized IGM of uniform density and uniform and constant temperature at all times.
A number of different definitions of the clumping factor exist in the literature, with different cuts in baryon overdensities, ionization and metallicity
levels, and including (or not) the temperature dependence of the recombination coefficient (see, e.g. Kohler et al., 2007; Pawlik et al., 2009; Finlator et al., 2012; Shull et al., 2012; Jeeson-Daniel et al., 2014; Kaurov & Gnedin, 2014; So et al., 2014).
Values of of order a few at are typically derived in recent work if the clumping factor is averaged only over gas with
overdensity (a threshold that excludes dense gas bound to halos as is done for the photoionization clumping factor ), but there are uncertanties.^{5}

A fitting formula to has been recently provided by So et al. (2014) that involves no ad hoc clumping factor,

(23) |

and is based on the analysis of a fully coupled radiation hydrodynamical realization of hydrogen reionization that begins at and completes at . This is the actual, appropriately-averaged recombination timescale in the simulation, again after applying a gas overdensity threshold of . The recombination time is a factor 1.6 shorter than at , and 1.9 times longer at . The two timescales are equal at , and are both shorter than the Hubble time at all redshifts . In an attempt to bracket the uncertainty in the recombination timescale of the IGM, we shall use both and in our numerical estimates below.

### 4.3. Reionization Revisited

Substituting Equations (20) and (21) into the ionization balance equation (2) for hydrogen, we finally obtain

(24) |

This is our revised reionization equation. Prior to overlap, when and , the source term is

(25) |

i.e. the volume-averaged photoionization rate becomes independent of the mean free path and equal to the ionizing photon emission rate into the IGM per hydrogen atom. In this limit, we recover the standard reionization Equation (1). After overlap, as and , the source term becomes

(26) |

i.e. the volume-averaged photoionization rate is reduced by the finite mean free path of LLSs and becomes smaller as the neutral fraction of the IGM decreases, as expected on physical grounds. This new formulation explicitly accounts for the presence of optically thick absorbers that cause the mean free path of LyC photons to remain small even after overlap, and never exceeds unity. The integration of Equation (24) therefore provides a link between the pre-overlap and post-overlap phases of the reionization process.

## 5. The IGM as a Photon Counting Device

We start by considering a simple model for the global ionization history of the universe, with the goal of providing an illustrative description of how reionization may proceed over cosmic time following Equation (24). Figure 1 shows theoretical curves for and obtained by numerically integrating Equation (24) from onwards, assuming a constant emission rate of ionizing photons per hydrogen atom, Gyr. This value corresponds to a comoving ionizing photon emissivity into the IGM of

(27) |

Becker & Bolton (2013) used measurements of the mean Ly and LyC opacity to estimate the ionizing emissivity in the post-reionization era, and found that this is in the range 2–14 Gyr at redshift . The fiducial constant value adopted here at is at the low end of this range, and leads to a reionization history (left panel) that is consistent with a number of observational constraints on the ionization state of the universe, from the redshift-dependent prevalence of LAEs in narrow band surveys at 7–8 (Schenker et al., 2014), to the damping wing absorption profiles measured in the spectra of (Schroeder et al., 2013) and quasars (Greig et al., 2017; Mortlock et al., 2011). The redshift of overlap corresponds to an integrated output into the IGM of about 2.5 ionizing photons per hydrogen atom. The redshift-asymmetric parameterization for the evolution of the ionized fraction adopted in Planck Collaboration et al. (2016),

(28) |

where is the redshift around which the first emitting sources form and is the redshift at which reionization “ends” (), provides a good fit to the reionization histories in the left panel of Figure 1 for and exponent . Denoting with the redshift at which , our models yield a duration of reionization – the pre-overlap phase of Gnedin (2000) – in the range , which is compatible with constraints on patchy reionization based on the kinetic Sunyaev-Zeldovich (kSZ) effect, (95% CL, George et al. 2015). Also plotted in the figure is the recent limit on the redshift of reionization (the redshift at which the hydrogen ionized fraction is 50%) extracted from the Planck CMB data (Planck Collaboration et al., 2016): (uniform prior, redshift-asymmetric parameterization).

The impact of the LLSs in shaping the end of the reionization process and keeping the mean hydrogen neutral fraction of the IGM above after overlap is clearly visible, as well as the sensitivity of the ionization state at these epochs to the recombination timescale. In this toy model, an IGM that recombines with the timescale provides a better fit to observations of the Gunn-Peterson optical depth in the spectra of Sloan Digital Sky Survey (SDSS) quasars at (Fan et al., 2006), as well as to estimates of the ionized fraction at and based on a combination of hydrodynamical simulations with measurements of the Ly forest opacity (Bolton & Haehnelt, 2007). While the pre-overlap stages extend over a considerable range of redshift, the phase of overlap, indicated by the sudden drop in the neutral gas fraction at , is clearly defined even in the presence of the LLSs: decreases by more than 2 orders of magnitude over a fraction of the then Hubble time. Such a dramatic transition marks the epoch when the photon mean free path becomes determined by the LLSs rather than by the typical size of H regions. Note how, in practice, the source term in our reionization equation becomes independent of the actual size of H bubbles in the pre-overlap stages provided this is much smaller than the spacing between LLSs, i.e. in the limit .

The right panel in Figure 1 shows the relative contribution of the source and sink terms to Equation (24). Photoionizations dominate over recombinations until just before overlap. As overlap is approached, the neutral fraction of the IGM dives below and the photoionization source term drops rapidly. Photoionizations come into balance with recombinations for the first time at overlap, where and

(29) |

and remain so thereafter. This relation modifies the original criterion by Madau et al. (1999),

(30) |

for the ionizing photon emissivity necessary to maintain a clumpy, recombining IGM in a highly ionized state in the absence of LLSs.
Just before overlap, however, the factor is very close to unity
and the number of ionizing photons emitted into the IGM in one recombination time is equal to the number of hydrogen atoms, as in the criterion by
Madau et al. (1999).^{6}

It has been recently argued by So et al. (2014) that the time-dependent differential Equation (24) and its quasi-instantaneous version, Equation (29), should not be used to predict the epoch of reionization completion. This is because, according to So et al. (2014), Equation (29) ignores history-dependent terms in the global ionization balance that are not ignorable, and both equations systematically overestimate the redshift of reionization completion compared to their simulations because the conversion of LyC photons into new ionized hydrogen atoms becomes inefficient at late times. The validity of the quasi-instantaneous approximation as predictor of reionization completion/maintenance is, however, clearly shown in the right panel of Figure 3. We believe instead that the disagreement is an artifact of the simulations underestimating the amount of UV absorption in the relatively small box adopted. According to So et al. (2014), the ratio of photoionizations to emitted UV photons decreases as , and is as low as 5% at overlap. By contrast, in our model, even as , the mean free path close to the Lyman edge is much smaller than the horizon because of the presence of LLSs, and the “ionization efficiency” remains 100%. It is the insufficient absorption of LyC radiation by the LLSs, unresolved in the So et al. (2014) simulations, which causes the early flattening in the curve observed by So et al. (2014) and not seen in our analytic calculations.

Together with evidence of patchy ionization in the IGM from an extreme Ly Gunn-Peterson trough at redshift 5.8 (Becker et al., 2015), the observational results listed above strongly suggest that cosmic reionization was “completed” (i.e. the universe became ionized at more than 99%) around redshift 6, and that above redshift 10 hydrogen was ionized to less than 10% levels. In particular, the latest Planck analysis disfavour an early beginning of reionization and yields a Thomson scattering optical depth (Planck Collaboration et al., 2016). The toy models in Figure 1 are consistent by design with the Planck determination, with

(31) |

where the Thomson cross section, and we have ignored the small contribution of photoelectrons from the double reionization of helium at late times. Also note how, taken at face values, the constraints from the Gunn-Peterson opacity, the dark pixel and gap/peak statistics, and the damping wing absorption profiles in the spectra of quasars all appear to require a very sharp transition around redshift 6 from a % neutral to a nearly fully ionized universe.

## 6. Discussion

It is widely believed that normal, star-forming galaxies dominate the production of LyC radiation in the pre-overlap, overlap, and immediate post-overlap stages of cosmic hydrogen reionization. The modeling in the previous sections offers an opportunity for a fresh look at this scenario.

### 6.1. Post-Overlap

In the post-overlap era, defined here as the epoch corresponding to , observations of the Gunn-Peterson optical depth in the spectra of SDSS quasars can be fit with a volume-average neutral fraction (Fan et al., 2006). Equations (11) and (13) imply then

(32) |

Inserting this expression into Equation (29), we derive from the photoionization balance condition after overlap the required comoving ionizing photon emissivity into the IGM,

(33) |

where the range of numerical values covers the redshift interval and the uncertanties in the recombination timescale (see the left panel of Fig. 2).

Over the redshift range , the galaxy UV (1500 Å) comoving luminosity density, , obtained by integrating the observed LF down to the absolute magnitude cut-off , evolves approximately as

(34) |

(“low model”) for mag, and as

(35) |

(“high model”) when the integration is carried down to (Bouwens et al., 2015).
A turn-over in the LF at these faint luminosities is suggested by combining abundance matching with detailed studies of the
color-magnitude diagram of low-luminosity dwarfs in the Local Group (Boylan-Kolchin et al., 2015). According to stellar population synthesis models, the LyC photon yield per
unit UV luminosity at 1500 Å is for a metal-poor stellar population with no dust
(e.g., Schaerer, 2003; Robertson et al., 2013; Madau & Dickinson, 2014).^{7}

(36) |

This is shown in the right panel of Figure 2, where the quantity is expressed in units of , The shadings reflect the uncertanties in the recombination timescale, while the upper and lower swaths correspond to magnitude cut-offs of mag and mag in the luminosity density, respectively. At these epochs, values of 0.04–0.11 can match the data at redshift 5, with a trend toward increasing ionizing photon returns with increasing redshifts. As shown below, the same trend may also hold true in the pre-overlap era.

### 6.2. Pre-Overlap

We can now use our revised reionization equation to bridge the transition from the pre-overlap to the post-overlap era, and ask, e.g., what values of the LyC photon return may, in combination with a given galaxy UV emissivity, reproduce all the observational constraints on the ionization state of the universe. We have therefore integrated Equation (24) numerically from redshift onwards, using the two expressions for the galaxy luminosity density in Equations (34) and (35) together with the two different recombination timescales in an attempt to bracket the uncertainties. Figure 3 shows the resulting reionization histories obtained by assuming a constant LyC photon return into the IGM, in the “high ” case, and in the “low ” case, with higher photon returns naturally leading to earlier epochs of overlap. All these histories are consistent within the errors with the Thomson scattering optical depth measured by the Planck satellite. In the “high ” solution, values of are required for overlap to occur at , while the constraint becomes in the “low ” case. Note how, even for photon production efficiencies as high as (Stanway et al., 2016; Matthee et al., 2017), relative photon leakages into the IGM of 10% appear to be needed. These escape fractions are in qualitative agreement with many recent studies of galaxy-dominated reionization (e.g, Ishigaki et al., 2017; Gnedin, 2016; Khaire et al., 2016; Mitra et al., 2015; Bouwens et al., 2015; Robertson et al., 2015).

As mentioned in the previous section, the kSZ effect provides an upper limit to the duration of reionization. All histories in Figure 3 with are characterized by , which is consistent with current constraints from Planck (Planck Collaboration et al., 2016). Reionization therefore proceeds rather quickly: the shortest durations () are found in the “low ” solutions with . Low luminosity density models typically have lower redshifts of reionization, with 6.9–7.3, compared to “high ” models, with 7.1–7.8.

### 6.3. Uncertainties and Limitations

Cosmic reionization involves a complex interplay between the abundance, clustering, spectrum, and leakage of LyC radiation of photoionizing sources, and the density, clumpiness, temperature, and spatial structure of intergalactic gas, and it is likely that some of the assumptions adopted in our numerical integration of the reionization equation may not be fulfilled in practice. In the previous section we postulated, for example, that normal, star-forming galaxies dominate the production of LyC photons in the pre-overlap, overlap, and immediate post-overlap stages of reionization with a constant ionizing photon return into the IGM. Figure 3 shows, however, that large, redshift-independent LyC photon returns may be difficult to reconcile with the neutral hydrogen fraction observed in the post-overlap IGM. In the “high ” model, in particular, star-forming galaxies can be the source of the cosmic ionizing emissivity at redshift 5 if (for ), a value that is too low for overlap to occur at .

Nearly all models with constant LyC photon returns into the IGM do not recombine rapidly enough at to reproduce the large Gunn-Peterson optical depths observed in the post-overlap universe. And while a number of poorly-known input parameters may affect our calculations at these epochs, from the opacity of the LLSs through the spectrum of ionizing sources to the temperature of the recombining IGM, it is likely that a time-varying photon return may be key for a detailed modeling of the reionization era (see, e.g., Haardt & Madau, 2012; Kuhlen & Faucher-Giguère, 2012; Robertson et al., 2013; Khaire et al., 2016; Sharma et al., 2017). This may be associated with an evolving globally-average escape fraction – with galaxies at being more porous than their lower redshift counterparts – and/or with an increasing LyC photon production efficiency with redshift. The latter has been shown to be as large as in stellar population synthesis models at low metallicity that include the impact of binary stars (Stanway et al., 2016) and in observations of luminous LAEs at 6-7 (Stark et al., 2015; Matthee et al., 2017). A rise in the effective photon production efficiency may also be caused by additional sources of UV radiation at early epochs other than massive stars, like, e.g., active galactic nuclei (AGNs) (e.g., Giallongo et al., 2015; Madau & Haardt, 2015; Chardin et al., 2017; Laporte et al., 2017). We also note that, at redshifts , our extrapolated galaxy UV emissivity is plagued by uncertainties, and does not consider the possibility of a rapidly declining luminosity density in the 260 Myr from to (e.g., Oesch et al., 2015). The inclusion of such a downturn would only strengthen the case for a late reionization epoch.

At a more basic level, we have integrated the revised reionization equation under the assumption that Equation (11) correctly represents the volume-averaged opacity of LLSs. In detail, however, direct measurements of the mean free path at may be biased towards higher values by the quasar proximity effect (D’Aloisio et al., 2016). Optically thin, highly-ionized absorbers with hydrogen columns below (“sub-LLSs”) may also play a non-negligible role in limiting the mean free path (e.g., Haardt & Madau, 2012). Rather than in the source term, we have included the contribution of these sub-LLSs to photon losses in the radiative recombination term.

We have extrapolated observations of the mean free path at to higher redshifts, which may be dangerous if the nature of the LLSs changes during reionization, as the ionizing background becomes weaker (see, e.g., Sobacchi & Mesinger, 2014; Muñoz et al., 2016, and references therein). In our model, LLSs only control the mean free path close to and after overlap, while it is the low-density neutral IGM that provides the bulk of the opacity at . We have checked, for example, that a more rapid decrease with redshift of the LLS mean free path at (vs. Eq. 11) would have little effect on the predicted volume-averaged neutral fraction at . Nevertheless, a firmer understanding of the nature and evolution of the LLSs at these epochs seems essential for an accurate modeling of the late stages of reionization.

Our semi-empirical approach to clumpy IGM absorption produces a sharp overlap phase, a sudden change in the neutral fraction when reionization is almost completed. This is clearly visible in Figures 1 and 3, and marks the epoch when the photon mean free path becomes determined by the LLSs rather than by the typical size of H bubbles. According to Equation (11), the mean separation of the LLSs exceeds 40 comoving Mpc (cMpc) at redshift 6. Therefore, if the typical size of H regions remain smaller than 40 cMpc during reionization, it is the diffuse IGM (with a “swiss-cheese” ionization topology) that will control the photon mean free path, and a well-defined overlap epoch should ensue. A dramatic rise in the mean free path (and in the corresponding UV background intensity) was already present in the first numerical hydrodynamics simulations of reionization (Gnedin, 2000), and is also seen in radiative transfer simulations that account for LLS absorption (Gnedin & Fan, 2006; Gnedin, 2014). One should be cautious, however, about numerical artifacts of studying reionization in simulations of finite box sizes that are comparable or smaller than the mean free path of LyC radiation at . A more gradual transition from bubble-dominated ionization at early times to a smoother “web-dominated” topology characteristic of the post-reionization universe has been argued for by, e.g., Furlanetto & Mesinger (2009). And while a rapid decrease in the volume-weighted neutral hydrogen fraction appear consistent with a number of recent observational constraints on the ionization state of the 5–9 universe, we should recognize that many of them are just educated guesses that rely on constantly changing information.

The resolution of many of the uncertanties discussed above feels overdue and may soon be achieved by the large wavelength coverage, unique sensitivity, and spectroscopic and imaging capabilities of the James Webb Space Telescope together with ongoing and future experiments aimed at measuring the redshifted 21-cm signal from neutral hydrogen during the epoch of reionization.

## Acknowledgements

We thank R. Bouwens, N. Gnedin, and F. Haardt for many useful discussions on the topics presented here, and the anonymous referee for valuable comments and suggestions. Support for this work was provided by NASA through grant HST-AR-13904.001-A (P.M.). P.M. also acknowledges a NASA contract supporting the WFIRST-EXPO Science Investigation Team (15-WFIRST15-0004), administered by GSFC, and thanks the Préfecture of the Ile-de-France Region for the award of a Blaise Pascal International Research Chair, managed by the Fondation de l’Ecole Normale Supérieure.

### Footnotes

- affiliationmark:
- affiliationmark:
- Such a description may not be accurate in the presence of hard-spectrum sources radiating copious penetrating X-ray photons. In the case of a galaxy-dominated reionization scenario, however, X-rays from compact binaries are expected to keep the electron fraction of the IGM between H cavities below 1 percent (Madau & Fragos, 2017).
- This is not the volume-averaged recombination time, which is defined as and weighs preferentially low-density regions.
- Photon losses by LLSs are also due to radiative recombinations. In analytical models it is standard practice, however, to include the effect of LLSs as a reduction in the source term through the finite mean free path of ionizing radiation. Three different quantities – the escape fraction, the clumping factor, and the mean free path – are then used to describe what are essentially radiative recombinations in the ISM, the IGM, and the LLSs. Numerical simulations of cosmic reionization show that such a separation into distinct regimes may indeed be reasonable (Kaurov & Gnedin, 2015).
- In our model, the factor () grows rapidly after overlap. In the case , for example, this factor is 1.007 at (neutral fraction 0.02), 2.15 at (neutral fraction ), and 2.87 at (neutral fraction ). The increase with decreasing redshift is even faster for .
- Note that, in this context, the symbol is sometimes used to express the LyC photon yield per unit rate of star formation instead (e.g. Gnedin, 2016).

### References

- Bauer, A., Springel, V., Vogelsberger, M., et al. 2015, MNRAS, 453, 3593
- Becker, G. D., & Bolton, J. S. 2013, MNRAS, 436, 1023
- Becker, G. D., Bolton, J. S., Madau, P., et al. 2015, MNRAS, 447, 3402
- Bolton, J. S., & Haehnelt, M. G. 2007, MNRAS, 382, 325
- Bouwens, R. J., Illingworth, G. D., Oesch, P. A., et al. 2015, ApJ, 811, 140
- Bouwens, R. J., Smit, R., Labbé, I., et al. 2016, ApJ, 831, 176
- Boylan-Kolchin, M., Weisz, D. R., Johnson, B. D., et al. 2015, MNRAS, 453, 1503
- Chardin, J., Puchwein, E., & Haehnelt, M. G. 2017, MNRAS, 465, 3429
- Ciardi, B., & Madau, P. 2003, ApJ, 596, 1
- D’Aloisio, A., McQuinn, M., Davies, F. B., & Furlanetto, S. R. 2016, ArXiv e-prints, arXiv:1611.02711
- Fan, X., Strauss, M. A., Becker, R. H., et al. 2006, AJ, 132, 117
- Faucher-Giguère, C.-A., Lidz, A., Hernquist, L., & Zaldarriaga, M. 2008, ApJ, 688, 85
- Finlator, K., Oh, S. P., Özel, F., & Davé, R. 2012, MNRAS, 427, 2464
- Finlator, K., Özel, F., Davé, R., & Oppenheimer, B. D. 2009, MNRAS, 400, 1049
- Friedrich, M. M., Mellema, G., Alvarez, M. A., Shapiro, P. R., & Iliev, I. T. 2011, MNRAS, 413, 1353
- Furlanetto, S. R., & Mesinger, A. 2009, MNRAS, 394, 1667
- Furlanetto, S. R., Zaldarriaga, M., & Hernquist, L. 2004, ApJ, 613, 1
- Gallerani, S., Ferrara, A., Fan, X., & Choudhury, T. R. 2008, MNRAS, 386, 359
- George, E. M., Reichardt, C. L., Aird, K. A., et al. 2015, ApJ, 799, 177
- Giallongo, E., Grazian, A., Fiore, F., et al. 2015, A&A, 578, A83
- Gnedin, N. Y. 2000, ApJ, 535, 530
- —. 2014, ApJ, 793, 29
- —. 2016, ApJ, 825, L17
- Gnedin, N. Y., & Fan, X. 2006, ApJ, 648, 1
- Gnedin, N. Y., & Ostriker, J. P. 1997, ApJ, 486, 581
- Greig, B., Mesinger, A., Haiman, Z., & Simcoe, R. A. 2017, MNRAS, 466, 4239
- Haardt, F., & Madau, P. 2012, ApJ, 746, 125
- Haiman, Z. 2016, Understanding the Epoch of Cosmic Reionization: Challenges and Progress, 423, 1
- Hui, L., & Haiman, Z. 2003, ApJ, 596, 9
- Iliev, I. T., Mellema, G., Pen, U.-L., et al. 2006, MNRAS, 369, 1625
- Ishigaki, M., Kawamata, R., Ouchi, M., Oguri, M., & Shimasaku, K. 2017, ArXiv e-prints, arXiv:1702.04867
- Jeeson-Daniel, A., Ciardi, B., & Graziani, L. 2014, MNRAS, 443, 2722
- Kaurov, A. A., & Gnedin, N. Y. 2013, ApJ, 771, 35
- —. 2014, ApJ, 787, 146
- —. 2015, ApJ, 810, 154
- Kewley, L. J., Dopita, M. A., Sutherland, R. S., Heisler, C. A., & Trevena, J. 2001, ApJ, 556, 121
- Khaire, V., Srianand, R., Choudhury, T. R., & Gaikwad, P. 2016, MNRAS, 457, 4051
- Kohler, K., Gnedin, N. Y., & Hamilton, A. J. S. 2007, ApJ, 657, 15
- Kuhlen, M., & Faucher-Giguère, C.-A. 2012, MNRAS, 423, 862
- Laporte, N., Nakajima, K., Ellis, R. S., et al. 2017, ArXiv e-prints, arXiv:1708.05173
- Lidz, A. 2016, in Astrophysics and Space Science Library, Vol. 423, Astrophysics and Space Science Library, ed. A. Mesinger, 23
- Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415
- Madau, P., & Fragos, T. 2017, ApJ, 840, 39
- Madau, P., & Haardt, F. 2015, ApJ, 813, L8
- Madau, P., Haardt, F., & Rees, M. J. 1999, ApJ, 514, 648
- Matthee, J., Sobral, D., Darvish, B., et al. 2017, ArXiv e-prints, arXiv:1706.06591
- McGreer, I. D., Mesinger, A., & D’Odorico, V. 2015, MNRAS, 447, 499
- Mesinger, A., Furlanetto, S., & Cen, R. 2011, MNRAS, 411, 955
- Miralda-Escudé, J., Haehnelt, M., & Rees, M. J. 2000, ApJ, 530, 1
- Mitra, S., Choudhury, T. R., & Ferrara, A. 2015, MNRAS, 454, L76
- Mortlock, D. J., Warren, S. J., Venemans, B. P., et al. 2011, Nature, 474, 616
- Muñoz, J. A., Oh, S. P., Davies, F. B., & Furlanetto, S. R. 2016, MNRAS, 455, 1385
- Oñorbe, J., Hennawi, J. F., & Lukić, Z. 2017, ApJ, 837, 106
- Oesch, P. A., Bouwens, R. J., Illingworth, G. D., et al. 2015, ApJ, 808, 104
- O’Meara, J. M., Prochaska, J. X., Worseck, G., Chen, H.-W., & Madau, P. 2013, ApJ, 765, 137
- Ouchi, M., Shimasaku, K., Furusawa, H., et al. 2010, ApJ, 723, 869
- Pawlik, A. H., Schaye, J., & van Scherpenzeel, E. 2009, MNRAS, 394, 1812
- Planck Collaboration, Adam, R., Aghanim, N., et al. 2016, A&A, 596, A108
- Prochaska, J. X., Madau, P., O’Meara, J. M., & Fumagalli, M. 2014, MNRAS, 438, 476
- Prochaska, J. X., Worseck, G., & O’Meara, J. M. 2009, ApJ, 705, L113
- Robertson, B. E., Ellis, R. S., Furlanetto, S. R., & Dunlop, J. S. 2015, ApJ, 802, L19
- Robertson, B. E., Furlanetto, S. R., Schneider, E., et al. 2013, ApJ, 768, 71
- Rudie, G. C., Steidel, C. C., Shapley, A. E., & Pettini, M. 2013, ApJ, 769, 146
- Schaerer, D. 2003, A&A, 397, 527
- Schenker, M. A., Ellis, R. S., Konidaris, N. P., & Stark, D. P. 2014, ApJ, 795, 20
- Schroeder, J., Mesinger, A., & Haiman, Z. 2013, MNRAS, 428, 3058
- Sharma, M., Theuns, T., Frenk, C., et al. 2017, MNRAS, 468, 2176
- Shull, J. M., Harness, A., Trenti, M., & Smith, B. D. 2012, ApJ, 747, 100
- So, G. C., Norman, M. L., Reynolds, D. R., & Wise, J. H. 2014, ApJ, 789, 149
- Sobacchi, E., & Mesinger, A. 2014, MNRAS, 440, 1662
- Songaila, A., & Cowie, L. L. 2010, ApJ, 721, 1448
- Stanway, E. R., Eldridge, J. J., & Becker, G. D. 2016, MNRAS, 456, 485
- Stark, D. P., Walth, G., Charlot, S., et al. 2015, MNRAS, 454, 1393
- Worseck, G., Prochaska, J. X., O’Meara, J. M., et al. 2014, MNRAS, 445, 1745
- Zahn, O., Mesinger, A., McQuinn, M., et al. 2011, MNRAS, 414, 727
- Zuo, L., & Phinney, E. S. 1993, Astrophysical Journal, 418, 28