Model-independent constraints on the hydrogen-ionizing emissivity at z>6

TL;DR

This paper derives model-independent constraints on the hydrogen-ionizing photon output during reionization at z>6, revealing a decline consistent with galaxy UV luminosity density and providing benchmarks for reionization models.

Contribution

It introduces a Bayesian, non-parametric method to infer ionizing emissivity from reionization data, independent of source properties.

Findings

Ionizing emissivity declines at z>6, consistent with galaxy UV luminosity density.

Steeper decline when including Lyα damping measurements.

Reionization midpoint constrained to z=6.93±0.14.

Abstract

Modelling reionization often requires significant assumptions about the properties of ionizing sources. Here, we infer the total output of hydrogen-ionizing photons (the ionizing emissivity, $\dot{N}_\textrm{ion}$) at $z=4-14$ from current reionization constraints, being maximally agnostic to the properties of ionizing sources. We use a Bayesian analysis to fit for a non-parametric form of $\dot{N}_\textrm{ion}$, allowing us to flexibly explore the entire prior volume. We infer a declining $\dot{N}_\textrm{ion}$ with redshift at $z>6$, which can be used as a benchmark for reionization models. Model-independent reionization constraints from the CMB optical depth and Ly$\alpha$ and Ly$\beta$ forest dark pixel fraction produce $\dot{N}_\textrm{ion}$ evolution ($d\log_{10}\dot{N}_\textrm{ion}/dz|_{z=6\rightarrow8} = -0.31\pm0.35$ dex) consistent with the declining UV luminosity density of galaxies, assuming constant ionizing photon escape fraction and efficiency. Including measurements from Ly$\alpha$ damping of galaxies and quasars produces a more rapid decline: $d\log_{10}\dot{N}_\textrm{ion}/dz|_{z=6\rightarrow8} =-0.44\pm0.22$ dex, steeper than the declining galaxy luminosity density (if extrapolated beyond $M_\mathrm{UV} \lesssim -13$), and constrains the mid-point of reionization to $z = 6.93\pm0.14$.

Equations

$$\frac{\mathrm{d}{\overline{x}}_{\text{hii}}}{\mathrm{d}t}=\frac{{\dot{N}}_{\text{ion}}}{⟨n_{\text{H}}⟩}-\frac{{\overline{x}}_{\text{hii}}}{t_{\text{rec}}}$$

$${\dot{N}}_{\text{ion,G}}=f_{\text{esc}}{\xi }_{\text{ion}}{\rho }_L$$

$${\dot{N}}_{\text{ion,Q}}=f_{\text{esc,Q}}\int d\nu \frac{{ϵ}_{\nu }}{h\nu }$$

$$p({\dot{\mathbf{N}}}_{\text{ion}}|\text{data})=\mathcal{L}(\text{data}|{\dot{\mathbf{N}}}_{\text{ion}})p({\dot{\mathbf{N}}}_{\text{ion}})$$

$${\tau }_{\mathrm{es}}(z)={\int }_0^{z}d{z}^{\prime }{\sigma }_{\text{textsc}t}{n}_e(1+z^{\prime }{)}^2{\overline{x}}_{\text{hii}}(z^{\prime })\frac{c}{H(z^{\prime })}$$

Full text

††pubyear: 2019††pagerange: Model-independent constraints on the hydrogen-ionizing emissivity at $z>6$–A

1 Introduction

The ultra-violet (UV) radiation of the first populations of stars and accreting black holes reionized hydrogen in the intergalactic medium (IGM) in the Universe’s first billion years. Understanding the reionization process enables us to learn about the properties of these first light sources (e.g., Loeb & Barkana, 2001; Robertson et al., 2010).

However, to model reionization and interpret observations, significant assumptions are often made about the properties of these sources. Were they galaxies or quasars? What was their number density? How clustered were they? How hard were their ionizing spectra? A common approach, assuming reionization is driven by galaxies, is to parameterise the comoving density of ionizing photons, $\dot{N}_{\textrm{ion}}$, as a product of: the galaxy UV luminosity density, $\rho_{L}$; a conversion from observed UV photons (rest-frame $\sim 1500$ Å) to unobserved hydrogen-ionizing photons ($<912$ Å), $\xi_{\mathrm{ion}}$; and the fraction of photons which escape galaxies to ionize the IGM, $f_{\textrm{esc}}$ (e.g., Robertson et al., 2010). However, in the reionization epoch, z\,\hbox to0.0pt{\lower 3.5pt\hbox{\mathchar 0\relax\sim}\hss}\raise 1.0pt\hbox{>}\,6, $f_{\textrm{esc}}$ cannot be directly measured due to the high opacity of the IGM to ionizing photons (e.g., Madau, 1995), and the most straightforward way to estimate $\xi_{\mathrm{ion}}$ requires dust-corrected H$\alpha$ emission (Leitherer & Heckman, 1995) which calls for \,\hbox to0.0pt{\lower 3.5pt\hbox{\mathchar 0\relax\sim}\hss}\raise 1.0pt\hbox{>}\,4.6\,\mum spectroscopy e.g. with JWST NIRSpec. Thus these parameters are often treated as constants in models (e.g., Robertson et al., 2015; Bouwens et al., 2015b). One may ask how valid such assumptions are in the context of current constraints on reionization.

In the past year, several independent measurements provided evidence that reionization was a late ($z<9$), moderately extended (\Delta z\,\hbox to0.0pt{\lower 3.5pt\hbox{\mathchar 0\relax\sim}\hss}\raise 1.0pt\hbox{<}\,4), process. The optical depth to cosmic microwave background (CMB) photons provides an integral constraint: the Planck Collaboration et al. (2018, hereafter P18) results suggest reionization’s midpoint is at $z=7.7\pm 0.7$. Observations of reduced Lyman-alpha (Ly$\alpha$) emission from high-redshift sources can be used to measure the neutral content of the IGM (e.g., McQuinn et al., 2007b; Dijkstra et al., 2011; Mason et al., 2018), as Ly$\alpha$ photons are absorbed by neutral hydrogen. Damped Ly$\alpha$ emission from galaxies and quasars currently provide the most precise constraints on the timeline of reionization, with most observations indicating a substantially neutral $z>7$ IGM (Bañados et al., 2017; Davies et al., 2018; Mason et al., 2018; Mason et al., 2019; Hoag et al., 2019; Greig et al., 2019).

