# Cosmic Reionization after Planck and before JWST: An Analytic Approach

Piero Madau¹,²

¹Department of Astronomy & Astrophysics, University of California, 1156 High Street, Santa Cruz, CA 95064, USA

²Institut d'Astrophysique de Paris, Sorbonne Universités, UPMC Univ Paris 6 et CNRS, UMR 7095, 98 bis bd Arago, F-75014 Paris, France
Received 2017 July 16; revised 2017 October 11; accepted 2017 October 26; published 2017 December 11

## Abstract

The reionization of cosmic hydrogen marks a critical juncture in the history of structure formation. Here we present a new formulation of the standard reionization equation for the evolution of the volume-averaged H II fraction that is more consistent with the accepted conceptual model of inhomogeneous intergalactic absorption. The revised equation 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 this equation provides a better characterization of the timing of reionization by smoothly linking the pre-overlap with the post-overlap phases of such a process. We confirm the validity of the quasi-instantaneous approximation as a 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 three Lyman continuum (LyC) photons per atom per gigayear leads to a reionization history that is consistent with a number of observational constraints on the ionization state of the z = 5–9 universe. While star-forming galaxies can dominate the reionization process if the luminosity-weighted fraction of LyC photons that escape into the IGM, $f_{esc}$, exceeds 15% (for a faint magnitude cut-off of the galaxy UV luminosity function of $M_{lim} = -13$ and a LyC photon yield per unit 1500 Å luminosity of $\xi_{ion} = 10^{25.3}$ erg$^{-1}$ Hz), simple models where the product of the two unknowns $f_{esc}$ and $\xi_{ion}$ is not evolving with redshift fail to reproduce the changing neutrality of the IGM observed at these epochs.

**Key words:** cosmology: theory – dark ages, reionization, first stars

## 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 z ~ 6 and fully completed by redshift 5.5 (e.g., Fan et al. 2006; Becker et al. 2015; McGreer et al. 2015). An early onset of reionization appears to be strongly disfavored by the latest Planck analysis: the combination of CMB temperature anisotropies with low-multipole polarization data yields an integrated Thomson scattering optical depth, $\tau_{es} = 0.058 \pm 0.012$ (Planck Collaboration et al. 2016), which is lower and more precise than previous measurements. The ionization state of the z = 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$\alpha$ emitting galaxies (LAEs; e.g., Ouchi et al. 2010; Schenker et al. 2014). 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)

$$ \frac{dQ}{dt} = \frac{\langle \dot{n}_{ion} \rangle}{\langle n_{H} \rangle} - \frac{Q}{\bar{t}_{rec}} \qquad (1) $$

for the time evolution of the volume-averaged hydrogen ionized fraction $Q$. Here, $\langle \dot{n}_{ion} \rangle$ is the emission rate into the IGM of ionizing photons per unit proper volume, $\langle n_{H} \rangle = 1.89 \times 10^{-7} (1+z)^3$ cm$^{-3}$ is the cosmological mean proper density of hydrogen, and $\bar{t}_{rec}$ 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 II 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; Madau & Haardt 2015; Robertson et al. 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 Q that are above unity when reionization is completed and H II regions overlap, which is physically impossible.

At $z < 6$, 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 I photoelectric absorption. The continuum opacity is instead dominated by the optically thick (to ionizing radiation) “Lyman-limit systems” (LLSs), high-density regions that trace nonlinear 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 Q 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 ($\Omega_M$, $\Omega_\Lambda$, $\Omega_b$) = (0.3, 0.7, 0.045) flat cosmology with $H_0 = 70$ km s$^{-1}$ Mpc$^{-1}$. The hydrogen and helium mass fractions are $X = 0.75$ and $Y = 0.25$, 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 II, He I, and He III,

$$ \frac{d\langle x_{H\text{II}} \rangle}{dt} = \langle x_{H\text{I}} \Gamma_{H\text{I}} \rangle - \langle n_e x_{H\text{II}} \alpha_{H\text{II}} \rangle, \quad (2) $$

$$ \frac{d\langle x_{\text{He I}} \rangle}{dt} = -\langle x_{\text{He I}} \Gamma_{\text{He I}} \rangle + \langle n_e x_{\text{He II}} \alpha_{\text{He II}} \rangle, \quad (3) $$

$$ \frac{d\langle x_{\text{He III}} \rangle}{dt} = \langle x_{\text{He II}} \Gamma_{\text{He III}} \rangle - \langle n_e x_{\text{He III}} \alpha_{\text{He III}} \rangle, \quad (4) $$

supplemented with three closure conditions for the conservation of charge and of the total abundances of hydrogen and helium. Here, $\Gamma_i$ is the photoionization rate of ion $i$, $\alpha_j(T)$ is the radiative recombination coefficient (in units of volume per unit time) of species $j$ into $i$, $n_e = n_{H\text{II}} + n_{\text{He II}} + 2n_{\text{He III}}$ is the proper electron density, and the angle brackets denote an average over all space. The evolution of, e.g., the quantity

$\langle x_{H\text{II}} \rangle$ with redshift fully specifies the global reionization history of cosmic hydrogen. The volume-averaged hydrogen photoionization rate can be written as

