Rapid late-time reionization: constraints and cosmological implications

TL;DR

This paper constrains the reionization optical depth using multiple independent methods, finding a rapid, late reionization history largely consistent across datasets but with some tension in combined constraints.

Contribution

It provides a comprehensive, model-insensitive analysis of reionization history and optical depth using Ly$eta$ constraints, CMB polarization, and small-scale anisotropies, highlighting the consistency and tensions among different data sets.

Findings

Reionization occurred rapidly around z=7.

Optical depth $ au$ is approximately 0.0492, largely model-independent.

Some tension exists between Ly$eta$ bounds and combined CMB and BAO data.

Abstract

We present constraints on the reionization optical depth, $\tau$, obtained using several independent methods. First, we perform a non-parametric reconstruction of the reionization history, using Lyman-$\alpha$ constraints on the evolution of the volume-averaged neutral hydrogen fraction, $x_\mathrm{HI}(z)$, including recent results from the James Webb Space Telescope. When combined with baryon acoustic oscillation (BAO) measurements from DESI and Big Bang nucleosynthesis constraints, these data imply a rapid reionization history ($z_\mathrm{mid}=7.00^{+0.12}_{-0.18}$ and $\Delta z_{50}=1.12^{+0.12}_{-0.29}$) and a value of $\tau=0.0492^{+0.0014}_{-0.0030}$, which is largely insensitive to the assumed cosmological model and independent of cosmic microwave background (CMB) data. The optical depth can also be measured from large-scale $(\ell<30)$ CMB polarization data, yielding constraints that are similarly model-insensitive and consistent with the Ly$\alpha$ bound. Third, $\tau$ may be constrained from the attenuation of small-scale $(\ell>30)$ CMB anisotropies, but the results are sensitive to the choice of cosmological model. Assuming $\Lambda$CDM and combining small-scale CMB data with CMB lensing and type 1a supernovae (SNe) yields tight constraints that are compatible with the Ly$\alpha$ bound. Adding galaxy clustering and lensing measurements brings the constraints further into agreement with the Ly$\alpha$ bound. These independent results reinforce a consensus picture in which reionization is rapid and late. However, the combination of small-scale CMB, CMB lensing, and BAO data yields $\tau=0.094\pm0.011$, which is in $4\sigma$ tension with our Ly$\alpha$ bound. Non-standard reionization scenarios can reconcile some but not all constraints. Concordance is restored in alternative cosmological models, such as models with dynamical dark energy favoured by BAO, CMB, and SNe data.

Tables

Table 1. Table 1: An overview of the astrophysical constraints on the neutral hydrogen fraction used in this paper.

Reference	$z$	$x_{\mathrm{HI}}\left(z\right)$
Jin et al. (2023)	$5.5\pm 0.2$
Jin et al. (2023)	$5.7\pm 0.2$
Umeda et al. (2025)	${5.8}_{-0.4}^{+0.5}$	${0.25}_{-0.20}^{+0.10}$
Jin et al. (2023)	$5.9\pm 0.2$
Greig et al. (2024)	$5.90\pm 0.05$
Greig et al. (2024)	$5.95\pm 0.05$
Greig et al. (2024)	$6.05\pm 0.05$
Greig et al. (2024)	$6.15\pm 0.05$	${0.20}_{-0.12}^{+0.14}$
Greig et al. (2024)	$6.35\pm 0.05$	${0.29}_{-0.13}^{+0.14}$
Greig et al. (2024)	$6.55\pm 0.05$
Mason et al. (2018)	$6.9\pm 0.5$	${0.59}_{-0.15}^{+0.11}$
Umeda et al. (2025)	${7.0}_{-0.4}^{+0.5}$	${0.65}_{-0.35}^{+0.27}$
Greig et al. (2022)	$7.29\pm 0.27$	${0.49}_{-0.14}^{+0.13}$
Hoag et al. (2019)	$7.6\pm 0.6$	${0.88}_{-0.10}^{+0.05}$
Mason et al. (2019a)	$7.9\pm 0.6$	$>0.76$
Umeda et al. (2025)	${8.6}_{-0.9}^{+1.0}$	${1.00}_{-0.20}^{+0.00}$
Umeda et al. (2025)	${10.4}_{-0.8}^{+2.0}$	${1.00}_{-0.40}^{+0.00}$
Bruton et al. (2023)	$10.6\pm 0.0013$
Curtis-Lake et al. (2023)	${11.49}_{-0.05}^{+0.04}$	${0.84}_{-0.27}^{+0.14}$

Table 2. Table 2: Constraints on the reionization optical depth, τ \tau , from small-scale CMB data combined with various low-redshift probes as indicated, from large-scale CMB polarization data (lowE) for three different likelihoods, and from x HI ​ ( z ) x_{\mathrm{HI}}(z) data. In all cases, we assume the Λ \Lambda CDM model.

Dataset	$\tau$
SNe + TTTEEE + $6\times 2$pt	$0.047\pm 0.015$
SNe + TTTEEE + Lensing	$0.062\pm 0.014$
BAO + TTTEEE + $6\times 2$pt	$0.073\pm 0.012$
BAO + TTTEEE + Lensing	$0.094\pm 0.011$
BAO + ACT-SPT + Lensing	$0.094\pm 0.011$
BAO + P-ACT-SPT + Lensing	$0.0914\pm 0.0094$
BAO + BBN + lowE (SimAll)	${0.058}_{-0.018}^{+0.010}$
BAO + BBN + lowE (SRoll2)	${0.0535}_{-0.011}^{+0.0066}$
BAO + BBN + lowE (LoLLiPoP)	${0.0494}_{-0.0082}^{+0.0059}$
BAO + BBN + $x_{\mathrm{HI}}\left(z\right)$	${0.0490}_{-0.0028}^{+0.0015}$

Table 3. Table 3: Constraints for different cosmological models on the reionization optical depth, τ \tau , and the parameter combination T = τ ​ Ω m ​ h 2 / 0.14 T=\tau\sqrt{\Omega_{\mathrm{m}}h^{2}/0.14} that is more directly measured by x HI ​ ( z ) x_{\mathrm{HI}}(z) constraints. The results are given as 68% credible intervals for four different data combinations labelled 1A through 3. The final column lists the tension in T T between combinations 1A ( x HI ​ ( z ) x_{\mathrm{HI}}(z) data with a BBN prior) and 2 (constraints from BAO, Planck TTTEEE, and CMB lensing).

	1A.	1B.	2.		3.		Tension in $T$
	BBN + $x_{\mathrm{HI}}\left(z\right)$	BAO + BBN + $x_{\mathrm{HI}}\left(z\right)$	BAO + TTTEEE + Lensing		BAO + BBN + lowE (SRoll2)		between
Cosmological model	$T$	$\tau$	$T$	$\tau$	$T$	$\tau$	$1A.\leftrightarrow 2$
$\Lambda$CDM	${0.0468}_{-0.0036}^{+0.0029}$	${0.0492}_{-0.0030}^{+0.0014}$	$0.094\pm 0.011$	$0.094\pm 0.011$	${0.0535}_{-0.011}^{+0.0066}$	${0.0535}_{-0.011}^{+0.0066}$	$3.7\sigma$
$\Lambda$CDM $+{\alpha }_s$	${0.0468}_{-0.0036}^{+0.0029}$	${0.0492}_{-0.0030}^{+0.0014}$	$0.093\pm 0.011$	$0.093\pm 0.011$	${0.0520}_{-0.013}^{+0.0078}$	${0.0520}_{-0.013}^{+0.0078}$	$4.3\sigma$
$\Lambda$CDM $+A_{\mathrm{lens}}$	${0.0468}_{-0.0036}^{+0.0029}$	${0.0492}_{-0.0030}^{+0.0014}$	${0.094}_{-0.055}^{+0.035}$	${0.094}_{-0.055}^{+0.036}$	${0.0546}_{-0.011}^{+0.0069}$	${0.0546}_{-0.011}^{+0.0069}$	$0.9\sigma$
$\Lambda$CDM $+{\Omega }_k$	${0.0481}_{-0.0060}^{+0.0073}$	${0.0497}_{-0.0032}^{+0.0019}$	$0.086\pm 0.012$	$0.085\pm 0.012$	${0.0537}_{-0.012}^{+0.0079}$	${0.0535}_{-0.011}^{+0.0074}$	$2.7\sigma$
$\Lambda$CDM $+N_{\mathrm{eff}}$	${0.0474}_{-0.0059}^{+0.0072}$	${0.0493}_{-0.0031}^{+0.0019}$	$0.091\pm 0.011$	$0.091\pm 0.012$	${0.0537}_{-0.011}^{+0.0070}$	${0.0538}_{-0.011}^{+0.0071}$	$3.3\sigma$
$\Lambda$CDM $+\sum m_{\nu ,\mathrm{eff}}$	${0.0471}_{-0.0036}^{+0.0028}$	${0.0493}_{-0.0047}^{+0.0032}$	${0.063}_{-0.021}^{+0.024}$	${0.063}_{-0.021}^{+0.024}$	${0.0542}_{-0.011}^{+0.0075}$	${0.0542}_{-0.011}^{+0.0071}$	$0.8\sigma$
$w$CDM	${0.0479}_{-0.0026}^{+0.0044}$	${0.0505}_{-0.0033}^{+0.0020}$	${0.091}_{-0.012}^{+0.014}$	$0.091\pm 0.013$	${0.0528}_{-0.011}^{+0.0070}$	${0.0545}_{-0.012}^{+0.0070}$	$3.0\sigma$
$w_0{w}_a$CDM	${0.0480}_{-0.0029}^{+0.0041}$	${0.0484}_{-0.0035}^{+0.0016}$	${0.069}_{-0.016}^{+0.014}$	${0.069}_{-0.016}^{+0.014}$	${0.0540}_{-0.011}^{+0.0072}$	${0.0524}_{-0.011}^{+0.0069}$	$1.2\sigma$

