Constraints on Reionization from Planck Low-$\ell$ EE Polarization and the Lyman-Series Forests

Introduction

Cosmic reionization marks the transition of the intergalactic medium (IGM) from a predominantly neutral state into the highly ionized medium observed at later times. The reionization history encodes information about the formation and evolution of the first ionizing sources, the escape of ionizing radiation into the IGM, and the density structure of the Universe. Extracting this information, however, has proven extremely challenging; it requires combining complementary constraints from a broad range of observables.

The large-scale Cosmic Microwave Background (CMB) $E$-mode polarization signal provides an important integral constraint on the reionization history. Thomson scattering by free electrons generates the “reionization bump” in the low-$\ell$ $EE$ power spectrum [1, 2], while suppressing the primary temperature and polarization anisotropies at smaller angular scales by $e^{-2\tau_{\mathrm{es}}}$, where $\tau_{\mathrm{es}}$ is the reionization Thomson scattering optical depth. In many CMB analyses, reionization is represented by a rapid and symmetric transition characterized by a single free parameter: the midpoint of reionization. Within this restricted family of histories, $\tau_{\mathrm{es}}$ maps almost uniquely onto the reionization midpoint and is therefore commonly described as primarily a constraint on that midpoint.

Recent analyses have emphasized, however, that the full low-$\ell$ polarization spectrum contains information beyond a single parameter. This information may be extracted using physics-based reionization models, flexible parametric histories, redshift bins, or a principal-component framework [3, 4, 5, 6, 7, 8, 9, 10]. Current Planck data robustly constrain $\tau_{\mathrm{es}}$ across a range of such modeling approaches, but the inferred reionization history remains considerably more sensitive to the adopted model and its implicit priors [10]. Joint analyses with astronomical observations have nevertheless begun to place constraints on additional parameters such as the duration and asymmetry of reionization [e.g. 3, 7, 9]. These results motivate constraining reionization models directly with CMB likelihoods, rather than compressing the CMB information into a $\tau_{\rm es}$ constraint under a fixed reionization prescription.

Meanwhile, substantial progress has been made in constraining the last phases of reionization using the Lyman-series forests of high-redshift quasars. The rapidly evolving mean transmission, dark gap statistics, and large scatter in the Ly$\alpha$ effective optical depth together suggest that the IGM was undergoing a transition at $5 \lesssim z \lesssim6$ [11, 12, 13, 14, 15, 16]. Additional evidence comes from the presence of damping-wing-like features in forest spectra [17, 18, 19], as well as the rapid evolution in the mean free path of ionizing photons during this epoch [20, 21, 22, 23]. Forward modeling of these observations favors a late completion of reionization, with a residual tail extending below $z = 6$ [16, 21, 24, 25, 26, 27, 28, 29] Thus, the forest should no longer be regarded simply as requiring that reionization be complete by approximately $z = 6$. Instead, it provides quantitative information about the final $\sim10$–$20\%$ of the transition.

The complementarity between the CMB and the Lyman-series forests has important implications for the interpretation of $\tau_{\rm es}$ measurements. Ref. [30] pointed out that once the endpoint of reionization is constrained to the relatively narrow interval preferred by forest data, $\tau_{\rm es}$ becomes sensitive to the duration of reionization as well as its midpoint. For an unconstrained endpoint, an increase in $\tau_{\rm es}$ can generally be produced by shifting a short reionization history to higher redshift. If the tail end of reionization is instead restricted by the forests, an earlier midpoint generally implies a more extended history and therefore a longer duration. In this sense, the forest changes the degeneracies and physical information carried by the CMB optical depth.

The analysis of [30] demonstrated this effect by examining the dependence of $\tau_{\rm es}$ on the midpoint, duration, and asymmetry of parametric reionization histories subject to an approximate forest constraint on the endpoint. It did not, however, perform a joint Bayesian analysis of the reionization and cosmological parameters using the full CMB temperature and polarization likelihoods. Such an analysis is important for at least two reasons. First, different reionization histories produce different low-$\ell$ EE spectra even when they have similar values of $\tau_{\rm es}$. Second, the high-$\ell$ CMB spectra constrain the combination $A_{\rm s}e^{-2\tau_{\rm es}}$, where $A_{\rm s}$ is the amplitude of primordial scalar fluctuations, so that the reionization parameters must be inferred jointly with $A_{\rm s}$ and the other cosmological parameters.

In this work, we carry out a joint inference using the full Planck temperature and polarization likelihoods. Our central goal is to determine how the different Lyman-series forest constraints alter reionization posteriors inferred from the CMB. To this end, we explore the use of different reionization models including the parametric model of [31] (used in [30]) and the non-parametric FlexKnot approach [32]. We will also apply the recently developed physical model of Ref. [33] which connects the ionizing source emissivity to the global ionized fraction and IGM photoionization rate self-consistently. This model was developed specifically for inference pipelines utilizing Lyman-series forest measurements. Here we will apply it for the first time. Several distinct forest observables can be used to extract information on the reionization history and they differ substantially in both constraining power and model dependence. For instance, the fraction of dark pixels in the Lyman-series forests provide a conservative and nearly model-independent upper limit on the volume-averaged neutral fraction [34, 35, 36]. On the other hand, a conservative but simulation-calibrated lower-limit on the neutral fraction at $z \approx5.8$ has been derived from Ly$\alpha$ damping-wing-like features on the boundaries of stacked forest dark gaps [18]. Among the most constraining statistics is the redshift-evolution of the mean forest transmission [12, 13]. Applying it, however, requires detailed forward-modeling. Another goal of this paper is to explore the information content within this range of forest statistics, and how different combinations alter the joint posteriors. Our analysis illustrates how information about reionization’s end complements the CMB polarization and temperature in constraining the reionization history.

This paper is organized as follows. In § 2, we describe our Bayesian inference method, data sets, and likelihoods. In § 3 we present the results of our joint inferences using CMB, dark pixel, and damping-wing-like forest measurements. We present posteriors for both parametric and non-parametric models. In § 4 we discuss how these constraints can be tightened considerably by implementing measurements of the forest transmission fraction. We conclude in § 5. Unless otherwise stated, all distances are reported in comoving units.

Methods

Bayesian Parameter Inference

We use the publicly available code, Cobaya, to perform Bayesian parameter inference [37]. Cobaya is widely used for obtaining constraints from cosmological data sets, especially CMB measurements. Conveniently, it is compatible with Boltzmann theory codes such as Code for Anisotropies in the Microwave Background (CAMB; [38]) and Cosmic Linear Anisotropy Solving System (CLASS; [39]). It also comes with a number of CMB likelihood functions, e.g. from Planck. For this work, we make two modifications to the Cobaya code base. First, we use a modified version of CAMB that accepts arbitrary reionization histories in the form of a numeric table of redshifts and ionized fraction values, $x_{\mathrm{HII}}$. This affords us flexibility in the range of models that we consider. Second, we construct new likelihoods for the forest constraints used in this work, which we describe below.

We sample the model parameters using Cobaya’s mcmc sampler and pass them to our external forest-likelihood functions. These functions calculate the reionization history from the parameter values and evaluate the log-likelihood of that history for the forest constraints considered. In what follows, we explore three different methods for modeling the reionization history: a parametric model (§ 3.1) in which a functional form is assumed, the more free-form FlexKnot approach of Ref. [32] (§ 3.2), and the physical model of Ref. [33] that connects the ionizing emissivity of the source population to the global ionized fraction and photoionization rate of intergalactic hydrogen (§ 4). The details of these methods will be described in their respective sections. In addition to the log-likelihood values, our external functions also return $x_{\mathrm{HII}}$ evaluated at 50 linearly-spaced points in the range $4 < z < 30$, which are then passed to CAMB to calculate the CMB power spectra and evaluate the CMB likelihoods.

Unless otherwise noted, we adopt the same uniform priors for cosmological parameters as in Ref. [40]: $\theta_{\mathrm{MC}} \in[0.0103, 0.0105]$, $\ln(10^{10} A_s) \in[2.6, 3.5]$, $n_s \in[0.9, 1.1]$, $\Omega_b h^2 \in[0.017, 0.027]$, $\Omega_c h^2 \in[0.09, 0.15]$. We also impose a conservative lower limit of $z = 4$ for the end of reionization. We use Cobaya’s generalized Gelman-Rubin statistic ($R - 1$) to evaluate the convergence of our mcmc chains. We consider our chains converged when $0.01 \lesssim R - 1 < 0.03$, where the upper limit of this range represents the most computationally intensive runs (generally the FlexKnot runs). To analyze the samples we use the GetDist code package [41].

Data and Likelihoods

CMB

We combine the following likelihoods in order to leverage the full range of scales measured by Planck:

• low-$\ell$ EE: We use the latest version of the Low-$\ell$ Likelihood on Polarised Power-spectra (LoLLiPoP) of Refs. [42, 43], based on Planck PR4 and covering the multipole range $2 \le\ell\le29$.