$$ \langle \Gamma_{H\text{I}} \rangle = \int d\nu \frac{4\pi J_\nu}{h\nu} \sigma_{H\text{I}}(\nu), \quad (5) $$

where $\sigma_{H\text{I}}$ is the photoionization cross section and $J_\nu(t)$ is the mean monochromatic UV radiation intensity, averaged over all space and directions,

$$ J_{\nu}(t) \equiv \lim_{V \to \infty} \frac{1}{4\pi V} \int_V d^3x \int d\Omega I_{\nu}(t, \mathbf{x}, \hat{n}). \quad (6) $$

Here, $I_\nu(t, \mathbf{x}, \hat{n})$ is the specific intensity of the radiation field at time $t$ and position $\mathbf{x}$, measured along the direction $\hat{n}$. The integral solution of the equation of radiative transfer is given by

$$ J_{\nu}(t) = \frac{c}{4\pi} \int_{0}^{t} dt' h\nu' \langle \dot{n}_{\nu'}(t') \rangle \left[ \frac{a(t')}{a(t)} \right]^3 e^{-\bar{\tau}(\nu, t, t')}, \quad (7) $$

where $\nu' = \nu a(t)/a(t')$, $\nu'' = \nu a(t)/a(t'')$, $\langle \dot{n}_\nu \rangle$ is the specific photon emission rate into the IGM per unit proper volume, $a(t)$ is the cosmic scale factor, and

$$ \bar{\tau}(\nu, t, t') = c \int_{t'}^{t} dt'' \bar{\kappa}_{\nu''}(t'') \quad (8) $$

is an “effective optical depth” for photons traveling from $t'$ to $t > t'$. A packet of photons will travel a proper mean-free path $1/\bar{\kappa}_\nu$ before suffering an $1/e$ attenuation.