Equations

$$\tau =n_{\mathrm{H}}c{\sigma }_{\mathrm{T}}{\int }_0^{z_{\max }}\mathrm{d}z{x}_{\mathrm{e}}(z)\frac{(1+z{)}^2}{H(z)},$$

$$x_{\mathrm{HI}}(z)=\frac{1}{2}\left[1-\tanh \left(\frac{(1+z_{\mathrm{mid}}{)}^{1.5}-(1+z{)}^{1.5}}{0.75(1+z_{\mathrm{mid}}{)}^{0.5}}\right)\right],$$

$$k(a_i,a_j;\sigma ,\ell )={\sigma }^2\exp \left(-\frac{(a_i-a_j{)}^2}{2{\ell }^2}\right),$$

$$T=0.0468_{-0.0036}^{+0.0029}(\text{BBN + }{x}_{\mathrm{HI}}).$$

$$\tau =0.0492_{-0.0030}^{+0.0014}(\text{Baseline; BAO + BBN + }{x}_{\mathrm{HI}}).$$

$$\displaystyle\left.\begin{array}[]{ll}\tau&=0.0544^{+0.0094}_{-0.014}\\ \Omega_{\mathrm{m}}h^{2}&=0.132^{+0.030}_{-0.080}\end{array}\right\}\quad(\text{BBN + }x_{\mathrm{HI}}\text{ + lowE ({SRoll2})}).$$

$$\tau =0.0481\pm 0.0018(\text{parametric; BAO + BBN + }{x}_{\mathrm{HI}}),$$

$$L{C}_L^{KK}=2.73\times 10^{-5}{\left[\frac{\Delta z_{90}}{4.0}\right]}^{1.74}{\left[\frac{z_{\mathrm{mid}}}{8.0}\right]}^{2.52},$$

$$\tau =0.0504_{-0.0044}^{+0.00033}(z_{\max }=30\text{; BAO + BBN + }{x}_{\mathrm{HI}}\text{ + kSZ}),$$

$$\tau =0.0526_{-0.015}^{+0.0081}(\text{lowE (}\text{SRoll2}\text{)}),$$

$$\tau =0.0535_{-0.011}^{+0.0066}(\text{BAO + BBN + lowE (}\text{SRoll2}\text{)}).$$

$$A_{\mathrm{lens}}=1.087\pm 0.021(\text{CMB + BAO + }{x}_{\mathrm{HI}}),$$

$$\displaystyle\textstyle\sum m_{\nu,\mathrm{eff}}=-0.120^{+0.034}_{-0.039}\,$\mathrm{eV}$\qquad(\text{CMB + BAO + }x_{\mathrm{HI}}).$$

$$w(a)=w_0+w_a(1-a),$$

$$\displaystyle\left.\begin{array}[]{ll}w_{0}&=-0.732\pm 0.056\\ w_{a}&=-0.99^{+0.22}_{-0.20}\end{array}\right\}\quad(\text{CMB + BAO + SNe + }x_{\mathrm{HI}}),$$

$$F(x)=\frac{1}{2}+\text{sign}(2x-1)\frac{1}{2\Gamma (1/\beta )}\gamma \left(1/\beta ,{\left|2x-1\right|}^{\beta }\right),$$

$$f(x|p,b)=\frac{b^{p/2}}{2K_p(\sqrt{b})}{x}^{p-1}{e}^{-(x+b/x)/2},x>0,$$

Full text

††pubyear: 2025††pagerange: Rapid late-time reionization: constraints and cosmological implications–A

1 Introduction

The epoch of reionization was a cosmic phase transition in which radiation from the first luminous objects ionized the neutral intergalactic medium (IGM) (Loeb & Barkana, 2001; Robertson et al., 2010; Choudhury, 2022). Establishing the timeline of reionization is a central challenge in cosmology, with profound implications for our understanding of the first stars and our interpretation of the cosmic microwave background (CMB). This primordial signal was perturbed at late times due to Thomson scattering by free electrons generated during reionization. The amount of scattering, characterized by the optical depth, $\tau$, is degenerate with other cosmological parameters of interest. Hence, reliable constraints on reionization are indispensable for cosmological model inference.

In recent years, a wide array of observations have made it possible to piece together the global history of reionization (e.g. Greig & Mesinger, 2017; Mason et al., 2019b; Hazra et al., 2020; Paoletti et al., 2025; Sims et al., 2025). Information on the timing of reionization can be derived from many different observables, such as damping wings in quasar and galaxy spectra (e.g. Mesinger & Haiman, 2007; Davies et al., 2018; Mason et al., 2018; Greig et al., 2022; Umeda et al., 2024), the fraction of dark pixels in the Lyman-$\alpha$ and Lyman-$\beta$ forests (e.g. McGreer et al., 2015; Jin et al., 2023), the clustering of Lyman-$\alpha$ emitters (e.g. Ouchi et al., 2010; Sobacchi & Mesinger, 2015), the kinematic Sunyaev-Zeldovich (kSZ) effect (e.g. Zahn et al., 2012; Smith & Ferraro, 2017; Raghunathan et al., 2024), the post-reionization Lyman-$\alpha$ forest (e.g. Cen et al., 2009), and the 21-cm power spectrum (e.g. Abdurashidova et al., 2022; Ghara et al., 2025). The arrival of the James Webb Space Telescope (JWST) has enabled observations of the onset of reionization at $z>10$ (e.g. Bruton et al., 2023; Curtis-Lake et al., 2023; Hsiao et al., 2024; Muñoz et al., 2024; Cohon et al., 2025; Llerena et al., 2025; Umeda et al., 2025).

The strongest CMB constraints on $\tau$ derive from large-scale polarization data. Thomson scattering of the CMB quadrupole results in linear polarization on the horizon scale, which produces a characteristic feature in the polarization power spectrum on large scales (Hu & White, 1997; Reichardt, 2016). Detection of this reionization feature was first achieved with the Wilkinson Microwave Anisotropy Probe (Kogut et al., 2003) and measurements were further refined with Planck (Planck Collaboration et al., 2014, 2016, 2020b), culminating in a precise determination of $\tau=0.051\pm 0.006$ (Planck Collaboration et al., 2020d). Recently, the polarization signature was detected with the CLASS telescope (Li et al., 2025), yielding $\tau=0.053_{-0.019}^{+0.018}$, in agreement with the Planck value, though with larger errors. On small scales, the first-order effect of reionization on the CMB is an overall attenuation of the signal, introducing a degeneracy between $\tau$ and the amplitude of the primordial power spectrum, $A_{\mathrm{s}}$. When small-scale CMB data are complemented with low-redshift probes to break this degeneracy, the resulting constraints on $\tau$ are competitive with other approaches (Giarè et al., 2024a).

Among low-redshift probes that can be combined with CMB data, measurements of baryon acoustic oscillations (BAO) are the most powerful (Eisenstein et al., 1998). They provide a robust geometric probe of the expansion history that effectively breaks parameter degeneracies. However, the latest BAO measurements by the Dark Energy Spectroscopic Instrument (DESI Collaboration et al., 2024, 2025c, 2025b, 2025a) are in tension with CMB experiments (Louis et al., 2025; Garcia-Quintero et al., 2025; Camphuis et al., 2025). This tension is at least partially responsible for anomalies, such as the preference for excess gravitational lensing (Calabrese et al., 2008; Planck Collaboration et al., 2020b; Mokeddem et al., 2023) and the tension between cosmological constraints on the sum of neutrino masses and bounds from neutrino oscillation experiments (Craig et al., 2024; Loverde & Weiner, 2024; Elbers et al., 2025b, a; Green & Meyers, 2025). Notably, it contributes to the $3\sigma$ preference for dynamical dark energy from the combination of DESI BAO and CMB data (DESI Collaboration et al., 2025a). To make progress on fundamental questions about neutrinos and the nature of dark energy, it is crucial to understand the origin of these anomalies. Recently, Sailer et al. (2025) and Jhaveri et al. (2025) pointed out that the tension between DESI and CMB data could be explained by unknown systematics in the large-scale polarization measurements of $\tau$, potentially allowing a greater optical depth of $\tau=0.09$. This value is much larger than predicted by standard reionization models, but could reconcile CMB measurements of the expansion history with those from DESI.