The low number density of quasars observed at $z>6$ (e.g., Parsa et al., 2018) has motivated recent work on galaxy-driven reionization. Fixing $f_{\textrm{esc}}=0.2$ and $\xi_{\mathrm{ion}}=10^{25.2}$ Hz/erg, Bouwens et al. (2015b) and Robertson et al. (2015) found $z<10$ galaxies could reionize the IGM by $z\sim 6$, providing undetected faint galaxies contribute significantly (extrapolating the galaxy luminosity function (LF) to $M_{\textsc{uv}}<-13$ with steep LF faint-end slopes, \alpha\,\hbox to0.0pt{\lower 3.5pt\hbox{\mathchar 0\relax\sim}\hss}\raise 1.0pt\hbox{<}\,-2). However, these studies used parametric models, which could rule out possible evolutionary pathways a priori: Bouwens et al. (2015b) parameterised $\log_{10}{\dot{N}_{\textrm{ion}}}\propto z$; and Robertson et al. (2015) modelled the cosmic star formation rate (SFR) density in a parametric form based on z\,\hbox to0.0pt{\lower 3.5pt\hbox{\mathchar 0\relax\sim}\hss}\raise 1.0pt\hbox{<}\,8 measurements.

Other works have explored models for $f_{\textrm{esc}}$ and $\xi_{\mathrm{ion}}$, allowing them to vary with galaxy properties. Finkelstein et al. (2019) used $f_{\textrm{esc}}$ as a function of halo mass, obtained in hydrodynamical simulations, where the lowest mass galaxies have the highest escape fractions, and modelled the galaxy UV luminosity density based on extrapolating $z\leq 10$ UV LF fits. By shifting the load of reionization to the lowest mass galaxies Finkelstein et al. (2019) find a more extended reionization process and relatively flat $\dot{N}_{\textrm{ion}}(z)$, and also require a contribution from quasars at the end of reionization. Duncan & Conselice (2015) modelled the evolution of $f_{\textrm{esc}}$ and $\xi_{\mathrm{ion}}$ as functions of UV luminosity and spectral slope $\beta$ based on stellar population modelling, aiming to reproduce the $\dot{N}_{\textrm{ion}}(z)$ inferred by Bouwens et al. (2015b). Duncan & Conselice (2015) predict that faint blue galaxies produce and emit the most ionizing photons, which softens the requirement to extrapolate the LF in their work. Naidu R. et al. in prep. (2019), based on the empirical galaxy evolution model of Tacchella et al. (2018), model $f_{\textrm{esc}}$ as a function of SFR surface density (see also, Sharma et al., 2017; Seiler et al., 2019), motivated by observations of Lyman Continuum leakers, and $\xi_{\textrm{ion}}$ using stellar population synthesis. With this model, Naidu R. et al. in prep. (2019) find $M_{\mathrm{UV}}<-18$ galaxies contribute $>80\%$ of the reionization budget (see also, Sharma et al., 2018). However, Greig & Mesinger (2017) demonstrated that significant degeneracies exist between reionization parameters, which could not be broken with current reionization constraints, motivating an integrated approach to modelling the timeline and sources of reionization.

In this work, we seek to update our knowledge of \dot{N}_{\textrm{ion}}(z\,\hbox to0.0pt{\lower 3.5pt\hbox{\mathchar 0\relax\sim}\hss}\raise 1.0pt\hbox{>}\,5), by constraining it from recent estimates of the reionization timeline, and ask how consistent simple parametric models are with the allowed form of $\dot{N}_{\textrm{ion}}(z)$. This method enables us to estimate $\dot{N}_{\textrm{ion}}(z)$ at redshifts higher than is possible with the Ly$\alpha$ forest, which becomes too heavily absorbed at $z>6$ (e.g., Fan et al., 2006; Becker & Bolton, 2013). We aim to be maximally agnostic to the evolution of $\dot{N}_{\textrm{ion}}$ and so in a novel step we fit for a non-parametric form of $\dot{N}_{\textrm{ion}}(z)$ allowing us to flexibly explore the allowed space. This offers an advantage over parametric models which rule out physically possible evolution a priori.

The paper is structured as follows: Section 2 describes the reionization history model; Section 3 describes the method for inferring the redshift evolution of $\dot{N}_{\textrm{ion}}$; Section 4 presents our inferred $\dot{N}_{\textrm{ion}}(z)$ and reionization history. We discuss our results in Section 5, and summarise in Section 6.

We use the Planck Collaboration et al. (2015) cosmology: $(\Omega_{\Lambda},\Omega_{\textrm{m}},\Omega_{\textrm{b}},n,\sigma_{8},H_{0})=$ (0.69, 0.31, 0.048, 0.97, 0.81, 68 km s*-1* Mpc*-1*), as implemented in astropy. Magnitudes are in the AB system. Distances, volumes, and densities are comoving. Confidence levels are 68% unless otherwise stated.

2 Modelling the reionization history

Reionization progresses as ionizations overtake recombinations in the IGM. The reionization history: the volume-averaged ionized hydrogen fraction as a function of redshift, ${\overline{x}_{\textrm{{{hii}}}}}(z)$, can be calculated by solving the following differential equation  (e.g., Madau et al., 1999):

[TABLE]

where $\dot{N}_{\textrm{ion}}$ is the number density of ionizing photons in Mpc*-3* s*-1* and $\langle n_{\textrm{H}}\rangle$ is the mean number density of hydrogen atoms. The recombination time of the IGM is $t_{\textrm{rec}}(z)=\left[C\alpha_{B}(T)n_{e}(1+z)^{3}\right]^{-1}$, where $\alpha_{B}(T)$ is the case B hydrogen recombination (i.e. opaque IGM) coefficient, $n_{e}$ is the number density of electrons (assuming singly ionized He), and $C=\langle n_{H}^{2}\rangle/\langle n_{\textrm{H}}\rangle^{2}$ is the ‘clumping factor’ which accounts for inhomogeneity in the IGM (Madau et al., 1999). For computational efficiency we fix $C=3$, motivated by IGM modelling and simulations (e.g., Shull et al., 2012; Finlator et al., 2012; Kaurov & Gnedin, 2015). Appendix A discusses how our results are insensitive to this assumption. Throughout most of this work we use the neutral fraction ${\overline{x}_{\textrm{{{hi}}}}}=1-{\overline{x}_{\textrm{{{hii}}}}}$.