• Low-$\ell$ TT: We use the Planck PR3 Commander likelihood in the multipole range $2 \le\ell\le29$ [44].

• high-$\ell$ TTTEEE: We use the nuisance-marginalized Planck PR3 Pliklite likelihood over $30 \le\ell\le2{,}508$ for TT, and $30 \le\ell\le1{,}996$ for TE and EE [44]. We choose this over the full (unmarginalized) Plik likelihood and the more recent HiLLiPoP for computational efficiency, as the small-scale CMB anisotropies are relatively less sensitive to the details of the reionization history.

The above combination mixes PR3 and PR4 likelihoods, pairing the PR4 low-$\ell$ EE likelihood with the PR3 high-$\ell$ TTTEEE likelihood (the PR4 analysis retained the PR3 Commander low-$\ell$ TT likelihood, since little improvement in large-scale temperature was expected). In Appendix § 7, we present tests verifying that this choice does not affect our main results. In summary, we ran an analysis where we replaced Pliklite with HiLLiPoP for the high-$\ell$ TTTEEE likelihood and found little difference in the inferred reionization histories.

Lyman Series Forest

We consider the combination of CMB measurements with three types of Lyman-series forest measurements: (1) upper limits on the volume-weighted mean neutral fraction, $\langle x_{\mathrm{HI}}\rangle_{\mathrm{V}}$, derived from dark pixels in the $z \approx4.9 - 6.7$ Ly$\alpha$, Ly$\beta$, and Ly$\gamma$ forests; (2) damping-wing-like features in the $z \approx5.5 - 6$ Ly$\alpha$ forest, which place a complementary lower-limit on $\langle x_{\mathrm{HI}}\rangle_{\mathrm{V}}$; (3) Ly$\alpha$ forest measurements of the ionizing radiation background strength at $4.9 \le z \le6$, as quantified by the hydrogen photoionization rate, $\Gamma_{\mathrm{HI}}$. The methodology for (3) is more detailed and the results that we obtain from it involve several important caveats. As such, we defer a description of those measurements to a dedicated section (§ 4). In this section, we describe the dark pixel and damping-wing-like measurements that form the backbone of our forest likelihoods.

The dark-pixel constraints on $\langle x_{\mathrm{HI}}\rangle_{\mathrm{V}}$ exploit the extreme opacity of the high-redshift IGM to Lyman-series absorption. At $z > 5$, a local residual neutral fraction of order $x_{\mathrm{HI}} \sim10^{-4}$ is enough to render the Ly$\alpha$ forest effectively dark; Ly$\beta$ absorption is also generally strongly saturated at this level. Thus any detected transmission identifies a highly ionized region, whereas a dark region may contain either neutral IGM or ionized gas with a sufficient level of residual neutral hydrogen. The dark-pixel fraction therefore provides a conservative upper limit on the volume-averaged neutral fraction. While this approach is generally model-independent, in practice one must still make choices on how to process the forest spectra and identify dark pixels. Following Ref. [34], the forest spectra are typically binned into 3.3 cMpc pixels.¹ Two methods have been employed to identify the dark pixels. The more conservative “threshold” approach simply flags as dark any pixel below some signal-to-noise ratio. The second approach, termed the negative-pixel method, assumes that truly dark pixels scatter symmetrically above or below the zero point, so the number of dark pixels can be approximated as twice the number of negative pixels. The negative pixel method provides tighter constraints on $\langle x_{\mathrm{HI}}\rangle_{\mathrm{V}}$ but is more vulnerable to biasing from imperfect subtraction of the sky signal. In what follows, we use limits obtained from the threshold method.

We use upper-limits from two sources. First, Ref. [36] provides the most up-to-date limits at $z = 4.85 - 6.25$ from 34 deep quasar spectra from the extended XQR-30 program [45]. They report multiple sets of limits by considering different combinations of Ly$\alpha$, Ly$\beta$ and Ly$\gamma$ forest measurements. We take the combination that yields their most stringent $\langle x_{\mathrm{HI}}\rangle_{\mathrm{V}}$ upper limits at each redshift bin (the bolded values in the 7th column of their Table 2).² To extend the redshift range of the dark pixel constraints, we also include the $z = 6.3$, $6.5$, and $6.7$ upper limits from Ref. [35], who similarly use dark pixels in the Ly$\alpha$ and Ly$\beta$ forests. The combination of Refs. [35, 36] represent the most up-to-date and model-independent upper limits on the volume-weighted mean neutral fraction of the IGM at $4.9 \lesssim z \lesssim6.7$. We approximate the likelihood at each redshift bin as one-sided Gaussian distributions. For a reported $\langle x_{\mathrm{HI}}\rangle_{\mathrm{V}}(z_i)$ upper limit at redshift $z_i$, if the predicted neutral fraction, $x_{\mathrm{HI}}(z_i) \ge\langle x_{\mathrm{HI}}\rangle_{\mathrm{V}}(z_i)$, we approximate the likelihood as a Gaussian centered at $\langle x_{\mathrm{HI}}\rangle_{\mathrm{V}}(z_i)$ with standard deviation equal to the reported $1\sigma$ uncertainty. On the other hand, if $x_{\mathrm{HI}}(z_i) \le\langle x_{\mathrm{HI}}\rangle_{\mathrm{V}}(z_i)$ then the likelihood returns a constant value chosen for continuity at $\langle x_{\mathrm{HI}}\rangle_{\mathrm{V}}(z_i)$. The total likelihood consists of the product of likelihoods at each redshift.

A complementary lower bound on $\langle x_{\mathrm{HI}}\rangle_{\mathrm{V}}$, at $z \simeq5.8$ has been inferred from a damping-wing-like absorption profile in the stacked Ly$\alpha$ forest adjacent to long Ly$\beta$ dark gaps.³ Following the method originally proposed by Ref. [46], Ref. [18] selected 24 Ly$\beta$ dark gaps with lengths $L_{\mathrm{gap}} \ge7\,h^{-1}\,\mathrm{cMpc}$, aligned the corresponding Ly$\alpha$ forest spectra at the red and blue edges of the gaps, and stacked the transmission (also see [19]). The resulting profile remains suppressed for approximately $1000\ \mathrm{km}\,\mathrm{s}^{-1}$ beyond the gap edges, substantially longer than expected for simulated dark gaps produced entirely by reionized gas, even in the presence of large ionizing background fluctuations. Applying the same stacking and fitting procedure to mock spectra from radiative transfer simulations, Ref. [18] found that only stacks with mean neutral fraction $x_{\mathrm{HI,gap}} \gtrsim0.4$ produced Ly$\alpha$ forest profiles consistent with the observed one within the simulated 95% interval. Combining this minimum neutral fraction with the measured fraction of the Ly$\alpha$ forest occupied by dark gaps longer than the mean length of their sample, $f_{L\ge15h^{-1}\mathrm{cMpc}} = 0.152 \pm0.098$, gives $\langle x_{\mathrm{HI}}\rangle_{\mathrm{V}} \gtrsim x_{\mathrm{HI}}^{\min} f_{L\ge15h^{-1}\mathrm{cMpc}} = 0.4(0.152 \pm0.098) = 0.061 \pm0.039$ at $z = 5.8$. This should be interpreted as a conservative, simulation-calibrated lower bound; it counts only the neutral hydrogen associated with the long-dark-gap population and neglects any neutral gas residing in shorter gaps or elsewhere in the IGM. Our likelihood for this limit is constructed similarly to the dark pixel constraints, but flipped to reflect it being a lower limit on $\langle x_{\mathrm{HI}}\rangle_{\mathrm{V}}$ rather than an upper limit.

In the ensuing sections, we explore four combinations of the observational constraints described above:

• CMB Only applies only the CMB likelihoods listed in 2.2.1.

• DP + CMB applies the dark pixel upper-limits on $\langle x_{\mathrm{HI}}\rangle_{\mathrm{V}}$ from Refs. [35, 36], together with the CMB likelihoods.

• DWlike + CMB includes the lower limit on $\langle x_{\mathrm{HI}}\rangle_{\mathrm{V}}$ derived from the damping-wing-like feature reported by Ref. [18] together with the CMB.

• DP + DWlike + CMB includes both the upper and lower limits on $\langle x_{\mathrm{HI}}\rangle_{\mathrm{V}}$ together with the CMB.

An important caveat is that the dark pixel and damping-wing-like measurements provide limits on the volume-weighted reionization history, whereas the CMB temperature and polarization anisotropies are sensitive to the mass-weighted history. Mapping these two quantities requires a model for the spatial structure of the reionization process (see e.g. [47]). For simplicity, we assume that the mass-weighted and volume-weighted ionized fractions are approximately the same in our joint CMB-forest analysis. In Appendix § 8, we use a radiative transfer simulation of reionization to test the effect of this assumption on our results. We find that correcting for the difference between mass-weighted and volume-weighted has a small effect on our constraints ($\lesssim5\%$) relative to the uncertainties in our derived reionization histories.