At the same time, new JWST observations suggest that galaxies were more efficient at producing ionizing radiation during the epoch of reionization (Atek et al., 2024; Simmonds et al., 2024; Llerena et al., 2025; Pahl et al., 2025; Papovich et al., 2025). How much of this radiation escapes into the IGM is crucial for reionization, but cannot be measured directly at high redshifts. The application of scaling relations, calibrated to low-redshift galaxy samples, suggests that the escape fraction was greater for reionization-era galaxies and increasing towards the faint end of the luminosity function (Chisholm et al., 2022; Mascia et al., 2024). As shown by Muñoz et al. (2024), the increased ionizing efficiencies and escape fractions, combined with incomplete observations at the faint end, could lead to an overabundance of ionizing photons, an early end to reionization, and a greater optical depth, $\tau=0.09$, in tension with CMB and Lyman-$\alpha$ observations (termed a ‘photon budget crisis’). However, more recent JWST studies have inferred the escape fractions of reionization-era galaxies by fitting their spectral energy distributions and found no relation between the escape fraction and UV brightness (Papovich et al., 2025; Giovinazzo et al., 2025), calling into question the use of relations derived from low-redshift samples. There is still considerable uncertainty regarding the escape fraction, but both studies agree that the average value is sufficiently low for a late end to reionization around $z\approx 5-6$, apparently resolving the photon budget crisis (see also Pahl et al. 2025).

Motivated by anomalies in the latest cosmological data and the new JWST observations, this paper presents constraints on $\tau$ employing several independent methods to minimize the dependence on any single dataset. First, we perform a non-parametric reconstruction of the reionization history using astrophysical constraints on the evolution of the volume-averaged neutral hydrogen fraction, $x_{\mathrm{HI}}(z)$, including the latest JWST measurements. A novelty of our analysis is that we simultaneously vary the reionization history and cosmological parameters, without relying on CMB data. The optical depth, $\tau$, is given by

[TABLE]

where $n_{\mathrm{H}}$ is the present-day density of hydrogen, $\sigma_{\mathrm{T}}$ the Thomson scattering cross section, $x_{\mathrm{e}}(z)$ the free electron fraction, and $H(z)$ the Hubble rate. We can determine $\tau$ independently of the CMB by combining constraints on $x_{\mathrm{HI}}(z)$ with BAO and Big Bang nucleosynthesis (BBN) measurements. Together, these data effectively constrain the hydrogen density, $n_{\mathrm{H}}$, and the expansion history, $H(z)$. We will compare the results of this reconstruction with constraints on $\tau$ obtained from small-scale and large-scale CMB data in combination with various low-redshift probes.

The remainder of the paper is structured as follows. We describe the data used in our analysis in Section 2. The non-parametric reconstruction of the reionization history is presented in Section 3. In Section 4, we compare the results with constraints on $\tau$ from the CMB and present results for extended cosmological models. Finally, we present our conclusions in Section 5.

2 Data

We make use of both astrophysical constraints on reionization and various other cosmological datasets as outlined in the two subsections below. These data are used for Bayesian cosmological parameter inference, employing a Markov Chain Monte Carlo algorithm implemented in the cobaya code (Torrado & Lewis, 2021). Cosmological calculations are performed with a modified version of the CLASS code (Lesgourgues, 2011). The chains are finally analyzed with the getdist package (Lewis, 2019).

2.1 Reionization data

The primary dataset used in our analysis is a compilation of recent constraints on the evolution of the volume-averaged neutral hydrogen fraction, $x_{\mathrm{HI}}(z)$, from the literature. The constraints are listed in Table 1 and derive mostly from damping wing features in quasar and galaxy spectra. We also use the Lyman-$\alpha$ and Lyman-$\beta$ forest dark pixel fraction constraints of Jin et al. (2023), which are similar to but more conservative than the commonly-used constraints from McGreer et al. (2015). In compiling this dataset, we gave priority to data points that were more recent. We also required that the probability density was publicly available (with one exception discussed below) and that the data were independent.

The values in Table 1 are given only for reference. In practice, we sample directly from the marginalized posterior probability density, $p(x_{\mathrm{HI}})$. Redshift errors are modelled by including one nuisance parameter for each data point, imposing a uniform prior, $\Delta z\in[-z_{\text{min}},z_{\text{max}}]$, with the values given in Table 1. This is a conservative choice, but only has a small impact on the results. Additionally, the following comments apply to specific data points:

The probability densities were not available for Jin et al. (2023). Following Greig & Mesinger (2017) and Mason et al. (2019b), we model the upper bounds $x_{\mathrm{HI}}<\mu\pm\sigma$ as one-sided Gaussians with $p(x_{\mathrm{HI}})=1$ for $x_{\mathrm{HI}}<\mu$ and $p(x_{\mathrm{HI}})=\phi_{\mathrm{Normal}}(x_{\mathrm{HI}};\mu,\sigma)$ for $x_{\mathrm{HI}}\geq\mu$. We tested an alternative approach in which the upper bounds were modelled as latent variables, but the results were similar. 2. 2.

The measurements of Mason et al. (2018, 2019a) and Hoag et al. (2019) assume that the IGM is fully ionized at $z=6$. Following Mason et al. (2019b), we interpret their posteriors, $p(x_{\mathrm{HI}})$, as posteriors for the difference $\Delta x_{\mathrm{HI}}=x_{\mathrm{HI}}(z)-x_{\mathrm{HI}}(6)$. 3. 3.

Following Greig et al. (2022), we average the posteriors from Greig et al. (2022) with those from Yang et al. (2020) for J1007+2115, Wang et al. (2020) for J0252-0503, and Davies et al. (2018) for J1120+0641 and J1342+0928 to obtain the $z=7.29$ data point.

2.2 Other cosmological data

Additional data are required to relate the ionization history, $x_{\mathrm{HI}}(z)$, to the optical depth, $\tau$, when varying all cosmological parameters. To do this, we make use of a BBN prior on the physical baryon density, $\Omega_{\mathrm{b}}h^{2}$, from Schöneberg (2024).111When we vary the effective number of relativistic species, $N_{\mathrm{eff}}$, we also incorporate the covariance between $N_{\mathrm{eff}}$ and $\Omega_{\mathrm{b}}h^{2}$ (Schöneberg, 2024). We constrain the cosmic expansion history using the latest BAO data from Data Release 2 (DR2) of the DESI Collaboration et al. (2025b, a).

We also consider a range of CMB datasets from the Planck Collaboration et al. (2020b, a). We specifically use the CamSpec likelihood for the small-scale ($30\leq\ell\leq 2500$) auto and cross power spectra of temperature and polarization anisotropies (Efstathiou & Gratton, 2019; Rosenberg et al., 2022), denoted as TTTEEE. When considering large-scale ($\ell<30$) polarization data from Planck, denoted as lowE, we compare three different likelihoods to determine the impact of differences in implementation. We consider the SimAll likelihood from the Planck Collaboration et al. (2020a), the SRoll2 likelihood from Delouis et al. (2019), and the LoLLiPoP likelihood from Tristram et al. (2021, 2024). We also use small-scale CMB data from the Atacama Cosmology Telescope (ACT) Data Release 6 (Louis et al., 2025; Calabrese et al., 2025) and the South Pole Telescope (SPT) ‘SPT-3G D1’ (Camphuis et al., 2025). The combination of ACT and SPT data is denoted as ACT-SPT, where we neglect their cross correlation due to the small overlap in area (Camphuis et al., 2025). Since ACT and SPT do not measure large-scale CMB anisotropies, we also employ a combination of Planck, ACT, and SPT denoted as P-ACT-SPT. Specifically, we combine the SPT likelihood with the P-ACT combination that relies on Planck on intermediate scales ($30<\ell<600$) and ACT on smaller scales ($\ell\geq 600$) (Louis et al., 2025).222The only difference with the P-ACT combination in Louis et al. (2025) is that we exclude large-scale temperature data (lowT), which can play a role when lowE is also excluded (e.g. Giarè et al., 2024a). However, we verified that including lowT only has a small impact once Planck or P-ACT is combined with DESI BAO, slightly increasing the tension with our Lyman-$\alpha$ bound. We also make use of CMB lensing data from Planck (Planck Collaboration et al., 2020c; Carron et al., 2022), ACT (Madhavacheril et al., 2024; Qu et al., 2024; MacCrann et al., 2024), and SPT (Ge et al., 2025), using the extended version of their joint likelihood (Qu et al., 2025), denoted simply as Lensing.

