Gaussian Process Reconstruction of Reionization History

In the late matter-dominated epoch of the universe, the neutral hydrogen and helium atoms in the intergalactic medium (IGM) were ionized by the photons emitted from the sources formed during structure formation. Observations of quasar spectroscopy have confirmed that the IGM was ionized before z ∼ 4 with a neutral hydrogen fraction of about 1 in 10,000 as we do not observe the Gunn & Peterson (1965) trough. Simulation of large-scale structure formation also reveals that the sources of high energetic photons needed for reionization could not have formed before z ∼ 30 (Shull et al. 2012). While these two bounds on the reionization history are well understood, we have yet to uncover the complete history. The observations that shed light on this history can be divided into two classes, direct and indirect probes. As a direct detection, spectroscopy of high redshift quasars or galaxies (Mortlock et al. 2011; Schenker et al. 2014) reveals the absorption by neutral hydrogen, which reveals its fraction as a function of redshift. Spectroscopy of high redshift quasars (Becker et al. 2001) shows evidence of IGM transition between its neutral and ionized state around redshift (z) 6. 21 cm observation (Pritchard & Loeb 2012; Loeb & Furlanetto 2013; Bowman et al. 2018), which is beyond the scope of this paper, and is the next promising avenue to constrain the epoch of reionization. Observation of luminosity from possible sources of reionization gives an indirect estimation of the reionization history. However, high energetic and fainter sources of reionization are difficult to observe. Until now we have obtained observations of luminosity from sources up to z ∼ 10. We have also been able to estimate the neutral hydrogen fractions until z ∼ 8, albeit with large uncertainties. Scattering of the Cosmic Microwave Background (CMB) photons by the electrons freed during the ionization process provides an alternate window to probe its history. The E-mode auto-correlation spectrum here constrains the integrated free electron fraction along the line of sight. Using these direct and indirect windows of observation, several efforts of constraining reionization history have been made (Hu & Holder 2003; Mitra et al. 2013; Douspis et al. 2015; Price et al. 2016; Hazra & Smoot 2017; Heinrich et al. 2017; Hazra et al. 2020; Paoletti et al. 2020; Qin et al. 2020; Chatterjee et al. 2021) in the last two decades. While comparing with data, the theory of ionization or the free electron fraction is modeled with parametric and nonparametric methods. CMB data releases have been providing crucial information about the line-of-sight integral of free electron fraction. With the Planck 2018 release the most precise information has been obtained about the optical depth with their low multipole E-mode polarization data from the high frequency instrument. In the last decade, in parallel, the detection of sources reached a milestone with Hubble Frontiers Fields (HFF). Here we make a joint analysis of these data in a new model-independent framework. Note that several simpler and popular parameterizations of reionization history in CMB also exist where the free electron fraction is parameterized (Lewis 2008; Hazra & Smoot 2017; Heinrich et al. 2017; Millea & Bouchet 2018). However, since these methods do not have source function parameterization, they cannot be used with the observation of UV luminosities and, in certain cases, they lead to unphysical histories of reionization.

In this paper, we explore two important aspects of reionization history. Using the model-dependent parametric formalism, we address the UV luminosity density data estimated from the observation of HFF with two different truncation magnitudes. With Gaussian process regression, we make model-independent reconstruction of sources and compare the reconstruction with parametric models. With this model-independent reconstruction of sources, we then reconstruct the reionization history by solving the ionization equations, using CMB and neutral hydrogen data along with the HFF data. The paper is organized as follows—in Section 2 we discuss the methodology of the parametric model and nonparametric Gaussian process regression. Following that we discuss the data sets used in this analysis in Section 3. The results are presented in Section 4. Finally, we conclude in Section 5.

The ionization equation describes the time evolution of the volume filling factor of ionized hydrogen in the intergalactic medium, ${Q}_{{{\rm{H}}}_{\mathrm{II}}}$, by a first-order ordinary differential equation (Madau et al. 1999; Wyithe & Loeb 2003):