Results

Parametric Reionization Model

We begin by considering a parametric model of the reionization history. We adopt the three-parameter Weibull function previously introduced in Ref. [31]. This functional form has been used to obtain constraints from the CMB optical depth and patchy kinematic Sunyaev–Zel’dovich effect [30, 48, 49], and to forecast optical depth constraints from 21-cm clustering measurements [50]. To our knowledge, this is the first analysis in which the Weibull midpoint, duration, and asymmetry parameters are jointly constrained using the full Planck temperature and polarization likelihoods. The reionization history is set by reionization’s midpoint ($z_{50}$), duration ($\Delta z_{90}$) and asymmetry ($A_z$). The midpoint is defined to be the redshift where the mass-weighted ionized fraction of hydrogen reaches 50%. The duration and asymmetry parameters are

where $z_{05}$ and $z_{95}$ are the redshifts of 5% and 95% ionization, respectively. Note that $A_z = 1$ corresponds to a symmetric history, where the first and second halves of reionization take place over the same redshift interval, and $A_z > 1$ corresponds to a longer first half. The mass-weighted reionization history is given by

where we solve for $a_{\mathrm{w}}$, $b_{\mathrm{w}}$, and $c_{\mathrm{w}}$ for given $z_{50}$, $\Delta z_{90}$ and $A_z$ following the prescription in Appendix A of Ref. [31]. We adopt wide uniform priors of $z_{50} \in[5,25]$, $\Delta z_{90} \in[0.5,25]$, and $A_z \in[0.5,9]$. As pointed out by Ref. [32] and emphasized in Ref. [10], taking uniform priors for the parameters describing a reionization history does not in general yield a uniform prior for $\tau_{\mathrm{es}}$. This can be an important consideration for interpreting $\tau_{\mathrm{es}}$, especially when comparing between different models of the reionization history. Our choice of priors is motivated by the goal of this paper; to understand how the combination of forest and CMB observations constrain the reionization history itself.

In Figure 1, we show the posteriors for $z_{50}$, $\Delta z_{90}$ and $A_z$ for different combinations of forest constraints (colored and filled contours), against the baseline CMB Only analysis (black contours). The top-left and bottom-right sets of panels correspond to the DP + CMB and DWlike+CMB analyses, respectively. The bottom-left set corresponds to the full DP+DWlike+CMB analysis. The inner and outer contours correspond to the 68% and 95% credible regions, respectively.

Image not extracted

We begin with the CMB Only constraints (which are the same in all sets of panels). The marginalized posterior distribution for the reionization midpoint ($z_{50}$; top-left) peaks near $\sim7.3$, consistent with the constraints from previous works [51]. The middle-left panel shows an anti-correlation between $z_{50}$ and $\Delta z_{90}$, which arises from a degeneracy in $\tau_{\mathrm{es}}$. Longer durations (higher $\Delta z_{90}$) can be offset by later midpoints (lower $z_{50}$), and conversely, to maintain approximately the same $\tau_{\mathrm{es}}$. We also see that low $z_{50}$ broadly favors a larger asymmetry parameter, $A_z$. This is related to our requirement that reionization be completed no later than $z = 4$. For models with later midpoints and longer durations, a larger asymmetry is the only way to finish reionization on time.

Consider now the combination of DP+CMB constraints (top-left set of panels). The inclusion of dark pixel information disfavors $z_{50} \lesssim6.5$, a result of the upper limit $\langle x_{\mathrm{HI}}\rangle_{\mathrm{V}} \le39.6\%^{+9.3\%}_{-?}(1\sigma)$ at $\bar{z} = 6.043$ from Ref [36]. The middle panel also shows a tightening of the upper limit on $\Delta z_{90}$. This results from two effects: (1) As noted above, a longer duration is correlated with a later midpoint owing to their degeneracy in $\tau_{\mathrm{es}}$. Therefore, longer durations are disfavored when the $z \sim6$ dark pixel constraints disfavor the latest midpoints; (2) The dark pixel constraints require $\langle x_{\mathrm{HI}}\rangle_{\mathrm{V}} \lesssim13\%$ at $z \lesssim5.5$. All else being equal, this has the effect of requiring the parametric reionization histories to complete earlier than our $z = 4$ hard cut-off, resulting in shorter durations.

The bottom-right panels of Figure 1 show results from the DWLIKE+CMB combination. Recall that the lower limit on the neutral fraction from Ref. [18] complements the upper limits obtained from dark pixel measurements, so it is interesting to consider how this combination differs from the DP+CMB analysis. We see that the $z_{50}$ posterior is shifted downward and $\Delta z_{90}$ now peaks at $\Delta z_{90} \sim10$. These extended but delayed histories are able to satisfy the constraint from [18] that $\langle x_{\mathrm{HI}}\rangle_{\mathrm{V}}(z = 5.8) > 6.1\% \pm3.8\%$ while also satisfying the CMB limits on $\tau_{\mathrm{es}}$. We note, however, that the effect of the $x_{\mathrm{HI}}$ limit is relatively weak. The likelihood saturates at $x_{\mathrm{HI}} = 0$ at a modest $\sim1.6\sigma$. Hence, it does not rule out reionization histories that complete earlier. Keeping in mind the fixed-$\tau_{\mathrm{es}}$ degeneracy, this accounts for the low-$\Delta z_{90}$ tail seen in the middle panel.

Finally, the lower-left panels show results from the full combination of constraints, DP+DWLIKE+CMB. Here we see a steep decline in the posterior for $z_{50} \lesssim6.5$ from the dark pixel constraints, as well as a tilt towards lower-$z_{50}$ from the damping-wing-like constraint. The $\Delta z_{90}$ posterior exhibits a peak as in the DWLIKE+CMB case, but is somewhat weaker. We note that $A_z$ is poorly constrained (bottom-right). In the $A_z$–$\Delta z_{90}$ plane (bottom-middle), the posterior now exhibits a multi-peaked structure. Following the discussion above, this structure arises from the competing preferences of the dark-pixel and damping-wing-like constraints. In combination with the CMB, the former favors comparatively rapid completion of reionization, while the latter favors a more extended history with significant late-time neutral fraction. As a result, the DP+DWlike+CMB analysis yields disjoint parameter combinations yielding similar posterior probabilities. In $\Delta z_{90} - A_z$ space (bottom-middle), the peaks at $\Delta z_{90} \sim1$ correspond to models in which reionization ends before $z \approx6.4$, satu-

rating the weak lower limit on the neutral fraction from the damping-wing-like measurement at $z = 5.8$ (see Figure 2). This region of the parameter space would be ruled out if future observations are able to strengthen lower limits on the neutral fraction.

Image not extracted

In Table 1 we summarize our constraints on $z_{50}$ and $\Delta_{90}$. We report the marginalized posterior mean and 68% credible intervals on the former and 95% upper limits on the latter. In the full joint CMB-forest analysis, we find

Parametric (Weibull)	$z_{50}$	$\Delta z_{90}$	$A_z$	$\tau_{\rm es}$	$\tau_{\rm es}^{z>15}$
CMB Only	$7.3 \pm0.9$	$<13.5$	-	$0.0607 \pm0.0063$	$<0.0075$
DWlike+CMB	$6.9^{+0.7}_{-1.1}$	$<14.2$	-	$0.0602 \pm0.0061$	$<0.0084$
DP+CMB	$7.6^{+0.6}_{-0.8}$	$<10.6$	-	$0.0612^{+0.0060}_{-0.0065}$	$<0.0048$
DP+DWlike+CMB	$7.5^{+0.5}_{-0.9}$	$<11.5$	-	$0.0611 \pm0.0062$	$<0.0058$
FlexKnot
CMB Only	$7.6^{+0.9}_{-0.7}$	$<12.3$	-	$0.0608^{+0.0058}_{-0.0065}$	$<0.0044$
DWlike+CMB	$7.4 \pm0.9$	$<12.8$	-	$0.0600 \pm0.0062$	$<0.0048$
DP+CMB	$7.8 \pm0.7$	$<10.6$	-	$0.0612^{+0.0058}_{-0.0064}$	$<0.0032$
DP+DWlike+CMB	$7.7^{+0.7}_{-0.8}$	$<11.4$	-	$0.0610^{+0.0058}_{-0.0065}$	$<0.0039$
Semi-analytic (fixed cosmology and IGM opacity model)
CMB Only	$7.4 \pm0.3$	$<7.8$	-	$0.0564 \pm0.0024$	$<0.0011$
DP+DWlike+CMB	$7.4 \pm0.3$	$<7.6$	-	$0.0562 \pm0.0024$	$<0.0011$
$\Gamma_{\rm HI}$ + CMB	$7.4^{+0.2}_{-0.3}$	$6.0^{+0.8}_{-1.0}$	$2.4^{+0.3}_{-0.7}$	$0.0564 \pm0.0024$	$<0.0011$