Assuming galaxies produce the bulk of ionizing photons during reionization $\dot{N}_{\textrm{ion}}$ can be expressed as:

[TABLE]

where $f_{\textrm{esc}}$ is the fraction of ionizing photons which escape galaxies to the IGM; $\xi_{\textrm{ion}}$ is the production rate of ionizing photons per UV luminosity, in Hz/ergs, which depends on the stellar populations’ initial mass function, metallicity, age and dust content; and $\rho_{L}$ is the dust-corrected UV luminosity density (e.g., Robertson et al., 2010). These parameters are likely functions of at least mass and redshift, but $f_{\textrm{esc}}$ and $\xi_{\textrm{ion}}$ are commonly treated as constant for simplicity.

The ionizing photon production rate from quasars can be derived from their spectral energy distribution (SED):

[TABLE]

where the integral limits are from $1-4$ Ry (the energy of hydrogen-ionizing photons) and $f_{\textrm{esc,Q}}$ is the ionizing photon escape fraction from quasars (usually assumed $\sim 1$). The UV SED of quasars is assumed to follow a double-power law, $\epsilon_{\nu}\propto\nu^{\alpha}$, with a pivot at 912 Å. For $<912$ Å we use $\alpha=-1.70$ (Madau & Haardt, 2015; Kulkarni et al., 2019). This value was derived from $M_{\textrm{UV}}\sim-27$ $z\sim 2$ quasars (Lusso et al., 2015). $\alpha$ may be steeper for lower luminosity quasars (e.g., -0.56 – -1.4, Stevans et al., 2014; Scott et al., 2004), which could increase the normalisation of $\dot{N}_{\textrm{ion}}$ by up to $\sim 0.5$ dex.

3 Non-parametric inference of $\dot{N}_{\textrm{ion}}$

To be maximally agnostic about the form of $\dot{N}_{\textrm{ion}}$ we fit for it as a non-parametric function of redshift, i.e. $\dot{N}_{\textrm{ion}}(z)$ can take any value at any redshift. We use redshift bins of width $\Delta z=1$, assuming smooth evolution on that scale, and fit for $\dot{\mathbf{N}}_{\textbf{ion}}(\mathbf{z})$ at $\mathbf{z}=\begin{bmatrix}4\;\dots\;14\end{bmatrix}$, i.e. 11 parameters. In the following we use $\dot{N}_{\textrm{ion}}$ to refer to the ionizing emissivity in general and the vector notation $\dot{\mathbf{N}}_{\textbf{ion}}$ to refer specifically to our inferred values of $\dot{N}_{\textrm{ion}}$ in each redshift bin. Using Bayes’ theorem, the posterior probability for $\dot{\mathbf{N}}_{\textbf{ion}}(\mathbf{z})$ is:

[TABLE]

where $\mathcal{L}(\textrm{data}\,|\,\dot{\mathbf{N}}_{\textbf{ion}})$ is the likelihood of obtaining observed data given our model $\dot{\mathbf{N}}_{\textbf{ion}}$ (described in Section 3.1). We set the prior, $p(\dot{\mathbf{N}}_{\textbf{ion}})$, such that:

• •

$\log_{10}\dot{\mathbf{N}}_{\textbf{ion}}(z=5)$ is uniformly distributed between 49 and 53 ($\dot{N}_{\textrm{ion}}$ in units of Mpc*-3* s*-1*), motivated by measurements of $\dot{N}_{\textrm{ion}}$ at z\,\hbox to0.0pt{\lower 3.5pt\hbox{\mathchar 0\relax\sim}\hss}\raise 1.0pt\hbox{<}\,5 (Becker & Bolton, 2013).

• •

The gradient between redshift steps, $d\log_{10}\dot{\mathbf{N}}_{\textbf{ion}}/dz$, is uniformly distributed between -1 and 1 dex. This is motivated by the SFR/luminosity density of galaxies and quasars which fall by $<1$ dex per redshift at these redshifts (e.g., Bouwens et al., 2015a; Finkelstein et al., 2015; Oesch et al., 2018; Kulkarni et al., 2019), even if bright galaxies dominate reionization (Sharma et al., 2018). The luminosity density of observable galaxies falls only $\sim 0.2$ dex per redshift at these redshift, thus we do not expect enormous jumps in $\dot{N}_{\textrm{ion}}(z)$.

Figure 1 shows the posterior is not restricted by this prior.

In each likelihood call we sample $\dot{\mathbf{N}}_{\textbf{ion}}$, linearly interpolate it over a redshift grid $\Delta z=0.2$ and solve the reionization history Equation 1. While e.g. $f_{\textrm{esc}}$ in individual galaxies could fluctuate on these timescales due to supernova feedback (Trebitsch et al., 2017), the average over the ensemble galaxy population should be smooth. The likelihood for the data, given $\dot{\mathbf{N}}_{\textbf{ion}}$, is described in the next section.

To obtain the posterior distribution for $\dot{N}_{\textrm{ion}}(z)$ we use dynamic nested sampling via dynesty111https://dynesty.readthedocs.io/ (Speagle, 2019), with the sampling settings optimised for posterior estimation.

3.1 Measurements used in $\dot{N}_{\textrm{ion}}$ likelihood

We use two sets of reionization constraints to create two likelihood functions. The total likelihood is the product of the individual likelihoods $\mathcal{L}_{i}$ for the below data given our model $\dot{\mathbf{N}}_{\textbf{ion}}$. Following Greig & Mesinger (2017), we define a ‘Gold’ sample of model-independent measurements, and an additional sample of more model-dependent measurements (described below). The ‘Gold’ sample constraints are:

CMB electron scattering optical depth. CMB photons scatter off a fraction $\tau_{\mathrm{es}}$ of free electrons created during reionization, suppressing CMB anisotropies by $\exp(-\tau_{\mathrm{es}})$, below angular scales corresponding to the size of the cosmological horizon at reionization.

[TABLE]

where $c$ is the speed of light, $\sigma_{\textsc{t}}$ is the Thomson scattering cross section, $n_{e}$ is the comoving number of free electrons (assuming doubly-ionized Helium at $z<4$, following Kuhlen & Faucher-Giguere, 2012; Bouwens et al., 2015b, the exact timing of Helium reionization has a negligible impact on our results), and $H(z)$ is the Hubble parameter. As $\tau_{\mathrm{es}}$ is measured at $z\sim 1100$ for the CMB it provides only an integral constraint on reionization and cannot give precise information about its timing.

In each call of our likelihood we solve for $\tau_{\mathrm{es}}(\dot{\mathbf{N}}_{\textbf{ion}},z=15)$ (assuming fully neutral IGM) and calculate the likelihood for obtaining the P18 optical depth: $\tau_{\mathrm{es}}=0.054\pm 0.007$, from our model, assuming Gaussian uncertainties. 2. 2.

Ly$\alpha$ forest dark fraction. The fraction of dark pixels in the Ly$\alpha$ and Ly$\beta$ forest of $z\sim 6$ quasars (hereafter, ‘dark fraction’) provides a model-independent constraint on ${\overline{x}_{\textrm{{{hi}}}}}$, as neutral hydrogen in the IGM produces fully saturated absorption in the Ly$\alpha$ and Ly$\beta$ forest (Mesinger, 2010). This is degenerate with absorption by self-shielded neutral gas post-reionization, so measurements provide upper limits. We use dark fraction measurements by McGreer et al. (2015): ${\overline{x}_{\textrm{{{hi}}}}}(z=5.6)<0.04+0.05\;(1\sigma)$, ${\overline{x}_{\textrm{{{hi}}}}}(z=5.9)<0.06+0.05$ and ${\overline{x}_{\textrm{{{hi}}}}}(z=6.1)<0.38+0.20$. Following Greig & Mesinger (2017) we implement these via a likelihood with a uniform probability for ${\overline{x}_{\textrm{{{hi}}}}}<{\overline{x}_{\textrm{{{hi}}}}}_{\mathrm{lim}}$ and a half-Gaussian distribution with $\mu={\overline{x}_{\textrm{{{hi}}}}}_{\mathrm{,lim}}$ and $\sigma=\sigma_{\mathrm{lim}}$ for higher values of ${\overline{x}_{\textrm{{{hi}}}}}$.

The second likelihood adds the following measurements of ${\overline{x}_{\textrm{{{hi}}}}}(z)$ from Ly$\alpha$ damping of high-redshift sources, which account for patchy reionization in statistically robust methods. ${\overline{x}_{\textrm{{{hi}}}}}$ can be measured from an observed reduction in Ly$\alpha$ as neutral hydrogen in the IGM absorbs Ly$\alpha$ photons. However, such measurements rely on modelling the sources observed as backlights to the neutral IGM:

Ly$\alpha$ equivalent width (EW) distribution. We use ${\overline{x}_{\textrm{{{hi}}}}}$ constraints at $z\sim 7,7.5$ and 8 from Mason et al. (2018); Hoag et al. (2019); Mason et al. (2019). These measurements infer ${\overline{x}_{\textrm{{{hi}}}}}$ from the redshift evolution of the Ly$\alpha$ EW distribution in Lyman-break galaxies via comparison with realistic IGM simulations. As these constraints were relative to $z=6$ (i.e. assuming the IGM was fully ionized at $z=6$) we use the relative measurement $\Delta_{6-z}{\overline{x}_{\textrm{{{hi}}}}}={\overline{x}_{\textrm{{{hi}}}}}(z)-{\overline{x}_{\textrm{{{hi}}}}}(6)$ in our likelihood. We use the posteriors $p({\overline{x}_{\textrm{{{hi}}}}}(z_{\mathrm{meas}}))$ from those works to directly calculate the likelihood $p(\Delta_{6-z}{\overline{x}_{\textrm{{{hi}}}}}\,|\,\dot{\mathbf{N}}_{\textbf{ion}})$ from our model. To account for the redshift uncertainty in these measurements due to using photometric samples ($\Delta z\sim 1$), we draw 1000 redshifts from the photometric redshift distribution for sources used in each measurement, maintaining detailed balance by using same drawn redshifts in every likelihood calculation, and use the median of the likelihoods calculated for the drawn redshifts. 2. 2.

Ly$\alpha$ emitter clustering. Reionization increases the clustering of Ly$\alpha$ emitting galaxies. We use the constraint of ${\overline{x}_{\textrm{{{hi}}}}}(z=6.6)\leq 0.5\,(1\sigma)$ by Sobacchi & Mesinger (2015) obtained from comparing the angular correlation function of Ly$\alpha$ emitters by Ouchi et al. (2010) to reionization simulations. We implement the likelihood as half-Gaussian distribution with $\mu=0$ and $\sigma=0.5$. 3. 3.

Ly$\alpha$ damping wings in quasar spectra. The neutral IGM creates smooth Ly$\alpha$ damping absorption in quasar spectra. We use recent constraints on ${\overline{x}_{\textrm{{{hi}}}}}$ from two $z>7$ quasars by Davies et al. (2018) and Greig et al. (2019). We use the ${\overline{x}_{\textrm{{{hi}}}}}$ posteriors from those works to calculate likelihoods for our model (using the ‘Intermediate HII regions’ IGM simulation from Greig et al., 2019). We note these works obtained slightly different constraints (though consistent within $1-1.5\sigma$) primarily due to differences in modelling the intrinsic quasar spectra. We include both constraints to ‘marginalise’ over differences in modelling.