$$\begin{eqnarray}&&\displaystyle \frac{{{dQ}}_{{{\rm{H}}}_{\mathrm{II}}}}{{dt}}=\displaystyle \frac{{\dot{n}}_{\mathrm{ion}}}{\langle {n}_{{\rm{H}}}\rangle }-\displaystyle \frac{{Q}_{{{\rm{H}}}_{\mathrm{II}}}}{{t}_{\mathrm{rec}}}.\end{eqnarray} \tag{ 1 }$$

The source term is characterized by the rate of production of ionizing photons ${\dot{n}}_{\mathrm{ion}}$, which depends on the UV luminosity density function ρ_UV, the efficiency of the source to produce ionizing photons ξ_ion, and the fraction of photons that actually escape into the IGM f_esc, and is defined as ${\dot{n}}_{\mathrm{ion}}={\rho }_{\mathrm{UV}}\langle {f}_{\mathrm{esc}}{\xi }_{\mathrm{ion}}\rangle$, where 〈f_esc ξ_ion〉 is a magnitude-averaged product. The sink term in the ionization equation accounts for the recombination process in the IGM. The recombination time t_rec is determined by the recombination coefficient α_B (T) and the clumping factor ${C}_{{{\rm{H}}}_{\mathrm{II}}};$ ${t}_{\mathrm{rec}}={\left[{C}_{{{\rm{H}}}_{\mathrm{II}}}{\alpha }_{B}(T)\left(1+\tfrac{{Y}_{p}}{4{X}_{p}}\right)\langle {n}_{{\rm{H}}}\rangle {\left(1+z\right)}^{3}\right]}^{-1}$. We keep the IGM temperature T fixed at 20,000 K in our analysis. X_p and Y_p are the primordial mass fractions of hydrogen and helium, respectively. n_H, n_He, and ${n}_{{{\rm{H}}}_{{\rm{II}}}}$ are the number densities of hydrogen, helium, and ionized hydrogen, respectively.

The evolution of the UV luminosity density with redshift can be obtained by parametric and nonparametric methods. In this article, as a parametric method we assume the density to be described by single power-law (Yu et al. 2012; Bouwens 2016) and double power-law (Ishigaki et al. 2015, 2018) forms,

$$\begin{eqnarray}&&{\rho }_{\mathrm{UV}}(z)={\rho }_{{\mathrm{UV}}_{1}}\cdot {10}^{-{az}}\,\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{\rho }_{\mathrm{UV}}(z)=\displaystyle \frac{2{\rho }_{\mathrm{UV},z={z}_{1}}}{{10}^{a(z-{z}_{1})}+{10}^{b(z-{z}_{1})}}.\end{eqnarray} \tag{ 3 }$$

While the first parameterization can be characterized completely by two parameters, namely the tilt (a) and the amplitude (${\rho }_{{\mathrm{UV}}_{1}}$), the second needs four parameters: the amplitude (${\rho }_{\mathrm{UV},z={z}_{1}}$), two tilts (a, b), and the redshift (z₁) at which the tilt in the power changes.

While parametric methods are useful, the functional form restricts their ability to address the data in several instances. Therefore model-independent reconstructions like Mason et al. (2019a) and Hazra et al. (2020) provide a conservative analysis owing to their flexibility. In this article, for a Bayesian, nonparametric reconstruction of luminosity density we make use of Gaussian process regression. A Gaussian process (GP) is a collection of random variables such that the joint distribution of any finite subset of it is described by a multivariate Gaussian (Rasmussen & Williams 2006). A GP is characterized by its mean function μ(x) and covariance function $k(x,x^{\prime} )$, where for a real process f(x), we have $\mu (x)={\mathbb{E}}[f(x)]$ and $k(x,x^{\prime} )={\mathbb{E}}[(f(x)-\mu (x))(f(x^{\prime} )-\mu (x^{\prime} ))]$. For a finite set of training points x = {x_i }, a function f(x) evaluated at each x_i can be represented by a random variable with a Gaussian distribution, such that the vector f = {f_i } has a multivariate Gaussian distribution given as $f\sim { \mathcal N }\left(\mu (x),C(x,x)\right)$, where C is the covariance matrix characterized by the covariance function k, which gives the covariance between two random variables. For the purpose of our analysis, we have chosen the covariance function to be the radial basis function represented as $k({x}_{i},{x}_{j})={\sigma }_{f}\exp \left(-\tfrac{{\left({x}_{i}-{x}_{j}\right)}^{2}}{2{l}^{2}}\right)$, with the correlation length l and amplitude σ_f .