The latter is complementary to recent constraints on $\Delta z_{90}$ from the patchy kinetic Sunyaev Zel’dovich (pkSZ) effect [52, 53, 54]. Our upper limit is weaker than the 95% limits of $\Delta z_{90} < 7.7$ obtained from ACT DR6 [52] and $\Delta z_{90} < 7.0$ from the recent SPT-3G analysis [53]. These rely on scaling relations calibrated using semi-numerical AMBER simulations [55]. Additionally, the reionization midpoint is fixed to $z_{50} = 8$ in the ACT analysis and $z_{50} = 7.68$ in the SPT-3G analysis, while the adopted AMBER scaling relation is calibrated at fixed asymmetry $A_z = 3$, rather than jointly marginalizing over $z_{50}$, $\Delta z_{90}$, and $A_z$.

We also tabulate our derived constraints on $\tau_{\mathrm{es}}$ and 95% upper limits on the $z > 15$ contribution to the optical depth, $\tau_{\mathrm{es}}^{z>15}$. We remind the reader of the implicit non-uniform prior for $\tau_{\mathrm{es}}$ that arises from adopting uniform priors on the Weibull function parameters. The inferred total optical depth is stable across the different likelihood combinations, with all central values agreeing to much better than $1\sigma$ and with essentially no change in the uncertainty. The forest constraints therefore primarily redistribute the allowed optical depth in redshift rather than altering its total value. The dark pixel data favor a slightly larger contribution from lower redshift while suppressing an extended high-redshift tail, whereas the damping-wing-like data permits a larger $z > 15$ contribution while shifting the total optical depth slightly downward. When both forest likelihoods are included, their opposing preferences largely cancel in $\tau_{\mathrm{es}}$, although the upper limit on $\tau_{\mathrm{es}}^{z>15}$ remains tighter than for the CMB alone.

FlexKnot Constraints

We now consider a more flexible model for the reionization history. In this section, we use the FlexKnot approach first applied in this context by Ref. [32]. In this approach, the reionization history is defined by a number of $(z_i, x_{\mathrm{HII}})$ pairs called knots. The reionization history between the knots is interpolated with piecewise cubic Hermite interpolating polynomials (PCHIP). We enforce monotonicity in the knots. In this case, the use of PCHIP ensures that the resultant reionization history will also be monotonic. Our fiducial FlexKnot model has four free parameters: $z_{\mathrm{beg}}$, $z_{\mathrm{end}}$, $z_1$ and $x_1$, where $z_{\mathrm{beg}}$ and $z_{\mathrm{end}}$ denote the redshifts at which $x_{\mathrm{HII}} = 0$ and $1$, respectively, and the $(z_1, x_1)$ pair determines the position of the free-floating knot between $z_{\mathrm{beg}}$ and $z_{\mathrm{end}}$.⁴ This model is similar to the 3-Flex model used in Ref. [10]. We adopt wide priors of $z_{\mathrm{end}}$, $z_1$, $z_{\mathrm{beg}} \in[4, 30]$, and $0 \leq x_1 \leq1$ but construct our likelihoods to ensure $z_{\mathrm{end}} < z_1 < z_{\mathrm{beg}}$. In Appendix § 7, we present our constraints on the cosmological parameters, which we find are stable across the different dataset combinations and reionization modeling methods applied in this paper.

Figure 2 shows the marginalized constraints on the reionization history for the different dataset combinations. From left to right, the columns correspond to CMB Only, DWlike+CMB, DP+CMB and DP+DWlike+CMB, respectively. The top row shows the entire reionization history, while the bottom row zooms in on the second half of reionization to highlight the impact of the forest constraints in that epoch. For reference, each panel also displays all of the dark pixel and damping-wing-like constraints used in our inference pipeline (although each column uses a different set of these). The inclusion of the damping-wing-like lower limit on $\langle x_{\mathrm{HI}}\rangle_{\mathrm{V}}$ has a modest effect on the reionization history contours, tightening the $1\sigma$ limits only slightly relative to the CMB Only case. On the other hand, the addition of the dark pixel constraints has a strong impact on the contours at $z \lesssim6.5$, requiring reionization to be more than $\sim80\%$ complete by $z \sim6$.

In Figure 3, we plot posteriors for the derived parameters $z_{50}$, $\Delta z_{90}$, and $A_z$. The middle section of Table 1 summarizes the constraints. Consistent with the results of the last section, the addition of the dark pixel constraints in DP+CMB disfavors $z_{50} \lesssim6.5$. This again leads to a preference for shorter $\Delta z_{90}$ relative to the CMB Only result in order to preserve a roughly constant $\tau_{\mathrm{es}}$. In the DWlike+CMB case, we also see a shift towards lower $z_{50}$ and longer $\Delta z_{90}$ relative to the CMB Only case. Comparing the posterior distributions of $z_{50}$, the main effects of the dark-pixel and damping-wing-like constraints are to pull the distributions toward earlier and latter midpoints, respectively. These constraints also pull the posterior of $\Delta z_{90}$ toward shorter and longer durations, respectively. When combined, the net result is a mild reduction in the widths of these distributions. In the full joint CMB-forest analysis, we find

Image not extracted

The posterior for $A_z$ exhibits two interesting features. First, the spike at $A_z < 1$ owes to the definition of $A_z$. When sampling over the uniform priors of the FlexKnot parameters, the samples are evenly split between $A_z > 1$ and $A_z < 1$. This means that half of the samples are constrained to $0 < A_z < 1$ and the other half to $A_z > 1$. Thus, the sample density with $A_z < 1$ is much higher, resulting in the observed spike. Secondly, the $A_z$ posterior is very different from the one obtained in the parametric analysis (compare to the bottom-right panels of each triangle plot if Figure 1). This owes to $A_z$ not being well constrained by the data and the different effective priors between the analyses. In the parametric analysis, we adopted flat priors on $z_{50}$, $\Delta z_{90}$, and $A_z$, while in the FlexKnot analysis we used flat priors for $z_{\mathrm{beg}}$, $z_{\mathrm{end}}$, $z_1$, and $x_1$.

In Figure 4, we explore the effects of implicit priors on our marginalized posteriors for $z_{50}$, $\Delta z_{90}$, and $A_z$ using importance sampling. We obtained the implicit prior distribution $\mathcal{P}(z_{50}, \Delta z_{90}, A_z)$, by sampling over the uniform prior distributions of $z_{\rm beg}$, $z_{\rm end}$, $z_1$, and $x_1$. We then assigned weights of $\mathcal{P}(z_{50}, \Delta z_{90}, A_z)^{-1}$ to each sample (see e.g. [10]). Figure 4 shows the resulting fully marginalized posteriors for the CMB Only (black) and DP+DWlike+CMB analyses. While there is little change in the $z_{50}$ posteriors, the $\Delta z_{50}$ and $A_z$ posteriors are qualitatively more similar to those obtained in the parametric model analyses (Figure 1). Notably, the peak in the $\Delta z_{90}$ posterior within the DP+DWlike+CMB analysis is significantly weaker, illustrating that the implicit prior on $\Delta z_{90}$ has an important effect on the shape of this distribution. The large difference in the $A_z$ posterior also confirms our assessment in the previous paragraph. After flattening the priors, we find $z_{50} = 7.4^{+1.1}_{-0.9}$, $\Delta z_{90} < 13.6$ for CMB Only, and $z_{50} = 7.5^{+0.6}_{-0.8}$, $\Delta z_{90} < 12.1$ for DP+DWlike+CMB.

Image not extracted

The top two sections of Table 1 show that the constraints on $\tau_{\rm es}$ are remarkably stable across not only the different combinations of data sets, but also across the parametric and FlexKnot modeling methods. The 95% limits on the $z > 15$ contribution to the optical depth, however, do depend on the choice of modeling method, with our implementation of FlexKnot yielding tighter upper limits than the parametric model. Again, at least part of these differences may owe to the different priors adopted between the two modeling methods. At a qualitative level, both analyses show that combining the CMB with information on the end of reionization tightens limits on the earliest stages of reionization.

Tightening Constraints with the Lyα Forest Transmission

In the preceding sections, we used dark-pixel and damping-wing-like absorption measurements to bracket the neutral fraction during the final stages of reionization. These relatively direct probes, however, do not exhaust the information contained in the Lyman-series forests. In particular, the redshift evolution and sightline-to-sightline fluctuations of the Lyα forest transmission are highly sensitive to the timing of reionization’s end. For example, forward-modeling the distribution of forest effective optical depths, $\tau_{\mathrm{eff}}$, can yield substantially tighter constraints on the reionization history [e.g., 29, 56]. Such an analysis requires a detailed model of the density, temperature, and ionizing background fluctuations that shape the forest, and lies beyond the scope of this paper. To illustrate the power of exploiting the additional information embedded within the forest, we instead use measurements of the hydrogen photoionization rate, $\Gamma_{\mathrm{HI}}$, which essentially encodes the evolution of the mean Lyα transmission. Although this approach does not capture the full constraining power of forest statistics, it demonstrates how strongly the mean transmission alone can sharpen the constraints. It also motivates the development of accurate and computationally efficient models for future analyses.