As an alternative to BAO data, we also rely on type 1a supernovae (SNe) to constrain the late-time expansion history. In this case, we employ the Dark Energy Survey (DES) Year 5 compilation of 1635 photometrically confirmed SNe combined with 194 low-$z$ external SNe (DES Collaboration et al., 2024). Finally, we make use of a joint analysis of 2-point correlation functions of galaxy positions, galaxy lensing, CMB lensing, and their cross correlations from DES, Planck, and SPT (Omori et al., 2023; Chang et al., 2023; Abbott et al., 2023a), which we denote as $6\times 2$pt. In the $6\times 2$pt likelihood, we adopt the same scale cuts for the dynamical dark energy model $w_{0}w_{a}$CDM as for $\Lambda$CDM, following Abbott et al. (2023b). We adapt the public likelihood code to a version compatible with cobaya using the CosmoSIS2cobaya package333https://github.com/JiangJQ2000/cosmosis2cobaya, following Ye et al. (2024).

3 Non-parametric reconstruction

The reionization history, $x_{\mathrm{HI}}(z)$, has been modelled using a wide range of methods, from simple functional forms to more complex and flexible approaches (e.g. Hazra & Smoot, 2017; Millea & Bouchet, 2018; Trac, 2018; Sailer et al., 2022; Montero-Camacho et al., 2024; Ilić et al., 2025; Paoletti et al., 2025). The most common choice, used by Planck, is to adopt a smooth step function given by

[TABLE]

where $z_{\mathrm{mid}}$ is the midpoint of reionization. While this approach is simple and useful for many applications, non-parametric methods offer greater flexibility with the distinct advantage of not excluding any reionization histories a priori (Mason et al., 2019b; Krishak & Hazra, 2021; Cheng et al., 2025). We adopt such an approach in this work, but also compare with fits to Eq. (2).

3.1 Methods

Here, we perform a non-parametric reconstruction of the reionization history using Gaussian process regression. We model the neutral hydrogen fraction by introducing a Gaussian process, $\{\chi_{a}\}$, where $a=(1+z)^{-1}$ is the scale factor and $X=(\chi_{a_{1}},\dots,\chi_{a_{n}})$ is a multivariate Gaussian random variable. The $\chi_{a}$ are unbounded latent variables that are mapped to the neutral hydrogen fraction, $x_{\mathrm{HI}}(a)$, by means of a wrapping function $F\colon\mathbb{R}\to(0,1)$.

This method allows us to generate smooth functions by sampling large numbers of points from a multivariate normal distribution with mean, $\mu(a)$, and a covariance matrix, $C(a_{i},a_{j})$, determined by the kernel function. We use a squared exponential kernel, given by

[TABLE]

where $\sigma$ is the typical size of fluctuations in $x_{\mathrm{HI}}$ around the mean. The correlation length is determined by the parameter $\ell$, which we include as a free hyper parameter. Further technical details of the Gaussian process are given in Appendix A.

We impose two boundary conditions at the endpoints of reionization. We require $x_{\mathrm{HI}}(z_{\mathrm{min}})=\epsilon$ at $z_{\mathrm{min}}=5.2$, motivated by observations that imply completion of reionization around $z\approx 5.2$ (e.g. Bosman et al., 2022; Gaikwad et al., 2023).444We use $\epsilon=10^{-3}$ because $x_{\mathrm{HI}}(z)$ is modelled on the open interval $(0,1)$. The beginning of reionization is much more uncertain. Evidence from large-scale CMB polarization and constraints on the duration of reionization from the kSZ effect (Zahn et al., 2012; Planck Collaboration et al., 2016; Reichardt et al., 2021; Smith & Ferraro, 2017; Raghunathan et al., 2024) suggest that there is little contribution to $\tau$ beyond $z>15$ (Planck Collaboration et al., 2016, 2020b), but see also e.g. Heinrich & Hu (2021); Cheng et al. (2025). Although our earliest data point is at $z=11.49$, we fix $x_{\mathrm{HI}}(z_{\mathrm{max}})=1-\epsilon$ at $z_{\mathrm{max}}=15$, to allow for early and gradual reionization scenarios (Finkelstein et al., 2019; Asthana et al., 2024; Lu et al., 2024; Muñoz et al., 2024; Cain et al., 2025b). We will relax this assumption in Section 3.3, where we consider an alternative analysis with $z_{\mathrm{max}}=30$ to allow for more extended reionization scenarios.

We fix the remaining two hyper parameters at $\sigma=0.3$ and $\mu(z)=\mu=0$. This choice is motivated by the behaviour at high redshifts. At low redshifts, the evolution of $x_{\mathrm{HI}}(z)$ is strongly constrained by the data. In the absence of constraints from the data or the boundary conditions, the results will be driven by the prior. For $z>12$, we reasonably assume that $x_{\mathrm{HI}}(z)\approx 0$ with an uncertainty of $\sigma=0.3$.

To allow for consistent comparisons with other analyses, we require a uniform prior on $\tau$. We accomplish this with the method of Millea & Bouchet (2018) and Handley & Millea (2019). Using Monte Carlo sampling, we estimate the implicit prior on $\tau$ determined by the choices outlined above. We then reweight samples accordingly in order to obtain a flat prior on $\tau$. This correction has a minimal effect on the results.

We assume that the first helium reionization happens simultaneously with hydrogen reionization and the second helium reioniziation is modelled as a smooth transition of the form (2) at $z=3.5$, following Planck Collaboration et al. (2020b). The impact of helium reionization on the optical depth is negligible compared to the uncertainty in $\tau$ from hydrogen reionization and this remains the case for the results presented in this paper.

3.2 Baseline results

In this section, we only report results for the $\Lambda$CDM model. Fitting the reionization model to the $x_{\mathrm{HI}}(z)$ constraints summarized in Table 1, we obtain the reconstructed reionization histories shown in Fig. 1. Overall, we find that the Gaussian process can describe the trends in the data well. The maximum a posteriori realization has $\chi^{2}=5.6$ for $19$ data points. The reconstruction favours rapid and late reionization histories that occur mostly between $6\lesssim z\lesssim 8$. Defining the midpoint of reionization, $z_{\mathrm{mid}}$, as the redshift at which $x_{\mathrm{HI}}=0.5$ and the duration as $\Delta z_{50}=z_{75}-z_{25}$, where $z_{75}$ and $z_{25}$ are the redshifts at which $x_{\mathrm{HI}}=0.75$ and $x_{\mathrm{HI}}=0.25$, we find $z_{\mathrm{mid}}=6.99^{+0.14}_{-0.17}$ and $\Delta z_{50}=1.10^{+0.14}_{-0.28}$.555Technically, $z_{25},z_{\mathrm{mid}}$, and $z_{75}$ are not necessarily unique, but this issue only affects a minimal fraction of reconstructed histories. We ignore the issue here, but mitigate it when we consider extended histories in Section 3.3.

To begin to constrain the optical depth, we first approximately factor out the dependence on $H(z)$ by defining the parameter $T\equiv\tau\sqrt{\Omega_{\mathrm{m}}h^{2}/0.14}$.666The value $0.14$ is motivated by the result $\Omega_{\mathrm{m}}h^{2}=0.1399\pm 0.0050$ from DESI DR2 BAO and BBN. With the addition of the BBN prior, we obtain a tight constraint,

[TABLE]

Fig. 2 shows that by adding DESI BAO data, we can break the degeneracy between $\tau$ and $\Omega_{\mathrm{m}}h^{2}$ and thereby constrain the optical depth itself. The rapid reionization histories preferred by the data then correspond to

[TABLE]

As an aside, we can also use $\tau$ as a “standard depth” and constrain $\Omega_{\mathrm{m}}h^{2}$ through the combintation of $x_{\mathrm{HI}}(z)$ data with an independent probe of $\tau$. We demonstrate this by adding large-scale CMB polarization data using the SRoll2 likelihood (see the next section) in Fig. 2. This combination yields

[TABLE]

The resulting constraint on $\Omega_{\mathrm{m}}h^{2}$ is compatible with standard CMB and BAO + BBN measurements, but the error bars are not competitive. Future CMB experiments could measure $\tau$ with an uncertainty of $0.002$ (LiteBIRD Collaboration et al., 2023), which may increase the utility of this method.

Including BAO and BBN constraints in the likelihood helps to break parameter degeneracies, but does not significantly impact the $x_{\mathrm{HI}}(z)$ fit. The timing of reionization is essentially unchanged ($z_{\mathrm{mid}}=7.00^{+0.12}_{-0.18}$ and $\Delta z_{50}=1.12^{+0.12}_{-0.29}$). In addition to the GP reconstruction, Fig. 1 also shows some representative reionization histories that are more gradual and achieve larger optical depths, $\tau=0.07$ or $\tau=0.09$, which could alleviate or resolve the tension between DESI BAO and CMB data. These histories are inconsistent with the Lyman-$\alpha$ constraints, particularly around $7\lesssim z\lesssim 8$.

