© The Authors 2025

1. Introduction

The epoch of reionisation (EoR) represents a critical phase in cosmic history during which the intergalactic medium transitioned from a neutral to ionised state due to the emergence of the earliest sources of ionising light. This process is estimated to have occurred between redshifts z ∼ 20 and z ∼ 6, corresponding to approximately 150 million to 1 billion years after the Big Bang (Barkana & Loeb 2001; Loeb & Furlanetto 2013), during which the Universe transitioned from the so-called Cosmic Dark Ages to the highly structured and luminous cosmos we observe today. Understanding the EoR is crucial for cosmic structure formation theories and offers insight into the properties of the first star and galaxy populations in the early Universe.

Observational constraints from the cosmic microwave background (CMB), particularly measurements of the optical depth, τ, from the Planck satellite, now yielding τ < 0.06, suggest that reionisation was a relatively recent process likely driven by the first galaxies and possibly by population III stars (Planck Collaboration Int. XLVII 2016). Additional insights come from high-redshift quasar spectra, which exhibit Gunn-Peterson troughs indicative of a fully ionised intergalactic medium by z ∼ 5.4 (Venemans et al. 2013; Becker et al. 2015; Bosman et al. 2022). The evolution of Lyman-alpha emitters also provides evidence of a patchy and inhomogeneous reionisation process (Mason et al. 2018). Recent observations from the James Webb Space Telescope (JWST; Gardner et al. 2006) have dramatically improved our understanding of galaxy populations during the EoR, revealing an unexpectedly high abundance of ultraviolet-luminous galaxies at redshifts z > 9 (Harikane et al. 2023; Bouwens et al. 2023). These observations extend the ultraviolet luminosity function (UVLF) out to z ∼ 13, particularly at the bright end. Moreover, JWST has uncovered surprisingly massive and evolved galaxies at z > 10, challenging existing models of early galaxy formation and suggesting a more rapid and efficient assembly of stellar mass than previously anticipated (Labbé et al. 2023). If these galaxies formed early and were efficient producers of ionising radiation, the reionisation of the Universe may have progressed more rapidly than inferred from CMB data, potentially implying a higher optical depth, up to τ ≳ 0.07 (Muñoz et al. 2024). Reconciling these findings with CMB-based constraints is essential for developing a coherent picture of the timeline and sources of cosmic reionisation.

In this context, the present work aims at providing a comprehensive study of the constraints on the EoR using CMB data, which encode imprints of the ionisation history through its polarisation and temperature anisotropies. We explore a broad range of reionisation models – from simple parametric forms to more flexible non-parametric descriptions – to evaluate the impact of model assumptions on the inferred constraints. Our analysis focuses primarily on the reionisation optical depth, τ, but also extends to the reconstruction of the ionisation history as a function of redshift, x_e(z). We pay particular attention to the implicit priors introduced by each model on the reionisation history and key parameters such as τ.

The Planck satellite has provided some of the mostdetailed CMB measurements to date on the full sky, enabling stringent constraints on the parameters of our cosmological model, including the reionisation history (Planck Collaboration Int. XLVII 2016). In this work, we use the latest data from the Planck mission – Public Release 4 (PR4) – to constrain various reionisation models. The PR4 dataset represents an improvement in data processing and error analysis (Planck Collaboration Int. LVII 2020), offering more precise measurements on large scales, especially in polarisation and tighter constraints on cosmological parameters (Rosenberg et al. 2022; Tristram et al. 2024).

This manuscript proceeds by presenting first in Sect. 2 the range of considered reionisation models to be confronted to our selected datasets. The latter are described in Sect. 3, alongside our methodology for deriving cosmological constraints. Our results on the reionisation optical depth are detailed in Sect. 4 and the interpretation in terms of the reionisation history (and other derived quantities) is discussed in Sect. 5. We summarise our findings and perspectives in Sect. 6.

2. Modelling and probing the reionisation history

Previous measurements of the CMB anisotropies by the Planck satellite (Planck Collaboration Int. XLVII 2016) have provided substantial advancements in our understanding of the EoR, improving our estimates of its timing and duration. Most notably, these observations have been key to constraining the optical depth due to Thomson scattering, usually denoted τ, a critical parameter in the ΛCDM cosmological model:

$$\begin{array}{c}\hfill \tau ={\displaystyle {\int }_0^{z_{\max }}}{\sigma }_{\mathrm{T}}{n}_{\mathrm{e}}(z)\frac{\mathrm{d}r}{\mathrm{d}z}\mathrm{d}z,\end{array}$$(1)

where σ_T is the Thomson cross-section, n_e(z) is the volume-averaged free electron (proper) number density, z_max is an arbitrary redshift upper bound (high enough to encompass the whole EoR), and dr/dz is the line-of-sight (proper) distance per unit redshift (throughout we assume that the Universe is flat, in which case dr/dz = c/(1 + z)/H(z), with H(z) the expansion rate).

Thomson scattering between the CMB photons and free electrons generates linear polarisation from the quadrupole moment of the CMB radiation field at the scattering epoch. Such scattering occurs whenever a significant amount of free electrons is present in the Universe, namely at recombination as well as during the EoR. For the latter, the probability of scattering is directly linked to τ (namely 1 − e^−τ). Such rescattering of the CMB photons during reionisation then generates additional polarisation anisotropies at large angular scales.

This signal, appearing primarily as a low-ℓ bump in the E-mode polarisation angular power spectrum, is linked to the horizon size at the “new last-rescattering surface” and thus depends on the redshift of reionisation. The amplitude of the bump scales roughly as the square of the reionisation optical depth, τ, but its exact shape encodes information about the evolution of the ionisation fraction as a function of redshift (Zaldarriaga et al. 1997; Kaplinghat et al. 2003). At the same time, the rescattering of CMB photons also reduces the overall amplitude of all scalar power spectra, scaling as e^−2τ. This introduces a correlation between τ and the primordial amplitude of the scalar perturbations, A_s, when inferring their values from observational data, as the latter is sensitive to the overall amplitude of the power spectra which scales as A_se^−2τ. Additionally, the kinetic Sunyaev-Zel’dovich (kSZ) effect, which arises from the scattering of CMB photons off moving electrons in ionised regions, provides complementary constraints on the duration and patchiness of reionisation through small-scale CMB anisotropies (McQuinn et al. 2005; Mesinger et al. 2012; Gorce et al. 2020). However, as its modelling requires additional assumptions and lies beyond the scope of this work, we do not include kSZ-based constraints in our analysis.

Accurate measurements of τ through the CMB are challenging due to foreground contamination and instrumental systematic uncertainties. These effects are particularly significant on large angular scales, from where most of the constraints on τ originate. The first estimates by WMAP suggested a high τ (Kogut et al. 2003) corresponding generally to an early start for reionisation processes, but subsequent measurements by Planck adjusted τ downwards (see Fig. 1), thanks to a re-assessment of the dust contamination in particular (Planck Collaboration Int. XLVII 2016; Planck Collaboration VI 2020). This shift restored consistency with other astrophysical data, suggesting a later reionisation not requiring exotic early-time ionising sources. Further refinements in the analysis of the Planck PR4 yielded an even better sensitivity in polarisation measurements, achieving a precision now reaching ∼2.5 times the cosmic-variance limited uncertainty (Tristram et al. 2024).