Methods

To connect the mean-transmission measurements to the reionization history, we use the semi-analytic model of Ref. [33] which is derived from the cosmological radiative transfer equation. Given a source ionizing emissivity history, the model jointly solves for the mass-weighted ionized fraction, $x_{\mathrm{HII}}(z)$, and $\Gamma_{\mathrm{HI}}(z)$, accounting for absorptions both within ionized regions and at the boundaries of the neutral IGM. In the model, these absorptions are accounted for with the total absorption coefficient such that the mean free path of ionizing photons can be written $\lambda^{-1} = \lambda_n^{-1} + \lambda_i^{-1}$. Here, $\lambda_n$ is the average distance ionizing photons travel to reach the edges of bubbles⁵, and $\lambda_i$ is the mean free path through ionized gas inside the bubbles. For the former, we measure the mean distance to neutral regions for randomly positioned skewers through the “Late start/Late end” reionization simulation of Ref. [57]. We find that this mean distance is well-fit by the following form

where $A = 0.278$ physical Mpc and $B = 2.060$. Absorptions within the reionized IGM are quantified by an effective clumping factor, termed the absorption rate clumping factor,

where $\dot{n}_{\mathrm{abs}}^{\mathrm{ionized}}$ is the ionizing photon absorption rate density in the gas, $\alpha_B(T_{\mathrm{ref}})$ is the case-B recombination coefficient evaluated at a reference temperature, $T_{\mathrm{ref}} = 10{,}000$ K, and the factor $1 + \chi= 1 + n_{\mathrm{H}}/n_{\mathrm{He}} \approx1.08$ (where $n_{\mathrm{H}}$ and $n_{\mathrm{He}}$ are the hydrogen and helium number densities, respectively) accounts for electrons from singly ionized helium.⁶ Motivated by results from high-resolution radiative hydrodynamics simulations [60, 61, 62], we adopt the power-law ansatz $C_A = 3 \times\left(\Gamma_{\mathrm{HI}}/3 \times10^{-13}\mathrm{s}^{-1}/x_{\mathrm{HII}}\right)^{1-\xi}$. Since $\lambda_i \propto\Gamma_{\mathrm{HI}}/C_A$, this scaling is equivalent to $\lambda_i \propto\Gamma_{\mathrm{HI}}^{\xi}$. Here we use the widely adopted value of $\xi= 2/3$ [see e.g. 16, 63–65].

In principle, the parameters describing $\lambda_n$ and $\lambda_i$ could be incorporated into the inference and marginalized over. To reduce computational cost for this illustrative analysis, we leave the parameters fixed to their simulation-motivated values. We parameterize the comoving ionizing emissivity with a double power law,

where, $N_*$, $z_*$, $\alpha_*$ and $\beta_*$ are the free parameters that we sample over. We take uniform priors of $\log(N_*) \in[-23, -20]$, $z_* \in[5, 12]$, $\alpha_* \in[-8, 0]$, $\beta_* \in[-6, 6]$, and require (without loss of generality) that $\alpha_* < \beta_*$. In this case, $\alpha_*$ sets how quickly the ionizing emissivity grows at $z > z_*$, and $\beta_*$ controls the evolution at $z < z_*$.

We use measurements of the spatially averaged hydrogen photoionization rate, $\langle\Gamma_{\mathrm{HI}}\rangle$, from Ref. [66]. These measurements are inferred from the Ly$\alpha$ forest transmission in a sample of 67 quasar spectra spanning $4.9 \lesssim z \lesssim6.0$. Ref. [66] forward-modeled the observations by post-processing hydrodynamic simulations with spatially fluctuating ionizing radiation fields and extracting mock Ly$\alpha$ forest spectra. They then compared the distribution of effective optical depths, $\tau_{\mathrm{eff}}$, between the observed and mock spectra to jointly constrain $\langle\Gamma_{\mathrm{HI}}\rangle$ and the ionizing-photon mean free path. The reported uncertainties account for the thermal state of the IGM, modeling assumptions, cosmic variance, and observational uncertainties such as noise, sky subtraction, and continuum modeling. In this work, we use only the marginalized measurements of $\langle\Gamma_{\mathrm{HI}}\rangle$. Our analysis therefore uses essentially the redshift evolution of the mean Ly$\alpha$ forest transmission, but does not attempt to model and leverage sightline-to-sightline fluctuations in the transmission.

For simplicity, we approximate the likelihood for the $\Gamma_{\mathrm{HI}}$ measurements with piecewise-asymmetric Gaussians and assume that the redshift bins are independent. For each of the 12 redshift bins, the width on either side of the mean is the relevant $1\sigma$ upper or lower measurement error. This approach neglects correlations between redshift bins arising from shared modeling assumptions and observational systematics in the inference of $\Gamma_{\mathrm{HI}}$. To assess the effects of neglecting such correlations, we performed two tests in which we retained only four and one of the $\Gamma_{\mathrm{HI}}$ measurements in our analysis (see Appendix § 9). We find that the posteriors in the four-redshift test are very similar to those in our fiducial analysis (with 12 redshift bins). Notably, the fully marginalized posterior for $\Delta z_{90}$ remains similar even in the one-redshift test run. Based on these findings, we argue that our results and conclusions would likely be similar in a full analysis that accounts for the covariance between $\Gamma_{\mathrm{HI}}$ redshift bins (all other details remaining the same). Because of the computational cost of solving differential equations in our semi-analytic model millions of times, and the approximate nature of our $\Gamma_{\mathrm{HI}}$ likelihood, in this section we keep the cosmological parameters fixed at: $A_s = 2.1 \times10^{-9}$, $n_s = 0.965$, $h = 0.674$, $\Omega_b h^2 = 0.022$, $\Omega_c h^2 = 0.12$, consistent with Planck measurements [40, 51].⁷

Results

Figure 5 shows the posterior constraints on the four parameters of the double-power-law emissivity model in Eq. (9). The values of $\alpha_*$ and $\beta_*$ determine the asymptotic slopes at $z > z_\ast$ and $z < z_\ast$, respectively, while $z_\ast$ sets the redshift at which the emissivity changes slope and $N_\ast$ determines its normalization.⁸ For reference, more negative values of $\alpha_\ast$ correspond to a more rapidly increasing emissivity (with cosmic time) before $z = z_\ast$, whereas positive values of $\beta_\ast$ imply that the emissivity declines toward lower redshift below the break.

The CMB Only posterior exhibits substantial degeneracies among the emissivity parameters. In particular, $N_\ast$ and $z_\ast$ are strongly anticorrelated: models in which the break occurs at lower redshift require a larger normalization to produce sufficient ionization at earlier times. These larger normalizations are in turn associated with more negative values of $\alpha_\ast$, which suppress the emissivity at $z > z_\ast$ and prevent excessive early ionization. Adding the dark-pixel and damping-wing-like constraints removes part of the allowed parameter space, but leaves this degeneracy broadly intact.

The measurements of $\Gamma_{\mathrm{HI}}$ substantially sharpen the constraints. Most notably, they tightly restrict the normalization of the emissivity and remove the high-$N_\ast$, low-$z_\ast$ tail that is allowed by the CMB plus neutral-fraction constraints. Models with high-$N_\ast$ and low-$z_\ast$ can stay consistent with $\tau_{\mathrm{es}}$ by combining a late break with a large normalization and a steep high-redshift slope, but they yield an ionizing background intensity near $z \sim5 - 6$ that is too large. The $\Gamma_{\mathrm{HI}}$+CMB posterior instead favors a break near $z_\ast\sim8$, together with an emissivity that declines toward lower redshift. However, significant correlations between $z_\ast$, $\alpha_\ast$, and $\beta_\ast$ remain.

The relative roles of $\alpha_\ast$ and $\beta_\ast$ depend on the location of the double power-law break. When $z_\ast\sim6$, the $\Gamma_{\mathrm{HI}}$ measurements cover a relatively short redshift interval below $z_\ast$ and therefore provide relatively little leverage on $\beta_\ast$. In this region of the parameter space, $\alpha_\ast$ mostly controls the preceding reionization history. For scenarios in which $z_\ast\sim8$–9, the forest measurements lie well below the break and therefore better constrain $\beta_\ast$. In this case, the value of $\alpha_\ast$ can vary over a broader range of negative values. These effects can be seen in the degeneracies between $\alpha_\ast$, $\beta_\ast$ and $z_\ast$ in Figure 5.

Image not extracted