The resulting one-dimensional marginalized posterior distribution of $\tau$ is shown as a thick solid line in Fig. 3. The figure also shows the impact of various modelling assumptions. Neglecting redshift errors has a negligible impact, as does the flat-prior correction. When comparing our Gaussian process reconstruction with a parametric fit based on Eq. (2), the results are in good agreement. For the parametric fit, we find

[TABLE]

compared to Eq. (5) for the GP reconstruction. Although the results are similar, the posterior distribution for the GP has a broader tail at large optical depths. This can be traced to uncertainty in the reconstruction at high $z$. Unlike the parametric form, the GP reconstruction permits an early drop and rise in the neutral hydrogen fraction before $z>9$, allowing larger values of $\tau$.

3.3 Extended reionization histories

The most important modelling assumption relates to the boundary condition at $z_{\mathrm{max}}=15$. Increasing the maximum starting redshift of reionization to $z_{\mathrm{max}}=30$ opens up the possibility of extended reionization histories. Due to the lack of observations beyond $z>12$, the possibility that there is a long tail with a non-negligible ionization fraction at high redshifts, possibly driven by Population III (Pop III) stars (Qin et al., 2020; Wu et al., 2021), cannot be excluded using current $x_{\mathrm{HI}}(z)$ constraints alone.

However, models with high-$z$ contributions to $\tau$ can be constrained when Lyman-$\alpha$ observations are combined with kSZ or large-scale CMB polarization data (e.g. Qin et al., 2020; Wu et al., 2021; Cain et al., 2025a). We illustrate this by imposing a prior on the midpoint and duration of reionization, based on the recent kSZ trispectrum analysis from SPT (Raghunathan et al., 2024). They find $\Delta z_{50}<4.5$ (95%), when imposing a prior that reionization ends by $z=6$. This is similar to the constraint $\Delta z_{50}<4.1$ (95%) from the kSZ power spectrum analysis of Reichardt et al. (2021).

Based on the amber simulations (Trac et al., 2022; Chen et al., 2023), the average reionization signal in the kSZ trispectrum, on scales $50\leq L\leq 300$, can be modelled as (Raghunathan et al., 2024)

[TABLE]

where $\Delta z_{90}=z_{95}-z_{5}$, is defined similarly as $\Delta z_{50}$. This relationship implicitly assumes a monotonic reionization history, which is not necessarily the case in our GP reconstruction. To allow as much freedom as possible, we only require that $z_{5}$, $z_{\mathrm{mid}}$, and $z_{95}$ are uniquely defined, but otherwise impose no monotonicity condition, although doing so would only tighten the constraints. We then include a joint prior on $p(\Delta z_{90},z_{\mathrm{mid}})=p(LC_{L}^{KK})$, sampling from the joint two-dimensional marginalized posterior distribution, which we reproduce from Fig. 5 in Raghunathan et al. (2024).

The outcome of the GP reconstruction, when setting $z_{\mathrm{max}}=30$ and imposing the kSZ prior, is shown in Fig. 4. We see that the high-redshift tail is strongly constrained beyond $z>15$. The posterior distribution of $\tau$ in this case, shown as a short dashed line in Fig. 3, agrees well with our baseline case. We find

[TABLE]

which is still consistent with the baseline result of Eq. (5). This confirms that our fiducial set-up with $z_{\mathrm{max}}=15$ is reasonable. Consistently with Cain et al. (2025a), we find that large optical depths, $\tau=0.09$, are ruled out at more than $2\sigma$ when Lyman-$\alpha$ constraints on the beginning of reionization are combined with kSZ constraints on its duration, even when allowing $z_{\mathrm{max}}=30$.

The uniqueness condition on $z_{95}$ does play a role in this result. Including the $x_{\mathrm{HI}}(z)$ data ensures that a crossing of $x_{\mathrm{HI}}=0.95$ occurs around $z\approx 9-10$ with high probability. Requiring this crossing to be unique then prevents $x_{\mathrm{HI}}(z)$ from dropping below this value at higher redshifts. However, out set-up still allows non-standard reionization scenarios in which the IGM is ionized multiple times, as in double reionization scenarios (Cen, 2003; Furlanetto & Loeb, 2005; Mason et al., 2019b; Ahn & Shapiro, 2021; Tan, 2025), provided that $x_{\mathrm{HI}}<0.95$ remains true during the middle recombination period. The dashed line in Fig. 4 corresponds to such a double reionization history, which achieves an optical depth of $\tau=0.09$ and fits the $x_{\mathrm{HI}}(z)$ data at $z<12$, but is excluded by the combination of Lyman-$\alpha$ and kSZ data.

By construction, we have excluded the subset of double reionization models that achieve near complete recombination in the middle period. We still expect that such models can be constrained by the combination of Lyman-$\alpha$ and kSZ data, since the kSZ signal is significantly larger for double reionization models (McQuinn et al., 2005), but more work is needed to obtain reliable constraints.

3.4 Comparison with the literature

Most previous works included a CMB optical depth constraint as input, whereas we have deliberately excluded it. However, anticipating the discussion in the next section, we find that astrophysical constraints on $x_{\mathrm{HI}}(z)$ lead to very tight constraints on the reionization history that agree with constraints from large-scale CMB polarization measurements. These findings are consistent with previous works (e.g. Greig & Mesinger, 2017; Mason et al., 2019b; Paoletti et al., 2025; Sims et al., 2025).

Our baseline constraint, $\tau=0.0492^{+0.0014}_{-0.0030}$, is consistent with most previous analyses, but towards the lower end of the range. For instance, Mason et al. (2019b) found $\tau=0.053\pm 0.004$, based on a non-parametric reconstruction from Lyman-$\alpha$ data combined with a CMB optical depth constraint. We can also compare with the recent analysis of Paoletti et al. (2025). They found $\tau=0.0515_{-0.0028}^{+0.0014}$ from CMB data combined with $x_{\mathrm{HI}}(z)$ constraints, which is fully consistent with our result. When also including UV luminosity function data, they obtain an increased optical depth of $\tau=0.0542_{-0.0028}^{+0.0017}$, although this result is more model dependent. Another recent analysis by Sims et al. (2025) found $\tau=0.052_{-0.0018}^{+0.0016}$ from a joint analysis of CMB, Lyman-$\alpha$ and Lyman-$\beta$, and 21-cm constraints. In all cases, the error bars are much smaller than the typical uncertainty from the CMB, in agreement with our findings.

The slight decrease in $\tau$ compared to previous works could be explained by the inclusion of the latest $x_{\mathrm{HI}}(z)$ constraints at $z>8$. In particular, the data points from Umeda et al. (2025) help to pinpoint the beginning of reionization, reinforcing the preference for a late start and a low optical depth. Another contributing factor could be that the large-scale CMB polarization data that we excluded prefer slightly larger values of $\tau$, though this depends on the implementation of the likelihood (see Section 4.1.2).

4 Constraints from the CMB

In this section, we will consider CMB constraints on $\tau$. Since these data are less sensitive to the shape of the reionization history, we adopt the simpler parametric form of Eq. (2) for all CMB-based analyses that exclude direct $x_{\mathrm{HI}}(z)$ constraints. This allows us to explore a greater variety of models at reduced computational cost.

4.1 Results for $\Lambda$CDM

Let us begin by discussing the constraints from small-scale CMB data and then move on to large-scale CMB polarization.

4.1.1 Small-scale measurements

It is interesting to consider what constraints on $\tau$ are possible using only small-scale CMB measurements (Giarè et al., 2024a; Sailer et al., 2025; Jhaveri et al., 2025). When analyzing such data within $\Lambda$CDM, there is a strong degeneracy between the parameters $\tau$ and $\Omega_{\mathrm{m}}$. This is especially true when including information from CMB lensing (Jhaveri et al., 2025), as may be gleaned from the left panel of Fig. 5. The solid gray contours represent the constraints from small-scale Planck data alone. The black dashed contours show that these constraints are significantly tightened, in the direction orthogonal to the $\Omega_{\mathrm{m}}$-$\tau$ degeneracy, by adding the CMB lensing likelihood of Qu et al. (2025), which is based on Planck, ACT, and SPT data. This can be explained by the fact that CMB lensing measurements are sensitive to $\sigma_{8}\Omega_{\mathrm{m}}^{0.25}$, the previously mentioned degeneracy between $\tau$ and $A_{\mathrm{s}}$, and the dependence of the clustering amplitude, $\sigma_{8}$, on the amplitude of the primordial power spectrum, $A_{\mathrm{s}}$.