The mean opacity $\bar{\kappa}_\nu$ in Equation (8) is defined as $\bar{\kappa}_\nu \equiv \langle \kappa_\nu I_\nu \rangle / J_\nu$, 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 $I_\nu$. The relation $\bar{\kappa}_\nu = \langle \kappa_\nu \rangle$ holds only in the limit of a uniform and isotropic ionizing background, or when $\kappa_\nu$ and $I_\nu$ are independent random variables. This is a reasonable approximation after overlap (the point at which H II 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” $C_f$, that accounts for the cross-correlation between the fractional abundance of H I and the radiation field, and write

$$ \bar{\kappa}_{\nu} = C_{f} \langle \kappa_{\nu} \rangle. \quad (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 $C_I$ is of the order of 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, $1/\bar{\kappa}_\nu \ll c/H$, cosmological effects such as source evolution and frequency shifts can ne neglected. In this local, $\bar{\tau} \to \infty$ approximation, the mean specific intensity relaxes to the source function $\langle \dot{n}_\nu \rangle / \bar{\kappa}_\nu$,

$$4\pi J_\nu \approx h\nu \frac{\langle \dot{n}_\nu \rangle}{C_I \langle \kappa_\nu \rangle}, \quad (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 ($\nu \gtrsim \nu_L$) at $z \sim 2$ (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 $J_\nu$ at the high redshifts, $z \gtrsim 5$, 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,

$$\langle \kappa_{\nu_L}^{\text{LLS}} \rangle^{-1} = (37 \pm 2) [(1+z)/5]^{-5.4 \pm 0.4} \text{ Mpc} \quad (11)$$

in the interval $2.3 < z < 5.5$ (Prochaska et al. 2009; O’Meara et al. 2013; Worseck et al. 2014), 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),

$$\langle \kappa_{\nu}^{\text{LLS}} \rangle^{-1} = \frac{cH^{-1}}{(1+z)} \left( \frac{dN_{\text{LLS}}}{dz} \right)^{-1}, \quad (12)$$

where $dN_{\text{LLS}}/dz$ is the mean number of absorbers per unit redshift with hydrogen columns $N_{\text{H}1} > \sigma_{\text{H}1}^{-1}(\nu)$, $H(z)$ is the Hubble parameter, and $cH^{-1}/(1+z) = c|dt/dz|$ 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 paths, carve out cosmological H II regions in the uniform component that expand in an otherwise largely neutral phase.<fnref n="3" /> Let us denote, as is customary, the volume filling factor of such cosmological H II regions as $Q$. Within the volume fraction $Q$, only the gas belonging to LLSs remains neutral because of its high density, whereas all the gas is neutral within the fractional volume $(1-Q)$, i.e., $\langle x_{\text{H}1} \rangle = (1-Q)$. The volume-averaged neutral hydrogen density in the uniform IGM is $\langle n_{\text{H}1} \rangle = \langle x_{\text{H}1} \rangle \langle n_{\text{H}} \rangle$, and the corresponding absorption coefficient for photons of frequency $\nu_L < \nu < 4\nu_L$ is

$$\langle \kappa_{\nu}^{\text{IGM}} \rangle = \langle n_{\text{H}} \rangle (1 - Q) \sigma_{\text{H}1}(\nu). \quad (13)$$

This quantity decreases rapidly during reionization as $Q \to 1$ 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 II regions as long as their mean spacing exceeds the typical ionized bubble radius. Once H II 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

$$\langle \kappa_{\nu} \rangle = \langle \kappa_{\nu}^{\text{IGM}} \rangle + \langle \kappa_{\nu}^{\text{LLS}} \rangle. \quad (14)$$

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

$$\frac{d\bar{\tau}}{dz}(\nu) = \int_{0}^{\infty} dN_{\text{H}1} \frac{\partial^2 N}{\partial N_{\text{H}1} \partial z} (1 - e^{-\tau}), \quad (15)$$

where $\partial^2 N/\partial N_{\text{H}1} \partial z$ is the joint frequency distribution in redshift and column density of Poisson-distributed absorbers along the line of sight, and $\tau = N_{\text{H}1}\sigma_{\text{H}1}(\nu)$ is the LyC optical depth through an individual absorber. Expanding the exponential in Equation (15) to first order when $\tau(\nu) < 1$ and

<footnote n="3">³ 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 II cavities below 1% (Madau & Fragos 2017).</footnote>
neglecting it compared to unity when $\tau(\nu) > 1$ gives

$$ \frac{d\bar{\tau}}{dz}(\nu) = \sigma_{H1}(\nu) \int_{0}^{\sigma_{H1}^{-1}(\nu)} N_{H1} dN_{H1} \frac{\partial^2 N}{\partial N_{H1} \partial z} + (dN_{LLS}/dz). \quad (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 $\sigma_{H1}(\nu)\langle x_{H1} \rangle \langle n_H \rangle cH^{-1}/(1+z)$. Using Equations (12) and (13), and writing the total average opacity as $\langle \kappa_{\nu} \rangle \equiv (d\bar{\tau}/dz)|dz/cdt|$, 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, $dN/dN_{H1} \propto N_{H1}^{-\beta}$, with breaks at pivot points. Around neutral hydrogen columns of $N_{H1} \sim 10^{17} \text{cm}^{-2}$, uncertainties are still large, with recently estimated slopes at $\langle z \rangle = 2.4$ ranging from $\beta \simeq 1.48$ (Rudie et al. 2013) to $\beta \simeq 2.11$ (Prochaska et al. 2014). At frequencies $\nu > \nu_L$, the absorption coefficient can be shown to decrease as (e.g., Madau et al. 1999)

$$ \langle \kappa_{\nu}^{\text{LLS}} \rangle = \langle \kappa_{\nu_L}^{\text{LLS}} \rangle (\nu/\nu_L)^{3-3\beta}, \quad (17) $$

where we have approximated the photoionization cross section as $\sigma_{H1} \propto (\nu/\nu_L)^{-3}$. When $\beta = 2$, this scaling matches that expected in a uniform absorbing medium, $\langle \kappa_{\nu}^{\text{IGM}} \rangle \propto \sigma_{H1}(\nu) \propto (\nu/\nu_L)^{-3}$. While we shall adopt $\beta = 2$ 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 photo-ionization 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, $P_V(\Delta_b)$ (where $\Delta_b = \rho_b/\bar{\rho}_b$), to describe the ionization state of an inhomogeneous universe. They assumed that reionization 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 $\Delta_i$ remains neutral, self-shielded from the ionizing background, while all lower density material is completely ionized in a fraction $Q$ of the volume. Additionally, their model postulates that the absorption mean-free path scales as $F_V(\Delta_i)^{-2/3}$, where $F_V(\Delta_i)$ is the volume filling factor of gas with $\Delta_b > \Delta_i$. 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 reionization 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 hydro-dynamics simulations.

## 4. Reionization Equation

### 4.1. Hydrogen Photoionization Rate

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

$$ \begin{aligned} \langle \Gamma_{H1} \rangle &= \int d\nu \frac{\langle \dot{n}_{\nu} \rangle \sigma_{H1}(\nu)}{\bar{\kappa}_{\nu}} \\ &= \int d\nu \frac{\langle \dot{n}_{\nu} \rangle \sigma_{H1}(\nu)}{C_I (\langle \kappa_{\nu}^{\text{IGM}} \rangle + \langle \kappa_{\nu}^{\text{LLS}} \rangle)}. \end{aligned} \quad (18) $$

The photon spectrum of star-forming galaxies between 1 and 4 Ryd can be approximated by a power law, $\dot{n}_{\nu} \propto \nu^{-2}$ (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 $\beta = 2$, and using the relation (e.g., Kohler et al. 2007; Oñorbe et al. 2017)

$$ \langle x_{H1} \Gamma_{H1} \rangle = C_I \langle x_{H1} \rangle \langle \Gamma_{H1} \rangle, \quad (19) $$

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

$$ \langle x_{H1} \Gamma_{H1} \rangle = \frac{\langle \dot{n}_{\text{ion}} \rangle}{\langle n_H \rangle} \frac{1}{(1 + \langle \kappa_{\nu_L}^{\text{LLS}} \rangle / \langle \kappa_{\nu_L}^{\text{IGM}} \rangle)}, \quad (20) $$

where $\langle \dot{n}_{\text{ion}} \rangle = \int_{\nu_L} d\nu \langle \dot{n}_{\nu} \rangle$. Note that the volume-averaged photoionization rate is independent of the photoionization clumping factor $C_I$, and can therefore be computed without any knowledge of the cross-correlation between the fractional abundance of H I and the radiation field.

### 4.2. Radiative Recombination Rate

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

$$ \langle n_e x_{H \text{ II}} \alpha_{H \text{ II}}(T) \rangle \equiv \frac{\langle x_{H \text{ II}} \rangle}{\bar{t}_{\text{rec}}} = \frac{Q}{\bar{t}_{\text{rec}}}, \quad (21) $$

where $\bar{t}_{\text{rec}}$ is the effective recombination timescale of a clumpy IGM.<fnref n="4" /> Because it scales quadratically with density, the volume-averaged recombination rate depends on the actual gas density distribution within the volume. It is often convenient to rewrite $\bar{t}_{\text{rec}}$ as (Madau et al. 1999)

$$ \bar{t}_{\text{rec,M}} = \frac{1}{(1 + \chi) \langle n_H \rangle \alpha(T_0) C_R}, \quad (22) $$

where $\chi = Y/4X = 0.083$ accounts for the presence of photoelectrons from He II, $\alpha$ is the recombination coefficient at temperature $T_0$, and $C_R$ 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, $C_R$ corrects for the enhanced recombination rate induced by structure formation relative to a highly ionized IGM of uniform density

<footnote n="4">⁴ This is not the volume-averaged recombination time, which is defined as $\langle t_{\text{rec}} \rangle = \langle [n_e \alpha_{H \text{ II}}(T)]^{-1} \rangle$ and weighs preferentially low-density regions.</footnote>
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 $C_R$ of the order of a few at $z \gtrsim 6$ are typically derived in recent work if the clumping factor is averaged only over gas with overdensity $\Delta_b < 100$ (a threshold that excludes dense gas bound to halos as is done for the photoionization clumping factor $C_I$), but there are uncertainties.<fnref n="5" /> In the following, we shall insert into $\tilde{t}_{\text{rec,M}}$ the expression $C_R = 2.9[(1+z)/6]^{-1.1}$ from Shull et al. (2012), which gives a clumping factor similar to the one derived by Finlator et al. (2012). We shall also use the Case-A recombination coefficient, $\alpha_A$, as recombination photons are typically redshifted below threshold before being absorbed (Kaurov & Gnedin 2014). The temperature of ionized gas will be fixed to $T_0 = 10^{4.3}$ K to account for the photoheating associated with the reionization process itself (Hui & Haiman 2003). Note that gas at this temperature recombines in Case-A at the same rate of gas in the oft-used Case-B situation at $T_0 = 10^4$ K. Average temperatures higher than $10^4$ K are found for the recombining IGM at $5 < z < 10$ in simulations (So et al. 2014).

A fitting formula to $\tilde{t}_{\text{rec}}$ has recently been provided by So et al. (2014) that involves no ad hoc clumping factor,

$$ \tilde{t}_{\text{rec,S}} = 2.3 [(1+z)/6]^{-4.35} \text{ Gyr}, \quad (23) $$

and is based on the analysis of a fully coupled radiation hydrodynamical realization of hydrogen reionization that begins at $z \simeq 10$ and completes at $z \simeq 5.8$. This is the actual, appropriately averaged recombination timescale in the simulation, again after applying a gas overdensity threshold of $\Delta_b < 100$. The recombination time $\tilde{t}_{\text{rec,M}}$ is a factor of 1.6 shorter than $\tilde{t}_{\text{rec,S}}$ at $z = 6$, and 1.9 times longer at $z = 10$. The two timescales are equal at $z = 7.5$, and are both shorter than the Hubble time at all redshifts $z \gtrsim 6$. In an attempt to bracket the uncertainty in the recombination timescale of the IGM, we shall use both $\tilde{t}_{\text{rec,M}}$ and $\tilde{t}_{\text{rec,S}}$ 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

$$ \frac{dQ}{dt} = \frac{\langle \dot{n}_{\text{ion}} \rangle}{\langle n_{\text{H}} \rangle (1 + \langle \kappa_{\nu_L}^{\text{LLS}} \rangle / \langle \kappa_{\nu_L}^{\text{IGM}} \rangle)} - \frac{Q}{\tilde{t}_{\text{rec}}}. \quad (24) $$

This is our revised reionization equation. Prior to overlap, when $Q < 1$ and $\langle \kappa_{\nu_L}^{\text{LLS}} \rangle \ll \langle \kappa_{\nu_L}^{\text{IGM}} \rangle$, the source term is

$$ \langle x_{\text{H}1} \Gamma_{\text{H}1} \rangle \sim \frac{\langle \dot{n}_{\text{ion}} \rangle}{\langle n_{\text{H}} \rangle}, \quad (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 $Q \to 1$ and $\langle \kappa_{\nu_L}^{\text{LLS}} \rangle \gg \langle \kappa_{\nu_L}^{\text{IGM}} \rangle$, the source term becomes

$$ \langle x_{\text{H}1} \Gamma_{\text{H}1} \rangle \sim \frac{\langle \dot{n}_{\text{ion}} \rangle \langle \kappa_{\nu_L}^{\text{IGM}} \rangle}{\langle n_{\text{H}} \rangle \langle \kappa_{\nu_L}^{\text{LLS}} \rangle} \propto (1 - Q), \quad (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 $Q$ 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 $(1 - Q)$ and $dQ/dt$ obtained by numerically integrating Equation (24) from $z = 12$ onwards, assuming a constant emission rate of ionizing photons per hydrogen atom, $\langle \dot{n}_{\text{ion}} \rangle / \langle n_{\text{H}} \rangle = 2.9 \text{ Gyr}^{-1}$. This value corresponds to a comoving ionizing photon emissivity into the IGM of

$$ \langle \dot{n}_{\text{ion}} \rangle \equiv \langle \dot{n}_{\text{ion}} \rangle / (1+z)^3 = 5.1 \times 10^{50} \text{ s}^{-1} \text{ Mpc}^{-3}. \quad (27) $$

Becker & Bolton (2013) used measurements of the mean Ly$\alpha$ and LyC opacity to estimate the ionizing emissivity in the post-reionization era, and found that this is in the range $\langle \dot{n}_{\text{ion}} \rangle / \langle n_{\text{H}} \rangle = 2-14 \text{ Gyr}^{-1}$ at redshift 4.75. The fiducial constant value adopted here at $5 < z < 12$ 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 $z > 5$ universe, from the redshift-dependent prevalence of LAEs in narrowband surveys at $z = 7-8$ (Schenker et al. 2014), to the damping wing absorption profiles measured in the spectra of $z = 6.3$ (Schroeder et al. 2013) and $z = 7.1$ quasars (Mortlock et al. 2011; Greig et al. 2017). 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),

$$ Q = \left( \frac{z_{\text{early}} - z}{z_{\text{early}} - z_{\text{end}}} \right)^{\alpha}, \quad (28) $$

<footnote n="5">⁵ 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).</footnote>
<figure><img src="image_1.png" /><figcaption>**Figure 1.** Reionization histories predicted by integrating Equation (24) with a constant emission rate of ionizing photons per hydrogen atom, $\langle \dot{n}_{\text{ion}} \rangle / \langle n_{\text{H}} \rangle = 2.9 \text{ Gyr}^{-1}$. Left panel: average neutral fraction ($1 - Q$) of an IGM that recombines according to Equation (22) (solid curve) or Equation (23) (dashed curve). In both models the photon mean-free path includes a contribution from the LLSs. The dotted line shows a model without LLSs and with $\bar{t}_{\text{rec}} = \bar{t}_{\text{rec,S}}$. The curves illustrate that the process of reionization is quite extended: even if it completes at redshift 6, the typical location in the IGM is 50% likely to be fully ionized already at $z \sim 8$. The data points represent constraints on the ionization state of the IGM from the Gunn-Peterson optical depth at $5 < z < 6.1$ (red open circles, Fan et al. 2006), measurements of the Ly$\alpha$ forest opacity combined with hydrodynamical simulations (green triangles, Bolton & Haehnelt 2007), the dark pixel statistics at $z = 5.6$ and $z = 5.9$ (blue squares, McGreer et al. 2015), the gap/peak statistics at $z = 6.32$ (magenta square, Gallerani et al. 2008), the damping wing absorption profiles in the spectra of quasars at $z = 6.3$ (orange pentagon, Schroeder et al. 2013) and $z = 7.1$ (turquoise hexagon, Greig et al. 2017; Mortlock et al. 2011), the redshift-dependent prevalence of LAEs in narrowband surveys at $z = 7$ and $z = 8$ (firebrick stars, Schenker et al. 2014), and their clustering properties at $z = 6.6$ (gold hexagon, Ouchi et al. 2010). The recent limit on the redshift at which $\langle x_{\text{H}} \rangle = 0.5$, extracted from the Planck CMB data (Planck Collaboration et al. 2016), is plotted as the green dot. Right panel: relative contribution of the photoionization source term, $\langle x_{\text{H}1} \Gamma_{\text{H}1} \rangle$ (red lines), and sink recombination term, $Q/\bar{t}_{\text{rec}}$ (blue lines), in the reionization equation. The solid and dashed curves correspond to the different expressions for the recombination timescale $\bar{t}_{\text{rec,M}}$ and $\bar{t}_{\text{rec,S}}$, respectively. Photoionizations dominate over recombinations until just before overlap; the two terms come swiftly into balance at overlap and remain in equilibrium thereafter.</figcaption></figure>

where $z_{\text{early}} = 20$ is the redshift around which the first emitting sources form and $z_{\text{end}}$ is the redshift at which reionization “ends” ($Q = 0.99$), provides a good fit to the reionization histories in the left panel of Figure 1 for $10 > z > z_{\text{end}} = 6$ and exponent $\alpha = 4$. Denoting with $z_{20}$ the redshift at which $Q = 0.2$, our models yield a duration of reionization—the pre-overlap phase of Gnedin (2000)—in the range $z_{20} - z_{\text{end}} = 4.0-4.3$, which is compatible with constraints on patchy reionization based on the kinetic Sunyaev-Zeldovich (kSZ) effect, $z_{20} - z_{\text{end}} < 5.4$ (95% CL, George et al. 2015). Also plotted in the figure is the recent limit on the redshift of reionization $z_{50}$ (the redshift at which the hydrogen ionized fraction is 50%) extracted from the Planck CMB data (Planck Collaboration et al. 2016): $z_{50} = 8.0_{-1.1}^{+0.9}$ (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 $10^{-5}$ 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 $\bar{t}_{\text{rec,M}}$ provides a better fit to observations of the Gunn-Peterson optical depth in the spectra of Sloan Digital Sky Survey (SDSS) quasars at $5 < z < 6.1$ (Fan et al. 2006), as well as to estimates of the ionized fraction at $z = 5$ and 6 based on a combination of hydrodynamical simulations with measurements of the Ly$\alpha$ 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 $z \simeq 5.9$, is clearly defined even in the presence of the LLSs: $\langle x_{\text{H}1} \rangle$ decreases by more than two 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 II regions. Note how, in practice, the source term in our reionization equation becomes *independent* of the actual size of H II bubbles in the pre-overlap stages provided this is much smaller than the spacing between LLSs, i.e., in the limit $\langle \kappa_{\nu_L}^{\text{IGM}} \rangle \gg \langle \kappa_{\nu_L}^{\text{LLS}} \rangle$.

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 $10^{-4}$ and the photoionization source term drops rapidly. Photoionizations come into balance with recombinations for the first time at overlap, where $dQ/dt \to 0$ and

$$ \langle \dot{n}_{\text{ion}} \rangle \bar{t}_{\text{rec}} = \langle n_{\text{H}} \rangle (1 + \langle \kappa_{\nu_L}^{\text{LLS}} \rangle / \langle \kappa_{\nu_L}^{\text{IGM}} \rangle), \quad (29) $$

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

$$ \langle \dot{n}_{\text{ion}} \rangle \bar{t}_{\text{rec}} = \langle n_{\text{H}} \rangle, \quad (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 $(1 + \langle \kappa_{\nu_L}^{\text{LLS}} \rangle / \langle \kappa_{\nu_L}^{\text{IGM}} \rangle)$ 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).<fnref n="6" />

<footnote n="6">⁶ In our model, the factor $(1 + \langle \kappa_{\nu_L}^{\text{LLS}} \rangle / \langle \kappa_{\nu_L}^{\text{IGM}} \rangle)$ grows rapidly after overlap. In the case $\bar{t}_{\text{rec}} = \bar{t}_{\text{rec,M}}$, for example, this factor is 1.007 at $z = 6.1$ (neutral fraction 0.02), 2.15 at $z = 6.0$ (neutral fraction $10^{-4}$), and 2.87 at $z = 5$ (neutral fraction $5 \times 10^{-5}$). The increase with decreasing redshift is even faster for $\bar{t}_{\text{rec}} = \bar{t}_{\text{rec,S}}$.</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 2.** Constraining the sources of photoionizations in the post-overlap era. Left panel: the comoving ionizing photon emissivity into the IGM that is required to reproduce the Gunn–Peterson optical depth measured at $5 < z < 6$ in the spectra of SDSS quasars (Fan et al. 2006). The shading highlights the impact of using different expressions for the IGM recombination timescale, $\bar{\tau}_{\text{rec}} = \bar{\tau}_{\text{rec,M}}$ (upper envelope), and $\bar{\tau}_{\text{rec}} = \bar{\tau}_{\text{rec,S}}$ (lower envelope). Right panel: limits on the “LyC photon return into the IGM” of star-forming galaxies, $f_{\text{esc}} \xi_{\text{ion}}$, set by the same data. The two models of the galaxy comoving UV luminosity density given in Equations (34) and (35) are shown with the blue and red bands, respectively. The quantity $\xi_{\text{ion},25.3}$ is expressed in units of $10^{25.3} \text{ erg}^{-1} \text{ Hz}$.</figcaption></figure>

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 $Q \to 1$, and is as low as 5% at overlap. By contrast, in our model, even as $Q \to 1$, 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, that causes the early flattening in the $Q(z)$ 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$\alpha$ 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 disfavors an early beginning of reionization and yields a Thomson scattering optical depth $\tau_{\text{es}} = 0.058 \pm 0.012$ (Planck Collaboration et al. 2016). The toy models in Figure 1 are consistent by design with the Planck determination, with

$$ \tau_{\text{es}} = c \sigma_{\text{T}} \langle n_{\text{H}} \rangle \int_{0}^{12} \frac{(1+z')^2 dz'}{H(z')} Q(1+\chi) \approx 0.06, \quad (31) $$

where $\sigma_{\text{T}}$ 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 value, 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 ~10% 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 $5 < z < 6$, observations of the Gunn–Peterson optical depth in the spectra of SDSS quasars can be fit with a volume-average neutral fraction $\langle x_{\text{H1}} \rangle \sim 7.2 \times 10^{-5} (z - 4.3)$ (Fan et al. 2006). Equations (11) and (13) imply then that

$$ \frac{\langle \kappa_{\nu_L}^{\text{LLS}} \rangle}{\langle \kappa_{\nu_L}^{\text{IGM}} \rangle} = 1.2 \times 10^{-6} \frac{(1+z)^{2.4}}{\langle x_{\text{H1}} \rangle} \sim 0.017 \frac{(1+z)^{2.4}}{(z-4.3)}. \quad (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,

$$ \begin{aligned} \langle \dot{N}_{\text{ion}} \rangle &= \frac{\langle n_{\text{H}} \rangle}{(1+z)^3 \bar{\tau}_{\text{rec}}} (1 + \langle \kappa_{\nu_L}^{\text{LLS}} \rangle / \langle \kappa_{\nu_L}^{\text{IGM}} \rangle) \\ &= 2 - 5 \times 10^{50} \text{ s}^{-1} \text{ Mpc}^{-3}, \end{aligned} \quad (33) $$

where the range of numerical values covers the redshift interval $5 < z < 6$ and the uncertainties in the recombination timescale (see the left panel of Figure 2).

Over the redshift range $4 < z < 10$, the galaxy UV (1500 Å) comoving luminosity density, $\rho_{\text{UV}}$, obtained by integrating the
observed LF down to the absolute magnitude cut-off $M_{lim}$, evolves approximately as

$$ \log \rho_{UV}(\text{erg s}^{-1} \text{Mpc}^{-3} \text{Hz}^{-1}) = 26.20 - 0.16(z - 6) \quad (34) $$

("low $\rho_{UV}$ model") for $M_{lim} = -16$ mag, and as

$$ \log \rho_{UV}(\text{erg s}^{-1} \text{Mpc}^{-3} \text{Hz}^{-1}) = 26.37 - 0.11(z - 6) \quad (35) $$

("high $\rho_{UV}$ model") when the integration is carried down to $M_{lim} = -13$ (Bouwens et al. 2015). A turn-over in the $z \sim 7$ 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 $\log \xi_{ion}/[\text{erg}^{-1} \text{Hz}] = 25.2-25.3$ for a metal-poor stellar population with no dust (e.g., Schaerer 2003; Robertson et al. 2013; Madau & Dickinson 2014).<fnref n="7" /> These values are consistent with those derived in a sample of sub-L* galaxies at $z = 4-5$ by Bouwens et al. (2016) and based on IRAC-inferred H$\alpha$ measurements. We can now write the comoving photon emission rate into the IGM as $\langle \dot{N}_{ion} \rangle = f_{esc} \xi_{ion} \rho_{UV}$, where $f_{esc}$ is the globally averaged relative fraction of LyC photons that escape from individual galaxies and make it into the IGM. The condition for star-forming galaxies to be the source of the cosmic ionizing emissivity after overlap then sets a constraint on the product of the two unknowns $f_{esc} \xi_{ion}$, the "ionizing photon return into the IGM,"

$$ f_{esc} \xi_{ion} = \frac{\langle n_H \rangle (1 + \langle \kappa_{vL}^{LLS} \rangle / \langle \kappa_{vL}^{IGM} \rangle)}{\rho_{UV} (1 + z)^3 \bar{t}_{rec}}. \quad (36) $$

This is shown in the right panel of Figure 2, where the quantity $\xi_{ion,25.3}$ is expressed in units of $10^{25.3} \text{ erg}^{-1} \text{ Hz}$, The shadings reflect the uncertainties in the recombination timescale, while the upper and lower swaths correspond to magnitude cut-offs of $M_{lim} = -16$ mag and $M_{lim} = -13$ mag in the luminosity density, respectively. At these epochs, values of $f_{esc} \xi_{ion,25.3} \simeq 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 $f_{esc} \xi_{ion}$ may, in combination with a given galaxy UV emissivity, reproduce all the observational constraints on the ionization state of the $z > 5$ universe. We have therefore integrated Equation (24) numerically from redshift $z = 12$ 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, $f_{esc} \xi_{ion,25.3} = 0.1, 0.15, 0.20$ in the "high $\rho_{UV}$" case, and $f_{esc} \xi_{ion,25.3} = 0.2, 0.25, 0.30$ in the "low $\rho_{UV}$" 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 $\rho_{UV}$" solution, values of $f_{esc} \xi_{ion,25.3} \gtrsim 0.15$ are required for overlap to occur at $z \gtrsim 5.9$, while the constraint becomes $f_{esc} \xi_{ion,25.3} \gtrsim 0.25$ in the "low $\rho_{UV}$" case. Note how, even for photon production efficiencies as high as $\xi_{ion} = 10^{25.5} \text{ erg}^{-1} \text{ Hz}$ (Stanway et al. 2016; Matthee et al. 2017), relative photon leakages into the IGM of $f_{esc} \gtrsim 10\%$ appear to be needed. These escape fractions are in qualitative agreement with many recent studies of galaxy-dominated reionization (e.g., Bouwens et al. 2015; Mitra et al. 2015; Robertson et al. 2015; Gnedin 2016; Khaire et al. 2016; Ishigaki et al. 2017).

As mentioned in the previous section, the kSZ effect provides an upper limit to the duration of reionization. All histories in Figure 3 with $z_{end} \ge 5.9$ are characterized by $\Delta z \equiv z_{10} - z_{end} < 4$, which is consistent with current constraints from Planck (Planck Collaboration et al. 2016). Reionization therefore proceeds rather quickly: the shortest durations ($\Delta z \lesssim 3.3$) are found in the "low $\rho_{UV}$" solutions with $\bar{t}_{rec} = \bar{t}_{rec,s}$. Low-luminosity density models typically have lower redshifts of reionization, with $z_{50} = 6.9-7.3$, compared to "high $\rho_{UV}$" models, with $z_{50} = 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 $\rho_{UV}$" model, in particular, star-forming galaxies can be the source of the cosmic ionizing emissivity at redshift 5 if $f_{esc} \xi_{ion,25.3} \simeq 0.04$ (for $\bar{t}_{rec} = \bar{t}_{rec,s}$), a value that is too low for overlap to occur at $z \gtrsim 5.9$.

Nearly all models with constant LyC photon returns into the IGM do not recombine rapidly enough at $5 < z < 6$ 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 $z \gtrsim 6$ 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 $\xi_{ion} = 10^{25.5} \text{ erg}^{-1} \text{ Hz}$ 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 $z \sim 6-7$ (Stark et al. 2015; Matthee

<footnote n="7">⁷ Note that, in this context, the symbol $\xi_{ion}$ is sometimes used to express the LyC photon yield per unit rate of star formation instead (e.g., Gnedin 2016).</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 3.** Star-forming galaxies as sources of reionization. The reionization histories predicted by numerically integrating Equation (24) from $z = 12$ onwards are shown for two models (“high” and “low”) of the galaxy comoving UV luminosity density and two different expressions of the IGM recombination timescale. The data points are the same as those in Figure 1. Top left panel: “high $\rho_{UV}$” model with $\bar{t}_{rec} = \bar{t}_{rec,M}$ and $f_{esc} \xi_{ion,25.3} = 0.2, 0.25, 0.3$. Top right panel: same as the top left with $\bar{t}_{rec} = \bar{t}_{rec,S}$. Bottom left panel: “low $\rho_{UV}$” model with $\bar{t}_{rec} = \bar{t}_{rec,M}$ and $f_{esc} \xi_{ion,25.3} = 0.1, 0.15, 0.2$. Bottom right panel: same as the bottom left with $\bar{t}_{rec} = \bar{t}_{rec,S}$. In all panels, the constant LyC photon returns into the IGM, $f_{esc} \xi_{ion,25.3}$, is expressed in units of $10^{25.3} \text{ erg}^{-1} \text{ Hz}$. All histories are consistent within the errors with the Thomson scattering optical depth measured by the Planck satellite (Planck Collaboration et al. 2016). Note how, in the “high $\rho_{UV}$” solution, values of $f_{esc} \xi_{ion,25.3} \gtrsim 0.15$ are required for overlap to occur at $z \gtrsim 5.9$. The constraint becomes $f_{esc} \xi_{ion,25.3} \gtrsim 0.25$ in the “low $\rho_{UV}$” case. Large LyC photon returns do not reproduce the neutral hydrogen fraction observed in the post-overlap IGM at $5 < z < 6$.</figcaption></figure>

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 (e.g., Giallongo et al. 2015; Madau & Haardt 2015; Chardin et al. 2017; Laporte et al. 2017). We also note that, at redshifts $8 < z < 12$, 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 $z = 8$ to $z = 12$ (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 $z \gtrsim 5$ may be biased toward higher values by the quasar proximity effect (D’Aloisio et al. 2016). Optically thin, highly

ionized absorbers with hydrogen columns below $N_{H1} < 1.6 \times 10^{17} \text{ cm}^{-2}$ (“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 $z < 5.5$ 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 $z \gtrsim 6$. We have checked, for example, that a more rapid decrease with redshift of the LLS mean-free path at $z > 5.5$ (versus Equation (11)) would have little effect on the predicted volume-averaged neutral fraction at
$z \ge 6$. 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 HII 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 HII 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 $z \lesssim 6$. 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 appears consistent with a number of recent observational constraints on the ionization state of the $z = 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 uncertainties 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.

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.

## ORCID iDs

Piero Madau https://orcid.org/0000-0002-6336-3293

## 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: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., Graziani, A., Fiore, F., et al. 2015, A&A, 578, A83

* Gnedin, N. Y. 2000, ApJ, 535, 530

* Gnedin, N. Y. 2014, ApJ, 793, 29

* Gnedin, N. Y. 2016, ApJL, 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, Vol. 423 (Dordrecht: Springer), 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:1702.04867

* Jeeson-Daniel, A., Ciardi, B., & Graziani, L. 2014, MNRAS, 443, 2722

* Kaurov, A. A., & Gnedin, N. Y. 2013, ApJ, 771, 35

* Kaurov, A. A., & Gnedin, N. Y. 2014, ApJ, 787, 146

* Kaurov, A. A., & Gnedin, N. Y. 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:1708.05173

* Lidz, A. 2016, in Astrophysics and Space Science Library, Vol. 423, ed. A. Mesinger (Dordrecht: Springer), 23

* Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415

* Madau, P., & Fragos, T. 2017, ApJ, 840, 39

* Madau, P., & Haardt, F. 2015, ApJL, 813, L8

* Madau, P., Haardt, F., & Rees, M. J. 1999, ApJ, 514, 648

* Matthee, J., Sobral, D., Darvish, B., et al. 2017, 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, Natur, 474, 616

* Muñoz, J. A., Oh, S. P., Davies, F. B., & Furlanetto, S. R. 2016, MNRAS, 455, 1385

* 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

* Oñorbe, J., Hennawi, J. F., & Lukić, Z. 2017, ApJ, 837, 106

* 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, ApJL, 705, L113

* Robertson, B. E., Ellis, R. S., Furlanetto, S. R., & Dunlop, J. S. 2015, ApJL, 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, ApJ, 418, 28