The CMB plays a key role in limiting the amount of early ionization. At fixed $N_\ast$ and $z_\ast$, making $\alpha_\ast$ less negative increases the emissivity at high redshift and can produce an excessively large $\tau_{\mathrm{es}}$ and low-$\ell$ $E$-mode polarization signal. On the other hand, the forest data anchor the ionizing background near the end of reionization and particularly constrain the evolution of the emissivity near this epoch. The two data sets therefore provide complementary constraints on different portions of the emissivity history.

Figure 6 illustrates how the constraints on the double-power-law emissivity parameters propagate into the ionizing emissivity, the global reionization history, and the hydrogen photoionization rate. The three columns correspond to the CMB Only, DP+DWLIKE+CMB, and $\Gamma_{\mathrm{HI}}$+CMB analyses, respectively. From top to bottom, the rows show the comoving ionizing emissivity ($\dot{n}_{\mathrm{ion}}$), the ionized fraction ($x_{\mathrm{HII}}$), a zoomed-in view of $x_{\mathrm{HII}}$ during the final stages of reionization, and the spatially averaged photoionization rate ($\Gamma_{\mathrm{HI}}$). The dark and light shaded regions denote the marginalized 68% and 95% credible intervals. In the panels showing $x_{\mathrm{HII}}$, for reference we overplot the observational constraints from Refs. [18, 35, 36], while the lower panels show the $\Gamma_{\mathrm{HI}}$ measurements from Ref. [66].

Image not extracted

The top row highlights the limited ability of the CMB and the direct neutral-fraction constraints to determine the ionizing emissivity near the end of reionization. Both the CMB Only and DP+DWLIKE+CMB posteriors allow a wide range of emissivities at $z \lesssim6$. The $\Gamma_{\mathrm{HI}}$ measurements substantially reduce this freedom. The $\Gamma_{\mathrm{HI}}$+CMB posterior favors an emissivity that rises toward a broad maximum at $z \sim7 - 8$ and subsequently declines. This is qualitatively similar to the behavior of the emissivities in reionization simulations that are calibrated to reproduce the evolution of the mean Ly$\alpha$ forest transmission [e.g. 24, 67, 68]. There has been some debate as to whether the rapid decline required in the simulations is an artifact of modeling deficiencies [see 68]. Our results provide additional evidence that the emissivity of the reionization source population may in fact decline after $z \sim7 - 8$. If confirmed, such a decline may owe to a decline in the escape fraction [69, 70, 71] and/or the onset of a feedback process, e.g. suppression of star formation by reionization itself [72, 73, 74, 75].

Although the forest measurements directly constrain $\Gamma_{\mathrm{HI}}$ only at $z \lesssim6$, the semi-analytic model connecting $\Gamma_{\mathrm{HI}}$ to the emissivity and reionization history allows this information to propagate to higher redshifts. This can be seen in the middle two rows, which show the corresponding effect on the inferred reionization history. As we saw in the previous sections, the CMB combined with direct neutral fraction constraints still allows a broad range of histories with different durations and completion redshifts. In contrast, the $\Gamma_{\mathrm{HI}}$+CMB analysis localizes the final stages of reionization to $z \sim6-8$. Notably, these histories are also broadly consistent with the neutral-fraction constraints from Refs. [18, 35, 36], even though those measurements are not included in the $\Gamma_{\mathrm{HI}}$+CMB likelihood.

The origin of this sensitivity can be understood from the dependence of the ionizing background on both the emissivity and the mean free path of ionizing photons. During reionization, $\Gamma_{\mathrm{HI}} \propto\dot{n}_{\mathrm{ion}}\lambda_{\mathrm{mfp}}$, up to factors describing the ionizing radiation spectrum and column density distribution of absorbers. Near the end of reionization, the mean free path grows rapidly as the remaining neutral regions disappear and the distribution of ionizing absorbers evolves and the last islands of neutral hydrogen disappear. Within our model, a history that completes reionization too early therefore tends to produce a mean free path, and hence a photoionization rate, that is too large at fixed redshift. Conversely, completing reionization too late suppresses $\Gamma_{\mathrm{HI}}$ below the observed values because the ionizing background is not able to grow enough when neutral islands are still present. Matching the measured amplitude and redshift evolution of $\Gamma_{\mathrm{HI}}$ therefore restricts both the timing of the end of reionization and the evolution of the emissivity. Indeed, Ref. [33] found that shifting the end of reionization by only $\Delta z \sim0.5$ can change the predicted $\Gamma_{\mathrm{HI}}$ by more than the quoted $1\sigma$ observational uncertainties.⁹

The bottom row Figure 6 highlights this connection. In the CMB Only and DP+DW Like+CMB runs, the allowed $\Gamma_{\mathrm{HI}}$ spans a few orders of magnitude. Incorporating the $\Gamma_{\mathrm{HI}}$ constraints, as expected, produces a much narrower band that follows both the normalization and redshift evolution of the observed photoionization rate. The key result, however, is the corresponding reduction in the allowed emissivity and reionization histories. The CMB limits the amount of ionization at high redshift, whereas the Ly$\alpha$ forest anchors the ionizing background near the end of reionization; together, they constrain substantially more of the reionization history than either data set alone.

Figure 7 shows the marginalized constraints on the midpoint and duration, which are also summarized in Table 1. The CMB Only and DP+DWlike+CMB posteriors both favor $z_{50} \simeq7.5$ and relatively short reionization histories, although their $\Delta z_{90}$ posteriors retain extended tails toward larger $\Delta z_{90}$. The secondary peak in the DP+DWlike+CMB posterior arises from the same competition between the CMB and the lower limit on $x_{\mathrm{HI}}$ from Ref. [18] discussed in § 3.1. Although the posteriors for $\Delta z_{90}$ and $A_z$ in the CMB Only and DP+DWlike+CMB analyses appear to close on both sides, we will find below that the shapes of these are highly sensitive to the implicit priors, i.e. that the data are not very constraining on these quantities. We therefore quote only a 95% upper limit for $\Delta z_{90}$ in Table 1 for these analyses.

Image not extracted

Including the $\Gamma_{\mathrm{HI}}$ measurements has only a modest effect on the $z_{50}$ posterior, but shifts and substantially narrows the posterior for $\Delta z_{90}$ toward longer durations. The posterior for $A_z$ is shifted toward more symmetric reionization histories. The origin of the longer duration can be understood from the inferred evolution of the ionized fraction. Matching the observed evolution of $\Gamma_{\mathrm{HI}}$ favors a later completion of reionization, with the IGM retaining a larger neutral fraction to lower redshift. This delayed completion is compensated by an earlier start to reionization in order to preserve roughly the same $\tau_{\mathrm{es}}$, as required by the CMB constraints. Most of this additional high-redshift ionization occurs at relatively small $x_{\mathrm{HII}}$, allowing $z_{50}$ to remain broadly similar to the CMB Only and DP+DWlike+CMB results. The resulting histories in the $\Gamma_{\mathrm{HI}}+$CMB analysis are therefore more extended and symmetric. The decreased asymmetry is explained by $z_{05}$ and $z_{95}$ increasing/decreasing by roughly similar amounts, which brings $A_z$ down. In the $\Gamma_{\mathrm{HI}}+$CMB analysis (with fixed IGM opacity and cosmological parameters), we find

The central conclusion here is that combining tight constraints on reionization’s end (provided here by the $\Gamma_{\mathrm{HI}}$ measurements) with the CMB turns the latter from primarily a constraint on the midpoint of reionization into a more sensitive probe of the reionization history’s shape. This confirms the complementarity first identified by Ref. [30] and extends that result beyond the optical-depth-level by using the full CMB temperature and polarization power-spectrum likelihoods and performing a Bayesian analysis using a combination of constraints.

How sensitive are these constraints to the choice of priors? We again address this question by flattening the prior distribution in $z_{50}$, $\Delta z_{90}$, and $A_z$, using the importance sampling methods described in the previous section. The resulting posteriors for the CMB Only and $\Gamma_{\mathrm{HI}}+$CMB analyses are shown in Figure 8. Comparing those results to Figure 7, we find that the CMB Only posteriors for all three quantities widen significantly upon adopting a uniform prior. This effect owes to the physical constraints imposed by the equations of the semi-analytic model and the limited constraining power of the CMB likelihoods alone. In our fiducial analysis, we assume broad uniform priors on the parameters of the double power-law emissivity, equation (9). Propagating these priors through the physics equations of the semi-analytic model yields an implicit prior distribution for $z_{50}$, $\Delta z_{90}$, and $A_z$. Flattening this distribution in post-processing effectively removes the physical constraints on these quantities, which accounts for the widening of the CMB Only posteriors in Figure 8. By contrast, the $\Gamma_{\mathrm{HI}}+$CMB posteriors are insensitive to the choice of priors, as evidenced by the similarity of those posteriors in Figure 7 and Figure 8. After prior flattening, we find $z_{50}=7.4^{+0.2}_{-0.3}$, $\Delta z_{90}=5.9^{+0.8}_{-1.1}$, and $A_z=2.3^{+0.3}_{-0.7}$, which are very similar to (4.4), (4.5), and (4.6). The similarity owes to the increased constraining power afforded by the forest $\Gamma_{\mathrm{HI}}$ constraints.