A similar conclusion holds when we add low-redshift galaxy clustering and lensing measurements, which are sensitive to $\sigma_{8}\Omega_{\mathrm{m}}^{0.5}$. We demonstrate this in the right panel of Fig. 5, where we add the $6\times 2$pt likelihood of Abbott et al. (2023a), which is based on DES, Planck, and SPT measurements. The degeneracy between $\Omega_{\mathrm{m}}$ and $\tau$ is again apparent from the black dashed contours, although with an expected rotation in the direction of the degeneracy.

To make progress, the matter density fraction can be further constrained using low-redshift probes. DESI BAO precisely determine $\Omega_{\mathrm{m}}=0.2975\pm 0.0086$ (DESI Collaboration et al., 2025a). The blue contours in the left panel of Fig. 5 show that the combination of small-scale Planck data, CMB lensing, and DESI BAO yields a large value of $\tau=0.094\pm 0.011$. These constraints are in tension with those from BAO, BBN, and $x_{\mathrm{HI}}(z)$ data, which were discussed extensively in Section 3 and are shown in purple in Fig. 5. The DES compilation of type 1a SNe prefers a significantly larger $\Omega_{\mathrm{m}}=0.352\pm 0.017$ (DES Collaboration et al., 2024). Moving along the CMB degeneracy, this higher $\Omega_{\mathrm{m}}$ gives rise to a smaller $\tau=0.062\pm 0.014$, as shown by the red contours. This constraint is compatible with the $x_{\mathrm{HI}}(z)$ measurement of $\tau$. However, the tension in $\Omega_{\mathrm{m}}$ between BAO and SNe is also apparent between the purple and red contours.

In the right panel of Fig. 5, we see better agreement among the different probes when including the $6\times 2$pt likelihood of Abbott et al. (2023a). This can be explained by the fact that the $6\times 2$pt analysis yields lower values of $S_{8}\equiv\sigma_{8}(\Omega_{\mathrm{m}}/0.3)^{0.5}=0.792\pm 0.012$ (Abbott et al., 2023a), compared to $S_{8}=0.828\pm 0.012$ for the CMB lensing analysis of Qu et al. (2025). As a result, the combinations with DESI BAO and DES SNe yield smaller values of $\tau$ that are more compatible with the $x_{\mathrm{HI}}(z)$ constraint. In the case of DES SNe, we find $\tau=0.047\pm 0.015$, with a central value close to that of the $x_{\mathrm{HI}}(z)$ analysis. It is important to note that, despite the widely reported tension in $S_{8}$ (e.g. Abdalla et al., 2022; Di Valentino et al., 2025), the results of the $6\times$2pt analysis of Abbott et al. (2023a) are consistent with the primary Planck CMB measurements (see also Chen et al. 2024; DES and KiDS et al. 2023; DES Collaboration et al. 2025; DESI Collaboration et al. 2024; Wright et al. 2025), justifying their combination.

In summary, while it is possible to obtain tight constraints on $\tau$ by combining small-scale CMB data with low-redshift probes, the degeneracy between $\tau$ and other cosmological parameters, particularly $\Omega_{\mathrm{m}}$ and $\sigma_{8}$, means that the results are sensitive to the choice of low-redshift probe. Later on, we will see that the results are also sensitive to the choice of cosmological model.

4.1.2 Large-scale measurements

Let us move on to the constraints from large-scale CMB polarization (lowE). We consider three different likelihoods based on Planck data: SimAll (Planck Collaboration et al., 2020a), SRoll2 (Delouis et al., 2019), and LoLLiPoP (Tristram et al., 2021, 2024). Compared to the SimAll likelihood used in the Planck 2018 analysis (Planck Collaboration et al., 2020b), the SRoll2 likelihood uses an improved version of the SRoll mapmaking algorithm that corrects for known instrumental effects. Both SimAll and SRoll2 use an empirical probability distribution constructed from simulated power spectra. LoLLiPoP uses a Gaussian likelihood function with a covariance matrix based on simulated spectra, handling non-Gaussianities with an approximation based on Hamimeche & Lewis (2008).

Although it is possible to constrain $\tau$ using lowE data alone, even when varying all $\Lambda$CDM parameters, we combine these likelihoods with DESI BAO and BBN constraints to break parameter degeneracies, which improves the constraint on $\tau$ by $40\%$. For instance, in the case of SRoll2, we find

[TABLE]

while adding BBN and BAO improves this to

[TABLE]

The results for the other two likelihoods are given in Table 2. The three likelihoods are mutually consistent and in good agreement with the Lyman-$\alpha$ constraint of Section 3. Compared to SimAll, the newer likelihoods produce tighter constraints that shift closer to the Lyman-$\alpha$ bound.

In Fig. 6, we present the marginalized posterior distributions of $\tau$ for the various CMB datasets we have considered. The three large-scale polarization likelihoods are shown in dark green, demonstrating the consistency between SimAll (dotted green), SRoll2 (solid green), and LoLLiPoP (dashed green). These measurements are also consistent with the $x_{\mathrm{HI}}(z)$ constraint, which is shown as a purple shaded band. The combination of DES SNe, small-scale Planck data, and the $6\times 2$pt likelihood of Abbott et al. (2023a), represented by the red long dashed line, also agrees well with the large-scale polarization and $x_{\mathrm{HI}}(z)$ bounds.

A second set of posteriors, centred on larger values of $\tau$, are shown in light blue. These constraints arise from the combination of DESI BAO, CMB lensing, and small-scale CMB data. Compared to the combination with primary CMB measurements from Planck (solid blue), we find very similar results with the alternative ACT-SPT dataset (dashed blue) and the most powerful P-ACT-SPT combination (dotted blue). See Table 2 for the numerical results. In all three cases, the constraints are in tension with the $x_{\mathrm{HI}}(z)$ bound (see below).

4.2 Results for extended cosmological models

In this section, we consider extensions of the $\Lambda$CDM model, including variations in the shape of the primordial power spectrum, the amplitude of the lensing potential, spatial curvature, neutrino content, and dark energy parameters. Table 3 provides constraints on $\tau$ for all cosmological models, obtained using three different approaches. The first is a non-parametric reconstruction of the reionization history, based primarily on $x_{\mathrm{HI}}(z)$ constraints. The second approach relies on constraints from small-scale Planck data, combined with CMB lensing and DESI BAO data. The third approach relies primarily on large-scale Planck polarization measurements.

We also provide constraints on $T=\tau\sqrt{\Omega_{\mathrm{m}}h^{2}/0.14}$, which can be precisely measured with $x_{\mathrm{HI}}(z)$ data, independently of probes of the expansion history. The final column in Table 3 lists the tension between the constraints on $T$ from BBN and $x_{\mathrm{HI}}(z)$ data and those from small-scale Planck, CMB lensing, and DESI BAO. The reason for computing the tension in $T$ is that this parameter can be constrained with $x_{\mathrm{HI}}(z)$ data, independently of DESI BAO. The resulting constraints are therefore uncorrelated with those from CMB and BAO data, which simplifies the calculations. We use the tensiometer package777https://github.com/mraveri/tensiometer to compute the level of tension, employing the parameter shift metric of Raveri et al. (2020) and Raveri & Doux (2021).

Before discussing each model in detail, we first make some general remarks. We observe that there is a $3.7\sigma$ tension in $T$ when assuming $\Lambda$CDM. The tension is reduced to $\sim 3\sigma$ for three models, for which the CMB and BAO constraints shift closer towards the Lyman-$\alpha$ bound. For three other models, the tension is reduced to below $2\sigma$. This is the case when the gravitational lensing potential is scaled by the $A_{\mathrm{lens}}$ parameter, when allowing for an effective cosmological neutrino mass parameter, $\sum m_{\nu,\mathrm{eff}}$, and when allowing for dynamical dark energy with an evolving equation of state (‘$w_{0}w_{a}$CDM’).

The most notable difference between the three approaches to constraining $\tau$ relates to their sensitivity to the choice of cosmological model. The distinctive signature of reionization on the large-scale polarization power spectrum means that the lowE constraints are largely insensitive to the choice of cosmological model. Similarly, the $x_{\mathrm{HI}}(z)$ constraints fix the reionization history, essentially independently of the cosmological model. By contrast, the constraints on $\tau$ from small-scale Planck, CMB lensing, and DESI BAO data do depend on the choice of cosmology, with shifts of up to $\sim 3\sigma$ from the $\Lambda$CDM value.

Some remarks are also in order about the validity of the $x_{\mathrm{HI}}(z)$ constraints in extended cosmological models. The dark pixel constraints are truly model-independent, but damping wing constraints typically rely on reionization simulations that assume a given cosmological model (e.g. Mesinger et al., 2016). In principle, there could be additional dependence on the cosmological model, not captured in our analysis. For instance, the coupling between cosmological parameters and galaxy formation processes could produce subtle effects (e.g. Elbers et al., 2025c; Pranjal et al., 2025). However, at least for the model variations considered in this work, these effects are expected to be subdominant compared to the uncertainty in subgrid reionization parameters. A more detailed investigation of such couplings will be left for future work.