These Ly$\alpha$ damping ${\overline{x}_{\textrm{{{hi}}}}}$ measurements are all obtained by comparison to large-scale reionizing IGM simulations. By construction these rely on parametric ‘sub-grid’ models for the properties of reionizing sources, which set the timing and morphology of reionization. However, currently these ${\overline{x}_{\textrm{{{hi}}}}}$ measurements are relatively insensitive to the properties of reionizing sources: there is little difference between ${\overline{x}_{\textrm{{{hi}}}}}$ results obtained using different simulation setups (Sobacchi & Mesinger, 2015; Mason et al., 2018; Greig et al., 2019), meaning our results should be robust to the simulation assumptions. This is due to the structure of the reionizing IGM depending primarily on ${\overline{x}_{\textrm{{{hi}}}}}$, with a lesser some dependence on the clustering scale of the primary reionizing sources, and only weak dependence on the redshift of reionization (as the luminosity-weighted power spectra of galaxies does not change much $z\sim 6-10$, McQuinn et al., 2007a). The impact of reionization morphology is mostly diluted in the measurements from galaxies, which span a redshift range greater than typical size of ionizing bubbles ($<\Delta z=0.1$), and in the case of quasars due to the large uncertainties in ${\overline{x}_{\textrm{{{hi}}}}}$ measurements from single objects.

Additional constraints on reionization come from the patchy kinetic Sunyaev-Zel’dovich (kSZ) effect, which requires large 3D simulations to model, and is beyond the scope of this work. Greig & Mesinger (2017) demonstrate current patchy kSZ results (George et al., 2015) favour late reionization, but do not provide more constraining power on ${\overline{x}_{\textrm{{{hi}}}}}(z)$ than $\tau_{\mathrm{es}}$ and the dark fraction.

4 Results

4.1 Evolution of $\dot{N}_{\textrm{ion}}(z)$

Figure 1 shows $\dot{\mathbf{N}}_{\textbf{ion}}(\mathbf{z})$ inferred from the measurements described in Section 3.1. Panel (a) shows draws from our posteriors and the prior, to demonstrate how $\dot{N}_{\textrm{ion}}$ evolution is constrained when the likelihood is included (Section 3.1). In both cases the inferred $\dot{\mathbf{N}}_{\textbf{ion}}$ is a smooth, declining function of redshift, with the tightest posteriors obtained at $z\sim 6-8$ where the majority of the reionization history constraints are found (Figure 2). When the Ly$\alpha$ damping ${\overline{x}_{\textrm{{{hi}}}}}$ constraints are included, the posterior for $\dot{\mathbf{N}}_{\textbf{ion}}$ is noticeably tighter and prefers a more rapid decline with redshift.

This is clearer in panel (b) which shows the distribution of derivatives of $\log_{10}\dot{\mathbf{N}}_{\textbf{ion}}$ between $z=6\rightarrow 8$ from our $\dot{\mathbf{N}}_{\textbf{ion}}$ posteriors. $d\log_{10}\dot{\mathbf{N}}_{\textbf{ion}}/dz|_{z=6\rightarrow 8}=-0.31\pm 0.35$ inferred from the $\tau_{\mathrm{es}}$ and dark fraction constraints, and $-0.44\pm 0.22$ when the Ly$\alpha$ damping constraints are added. We also compare to models by Bouwens et al. (2015b), Finkelstein et al. (2019) and Naidu R. et al. in prep. (2019), also shown in panel (d) and described in more detail below.

Panel (c) compares our inferred $\dot{\mathbf{N}}_{\textbf{ion}}$ to simple ‘Galaxies’ and ‘Quasars’ models. For the ‘Galaxies’ model we take the galaxy UV luminosity density from the Mason et al. (2015) UV luminosity function (LF) model (consistent with observations from $0<z<10$, see also, Tacchella et al., 2013), and calculate $\dot{N}_{\textrm{ion,G}}$ (Equation 2), fixing $f_{\textrm{esc}}=0.2$ and $\log_{10}(\xi_{\mathrm{ion}})=25.2$ for comparison with previous works (Bouwens et al., 2015b; Robertson et al., 2015). To encompass uncertainty in the faint-end cut-off of the galaxy LF we show the range of $\dot{N}_{\textrm{ion,G}}$ allowed by integrating the LFs to $M_{\textsc{uv}}<-13$ (comparable to the atomic cooling limit at $z\sim 8-10$, M_{h}\,\hbox to0.0pt{\lower 3.5pt\hbox{\mathchar 0\relax\sim}\hss}\raise 1.0pt\hbox{>}\,10^{8}\,M_{\odot}, Bromm & Yoshida, 2011, resulting in a flatter luminosity density and $\dot{N}_{\textrm{ion}}$), $M_{\textsc{uv}}<-15$ (the depth reached in deep HST imaging of galaxy cluster lensed fields, e.g., Livermore et al., 2017; Bouwens et al., 2017) and $M_{\textsc{uv}}<-18$. For the quasar model we take the quasar ionizing emissivity at 912 Å from a homogenised compilation by Kulkarni et al. (2019) and assume a power-law quasar UV SED to calculate the hydrogen-ionizing emissivity (Equation 3). We assume the ionizing photon escape fraction from quasars, $f_{\textrm{esc,Q}}=1$, though this may be optimistic (e.g., Micheva et al., 2017), thus our model shows the upper limit of the quasar contribution. To encompass the uncertainty in the quasar LF faint-end cut-off we show the range of $\dot{N}_{\textrm{ion,Q}}$ obtained by integrating to $M_{\textsc{uv}}<-18$ or $M_{\textsc{uv}}<-21$. The galaxy model is consistent with $\dot{\mathbf{N}}_{\textbf{ion}}$ inferred from $\tau_{\mathrm{es}}$ and the dark fraction constraints. However, including the Ly$\alpha$ damping constraints requires a steeper evolution in $\dot{N}_{\textrm{ion}}$ between $z\sim 6-8$, as described above and shown in panel (b). The majority of the posterior inferred for the slope $d\log_{10}\dot{\mathbf{N}}_{\textbf{ion}}/dz|_{z=6\rightarrow 8}$ favours steeper values than the galaxy luminosity density: 72% of the posterior has $d\log_{10}\dot{\mathbf{N}}_{\textbf{ion}}/dz|_{z=6\rightarrow 8}<-0.31$ dex (steeper than the slope of the galaxy luminosity density if $M_{\textsc{uv}}<-18$), 80% of the posterior has $d\log_{10}\dot{\mathbf{N}}_{\textbf{ion}}/dz|_{z=6\rightarrow 8}<-0.25$ dex, steeper than the galaxy luminosity density for $M_{\textsc{uv}}<-15$, and 83% has $d\log_{10}\dot{\mathbf{N}}_{\textbf{ion}}/dz|_{z=6\rightarrow 8}<-0.23$ dex (steeper than the galaxy luminosity density for $M_{\textsc{uv}}<-12$). The inferred steep $\dot{N}_{\textrm{ion}}(z)$ gradient thus makes a model with constant $f_{\textrm{esc}}$ and $\xi_{\mathrm{ion}}$ and/or abundant faint galaxies less likely.