Image not extracted

We conclude this section with two important caveats. First, we have used an approximate form of the $\Gamma_{\mathrm{HI}}$ likelihood for this analysis. Secondly, we have held the cosmological parameters, the clumping factor normalization, and its scaling with $\Gamma_{\mathrm{HI}}$ fixed. We have also adopted a simulation-calibrated form for the ionizing bubble size distribution in our reionization modeling (see eqs. 4.1, 4.2 and surrounding text). A more detailed analysis will require: (1) A more careful treatment of possible non-Gaussianities and covariances in the $\Gamma_{\mathrm{HI}}$ measurements; (2) Marginalizing over cosmology as well as the parameters describing ionizing photon absorption. The preceding analysis nevertheless illustrates the constraining power of the Ly$\alpha$ forest transmission when combined with the CMB.

Conclusions

We have investigated how Lyman-series forest measurements complement the CMB temperature and polarization anisotropies in constraining reionization. A CMB-only analysis permits a broad range of degenerate reionization histories, despite yielding tight constraints on the integrated Thomson optical depth. For instance, a relatively late midpoint can be compensated by an extended high-redshift tail, producing nearly the same optical depth and large-scale E-mode polarization signal. Incorporating Lyman-series forest constraints on the end of reionization can alleviate the degeneracy and narrow the range of allowed histories.

To explore the information content of different Lyman-series forest statistics, we have considered 3 forms of constraints: (1) Upper limits on the IGM neutral fraction at $4.9 \lesssim z \lesssim6.9$ from dark pixel counts [35, 36]; (2) A lower limit on the neutral fraction at $z=5.8$ based on damping-wing-like features adjacent to Ly$\beta$ dark gaps in quasar spectrum stacks [18]; (3) Measurements of the mean IGM photoionization rate at $4.9 \lesssim z \lesssim6$ derived from Ly$\alpha$ forest transmission measurements [66]. To explore model dependencies, we also compared 3 reionization modeling approaches. We first considered a parametric model in which the reionization history is described by a Weibull function, characterized by a midpoint, duration, and asymmetry. We also considered the more free-form FlexKnots approach. Lastly, we considered a new physical model from Ref. [33] based on the cosmological radiative transfer equation.

Our main conclusions can be summarized as follows:

• The dark-pixel and damping-wing-like constraints provide complementary information on the reionization history. Incorporating the dark pixel upper limits on $\langle x_{\mathrm{HI}}\rangle_{\mathrm{V}}$ into a joint CMB-forest analysis imposes a lower-bound on the midpoint of reionization and favors shorter reionization durations relative to a CMB-only analysis. On the other hand, incorporating the damping-wing-like lower limit on $\langle x_{\mathrm{HI}}\rangle_{\mathrm{V}}$ favors later midpoints and more extended reionization histories.

• Combining the dark pixel and damping-wing-like measurements with the CMB yields tighter constraints on reionization’s midpoint and duration. In our DP+DW Like+CMB analysis using the parametric reionization model, we find $z_{50}=7.5^{+0.5}_{-0.9}$ and $\Delta z_{90}<11.5$ (95% upper limit). For the same analysis using FlexKnots, we find $z_{50}=7.7^{+0.7}_{-0.8}$ and $\Delta z_{90}<11.4$. The asymmetry of the reionization history, however, is poorly constrained by these data.

• Neither the choice between our parametric and non-parametric descriptions, nor the choice among different forest data combinations, produces a significant shift in $\tau_{\mathrm{es}}$ or the other cosmological parameters considered here. Instead, the forest measurements primarily redistribute the optical depth in redshift. Constraints on the high-redshift ($z>15$) contribution to the optical depth are more sensitive to these choices.

• Stronger constraints can be obtained by leveraging measurements of the forest transmission at $z\lesssim6$. In this work, we have used measurements of the IGM photoionization rate, $\Gamma_{\mathrm{HI}}$, as proxy for the mean transmission. Using the semi-analytic reionization model of Ref. [33], we demonstrated that the $\Gamma_{\mathrm{HI}}$ measurements strongly restrict the allowed emissivity and reionization histories compared to the corresponding CMB-only analysis. In the $\Gamma_{\mathrm{HI}}+$CMB analysis (with fixed IGM opacity model and fixed cosmological parameters), we find $z_{50}=7.4^{+0.2}_{-0.3}$, $\Delta z_{90}=6.0^{+0.8}_{-1.0}$, and $A_z=2.4^{+0.3}_{-0.7}$. In qualitative agreement with Ref. [30], our analysis confirms that the CMB becomes a more sensitive probe of the reionization history once the last phases of reionization are tightly constrained by forest measurements.

• Our joint analysis with $\Gamma_{\mathrm{HI}}$ and CMB measurements favors models in which the ionizing emissivity of reionization sources begins declining after $z\approx7-8$. This evolution may owe to a decline in the escape fraction of ionizing photons and/or the onset of a feedback process regulating star formation.

Taken together, our results show that combining the Lyman-series forest with CMB power spectra provides access to the chronology of reionization, rather than merely its integrated optical depth. Improvements in high-redshift forest measurements, their associated likelihood products, and large-scale CMB polarization measurements should therefore translate directly into stronger joint constraints on the full reionization history and on the evolution of the ionizing sources that drove it. A further challenge to this goal will be to develop computationally efficient frameworks that allow the Ly$\alpha$ forest transmission itself to be forward-modeled within a Bayesian MCMC inference pipeline. At present, predicting the forest transmission with sufficient fidelity generally requires expensive radiative hydrodynamic simulations, making direct incorporation into inference pipelines impractical. Our results provide additional motivation for developing fast and accurate forward-modeling techniques for the forest.

Acknowledgments

The authors thank Shangqi Liu, Simeon Bird, and George Becker for helpful conversations. AD used ChatGPT GPT-5.6 Sol to edit/sharpen some passages of text in this manuscript. AD acknowledges support from grant NSF AST-2045600. CC acknowledges support from the Beus Center for Cosmic Foundations at Arizona State University. Computations were run using computational resources distributed under NSF ACCESS allocations TG-PHYS210041 and TG-PHY240332.

Cosmological parameters and testing CMB likelihood choices

In Figure 9, we show a corner plot for the 6 $\Lambda$CDM parameters obtained using the FlexKnot analyses of § 3.2. The plot compares all of the different dataset combinations as indicated in the legend. We find that the cosmological parameters are stable to different combinations of the Lyman-series forest likelihoods. Although not shown in the figure, we have also compared these results to the cosmological parameters obtained using the parametric model in § 3.1, again finding no significant difference between them. In Table 2, we provide a numerical comparison between our obtained values and the Planck PR3 and PR4 analyses (see discussion below).

Image not extracted

Parameter	Planck PR3 [51]	Ref. [40] (PR4-based)	DP+DWlike+CMB (Plik)	DP+DWlike+CMB (HiLLiPoP)
$\log(10^{10}A_{\mathrm{s}})$	$3.044 \pm0.014$	$3.040 \pm0.014$	$3.057 \pm0.013$	$3.041 \pm0.014$
$n_{\mathrm{s}}$	$0.9649 \pm0.0042$	$0.9681 \pm0.0039$	$0.9646 \pm0.0043$	$0.9676 \pm0.0040$
$\tau_{\mathrm{es}}$	$0.0544 \pm0.0073$	$0.0580 \pm0.0062$	$0.0610^{+0.0058}_{-0.0065}$	$0.0594 \pm0.0061$
$\Omega_{\mathrm{b}}h^2$	$0.02237 \pm0.0015$	$0.02226 \pm0.00013$	$0.02237 \pm0.00015$	$0.02225 \pm0.000126$
$\Omega_{\mathrm{c}}h^2$	$0.1200 \pm0.0012$	$0.1188 \pm0.0012$	$0.1202 \pm0.0014$	$0.1188 \pm0.0012$
$100\theta_*$	$1.04110 \pm0.00031$	$1.04108 \pm0.00026$	$1.04109 \pm0.00031$	$1.0411 \pm0.00025$
$H_0$	$67.36 \pm0.54$	$67.64 \pm0.52$	$67.26 \pm0.60$	$67.64 \pm0.53$
$\sigma_8$	$0.8111 \pm0.0060$	$0.8070 \pm0.0065$	$0.8169 \pm0.0067$	$0.8073 \pm0.0066$