In the following subsections, we will briefly comment on each cosmological model and assess the tension in $T$. In cases where the tension between CMB, BAO, and $x_{\mathrm{HI}}(z)$ data is resolved, we will combine the most powerful CMB dataset (P-ACT-SPT with lensing) and the most powerful BAO dataset (DESI DR2) with a dataset that tightly constrains the optical depth (Table 1) to obtain the strongest possible constraints on extensions of $\Lambda$CDM.

4.2.1 Running of the spectral index

In $\Lambda$CDM, the spectral index, $n_{\mathrm{s}}$ is independent of scale, but inflationary models predict departures from pure scale invariance (Escudero et al., 2016). The running of the spectral index is defined as $\alpha_{\mathrm{s}}=\mathrm{d}n_{\mathrm{s}}/\mathrm{d}\!\log k$. Table 3 shows that allowing for non-zero running does not significantly impact the constraints on $\tau$. The tension between the Lyman-$\alpha$ constraint and the CMB + BAO combination remains at the $4\sigma$-level.

4.2.2 Lensing excess

The CMB lensing anomaly refers to the preference of CMB data for a greater amount of gravitational lensing, either measured through the lensing smoothing effect on the primary CMB measurements or through the lensing reconstruction from the four-point correlation function, than the level that would be self-consistently predicted within a given cosmological model. It is commonly quantified using the $A_{\mathrm{lens}}$ parameter that scales the lensing potential (Calabrese et al., 2008; Renzi et al., 2018; Mokeddem et al., 2023).

By constraining the matter density fraction, BAO data play an important role in the $A_{\mathrm{lens}}$ anomaly. As discussed in Section 4.1.1, CMB lensing measurements impose a strong degeneracy between $\tau$ and $\Omega_{\mathrm{m}}$, such that the lower $\Omega_{\mathrm{m}}$ measured by DESI leads to a higher $\tau$ than the one preferred by Lyman-$\alpha$ constraints or large-scale CMB polarization measurements. Introducing the $A_{\mathrm{lens}}$ factor breaks the link between $\Omega_{\mathrm{m}}$ and $\tau$, making it possible to achieve a lower $\tau$ for a given $\Omega_{\mathrm{m}}$ through a value $A_{\mathrm{lens}}$ greater than unity.

Consequently, the tension in the optical depth can be entirely absorbed by adopting $A_{\mathrm{lens}}>1$. The tension in $T$ is reduced to $0.9\sigma$. However, combining the CMB dataset from P-ACT-SPT with DESI BAO and the Lyman-$\alpha$ constraints then leads to

[TABLE]

corresponding to a $3.7\sigma$ preference for $A_{\mathrm{lens}}>1$.

4.2.3 Curvature

The lensing anomaly is related to the preference of CMB data (including large-scale polarization) for a closed Universe with $\Omega_{\mathrm{k}}<0$, due to the correlation between $A_{\mathrm{lens}}$ and $\Omega_{\mathrm{k}}$ (Di Valentino et al., 2020). This can be seen as another expression of the preference of small-scale CMB data for an optical depth greater than the one allowed by lowE or Lyman-$\alpha$ constraints.

Allowing for non-zero spatial curvature has some impact on the CMB + BAO constraint on $\tau$, shifting the best-fitting value by more than $1\sigma$ towards lower optical depths (see Table 3). However, $\Omega_{\mathrm{k}}<0$ is no longer preferred once BAO data are included (Chen & Zaldarriaga, 2025). Indeed, the combination of DESI and CMB data is consistent with a flat Universe (DESI Collaboration et al., 2025c). Allowing for $\Omega_{\mathrm{k}}\neq 0$ therefore only reduces the tension with the Lyman-$\alpha$ bound to $2.7\sigma$.

4.2.4 New relativistic species

We also consider the possibility of additional light degrees of freedom, as parametrized by the effective number of neutrino species, $N_{\mathrm{eff}}$. The Standard Model prediction for the non-instantaneous decoupling of three neutrinos is $N_{\mathrm{eff}}=3.044$ (Froustey et al., 2020; Bennett et al., 2021; Binder et al., 2024). The data are fully consistent with the Standard Model prediction assumed in our baseline analysis. Consequently, allowing variations in $N_{\mathrm{eff}}$ does not significantly impact the constraints on $\tau$ and the tension in $T$ remains at $3.3\sigma$.

4.2.5 Neutrino masses

One of the most striking anomalies arising from the tension between CMB and BAO data is the cosmological preference for negative effective neutrino masses (Craig et al., 2024; Green & Meyers, 2025; Elbers et al., 2025b, a). As pointed out in these and other works (e.g. Loverde & Weiner, 2024; Lynch & Knox, 2025), this preference is also related to the $A_{\mathrm{lens}}$ and $\tau$ problems.

To investigate this anomaly, we adopt the effective neutrino mass parametrization of Elbers et al. (2025b), which allows for negative neutrino contributions to the energy density at late times, implemented by extending standard perturbative calculations. This mathematical extension allows us to test for signs of new physics or systematics that correlate with the cosmological neutrino mass parameter, $\sum m_{\nu,\mathrm{eff}}$. When we allow $\sum m_{\nu,\mathrm{eff}}<0$, the tension in $T$ is greatly reduced from $3.7\sigma$ to $0.8\sigma$. A significant part of this reduction is due to a $3\sigma$ shift of the best-fitting value of $\tau$ for the CMB + BAO combination towards lower optical depths. This justifies the combination of CMB, BAO, and Lyman-$\alpha$ constraints on reionization. For the most constraining CMB dataset that combines small-scale P-ACT-SPT data with CMB lensing, we find

[TABLE]