In our analysis, we define the UV luminosity density in a few redshift nodes, and the values of the UV luminosity density at these nodes are taken as free parameters. We use the sklearn (Pedregosa et al. 2011) package in Python to implement Gaussian process regression. We have also developed independent Fortran codes that are used to obtain the bounds on hyperparameters. Gaussian processes have been useful in the reconstruction of several estimators related to cosmological processes such as dark energy (Holsclaw et al. 2010; Seikel et al. 2012; Shafieloo et al. 2012, 2013), test of the standard model (Aghamousa et al. 2017) with CMB, etc. While we introduce the use of GP in obtaining reionization history by reconstructing the UV luminosity densities, for recombination time, we use the aforementioned analytical formalism as has been demonstrated in Hazra et al. (2020), Mason et al. (2019a), and Paoletti et al. (2021) that apart from strict monotonic histories, t_rec cannot be constrained due to degeneracies with source terms.

Using the solutions of the ionization Equation (1), the Thomson scattering optical depth is defined as:

$$\begin{eqnarray}&&\tau =\int \displaystyle \frac{c{\left(1+z\right)}^{2}}{H(z)}{Q}_{{{\rm{H}}}_{\mathrm{II}}}(z)\langle {n}_{{\rm{H}}}\rangle {\sigma }_{{\rm{T}}}\left(1+\eta \displaystyle \frac{{Y}_{p}}{4{X}_{p}}\right){dz},\end{eqnarray} \tag{ 4 }$$

where helium is assumed to be singly ionized for z > 4 (η = 1) and doubly ionized for z < 4 (η = 2) (Kuhlen & Faucher-Giguère 2012).