In this work we use the Pliklite and Commander likelihoods for high-$\ell$ TTTEEE and low-$\ell$ TT, respectively. Both are based on Planck PR3. For low-$\ell$ EE, we use the more recent Lollipop likelihood (PR4). Given changes to the data reduction pipeline between PR3 and PR4, in this section we also verify that the choice to mix PR3 and PR4 likelihoods does not significantly affect our main results. We perform a separate DP+DWlike+CMB analysis under our FlexKnot model where we replace Pliklite with HiLLiPoP as the high-$\ell$ TTTEEE likelihood. In this case, the likelihoods are fully up-to-date with PR4 since there were no updates to the Commander low-$\ell$ TT likelihood in the PR4 analysis. We find little difference with our fiducial analysis in the reionization history posteriors. In Figure, 10, we compare our fiducial constraints (red) on $z_{50}$, $\Delta z_{90}$, and $A_z$, against the case with HiLLiPoP (purple).

Table 2 compares the cosmological parameters inferred from our fiducial analysis (4th column) and HiLLiPoP test (5th column) against Planck PR3 [51] (2nd column) and the PR4-based results of [40] (3rd column). Relative to the PR3 baseline, our fiducial Plik-based results exhibit a modest shift toward larger values of $\tau_{\mathrm{es}}$, $A_s$, and $\sigma_8$, with each difference amounting to approximately 0.6 – 0.7$\sigma$ when the quoted uncertainties are combined in quadrature. These shifts are not independent, but instead follow the $\tau_{\mathrm{es}}$–$A_s$–$\sigma_8$ degeneracy. The remaining cosmological parameters are effectively unchanged. In contrast, the HiLLiPoP-based results closely reproduce the PR4-based constraints from [40] across all parameters. Based on the results of this section, we conclude that our choice of CMB likelihoods does not affect the conclusions of our work.

Ref. [10] performed a closely related CMB-only analysis using a 3-Flex model, similar to our FlexKnot parameterization, obtaining $\tau_{\mathrm{es}} = 0.0593^{+0.0057}_{-0.0062}$. Reweighting their posterior to impose a flat prior on $\tau_{\mathrm{es}}$ reduced this to $0.0579 \pm0.0061$, demonstrating a modest effect due to choice of prior. The raw posterior, however, agrees remarkably well with our HiLLiPoP-based result, $\tau_{\mathrm{es}} = 0.05936 \pm0.0062$, as expected given the similar combination of CMB likelihoods. The somewhat larger values of $\tau_{\mathrm{es}}$, $A_{\mathrm{s}}$, and $\sigma_{8}$ in our fiducial Plik-based analysis therefore appear to arise primarily from the choice of CMB likelihoods. Consistent with this interpretation, replacing LoLLiPoP with SimAll in our CMB-only analysis gives $\tau_{\mathrm{es}} \simeq0.0572$ and $\ln(10^{10} A_{\mathrm{s}}) \simeq3.050$, closer to the Planck 2018 constraints. Overall, the differences are small compared with the uncertainties. Any such differences between our constraints on $\tau_{\mathrm{es}}$ and other works likely owe to both the implicit $\tau_{\mathrm{es}}$ priors and different choices in the CMB likelihoods.

Image not extracted

Volume-weighted versus mass-weighted reionization histories

The dark pixel and damping-wing-like measurements described in §2.2.2 place limits on the volume-weighted mean ionized fraction of the IGM, $\langle x_{\mathrm{HII}}\rangle_{\mathrm{V}}$. The CMB temperature and polarization anisotropies, however, are sensitive to the mass-weighted mean ionized fraction, $\langle x_{\mathrm{HII}}\rangle_{\mathrm{M}}$. In principle, this difference should be accounted for when combining constraints in a joint analysis. Broadly, simulations find that $\langle x_{\mathrm{HII}}\rangle_{\mathrm{M}}$ leads $\langle x_{\mathrm{HII}}\rangle_{\mathrm{V}}$ for most of reionization because over-dense regions are the first to reionize. The difference can be as large as $\sim10$–$15\%$ near the midpoint of reionization, reducing to less than $5\%$ for $\langle x_{\mathrm{HII}}\rangle_{\mathrm{M}} > 90\%$. Near the end of reionization, these differences are comparable to the $1\sigma$ uncertainties reported in the dark pixel and damping-wing-like measurements. Applying both $\langle x_{\mathrm{HII}}\rangle_{\mathrm{M}}$ and $\langle x_{\mathrm{HII}}\rangle_{\mathrm{V}}$ in a joint analysis introduces an additional model dependence in order to convert one quantity into the other (see e.g. [47]). For simplicity, we have assumed $\langle x_{\mathrm{HII}}\rangle_{\mathrm{M}} \approx\langle x_{\mathrm{HII}}\rangle_{\mathrm{V}}$ throughout this paper. Here we test the effect of this assumption.

We construct test set of likelihoods where we convert $\langle x_{\mathrm{HII}}\rangle_{\mathrm{M}} \rightarrow\langle x_{\mathrm{HII}}\rangle_{\mathrm{V}}$ using a mapping obtained from a radiative transfer simulation of reionization. We construct an interpolation function using the Late Start/Late End simulation of Ref. [57], chosen among the three simulations in that work because it maximizes differences between $\langle x_{\mathrm{HII}}\rangle_{\mathrm{M}}$ and $\langle x_{\mathrm{HII}}\rangle_{\mathrm{V}}$. We use the interpolation function to convert the $\langle x_{\mathrm{HII}}\rangle_{\mathrm{M}}$ from the CMB analysis into $\langle x_{\mathrm{HII}}\rangle_{\mathrm{V}}$ before computing the likelihoods for the dark pixel and damping-wing-like constraints. As an example, our conversion gives $\langle x_{\mathrm{HII}}\rangle_{\mathrm{V}} \approx0.94$ for $\langle x_{\mathrm{HII}}\rangle_{\mathrm{M}} = 0.96$. We apply this model to our DP+DWLIKE+CMB analysis using the parametric reionization model of § 3.1. In Figure 11, we compare the posterior constraints from our fiducial analysis (red) against the “corrected” analysis (blue). Table 3 summarizes the marginalized constraints on $z_{50}$, $\Delta z_{90}$, and $A_z$. Note that these quantities are still defined with respect to the mass-weighted reionization history. The constraints change only modestly, with the corrected results preferring slightly higher $z_{50}$, shorter $\Delta z_{90}$, and slightly lower $A_z$. These differences, however, are well within the derived uncertainties. From these results we conclude that our main results and conclusions are not sensitive to the approximation $\langle x_{\mathrm{HII}}\rangle_{\mathrm{V}} \approx\langle x_{\mathrm{HII}}\rangle_{\mathrm{M}}$ adopted throughout this paper.

Image not extracted

Analysis	$z_{50}$	$\Delta z_{90}$, (95% UL)
Fiducial	$7.5^{+0.5}_{-0.9}$	$< 11.5$
Corrected	$7.6^{+0.5}_{-0.8}$	$< 10.9$

Covariance in $\Gamma_{\mathrm{HI}}$ measurements

In § 4, our approximate $\Gamma_{\mathrm{HI}}$ likelihood treats each redshift bin as independent. In this section, we present a simple test that help gauge the impact of this approximation on our results. Our fiducial $\Gamma_{\mathrm{HI}}+$CMB analysis in § 4 uses measurements at 12 redshifts. For our two test runs, we perform the same analysis but using only four (spaced by $\Delta z \approx0.4$) and one of those measurements. Specifically, the four-redshift run uses the bins with $\bar{z} = \{4.9, 5.2, 5.6, 6.0\}$, while the one-redshift run uses the bin with $\bar{z} = 5.5$. These runs crudely mimic the effects of correlations among the 12 redshift bins, with the latter bracketing the loss of constraining power that could result from strong correlations between redshift bins.

In Figure 12 we compare constraints on emissivity parameters from the four-redshift and one-redshift analyses against our fiducial analysis. The four-redshift results are quite similar to the fiducial results. Even the extreme one-redshift analysis recovers the core qualitative features of the fiducial analysis, but with larger uncertainties. A noteworthy difference in the one-redshift analysis is the diminished preference for $\beta_* > 0$ seen in the bottom-right panel of Figure 12, i.e. a diminished preference for a declining $\dot{n}_{\mathrm{em}}$ near the end of reionization. Naturally, constraining $\Gamma_{\mathrm{HI}}$ at only one redshift severely limits information about the $z < z_*$ slope of $\dot{n}_{\mathrm{em}}$. In Figure 13 we show the fully marginalized posteriors for $z_{50}$, $\Delta z_{90}$, $A_z$. Again, the four-redshift analysis gives results very similar to the fiducial case. Notably, the posterior for $\Delta z_{90}$ remains broadly intact even in the one-redshift analysis.

Image not extracted

Image not extracted

These simple tests suggest that our main conclusions in § 4 would be broadly the same in a full analysis that accounts for the covariance between $\Gamma_{\mathrm{HI}}$ redshift bins.