The marginalized error is $\sigma(\sum m_{\nu,\mathrm{eff}})=0.038\,$\mathrm{eV}. To our knowledge, this is the strongest cosmological bound on the neutrino mass, representing a $30\%$ improvement in precision compared to Elbers et al. ([2025a](https://arxiv.org/html/2508.21069v1#bib.bib46)). However, much like the extension with $A_{\mathrm{lens}}$, we have absorbed the tension between the underlying datasets into a new parameter that takes unphysical values. The tension with the minimal neutrino mass under the normal ordering, $\sum m_{\nu}>$0.059\text{\,}\mathrm{eV} (Esteban et al., 2024), exceeds $4\sigma$. Physically, the finding that $\sum m_{\nu,\mathrm{eff}}<0$ indicates a preference in the data for enhanced small-scale clustering and a late-time reduction in the matter density fraction, $\Omega_{\mathrm{m}}$, which represent a reversal of the cosmological effects of standard massive neutrinos. The connection with $\tau$ is then apparent from Fig. 5.

4.2.6 Dark energy

As a final extension, we consider variations of the dark energy model. Generalizing from $\Lambda$CDM, we first consider dynamical dark energy with a constant equation of state, $w=P/\rho$, where $w$ is a free parameter. Adopting this ‘$w$CDM’ model does not significantly shift the constraints on $\tau$. The tension between the Lyman-$\alpha$ bound and the BAO and CMB constraints remains at $3\sigma$. This is consistent with expectations, given that current BAO and CMB data do not prefer $w$CDM over $\Lambda$CDM (e.g. DESI Collaboration et al., 2025a).

Next, we consider models of dynamical dark energy with an evolving equation of state, described by the Chevallier-Polarski-Linder parametrization (Chevallier & Polarski, 2001; Linder, 2003), in which

[TABLE]

where $w_{0}$ and $w_{a}$ are free parameters. BAO measurements from the first data release of DESI (DESI Collaboration et al., 2025c) were consistent with $\Lambda$CDM, but preferred $w_{0}w_{a}$CDM when combined with CMB and SNe data. With the second data release of DESI, this preference increased to $2.8\sigma$-$4.2\sigma$, depending on the SNe dataset (DESI Collaboration et al., 2025a). See Lodha et al. (2025) for an extended dark energy analysis of the DESI results and Giarè et al. (2024b) for references to the recent literature.

Allowing for an evolving equation of state of dark energy allows the BAO and CMB combination to be more compatible with lower optical depths. The tension in $T$ is reduced to $1.2\sigma$. Fig. 7 shows the marginalized posterior distributions of $\tau$ for various data combinations, assuming $w_{0}w_{a}$CDM. It is fruitful to compare this figure with its $\Lambda$CDM counterpart (Fig. 6). First of all, we observe that the posteriors for the large-scale polarization measurements (dark green lines) are very similar to ones obtained under $\Lambda$CDM, demonstrating their robustness to changes in the assumed cosmological model. Both the Lyman-$\alpha$ bound (purple band) and the posterior from the combination of CMB, SNe, and $6\times 2$pt measurements (red long dashed line) are also similar to their $\Lambda$CDM counterparts, but with larger uncertainties. By contrast, the posteriors from the combinations of DESI BAO, CMB, and CMB lensing data (light blue lines) are shifted towards lower $\tau$, such that all constraints are compatible with a rapid and late reionization history.

Given that the tension in $T$ is significantly reduced, we combine CMB, BAO, and $x_{\mathrm{HI}}(z)$ data to constrain the dark energy parameters. We also include the DES SNe dataset, since the constraints from SNe and BAO are compatible within $w_{0}w_{a}$CDM (see Fig. 8). Using the CMB dataset from P-ACT-SPT, we obtain the marginalized constraints

[TABLE]

which are close to the values reported in DESI Collaboration et al. (2025c) for the combination of DESI BAO, CMB, and DES SNe. Fig. 8 illustrates how the addition of $x_{\mathrm{HI}}(z)$ data pulls the constraints away from $\Lambda$CDM. We estimate the preference for $w_{0}w_{a}$CDM to be $4.5\sigma$ when combining P-ACT-SPT, DESI BAO, DES SNe, and $x_{\mathrm{HI}}(z)$ data.888We employ the $\Delta\chi^{2}_{\mathrm{MAP}}$ metric used by DESI Collaboration et al. (2025c, a). This is an increase compared to the $4.2\sigma$ preference reported in DESI Collaboration et al. (2025a). Compared to that analysis, we have replaced their CMB dataset with the more constraining P-ACT-SPT dataset (including lensing), but excluded large-scale CMB temperature and polarization data, which we replaced with our $x_{\mathrm{HI}}(z)$ likelihood. Including large-scale temperature and polarization data would likely further increase the preference for $w_{0}w_{a}$CDM.

Since we deliberately excluded lowE data, these results are invulnerable to any large-scale CMB systematics. On the other hand, one must now contend with the possibility of systematics in the analysis of Lyman-$\alpha$ damping wings, particularly at the highest redshifts (Mesinger & Furlanetto, 2008; Keating et al., 2024; Umeda et al., 2024; Huberty et al., 2025).

5 Conclusion

Motivated by anomalies in the latest cosmological data, we have presented new constraints on the reionization optical depth, $\tau$, using several independent methods. Our strongest bound derives from a non-parametric reconstruction of the reionization history from Lyman-$\alpha$ and Lyman-$\beta$ constraints on the evolution of the volume-averaged neutral hydrogen fraction, $x_{\mathrm{HI}}(z)$, between $5<z<12$. The constraints mainly derive from dark pixel fraction constraints at low redshifts and damping wing constraints at high redshifts (see Table 1 for references).

Our results point to a late and rapid reionization history, with a midpoint at $z_{\mathrm{mid}}=7.00^{+0.12}_{-0.18}$ and a duration of $\Delta z_{50}=1.12^{+0.12}_{-0.29}$. We obtain an optical depth, $\tau$, that is perfectly consistent with Planck measurements of the large-scale CMB polarization power spectrum. The strong constraining power of quasar damping wings and the consistency between damping wing, dark pixel fraction, and CMB constraints on the optical depth agree with the findings of earlier works (e.g. Greig & Mesinger, 2017). Our baseline constraint, $\tau=0.0492^{+0.0014}_{-0.0030}$, is consistent with previous analyses (Mason et al., 2019b; Paoletti et al., 2025; Sims et al., 2025), but incorporates the latest data and is independent of the CMB.

Under the assumption of the cosmological $\Lambda$CDM model, BAO data from DESI, combined with CMB lensing and small-scale CMB data from Planck (or the independent ACT and SPT datasets) imply a larger optical depth, $\tau=0.094\pm 0.011$. Such a large optical depth provides a solution to the tension between DESI BAO and CMB constraints that preserves the $\Lambda$CDM model (Sailer et al., 2025; Jhaveri et al., 2025). However, it would have a number of consequences that we summarize below:

• •

The large optical depth is in $3.7\sigma$ tension with our Lyman-$\alpha$ bound, unless we relax the assumption about the maximum starting redshift of reionization and allow for reionization scenarios with a significant contribution to $\tau$ from $z>15$, as in models in which the IGM is reionized multiple times (Cen, 2003; Furlanetto & Loeb, 2005; Mason et al., 2019b; Ahn & Shapiro, 2021; Tan, 2025) or with a long high-redshift tail from Pop III stars (Qin et al., 2020; Wu et al., 2021).

• •

We have shown that such early reionization models are still in tension with constraints on the midpoint and duration of reionization from the kSZ effect (Raghunathan et al., 2024; Cain et al., 2025a), but further work is needed to study the kSZ effect in non-monotonic models.

• •

It would also conflict with Planck lowE measurements. There is no evidence for systematic errors in these measurements, but Giarè et al. (2024a) give some reasons why the possibility cannot be dismissed. Jhaveri et al. (2025) offer another way of relaxing the lowE constraint on $\tau$, through a change in the curvature fluctuation power spectrum on large scales. However, this would no longer be a standard $\Lambda$CDM scenario.

• •

A large optical depth would also be in tension with the constraint from the combination of small-scale CMB data with type 1a supernovae and $6\times 2$pt galaxy lensing, galaxy clustering, and CMB lensing measurements from DES, Planck, and SPT. Although the DES SNe results have been challenged based on differences with the Pantheon+ compilation (Efstathiou, 2025a), these differences can be largely explained by well-motivated methodological choices (Vincenzi et al., 2025). Galaxy lensing and clustering measurements have long been associated with a tension in $S_{8}$ (e.g. Abdalla et al., 2022; Di Valentino et al., 2025), but some recent results (including the ones used in this paper) show consistency with the CMB (Abbott et al., 2023a; Chen et al., 2024; DES and KiDS et al., 2023; DES Collaboration et al., 2025; DESI Collaboration et al., 2024; Wright et al., 2025). Future constraints on $S_{8}$ and $\Omega_{\mathrm{m}}$ from DESI, Euclid, and LSST could adjudicate this issue.

On the other hand, with an optical depth as low as $\tau=0.0492^{+0.0014}_{-0.0030}$, the constraints from DESI BAO and the CMB are in tension. This gives rise to several anomalies or apparent departures from $\Lambda$CDM. We have shown that the Lyman-$\alpha$ constraint on $\tau$ strengthens the evidence for the gravitational lensing excess, unphysical neutrino masses, and dynamical dark energy with an evolving equation of state. It remains an open questions whether these hints truly point to a new cosmological paradigm (Ishak & Medina-Varela, 2025; Leauthaud & Riess, 2025; Efstathiou, 2025b). Future CMB experiments could significantly reduce the uncertainty in $\tau$ (LiteBIRD Collaboration et al., 2023). Analyses of the 21-cm signal, the patchy kSZ effect, and the Lyman-$\alpha$ forest provide other promising avenues. Ongoing observations probing the early stages of reionization will bolster the robustness of astrophysical constraints on $\tau$, which may prove pivotal in testing our currently favoured cosmological model.

Acknowledgements

I thank Carlos Frenk, Adrian Jenkins, Baojiu Li, and Silvia Pascoli for encouragement and useful discussions. I acknowledge STFC Consolidated Grant ST/X001075/1. This work used the DiRAC@Durham facility managed by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The equipment was funded by BEIS capital funding via STFC capital grants ST/K00042X/1, ST/P002293/1 and ST/R002371/1, Durham University and STFC operations grant ST/R000832/1. DiRAC is part of the National e-Infrastructure.

Data Availability

We publicly release our likelihood code at https://github.com/wullm/reionlik. The data will be made available upon request.

Appendix A Technical details

To map the Gaussian process to the interval $(0,1)$, we use a wrapping function $F\colon\mathbb{R}\to(0,1)$. For this purpose, we choose the cumulative distribution function of the standard generalized normal distribution with shape parameter $\beta=8$,

[TABLE]

where $\gamma(s,x)=\int_{0}^{x}t^{s-1}e^{-t}\mathrm{d}t$ is the unnormalized incomplete lower gamma function. This function is very close to linear over most of its range, ensuring that $\chi(z)\approx x_{\mathrm{HI}}(z)$.

Inspired by the Gaussian process reconstruction of the dark energy equation of state, $w(z)$, of Lodha et al. (2025), we impose a generalized inverse Gaussian prior on the quantity $100\ell$,

[TABLE]

where $p=3$ and $b=1.5$ and where $K_{p}$ is the modified Bessel function of the second kind. This choice penalizes functions with too much freedom at low values of $\ell$ without oversampling nearly linear functions with high values of $\ell$.

Similar to Lodha et al. (2025), we introduce a latent variable corresponding to the value of $x_{\mathrm{HI}}(z_{\mathrm{lat}})$ at $z_{\mathrm{lat}}=7$. We then generate the Gaussian process conditional on this variable, which accelerates convergence. We perform our Gaussian process regression using a modified version of the CLASS code (Lesgourgues, 2011) that allows for a dense sampling of the free electron fraction, $x_{\mathrm{e}}(z)$, over the redshift range $[z_{\mathrm{min}},z_{\mathrm{max}}]$.