We have used data sets from observations of three different types. First, we use the optical depth constraints from the Planck 2018 release of cosmic microwave background observation. With the low multipole polarization likelihood using data from high frequency instruments, the best CMB constraint on the optical depth is obtained as τ = 0.054 ± 0.007 (Aghanim et al. 2020). Since it has been demonstrated that the degeneracy between the primordial spectral tilt and the optical depth has reduced significantly in the new CMB data (Hazra et al. 2020; Paoletti et al. 2020), instead of using the full data we use this constraint on optical depth as the summary statistics. Second we use the derived UV luminosity density data (Bouwens et al. 2015; Ishigaki et al. 2018) in the compilation used in the Ishigaki et al. (2018) analysis, from HFF observations (Lotz et al. 2017, http://www.stsci.edu/hst/campaigns/frontier-fields/). Here we use the luminosity density measurements with truncation magnitude M_trunc of −17 and −15. Finally we use the measurements of neutral hydrogen fractions from the Lyα emission from galaxies (Ono et al. 2012; Schenker et al. 2014; Tilvi et al. 2014; Mason et al. 2019b), damping wings of gamma-ray bursts (Totani et al. 2006; McQuinn et al. 2008) and quasar spectra (Schroeder et al. 2013; Greig et al. 2017; Davies et al. 2018), and the dark gap in quasar spectra (McGreer et al. 2015) ionized zones near high redshift quasars (Bolton et al. 2011; Mortlock et al. 2011).

4.1. Can a Power Law Explain the UV Luminosity Density Data?

We perform a Gaussian process regression (GPR) for the UV luminosity density data for both UV17(M_trunc = −17) and UV15(M_trunc = −15). In Figure 1 we plot the log of GP likelihood as a function of the hyperparameters. We present the results for three different choices of the mean function: constant UV luminosity density (top), single power-law model (middle) and logarithmic double power-law model (bottom). The mean functions are chosen as the best-fitted function to the corresponding data sets. 1σ–3σ confidence contours are plotted in each case. If amplitude parameter σ_f is well-bounded for a particular mean function, it implies the rejection of that mean function with corresponding significance. As expected, constant mean functions are rejected with high significance. The single power-law model is rejected by the UV17 data at more than 2σ C.L. The UV15 data, however, does not rule out the single power-law model for higher uncertainties and larger values of UV luminosity density at high redshifts (this does not require a drop in power after a particular redshift). The logarithmic double power-law model is completely consistent with both the data sets.

Zoom In Zoom Out Reset image size

We now compare the power-law best fits within the uncertainty bands of this reconstruction as a function of redshifts. Figure 2 shows the model-independent reconstruction of the UV luminosity density using GPR. The training points here are the ${\mathrm{log}}_{10}{\rho }_{\mathrm{UV}}$ data from UV15 and UV17, which are used to train the kernel hyperparameters. Then, the optimized hyperparameters are used to interpolate between the training points. The top panel shows the results for UV17 and the bottom panel for UV15. The plots on the left show the interpolation of ${\mathrm{log}}_{10}{\rho }_{\mathrm{UV}}$ obtained using GPR, overlaid with the corresponding data points used as training points, along with the power-law best fits.

Zoom In Zoom Out Reset image size

We perform χ² minimization to obtain the best-fit parameters for the two power-law (single and logarithmic double power law) forms. The plots on the right show the residuals of power-law, GPR construction, and data points with respect to the single power law. Therefore, the zero line represents the single power law. In the top panel, for UV17, we notice that for higher redshifts the single power law does not fall within the 95% confidence interval of the model-independent reconstruction and therefore it is ruled out (as expected, this is consistent with the results from the hyperparameter contours in Figure 1). The double power law seems to be agreeing with the reconstruction at all redshifts in this case. In the case of UV15 (bottom panel), it can be observed that power-law models agree with the model-independent reconstruction with small deviations. In our analysis with constant, single power law and logarithmic double power law as mean functions, we find that the double power law is completely consistent with both UV15 and UV17 data. While the single power law is consistent with UV15 data, it is ruled out at more than 95% C.L. by UV17 data.

4.2. Constraints on Reionization History

To solve the ionization equation (Equation (1)), instead of assuming a specific function for the UV luminosity density, we follow a model-independent approach. Compared to earlier model-independent approaches in Hazra et al. (2020) and Mason et al. (2019a), we use Gaussian process regression to construct the UV luminosity density at different redshifts and find constraints on reionization history by solving the ionization equation using this model-independent form.

We follow a similar approach to Gerardi et al. (2019). Here we use four equidistant nodes between redshifts 4–10 to define the UV luminosity densities. The redshift range is chosen in such a way that it falls completely inside the range of redshifts of available data. It should be noted that to capture this correlation the distances between the nodes have to be smaller than the correlation length and we have tested the robustness of the analysis with varying number of nodes. The values of UV luminosity density at the redshift nodes are taken as free parameters for Markov Chain Monte Carlo (MCMC) sampling (using emcee Foreman-Mackey et al. 2013), and at each step these points are used as training points for Gaussian process regression. Therefore, here the Gaussian process reconstruction provides the samples of the history of UV luminosity densities for different choices of training points. Using the source function obtained from the sample of UV luminosity density evolution, we then solve Equation (1) to get the reionization history. For the integration we need ρ_UV(z) to extend beyond the available data. Since outside the range of data, GP reconstruction simply merges with the mean function, the choice of mean function is important here. We have used the logarithmic double power law as the mean function for the GP here as both data sets agree with this functional form. While for UV15 the single power law would also act as an appropriate choice of mean function, to keep the uniformity, here too we use the double power-law mean function that includes the single power law within its parameters.

We now obtain joint constraints on reionization history using all three data sets described earlier: Planck optical depth, ${Q}_{{{\rm{H}}}_{\mathrm{II}}}$ data, and UV luminosity density data for M_trunc = −17 and −15. In the ionization equation, we vary 〈f_esc ξ_ion〉 as one parameter by absorbing f_esc into the parameter ξ_ion (see the review by Dayal & Ferrara 2018 for bounds on the escape fraction). We keep a uniform prior on ${\mathrm{log}}_{10}{\xi }_{\mathrm{ion}}$ between 23.5 and 27.5 $[{\mathrm{log}}_{10}$ Hz erg⁻¹] as suggested in Price et al. (2016). We perform analyses with the clumping factor fixed as well as free with a uniform prior between 1 and 10. The other free parameters in the analysis are the values of UV luminosity density at the four redshift nodes: ρ₁₋₄. For the case with fixed ${C}_{{{\rm{H}}}_{\mathrm{II}}}$, the left panel of Figure 3 shows the MCMC constraints obtained for the five free parameters using both combinations of data sets: CMB+QHII+UV17 (in blue) and CMB+QHII+UV15 (in red). The constraints on UV luminosity densities at different redshift nodes are significantly different in these two data combinations, and since UV15 integrates luminosity function until larger magnitudes, we find constraints on the density to be shifted at larger values consistently at all redshift nodes. At higher redshifts, the constraints degrade as the data has larger uncertainties. Since the source term depends on both the UV luminosity density and the ionizing photon production efficiency, they must be anticorrelated in the CMB+UV+QHII data combination. This is easy to appreciate, as higher UV luminosity density and higher efficiency will lead to a larger optical depth, which will be ruled out by CMB optical depth constraints. We notice this anticorrelation in the last two redshift nodes as they mainly fall within the epoch of reionization. As expected, we also notice that the marginalized posteriors of ${\mathrm{log}}_{10}{\xi }_{\mathrm{ion}}$ are shifted in opposite directions compared to the UV luminosity density posterior distributions.

Zoom In Zoom Out Reset image size

The resulting constraints on the redshift of reionization z_re, duration of reionization Δz, and the reionization optical depth τ are listed in Table 1. The optical depth constraints obtained in both cases agree with 1σ optical depth from Planck results. We find that the optical depth is larger in the analysis where UV15 is used compared to UV17 data. Our constraint on the duration of reionization (Δz) indicates that the significant part of the reionization (10%–90% ionization) has occurred in ∼2 (for CMB+UV17+QHII) and ∼3 redshifts (for CMB+UV15+QHII). Since UV15 data show flattening of the luminosity densities at higher redshifts, it allows for a longer and earlier history of reionization compared to the data combination where UV17 data is used. The 68% and 95% bounds indicate that the marginalized posterior of Δz is significantly skewed. The redshift where the reionization is 50% completed is denoted as the redshift of reionization (z_re) and for both data combinations we find ${z}_{\mathrm{re}}\sim 7$.

Table 1. Joint Constraints (with ${C}_{{{\rm{H}}}_{\mathrm{II}}}$ Fixed) on the Reionization Optical Depth τ, Duration of Reionization Δz, and the Redshift of Reionization z_re with 68% and 95% Intervals

	CMB+UV17+QHII	CMB+UV15+QHII
τ	0.052 ± 0.001 ± 0.002	${0.056}_{-0.004-0.005}^{+0.002+0.006}$
Δz	${2.05}_{-0.21-0.30}^{+0.11+0.37}$	${2.98}_{-1.0-1.1}^{+0.30+1.8}$
zre	${6.93}_{-0.04-0.13}^{+0.08+0.11}$	${7.00}_{-0.10-0.22}^{+0.11+0.22}$

Download table as:  ASCIITypeset image

For the case with ${C}_{{{\rm{H}}}_{\mathrm{II}}}$ assumed as a free parameter, we plot the correlation between the ${\mathrm{log}}_{10}{\xi }_{\mathrm{ion}}$ and ${C}_{{{\rm{H}}}_{\mathrm{II}}}$ in the right panel of Figure 3. We find that the efficiency and the clumping factors are positively correlated and both cannot be constrained with these data combinations. This degeneracy is also discussed in Gorce et al. (2018), Mason et al. (2019a), and Paoletti et al. (2021). This indicates anticorrelation with the recombination time. A decrease in the ionizing photon production efficiency decreases the source term, which is only compensated by increased time for recombination to sustain ionization.

From the samples obtained in CMB+UV17+QHII and CMB+UV15+QHII data analysis (with fixed ${C}_{{{\rm{H}}}_{\mathrm{II}}}$), the reconstructed histories of reionization are plotted in the top and bottom panels, respectively, of Figure 4. We use the fgivenx (Handley 2018) package to compute the confidence band on the reconstructed ionization history. Different confidence limits are indicated in the color bar on the right. On top of the band, we have plotted the QHII constraints we have used in our analysis. Note that when UV17 data is used in the combination, the Gaussian process reconstruction does not allow early onset of ionization as it constrains the ${Q}_{{{\rm{H}}}_{\mathrm{II}}}\lt 0.1$ at z ∼ 10 with 95% C.L. When UV15 data is used we find the constraints are much more relaxed, with ${Q}_{{{\rm{H}}}_{\mathrm{II}}}\lt 0.3$ at the same redshift. Here we also find a longer tail of reionization between z = 8 and 15. It is important to note that with the present observation at high redshifts the constraints cannot be made tighter unless we restrict our model to a specific redshift symmetric form.

Zoom In Zoom Out Reset image size

In this paper, using Gaussian process regression as a nonparametric method for model-independent reconstruction of the evolution of UV luminosity density we address two questions. First, within the confidence band of reconstruction, we show that while the commonly used logarithmic double power-law model of UV luminosity density evolution agrees well with the data, the single power-law is ruled out. Second, using the reconstructed UV luminosity density evolution we reconstruct the reionization history using the optical depth from the cosmic microwave background observation, UV luminosity data from Hubble Frontiers Field observation, and neutral hydrogen fraction data from galaxy, quasar, and gamma-ray burst observations. Gaussian process regression allows a free form reconstruction of reionization history and thereby provides the platform for a conservative analysis. The use of three types of data sets, even in this conservative framework, breaks degeneracies between parameters and provides stringent constraints. From the joint constraints (when UV luminosity data with truncation magnitude 17 is used in the combination) we highlight that at 68%–95% C.L. we find the optical depth τ = 0.052 ± 0.001 ± 0.002 and the duration between 10% and 90% ionization is ${2.05}_{-0.21-0.30}^{+0.11+0.37}$ in redshifts. When the joint analysis uses UV luminosity data assuming truncation magnitude of −15, a longer tail in the ionization history is observed between redshift 10–14 with less than 30% ionization fraction (at 95% C.L.) at redshifts higher than 10.

Since we find that the use of UV luminosity density obtained by assuming higher (−15) truncation magnitude results in higher optical depth and larger duration of ionization, high redshift observations with JWST (Gardner et al. 2006; Madau 2017) and THESEUS (Amati et al. 2021) will definitely be helpful in providing better constraints, particularly around the tail. Cosmic variance limited CMB polarization data (Hazumi et al. 2020) will also help in constraining the optical depth better and therefore can also indirectly improve the joint constraints.

We would like to acknowledge the use of the HPC system Nandadevi (hpc.imsc.res.in) at the Institute of Mathematical Sciences, Chennai. We would like to thank Abraham Loeb, Sourav Mitra, Daniela Paoletti, Tirthankar Roy Choudhury, and Arman Shafieloo for their comments on the manuscript. D.K.H. thanks Masami Ouchi and Masafumi Ishigaki for providing their UV luminosity density data compilation.