Panel (d) compares our inferred $\dot{\mathbf{N}}_{\textbf{ion}}$ to models by Bouwens et al. (2015b), Finkelstein et al. (2019) and Naidu R. et al. in prep. (2019). The Bouwens et al. (2015b) model was derived using older constraints on the reionization history, and used tighter priors on $\dot{N}_{\textrm{ion}}$, but is mostly consistent with our $\dot{\mathbf{N}}_{\textbf{ion}}$ inferred from the P18 $\tau_{\mathrm{es}}$ and dark fraction. This model is similar to the ‘Galaxies’ model in panel (c), with its evolution consistent with the galaxy luminosity density, dropping $\sim 0.2$ dex per redshift. It is difficult for such a model, where $f_{\textrm{esc}}$ and $\xi_{\mathrm{ion}}$ are constant, to match the $\dot{N}_{\textrm{ion}}(z)$ evolution inferred from the full set of ${\overline{x}_{\textrm{{{hi}}}}}$ constraints which show later, more rapid, reionization – $d\log_{10}\dot{\mathbf{N}}_{\textbf{ion}}/dz|_{z=6\rightarrow 8}=-0.44\pm 0.22$ dex.

The models by Finkelstein et al. (2019) and Naidu R. et al. in prep. (2019) model $f_{\textrm{esc}}$ as functions: of halo mass (with the highest $f_{\textrm{esc}}$ in low mass halos); and SFR surface density (highest $f_{\textrm{esc}}$ in high SFR surface density objects, see also Sharma et al., 2017; Seiler et al., 2019) respectively. The Finkelstein et al. (2019) $\dot{N}_{\textrm{ion}}$ model slightly increases $z\sim 4\rightarrow 10$, due to low mass halos being more prevalent at higher redshifts with a steep LF faint-end slope ($\alpha<-2$). As demonstrated in panel (b), with the $\tau_{\mathrm{es}}$ and dark fraction constraints, 21% of our $\dot{\mathbf{N}}_{\textbf{ion}}$ posterior increases over $z=6\rightarrow 8$ ($d\log_{10}\dot{\mathbf{N}}_{\textbf{ion}}/dz|_{z=6\rightarrow 8}>0$), and only 3% when the Ly$\alpha$ damping constraints are included, making an increasing $\dot{N}_{\textrm{ion}}(z)$ much more unlikely. The Naidu R. et al. in prep. (2019) SFR surface-density $f_{\textrm{esc}}$ model shifts the burden of reionization onto brighter galaxies, producing a steeper rise in $\dot{N}_{\textrm{ion}}$ with decreasing redshift (see also, Sharma et al., 2018; Seiler et al., 2019). We note that the Naidu R. et al. in prep. (2019) model was fit to most of the same constraints as presented here, so is likely to match by construction, but demonstrates a physically-motivated model which can match the inferred evolution in $\dot{N}_{\textrm{ion}}$.

The right axes of panels (c) and (d) show the number of ionizing photons per hydrogen atom released in a Hubble time: $\dot{N}_{\textrm{ion}}/\langle n_{\mathrm{H}}\rangle/H(z)$, calculated at $z=6$. We infer $\sim 5-23$ (68% range) ionizing photons per hydrogen atom from our $\dot{\mathbf{N}}_{\textbf{ion}}$ posteriors including the Ly$\alpha$ damping constraints ($\sim 3-36$ with only $\tau_{\mathrm{es}}$ and the dark fraction), supporting a scenario where ionizing photons are relatively abundant in the late stages of reionization (Becker & Bolton, 2013).

4.2 Reionization history and CMB optical depth

Figure 2 shows the reionization history obtained in our inference and the observational constraints used in the likelihoods (Section 3.1). As required by the likelihoods our model fits the observations by construction, but it is instructive to observe the allowed reionization histories. While only using the most model-independent constraints allows a broad range of reionization histories (see also Greig & Mesinger, 2017), the Ly$\alpha$ damping constraints prefer a later reionization.

The reionization timeline we infer is rapid between $z\sim 6-8$, driven by the Ly$\alpha$ damping measurements, and thus does not overlap significantly with an extended timeline, such as the model by Finkelstein et al. (2019). This is primarily due to our sharply decreasing $\dot{N}_{\textrm{ion}}(z)$. In the Finkelstein et al. (2019) model low mass galaxies with extremely high escape fractions produce flatter $\dot{N}_{\textrm{ion}}(z)$ at early times than we infer (see panel (d) of Figure 1), enabling reionization to start early, while more massive galaxies with lower escape fractions dominate the SFR density at later times, as the galaxy LF faint-end slope flattens, extending reionization. Our results imply that the latest reionization constraints favour an alternative balancing of the ionizing photon budget, likely reducing the dependence on very low mass galaxies with high $f_{\textrm{esc}}$ to complete reionization.

Double reionization (where the neutral fraction dips again at high redshift) is not ruled out a priori in our model and can be seen in a few posterior draws. In most physical models the luminosity density of galaxies steadily declines at $z>10$ (e.g., Mason et al., 2015; Tacchella et al., 2018), making double reionization unlikely. Within our framework, without placing model-dependent priors on the evolution of $\dot{N}_{\textrm{ion}}$, ruling out this scenario requires higher redshift measurements of the reionization timeline.