Fig. 1. Mean values (filled circles) and 68% credible intervals (horizontal bars) obtained from posterior distributions of τ, derived from CMB power spectra measurements including: WMAP 2013 (Hinshaw et al. 2013), Planck PR1 + WMAP 2014 (Planck Collaboration XV 2014), Planck PR2 + WMAP 2016 (Planck Collaboration XIII 2016), Planck-LFI 2016 (Planck Collaboration XIII 2016), Planck-HFI 2016 (Planck Collaboration Int. XLVII 2016), WMAP+LFI 2018 (Weiland et al. 2018), Planck PR3 2020 (Planck Collaboration VI 2020), Planck 2020 (Pagano et al. 2020), WMAP+LFI 2022 (Paradiso et al. 2023), Planck PR4 2022 (Planck Collaboration Int. LVII 2020). The shaded grey region highlights the cosmic variance limit on the measurement of τ, centred on the value from Planck PR4 2024 (Tristram et al. 2024).

However, it is important to highlight that the aforementioned CMB measurements were derived under the assumption of a simple reionisation history, modelled as a single-step transition in the ionisation fraction, x_e(z), using a hyperbolic tangent function. While alternative approaches – such as principal component analysis (Heinrich & Hu 2021), FlexKnot (Planck Collaboration VI 2020), and Poly-reion (Paoletti et al. 2020) – have been explored using Planck 2018 data, it should be emphasised that the low-ℓ EE likelihood presented in Planck Collaboration V (2020) is an empirical likelihood based on simulations constructed within a fixed ΛCDM framework assuming a hyperbolic tangent reionisation. As such, it is not suitable for constraining more extended or flexible reionisation scenarios.

In the present work, we set out to test various reionisation scenarios by constraining the evolution of the fraction of ionised matter as a function of redshift, thereafter x_e(z), defined as the ratio of n_e(z) over the hydrogen number density n_H(z) and ranging between 0 ≤ x_e(z)≤(1 + 2f_He) (where f_He = n_He/n_H is the primordial helium-to-hydrogen nucleon ratio). This process involves confronting observational data to models with various levels of sophistication, ranging from instantaneous to extended, from symmetric to asymmetric, and from physically motivated to model-independent. The variety of reionisation models comes from the recognition that different, still mostly unknown mechanisms – ranging from the first stars and galaxies to quasars – may have their own distinct signatures, at distinct moments on the reionisation timeline.

In the following, we briefly describe the suite of parametric and non-parametric models that we examine, and provide a visual summary in Fig. 2. A list of all model parameters and the corresponding prior choices are summarised in Table 1. We note that, unless specified otherwise, we include in all models a second reionisation of helium between redshift 3 and 4 (constructed using a simple tanh function of width $\mathrm{\Delta }{z}_{\mathrm{re}}^{\mathrm{He}}=0.5$ and central redshift $z_{\mathrm{re}}^{\mathrm{He}}=3.5$)¹, which is not mentioned in the model descriptions. Finally, in all our considered models, we set the beginning of reionisation at redshift 25 at the earliest (sufficiently deep to allow for an early start) and ensure that the hydrogen and first helium reionisation are ended by redshift 5 in accordance with astrophysical constraints (see references in Sect. 1).

Fig. 2. Graphical summary of the reionisation models considered in the present study, illustrating with coloured arrows and lines their respective degrees of freedom in modelling the history of the ionising fraction. The models are grouped by their general approach: hyperbolic tangents, power-laws, step functions, interpolations, bins, and PCA. The models are described in detail in Sect. 2.

2.1. Hyperbolic tangent models

The classic hyperbolic tangent model features a rapid transition of the Universe from a neutral to ionised state, and has been extensively utilised in CMB analyses as the default model implemented most notably in widespread Boltzmann codes (e.g. CAMB, Lewis 2008). It is parametrised by a midpoint redshift, z_re, and a duration, Δz_re, for the transition, as

$$\begin{array}{c}\hfill x_{\mathrm{e}}(z)\propto \frac{1}{2}[1+tanh\left(\frac{1}{\mathrm{\Delta }{z}_{\mathrm{re}}}\frac{{(1+z_{\mathrm{re}})}^{\gamma }-{(1+z)}^{\gamma }}{\gamma {(1+z_{\mathrm{re}})}^{\gamma -1}}\right)],\end{array}$$(2)

where the exponent, γ, is typically set to 1.5. We consider two forms of this model: one with a fixed duration of Δz_re = 0.5 and the other where the duration is allowed to vary (see the upper-left panel in Fig. 2), enhancing the flexibility of the model to account for more or less extended reionisations, although those rarely reflect actual expected physical evolutions. We denote these two models tanh and tanhdz, respectively. We note that tanhdz is unique in our suite of models, as it is the only one that allows for unfinished reionisations at z = 0.

2.2. Power-law models

A significant concern regarding the use of the hyperbolic tangent as a model relates to its lack of agreement with state-of-the-art hydrodynamical simulations of the EoR (Aubert et al. 2015; Seiler et al. 2019; Puchwein et al. 2023; Meriot & Semelin 2024). Indeed, these simulations tend to reveal a more intricate picture, with for instance a reionisation that starts slowly with the formation of the first stars and culminates abruptly. A more physically motivated parameterisation of the EoR was first proposed and explored in Douspis et al. (2015) and consequently used in the EoR analysis of the Planck data (Planck Collaboration Int. XLVII 2016), albeit in a simplified form, which we adopt here:

$$\begin{array}{c}\hfill x_{\mathrm{e}}(z)=\{\begin{array}{cc}(1+f_{\mathrm{He}})\hfill & \text{for}z\le z_{\mathrm{end}},\hfill \\ (1+f_{\mathrm{He}}){\left(\frac{z_{\mathrm{beg}}-z}{z_{\mathrm{beg}}-z_{\mathrm{end}}}\right)}^{\alpha }\hfill & \text{for}z>z_{\mathrm{end}}.\hfill \end{array}\end{array}$$(3)

This model, by design, incorporates an element of asymmetry absent from the hyperbolic tangent via the exponent, α, and the choice of the redshift of the start (z_beg) and end (z_end) of reionisation. It can accommodate a diverse range of astrophysically-plausible scenarios, which may include either a protracted onset prior to a swift completion or the inverse sequence. Two variants of the model are also considered, where the start redshift is either fixed to z_beg = 20 (as in Planck Collaboration Int. XLVII 2016) or free (see the upper-central panel in Fig. 2). We denote these models 2-pow and 3-pow, respectively.

2.3. Step-based models

We also explore in parallel a set of models based on step-like, sharp (Δz = 0.1) hyperbolic tangent transitions (see the upper-right panel in Fig. 2). More specifically, we adopt two models constructed out of 2 and 3 transitions, where the position of each transition is allowed to vary, as well as the height of the plateaus between transitions. This results in two models with 3 and 5 degrees of freedom, respectively, which we dub 2-steps (as used in Watts et al. 2020) and 3-steps in the following.

2.4. Interpolation-based models

Another approach we consider here involves flexible histories, namely via interpolation-based models in which the ionisation fraction is interpolated across a series of user-specified {z_i, x_e, i} points, called “knots” in the following. Such a model is simple yet versatile, allowing for the construction of a wide variety of reionisation trajectories by adjusting the specified points. We consider two types of interpolation: a simple linear interpolation model, and the “flexknot” model of Millea & Bouchet (2018). The latter employs Piecewise Cubic Hermite Interpolating Polynomials (PCHIP), which ensure monotonic transitions between nodes and provide a smooth reionisation history (see lower-left panel in Fig. 2). We consider models with 2 and 3 knots, and we note that the ionised fraction of the lowest and highest redshift nodes are always fixed, respectively, to 1 + f_He (to ensure reionisation is finished) and to its residual value from recombination (for continuity). We denote these models 2-knots and 3-knots for the linear interpolation, and 2-flex and 3-flex for the PCHIP interpolation, respectively.

2.5. Binned models

One of the most general model-independent approaches to study the reionisation history is to only consider the average ionisation fraction over a succession of redshift intervals. In practice, this can be achieved by adopting a binned history of the ionisation fraction, which we implement through a series of sharp (Δz = 0.1) consecutive hyperbolic tangent transitions with constant plateaus in-between, each corresponding to distinct, fixed redshift intervals. This discretised approach enables notably the investigation of non-monotonic reionisation histories, allowing for periods of both increased and decreased ionisation fraction within the EoR, but also more generally allows for a model-independent exploration, mostly independent of any assumption or model prior.

In the present study, we settle on a 10-bin model, spanning z = 5 to z = 25 and with a bin of width 2 in redshift – our model with the seemingly highest flexibility in our analysis, which we dub 10-bins in the following.

2.6. Principal component analysis model

The principal component analysis (PCA) model offers a generic approach to reconstruct the reionisation history from CMB data, without assuming any functional form but instead using a set of complex, data-driven functions. It involves decomposing the ionisation fraction history, x_e(z), into a set of orthonormal eigenfunctions or principal components (Hu & Holder 2003), which are specifically designed to maximise the sensitivity of the model to the underlying data.

In the PCA model, the ionisation fraction is then written as

$$\begin{array}{c}\hfill x_{\mathrm{e}}(z)=x_{\mathrm{e}}^{\mathrm{fid}}(z)+{\displaystyle \sum _i}{m}_i{S}_i(z),\end{array}$$(4)

where $x_{\mathrm{e}}^{\mathrm{fid}}(z)$ denotes an arbitrary chosen fiducial ionisation history, m_i represents the amplitudes of the principal components, and S_i(z) are the principal component templates. These templates are constructed to satisfy several properties, namely: forming a special orthonormal basis over the redshift range of interest, minimising the covariance of the amplitude parameters given cosmic-variance-limited CMB data, and being ordered by increasing uncertainty, ensuring that the most constrained modes are considered first. This construction allows the PCA model to flexibly adapt to various possible reionisation histories without being overly prescriptive. This model has been applied to WMAP (Mortonson & Hu 2008) in the past, as well as to Planck data (Heinrich et al. 2017; Heinrich & Hu 2021).

In the present work, we built a PCA model with a set of basis functions spanning the redshift range z = 5 to z = 25, and a fiducial model $x_{\mathrm{e}}^{\mathrm{fid}}=0.15$ (corresponding to τ_fid = 0.066, reasonably close to current estimates). The components are derived from a Fisher matrix-based approach, assuming a Planck best-fit cosmology for an all-sky, cosmic variance-limited experiment. By construction, due to the finite number of components used, the ionisation fraction, x_e(z), reconstructed from the PCA is not inherently restricted to physical bounds (i.e. between 0 and 1 + f_He). We keep only the first 5 components which according to Mortonson & Hu (2008) are sufficient to describe all the information on x_e carried by C_ℓ^EE up to the cosmic variance limit. During Markov chains Monte Carlo (MCMC) exploration, we impose priors to the amplitude of each mode following the prescription of Mortonson & Hu (2008), ensuring that reionisation histories with unphysical optical depths over a broad redshift range are treated as unphysical models². We denote this model pca-5 in the following.

3. Data and methodology

The following section outlines the specific datasets and the computational tools used in the present analysis to test the reionisation histories described in the previous section.

3.1. Data

The primary dataset used in this analysis is derived from the latest and final Planck PR4 maps. These maps have been processed through the NPIPE pipeline (Planck Collaboration Int. LVII 2020) which includes data from the Planck Low-Frequency Instrument (LFI) and the High-Frequency Instrument (HFI) in both temperature (T) and polarisation (P). The NPIPE processing pipeline improves upon previous releases by incorporating data from the repointing periods and by providing lower noise levels and reduced systematics across all angular scales, enhancing the internal consistency between frequencies. Additionally, the PR4 includes a comprehensive set of “End-to-End” Monte Carlo simulations aiming to help characterise potential biases and uncertainties associated with the data processing.

To ensure unbiased estimates of the angular power spectra, cross-correlations of two independent splits of the data – referred to as detsets – were performed, both across frequencies and across T and P signals. Detset maps comprise subsets of detectors that are combined to produce maps with nearly independent noise properties, crucial for robust cosmological parameter estimation.

In order to extract cosmological information from the Planck PR4 data, we split the dataset in two likelihoods: one for the large angular scales (low ℓ) and the other for small angular scales (high ℓ). The details of these respective likelihoods are briefly described in the following.

3.1.1. LOLLIPOP

LOLLIPOP (Low-ℓ Likelihood on Polarised Power-spectra) is a likelihood focused on the low multipoles of the CMB polarisation, including the EE, BB, and EB spectra. Initially applied to Planck PR3 data for reionisation studies, LOLLIPOP has been upgraded to PR4 as described in Tristram et al. (2021, 2022). LOLLIPOP uses cross-power spectra calculated from component-separated CMB detset maps processed by the Commander code (Eriksen et al. 2008) using all available polarised frequencies from 30 to 353 GHz.

This likelihood utilises the Hammimeche–Lewis approximation, modified to apply to cross-power spectra (Hamimeche & Lewis 2008; Mangilli et al. 2015). This method ensures a nearly Gaussian distribution of a transformed variable, derived from the measured and model power spectra. The likelihood function incorporates the covariance matrix derived from Monte Carlo simulations, ensuring that uncertainties from the CMB sample variance, statistical noise, and systematics are accurately propagated. Finally, the likelihood is marginalised over the unknown true covariance matrix (as proposed in Sellentin & Heavens 2016) in order to propagate the uncertainty in the estimation of the covariance matrix caused by a limited number of simulations.

In this work, we use LOLLIPOP to cover the multipole range from ℓ = 2 to ℓ = 30 for the EE spectrum, providing crucial constraints on the reionisation optical depth, τ.

3.1.2. HILLIPOP

HILLIPOP (High-ℓ Likelihood on Polarised Power-spectra) is designed for analysing the small-scale angular power spectra (ℓ > 30) from the Planck PR4 data. HILLIPOP uses a multi-frequency Gaussian likelihood approximation, incorporating sky maps at 100, 143, and 217 GHz, with the 353 GHz channel used to model dust contamination. This likelihood has been updated to V4.2 for the PR4 data release, utilising a larger sky fraction (up to 75%) and refined models for foreground components such as dust and point sources (Tristram et al. 2024).

Small-scale power spectra in TT, TE, and EE are computed using Xpol (an extension to polarisation of the Xspect pseudo-C_ℓ method, described in Tristram et al. 2005), which corrects for beam effects and normalises the power spectra with a mode-mixing matrix. HILLIPOP incorporates a comprehensive semi-analytic covariance matrix that includes contributions from cosmic variance, instrumental noise, and the residualforegrounds.

3.1.3. Combining likelihoods

By combining LOLLIPOP and HILLIPOP, we leverage the full range of Planck data and make full use of the PR4 improvements. We also incorporate the Planck Commander low-ℓ temperature likelihood in our analysis, which remains based on PR3 as no significant improvements were expected in the PR4 intensity maps for this range of multipoles. The synergy between these likelihoods enables a robust estimation of both the reionisation parameter, τ, from the large scales and the other cosmological parameters, in both temperature and polarisation. On the other hand, Appendix A demonstrates that incorporating information from the CMB-derived gravitational lensing power spectrum, while tightening the constraints on the amplitude of primordial scalar fluctuations (A_s), has no effect on the reconstruction of the reionisation optical depth despite the known correlation between the two parameters. Consequently, the Planck CMB lensing likelihood is not included in our analysis.

Unless explicitly mentioned otherwise, we sampled the likelihoods for the free reionisation parameters along with the five parameters of the ΛCDM cosmological model (the physical baryon density, Ω_bh², the physical density of cold dark matter, Ω_cdmh², the amplitude of primordial scalar fluctuations, A_s, the scalar index, n_s, and the Hubble constant, H₀), as well as the nuisance parameters associated with each likelihood. The latter are related to instrumental calibration and residual foreground modelling. We adopted the same priors on those nuisance parameters as described in Tristram et al. (2024). Additionally, we assumed three degenerate neutrinos of fixed total mass 0.06 eV.

3.2. Tools

Throughout our analysis, theoretical predictions for the CMB anisotropy power spectra were computed using the CLASS Boltzmann code (Blas et al. 2011). We developed a custom-modified version of CLASS³, which incorporates all the reionisation scenarios considered in this study. The only exception is the pca-5 model, for which we use a modified version of CAMB (Lewis et al. 2000) as it allows for the computation of CMB power spectra even in cases where the ionisation fraction temporarily goes negative.

To perform cosmological inference, we used the ECLAIR suite of codes⁴, designed as a versatile tool for sampling the posterior distribution of model parameters. ECLAIR (first introduced in Ilić et al. 2021) integrates cutting-edge datasets, efficient affine-invariant ensemble MCMC algorithms (via emcee and zeus, Foreman-Mackey et al. 2013; Karamanis et al. 2021; Karamanis & Beutler 2021), and interfaces with the CAMB and CLASS Boltzmann solvers, or any custom modifications of them. The pipeline is highly parallelised, leveraging multiprocessing on a single machine or distributing computation with MPI, proving advantageous for the computationally demanding tasks associated with parameter space exploration. ECLAIR also includes a dedicated minimiser that combines ensemble sampling with simulated annealing, efficiently targeting the global maximum of any posterior. The convergence and robustness of the analyses are further supported by integrated scripts that facilitate validation of the estimated marginalised posteriordistributions and credibility intervals, and take advantage of the capabilities of the getdist package (Lewis 2019) for outputting results.

We paid close attention to the convergence of the chains, especially for the more complex, weakly constrained models that exhibit very long correlation lengths (up to hundreds of steps). We ensured effective sample sizes of at least a few thousands (after burn-in), and Gelman-Rubin R − 1 values below 0.01 for all cosmological parameters.

3.3. Handling of implicit priors

Before diving into our results, it is essential to acknowledge that all reionisation models inherently carry a prior in the space of reionisation histories, x_e(z). Indeed, while we use uniform priors for all parameters of each reionisation model, this choice does not guarantee a uniform distribution of x_e values at any given z, and even less so across an interval of z values. Consequently, this also imposes an implicit prior on the derived value of the optical depth, τ, which can significantly influence its inferred posterior distribution and other reionisation-related quantities. As this τ prior is unique to each reionisation model, not accounting for it would bias any comparison between models since we would not be able to disentangle the effect of the prior versus what model the data actually prefer. We addressed this potential issue via two different approaches, described in the following sections.

3.3.1. Prior flattening

One method consists in precisely characterising this implicit prior on τ in order to cancel its influence on the final posterior distribution. To do so, we consider the methodology proposed by Millea & Bouchet (2018), which involves two key steps. First, for each reionisation model considered we sampled the prior distribution, 𝒫(θ), of its parameters, θ, and computed the corresponding marginalised distribution of derived τ values, 𝒫(τ(θ)). This step allowed us to quantify any intrinsic bias on τ imposed by the chosen reionisation history priors.

Next, we corrected for this prior-induced bias by importance sampling the points in our MCMC chains. More specifically, each MCMC sample was weighted by the inverse of the previously obtained distribution, 𝒫(τ(θ)). This procedure ensures that our posterior distributions of τ are not artificially skewed by the initial choice of priors in the reionisation models. Written explicitly, we can create a posterior where we have assumed a flat τ prior by computing

$$\begin{array}{c}\hfill {\mathcal{P}}^{\mathrm{flat}-\tau }(\theta |d)=\frac{\mathcal{P}(\theta |d)}{\mathcal{P}(\tau (\theta ))}\propto \mathcal{L}(\theta )\underset{{\mathcal{P}}^{\mathrm{flat}-\tau }(\theta )}{\underbrace{\left[\frac{\mathcal{P}(\theta )}{\mathcal{P}(\tau (\theta ))}\right]}},\end{array}$$(5)

where the posterior given data d, 𝒫(θ | d), is the original posterior that did not assume a flat τ prior; 𝒫^flat-τ(θ | d) is the “τ-flattened” posterior, and the likelihood, ℒ(θ), is the probability of the data given parameters, 𝒫(d | θ). Figure 3 illustrates the “flattening" process for a reionisation model whose τ posterior is noticeably affected, namely our 3-steps model.

Fig. 3. Illustration of the flattening process of the τ prior distribution, applied to the 3-steps reionisation model. The original posterior on τ is shown in blue, while the posterior after flattening and the τ prior are shown, respectively, in orange and green.

We note that there are an infinite number of ways to adjust the posterior distribution to achieve a flat prior on τ. However, the aforementioned method ensures that the minimum amount of additional information is introduced into the analysis, or equivalently that the entropy of the prior-adjusted τ posterior ismaximised (see Appendix A of Millea & Bouchet 2018 for details). We further note that this method can be extended to any derived parameter of interest.

3.3.2. Profile likelihoods

While the previously described method effectively removes the effects of prior assumptions on the final τ constraints, it does not address the further possibility of volume effects in the parameter space. Indeed, whole swaths of the oddly shaped parts of the parameter space can have similarly acceptable likelihood – and τ values, which can then affect the posterior distribution of τ after marginalisation. This can be particularly pronounced for models with a large number of degrees of freedom. While such an effect is inherent to – and expected from – a Bayesian formalism, one may want to infer constraints without being sensitive to either the choice of priors or volume effects in the parameter space. To mitigate those effects, we consider an alternative approach based on profile likelihoods.

Using τ as an example, this method consists in successively fixing τ to a grid of predefined values and then minimising the chi-square for each distinct value of τ, varying all other parameters in the analysis (both cosmological and nuisance-related). Selecting the resulting best-fit value and applying a standard chi-square threshold on the resulting likelihood profile then provides, respectively, a robust estimate of τ and its uncertainty, without the influence of any explicit or implicit prior (see e.g. Planck Collaboration Int. XVI 2014, for an application in the context of the cosmological analysis of the Planck CMB data). We apply this method to the reionisation models considered in this work and compare the results to those obtained with the prior-flattening technique.

As a technical detail, we note that to perform the minimisations required by the profile likelihood, we had to reparameterise each of our models to include τ explicitly as an input parameter. For most models, this was simply achieved by introducing a parameter that can shift all redshift-related parameters at once, thus allowing for a direct control of τ. For the pca-5 and 10-bins models, for the vector of model parameters, we switched from cartesian to spherical coordinates, and use τ to scale the norm of the aforementioned vector. Finally, we also had to increase the precision settings of the Boltzmann solvers to facilitate the minimisation procedure. Specifically, for CLASS we identified the most relevant parameters as reionization_sampling (set to 0.01) and the perturbations_sampling_stepsize (set to 0.05), while for CAMB we set the AccuracyBoost parameter to 3.

4. Results on the reionisation optical depth

In this section, we present our measurements of the reionisation optical depth, τ, obtained as a derived parameter from our various reionisation models, as constrained by the Planck PR4 data. In our work we define τ as the integral defined in Eq. (1) with bounds set from z = 0 to z = 30⁵. We emphasize that all τ measurements reported below include an identical contribution (about 0.03) from z < 5 resulting from the fixed timing and duration of the second helium reionisation.

4.1. Posterior distributions

For each of the 12 models considered in this work, we derived the posterior distribution of the reionisation optical depth from our corresponding MCMC chains. We summarise our results in Fig. 4, showing the obtained means and maxima of the posteriors as well as the 68% credible intervals. The left panel of the figure displays the raw posteriors, while the right panel shows the posteriors after applying the prior-flattening correction described in Sect. 3.3.

Fig. 4. Constraints on the reionisation optical depth from Planck PR4 for the reionisation models considered. The plotted segments represent the 68% credible intervals, where the crosses and circles mark, respectively, the maximum and mean of each posterior. The left panel shows the raw posterior results, while the right panel includes the prior-flattening correction described in Sect. 3.3.

We first notice that the measurements of τ are all consistent between the different models. The raw posteriors of τ (left panel of Fig. 4) exhibit some relatively small dispersion with the exception of the 10-bins model that has a higher value compared to the others by approximately 1σ. This aligns with expectations regarding the implicit prior effect, which is anticipated to be more pronounced in higher-dimensional parameter spaces, especially in cases such as ours where the data are not particularly constraining.

This is further confirmed by the results obtained from the posteriors after applying the prior-flattening correction (right panel of Fig. 4), which show better agreement between the different models. This alignment indicates that the correction of the prior effect effectively mitigates the bias induced by model-specific priors, providing a more consistent estimate of τ.

This consistency not only reinforces the credibility of our τ estimates but also reinforces the reliability of other cosmological parameters. Indeed, our findings indicate that the modelling of the reionisation history has no impact on the constraints for the other ΛCDM parameters after applying the τ prior-flattening procedure (see Appendix C).

We also note that the width of the posteriors is quite stable across models, despite significant differences in numbers of degrees of freedom (dof). The Planck PR4 data thus appear to be able to constrain τ with a similar precision, regardless of the complexity of the reionisation model. We may hypothesise about the origin of this observation: as briefly mentioned in Sect. 2, the constraining power of the CMB on reionisation is two-fold: scalar spectra are damped by e^−2τ and the low-ℓ polarisation signal scales as τ² (at first order). Together, these two effects represent the bulk of the constraint on τ from the CMB, and are likely to dominate the total variance budget. However, the large-scale polarisation signal also exhibits additional multipole-dependent features depending on the precise details of the reionisation history, x_e(z), which translates into an additional constraint on τ but is likely to be of secondary order. This could explain why the width of the posteriors for τ is so stable and is not significantly affected by the details and complexity of the chosen reionisation model.

The only notable exception in terms of consistency in the constraints comes from the 10-bins model that, despite the prior-flattening correction, still exhibits a different, slightly smaller (by ∼20%) τ uncertainty compared to the other models. This may seem paradoxical, as the 10-bins model has the largest number of degrees of freedom. However, even though the prior-flattening procedure removes any volume effect affecting τ, it does not remove potential volume effects affecting the likelihood itself. If present, the profile likelihood method (see next section) should be able to identify such remaining effects.

We summarise the quantitative constraints on τ in all the cases considered in Table 2 for each model. Ultimately, from the posterior distributions obtained with Planck CMB data, we measure the following averaged reionisation optical depth:

$$\begin{array}{c}\hfill \tau =0.0576\pm 0.0060\pm 0.0006,\end{array}$$(6)

where the first uncertainty accounts for statistical errors (average of the 68% τ uncertainties across models) while the second is given by the dispersion over the reionisation history models (standard deviation of the mean τ value across models).

4.2. Profile likelihoods

We also computed the profile likelihoods for τ for each reionisation model considered in this work. The results are shown in Fig. 5 and reported in Table 2. The minimisation procedure in ECLAIR enables the reconstruction of the full shape of the profile likelihood, from which we extract the global minimum and define the confidence interval by applying a Δχ² = 1 threshold. Most profiles are observed to be roughly parabolic in shape. The simpler models in particular form a remarkably consistent group with nearly identical profiles. The more complex models allow marginally better chi-square values (differing by at most six units)⁶, while some exhibit deviations from quadratic behaviour in their shape. Indeed, complex models (including tanhdz, pca-5, 10-bins and 3-steps) introduce additional degrees of freedom, subtle dependencies on the parameter, non-linear parameter effects and potential multimodality, all of which can lead to more intricate profile likelihoods.

Fig. 5. Left panel: profile likelihoods of the reionisation optical depth, τ, for the reionisation models considered, where crosses and small vertical bars, respectively, mark the maximum-likelihood values and the Δχ2 < 1 intervals. Right panel: summary of the maximum likelihood values and Δχ2 < 1 intervals derived from the profiles, plotted as segments.

The τ constraints from the profile likelihoods display the same consistency across different reionisation models as observed in the MCMC-derived posteriors. The profile likelihoods also exhibit a stable width across the various models, confirming the robustness of our results. Interestingly, the profile likelihoods make the results on τ even more coherent, including for the 10-bins model that was preferring slightly lower values of τ and tighter errors as noted in the previous section, thus confirming our hypothesis regarding the remaining volume effects.

Following the same prescription as for Eq. (6), we find as constraints on τ from profile likelihoods

$$\begin{array}{c}\hfill \tau =0.0581\pm 0.0061\pm 0.0005.\end{array}$$(7)

As was previously discussed, profile likelihoods are not affected by volume effects, thus any visible discrepancy between them and MCMC results are likely a consequence of these influences. The most notable impact is observed in the 10-bins model (see Table 2). The average value of τ from profile likelihoods is slightly higher (albeit by only ∼0.1σ), while the uncertainties become noticeably more consistent.

5. Results on the history of reionisation

In this section, we examine in more detail the constraints on the evolution of the ionised fraction of matter itself as a function of redshift as well as associated derived quantities.

5.1. Redshift evolution of the ionisation fraction

We explore the reconstructed reionisation histories for all the models considered in Sect. 2. In each case, we plot x_e(z), providing a visual representation of the constraints on the evolution of the ionisation fraction over cosmic time. Each reionisation model, from the classic hyperbolic tangent to the flexknot approach, yields distinct x_e(z) curves that reflect their underlying assumptions and parametrisations. These reconstructions allow us to assess the flexibility and accuracy of each model in fitting the Planck PR4 data, thereby offering insights into the nature and duration of the reionisation process.

We summarise the reconstruction results in Fig. 6, which displays in blue the MCMC distributions of histories of the ionisation fraction as a function of redshift, x_e(z), for all models considered in this work. The blue curves envelop the 68% credible intervals at each z. Given that the sampling of the reionisation models is performed using flat priors on the parameters, this consequently induces non-flat priors on x_e(z). As a result, their posterior distributions may be significantly influenced by prior effects, which can vary substantially between models. Similarly to the previous discussion on the reionisation optical depth (Sect. 4), we also derived the best-fit history for each model, shown as a solid red line. The two dashed red lines mark the boundaries of the region containing models that fall within Δχ² = 1 of the best-fit value, illustrating the confidence interval associated with the best-fit. We note that, as a consequence of the prior effects, combined with complex likelihood shapes associated with the most sophisticated models, the best-fit x_e(z) can fall outside of the 68% credible intervals.

Fig. 6. Ionisation history as a function of the redshift, xe(z), for the reionisation models considered. A few hundred samples from the MCMC distributions are shown in blue, while dashed black lines represent the 68th percentile of the distribution of xe(z) for each redshift z. The solid and dashed red lines, respectively, correspond to the best-fit solutions and the frequentist confidence intervals determined by Δχ2 < 1.

We observe that the simpler models, with less than 3 degrees of freedom (d.o.f.), all exhibit a similar behaviour. These models include tanh, 2-pow, 2-knots, 2-flex, 3-pow, and to a lesser extent tanhdz. They are largely consistent with an instantaneous reionisation scenario, resembling the classic hyperbolic tangent model. This shows that when the model flexibility is limited, both the MCMC exploration and best-fit solutions tend to favour a rapid reionisation history.

In contrast, the intermediate complexity models, with up to 5 d.o.f., show a more nuanced picture. These models include 2-steps, 3-knots, and 3-flex. They consistently exhibit a tail towards higher redshifts, indicating the possibility of a moderate contribution to reionisation at early times (z > 10). This tail is coherent across different models, suggesting that the data can support the presence of a limited ionisation fraction at higher redshifts.

Finally, the highest complexity models (such as 3-steps, pca-5, and 10-bins), offer an even more detailed insight. These models are consistent with the intermediate-complexity models, but we note that their best-fit consistently suggests a non-zero ionisation fraction around redshift z ∼ 15. This peak, while moderate in significance, suggests that reionisation history may have been more complex than expected. We note that, by construction, the pca-5 model allows the ionisation fraction to reach values beyond the physical range [0, 1 + f_He]. This occurs because the completeness of the orthogonal basis is broken when limiting the number of modes, even more so in our case where we only use five (see Sect. 2.6).

As outlined in Sect. 4, the constraints on the ionisation fraction derived from the CMB power spectra are primarily driven by the optical depth. In practice, the reconstructed ionisation histories all conform to the same τ, within the flexibility permitted by each model. More complex models can extract additional information from the detailed shape of the C_ℓ^EE spectrum, providing a further constraining power – though its statistical significance remains limited. We present in Appendix B the E-mode CMB polarisation angular power spectra (C_ℓ^EE) corresponding to the best-fit of each reionisation history model. Notably, models that exhibit a non-zero ionisation fraction at z ∼ 15 show an increase in power at scales around ℓ = 6 − 25 (compared to the more instantaneous models), which corresponds to the region just after the reionisation bump. While this additional power improves the fit to the data, its statistical significance is not high enough to determine whether it originates from a feature of thereionisation history. Indeed, in this multipole range, the signal-to-noise ratio of Planck is below one, and the data may be affected by systematic residuals (arising from foreground or instrumental contamination).

5.2. Constraints on the end of reionisation

In this section, we explore to what extent the CMB could place constraints on the end of reionisation. To investigate this, we extend the 10-bins model by introducing two additional bins at very low redshifts: one spanning z = 0 to z = 2.5, and the other from z = 2.5 to z = 5. Unlike previous models, we do not impose a second reionisation of helium but instead allow the first two bins to vary freely between 0 and 1 + 2f_He.

The resulting ionisation fraction as a function of redshift is shown in Fig. 7. The recovered reionisation history closely resembles that of the 10-bins model, with a global transition occurring between z = 7.5 and z = 5. As in the 10-bins case, the best-fit solution for the 12-bins model also shows a non-zero ionisation fraction around z = 15. Notably, in contrast to the higher redshift bins, we find a preference for a non-zero ionisation fraction in the first two bins, with both a peak of the x_e posterior distribution and a best-fit x_e value close to 1 + 2f_He.

Fig. 7. Ionisation history as a function of the redshift, xe(z), for the 12-bins model. The dashed black lines represent the 68th percentile of the distribution of xe(z) for each redshift z. The red lines correspond to the best-fit solution with the corresponding 1-σ confidence interval given by the dashed red lines.

The optical depth values recovered for this 12-bins model are τ = 0.0548 ± 0.0069 (posterior with flat-τ prior) and $\tau =0.{0575}_{-0.0072}^{+0.0076}$ (profile likelihood). In this case, the mean posterior distribution of τ is slightly lower than in the other models, as this model allows for incomplete reionisation at z = 0, a scenario that was excluded in previous models except for the tanhdz one (which also exhibited such slight shift towards lower τ). However, the best-fit value of τ remains fully consistent with the one obtained from the 10-bins model.

5.3. Constraints on early sources of reionisation

In addition to constraining the ionisation fraction, x_e(z), we also derived estimates for redshift-integrated quantities such as the cumulative reionisation optical depth from z to z_max, τ_[z, zmax] (see Appendix D). As noted by Heinrich et al. (2017), this quantity is more directly connected to the polarisation observables than the instantaneous ionisation fraction, x_e(z). In our case, the shape of τ_[z, zmax] below z = 5 is fixed by construction. Similar to x_e(z), our results show that the contribution to optical depth at z > 10 are highly model-dependent.

Specifically, we consider the optical depth integrated over the redshift range z > 15. Indeed, this quantity is of particular interest for probing early reionisation signatures, as it highlights the contribution to τ from early ionising sources, prior to redshift 15. Figure 8 displays the Planck PR4 constraints on τ_{> 15} for the various models discussed in this work. The left panel presents the posterior distributions, while the right panel depicts the profile likelihoods estimated for τ_{> 15}. As with τ, the posterior distribution for τ_{> 15} is similarly influenced by prior choices. We switch to a flat prior on τ_{> 15} after applying the same type of correction as for τ (see Sect. 3.3).

Fig. 8. Constraints on τ> 15 for the reionisation models considered. Left panel: full posterior distributions, with a flattening procedure applied to the τ> 15 prior. Right panel: profile likelihoods for τ> 15.

The results for the simplest models (with less than 3 degrees of freedom) yield almost identical and tightly constrained values of τ_{> 15}, centred around ∼4.75 × 10⁻⁵, consistent with the minimal residual ionisation remaining from the recombination era. This outcome is within expectations, as these models are designed to capture only the primary transition from a neutral to an ionised intergalactic medium, lacking the flexibility to model an additional, independent contribution to ionisation at earlier times. Indeed, for such models, higher values of τ_{> 15} can only be achieved for sets of parameters that lead to a correspondingly larger total optical depths, which are disfavoured by the data.

The models that allow for a more complex transition with an early contribution at high redshift are also compatible with zero early contribution to reionisation. Only the posteriors for the most complex models, namely the 10-bins and pca-5, exhibit a preference for a non-zero value even after prior flattening. However, the profile likelihoods on τ_{> 15} for these complex models are less decisive and show a large region of equivalent statistical significance between τ_{> 15} = 0 and 0.02.

Based on the profile likelihoods, we derive 95% confidence level upper limits for τ_{> 15} by applying a purely data-based threshold (independent of priors or parameterisations) of Δχ² < 3.841, which corresponds to the critical value of the chi-square distribution for 95%, resulting in the following upper-bounds:

$$\begin{array}{cc}\hfill {\tau }_{>15}<0.028& \text{(10-bins)},\hfill \end{array}$$(8)

$$\begin{array}{cc}\hfill {\tau }_{>15}<0.024& \text{(3-steps)},\hfill \end{array}$$(9)

$$\begin{array}{cc}\hfill {\tau }_{>15}<0.024& \text{(pca-5)},\hfill \end{array}$$(10)

$$\begin{array}{cc}\hfill {\tau }_{>15}<0.018& \text{(2-steps)},\hfill \end{array}$$(11)

$$\begin{array}{cc}\hfill {\tau }_{>15}<0.011& \text{(3-knots)},\hfill \end{array}$$(12)

while the upper-bound is lower than 0.002 for the other models, including tanh, tanhdz, 2-pow, 3-pow, 2-knots, 2-flex, 3-flex. This highlights the challenges in measuring potential contributions from early sources of reionisation with CMB data. Our results show that the constraints on τ_{> 15} are strongly model-dependent. Our upper-bounds are compatible with the results derived from previous Planck data sets with specific models (cf. τ_{> 15} < 0.020 with pca-5, Heinrich & Hu 2021; τ_{> 15} < 0.018 with FlexKnot, Planck Collaboration VI 2021).

6. Conclusions

In this work, we have reconstructed constraints on the history of cosmic reionisation using the measurements of CMB anisotropies, systematically exploring a wide range of reionisation models, from simple parametric forms to more flexible, non-parametric approaches. By employing both Bayesian and frequentist methods, we have carefully examined the impact of model-dependent priors and volume effects on the inferred constraints.

Our work highlights the importance of accounting for implicit priors when interpreting CMB-based constraints on reionisation. Different modelling approaches can introduce substantial biases, particularly for models with higher degrees of freedom, where prior effects become more pronounced. By correcting for these biases using prior-flattening techniques and complementing MCMC posteriors with profile likelihoods, we have ensured a more robust and unbiased interpretation of the results.

Our analysis confirms that Planck PR4 data provide robust constraints on the integrated optical depth, τ, reaching an accuracy of approximatively 10%. We find consistent τ estimates across the various reionisation models considered. Assuming a flat prior on τ, we obtain τ = 0.0576 ± 0.0060, to which we add an additional uncertainty of around 0.0006 to account for reionisation history modelling. In addition to prior effects, we clearly identify the presence of volume effects in the distribution of reionisation parameters. While their impact remains limited with Planck data, such effects will require careful consideration in future CMB surveys with higher precision, where they could play a more significant role in shaping parameter constraints.

Despite consistent estimates of the optical depth, the reconstruction of the ionisation history, x_e(z), remains strongly model-dependent, as the CMB primarily constrains the integrated optical depth rather than the detailed redshift evolution of reionisation. Simpler models, which impose tighter assumptions, typically favour a rapid and relatively late (z < 10) reionisation process, closely matching the standard hyperbolic tangent scenario. In contrast, more flexible models accommodate a wider variety of histories, including some early ionisation at z ∼ 15. However, this feature does not appear in the Bayesian posterior distributions and its significance remains below the 1-σ level. As such, it cannot be determined from current data whether this early ionisation signal is physically motivated or simply a result of model flexibility.

Building upon the model-dependent nature of the reconstructed x_e(z), we further explored the potential for early ionisation by examining the parameter τ_{> 15}, which isolates the contribution to the optical depth from redshifts z > 15. Among all models considered, only the more flexible 10-bins model shows a slight preference for non-zero values in the posterior distribution. However, the corresponding profile likelihood remains broad, reflecting significant degeneracies and reaffirming that current CMB data alone lack the constraining power to robustly determine the presence of early ionising sources. From the profile likelihoods, we derived upper-bounds on τ_{> 15}, which are highly model-dependent, ranging from negligible up to τ_{> 15} < 0.028 depending on the model considered.

Overall, our results show that the low value of τ measured from CMB data could be naturally reconciled with a high-redshift onset of reionisation suggested by recent JWST observations, by adopting more complex reionisation histories that allow for a slow and extended initial phase of ionisation.

While Planck PR4 data represent the most precise CMB measurements to date, further progress in understanding the EoR will require complementary observations. Next-generation CMB experiments, such as LiteBIRD (LiteBIRD Collaboration 2023) or CLASS (Watts et al. 2018), will significantly improve the precision of optical depth measurements through their enhanced sensitivity to large-scale polarisation anisotropies. This will help to break parameter degeneracies – such as with the sum of neutrino masses – and provide a more accurate anchor for the reionisation timeline. In particular, more accurate measurements of the large-scale CMB anisotropies in the ℓ = 10 − 25 range will help constraining possible ionisation onset at high redshift (Watts et al. 2020). However, our study shows that while current CMB data can robustly constrain the overall integrated optical depth, achieving a complete understanding of the redshift-dependent evolution of reionisation remains very challenging. This difficulty arises because the CMB spectrum is only weakly sensitive to the detailed shape of the reionisation history, leading to strong parameter degeneracies when fitting the latter to CMB data alone. Future observations of the 21 cm power spectrum, such as HERA or SKA, will be crucial in directly mapping the spatial distribution of neutral hydrogen at high redshifts (Morales & Wyithe 2010; Pritchard & Loeb 2012; Mellema et al. 2013; DeBoer et al. 2017). The combination of upcoming CMB with 21 cm and UV luminosity observations will enable tighter constraints on the timing, duration, and morphology of reionisation, helping to identify the dominant ionising sources that shaped this pivotal phase of cosmic history (Qin et al. 2020; Paoletti et al. 2021).

Data availability

All data products generated in the course of this analysis are publicly available on Zenodo at 10.5281/zenodo.15706271

Acknowledgments

The authors acknowledge the support of the French Agence Nationale de la Recherche (ANR), under grant ANR-22-CE31-0010 (project BATMAN). This work was supported by the French national space agency (CNES). This work was performed using computational resources from the “Mésocentre” computing centre of Université Paris-Saclay, CentraleSupélec and École Normale Supérieure Paris-Saclay supported by CNRS and Région Île-de-France(https://mesocentre.universite-paris-saclay.fr).

References

Appendix A: Effects of adding lensing likelihood

In Fig. A.1, we show the impact of the inclusion of the latest lensing likelihoods from Planck PR4 (Carron et al. 2022) and ACT DR6 (Madhavacheril et al. 2024; Qu et al. 2024) on constraints for the tanh reionisation model. As expected, adding lensing tightens the constraints on the amplitude of the primordial power spectrum, A_s, reflecting its sensitivity to the matter distribution at intermediate redshifts (z∼ 1 – 5). However, it has negligible effect on the marginalised constraints on the reionisation optical depth, τ, with nearly identical 68% confidence intervals compared to analyses without lensing. While the lensing power spectrum, C_ℓ^ϕϕ, is almost independent of τ, the lensing likelihood depends on the lensed CMB spectra through the bias terms (N⁽⁰⁾, N⁽¹⁾), mean-field, and estimator normalisation. In practice, with the current noise level, the τ–A_s degeneracy remains mostly aligned with that from primary CMB anisotropies, so lensing does not significantly help in breaking it.

Fig. A.1. Effects of adding the Planck CMB lensing likelihood to the tanh model constraints on τ and As. The contours show the 68% and 95% credible regions for the parameters.

Appendix B: Angular power spectra

In this appendix, we discuss the impact of the reionisation history on the CMB E-mode polarisation angular power spectra. Figure B.1 shows D_ℓ = ℓ(ℓ + 1)C_ℓ^EE/2π for the best-fit corresponding to each reionisation model described in the main text. We observe that models with a greater number of degrees of freedom, which allow for a more extended reionisation history, consistently exhibit increased power at multipole scales between ℓ = 5 and ℓ = 25. On the contrary, in these models, the reionisation bump appears lower in amplitude, maintaining a consistent total reionisation optical depth across all models.

We note on the other hand that the TT and TE power spectra do not change significantly across the various best-fit reionisation history, as they are primarily sensitive to the amplitude of the primordial power spectrum and the reionisation optical depth.

Fig. B.1. Angular power spectra of the CMB E-mode polarisation for the best-fit set of parameters for each reionisation model considered. The grey band indicates the Planck PR4 data full error bars (including sample variance, instrumental noise and systematics), centred on the tanh best-fit.

Appendix C: Cosmological parameters

Figure C.1 presents the posterior distributions for the ΛCDM parameters obtained when fitting the different reionisation models explored in this study, after the flattening procedure has been applied to the τ prior. These parameters include the physical baryon density, Ω_bh², the physical density of cold dark matter, Ω_cdmh², the amplitude of primordial scalar fluctuation, A_s, the scalar index, n_s, and the Hubble constant, H₀. We also include the reionisation optical depth, τ, derived from the reconstructed reionisation history. Our results indicate that the inferred constraints on these parameters are robust and largely unaffected by the specific reionisation model assumed. The only exception is n_s, which shows a slight shift towards higher values in the case of the pca-5 model. This may be due to minor differences between CLASS and CAMB, as the latter is used to compute the spectra for the pca-5 model. However, we were able to reproduce identical parameter distributions – including n_s – for the tanh model using both codes, suggesting that the shift is specific to the pca-5 implementation rather than the choice of Boltzmann solver.

Fig. C.1. Posterior distributions for the ΛCDM cosmological parameters, including the reionisation optical depth, τ, as a derived parameter, obtained from Planck PR4 data for the reionisation models considered.

Appendix D: Cumulative optical depth

The cumulative optical depth quantifies the integral of the electron density weighted by the Thomson cross section along the line of sight, and is calculated between redshifts z and z_max as:

$$\begin{array}{c}\hfill \tau (z,z_{\max })={\displaystyle {\int }_z^{z_{\max }}}{\sigma }_T{n}_{\mathrm{e}}(z)\frac{\mathit{dr}}{\mathit{dz}}dz.\end{array}$$(D.1)

By definition, the total optical depth is given by τ = τ(0, z_max). Figure D.1 displays the reconstructed τ(z, z_max) for each reionisation model considered in this work. Because we fixed the ionisation fraction, x_e(z), up to z = 5, the shape of τ(z, z_max) below z = 5 is fixed by construction. Our results indicate that the contribution to optical depth at z > 10 are highly model-dependent. Among models with enough flexibility to allow some contribution to τ at high redshifts, the 2-steps, 3-steps, 3-knots, 10-bins and pca-5 show non-zero contribution at z > 15.

Fig. D.1. Cumulative optical depth (68% and 95% C.L. contours) for the reionisation models considered. The dashed black lines represent the 68th percentile of the distribution of τ(z) for each redshift z. The red lines correspond to the best-fit solution (solid) with the corresponding 1-σ confidence interval (dashed).

All Tables

All Figures

Fig. 1. Mean values (filled circles) and 68% credible intervals (horizontal bars) obtained from posterior distributions of τ, derived from CMB power spectra measurements including: WMAP 2013 (Hinshaw et al. 2013), Planck PR1 + WMAP 2014 (Planck Collaboration XV 2014), Planck PR2 + WMAP 2016 (Planck Collaboration XIII 2016), Planck-LFI 2016 (Planck Collaboration XIII 2016), Planck-HFI 2016 (Planck Collaboration Int. XLVII 2016), WMAP+LFI 2018 (Weiland et al. 2018), Planck PR3 2020 (Planck Collaboration VI 2020), Planck 2020 (Pagano et al. 2020), WMAP+LFI 2022 (Paradiso et al. 2023), Planck PR4 2022 (Planck Collaboration Int. LVII 2020). The shaded grey region highlights the cosmic variance limit on the measurement of τ, centred on the value from Planck PR4 2024 (Tristram et al. 2024).

In the text

	Fig. 2. Graphical summary of the reionisation models considered in the present study, illustrating with coloured arrows and lines their respective degrees of freedom in modelling the history of the ionising fraction. The models are grouped by their general approach: hyperbolic tangents, power-laws, step functions, interpolations, bins, and PCA. The models are described in detail in Sect. 2.
In the text

	Fig. 3. Illustration of the flattening process of the τ prior distribution, applied to the 3-steps reionisation model. The original posterior on τ is shown in blue, while the posterior after flattening and the τ prior are shown, respectively, in orange and green.
In the text

	Fig. 4. Constraints on the reionisation optical depth from Planck PR4 for the reionisation models considered. The plotted segments represent the 68% credible intervals, where the crosses and circles mark, respectively, the maximum and mean of each posterior. The left panel shows the raw posterior results, while the right panel includes the prior-flattening correction described in Sect. 3.3.
In the text

	Fig. 5. Left panel: profile likelihoods of the reionisation optical depth, τ, for the reionisation models considered, where crosses and small vertical bars, respectively, mark the maximum-likelihood values and the Δχ2 < 1 intervals. Right panel: summary of the maximum likelihood values and Δχ2 < 1 intervals derived from the profiles, plotted as segments.
In the text

Fig. 6. Ionisation history as a function of the redshift, xe(z), for the reionisation models considered. A few hundred samples from the MCMC distributions are shown in blue, while dashed black lines represent the 68th percentile of the distribution of xe(z) for each redshift z. The solid and dashed red lines, respectively, correspond to the best-fit solutions and the frequentist confidence intervals determined by Δχ2 < 1.

In the text

	Fig. 7. Ionisation history as a function of the redshift, xe(z), for the 12-bins model. The dashed black lines represent the 68th percentile of the distribution of xe(z) for each redshift z. The red lines correspond to the best-fit solution with the corresponding 1-σ confidence interval given by the dashed red lines.
In the text

	Fig. 8. Constraints on τ> 15 for the reionisation models considered. Left panel: full posterior distributions, with a flattening procedure applied to the τ> 15 prior. Right panel: profile likelihoods for τ> 15.
In the text

	Fig. A.1. Effects of adding the Planck CMB lensing likelihood to the tanh model constraints on τ and As. The contours show the 68% and 95% credible regions for the parameters.
In the text

	Fig. B.1. Angular power spectra of the CMB E-mode polarisation for the best-fit set of parameters for each reionisation model considered. The grey band indicates the Planck PR4 data full error bars (including sample variance, instrumental noise and systematics), centred on the tanh best-fit.
In the text

	Fig. C.1. Posterior distributions for the ΛCDM cosmological parameters, including the reionisation optical depth, τ, as a derived parameter, obtained from Planck PR4 data for the reionisation models considered.
In the text

	Fig. D.1. Cumulative optical depth (68% and 95% C.L. contours) for the reionisation models considered. The dashed black lines represent the 68th percentile of the distribution of τ(z) for each redshift z. The red lines correspond to the best-fit solution (solid) with the corresponding 1-σ confidence interval (dashed).
In the text