Figure 3 shows the distributions of reionization’s midpoint and duration obtained from our analysis. We define the midpoint: $z_{\textrm{0.5}}=z({\overline{x}_{\textrm{{{hi}}}}}=0.5)$, and the duration: $\Delta z_{\mathrm{reion}}=z({\overline{x}_{\textrm{{{hi}}}}}=0.9)-z({\overline{x}_{\textrm{{{hi}}}}}=0.1)$. With just the P18 and dark fraction constraints, we obtain $z_{\textrm{0.5}}=7.49_{-0.57}^{+0.71}\,(68\%)$ and $\Delta z_{\mathrm{reion}}=2.26_{-0.99}^{+2.00}$, comparable to the P18 analysis. Including the Ly$\alpha$ damping constraints more tightly constrains the reionization midpoint to later times but does not significantly update its duration: $z_{\textrm{0.5}}=6.93\pm 0.14$ and $\Delta z_{\mathrm{reion}}=2.83_{-0.99}^{+2.40}$. Again, this is due to the current lack of constraints on reionization’s early stages ($z>8$).

Figure 4 shows the electron scattering optical depth to the CMB obtained in our inference. As required in our likelihood, both models are consistent with the P18 $\tau_{\mathrm{es}}$ at $z=14$. Using the full ${\overline{x}_{\textrm{{{hi}}}}}$ constraints we infer $\tau_{\mathrm{es}}=0.053\pm 0.004$.

5 Discussion

The CMB optical depth and dark fraction provide the most model-independent constraint on $\dot{N}_{\textrm{ion}}(z)$. Any model of reionization should produce $\dot{N}_{\textrm{ion}}(z)$ and reionization histories which are consistent with these results. As noted in Section 4.1, the $\dot{N}_{\textrm{ion}}(z)$ we infer from these constraints is consistent with the decline in the galaxy UV luminosity density, with fixed $f_{\textrm{esc}}$ and $\xi_{\textrm{ion}}$, as shown in previous work (e.g., Bouwens et al., 2015b; Robertson et al., 2015). As measurements of the reionization history improve, with further Ly$\alpha$ observations and 21cm experiments, our knowledge of $\dot{N}_{\textrm{ion}}$ will increase and can be more directly decomposed into its constituent sources.

Figure 1 shows the Ly$\alpha$ absorption constraints prefer a more rapid decline of $\dot{N}_{\textrm{ion}}$ over $z\sim 6\rightarrow 8$ to produce later reionization. While a contribution from quasars could be invoked to produce the build-up of $\dot{N}_{\textrm{ion}}$ from $z\sim 8\rightarrow 6$, using the model described in Section 2 the contribution of quasars to $\dot{N}_{\textrm{ion}}$ is negligible at $z>5$: the combined galaxies and quasar models shown in Figure 1 do not produce an evolution as steep as our inferred $\dot{N}_{\textrm{ion}}$. This is true even assuming quasars exist to $M_{\textsc{uv}}<-18$ and have hard ionizing spectra, $\alpha=-0.56$ (Scott et al., 2004). However, the faint-end slope of the quasar LF is still uncertain at $z>7$ (e.g., Manti et al., 2016): a large population of faint quasars at $z\sim 8-10$ could contribute to the increase in $\dot{N}_{\textrm{ion}}$ at $z<8$.

The inferred rapid decline in $\dot{N}_{\textrm{ion}}$ is steeper than the evolution in the galaxy luminosity density (if integrated down to $M_{\textsc{uv}}\geq-13$). While many different effects contribute to $\dot{N}_{\textrm{ion}}$, using this luminosity density and constant $f_{\textrm{esc}}$ and $\xi_{\mathrm{ion}}$ it thus is difficult to describe $\dot{N}_{\textrm{ion}}$ evolution. Our inferred $\dot{N}_{\textrm{ion}}$ therefore provides tentative a posteriori evidence that galaxies’ ionizing photon emission properties evolve with redshift and/or that frequently-solicited, but undetected, faint galaxies contribute less to reionization than previously required (e.g., Bouwens et al., 2015b; Robertson et al., 2015). For example, models where $f_{\textrm{esc}}$ is higher in more massive galaxies will naturally produce a steeper $\dot{N}_{\textrm{ion}}$ evolution, due to the more rapid build-up of high mass galaxies at high redshifts (Sharma et al., 2018; Seiler et al., 2019; Naidu R. et al. in prep., 2019). A steep $\dot{N}_{\textrm{ion}}$ can also be produced and/or enhanced if the galaxy luminosity density is dominated by ‘bright’ galaxies (M_{\textsc{uv}}\,\hbox to0.0pt{\lower 3.5pt\hbox{\mathchar 0\relax\sim}\hss}\raise 1.0pt\hbox{<}\,-18), which is the case if the UV LF faint-end slope is shallow, $\alpha>-2$. Accurate measurements of the faint-end slope at $z>4$ are therefore important for understanding the ionizing photon budget, and will be improved in the next decade with JWST deep imaging observations. The dominant sources of reionization impact the topology of HII regions during reionization, measurable by future 21cm intensity mapping experiments: e.g. bright, massive sources produce more biased and larger ionized bubbles compared to if extremely faint, low mass sources dominate the ionizing budget (e.g., McQuinn et al., 2007a; Seiler et al., 2019).

As demonstrated in Figure 2 the Ly$\alpha$ emitter clustering and EW evolution measurements most tightly constrain the $\dot{\mathbf{N}}_{\textbf{ion}}$ posterior at $z\sim 7$, due to their smaller uncertainties compared to the ${\overline{x}_{\textrm{{{hi}}}}}$ measurements from individual quasars. As discussed in Section 3.1 the Ly$\alpha$ damping measurements may introduce systematics due to modelling the intrinsic Ly$\alpha$ emission (both in galaxies and quasars). In particular, an increase in circumgalactic medium absorption at $z>6$ could also play a role in absorbing Ly$\alpha$ emission from galaxies (e.g., from an increase in self-shielding systems Bolton & Haehnelt, 2013), but is likely subdominant to reionization as explaining the Ly$\alpha$ damping without reionization requires the ionizing background to drop by at least a factor of 20 from $z\sim 6\rightarrow 7$ (Mesinger et al., 2015).

We find constraints from Ly$\alpha$ damping provide increasing evidence that the mid-point of reionization was relatively late, $z\sim 7$. However, the duration of reionization is still mostly unconstrained due to the lack of observations at $z>8$. Physical models of the $z>10$ galaxy population (e.g., Mason et al., 2015; Tacchella et al., 2018) predict a declining luminosity density, due to the lower abundance of halos to host star formation, which should provide continually declining $\dot{N}_{\textrm{ion}}(z)$, but our work demonstrates that this is not yet confirmed observationally. $\dot{N}_{\textrm{ion}}$ is best constrained $z\sim 6-8$, but future $z>8$ Ly$\alpha$ surveys and 21cm experiments will measure ${\overline{x}_{\textrm{{{hi}}}}}$, and thus $\dot{N}_{\textrm{ion}}$, at higher redshifts.

6 Conclusions

We have fit for the non-parametric evolution of the hydrogen-ionizing emissivity, $\dot{N}_{\textrm{ion}}$, at $z=4-14$, being maximally agnostic about the ionizing sources. We use the most recent constraints on hydrogen reionization to constrain our model. This method enables inference of $\dot{N}_{\textrm{ion}}(z)$ at redshifts higher than is possible with the Ly$\alpha$ forest, which becomes too heavily absorbed at $z>6$. Our main conclusions can be summarised as follows:

Current constraints on reionization favour a declining $\dot{N}_{\textrm{ion}}$ with redshift at $z>6$. This is moderately favoured by model-independent $\tau_{\mathrm{es}}$ and dark fraction constraints, $d\log_{10}\dot{\mathbf{N}}_{\textbf{ion}}/dz|_{z=6\rightarrow 8}=-0.31\pm 0.35$ and more strongly favoured when Ly$\alpha$ damping constraints are included ($d\log_{10}\dot{\mathbf{N}}_{\textbf{ion}}/dz|_{z=6\rightarrow 8}=-0.44\pm 0.22$). 2. 2.

$\dot{N}_{\textrm{ion}}(z)$ inferred from model-independent reionization constraints is consistent with the declining UV luminosity density of galaxies, with constant $f_{\textrm{esc}}$ and $\log_{10}\xi_{\mathrm{ion}}$, as found by previous studies. 3. 3.

When reionization constraints from Ly$\alpha$ damping in quasars and galaxies are included, a more rapid decline in $\dot{N}_{\textrm{ion}}(z)$ is inferred, which is less likely to be explained purely by the declining galaxy UV luminosity density (integrated to $M_{\textsc{uv}}>-13$) with constant ionizing photon output of galaxies, relative to their non-ionizing UV emission, $z=8\rightarrow 6$. 4. 4.

Including the Ly$\alpha$ damping measurements significantly constrains the midpoint of reionization to $z_{\textrm{0.5}}=6.93\pm 0.14$ (compared with $7.49_{-0.57}^{+0.71}$ from the P18 optical depth and dark fraction alone), but does not provide more information on its duration due to the lack of observational constraints in reionization’s earliest stages, $z>8$.

Without making assumptions about galaxy evolution, our analysis demonstrates current reionization measurements are broadly consistent with the evolution in the galaxy luminosity density, but hint at evolution in the ionizing photon output at $z>6$ and/or a lessened requirement for undetected extremely faint galaxies to dominate reionization. Future measurements of the reionization history, from galaxy surveys and 21cm experiments, can be used in this framework to more tightly constrain the evolution of $\dot{N}_{\textrm{ion}}$, and combined with statistics of the $z>8$ galaxy population, observable with JWST, and measurements of the reionization topology, expose the sources of reionization.

Acknowledgements

The authors thank Charlie Conroy for useful discussions and for providing comments on a draft of this paper. We thank Steve Finkelstein, Fred Davies and Brad Greig for useful discussions and for providing: the $\dot{N}_{\textrm{ion}}$ model from Finkelstein et al. (2019), and ${\overline{x}_{\textrm{{{hi}}}}}$ posteriors from Davies et al. (2018) and Greig et al. (2019) respectively. We thank Austin Hoag for sharing the ${\overline{x}_{\textrm{{{hi}}}}}$ posterior from Hoag et al. (2019). CAM acknowledges support by NASA Headquarters through the NASA Hubble Fellowship grant HST-HF2-51413.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. RPN gratefully acknowledges an Ashford Fellowship and Peirce Fellowship granted by Harvard University. ST is supported by the Smithsonian Astrophysical Observatory through the CfA Fellowship. JL is supported by an NSF Astronomy and Astrophysics Postdoctoral Fellowship under award AST-1701487.

Software: IPython (Pérez & Granger, 2007), matplotlib (Hunter, 2007), NumPy (Van Der Walt et al., 2011), SciPy (Oliphant, 2007), Astropy (Robitaille et al., 2013), dynesty (Speagle, 2019).

Appendix A Impact of clumping factor

The only parameter in Equation 1 with significant uncertainty is the IGM clumping factor, $C$, which could impact our results. However, as noted by Bouwens et al. (2015b), $\dot{N}_{\textrm{ion}}$ is remarkably insensitive to $C$. We tested our inference with $C=1-30$ (where $C=2-6$ is the physically motivated range based on IGM simulations, e.g., Finlator et al., 2012; Kaurov & Gnedin, 2015), shown in Figure 5.

Increasing the clumping factor shifts the normalisation of $\dot{N}_{\textrm{ion}}$ to higher values, in order to balance increased recombinations in the IGM and complete reionization on time, but does not change its evolution with redshift. All of the resulting $\dot{\mathbf{N}}_{\textbf{ion}}$ posteriors are consistent within the 68% confidence region, and negligibly different for $C\sim 1-10$. We also tested an evolving $C(z)=2.9(1+z/6)^{-1.1}$ (Shull et al., 2012) and found no significant difference from the non-evolving cases. Therefore, our inferred $\dot{\mathbf{N}}_{\textbf{ion}}$ is robust to assumptions about the clumping factor, suggesting the impact of the clumping factor is less than the uncertainties in inferring $\dot{N}_{\textrm{ion}}$ from reionization constraints. All our results on the reionization timeline are insensitive to changes in the clumping factor as $\dot{N}_{\textrm{ion}}$ is modified self-consistently with $C$ to produce ${\overline{x}_{\textrm{{{hi}}}}}(z)$ consistent with observations.