# Constraining cosmic reionization by combining the kinetic Sunyaev-Zel'dovich and the 21 cm power spectra

Ivelin Georgiev, Adélie Gorce, Garrelt Mellema

► To cite this version:

Ivelin Georgiev, Adélie Gorce, Garrelt Mellema. Constraining cosmic reionization by combining the kinetic Sunyaev-Zel'dovich and the 21 cm power spectra. Monthly Notices of the Royal Astronomical Society, 2024, 528 (4), pp.7218-7235. 10.1093/mnras/stae506. hal-04356160

HAL Id: hal-04356160

https://hal.science/hal-04356160

Submitted on 20 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Constraining cosmic reionization by combining the kinetic Sunyaev–Zel’dovich and the 21 cm power spectra

Ivelin Georgiev<fnref n="1" />,<fnref n="1" />★ Adélie Gorce<fnref n="2" />,<fnref n="3" /> and Garrelt Mellema<fnref n="1" />

Accepted 2024 February 14. Received 2024 February 12; in original form 2023 December 7

## ABSTRACT

During the Epoch of Reionization (EoR), the ultraviolet radiation from the first stars and galaxies ionized the neutral hydrogen of the intergalactic medium, which can emit radiation through its 21 cm hyperfine transition. Measuring the 21 cm power spectrum is a key science goal for the future Square Kilometre Array (SKA); however, observing and interpreting it is a challenging task. Another high-potential probe of the EoR is the patchy kinetic Sunyaev–Zel’dovich (pkSZ) effect, observed as a foreground to the cosmic microwave background temperature anisotropies on small scales. Despite recent promising measurements, placing constraints on reionization from pkSZ observations is a non-trivial task, subject to strong model dependence. We propose to alleviate the difficulties in observing and interpreting the 21 cm and pkSZ power spectra by combining them. With a simple yet effective parametric model that establishes a formal connection between them, we can jointly fit mock 21 cm and pkSZ data points. We confirm that these observables provide complementary information on reionization, leading to significantly improved constraints when combined. We demonstrate that with as few as two measurements of the 21 cm power spectrum with 100 h of observations with the SKA, as well as a single $\ell = 3000$ pkSZ data point, we can reconstruct the reionization history of the universe and its morphology. We find that the reionization history (morphology) is better constrained with two 21 cm measurements at different redshifts (scales). Therefore, a combined analysis of the two probes will give access to tighter constraints on cosmic reionization even in the early stages of 21 cm detections.

**Key words:** cosmological parameters – dark ages, reionization, first stars – observations.

## 1 INTRODUCTION

The Epoch of Reionization (EoR) is a large-scale phase transition during which the Universe transitioned from a cold and neutral to a hot and ionized state (for a review, see e.g. Wise 2019). During the EoR, the first stars and galaxies formed in the densest regions of the Universe due to the accretion of baryonic matter on to dark matter (DM) haloes. The radiation produced by these young stars and galaxies ionized the neutral hydrogen in the intergalactic medium (IGM), forming H II ‘bubbles’ that progressively grew and overlapped as new ionizing sources formed. Because of this, the properties of cosmic reionization contain information on both cosmology and astrophysics.

A powerful probe of the EoR comes from the measurement of the Thomson scattering optical depth $\tau$ from the cosmic microwave background (CMB). As CMB photons emitted during the recombination epoch travel through the IGM, they scatter off the free electrons produced during reionization. The Thomson scattering generates a polarization signal on scales larger than the horizon scale during reionization while it suppresses the temperature anisotropies on angular scales lower than the horizon size during reionization.

Assuming an instantaneous reionization history, measurements of $\tau = 0.051 \pm 0.006$ by the Planck Collaboration I (2020) indicate that the mid-point of EoR lies around a redshift of $z_{re} \sim 8$. On the other hand, observations of the fluctuations of the Ly $\alpha$ optical depth caused by the Gunn–Peterson effect in high-$z$ quasar spectra (Bosman et al. 2022) and the inferred low mean free path of ionizing photons ($\lambda_{MFP}$; Gaikwad et al. 2023) hint that reionization may extend past redshift six and complete by $z_{end} \sim 5.2$.

Another promising probe of the reionization process is the kinetic Sunyaev–Zel’dovich (SZ) effect where rest-frame CMB photons scatter off the free electrons along the line of sight. As the free electrons from the EoR have a non-zero bulk velocity relative to that of the CMB photons, the latter gain or lose energy, producing secondary temperature anisotropies in the observed CMB (Sunyaev & Zeldovich 1980). This effect occurs during and after the EoR so that it can be divided into two stages. The ‘homogeneous’ kSZ effect is related to the ionized IGM of the post-EoR Universe (Shaw, Rudd & Nagai 2012), while the patchy kinetic Sunyaev–Zel’dovich (pkSZ) effect has been shown to depend strongly on the morphology of H II bubbles during reionization (McQuinn et al. 2005; Mesinger, McQuinn & Spergel 2012; Alvarez 2016; Chen et al. 2023). The pkSZ signal is an integrated observable and the amplitude and peak of its angular power spectrum contains information about, e.g. the duration of reionization and the characteristic sizes of H II

© 2024 The Author(s).

Published by Oxford University Press on behalf of Royal Astronomical Society. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.

<footnote n="1">¹The Oskar Klein Centre, Department of Astronomy, Stockholm University, AlbaNova, SE-10691 Stockholm, Sweden</footnote>

<footnote n="2">²Department of Physics and Trottier Space Institute, McGill University, 3600 University Street, Montreal, QC H3A 2T8, Canada</footnote>

<footnote n="3">³Institut d'Astrophysique Spatiale, CNRS, Université Paris-Saclay, F-91405 Orsay, France</footnote>

<footnote>*E-mail: ivelin.georgiev@astro.su.se</footnote>
<figure><img src="image_2.png" /></figure>

# Constraining cosmic reionization by combining the kinetic Sunyaev–Zel’dovich and the 21 cm power spectra

Ivelin Georgiev<fnref n="1" />,<fnref n="1" />★ Adélie Gorce<fnref n="2" />,<fnref n="3" /> and Garrelt Mellema<fnref n="1" />

Accepted 2024 February 14. Received 2024 February 12; in original form 2023 December 7

## ABSTRACT

During the Epoch of Reionization (EoR), the ultraviolet radiation from the first stars and galaxies ionized the neutral hydrogen of the intergalactic medium, which can emit radiation through its 21 cm hyperfine transition. Measuring the 21 cm power spectrum is a key science goal for the future Square Kilometre Array (SKA); however, observing and interpreting it is a challenging task. Another high-potential probe of the EoR is the patchy kinetic Sunyaev–Zel’dovich (pkSZ) effect, observed as a foreground to the cosmic microwave background temperature anisotropies on small scales. Despite recent promising measurements, placing constraints on reionization from pkSZ observations is a non-trivial task, subject to strong model dependence. We propose to alleviate the difficulties in observing and interpreting the 21 cm and pkSZ power spectra by combining them. With a simple yet effective parametric model that establishes a formal connection between them, we can jointly fit mock 21 cm and pkSZ data points. We confirm that these observables provide complementary information on reionization, leading to significantly improved constraints when combined. We demonstrate that with as few as two measurements of the 21 cm power spectrum with 100 h of observations with the SKA, as well as a single $\ell = 3000$ pkSZ data point, we can reconstruct the reionization history of the universe and its morphology. We find that the reionization history (morphology) is better constrained with two 21 cm measurements at different redshifts (scales). Therefore, a combined analysis of the two probes will give access to tighter constraints on cosmic reionization even in the early stages of 21 cm detections.

**Key words:** cosmological parameters – dark ages, reionization, first stars – observations.

## 1 INTRODUCTION

The Epoch of Reionization (EoR) is a large-scale phase transition during which the Universe transitioned from a cold and neutral to a hot and ionized state (for a review, see e.g. Wise 2019). During the EoR, the first stars and galaxies formed in the densest regions of the Universe due to the accretion of baryonic matter on to dark matter (DM) haloes. The radiation produced by these young stars and galaxies ionized the neutral hydrogen in the intergalactic medium (IGM), forming H II ‘bubbles’ that progressively grew and overlapped as new ionizing sources formed. Because of this, the properties of cosmic reionization contain information on both cosmology and astrophysics.

A powerful probe of the EoR comes from the measurement of the Thomson scattering optical depth $\tau$ from the cosmic microwave background (CMB). As CMB photons emitted during the recombination epoch travel through the IGM, they scatter off the free electrons produced during reionization. The Thomson scattering generates a polarization signal on scales larger than the horizon scale during reionization while it suppresses the temperature anisotropies on angular scales lower than the horizon size during reionization.

Assuming an instantaneous reionization history, measurements of $\tau = 0.051 \pm 0.006$ by the Planck Collaboration I (2020) indicate that the mid-point of EoR lies around a redshift of $z_{re} \sim 8$. On the other hand, observations of the fluctuations of the Ly $\alpha$ optical depth caused by the Gunn–Peterson effect in high-$z$ quasar spectra (Bosman et al. 2022) and the inferred low mean free path of ionizing photons ($\lambda_{MFP}$; Gaikwad et al. 2023) hint that reionization may extend past redshift six and complete by $z_{end} \sim 5.2$.

Another promising probe of the reionization process is the kinetic Sunyaev–Zel’dovich (SZ) effect where rest-frame CMB photons scatter off the free electrons along the line of sight. As the free electrons from the EoR have a non-zero bulk velocity relative to that of the CMB photons, the latter gain or lose energy, producing secondary temperature anisotropies in the observed CMB (Sunyaev & Zeldovich 1980). This effect occurs during and after the EoR so that it can be divided into two stages. The ‘homogeneous’ kSZ effect is related to the ionized IGM of the post-EoR Universe (Shaw, Rudd & Nagai 2012), while the patchy kinetic Sunyaev–Zel’dovich (pkSZ) effect has been shown to depend strongly on the morphology of H II bubbles during reionization (McQuinn et al. 2005; Mesinger, McQuinn & Spergel 2012; Alvarez 2016; Chen et al. 2023). The pkSZ signal is an integrated observable and the amplitude and peak of its angular power spectrum contains information about, e.g. the duration of reionization and the characteristic sizes of H II

© 2024 The Author(s).

Published by Oxford University Press on behalf of Royal Astronomical Society. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.

<footnote n="1">¹The Oskar Klein Centre, Department of Astronomy, Stockholm University, AlbaNova, SE-10691 Stockholm, Sweden</footnote>

<footnote n="2">²Department of Physics and Trottier Space Institute, McGill University, 3600 University Street, Montreal, QC H3A 2T8, Canada</footnote>

<footnote n="3">³Institut d'Astrophysique Spatiale, CNRS, Université Paris-Saclay, F-91405 Orsay, France</footnote>

<footnote>*E-mail: ivelin.georgiev@astro.su.se</footnote>
bubbles, respectively (Zahn et al. 2005; Iliev et al. 2007; Gorce et al. 2020). By combining data from the South Pole Telescope<fnref n="1" /> (Ruhl et al. 2004) and the Planck satellite<fnref n="2" /> (Planck Collaboration I 2020), George et al. (2015) and, later, Reichardt et al. (2021) have constrained the amplitude of the pkSZ angular power spectrum to $D_{3000}^{\text{pkSZ}} = 3.0 \pm 1.0 \mu\text{K}^2$ using the post-reionization models described in Shaw, Rudd & Nagai (2012) for the homogeneous part and Battaglia et al. (2013) for the patchy component. They deduce a $2\sigma$ upper limit on duration of reionization from 25 to 75 per cent of $\Delta z < 4.1$. However, these upper limits are loosened when accounting for the angular correlation between the cosmic infrared background and thermal SZ power spectrum (Reichardt et al. 2021). Moreover, the equations relating the pkSZ amplitude to reionization parameters used to derive such constraints are model dependent (Park et al. 2013). Zahn et al. (2012) show the pkSZ amplitude can be suppressed by up to $1.0 \mu\text{K}^2$ due to radiative cooling and depending on the star formation models considered, as they affect the mean gas density within clusters. In addition, in order to properly model the pkSZ power spectrum, the required simulations must be large as well as highly resolved (see e.g. Shaw, Rudd & Nagai 2012). One way to get around this computational challenge is using a parametrized model calibrated on hydrodynamical simulations (see Gorce et al. 2020, for an example).

On the other hand, an extremely promising probe of cosmic reionization comes from the 21 cm signal emitted by neutral hydrogen within the IGM. One prospect in detecting this signal is measuring the spherically averaged power spectrum of its spatial fluctuations. This power spectrum contains information about, e.g. the global neutral fraction of the IGM and the growth of the ionizing bubbles during reionization (Furlanetto, Zaldarriaga & Hernquist 2004; Furlanetto, McQuinn & Hernquist 2006; McQuinn et al. 2006). For example, Georgiev et al. (2022) find a transition scale within the 21 cm power spectrum which can be directly related to the value of the mean-free path of ionizing photons $\lambda_{\text{MFP}}$ through an empirical formula $k_{\text{trans}} \approx 2/\lambda_{\text{MFP}}$. Examples of low-frequency radio interferometers and the large variety of upper limit values on the 21 cm power spectrum have been reported at various redshifts and scales by, e.g. the Low-Frequency Array (LOFAR<fnref n="3" />; Mertens et al. 2020), the Hydrogen Epoch of Reionization Array (HERA<fnref n="4" />; The HERA Collaboration 2022; HERA Collaboration 2023), the Murchison Widefield Array (MWA<fnref n="5" />; Trott et al. 2020; Yoshiura et al. 2021), the Giant Metrewave Radio Telescope (GMRT<fnref n="6" />; Paciga et al. 2013), as well as the forthcoming Square Kilometre Array (SKA<fnref n="7" />; Koopmans et al. 2015). However, no detection has been made and the derived upper limits set only weak constraints on astrophysical quantities (The HERA Collaboration 2022; HERA Collaboration 2023).

Indeed, the 21 cm power spectrum is affected by extra-galactic foregrounds from radio bright sources, radio frequency interference, and ionospheric activity, which complicate the calibration of the antennas of the radio interferometers and, in turn, its measurement. However, the foreground signal is anticipated to mainly affect the lower k region of the spectrum in Fourier space and techniques

of foreground avoidance, suppression, and subtraction have been considered (see e.g. Chapman et al. 2013; Liu, Parsons & Trott 2014; Mertens, Ghosh & Koopmans 2018).

In this work, we leverage the complementarity of 21 cm and kSZ observations in probing cosmic reionization in order to obtain significant constraints before a complete measurement of the 21 cm power spectrum is achieved. In Bégin, Liu & Gorce (2022), the authors already demonstrate this complementarity at the level of the global 21 cm signal, showing that a combined analysis makes it possible to reconstruct a model-independent reionization history, with no assumed parametrization of the redshift-evolution. Here, we push this analysis a step further by including second-order statistics in the assessment and considering the 21 cm power spectrum. We place constraints on cosmic reionization by utilizing the fundamental relationship between the power spectra of the 21 cm and the pkSZ signal from the EoR, which we formalize through a simple yet effective parametric method.

This paper is organized as follows. In Section 2, we describe the physics behind the 21 cm signal from the EoR, the derivation of the pkSZ angular power spectrum, and the formalism used to formally link the 21 cm signal and the pkSZ effect, and, in turn, their power spectra. We also introduce the methodology of the Markov chain Monte Carlo (MCMC) sampling used in our forecast. Section 3 briefly outlines the simulation used to check the accuracy of our forecast. In Section 4, we present our findings for an ideal case, with the detection of both the 21 cm with SKA-Low for two redshifts as well as the pkSZ power spectra. We discuss how the forecast behaves and how constraining it is on the parameters of reionization. In Section 5, we outline the caveats of the forecast and explore special cases. We summarize our conclusion and elaborate on future improvements in Section 6. Unless stated otherwise, we address comoving megaparsec as Mpc.

## 2 METHODS

In this section, we study the connection between the neutral and ionized density components of the IGM during the EoR by investigating the relation of the 21 cm and pkSZ power spectra. We go through the derivation of the 21 cm power spectrum, where we construct our model based on its relation to the electron autocorrelation power spectrum. Using a parametrization of the latter, we present how to analytically reconstruct both the 21 cm and pkSZ signals and generate mock data. Lastly, we describe the methodology of the forecast analysis based on the mock data.

### 2.1 Relation between the 21 cm and pkSZ signals and the electron overdensity field

#### 2.1.1 Derivation of the 21 cm power spectrum

The spin-flip transition between the hyperfine states of a neutral hydrogen H1 atom results in the emission or absorption of photons with a 21 cm wavelength (Field 1957). The abundance of H1 in the IGM makes this radiation an appealing probe of the large-scale structure of the Universe and a scientific goal for radio interferometry telescopes such as the SKA-Low.

The 21 cm signal from the EoR is seen against the Rayleigh-Jeans tail of the CMB. Assuming small optical depth and neglecting redshift space distortions, the brightness temperature of the signal can be written as

$$ \delta T_{21}(\mathbf{r}, z) = T_0(\mathbf{r}, z) x_{\text{H1}}(\mathbf{r}, z) [1 + \delta_b(\mathbf{r}, z)], \quad (1) $$

<footnote n="1">¹https://pole.uchicago.edu</footnote>

<footnote n="2">²https://www.esa.int/planck</footnote>

<footnote n="3">³www.lofar.org</footnote>

<footnote n="4">⁴https://reionization.org</footnote>

<footnote n="5">⁵www.mwatelescope.org</footnote>

<footnote n="6">⁶www.gmrt.ncra.tifr.res.in</footnote>

<footnote n="7">⁷www.skatelescope.org</footnote>
bubbles, respectively (Zahn et al. 2005; Iliev et al. 2007; Gorce et al. 2020). By combining data from the South Pole Telescope<fnref n="1" /> (Ruhl et al. 2004) and the Planck satellite<fnref n="2" /> (Planck Collaboration I 2020), George et al. (2015) and, later, Reichardt et al. (2021) have constrained the amplitude of the pkSZ angular power spectrum to $D_{3000}^{\text{pkSZ}} = 3.0 \pm 1.0 \mu\text{K}^2$ using the post-reionization models described in Shaw, Rudd & Nagai (2012) for the homogeneous part and Battaglia et al. (2013) for the patchy component. They deduce a $2\sigma$ upper limit on duration of reionization from 25 to 75 per cent of $\Delta z < 4.1$. However, these upper limits are loosened when accounting for the angular correlation between the cosmic infrared background and thermal SZ power spectrum (Reichardt et al. 2021). Moreover, the equations relating the pkSZ amplitude to reionization parameters used to derive such constraints are model dependent (Park et al. 2013). Zahn et al. (2012) show the pkSZ amplitude can be suppressed by up to $1.0 \mu\text{K}^2$ due to radiative cooling and depending on the star formation models considered, as they affect the mean gas density within clusters. In addition, in order to properly model the pkSZ power spectrum, the required simulations must be large as well as highly resolved (see e.g. Shaw, Rudd & Nagai 2012). One way to get around this computational challenge is using a parametrized model calibrated on hydrodynamical simulations (see Gorce et al. 2020, for an example).

On the other hand, an extremely promising probe of cosmic reionization comes from the 21 cm signal emitted by neutral hydrogen within the IGM. One prospect in detecting this signal is measuring the spherically averaged power spectrum of its spatial fluctuations. This power spectrum contains information about, e.g. the global neutral fraction of the IGM and the growth of the ionizing bubbles during reionization (Furlanetto, Zaldarriaga & Hernquist 2004; Furlanetto, McQuinn & Hernquist 2006; McQuinn et al. 2006). For example, Georgiev et al. (2022) find a transition scale within the 21 cm power spectrum which can be directly related to the value of the mean-free path of ionizing photons $\lambda_{\text{MFP}}$ through an empirical formula $k_{\text{trans}} \approx 2/\lambda_{\text{MFP}}$. Examples of low-frequency radio interferometers and the large variety of upper limit values on the 21 cm power spectrum have been reported at various redshifts and scales by, e.g. the Low-Frequency Array (LOFAR<fnref n="3" />; Mertens et al. 2020), the Hydrogen Epoch of Reionization Array (HERA<fnref n="4" />; The HERA Collaboration 2022; HERA Collaboration 2023), the Murchison Widefield Array (MWA<fnref n="5" />; Trott et al. 2020; Yoshiura et al. 2021), the Giant Metrewave Radio Telescope (GMRT<fnref n="6" />; Paciga et al. 2013), as well as the forthcoming Square Kilometre Array (SKA<fnref n="7" />; Koopmans et al. 2015). However, no detection has been made and the derived upper limits set only weak constraints on astrophysical quantities (The HERA Collaboration 2022; HERA Collaboration 2023).

Indeed, the 21 cm power spectrum is affected by extra-galactic foregrounds from radio bright sources, radio frequency interference, and ionospheric activity, which complicate the calibration of the antennas of the radio interferometers and, in turn, its measurement. However, the foreground signal is anticipated to mainly affect the lower k region of the spectrum in Fourier space and techniques

of foreground avoidance, suppression, and subtraction have been considered (see e.g. Chapman et al. 2013; Liu, Parsons & Trott 2014; Mertens, Ghosh & Koopmans 2018).

In this work, we leverage the complementarity of 21 cm and kSZ observations in probing cosmic reionization in order to obtain significant constraints before a complete measurement of the 21 cm power spectrum is achieved. In Bégin, Liu & Gorce (2022), the authors already demonstrate this complementarity at the level of the global 21 cm signal, showing that a combined analysis makes it possible to reconstruct a model-independent reionization history, with no assumed parametrization of the redshift-evolution. Here, we push this analysis a step further by including second-order statistics in the assessment and considering the 21 cm power spectrum. We place constraints on cosmic reionization by utilizing the fundamental relationship between the power spectra of the 21 cm and the pkSZ signal from the EoR, which we formalize through a simple yet effective parametric method.

This paper is organized as follows. In Section 2, we describe the physics behind the 21 cm signal from the EoR, the derivation of the pkSZ angular power spectrum, and the formalism used to formally link the 21 cm signal and the pkSZ effect, and, in turn, their power spectra. We also introduce the methodology of the Markov chain Monte Carlo (MCMC) sampling used in our forecast. Section 3 briefly outlines the simulation used to check the accuracy of our forecast. In Section 4, we present our findings for an ideal case, with the detection of both the 21 cm with SKA-Low for two redshifts as well as the pkSZ power spectra. We discuss how the forecast behaves and how constraining it is on the parameters of reionization. In Section 5, we outline the caveats of the forecast and explore special cases. We summarize our conclusion and elaborate on future improvements in Section 6. Unless stated otherwise, we address comoving megaparsec as Mpc.

## 2 METHODS

In this section, we study the connection between the neutral and ionized density components of the IGM during the EoR by investigating the relation of the 21 cm and pkSZ power spectra. We go through the derivation of the 21 cm power spectrum, where we construct our model based on its relation to the electron autocorrelation power spectrum. Using a parametrization of the latter, we present how to analytically reconstruct both the 21 cm and pkSZ signals and generate mock data. Lastly, we describe the methodology of the forecast analysis based on the mock data.

### 2.1 Relation between the 21 cm and pkSZ signals and the electron overdensity field

#### 2.1.1 Derivation of the 21 cm power spectrum

The spin-flip transition between the hyperfine states of a neutral hydrogen H1 atom results in the emission or absorption of photons with a 21 cm wavelength (Field 1957). The abundance of H1 in the IGM makes this radiation an appealing probe of the large-scale structure of the Universe and a scientific goal for radio interferometry telescopes such as the SKA-Low.

The 21 cm signal from the EoR is seen against the Rayleigh-Jeans tail of the CMB. Assuming small optical depth and neglecting redshift space distortions, the brightness temperature of the signal can be written as

$$ \delta T_{21}(\mathbf{r}, z) = T_0(\mathbf{r}, z) x_{\text{H1}}(\mathbf{r}, z) [1 + \delta_b(\mathbf{r}, z)], \quad (1) $$

<footnote n="1">¹https://pole.uchicago.edu</footnote>

<footnote n="2">²https://www.esa.int/planck</footnote>

<footnote n="3">³www.lofar.org</footnote>

<footnote n="4">⁴https://reionization.org</footnote>

<footnote n="5">⁵www.mwatelescope.org</footnote>

<footnote n="6">⁶www.gmrt.ncra.tifr.res.in</footnote>

<footnote n="7">⁷www.skatelescope.org</footnote>
with (Pritchard & Loeb 2012)

$$T_0(r, z) \approx 27 \, \text{mK} \left(\frac{1+z}{10}\right)^{1/2} \left(\frac{T_s(r, z) - T_{\text{CMB}}(z)}{T_s(r, z)}\right) \times \left(\frac{\Omega_b}{0.044} \frac{h}{0.7}\right) \times \frac{\Omega_b h \sqrt{0.15(1+z)}}{0.023\sqrt{10\Omega_m}}. \quad (2)$$

Here, $\delta_b = \rho_b/\bar{\rho}_b - 1$ is the baryonic density fluctuation and the bar denotes a spatial average and $\delta_{H1}$ is the fluctuation of the neutral fraction $x_{H1}$ field, i.e. the fraction of the IGM baryons that are neutral. Naturally, $x_{H1} = 1 - x_{HII}$, where $x_{HII}$ is the ionization field. The pre-factor $T_0$ includes most of the spatially averaged information, including cosmological parameters, and the spin temperature information $T_s$. We assume the limit where during cosmic reionization the gas in the IGM is sufficiently heated by the first ionizing sources so $T_s \gg T_{\text{CMB}}$ and the spin-temperature dependence in $T_0$ is dropped.

Let us consider the electron fraction field

$$x_e \equiv f_H x_{HII} (1 + \delta_b), \quad (3)$$

where $f_H$ is the fractional quantity of electrons per hydrogen atom. If the first reionization of helium is considered, $f_H \simeq 1.08$. We distinguish the mass-weighted $x_m$ and volume-weighted $x_v$ average ionized fractions of H I atoms. By definition, $\langle x_e \rangle \equiv f_H x_m$, and, if we take the fluctuations of each element in equation (3), we have

$$x_m (1 + \delta_e) = x_v (1 + \delta_i) (1 + \delta_b), \quad (4)$$

where we refer to the density and ionization field perturbations as ‘b’ and ‘i’, respectively. Including these definitions in equation (1) and using the fact that $x_{H1} = 1 - x_{HII}$, we have

$$\begin{aligned} \frac{\delta T_{21}}{T_0(z)} &= (1 + \delta_b) - x_m(1 + \delta_e) \\ &= (1 + \delta_b)[1 - x_v(1 + \delta_i)], \end{aligned} \quad (5)$$

where the spatial- and time-dependence have been omitted for simplicity. Taking the Fourier transform and averaging, we find

$$\begin{aligned} \frac{P_{21}(k, z)}{T_0(z)^2} &= P_{bb}(k, z) + x_m(z)^2 P_{ee}(k, z) \\ &\quad - 2x_v(z) (P_{bb}(k, z) + P_{bi}(k, z) + P_{bi,b}(k, z)), \end{aligned} \quad (6)$$

where $P_{21}$ is the power spectrum of the 21 cm signal fluctuations, $P_{ee}$ of the electron density fluctuations, $P_{bi}$ of the ionization fluctuations, $P_{bb}$ of the baryon overdensity, and their respective cross-terms. Note that the cross-terms in the brackets can also be re-expressed as the cross-correlation between the 21 cm signal and the baryonic density field following Georgiev et al. (2022).

To avoid modelling the cross-terms and the three-point correlations, we simplify equation (6) as

$$\frac{P_{21}(k, z)}{T_0(z)^2} = [1 - 2x_v(z)] P_{bb}(k, z) + x_m(z)^2 P_{ee}(k, z). \quad (7)$$

Some further simplifications of equation (7) are assumed. We identify the mass-weighted $x_m$ and the volume-weighted $x_v$ ionized fractions and approximate the baryon power spectrum as a biased DM power spectrum $P_{bb}(k, z) = b_{\delta b}^2(k) P_{\delta\delta}(k, z)$, such that the fully simplified expression used in this work, unless stated otherwise, is

$$\frac{P_{21}(k, z)}{T_0(z)^2} = [1 - 2x_v(z)] b_{\delta b}(k)^2 P_{\delta\delta}(k, z) + x_v(z)^2 P_{ee}(k, z). \quad (8)$$

We will investigate the limitations of these assumptions in Section 4.1. Equation (8) yields a relation between the 21 cm power spectrum and the electron power spectrum $P_{ee}$. We will now look for a similar relation for the kSZ angular power spectrum.

### 2.1.2 Derivation of the patchy kSZ angular power spectrum

The angular power spectrum of the pkSZ effect at multipole $\ell$ can be derived from the electron density power spectrum $P_{ee}$ under some assumptions,<fnref n="8" /> such that (Mesinger, McQuinn & Spergel 2012; Gorce et al. 2020)

$$C_{\ell} = \frac{8\pi^2}{(2\ell + 1)^3} \frac{\sigma_T^2}{c^2} \int \frac{\bar{n}_e(z)^2}{(1+z)^2} \Delta_{B,e}^2 (\ell/\eta, z) e^{-2\tau(z)} \eta(z) \frac{d\eta}{dz} dz, \quad (9)$$

with $\sigma_T$ being the Thomson scattering cross-section, $\bar{n}_e(z)$ and $\eta(z)$ are the mean electron density and the comoving distance for redshift $z$, respectively. Moreover, $P_{B,e}$ is the power spectrum of the curl component of the momentum field $q_{B,e}$ such that $(2\pi)^3 P_{B,e} \delta_D(k - k') = \langle \tilde{q}_{B,e}(k) \tilde{q}_{B,e}^*(k') \rangle$, where $\delta_D$ is the Dirac delta function, the tilde denotes a Fourier transform, the asterisk a complex conjugate, and we define the dimensionless power spectrum, for a given field $a$ and $b$ at a certain wave number $k$, as $\Delta_{a,b}^2(k) = k^3 P_{a,b}(k)/(2\pi^2)$. We have

$$\begin{aligned} \frac{\langle \tilde{q}_{B,e}(k) \tilde{q}_{B,e}^*(k') \rangle}{(2\pi\nu)^3 \delta_D(|k-k'|)} &= \frac{2\pi^2}{k^3} \Delta_{B,e}^2(k, z) \\ &= \frac{1}{(2\pi\nu)^3} \int d^3k' \left[ (1-\mu^2) P_{ee}(|k-k'|) P_{vv}(k') - \frac{(1-\mu^2)k'}{|k-k'|} P_{ev}(|k-k'|) P_{ee}(k') \right]. \end{aligned} \quad (10)$$

where $\mu = \hat{k} \cdot \hat{k}'$ and the z-dependencies have been omitted for simplicity. $P_{ee}(k, z)$ is the power spectrum of the free electrons density fluctuations and $P_{ev}$ is the free electrons density–velocity cross-spectrum. The latter is computed as

$$P_{ee}(k, z) = \frac{f\dot{a}(z)}{k} b_{\delta e}(k, z) P_{\delta\delta}^{\text{lin}}(k, z), \quad (11)$$

where $a$ is the scale factor, $f$ the linear growth rate, and the bias is defined by the ratio $b_{\delta e}(k, z)^2 \equiv P_{ee}(k, z)/P_{\delta\delta}(k, z)$.

For a given electron power spectrum, both the angular pkSZ and the spherical 21 cm power spectra can thus be derived and it is this relationship we will explore in this paper. Note that, reciprocally, the pkSZ is effectively an integral of the 21 cm power spectrum, equation (7) can be used to reconstruct the former from the latter, a potential that we investigate in Appendix A. In Section 4.2, we will use this relationship to perform a joined fit of 21 cm and pkSZ power-spectrum measurements. We now turn to the model used to relate the electron power spectrum to reionization.

### 2.1.3 Electron power spectrum

We assume the $P_{ee}(k, z)$ parametrization introduced in Gorce et al. (2020), that is

$$P_{ee}(k, z) = f_H [1 - x_v(z)] \times \frac{\alpha_0 x_v(z)^{-1/5}}{1 + [k/\kappa]^3 x_v(z)} + x_v(z) b_{\delta b}(k, z)^2 P_{\delta\delta}(k, z), \quad (12)$$

where $\kappa$ is the electron drop-off frequency, which can be related to the typical size of ionized bubbles during the EoR, $\log_{10}(\alpha_0)$ is the large-scale $P_{ee}$ amplitude, related to the variance of the electron field during reionization. The redshift-independent baryon–DM bias $b_{\delta b}(k)$ is

<footnote n="8">⁸These assumptions include the Limber approximation, the assumption that the velocity power spectrum is a biased linear matter power spectrum, and the omission of third- and fourth-order correlation terms (Gorce et al. 2020). The latter implies variations of the order of 10 per cent (~0.05 µK$^2$) in the patchy kSZ amplitude (see Alvarez 2016 for details).</footnote>
with (Pritchard & Loeb 2012)

$$T_0(r, z) \approx 27 \, \text{mK} \left(\frac{1+z}{10}\right)^{1/2} \left(\frac{T_s(r, z) - T_{\text{CMB}}(z)}{T_s(r, z)}\right) \times \left(\frac{\Omega_b}{0.044} \frac{h}{0.7}\right) \times \frac{\Omega_b h \sqrt{0.15(1+z)}}{0.023\sqrt{10\Omega_m}}. \quad (2)$$

Here, $\delta_b = \rho_b/\bar{\rho}_b - 1$ is the baryonic density fluctuation and the bar denotes a spatial average and $\delta_{H1}$ is the fluctuation of the neutral fraction $x_{H1}$ field, i.e. the fraction of the IGM baryons that are neutral. Naturally, $x_{H1} = 1 - x_{HII}$, where $x_{HII}$ is the ionization field. The pre-factor $T_0$ includes most of the spatially averaged information, including cosmological parameters, and the spin temperature information $T_s$. We assume the limit where during cosmic reionization the gas in the IGM is sufficiently heated by the first ionizing sources so $T_s \gg T_{\text{CMB}}$ and the spin-temperature dependence in $T_0$ is dropped.

Let us consider the electron fraction field

$$x_e \equiv f_H x_{HII} (1 + \delta_b), \quad (3)$$

where $f_H$ is the fractional quantity of electrons per hydrogen atom. If the first reionization of helium is considered, $f_H \simeq 1.08$. We distinguish the mass-weighted $x_m$ and volume-weighted $x_v$ average ionized fractions of H I atoms. By definition, $\langle x_e \rangle \equiv f_H x_m$, and, if we take the fluctuations of each element in equation (3), we have

$$x_m (1 + \delta_e) = x_v (1 + \delta_i) (1 + \delta_b), \quad (4)$$

where we refer to the density and ionization field perturbations as ‘b’ and ‘i’, respectively. Including these definitions in equation (1) and using the fact that $x_{H1} = 1 - x_{HII}$, we have

$$\frac{\delta T_{21}}{T_0(z)} = (1 + \delta_b) - x_m(1 + \delta_e) \\ = (1 + \delta_b)[1 - x_v(1 + \delta_i)], \quad (5)$$

where the spatial- and time-dependence have been omitted for simplicity. Taking the Fourier transform and averaging, we find

$$\frac{P_{21}(k, z)}{T_0(z)^2} = P_{bb}(k, z) + x_m(z)^2 P_{ee}(k, z) - 2x_v(z) (P_{bb}(k, z) + P_{bi}(k, z) + P_{bi,b}(k, z)), \quad (6)$$

where $P_{21}$ is the power spectrum of the 21 cm signal fluctuations, $P_{ee}$ of the electron density fluctuations, $P_{bi}$ of the ionization fluctuations, $P_{bb}$ of the baryon overdensity, and their respective cross-terms. Note that the cross-terms in the brackets can also be re-expressed as the cross-correlation between the 21 cm signal and the baryonic density field following Georgiev et al. (2022).

To avoid modelling the cross-terms and the three-point correlations, we simplify equation (6) as

$$\frac{P_{21}(k, z)}{T_0(z)^2} = [1 - 2x_v(z)] P_{bb}(k, z) + x_m(z)^2 P_{ee}(k, z). \quad (7)$$

Some further simplifications of equation (7) are assumed. We identify the mass-weighted $x_m$ and the volume-weighted $x_v$ ionized fractions and approximate the baryon power spectrum as a biased DM power spectrum $P_{bb}(k, z) = b_{\delta b}^2(k) P_{\delta\delta}(k, z)$, such that the fully simplified expression used in this work, unless stated otherwise, is

$$\frac{P_{21}(k, z)}{T_0(z)^2} = [1 - 2x_v(z)] b_{\delta b}(k)^2 P_{\delta\delta}(k, z) + x_v(z)^2 P_{ee}(k, z). \quad (8)$$

We will investigate the limitations of these assumptions in Section 4.1. Equation (8) yields a relation between the 21 cm power spectrum and the electron power spectrum $P_{ee}$. We will now look for a similar relation for the kSZ angular power spectrum.

### 2.1.2 Derivation of the patchy kSZ angular power spectrum

The angular power spectrum of the pkSZ effect at multipole $l$ can be derived from the electron density power spectrum $P_{ee}$ under some assumptions,<fnref n="8" /> such that (Mesinger, McQuinn & Spergel 2012; Gorce et al. 2020)

$$C_l = \frac{8\pi^2}{(2\ell + 1)^3} \frac{\sigma_T^2}{c^2} \int \frac{\bar{n}_e(z)^2}{(1+z)^2} \Delta_{B,e}^2 (\ell/\eta, z) e^{-2\tau(z)} \eta(z) \frac{d\eta}{dz} dz, \quad (9)$$

with $\sigma_T$ being the Thomson scattering cross-section, $\bar{n}_e(z)$ and $\eta(z)$ are the mean electron density and the comoving distance for redshift $z$, respectively. Moreover, $P_{B,e}$ is the power spectrum of the curl component of the momentum field $q_{B,e}$ such that $(2\pi)^3 P_{B,e} \delta_D(k - k') = \langle \tilde{q}_{B,e}(k) \tilde{q}_{B,e}^*(k') \rangle$, where $\delta_D$ is the Dirac delta function, the tilde denotes a Fourier transform, the asterisk a complex conjugate, and we define the dimensionless power spectrum, for a given field $a$ and $b$ at a certain wave number $k$, as $\Delta_{a,b}^2(k) = k^3 P_{a,b}(k)/(2\pi^2)$. We have

$$\frac{\langle \tilde{q}_{B,e}(k) \tilde{q}_{B,e}^*(k') \rangle}{(2\pi\nu)^3 \delta_D(|k-k'|)} \equiv \frac{2\pi^2}{k^3} \Delta_{B,e}^2(k, z) = \frac{1}{(2\pi\nu)^3} \int d^3k' \left[ (1-\mu^2) P_{ee}(|k-k'|) P_{vv}(k') - \frac{(1-\mu^2)k'}{|k-k'|} P_{ev}(|k-k'|) P_{ee}(k') \right], \quad (10)$$

where $\mu = \hat{k} \cdot \hat{k}'$ and the $z$-dependencies have been omitted for simplicity. $P_{ee}(k, z)$ is the power spectrum of the free electrons density fluctuations and $P_{ev}$ is the free electrons density-velocity cross-spectrum. The latter is computed as

$$P_{ev}(k, z) = \frac{f \dot{a}(z)}{k} b_{\delta e}(k, z) P_{\delta\delta}^{\text{lin}}(k, z), \quad (11)$$

where $a$ is the scale factor, $f$ the linear growth rate, and the bias is defined by the ratio $b_{\delta e}(k, z)^2 \equiv P_{ee}(k, z)/P_{\delta\delta}(k, z)$.

For a given electron power spectrum, both the angular pkSZ and the spherical 21 cm power spectra can thus be derived and it is this relationship we will explore in this paper. Note that, reciprocally, the pkSZ is effectively an integral of the 21 cm power spectrum, equation (7) can be used to reconstruct the former from the latter, a potential that we investigate in Appendix A. In Section 4.2, we will use this relationship to perform a joined fit of 21 cm and pkSZ power-spectrum measurements. We now turn to the model used to relate the electron power spectrum to reionization.

### 2.1.3 Electron power spectrum

We assume the $P_{ee}(k, z)$ parametrization introduced in Gorce et al. (2020), that is

$$P_{ee}(k, z) = f_H [1 - x_v(z)] \times \frac{\alpha_0 x_v(z)^{-1/5}}{1 + [k/\kappa]^3 x_v(z)} + x_v(z) b_{\delta b}(k, z)^2 P_{\delta\delta}(k, z), \quad (12)$$

where $\kappa$ is the electron drop-off frequency, which can be related to the typical size of ionized bubbles during the EoR, $\log_{10}(\alpha_0)$ is the large-scale $P_{ee}$ amplitude, related to the variance of the electron field during reionization. The redshift-independent baryon–DM bias $b_{\delta b}(k)$ is

<footnote n="8">⁸These assumptions include the Limber approximation, the assumption that the velocity power spectrum is a biased linear matter power spectrum, and the omission of third- and fourth-order correlation terms (Gorce et al. 2020). The latter implies variations of the order of 10 per cent (~0.05 µK²) in the patchy kSZ amplitude (see Alvarez 2016 for details).</footnote>
<table><caption>**Table 1.** Reference cosmological and reionization parameters used to generate the mock data, based on the tuned to best fit the RSAGE simulation. All parameters apart from $z_{re}$, $z_{end}$, $\log_{10}(\alpha_0)$, and $\kappa$ are fixed.</caption><thead><tr><th>Parameter</th><th>Ref. value</th><th>Prior</th></tr></thead><tbody><tr><td>$h$</td><td>0.681</td><td>NA</td></tr><tr><td>$\Omega_b$</td><td>0.050</td><td>NA</td></tr><tr><td>$\Omega_m$</td><td>0.302</td><td>NA</td></tr><tr><td>$n_s$</td><td>0.96</td><td>NA</td></tr><tr><td>$A_s$</td><td>2.10 × 10<sup>-9</sup></td><td>NA</td></tr><tr><td>$z_{re}$</td><td>7.37</td><td>[5.0, 10.0]</td></tr><tr><td>$z_{end}$</td><td>6.15</td><td>[4.5, 9.0]</td></tr><tr><td>$d\bar{z} = z_{re} - z_{end}$</td><td>1.22</td><td>[0.5, 5.5]</td></tr><tr><td>$\log(\alpha_0 \text{ Mpc}^{-3})$</td><td>3.12</td><td>[2.5, 4.5]</td></tr><tr><td>$\kappa/\text{Mpc}^{-1}$</td><td>0.145</td><td>[0.05, 0.25]</td></tr></tbody></table>

given by the adapted Shaw, Rudd & Nagai (2012) parametrization, that is

$$b_{\delta b}(k)^2 = \frac{1}{2} \left[ e^{-k/k_f} + \frac{1}{1 + (gk/k_f)^2} \right], \quad (13)$$

where $k_f = 9.4 \text{ Mpc}^{-1}$ and $g = 0.5$ are constant with redshift and calibrated on the EMMA simulation (Aubert, Deparis & Ocvirk 2015), which includes coupled radiative transfer and hydrodynamics and is therefore sensitive to the thermal and reionization history. The global reionization history $x_v(z)$ in equation (12) is defined according to the following parametrization (Douspis et al. 2015; Gorce, Douspis & Salvati 2022):

$$z < z_{\text{end}}, \quad x_v(z) = 1, \quad (14)$$

$$z \ge z_{\text{end}}, \quad x_v(z) = \left( \frac{z_{\text{early}} - z}{z_{\text{early}} - z_{\text{end}}} \right)^{\alpha}, \quad (15)$$

where $z_{\text{early}} = 20$ is the redshift for which $x_v(z_{\text{early}}) \approx 10^{-4}$, the ionization leftover from recombination. Random forests are used to speed up computations and predict the pkSZ power spectrum given the parameter set; see Gorce, Douspis & Salvati (2022) for more details on the training and testing of the random forests. Note that for both kSZ and 21 cm derivations, the required matter power spectrum is obtained using the Boltzmann integrator CAMB<fnref n="9" /> (Lewis, Challinor & Lasenby 2000; Howlett et al. 2012).

We now have the full framework enabling us to build the spherical 21 cm and the angular kSZ power spectra given a reionization model, which we will use to jointly fit reionization parameters to mock measurements of these power spectra.

## 2.2 Fits to mock data

### 2.2.1 MCMC sampling

In Section 4, we perform a joint fit of mock pkSZ and 21 cm data points, derived as described in Section 2.1. To do so, we use a version of the COSMOMC MCMC sampler (Lewis & Bridle 2002; Lewis 2013), modified as described in Gorce, Douspis & Salvati (2022). Rather than sampling the Thomson optical depth $\tau$, we fit for the reionization mid- and endpoint, $z_{re}$ and $z_{end}$, as well as for the reionization morphology parameters $\log_{10}(\alpha_0)$ and $\kappa$, defined in equation (12). The model parameters are listed in Table 1. The assumed (flat) prior range is given for sampled parameters and the

<table><caption>**Table 2.** Parameters describing the mock observations by different telescopes used to derive our 21 cm power spectrum errors (equation 17).</caption><thead><tr><th>Parameters</th><th>SKA-Low</th><th>MWA</th><th>LOFAR</th><th>HERA</th></tr></thead><tbody><tr><td>$B$ (MHz)</td><td>10</td><td>10</td><td>10</td><td>10</td></tr><tr><td>$N_{\text{ant}}$</td><td>224</td><td>128</td><td>12</td><td>350</td></tr><tr><td>$A_{\text{eff}}$ (m<sup>2</sup>)</td><td>600</td><td>21.5</td><td>804</td><td>150</td></tr><tr><td>$R_{\text{core}}$ (m)</td><td>500</td><td>150</td><td>150</td><td>150</td></tr><tr><td>$\theta$ (deg<sup>2</sup>)</td><td>3<sup>2</sup></td><td>24.7<sup>2</sup></td><td>3.7<sup>2</sup></td><td>8<sup>2</sup></td></tr></tbody></table>

value of fixed parameters is listed. For each set of model parameters, we generate the corresponding reionization history and electron power spectrum $P_{ee}$. With the latter, we obtain the dimensionless 21 cm power spectrum $\Delta_{21}^2$ applying equation (7), as well as the pkSZ angular power spectrum according to equation (9), at each iteration of the sampler. We evaluate the agreement of our model with (mock) data by assuming two independent Gaussian likelihoods with uncorrelated errors for each data set:

$$\log \mathcal{L}_{\text{tot}} = - \frac{(\Delta_{21, \text{true}}^2 - \Delta_{21, \text{model}}^2)^2}{2(\Delta_{\text{noise}}^2 + \Delta_{\text{var}}^2)} - \frac{(C_{\ell, \text{true}}^{\text{pkSZ}} - C_{\ell, \text{model}}^{\text{pkSZ}})^2}{2\sigma_{\text{pkSZ}}^2}, \quad (16)$$

where uncertainty in the measurements $\Delta_{\text{noise}}^2$, $\Delta_{\text{var}}^2$, and $\sigma_{\text{pkSZ}}$ are calculated in Section 2.2.2. The true model, used to obtain the mock data points, is generated from equations (7) and (9) for the parameters in Table 1.

We employ the Gelman–Rubin test (Gelman & Rubin 1992) to assess the convergence of the MCMC. The *R* parameter used in this method represents the variance of the chain means compared to the mean of chain variances. In the cases where *R* ≲ 10<sup>-2</sup>, we consider the chains to have converged. All parameter values in this work are reported as the marginalized posterior probability’s maximum, which is more suitable for skewed distributions. Confidence intervals correspond to intervals with the highest probability density at 68 per cent.

### 2.2.2 Error estimation

Following equation (11) from Mellema et al. (2013), also used in Koopmans et al. (2015), we estimate the dimensionless noise power spectrum for different experiments and observation strategies with

$$\Delta_{\text{noise}}^2(k, z) = k^{3/2} \lambda_{21}(z) \times \frac{2}{\pi} \sqrt{\frac{D_c^2(z) \Delta D_c}{A_{\text{eff}}}} \frac{T_{\text{sys}}(z)^2}{B} \frac{A_{\text{core}} A_{\text{eff}}}{A_{\text{coll}}^2}, \quad (17)$$

where $\lambda_{21}(z)$ is the redshifted 21 cm wavelength, *B* is the bandwidth centred on redshift *z*, and $\Delta D_c$ is its length in comoving Mpc. We write the system noise $T_{\text{sys}} = 100 + 300(\nu_{\text{obs}}/\nu_0)^{-2.55}$ K, for $\nu_0 = 150$ MHz and $\nu_{\text{obs}}(z) = c/\lambda_{21}(z)$ is the observed 21 cm frequency. The total integration time is $t_{\text{int}}$, $D_c(z)$ is the comoving distance to redshift *z*, and $n_{\text{base}}$ the number of baselines, roughly equal to the number of antennas squared $N_a^2$. The total collecting area of the telescope is $A_{\text{coll}}$, such that $A_{\text{coll}} = N_a \pi R_a^2$ where $R_a$ is the radius of the antenna. We have $A_{\text{core}}$ the core area of the array and $A_{\text{eff}}$ the effective collecting area of each antenna. We use the values presented in Table 2 to estimate the noise for four different experiments: MWA, HERA, LOFAR, and SKA-Low. Note that for the latter and LOFAR, we consider each station as a single antenna.

We also include the contribution of sample variance to the uncertainty of the measurement following the relation (see equations 9 and

<footnote n="9">⁹Available at https://github.com/cmbant/CAMB.</footnote>
<table><caption>**Table 1.** Reference cosmological and reionization parameters used to generate the mock data, based on the tuned to best fit the RSAGE simulation. All parameters apart from $z_{re}$, $z_{end}$, $\log_{10}(\alpha_0)$, and $\kappa$ are fixed.</caption><thead><tr><th>Parameter</th><th>Ref. value</th><th>Prior</th></tr></thead><tbody><tr><td>$h$</td><td>0.681</td><td>NA</td></tr><tr><td>$\Omega_b$</td><td>0.050</td><td>NA</td></tr><tr><td>$\Omega_m$</td><td>0.302</td><td>NA</td></tr><tr><td>$n_s$</td><td>0.96</td><td>NA</td></tr><tr><td>$A_s$</td><td>2.10 × 10<sup>-9</sup></td><td>NA</td></tr><tr><td>$z_{re}$</td><td>7.37</td><td>[5.0, 10.0]</td></tr><tr><td>$z_{end}$</td><td>6.15</td><td>[4.5, 9.0]</td></tr><tr><td>$d\bar{z} = z_{re} - z_{end}$</td><td>1.22</td><td>[0.5, 5.5]</td></tr><tr><td>$\log(\alpha_0 \text{ Mpc}^{-3})$</td><td>3.12</td><td>[2.5, 4.5]</td></tr><tr><td>$\kappa/\text{Mpc}^{-1}$</td><td>0.145</td><td>[0.05, 0.25]</td></tr></tbody></table>

given by the adapted Shaw, Rudd & Nagai (2012) parametrization, that is

$$b_{\delta b}(k)^2 = \frac{1}{2} \left[ e^{-k/k_f} + \frac{1}{1 + (gk/k_f)^2} \right], \quad (13)$$

where $k_f = 9.4 \text{ Mpc}^{-1}$ and $g = 0.5$ are constant with redshift and calibrated on the EMMA simulation (Aubert, Deparis & Ocvirk 2015), which includes coupled radiative transfer and hydrodynamics and is therefore sensitive to the thermal and reionization history. The global reionization history $x_v(z)$ in equation (12) is defined according to the following parametrization (Douspis et al. 2015; Gorce, Douspis & Salvati 2022):

$$z < z_{\text{end}}, \quad x_v(z) = 1, \quad (14)$$

$$z \ge z_{\text{end}}, \quad x_v(z) = \left( \frac{z_{\text{early}} - z}{z_{\text{early}} - z_{\text{end}}} \right)^{\alpha}, \quad (15)$$

where $z_{\text{early}} = 20$ is the redshift for which $x_v(z_{\text{early}}) \approx 10^{-4}$, the ionization leftover from recombination. Random forests are used to speed up computations and predict the pkSZ power spectrum given the parameter set; see Gorce, Douspis & Salvati (2022) for more details on the training and testing of the random forests. Note that for both kSZ and 21 cm derivations, the required matter power spectrum is obtained using the Boltzmann integrator CAMB<fnref n="9" /> (Lewis, Challinor & Lasenby 2000; Howlett et al. 2012).

We now have the full framework enabling us to build the spherical 21 cm and the angular kSZ power spectra given a reionization model, which we will use to jointly fit reionization parameters to mock measurements of these power spectra.

## 2.2 Fits to mock data

### 2.2.1 MCMC sampling

In Section 4, we perform a joint fit of mock pkSZ and 21 cm data points, derived as described in Section 2.1. To do so, we use a version of the COSMOMC MCMC sampler (Lewis & Bridle 2002; Lewis 2013), modified as described in Gorce, Douspis & Salvati (2022). Rather than sampling the Thomson optical depth $\tau$, we fit for the reionization mid- and endpoint, $z_{re}$ and $z_{end}$, as well as for the reionization morphology parameters $\log_{10}(\alpha_0)$ and $\kappa$, defined in equation (12). The model parameters are listed in Table 1. The assumed (flat) prior range is given for sampled parameters and the

<table><caption>**Table 2.** Parameters describing the mock observations by different telescopes used to derive our 21 cm power spectrum errors (equation 17).</caption><thead><tr><th>Parameters</th><th>SKA-Low</th><th>MWA</th><th>LOFAR</th><th>HERA</th></tr></thead><tbody><tr><td>$B$ (MHz)</td><td>10</td><td>10</td><td>10</td><td>10</td></tr><tr><td>$N_{\text{ant}}$</td><td>224</td><td>128</td><td>12</td><td>350</td></tr><tr><td>$A_{\text{eff}}$ (m<sup>2</sup>)</td><td>600</td><td>21.5</td><td>804</td><td>150</td></tr><tr><td>$R_{\text{core}}$ (m)</td><td>500</td><td>150</td><td>150</td><td>150</td></tr><tr><td>$\theta$ (deg<sup>2</sup>)</td><td>3<sup>2</sup></td><td>24.7<sup>2</sup></td><td>3.7<sup>2</sup></td><td>8<sup>2</sup></td></tr></tbody></table>

value of fixed parameters is listed. For each set of model parameters, we generate the corresponding reionization history and electron power spectrum $P_{ee}$. With the latter, we obtain the dimensionless 21 cm power spectrum $\Delta_{21}^2$ applying equation (7), as well as the pkSZ angular power spectrum according to equation (9), at each iteration of the sampler. We evaluate the agreement of our model with (mock) data by assuming two independent Gaussian likelihoods with uncorrelated errors for each data set:

$$\log \mathcal{L}_{\text{tot}} = - \frac{(\Delta_{21, \text{true}}^2 - \Delta_{21, \text{model}}^2)^2}{2(\Delta_{\text{noise}}^2 + \Delta_{\text{var}}^2)} - \frac{(C_{\ell, \text{true}}^{\text{pkSZ}} - C_{\ell, \text{model}}^{\text{pkSZ}})^2}{2\sigma_{\text{pkSZ}}^2}, \quad (16)$$

where uncertainty in the measurements $\Delta_{\text{noise}}^2$, $\Delta_{\text{var}}^2$, and $\sigma_{\text{pkSZ}}$ are calculated in Section 2.2.2. The true model, used to obtain the mock data points, is generated from equations (7) and (9) for the parameters in Table 1.

We employ the Gelman–Rubin test (Gelman & Rubin 1992) to assess the convergence of the MCMC. The *R* parameter used in this method represents the variance of the chain means compared to the mean of chain variances. In the cases where $R \lesssim 10^{-2}$, we consider the chains to have converged. All parameter values in this work are reported as the marginalized posterior probability’s maximum, which is more suitable for skewed distributions. Confidence intervals correspond to intervals with the highest probability density at 68 per cent.

### 2.2.2 Error estimation

Following equation (11) from Mellema et al. (2013), also used in Koopmans et al. (2015), we estimate the dimensionless noise power spectrum for different experiments and observation strategies with

$$\Delta_{\text{noise}}^2(k, z) = k^{3/2} \lambda_{21}(z) \times \frac{2}{\pi} \sqrt{\frac{D_c^2(z) \Delta D_c}{A_{\text{eff}}}} \frac{T_{\text{sys}}(z)^2}{B} \frac{A_{\text{core}} A_{\text{eff}}}{A_{\text{coll}}^2}, \quad (17)$$

where $\lambda_{21}(z)$ is the redshifted 21 cm wavelength, *B* is the bandwidth centred on redshift *z*, and $\Delta D_c$ is its length in comoving Mpc. We write the system noise $T_{\text{sys}} = 100 + 300(\nu_{\text{obs}}/\nu_0)^{-2.55}$ K, for $\nu_0 = 150$ MHz and $\nu_{\text{obs}}(z) = c/\lambda_{21}(z)$ is the observed 21 cm frequency. The total integration time is $t_{\text{int}}$, $D_c(z)$ is the comoving distance to redshift *z*, and $n_{\text{base}}$ the number of baselines, roughly equal to the number of antennas squared $N_a^2$. The total collecting area of the telescope is $A_{\text{coll}}$, such that $A_{\text{coll}} = N_a \pi R_a^2$ where $R_a$ is the radius of the antenna. We have $A_{\text{core}}$ the core area of the array and $A_{\text{eff}}$ the effective collecting area of each antenna. We use the values presented in Table 2 to estimate the noise for four different experiments: MWA, HERA, LOFAR, and SKA-Low. Note that for the latter and LOFAR, we consider each station as a single antenna.

We also include the contribution of sample variance to the uncertainty of the measurement following the relation (see equations 9 and

<footnote n="9">⁹Available at https://github.com/cmbant/CAMB.</footnote>
<figure><img src="image_1.png" /><figcaption>**Figure 1.** Measurement errors on the 21 cm power spectrum from RSAGE (seen in black) for 1000 h of integration with different interferometers, corresponding to different colours. These include noise and sample variance errors. The noise and sample variance errors across different k-scales are presented for redshifts $z = 6.5, 7.8$ as thick solid and dashed lines, respectively.</figcaption></figure>

10 of Mellema et al. 2013):

$$ \Delta_{\text{var}}^{2}(k, z) = 0.01 \left( \frac{k}{0.1 \text{ Mpc}^{-1}} \right)^{-3/2} \left( \frac{V}{1 \text{ eGpc}^{3}} \right)^{-1/2} \left( \frac{\epsilon}{0.5} \right)^{-1/2} \Delta_{21}^{2}, \quad (18) $$

where $\epsilon$ is the logarithmic binning of $\Delta k = \epsilon k$ and $V$ is the survey size, which can be expressed as

$$ V = 0.1 \text{ Gpc}^3 \times \left( \frac{\theta}{5^{\circ}} \right)^2 \left( \frac{B}{12 \text{ MHz}} \right) [(1+z)^{1/2} - 2], \quad (19) $$

and $\theta$ is the field of view.

An example of measurement errors for an integration time of 1000 h is presented in Fig. 1 for each radio telescope considered for two redshifts, chosen following the upper limits of Trott et al. (2020). The 21 cm power spectrum for each redshift is also included. Note that the sample variance dominates the total error budget for the SKA-Low at $k \lesssim 0.3 \text{ Mpc}^{-1}$ because of its strong k-dependency ($\Delta_{\text{var}}^2 \propto k^{-3/2}$). The opposite is true of the noise estimate for MWA, which has a larger field of view and as $\Delta_{\text{var}}^2 \propto 1/\theta$ has a lower error due to sample variance but a higher uncertainty on the noise due to its configuration.

Regarding the measurement errors on the pkSZ angular power spectrum, we choose a 10 per cent uncertainty on the data point throughout this work. We base this on the current constraint of the total kSZ amplitude $D_{3000}^{\text{pkSZ}} = 3.0 \pm 1.0 \text{ µK}^2$ (Reichardt et al. 2021), assuming that in future observations longer observation times will allow to reduce the noise, while numerous frequency channels and improved modelling will help characterizing and removing other CMB foregrounds (Maniyar, Béthermin & Lagache 2021; Douspis et al. 2022; Raghunathan & Omori 2023).

## 3 DATA

We briefly describe the RSAGE simulation used in this work (see Seiler et al. 2019, for a detailed description) to validate our simplified parametrization of the 21 cm power spectrum given in equation (7). The N-body data underlying RSAGE has been generated with the hydrodynamic KALI code (Seiler et al. 2018) with $2400^3$

<figure><img src="image_2.png" /><figcaption>**Figure 2.** Slices of the differential surface brightness temperature $\delta T_{21}$ for the RSAGE simulation. The upper panel is at a redshift of $z = 7.83$ and has a volume-averaged ionization fraction of $x_v = 0.37$, while for the bottom panel, the same values are $z = 6.52$ and $x_v = 0.87$. Both redshifts correspond to the mock data points used in Section 4.2.</figcaption></figure>

DM particles and a volume of $(160 \text{ Mpc})^3$. The radiative transfer is conducted using the seminumerical code CIFOG (Hutter 2018), accounting for star formation and feedback processes. We use the const version of the simulation, which assumes a constant escape fraction of ionizing photons $f_{\text{esc}} = 0.2$. In Fig. 2, we show simulation slices of the differential surface brightness temperature of the 21 cm signal. The cosmology used is consistent with Planck Collaboration I (2020) and summarized in Table 1. The reionization mid- and endpoint values given in Table 1 are obtained by fitting the asymmetric parametrization from equation (14) to the reionization history of the simulation. Similarly, the values of the morphology parameters $\log_{10}(\alpha_0)$ and $\kappa$ from Table 1 are derived by fitting the RSAGE electron density power spectrum with equation (12) (see Gorce et al. 2020, for details of the fit). We have checked that the baryon power spectrum from RSAGE is well described using equation (13) for the given values of $k_f$ and $g$. Overall, we find the model of the electron density power spectrum in equation (12) to be a good fit for the redshift and k-scales explored in this work. In the following section, we therefore limit our inspections to the accuracy of the 21 cm power-spectrum reconstruction.
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Measurement errors on the 21 cm power spectrum from RSAGE (seen in black) for 1000 h of integration with different interferometers, corresponding to different colours. These include noise and sample variance errors. The noise and sample variance errors across different k-scales are presented for redshifts $z = 6.5, 7.8$ as thick solid and dashed lines, respectively.</figcaption></figure>

10 of Mellema et al. 2013):

$$ \Delta_{\text{var}}^{2}(k, z) = 0.01 \left( \frac{k}{0.1 \text{ Mpc}^{-1}} \right)^{-3/2} \left( \frac{V}{1 \text{ eGpc}^{3}} \right)^{-1/2} \left( \frac{\epsilon}{0.5} \right)^{-1/2} \Delta_{21}^{2}, \quad (18) $$

where $\epsilon$ is the logarithmic binning of $\Delta k = \epsilon k$ and $V$ is the survey size, which can be expressed as

$$ V = 0.1 \text{ Gpc}^3 \times \left( \frac{\theta}{5^{\circ}} \right)^2 \left( \frac{B}{12 \text{ MHz}} \right) [(1+z)^{1/2} - 2], \quad (19) $$

and $\theta$ is the field of view.

An example of measurement errors for an integration time of 1000 h is presented in Fig. 1 for each radio telescope considered for two redshifts, chosen following the upper limits of Trott et al. (2020). The 21 cm power spectrum for each redshift is also included. Note that the sample variance dominates the total error budget for the SKA-Low at $k \lesssim 0.3 \text{ Mpc}^{-1}$ because of its strong k-dependency ($\Delta_{\text{var}}^2 \propto k^{-3/2}$). The opposite is true of the noise estimate for MWA, which has a larger field of view and as $\Delta_{\text{var}}^2 \propto 1/\theta$ has a lower error due to sample variance but a higher uncertainty on the noise due to its configuration.

Regarding the measurement errors on the pkSZ angular power spectrum, we choose a 10 per cent uncertainty on the data point throughout this work. We base this on the current constraint of the total kSZ amplitude $D_{3000}^{\text{pkSZ}} = 3.0 \pm 1.0 \text{ µK}^2$ (Reichardt et al. 2021), assuming that in future observations longer observation times will allow to reduce the noise, while numerous frequency channels and improved modelling will help characterizing and removing other CMB foregrounds (Maniyar, Béthermin & Lagache 2021; Douspis et al. 2022; Raghunathan & Omori 2023).

## 3 DATA

We briefly describe the RSAGE simulation used in this work (see Seiler et al. 2019, for a detailed description) to validate our simplified parametrization of the 21 cm power spectrum given in equation (7). The N-body data underlying RSAGE has been generated with the hydrodynamic KALI code (Seiler et al. 2018) with $2400^3$

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Slices of the differential surface brightness temperature $\delta T_{21}$ for the RSAGE simulation. The upper panel is at a redshift of $z = 7.83$ and has a volume-averaged ionization fraction of $x_v = 0.37$, while for the bottom panel, the same values are $z = 6.52$ and $x_v = 0.87$. Both redshifts correspond to the mock data points used in Section 4.2.</figcaption></figure>

DM particles and a volume of $(160 \text{ Mpc})^3$. The radiative transfer is conducted using the seminumerical code CIFOG (Hutter 2018), accounting for star formation and feedback processes. We use the const version of the simulation, which assumes a constant escape fraction of ionizing photons $f_{\text{esc}} = 0.2$. In Fig. 2, we show simulation slices of the differential surface brightness temperature of the 21 cm signal. The cosmology used is consistent with Planck Collaboration I (2020) and summarized in Table 1. The reionization mid- and endpoint values given in Table 1 are obtained by fitting the asymmetric parametrization from equation (14) to the reionization history of the simulation. Similarly, the values of the morphology parameters $\log_{10}(\alpha_0)$ and $\kappa$ from Table 1 are derived by fitting the RSAGE electron density power spectrum with equation (12) (see Gorce et al. 2020, for details of the fit). We have checked that the baryon power spectrum from RSAGE is well described using equation (13) for the given values of $k_f$ and $g$. Overall, we find the model of the electron density power spectrum in equation (12) to be a good fit for the redshift and k-scales explored in this work. In the following section, we therefore limit our inspections to the accuracy of the 21 cm power-spectrum reconstruction.
<figure><img src="image_3.png" /><figcaption>**Figure 3.** *Upper Panel:* The dimensionless 21 cm power spectrum at $k = 0.141 \text{ Mpc}^{-1}$ as a function of redshift. The power spectrum computed directly from the RSAGE simulation is represented by the black dashed line. The coloured lines correspond to different reconstruction levels (see the text for details). *recons no cross* corresponds to equation (8), which is used for the forecast in Section 4.2. Additionally, *recons no cross* (*with $x_m$*) corresponds to equation (8) but for $x_m \neq x_v$. In *recons with cross*, the cross-power spectrum is included: $-2x_v(z)P_{bi}(k, z)$. Lastly, *recons with third order* adds three-point power spectrum contribution to *recons with cross*: $-2x_v(z)P_{bi, b}(k, z)$, which corresponds to the full form of the 21 cm power spectrum in equation (6). *Middle Panel:* Ratio of each component in equation (6) against the total 21 cm power spectrum. In this panel, solid (dashed) lines represent a positive (negative) contribution. *Bottom Panel:* Reionization history of the RSAGE simulation.</figcaption></figure>

# 4 RESULTS

In this section, we first look at the model uncertainties associated with our simplified derivation of the 21 cm power spectrum from equation (8), given an electron power spectrum by comparing it to the full spectrum computed from the RSAGE simulation. In Section 4.2, we derive forecast constraints on reionization by applying this derivation to fit mock 21 cm and kSZ data points, using equations (8) and (12), to a reionization model with the parameters from Table 1. We showcase the advantages of combining the two data sets compared to a separate analysis. In Appendix C, we extend this analysis with 21 cm and kSZ data for the from RSAGE simulation to study the effect of the model assumptions.

## 4.1 Reconstructing the power spectra

We want to understand the role of each contributing term in equation (6), especially the ones we have discarded in our simplified $P_{21}$ reconstruction (equation 8), in order to assess its accuracy. We reconstruct and examine the 21 cm power spectrum from the RSAGE simulation described in Section 3. The top panel of Fig. 3 showcases the true $P_{21}$ plotted against different reconstruction precision levels:

(1) The *recons no cross* case corresponds to the highest simplification level, corresponding to equation (8), and the approximation

we will use in the remaining of this work. Here, it is computed with the true $P_{ee}(k, z)$ spectrum obtained from the simulation, and not equation (12); therefore, we only assess the accuracy of the 21 cm reconstruction. The contributions from both third-order terms and cross-terms are discarded. We also include an additional case, *recons no cross* (*using $x_m$*), where we lift the assumption of equating the mass-weighted and volume-weighted ionization fractions.

(2) For *recons with cross*, we add the cross-power spectra contributions: $-2x_v(z)P_{bi}(k, z)$,<fnref n="10" />

(3) For *recons with third order*, we further add the three-point power spectrum contribution: $-2x_v(z)P_{bi, b}(k, z)$.<fnref n="10" /> This corresponds to the full form of the 21 cm power spectrum, without simplification.

We look at the relative contribution of each of the terms in equation (6) to the 21 cm power spectrum in the RSAGE simulation. This is represented in the middle panel of Fig. 3 for $k = 0.141 \text{ Mpc}^{-1}$. We see that for $z > 10$, in the very early stages of EoR, most of the 21 cm power spectrum amplitude comes from the matter density. In this period, despite the low ionization level ($x_v \lesssim 10^{-4}$), the ionization and density fields are weakly correlated, and the role of higher order terms is negligible (Lidz et al. 2007): All reconstruction levels of $P_{21}$ look identical in the upper panel. Beyond $x_v \approx 10$ per cent, when the first luminous sources have reionized their local overdensities and the ionization field is correlated with the density field, $P_{21}$ is primarily determined by the difference between the electron density power spectrum and the $P_{bi}$ cross-spectrum. The third-order term $P_{bib}$, in red, plays a minor role during most of reionization, only gaining significance in the final stages. This term is initially negative and changes sign at the midpoint of reionization because of the change in correlation between the density and ionization fields (Georgiev et al. 2022). Conversely, for the density term, the change in sign is due to the $[1 - 2x_v(z)]$ factor present in equation (7).

We can form a clearer understanding of the impact of each term by examining the top panel of Fig. 3, where the power spectrum of the 21 cm signal is calculated according to equation (8), in blue, while the full model presented in equation (6), in green. We see that the simplified expression in the *recons no cross* model overestimates the true power throughout reionization. Early on, *recons no cross* follows the true *recons with third order* model. With the onset of reionization, a positive bias of an order of two appears between the true spectrum and its model, reaching a global maximum at the midpoint of reionization. Additionally, in the simplified model, reionization is delayed by four per cent compared to the true model. Hence, our simplified model will tend to overestimate $P_{21}$ and will be biased towards late reionization scenarios, the implications of which are discussed in Section 4.2. Additionally, we note that by removing the $x_v \equiv x_m$ assumption and replacing $x_v^2 P_{ee}$ by $x_m^2 P_{ee}$ in equation (7) results in a higher amplitude (seen in purple) compared to the *recons no cross* model (light blue) using the volume-weighted fraction. The increased bias of the model is due to the mass-weighted ionization fraction being larger than the volume-weighted ionization fraction $x_v > x_m$ at all redshifts in inside-out models of reionization (see e.g. Dixon et al. 2016). We have compared this across the k-scales accessible within the RSAGE simulation volume, and we find the $x_v \equiv x_m$ simplification leads to a $P_{21}$ model closer to the truth when excluding the higher order contributions.

<footnote n="10">¹⁰Note that this contribution can be positive or negative, as illustrated in Fig. 3.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 3.** *Upper Panel:* The dimensionless 21 cm power spectrum at $k = 0.141 \text{ Mpc}^{-1}$ as a function of redshift. The power spectrum computed directly from the RSAGE simulation is represented by the black dashed line. The coloured lines correspond to different reconstruction levels (see the text for details). *recons no cross* corresponds to equation (8), which is used for the forecast in Section 4.2. Additionally, *recons no cross* (*with $x_m$*) corresponds to equation (8) but for $x_m \neq x_v$. In *recons with cross*, the cross-power spectrum is included: $-2x_v(z)P_{bi}(k, z)$. Lastly, *recons with third order* adds three-point power spectrum contribution to *recons with cross*: $-2x_v(z)P_{bi, b}(k, z)$, which corresponds to the full form of the 21 cm power spectrum in equation (6). *Middle Panel:* Ratio of each component in equation (6) against the total 21 cm power spectrum. In this panel, solid (dashed) lines represent a positive (negative) contribution. *Bottom Panel:* Reionization history of the RSAGE simulation.</figcaption></figure>

# 4 RESULTS

In this section, we first look at the model uncertainties associated with our simplified derivation of the 21 cm power spectrum from equation (8), given an electron power spectrum by comparing it to the full spectrum computed from the RSAGE simulation. In Section 4.2, we derive forecast constraints on reionization by applying this derivation to fit mock 21 cm and kSZ data points, using equations (8) and (12), to a reionization model with the parameters from Table 1. We showcase the advantages of combining the two data sets compared to a separate analysis. In Appendix C, we extend this analysis with 21 cm and kSZ data for the from RSAGE simulation to study the effect of the model assumptions.

## 4.1 Reconstructing the power spectra

We want to understand the role of each contributing term in equation (6), especially the ones we have discarded in our simplified $P_{21}$ reconstruction (equation 8), in order to assess its accuracy. We reconstruct and examine the 21 cm power spectrum from the RSAGE simulation described in Section 3. The top panel of Fig. 3 showcases the true $P_{21}$ plotted against different reconstruction precision levels:

(1) The *recons no cross* case corresponds to the highest simplification level, corresponding to equation (8), and the approximation

we will use in the remaining of this work. Here, it is computed with the true $P_{ee}(k, z)$ spectrum obtained from the simulation, and not equation (12); therefore, we only assess the accuracy of the 21 cm reconstruction. The contributions from both third-order terms and cross-terms are discarded. We also include an additional case, *recons no cross* (*using $x_m$*), where we lift the assumption of equating the mass-weighted and volume-weighted ionization fractions.

(2) For *recons with cross*, we add the cross-power spectra contributions: $-2x_v(z)P_{bi}(k, z)$,<fnref n="10" />

(3) For *recons with third order*, we further add the three-point power spectrum contribution: $-2x_v(z)P_{bi, b}(k, z)$.<fnref n="10" /> This corresponds to the full form of the 21 cm power spectrum, without simplification.

We look at the relative contribution of each of the terms in equation (6) to the 21 cm power spectrum in the RSAGE simulation. This is represented in the middle panel of Fig. 3 for $k = 0.141 \text{ Mpc}^{-1}$. We see that for $z > 10$, in the very early stages of EoR, most of the 21 cm power spectrum amplitude comes from the matter density. In this period, despite the low ionization level ($x_v \lesssim 10^{-4}$), the ionization and density fields are weakly correlated, and the role of higher order terms is negligible (Lidz et al. 2007): All reconstruction levels of $P_{21}$ look identical in the upper panel. Beyond $x_v \approx 10$ per cent, when the first luminous sources have reionized their local overdensities and the ionization field is correlated with the density field, $P_{21}$ is primarily determined by the difference between the electron density power spectrum and the $P_{bi}$ cross-spectrum. The third-order term $P_{bib}$, in red, plays a minor role during most of reionization, only gaining significance in the final stages. This term is initially negative and changes sign at the midpoint of reionization because of the change in correlation between the density and ionization fields (Georgiev et al. 2022). Conversely, for the density term, the change in sign is due to the $[1 - 2x_v(z)]$ factor present in equation (7).

We can form a clearer understanding of the impact of each term by examining the top panel of Fig. 3, where the power spectrum of the 21 cm signal is calculated according to equation (8), in blue, while the full model presented in equation (6), in green. We see that the simplified expression in the *recons no cross* model overestimates the true power throughout reionization. Early on, *recons no cross* follows the true *recons with third order* model. With the onset of reionization, a positive bias of an order of two appears between the true spectrum and its model, reaching a global maximum at the midpoint of reionization. Additionally, in the simplified model, reionization is delayed by four per cent compared to the true model. Hence, our simplified model will tend to overestimate $P_{21}$ and will be biased towards late reionization scenarios, the implications of which are discussed in Section 4.2. Additionally, we note that by removing the $x_v \equiv x_m$ assumption and replacing $x_v^2 P_{ee}$ by $x_m^2 P_{ee}$ in equation (7) results in a higher amplitude (seen in purple) compared to the *recons no cross* model (light blue) using the volume-weighted fraction. The increased bias of the model is due to the mass-weighted ionization fraction being larger than the volume-weighted ionization fraction $x_v > x_m$ at all redshifts in inside-out models of reionization (see e.g. Dixon et al. 2016). We have compared this across the k-scales accessible within the RSAGE simulation volume, and we find the $x_v \equiv x_m$ simplification leads to a $P_{21}$ model closer to the truth when excluding the higher order contributions.

<footnote n="10">¹⁰Note that this contribution can be positive or negative, as illustrated in Fig. 3.</footnote>
Adding only the contribution of the cross-power spectrum in the calculation (*recons with cross*) does not help the reconstruction as much as one would hope. While the additional negative power from the cross-correlation decreases the bias of the estimate to that of the true power spectrum, the exclusion of the higher order power spectrum shifts the midpoint redshift and biases the result to an early end of reionization.

The role of the higher order terms in the derivation of the 21 cm power spectrum is non-trivial and varies across the *k*-scales and redshifts. Generally, the negative cross-correlation between the density and ionization fields weakens at higher *k*-scales as well as with the progression of reionization. We see that by excluding the higher order power spectrum and varying the *k*-scale, the *recons with cross* model converges to the *recons with no cross* model on large *k*-scales. In Appendix A, we look at the possibility of using equation (7) to reconstruct the pkSZ angular power spectrum from an interferometric measurement of the 21 cm power spectrum.

In summary, the different terms comprising our expression of the 21 cm power spectrum, described in equation (6), are non-trivial. The importance of each term presented in Fig. 3 is shown to evolve with redshift and *k*-scale. The model corresponding to the highest degree of simplification is the most promising for this work, as the amplitude of the simplified 21 cm power spectrum is positively biased on all *k*-scales, resulting in a positive bias on the end of reionization, which we can easily quantify. In contrast, including the cross-power spectrum $P_{bi}$ contribution to the model introduces an evolving amplitude bias with *k*-scale. Hence, we choose to perform our forecast with the *recons no cross* model while keeping its limitations in mind. We use this model to derive our true mock data points and each sampled model.

For the sake of clarity, we list below the different layers of assumptions for our model in equation (8) and how we address them:

(i) The baryon power spectrum follows equation (13) when using the CAMB-derived theory non-linear matter power spectrum. Comparing to the spectrum directly measured from RSAGE, we find a good match with the model across all redshifts and *k*-scales covered in this work.

(ii) The global reionization history follows equation (14). We have compared the evolution obtained with equation (14) and the parameters of Table 1 to the RSAGE history seen at the bottom panel of Fig. 3 and found to be a good match for the redshifts considered in this paper.

(iii) The mass- and volume-weighted ionization fractions $x_m$ and $x_v$ are identical. As seen in the upper panel of Fig. 3, the positive bias induced by excluding the higher order terms is slightly alleviated for the *recons no cross* model for the inside-out reionization models considered in this work, for which $x_v < x_m$.

(iv) The true electron power spectrum in RSAGE follows equation (12). This model has been compared against the RSAGE simulation within this work while the best-fitting morphology parameters are reported in Table 1 (see appendix B.2 of Gorce et al. 2020, for a detailed discussion).

We conclude this assessment of our model uncertainties by comparing, in Fig. 4, the true 21 cm power spectrum from RSAGE against the one generated by our fully simplified parametrization, given in equation (8). Overall, the parametrization performs remarkably well on the *k*-scales considered in this work ($k \sim 0.1$ Mpc$^{-1}$). At high *k*-scales ($k > 1.0$ Mpc$^{-1}$), the model in equation (12) scales as $k^{-3}$ and quickly approaches the baryon power spectrum. The result is a negative 21 cm power spectrum as can be observed in Fig. 4 at $z = 6.5$. Therefore, the model should not be used as is at very small scales.

<figure><img src="image_4.png" /><figcaption>**Figure 4.** The 21 cm power spectrum RSAGE (seen in black) compared against our model (seen in fuchsia) for redshifts $z = 6.5, 7.8$ ($x_v = 0.87, 0.37$) as thick solid and dashed lines, respectively. Our model 21 cm power spectrum has been generated using equation (8) for the reference parameters in Table 1 and the CAMB non-linear matter power spectrum.</figcaption></figure>

At low-*k*, differences between the model and the RSAGE power spectrum will be sensitive to sample variance due to the simulation box size and hence have not been studied within this work.

We discuss the effect of the aforementioned bias on our results in Section C, where we derive the mock data using the RSAGE simulation and no approximation. In Section 6, we discuss future improvements and possible avenues in modelling the higher order terms in our model.

## 4.2 Efficient probe combination

In this section, we fit the mid-point, endpoint, and morphology parameters of reionization to three mock data points: one measurement of the pkSZ angular power spectrum at multipole $\ell = 3000$ and two measurements of the 21 cm power spectrum. The choice of the data points is motivated by current upper limits on both observables (e.g. Trott et al. 2020; Reichardt et al. 2021; Gorce, Douspis & Salvati 2022), as well as our estimate of measurement errors from Section 2.2.2. The number of data points used is chosen to efficiently retrieve information on the EoR while intuitively illustrating the role of each measurement. Limiting the number of data points will also help us understand how many observations and of what quality are necessary to begin constraining the properties of reionization, in a context of gradual improvement of upper limits on the 21 cm power spectrum. In Section 5, we investigate how the quality of constraints depends on the redshift and scale chosen for our data points. We first forecast results obtained with 1000 h of SKA-Low observations before turning to mock observations by MWA and LOFAR.

### 4.2.1 Forecast with SKA-Low

Our primary scenario is based on detection of the 21 cm power spectrum at redshifts $z = 6.5$ and 7.8, both at $k = 0.50$ Mpc$^{-1}$, that is outside the foreground wedge (Liu, Parsons & Trott 2014), considering noise levels associated with 100 and 1000 h of observations with the SKA-Low. Additionally, a measurement of the pkSZ power spectrum is assumed at $\ell = 3000$ with a ten per cent error bar such that $D_{\ell=3000}^{\text{pkSZ}} = 0.86 \pm 0.09$ µK$^2$. For simplicity, we will refer to these cases as 1k2z, based on the choice of two observations of the 21 cm power spectrum at the same *k*-scale.
Adding only the contribution of the cross-power spectrum in the calculation (*recons with cross*) does not help the reconstruction as much as one would hope. While the additional negative power from the cross-correlation decreases the bias of the estimate to that of the true power spectrum, the exclusion of the higher order power spectrum shifts the midpoint redshift and biases the result to an early end of reionization.

The role of the higher order terms in the derivation of the 21 cm power spectrum is non-trivial and varies across the *k*-scales and redshifts. Generally, the negative cross-correlation between the density and ionization fields weakens at higher *k*-scales as well as with the progression of reionization. We see that by excluding the higher order power spectrum and varying the *k*-scale, the *recons with cross* model converges to the *recons with no cross* model on large *k*-scales. In Appendix A, we look at the possibility of using equation (7) to reconstruct the pkSZ angular power spectrum from an interferometric measurement of the 21 cm power spectrum.

In summary, the different terms comprising our expression of the 21 cm power spectrum, described in equation (6), are non-trivial. The importance of each term presented in Fig. 3 is shown to evolve with redshift and *k*-scale. The model corresponding to the highest degree of simplification is the most promising for this work, as the amplitude of the simplified 21 cm power spectrum is positively biased on all *k*-scales, resulting in a positive bias on the end of reionization, which we can easily quantify. In contrast, including the cross-power spectrum $P_{bi}$ contribution to the model introduces an evolving amplitude bias with *k*-scale. Hence, we choose to perform our forecast with the *recons no cross* model while keeping its limitations in mind. We use this model to derive our true mock data points and each sampled model.

For the sake of clarity, we list below the different layers of assumptions for our model in equation (8) and how we address them:

(i) The baryon power spectrum follows equation (13) when using the CAMB-derived theory non-linear matter power spectrum. Comparing to the spectrum directly measured from RSAGE, we find a good match with the model across all redshifts and *k*-scales covered in this work.

(ii) The global reionization history follows equation (14). We have compared the evolution obtained with equation (14) and the parameters of Table 1 to the RSAGE history seen at the bottom panel of Fig. 3 and found to be a good match for the redshifts considered in this paper.

(iii) The mass- and volume-weighted ionization fractions $x_m$ and $x_v$ are identical. As seen in the upper panel of Fig. 3, the positive bias induced by excluding the higher order terms is slightly alleviated for the *recons no cross* model for the inside-out reionization models considered in this work, for which $x_v < x_m$.

(iv) The true electron power spectrum in RSAGE follows equation (12). This model has been compared against the RSAGE simulation within this work while the best-fitting morphology parameters are reported in Table 1 (see appendix B.2 of Gorce et al. 2020, for a detailed discussion).

We conclude this assessment of our model uncertainties by comparing, in Fig. 4, the true 21 cm power spectrum from RSAGE against the one generated by our fully simplified parametrization, given in equation (8). Overall, the parametrization performs remarkably well on the *k*-scales considered in this work ($k \sim 0.1$ Mpc$^{-1}$). At high *k*-scales ($k > 1.0$ Mpc$^{-1}$), the model in equation (12) scales as $k^{-3}$ and quickly approaches the baryon power spectrum. The result is a negative 21 cm power spectrum as can be observed in Fig. 4 at $z=6.5$. Therefore, the model should not be used as is at very small scales.

<figure><img src="image_6.png" /><figcaption>**Figure 4.** The 21 cm power spectrum RSAGE (seen in black) compared against our model (seen in fuchsia) for redshifts $z = 6.5, 7.8$ ($x_v = 0.87, 0.37$) as thick solid and dashed lines, respectively. Our model 21 cm power spectrum has been generated using equation (8) for the reference parameters in Table 1 and the CAMB non-linear matter power spectrum.</figcaption></figure>

At low-*k*, differences between the model and the RSAGE power spectrum will be sensitive to sample variance due to the simulation box size and hence have not been studied within this work.

We discuss the effect of the aforementioned bias on our results in Section C, where we derive the mock data using the RSAGE simulation and no approximation. In Section 6, we discuss future improvements and possible avenues in modelling the higher order terms in our model.

## 4.2 Efficient probe combination

In this section, we fit the mid-point, endpoint, and morphology parameters of reionization to three mock data points: one measurement of the pkSZ angular power spectrum at multipole $\ell = 3000$ and two measurements of the 21 cm power spectrum. The choice of the data points is motivated by current upper limits on both observables (e.g. Trott et al. 2020; Reichardt et al. 2021; Gorce, Douspis & Salvati 2022), as well as our estimate of measurement errors from Section 2.2.2. The number of data points used is chosen to efficiently retrieve information on the EoR while intuitively illustrating the role of each measurement. Limiting the number of data points will also help us understand how many observations and of what quality are necessary to begin constraining the properties of reionization, in a context of gradual improvement of upper limits on the 21 cm power spectrum. In Section 5, we investigate how the quality of constraints depends on the redshift and scale chosen for our data points. We first forecast results obtained with 1000 h of SKA-Low observations before turning to mock observations by MWA and LOFAR.

### 4.2.1 Forecast with SKA-Low

Our primary scenario is based on detection of the 21 cm power spectrum at redshifts $z = 6.5$ and $7.8$, both at $k = 0.50$ Mpc$^{-1}$, that is outside the foreground wedge (Liu, Parsons & Trott 2014), considering noise levels associated with 100 and 1000 h of observations with the SKA-Low. Additionally, a measurement of the pkSZ power spectrum is assumed at $\ell = 3000$ with a ten per cent error bar such that $D_{\ell=3000}^{\text{pkSZ}} = 0.86 \pm 0.09$ µK$^2$. For simplicity, we will refer to these cases as 1k2z, based on the choice of two observations of the 21 cm power spectrum at the same *k*-scale.
<figure><img src="image_5.png" /><figcaption>**Figure 5.** Posterior distributions on reionization mid- and endpoint (upper panel) and the morphology parameters $\log_{10}(\alpha_0)$ and $\kappa$ (bottom panel), considering a detection of the pkSZ power spectrum at $\ell = 3000$ and of the 21 cm power spectrum at redshifts $z = 6.5, 7.8$, and $k = 0.50 \text{ Mpc}^{-1}$ for either 100 or 1000 h of integration with SKA-Low, in purple and in blue, respectively. Note that the axis ranges in both panels are smaller than the priors in Table 1 or the ones presented in following figures.</figcaption></figure>

In Fig. 5, we show the joint probability distributions for our four sampled parameters (see appendix B for a discussion on parameter degeneracies). The best-fitting parameters, their $1\sigma$ error, and Gelman–Rubin convergence parameters are presented in Table 3.<fnref n="11" />

Our results indicate that we can recover the true values of $z_{re}$ and $z_{end}$ in both the 1k2z cases, despite the order of magnitude difference in integration time between them. Moreover, we independently constrain the value of the Thomson optical depth with a deviation of less than one per cent of the true value, $\tau = 0.066 \pm 0.002$ (see Section D for a more detailed analysis). Therefore, even early on in the operation of SKA-Low, we will gain insight into the cosmological properties of reionization with minimal knowledge of the properties of ionizing sources (marginalizing over them). By fixing the k-scale and varying the redshift, the forecast is naturally sensitive to the variance of the 21 cm field. This is partly due to the choice of redshift,

as once significant overlap of H II bubbles has taken place, that is past the mid-point of the EoR, the 21 cm power spectrum decreases in amplitude.

For both cases presented in Fig. 5, we can place loose lower limits on $\log_{10}(\alpha_0)$, although the 1000 h case is more constraining on both distributions: Adding integration time additionally provides a lower limit on $\kappa$. Note that the ranges of the axes in the figure are smaller than those of the priors from Table 1 (while the limit on $\kappa$ is extended), confirming that these limits are not only informed by the priors. We further investigate whether the choice of k-scale has an effect on these results in Section 5.1.2.

To get more insight into the results above, we look in more detail at the dependence of our two observables on the model parameters. Fig. 6 presents the reionization history, the pkSZ power spectrum, and the 21 cm power spectrum at $z = 6.5$ and $7.8$, corresponding to the two redshifts considered above for the mock data points (in each column of the figure, respectively), for different values of the model parameters. All parameters are varied within the value ranges used as flat priors in the analysis. Regarding the global reionization parameters, we see that varying the mid- and endpoint leads to similar, although opposite, variations in both observables. Increasing $z_{re}$ extends the duration of the EoR, which leads to a boosted amplitude of the pkSZ power spectrum while having a smaller effect on the amplitude of the 21 cm power spectrum. Conversely, increasing $z_{end}$ decreases the duration of the EoR, suppressing the pkSZ amplitude. This explains the strong anticorrelation observed between the two parameters in Fig. 5. Note that for the model where $z_{end} = 7$, naturally, the 21 cm power spectrum is zero at $z = 6.5$.

On the other hand, $\log_{10}(\alpha_0)$ and $\kappa$ are also strongly correlated. The morphology parameters have no influence on the global reionization history but impact the shape of the electron power spectrum and, in turn, the shape and amplitude of the 21 cm and pkSZ power spectra. However, compared to the global history parameters, increasing $\log_{10}(\alpha_0)$ or $\kappa$, as seen in the third and fourth columns of Fig. 6, results in a similar impact on both observables. The $\log_{10}(\alpha_0)$ parameter is effectively a measurement of the large-scale amplitude of $P_{ee}$ at high redshift (before the reionization mid-point), while the $P_{ee}$ power breaks as $k^{-3}$ for scales $k > \kappa$, such that $\kappa$ can be interpreted as the minimal size of ionized regions during reionization. Hence, varying $\log_{10}(\alpha_0)$ directly impacts the amplitude of both observables, while varying $\kappa$ changes the shapes of both power spectra. The maximum of the pkSZ power and the shoulder in the dimensionless 21 cm power spectrum shift towards smaller scales as $\kappa$ increases. This dependence on the observables on these two parameters explains the degeneracy between the morphology parameters seen in Fig. 5. For the 1k2z case, where both data points are at a fixed k-scale, a lower value $\log_{10}(\alpha_0)$ can be compensated by a higher $\kappa$ value, within the measurement error of the data. However, this degeneracy can be reduced in the presence of 21 cm data points at different k-scales or by kSZ data points at different multipoles, since $\kappa$ impacts the shape of both power spectra. We explore the k-scale dependency of our chosen data in Section 5.1.2.

### 4.2.2 Constraining power of each data set

In this subsection, we examine each data set's ability to individually impose constraints on reionization. In Fig. 7, we present the joint posterior distribution of the sampled parameters for three different cases assuming an SKA-Low detection after integrating over 1000 h. We compare a case where the forecast is run only with the 21 cm data points (○21), a case with the pkSZ measurement only (○pkSZ),

<footnote n="11">¹¹The Gelman–Rubin parameters for morphology parameters of the 1000 h case, given in Table 3, are of order $R - 1 \approx 1$, so that the chains are not fully converged despite running for over a million iterations. However, the worst values are obtained for parameters we marginalize over, and we noticed that the posteriors of relevant parameters do not change significantly as iterations increase, confirming that the results still hold good qualitative significance.</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Posterior distributions on reionization mid- and endpoint (upper panel) and the morphology parameters $\log_{10}(\alpha_0)$ and $\kappa$ (bottom panel), considering a detection of the pkSZ power spectrum at $\ell = 3000$ and of the 21 cm power spectrum at redshifts $z = 6.5, 7.8$, and $k = 0.50 \text{ Mpc}^{-1}$ for either 100 or 1000 h of integration with SKA-Low, in purple and in blue, respectively. Note that the axis ranges in both panels are smaller than the priors in Table 1 or the ones presented in following figures.</figcaption></figure>

In Fig. 5, we show the joint probability distributions for our four sampled parameters (see appendix B for a discussion on parameter degeneracies). The best-fitting parameters, their $1\sigma$ error, and Gelman–Rubin convergence parameters are presented in Table 3.<fnref n="11" />

Our results indicate that we can recover the true values of $z_{re}$ and $z_{end}$ in both the 1k2z cases, despite the order of magnitude difference in integration time between them. Moreover, we independently constrain the value of the Thomson optical depth with a deviation of less than one per cent of the true value, $\tau = 0.066 \pm 0.002$ (see Section D for a more detailed analysis). Therefore, even early on in the operation of SKA-Low, we will gain insight into the cosmological properties of reionization with minimal knowledge of the properties of ionizing sources (marginalizing over them). By fixing the k-scale and varying the redshift, the forecast is naturally sensitive to the variance of the 21 cm field. This is partly due to the choice of redshift,

as once significant overlap of H II bubbles has taken place, that is past the mid-point of the EoR, the 21 cm power spectrum decreases in amplitude.

For both cases presented in Fig. 5, we can place loose lower limits on $\log_{10}(\alpha_0)$, although the 1000 h case is more constraining on both distributions: Adding integration time additionally provides a lower limit on $\kappa$. Note that the ranges of the axes in the figure are smaller than those of the priors from Table 1 (while the limit on $\kappa$ is extended), confirming that these limits are not only informed by the priors. We further investigate whether the choice of k-scale has an effect on these results in Section 5.1.2.

To get more insight into the results above, we look in more detail at the dependence of our two observables on the model parameters. Fig. 6 presents the reionization history, the pkSZ power spectrum, and the 21 cm power spectrum at $z = 6.5$ and $7.8$, corresponding to the two redshifts considered above for the mock data points (in each column of the figure, respectively), for different values of the model parameters. All parameters are varied within the value ranges used as flat priors in the analysis. Regarding the global reionization parameters, we see that varying the mid- and endpoint leads to similar, although opposite, variations in both observables. Increasing $z_{re}$ extends the duration of the EoR, which leads to a boosted amplitude of the pkSZ power spectrum while having a smaller effect on the amplitude of the 21 cm power spectrum. Conversely, increasing $z_{end}$ decreases the duration of the EoR, suppressing the pkSZ amplitude. This explains the strong anticorrelation observed between the two parameters in Fig. 5. Note that for the model where $z_{end} = 7$, naturally, the 21 cm power spectrum is zero at $z = 6.5$.

On the other hand, $\log_{10}(\alpha_0)$ and $\kappa$ are also strongly correlated. The morphology parameters have no influence on the global reionization history but impact the shape of the electron power spectrum and, in turn, the shape and amplitude of the 21 cm and pkSZ power spectra. However, compared to the global history parameters, increasing $\log_{10}(\alpha_0)$ or $\kappa$, as seen in the third and fourth columns of Fig. 6, results in a similar impact on both observables. The $\log_{10}(\alpha_0)$ parameter is effectively a measurement of the large-scale amplitude of $P_{ee}$ at high redshift (before the reionization mid-point), while the $P_{ee}$ power breaks as $k^{-3}$ for scales $k > \kappa$, such that $\kappa$ can be interpreted as the minimal size of ionized regions during reionization. Hence, varying $\log_{10}(\alpha_0)$ directly impacts the amplitude of both observables, while varying $\kappa$ changes the shapes of both power spectra. The maximum of the pkSZ power and the shoulder in the dimensionless 21 cm power spectrum shift towards smaller scales as $\kappa$ increases. This dependence on the observables on these two parameters explains the degeneracy between the morphology parameters seen in Fig. 5. For the 1k2z case, where both data points are at a fixed k-scale, a lower value $\log_{10}(\alpha_0)$ can be compensated by a higher $\kappa$ value, within the measurement error of the data. However, this degeneracy can be reduced in the presence of 21 cm data points at different k-scales or by kSZ data points at different multipoles, since $\kappa$ impacts the shape of both power spectra. We explore the k-scale dependency of our chosen data in Section 5.1.2.

### 4.2.2 Constraining power of each data set

In this subsection, we examine each data set's ability to individually impose constraints on reionization. In Fig. 7, we present the joint posterior distribution of the sampled parameters for three different cases assuming an SKA-Low detection after integrating over 1000 h. We compare a case where the forecast is run only with the 21 cm data points (○21), a case with the pkSZ measurement only (○pkSZ),

<footnote n="11">¹¹The Gelman–Rubin parameters for morphology parameters of the 1000 h case, given in Table 3, are of order $R - 1 \approx 1$, so that the chains are not fully converged despite running for over a million iterations. However, the worst values are obtained for parameters we marginalize over, and we noticed that the posteriors of relevant parameters do not change significantly as iterations increase, confirming that the results still hold good qualitative significance.</footnote>
<table><caption>**Table 3.** Best-fitting values of the model parameters for the 1k2z cases, their corresponding 1σ uncertainty and Gelman–Rubin parameter. Best-fitting values and errors are also given for derived parameters such as τ.</caption><thead><tr><th>Models</th><th></th><th></th><th>z<sub>re</sub></th><th>z<sub>end</sub></th><th>log<sub>10</sub>(α<sub>0</sub>)</th><th>κ</th><th>τ</th><th>dz</th></tr></thead><tbody><tr><td rowspan="2">Label</td><td>Data</td><td>True</td><td>7.37</td><td>6.15</td><td>3.12</td><td>0.16</td><td>0.0649</td><td>1.22</td></tr><tr><td>1k2z</td><td>z = 6.5, 7.8<br/>k = 0.5 Mpc<sup>-1</sup></td><td>1000 h<br/>R - 1</td><td>7.39 ± 0.14<br/>0.3</td><td>6.15 ± 0.04<br/>1.5</td><td>3.04 ± 0.32<br/>1.6</td><td>0.18 ± 0.04<br/>0.5</td><td>0.0651 ± 0.0018<br/>N/A</td><td>1.24 ± 0.17<br/>N/A</td></tr><tr><td rowspan="2">Label</td><td>Data</td><td>True</td><td>7.42 ± 0.13</td><td>6.15 ± 0.04</td><td>2.95 ± 0.31</td><td>0.19 ± 0.04</td><td>0.0655 ± 0.0017</td><td>1.27 ± 0.16</td></tr><tr><td>1k2z</td><td>z = 6.5, 7.8<br/>k = 0.5 Mpc<sup>-1</sup></td><td>100h<br/>R - 1</td><td>0.1</td><td>0.1</td><td>0.3</td><td>0.3</td><td>N/A</td><td>N/A</td></tr></tbody></table>

<table><caption>**Table 4.** Best-fitting values of the distributions in Fig. 9 and their corresponding 1σ uncertainty as well as the Gelman–Rubin convergence diagnostic for each parameter.</caption><thead><tr><th>Models</th><th></th><th></th><th>z<sub>re</sub></th><th>z<sub>end</sub></th><th>log<sub>10</sub>(α<sub>0</sub>)</th><th>κ</th><th>τ</th><th>dz</th></tr></thead><tbody><tr><td rowspan="2">Label</td><td>Data</td><td>True</td><td>7.37</td><td>6.15</td><td>3.12</td><td>0.16</td><td>0.0649</td><td>1.22</td></tr><tr><td>MWA</td><td>z = 6.5, 7.8<br/>k = 0.1 Mpc<sup>-1</sup></td><td>1000 h<br/>R - 1</td><td>7.02 ± 0.49<br/>0.01</td><td>5.39 ± 0.54<br/>0.008</td><td>3.25 ± 0.45<br/>0.03</td><td>0.13 ± 0.04<br/>0.03</td><td>0.0623 ± 0.0051<br/>N/A</td><td>1.53 ± 0.49<br/>N/A</td></tr><tr><td rowspan="2">LOFAR</td><td>Data</td><td>True</td><td>7.30 ± 0.87</td><td>5.74 ± 0.80</td><td>3.15 ± 0.33</td><td>0.14 ± 0.03</td><td>0.0643 ± 0.0094</td><td>1.38 ± 0.83</td></tr><tr><td>LOFAR</td><td>z = 8.3, 9.1<br/>k = 0.1 Mpc<sup>-1</sup></td><td>1000 h<br/>R - 1</td><td>0.06</td><td>0.03</td><td>0.1</td><td>0.04</td><td>N/A</td><td>N/A</td></tr></tbody></table>

<table><caption>**Table 5.** Best-fitting values of the distributions in and their corresponding 1σ uncertainty as well as the Gelman–Rubin convergence diagnostic for each parameter.</caption><thead><tr><th>Models</th><th></th><th></th><th>z<sub>re</sub></th><th>z<sub>end</sub></th><th>log<sub>10</sub>(α<sub>0</sub>)</th><th>κ</th><th>τ</th><th>dz</th></tr></thead><tbody><tr><td rowspan="2">Label</td><td>Data</td><td>True</td><td>7.37</td><td>6.15</td><td>3.12</td><td>0.16</td><td>0.0649</td><td>1.22</td></tr><tr><td>2k1z</td><td>z = 6.5, 6.5<br/>k = 0.1, 0.5 Mpc<sup>-1</sup></td><td>1000 h<br/>R - 1</td><td>7.22 ± 0.21<br/>0.3</td><td>5.29 ± 0.55<br/>0.9</td><td>2.97 ± 0.08<br/>4.0</td><td>0.15 ± 0.01<br/>1.2</td><td>0.0649 ± 0.0016<br/>N/A</td><td>1.90 ± 0.40<br/>N/A</td></tr><tr><td rowspan="2">Label</td><td>Data</td><td>True</td><td>7.28 ± 0.33</td><td>6.02 ± 0.52</td><td>3.09 ± 0.15</td><td>0.17 ± 0.02</td><td>0.0642 ± 0.0025</td><td>1.27 ± 0.19</td></tr><tr><td>hiz</td><td>z = 7.8, 10.4<br/>k = 0.1, 0.1 Mpc<sup>-1</sup></td><td>1000 h<br/>R - 1</td><td>0.7</td><td>0.7</td><td>0.5</td><td>0.5</td><td>N/A</td><td>N/A</td></tr><tr><td rowspan="2">Label</td><td>Data</td><td>True</td><td>7.37 ± 0.07</td><td>6.15 ± 0.03</td><td>3.12 ± 0.01</td><td>0.164 ± 0.001</td><td>0.0649 ± 0.0010</td><td>1.22 ± 0.10</td></tr><tr><td>3k2z</td><td>z = 6.5, 6.5, 7.8<br/>k = 0.1, 0.5, 0.1 Mpc<sup>-1</sup></td><td>1000 h<br/>R - 1</td><td>0.02</td><td>0.02</td><td>0.005</td><td>0.004</td><td>N/A</td><td>N/A</td></tr></tbody></table>

and, lastly, the joint forecast previously illustrated in Fig. 5. To better understand the role of each parameter in each of the aforementioned cases, we firstly group our analysis by examining the global reionization parameters $z_{re}$ and $z_{end}$, then the astrophysical $log_{10}(\alpha_0)$ and $\kappa$, in Fig. 7.

Examining the two-dimensional $z_{re}$ and $z_{end}$ posterior distributions for the opkSZ case, we find that the parameters are positively correlated. As discussed in Section 4.2, this occurs because the pkSZ effect is sensitive to the duration of reionization $dz \equiv z_{re} - z_{end}$ as the rate of the growth of ionizing bubbles impacts the amplitude of the signal (see also the first and second column of Fig. 6). Models with an extended reionization period result in a larger amplitude of the pkSZ power spectrum as more free electrons interact with the ionized medium, and reciprocally (McQuinn et al. 2005; Mesinger, McQuinn & Spergel 2012; Battaglia et al. 2013). Note that the posterior distribution for the opkSZ case covers a region of parameter space where the universe has reionized earlier than the redshift at which our 21 cm data points are measured. These models are naturally excluded when combining both data sets. For the 021 case, the data constraints tightly constrain the end of reionization but only place an upper limit on the mid-point. Compared to the pkSZ-only case, the 021 case is not sensitive to the duration of reionization and favours extended models of the EoR, as previously noted by Bégin, Liu & Gorce (2022). Conversely, a detection of the 21 cm power spectrum at $z = 6.5$ naturally excludes models where reionization finishes earlier than that redshift, hence, the forecast gives a lower limit on $z_{end}$. Lastly, models at the lower left of the

two-dimensional posterior ($z_{re} \le 7$ and $z_{end} \le 6$) are excluded by the 21 cm data at $z = 6.5$ and 7.7. For these models, the 21 cm power spectrum at these redshifts would trace the density power spectrum. Such a feature is inconsistent with the amplitude of the 21 cm data points and their associated error bars. Therefore, the complementary nature of both probes breaks the degeneracy inherent to each data set and allows us to constrain the parameter space and recover the true reionization history.

Regarding $log_{10}(\alpha_0)$ and $\kappa$, we see that both morphology parameters are anticorrelated and the pkSZ data alone only provides a biased measurement. This can be better understood by examining the two lower rows of Fig. 6: Increasing either of the morphology parameters results in a boosting of the pkSZ and 21 cm power spectra, resulting in a degeneracy between them. The pkSZ data seems to favour a low $\kappa$, which in order to match the measurement is then compensated by a low $log_{10}(\alpha_0)$. Consequently, while the true pkSZ power spectra peaks at $\ell \approx 3000$, here there are multiple models for which the pkSZ reaches its maximum amplitude at $\ell \approx 2000$. Part of this bias stems from the simplifications in equation (7), which does not fully account for the cross-correlation between the ionization and density fields (see Section 4.1) and results in a boosted amplitude of the 21 cm power spectrum. As the opkSZ case does not contain information from 21 cm data points, these models are not constrained, resulting in a biased distribution. We confirm that this effect does not occur in the 021 case, where we only use the 21 cm data points. However, as both data points are for the same $k$ value, the parameter values are not well constrained (see Section 5.1.2).
<table><caption>**Table 3.** Best-fitting values of the model parameters for the 1k2z cases, their corresponding 1σ uncertainty and Gelman–Rubin parameter. Best-fitting values and errors are also given for derived parameters such as τ.</caption><thead><tr><th>Models</th><th></th><th></th><th>z<sub>re</sub></th><th>z<sub>end</sub></th><th>log<sub>10</sub>(α<sub>0</sub>)</th><th>κ</th><th>τ</th><th>dz</th></tr></thead><tbody><tr><td rowspan="2">Label</td><td>Data</td><td>True</td><td>7.37</td><td>6.15</td><td>3.12</td><td>0.16</td><td>0.0649</td><td>1.22</td></tr><tr><td>1k2z</td><td>z = 6.5, 7.8<br/>k = 0.5 Mpc<sup>-1</sup></td><td>1000 h<br/>R - 1</td><td>7.39 ± 0.14<br/>0.3</td><td>6.15 ± 0.04<br/>1.5</td><td>3.04 ± 0.32<br/>1.6</td><td>0.18 ± 0.04<br/>0.5</td><td>0.0651 ± 0.0018<br/>N/A</td><td>1.24 ± 0.17<br/>N/A</td></tr><tr><td rowspan="2">Label</td><td>Data</td><td>True</td><td>7.42 ± 0.13</td><td>6.15 ± 0.04</td><td>2.95 ± 0.31</td><td>0.19 ± 0.04</td><td>0.0655 ± 0.0017</td><td>1.27 ± 0.16</td></tr><tr><td>1k2z</td><td>z = 6.5, 7.8<br/>k = 0.5 Mpc<sup>-1</sup></td><td>100h<br/>R - 1</td><td>0.1</td><td>0.1</td><td>0.3</td><td>0.3</td><td>N/A</td><td>N/A</td></tr></tbody></table>

<table><caption>**Table 4.** Best-fitting values of the distributions in Fig. 9 and their corresponding 1σ uncertainty as well as the Gelman–Rubin convergence diagnostic for each parameter.</caption><thead><tr><th>Models</th><th></th><th></th><th>z<sub>re</sub></th><th>z<sub>end</sub></th><th>log<sub>10</sub>(α<sub>0</sub>)</th><th>κ</th><th>τ</th><th>dz</th></tr></thead><tbody><tr><td rowspan="2">Label</td><td>Data</td><td>True</td><td>7.37</td><td>6.15</td><td>3.12</td><td>0.16</td><td>0.0649</td><td>1.22</td></tr><tr><td>MWA</td><td>z = 6.5, 7.8<br/>k = 0.1 Mpc<sup>-1</sup></td><td>1000 h<br/>R - 1</td><td>7.02 ± 0.49<br/>0.01</td><td>5.39 ± 0.54<br/>0.008</td><td>3.25 ± 0.45<br/>0.03</td><td>0.13 ± 0.04<br/>0.03</td><td>0.0623 ± 0.0051<br/>N/A</td><td>1.53 ± 0.49<br/>N/A</td></tr><tr><td rowspan="2">LOFAR</td><td>Data</td><td>True</td><td>7.30 ± 0.87</td><td>5.74 ± 0.80</td><td>3.15 ± 0.33</td><td>0.14 ± 0.03</td><td>0.0643 ± 0.0094</td><td>1.38 ± 0.83</td></tr><tr><td>LOFAR</td><td>z = 8.3, 9.1<br/>k = 0.1 Mpc<sup>-1</sup></td><td>1000 h<br/>R - 1</td><td>0.06</td><td>0.03</td><td>0.1</td><td>0.04</td><td>N/A</td><td>N/A</td></tr></tbody></table>

<table><caption>**Table 5.** Best-fitting values of the distributions in and their corresponding 1σ uncertainty as well as the Gelman–Rubin convergence diagnostic for each parameter.</caption><thead><tr><th>Models</th><th></th><th></th><th>z<sub>re</sub></th><th>z<sub>end</sub></th><th>log<sub>10</sub>(α<sub>0</sub>)</th><th>κ</th><th>τ</th><th>dz</th></tr></thead><tbody><tr><td rowspan="2">Label</td><td>Data</td><td>True</td><td>7.37</td><td>6.15</td><td>3.12</td><td>0.16</td><td>0.0649</td><td>1.22</td></tr><tr><td>2k1z</td><td>z = 6.5, 6.5<br/>k = 0.1, 0.5 Mpc<sup>-1</sup></td><td>1000 h<br/>R - 1</td><td>7.22 ± 0.21<br/>0.3</td><td>5.29 ± 0.55<br/>0.9</td><td>2.97 ± 0.08<br/>4.0</td><td>0.15 ± 0.01<br/>1.2</td><td>0.0649 ± 0.0016<br/>N/A</td><td>1.90 ± 0.40<br/>N/A</td></tr><tr><td rowspan="2">Label</td><td>Data</td><td>True</td><td>7.28 ± 0.33</td><td>6.02 ± 0.52</td><td>3.09 ± 0.15</td><td>0.17 ± 0.02</td><td>0.0642 ± 0.0025</td><td>1.27 ± 0.19</td></tr><tr><td>hiz</td><td>z = 7.8, 10.4<br/>k = 0.1, 0.1 Mpc<sup>-1</sup></td><td>1000 h<br/>R - 1</td><td>0.7</td><td>0.7</td><td>0.5</td><td>0.5</td><td>N/A</td><td>N/A</td></tr><tr><td rowspan="2">Label</td><td>Data</td><td>True</td><td>7.37 ± 0.07</td><td>6.15 ± 0.03</td><td>3.12 ± 0.01</td><td>0.164 ± 0.001</td><td>0.0649 ± 0.0010</td><td>1.22 ± 0.10</td></tr><tr><td>3k2z</td><td>z = 6.5, 6.5, 7.8<br/>k = 0.1, 0.5, 0.1 Mpc<sup>-1</sup></td><td>1000 h<br/>R - 1</td><td>0.02</td><td>0.02</td><td>0.005</td><td>0.004</td><td>N/A</td><td>N/A</td></tr></tbody></table>

and, lastly, the joint forecast previously illustrated in Fig. 5. To better understand the role of each parameter in each of the aforementioned cases, we firstly group our analysis by examining the global reionization parameters $z_{re}$ and $z_{end}$, then the astrophysical $log_{10}(\alpha_0)$ and $\kappa$, in Fig. 7.

Examining the two-dimensional $z_{re}$ and $z_{end}$ posterior distributions for the opkSZ case, we find that the parameters are positively correlated. As discussed in Section 4.2, this occurs because the pkSZ effect is sensitive to the duration of reionization $dz \equiv z_{re} - z_{end}$ as the rate of the growth of ionizing bubbles impacts the amplitude of the signal (see also the first and second column of Fig. 6). Models with an extended reionization period result in a larger amplitude of the pkSZ power spectrum as more free electrons interact with the ionized medium, and reciprocally (McQuinn et al. 2005; Mesinger, McQuinn & Spergel 2012; Battaglia et al. 2013). Note that the posterior distribution for the opkSZ case covers a region of parameter space where the universe has reionized earlier than the redshift at which our 21 cm data points are measured. These models are naturally excluded when combining both data sets. For the 021 case, the data constraints tightly constrain the end of reionization but only place an upper limit on the mid-point. Compared to the pkSZ-only case, the 021 case is not sensitive to the duration of reionization and favours extended models of the EoR, as previously noted by Bégin, Liu & Gorce (2022). Conversely, a detection of the 21 cm power spectrum at $z = 6.5$ naturally excludes models where reionization finishes earlier than that redshift, hence, the forecast gives a lower limit on $z_{end}$. Lastly, models at the lower left of the

two-dimensional posterior ($z_{re} \le 7$ and $z_{end} \le 6$) are excluded by the 21 cm data at $z = 6.5$ and 7.7. For these models, the 21 cm power spectrum at these redshifts would trace the density power spectrum. Such a feature is inconsistent with the amplitude of the 21 cm data points and their associated error bars. Therefore, the complementary nature of both probes breaks the degeneracy inherent to each data set and allows us to constrain the parameter space and recover the true reionization history.

Regarding $log_{10}(\alpha_0)$ and $\kappa$, we see that both morphology parameters are anticorrelated and the pkSZ data alone only provides a biased measurement. This can be better understood by examining the two lower rows of Fig. 6: Increasing either of the morphology parameters results in a boosting of the pkSZ and 21 cm power spectra, resulting in a degeneracy between them. The pkSZ data seems to favour a low $\kappa$, which in order to match the measurement is then compensated by a low $log_{10}(\alpha_0)$. Consequently, while the true pkSZ power spectra peaks at $\ell \approx 3000$, here there are multiple models for which the pkSZ reaches its maximum amplitude at $\ell \approx 2000$. Part of this bias stems from the simplifications in equation (7), which does not fully account for the cross-correlation between the ionization and density fields (see Section 4.1) and results in a boosted amplitude of the 21 cm power spectrum. As the opkSZ case does not contain information from 21 cm data points, these models are not constrained, resulting in a biased distribution. We confirm that this effect does not occur in the 021 case, where we only use the 21 cm data points. However, as both data points are for the same $k$ value, the parameter values are not well constrained (see Section 5.1.2).
<figure><img src="image_6.png" /><figcaption>**Figure 6.** Evolution of the reionization history (first column), the pkSZ angular power spectrum (second column), and the dimensionless 21 cm power spectrum at redshifts $z = 6.5, 7.8$, considered in Section 4.2 (last two columns, respectively) as a function of our four model parameters. The values of $\kappa$ are shown as vertical dotted lines in the final row with the corresponding colour. Note the difference in the shape of the 21 cm power spectrum at high-$k$ compared to Fig. 1 is discussed in Section 4.1.</figcaption></figure>

To summarize, we find that with as few as 100 h of integration with SKA-Low, we can successfully constrain the global reionization parameters $z_{re}$ and $z_{end}$, the data at different redshifts giving access to the evolution of the 21 cm signal. By conducting our analysis with either only the pkSZ or only the 21 cm data, we find that, while the 21 cm signal can limit the range of values of $z_{end}$, it cannot constrain $z_{re}$. Conversely, we find that the pkSZ measurement mainly constrains the duration of reionization ($dz = z_{re} - z_{end}$). Naturally, the combined analysis benefits from the complementarity of probes. While 100 h of integration are sufficient to obtain upper limits on the reionization morphology parameters, it is only with 1000 h that we are capable of constraining them. The morphology parameters have a direct impact on the shape of the 21 cm power spectrum across different $k$-scales and we discuss how we can further constrain them in the following section.

# 5 DISCUSSION AND LIMITATIONS

No estimator is without its limitations and it is essential to consider its biases to better grasp the scope of the results. In this section, we delve into a series of cases designed to examine the capabilities of the forecast. We vary the redshifts and $k$-ranges at which the data points were selected in Section 4.2.1, as well as the noise models and study the impact on our results. We also explore how the choice of

input parameters affects the forecast. Note that for all cases in this discussion, we limit ourselves to the fixed number of data points as in Section 4.2.1 to quantify the role of each data point and its driving physical processes. The interested reader can refer to Section 5.3 for results on how an additional 21 cm data point can significantly improve the forecast.

## 5.1 Importance of the choice of data points

### 5.1.1 Choice of redshift

We explore the constraining ability of the forecast based on the choice of redshifts for the 21 cm observation. To do so, we examine the SKA-Low 1000h case (hiz) and vary the redshift of the observed 21 cm power spectrum from $z = 6.5$ and $7.8$ to $z = 7.8$ and $10.4$, motivated by the recent HERA upper limits (The HERA Collaboration 2022). The resulting two-dimensional posterior distributions for $z_{re}-z_{end}$ and $\log \alpha_0 - \kappa$ are presented in Fig. 8 in purple. Compared to the fiducial case (in blue), two main differences are apparent. First, the data are now compatible with models whose mid- and endpoint are generally at higher redshifts (early and rapid reionization scenarios). For comparison, for the 1k2z cases in Section 4.2.1, the fact that the 21 cm signal was non-zero at $z = 6.5$ naturally excludes models where reionization was completed by then. Additionally, the high-redshift
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Evolution of the reionization history (first column), the pkSZ angular power spectrum (second column), and the dimensionless 21 cm power spectrum at redshifts $z = 6.5, 7.8$, considered in Section 4.2 (last two columns, respectively) as a function of our four model parameters. The values of $\kappa$ are shown as vertical dotted lines in the final row with the corresponding colour. Note the difference in the shape of the 21 cm power spectrum at high-$k$ compared to Fig. 1 is discussed in Section 4.1.</figcaption></figure>

To summarize, we find that with as few as 100 h of integration with SKA-Low, we can successfully constrain the global reionization parameters $z_{re}$ and $z_{end}$, the data at different redshifts giving access to the evolution of the 21 cm signal. By conducting our analysis with either only the pkSZ or only the 21 cm data, we find that, while the 21 cm signal can limit the range of values of $z_{end}$, it cannot constrain $z_{re}$. Conversely, we find that the pkSZ measurement mainly constrains the duration of reionization ($dz = z_{re} - z_{end}$). Naturally, the combined analysis benefits from the complementarity of probes. While 100 h of integration are sufficient to obtain upper limits on the reionization morphology parameters, it is only with 1000 h that we are capable of constraining them. The morphology parameters have a direct impact on the shape of the 21 cm power spectrum across different $k$-scales and we discuss how we can further constrain them in the following section.

# 5 DISCUSSION AND LIMITATIONS

No estimator is without its limitations and it is essential to consider its biases to better grasp the scope of the results. In this section, we delve into a series of cases designed to examine the capabilities of the forecast. We vary the redshifts and $k$-ranges at which the data points were selected in Section 4.2.1, as well as the noise models and study the impact on our results. We also explore how the choice of

input parameters affects the forecast. Note that for all cases in this discussion, we limit ourselves to the fixed number of data points as in Section 4.2.1 to quantify the role of each data point and its driving physical processes. The interested reader can refer to Section 5.3 for results on how an additional 21 cm data point can significantly improve the forecast.

## 5.1 Importance of the choice of data points

### 5.1.1 Choice of redshift

We explore the constraining ability of the forecast based on the choice of redshifts for the 21 cm observation. To do so, we examine the SKA-Low 1000h case (hiz) and vary the redshift of the observed 21 cm power spectrum from $z = 6.5$ and $7.8$ to $z = 7.8$ and $10.4$, motivated by the recent HERA upper limits (The HERA Collaboration 2022). The resulting two-dimensional posterior distributions for $z_{re}-z_{end}$ and $\log \alpha_0 - \kappa$ are presented in Fig. 8 in purple. Compared to the fiducial case (in blue), two main differences are apparent. First, the data are now compatible with models whose mid- and endpoint are generally at higher redshifts (early and rapid reionization scenarios). For comparison, for the 1k2z cases in Section 4.2.1, the fact that the 21 cm signal was non-zero at $z = 6.5$ naturally excludes models where reionization was completed by then. Additionally, the high-redshift
<figure><img src="image_7.png" /><figcaption>**Figure 7.** Posterior distributions on reionization mid- and endpoint (upper panel) and the morphology parameters $\log_{10}(\alpha_0)$ and $\kappa$ (bottom panel) when fitting 21 cm and pkSZ power spectra individually (green and magenta) or jointly (in blue), for 1000 h of SKA-Low observations. The vertical and horizontal lines show the ‘true’ values used to generate the mock data.</figcaption></figure>

data allow for later reionization histories than the truth (lower left quadrant of the figure). Indeed, the data point at $z = 10.4$ has limited constraining power on the global history parameters but provides better constraints on morphology parameters, as illustrated in the bottom panel of Fig. 8, in orange. The reason for this is that at that redshift, $x_v \ll 1$, such that $P_{21} \sim P_{bb} + x_m^2 P_{ee}$ and we do not vary the baryon density power spectrum in our analysis. Meanwhile, the shape of the electron density power spectrum is governed by the power law in equation (12).

### 5.1.2 Choice of k-scales

Next, we investigate the constraining power of the forecast using two detections of the 21 cm power spectrum at $z = 6.5$ and $k = 0.1$ and $0.5 \text{ Mpc}^{-1}$. We assume thermal noise errors corresponding to 1000 h of observations with SKA-Low. As before, a ten per cent error bar is assumed on the pkSZ power spectrum. We refer to this case as 2k1z. Previously, we utilized the fact that reionization is an evolving process, while here, we make use of its multiscale nature. Cosmic reionization is thought to be an inside-out phase transition, where the densest parts of the Universe are the first to be re-ionized and cosmic voids are the last (Iliev et al. 2014). This will impact the shape of the 21 cm power spectrum, specifically, the shape of the high-k ‘tail’ that encodes information about the IGM on small scales.

<figure><img src="image_8.png" /><figcaption>**Figure 8.** Posterior distributions on reionization mid- and endpoint (upper panel) and the morphology parameters $\log_{10}(\alpha_0)$ and $\kappa$ (bottom panel) using 21 cm power spectra measurements from redshifts 7.8 and 10.4 at $k = 0.5 \text{ Mpc}^{-1}$ and kSZ power spectrum measurement at $\ell = 3000$ for 1000 h of SKA-Low observations compared to the fiducial case from Section 4.2 (in blue). The vertical and horizontal lines show the ‘true’ values used to generate the mock data.</figcaption></figure>

Fig. 8 illustrates the constraining power of this approach and results for 2k1z are shown in yellow, in contrast to the fiducial 1k2z case, shown in blue. While the data have clear constraining power on the $z_{re}$ and $z_{end}$ values, the two-dimensional posterior distribution exhibits a tail and is biased towards late reionization models ($z_{end} = 5.29 \pm 0.55$, while the true value is at redshift 6.15). Observing a non-zero 21 cm data point close to the true $z_{end}$ naturally excludes early reionization models. Despite the increased measurement error, the low-k (large-scale) amplitude of the 21 cm power spectrum constrains the upper limit on the midpoint. Indeed, models with an early mid-point of reionization would either re-ionize earlier or have a lower large-scale power, which is inconsistent with the measurement. Additionally, as can be inferred from Fig. 1, the low-k 21 cm measurement has a higher uncertainty, primarily due to a larger sample variance. This results in a poorly constrained tail of the 21 cm power spectrum, which undermines our ability to probe the state of the IGM with a single redshift measurement. Compared to the fiducial 1k2z case, the choice of probing two k-scales results in a tighter constraining power on $\log_{10}(\alpha_0)$ and $\kappa$, although the two-dimensional posterior is biased towards lower values of the parameters, corresponding to later reionization models. Examining and comparing the 21 cm power spectra in Fig. 6 reveals that the differences in k-scale are most prominent for $k > 1.0 \text{ Mpc}^{-1}$ scales.
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Posterior distributions on reionization mid- and endpoint (upper panel) and the morphology parameters $\log_{10}(\alpha_0)$ and $\kappa$ (bottom panel) when fitting 21 cm and pkSZ power spectra individually (green and magenta) or jointly (in blue), for 1000 h of SKA-Low observations. The vertical and horizontal lines show the ‘true’ values used to generate the mock data.</figcaption></figure>

data allow for later reionization histories than the truth (lower left quadrant of the figure). Indeed, the data point at $z = 10.4$ has limited constraining power on the global history parameters but provides better constraints on morphology parameters, as illustrated in the bottom panel of Fig. 8, in orange. The reason for this is that at that redshift, $x_v \ll 1$, such that $P_{21} \sim P_{bb} + x_m^2 P_{ee}$ and we do not vary the baryon density power spectrum in our analysis. Meanwhile, the shape of the electron density power spectrum is governed by the power law in equation (12).

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Posterior distributions on reionization mid- and endpoint (upper panel) and the morphology parameters $\log_{10}(\alpha_0)$ and $\kappa$ (bottom panel) using 21 cm power spectra measurements from redshifts 7.8 and 10.4 at $k = 0.5$ Mpc$^{-1}$ and kSZ power spectrum measurement at $\ell = 3000$ for 1000 h of SKA-Low observations compared to the fiducial case from Section 4.2 (in blue). The vertical and horizontal lines show the ‘true’ values used to generate the mock data.</figcaption></figure>

### 5.1.2 Choice of k-scales

Next, we investigate the constraining power of the forecast using two detections of the 21 cm power spectrum at $z = 6.5$ and $k = 0.1$ and $0.5$ Mpc$^{-1}$. We assume thermal noise errors corresponding to 1000 h of observations with SKA-Low. As before, a ten per cent error bar is assumed on the pkSZ power spectrum. We refer to this case as 2k1z. Previously, we utilized the fact that reionization is an evolving process, while here, we make use of its multiscale nature. Cosmic reionization is thought to be an inside-out phase transition, where the densest parts of the Universe are the first to be re-ionized and cosmic voids are the last (Iliev et al. 2014). This will impact the shape of the 21 cm power spectrum, specifically, the shape of the high-$k$ ‘tail’ that encodes information about the IGM on small scales.

Fig. 8 illustrates the constraining power of this approach and results for 2k1z are shown in yellow, in contrast to the fiducial 1k2z case, shown in blue. While the data have clear constraining power on the $z_{re}$ and $z_{end}$ values, the two-dimensional posterior distribution exhibits a tail and is biased towards late reionization models ($z_{end} = 5.29 \pm 0.55$, while the true value is at redshift 6.15). Observing a non-zero 21 cm data point close to the true $z_{end}$ naturally excludes early reionization models. Despite the increased measurement error, the low-$k$ (large-scale) amplitude of the 21 cm power spectrum constrains the upper limit on the midpoint. Indeed, models with an early mid-point of reionization would either re-ionize earlier or have a lower large-scale power, which is inconsistent with the measurement. Additionally, as can be inferred from Fig. 1, the low-$k$ 21 cm measurement has a higher uncertainty, primarily due to a larger sample variance. This results in a poorly constrained tail of the 21 cm power spectrum, which undermines our ability to probe the state of the IGM with a single redshift measurement. Compared to the fiducial 1k2z case, the choice of probing two $k$-scales results in a tighter constraining power on $\log_{10}(\alpha_0)$ and $\kappa$, although the two-dimensional posterior is biased towards lower values of the parameters, corresponding to later reionization models. Examining and comparing the 21 cm power spectra in Fig. 6 reveals that the differences in $k$-scale are most prominent for $k > 1.0$ Mpc$^{-1}$ scales.
Since the high *k*-scale measurements of the 21 cm power spectrum exhibit a higher level of uncertainty due to larger instrumental noise (see Fig. 1), models are less tightly constrained on such scales.

In summary, with two 21 cm data points at different *k*-scales, the morphology parameters are better constrained as the scale-evolution of the 21 cm power spectrum is included in the forecast. This result is, however, mitigated by the fact that, by fixing the redshift, the uncertainty on the middle and end of reionization is increased, resulting in a tail of the $z_{re}-z_{end}$ posterior as well as biased morphology parameters. One way to improve constraints on $z_{end}$ would be to consider 21 cm data for $k > 1.0 \text{ Mpc}^{-1}$ scales, or to limit the prior on $z_{end}$ to values allowed by measurements of Ly $\alpha$ absorption in quasar spectra (Bosman et al. 2022).

## 5.2 Current telescope capabilities

Currently operating radio telescopes are getting closer and closer to a detection of the 21 cm power spectrum during the EoR. In this section, we examine the constraining ability of the forecast, assuming a potential measurement of the 21 cm power spectrum by the MWA and LOFAR radio telescopes, i.e. with higher noise levels than previously considered.<fnref n="12" /> Because of its specific characteristics, each telescope is sensitive to different redshifts and scales and, in turn, to different reionization parameters. For this reason, we choose LOFAR data points corresponding to a fixed *k*-scale of $k = 0.1 \text{ Mpc}^{-1}$ for redshifts $z = 8.2$ and 9.1 (based on Patil et al. 2017; Mertens et al. 2020). The MWA data is consistent with the data points chosen in Section 4.2.1. For both cases, we assume an integration time of 1000 h. The resulting two-dimensional posterior probability distributions of the cosmological and morphology parameters are presented in Fig. 9 and Table 4.

Results for the MWA case are overall consistent with Section 4.2.1, with key differences emerging due to the higher uncertainty on the 21 cm power spectrum, which is approximately 1000 times larger than that of the SKA (see Fig. 1). As shown in orange in the upper panel of Fig. 9, we are still able to constrain $z_{re} = 7.02 \pm 0.49$. However, the spread in the distribution has increased compared to the SKA-Low case in blue (for which $z_{re} = 7.39 \pm 0.14$), and the result is biased to low values (see Table 1). On the other hand, the constraint on the end of reionization is now limited to an upper limit $z_{end} < 7.5$. A lower constraining power is also observed for the morphology parameters. We can primarily obtain upper limits for $\log_{10}(\alpha_0) < 3.75$ and $\kappa > 0.08$, showcasing the ability of the mock MWA data to exclude outlying models of the EoR. The decrease in the overall constraining power of the forecast is linked to the increase in the uncertainty of the 21 cm power spectrum points. This increase results in favouring late-time reionization models where the 21 cm power spectrum at $z = 6.5$ has a higher amplitude and a lower value of $\kappa$, and the pkSZ power spectrum peaks at $\ell \approx 2000$, similar to that discussed in Section 4.2.2 for the opKSZ case. Such models could potentially be excluded with the addition of a pkSZ data point at $\ell = 2000$ or more data on the 21 cm power spectrum at lower *k* values.

Regarding the LOFAR case, results are roughly similar to what was observed in Section 5.1.2 where the 21 cm mock data is taken at a fixed *k*-scale but different redshifts. Naturally, increasing the uncertainty on the measurements by approximately two orders of magnitude compared to that of the SKA (see Fig. 1), the constraining power of the data is decreased. The choice of a higher *z* mock 21 cm

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Posterior distributions on reionization mid- and endpoint (upper panel) and the morphology parameters $\log_{10}(\alpha_0)$ and $\kappa$ (bottom panel) when fitting 21 cm and pkSZ power spectra jointly, for 1000 h of MWA, LOFAR, & SKA observations, seen in grey, orange, and blue. The vertical and horizontal lines show the ‘true’ values used to generate the mock data.</figcaption></figure>

data point further reduces the ability of the forecast to place upper limits on the global reionization parameters (see Section 5.1). For example, as seen in grey in the upper panel of Fig. 9, we can constrain $z_{re} < 8.0$ and can only exclude only late reionization models for which $z_{end} > 5.5$. Compared to the MWA case, having higher $z$ 21 cm data leads to higher constraining power on the morphology parameters. While we can only place an upper limit on $\kappa > 0.1$, the constraint on $\log_{10}(\alpha_0) = 3.15 \pm 0.33$ is fairly comparable to that of 1k2z 1000 h (where $\log_{10}(\alpha_0) = 3.04 \pm 0.32$). While this is mostly related to the choice of data, for example, the hiz case can constrain $\log_{10}(\alpha_0) = 3.09 \pm 0.15$, our results indicate the potential of high-redshift LOFAR measurements on constraining on the morphology parameters of the EoR.

A possible improvement is the addition of an informed prior on the Thomson optical depth from the Planck Collaboration I (2020) measurement, which could potentially aid in increasing the precision of the forecast and constraining the parameter space of the global history parameters.

## 5.3 Increasing the number of data points

In previous sections, we used a minimalist approach by limiting the number of detected modes to test how much information would be obtained from early, partial 21 cm data, as a detection will only be achieved incrementally. On the one hand, with two data points at

<footnote n="12">¹²Note that we do not consider HERA in this section, since its ideal noise levels are lower than the ones for the SKA (Fig. 1).</footnote>
Since the high *k*-scale measurements of the 21 cm power spectrum exhibit a higher level of uncertainty due to larger instrumental noise (see Fig. 1), models are less tightly constrained on such scales.

In summary, with two 21 cm data points at different *k*-scales, the morphology parameters are better constrained as the scale-evolution of the 21 cm power spectrum is included in the forecast. This result is, however, mitigated by the fact that, by fixing the redshift, the uncertainty on the middle and end of reionization is increased, resulting in a tail of the $z_{re}-z_{end}$ posterior as well as biased morphology parameters. One way to improve constraints on $z_{end}$ would be to consider 21 cm data for $k > 1.0 \ \text{Mpc}^{-1}$ scales, or to limit the prior on $z_{end}$ to values allowed by measurements of Ly $\alpha$ absorption in quasar spectra (Bosman et al. 2022).

## 5.2 Current telescope capabilities

Currently operating radio telescopes are getting closer and closer to a detection of the 21 cm power spectrum during the EoR. In this section, we examine the constraining ability of the forecast, assuming a potential measurement of the 21 cm power spectrum by the MWA and LOFAR radio telescopes, i.e. with higher noise levels than previously considered.<fnref n="12" /> Because of its specific characteristics, each telescope is sensitive to different redshifts and scales and, in turn, to different reionization parameters. For this reason, we choose LOFAR data points corresponding to a fixed *k*-scale of $k = 0.1 \ \text{Mpc}^{-1}$ for redshifts $z = 8.2$ and 9.1 (based on Patil et al. 2017; Mertens et al. 2020). The MWA data is consistent with the data points chosen in Section 4.2.1. For both cases, we assume an integration time of 1000 h. The resulting two-dimensional posterior probability distributions of the cosmological and morphology parameters are presented in Fig. 9 and Table 4.

Results for the MWA case are overall consistent with Section 4.2.1, with key differences emerging due to the higher uncertainty on the 21 cm power spectrum, which is approximately 1000 times larger than that of the SKA (see Fig. 1). As shown in orange in the upper panel of Fig. 9, we are still able to constrain $z_{re} = 7.02 \pm 0.49$. However, the spread in the distribution has increased compared to the SKA-Low case in blue (for which $z_{re} = 7.39 \pm 0.14$), and the result is biased to low values (see Table 1). On the other hand, the constraint on the end of reionization is now limited to an upper limit $z_{end} < 7.5$. A lower constraining power is also observed for the morphology parameters. We can primarily obtain upper limits for $\log_{10}(\alpha_0) < 3.75$ and $\kappa > 0.08$, showcasing the ability of the mock MWA data to exclude outlying models of the EoR. The decrease in the overall constraining power of the forecast is linked to the increase in the uncertainty of the 21 cm power spectrum points. This increase results in favouring late-time reionization models where the 21 cm power spectrum at $z = 6.5$ has a higher amplitude and a lower value of $\kappa$, and the pkSZ power spectrum peaks at $\ell \approx 2000$, similar to that discussed in Section 4.2.2 for the opKSZ case. Such models could potentially be excluded with the addition of a pkSZ data point at $\ell = 2000$ or more data on the 21 cm power spectrum at lower *k* values.

Regarding the LOFAR case, results are roughly similar to what was observed in Section 5.1.2 where the 21 cm mock data is taken at a fixed *k*-scale but different redshifts. Naturally, increasing the uncertainty on the measurements by approximately two orders of magnitude compared to that of the SKA (see Fig. 1), the constraining power of the data is decreased. The choice of a higher *z* mock 21 cm

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 9.** Posterior distributions on reionization mid- and endpoint (upper panel) and the morphology parameters $\log_{10}(\alpha_0)$ and $\kappa$ (bottom panel) when fitting 21 cm and pkSZ power spectra jointly, for 1000 h of MWA, LOFAR, & SKA observations, seen in grey, orange, and blue. The vertical and horizontal lines show the ‘true’ values used to generate the mock data.</figcaption></figure>

data point further reduces the ability of the forecast to place upper limits on the global reionization parameters (see Section 5.1). For example, as seen in grey in the upper panel of Fig. 9, we can constrain $z_{re} < 8.0$ and can only exclude only late reionization models for which $z_{end} > 5.5$. Compared to the MWA case, having higher $z$ 21 cm data leads to higher constraining power on the morphology parameters. While we can only place an upper limit on $\kappa > 0.1$, the constraint on $\log_{10}(\alpha_0) = 3.15 \pm 0.33$ is fairly comparable to that of 1k2z 1000 h (where $\log_{10}(\alpha_0) = 3.04 \pm 0.32$). While this is mostly related to the choice of data, for example, the hiz case can constrain $\log_{10}(\alpha_0) = 3.09 \pm 0.15$, our results indicate the potential of high-redshift LOFAR measurements on constraining on the morphology parameters of the EoR.

A possible improvement is the addition of an informed prior on the Thomson optical depth from the Planck Collaboration I (2020) measurement, which could potentially aid in increasing the precision of the forecast and constraining the parameter space of the global history parameters.

## 5.3 Increasing the number of data points

In previous sections, we used a minimalist approach by limiting the number of detected modes to test how much information would be obtained from early, partial 21 cm data, as a detection will only be achieved incrementally. On the one hand, with two data points at

<footnote n="12">¹²Note that we do not consider HERA in this section, since its ideal noise levels are lower than the ones for the SKA (Fig. 1).</footnote>
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Posterior distributions on reionization mid- and endpoint (upper panel) and the morphology parameters $log_{10}(\alpha_0)$ and $\kappa$ (bottom panel), considering a detection of the pkSZ power spectrum measurement at $\ell = 3000$ and of the 21 cm power spectrum at redshifts of $z = 6.5, 6.5, 7.8$ at $k = 0.1, 0.50, 0.5$ Mpc$^{-1}$ for a 1000 h of integration time with SKA, seen in brown and compared to the 1k2z 1000 h fiducial case shown in blue (see Fig. 5).</figcaption></figure>

a fixed $k$-scale ($k = 0.5$ Mpc$^{-1}$ for $z = 6.5, 7.8, 1k2z$ case), the forecast is sensitive to the redshift evolution of the 21 cm signal and we are able to well constrain the global reionization parameters $z_{re}$ and $z_{end}$. On the other hand, with two data points at fixed redshift but different scales ($z = 6.5$ for $k = 0.1, 0.5$ Mpc$^{-1}, 2k1z$ case), we can constrain the morphology parameters $log_{10}(\alpha_0)$ and $\kappa$. A natural way to extend these results is to venture into the multiscale regime (3k2z case, seen in brown in Fig. 10) by combining both of the aforementioned data sets such that the 21 cm power spectrum is measured at $z = 6, 5, 6.5, 7.8$ and $k = 0.1, 0.5, 0.5$ Mpc$^{-1}$ for 1000 h of integration with the SKA-Low.

We find that the forecast inherits the merits of the 1k2z and 2k1z cases as we access information on both the amplitude and the shape of the 21 cm power spectrum. Compared to our primary 1k2z case (seen in blue in Fig. 10), we improve the accuracy of our constraints of the reionization parameters: We get $z_{re} = 7.37 \pm 0.07$ and $z_{end} = 6.15 \pm 0.03$, significantly reducing the $1\sigma$ standard deviation (especially for $z_{end}$). Consequently, we can tightly constrain the parameter space around the true value and retrieve the correct Thomson optical depth with a precision of $\Delta\tau = \pm 0.001$ (see Appendix D). Hence, our method could provide an independent measurement of the optical depth, separate from the analysis of large-scale CMB polarization anisotropies (Planck Collaboration XLVII

2016), and potentially break well-known degeneracies with other cosmological parameters (Liu et al. 2016). Conversely, large-scale CMB data could be included in the analysis in order to improve constraints on the reionization parameters (Gorce, Douspis & Salvati 2022).

Constraints on the morphology parameters are significantly tighter. We measure $log_{10}(\alpha_0) = 3.12 \pm 0.01$ and $\kappa = 0.164 \pm 0.001$. This improvement is partly explained by the fact that adding data at small scales breaks the degeneracy between the morphology parameters. We can intuitively explain this by referring to the results in Section 5.1.2. The 3k2z case can be understood as the combination of the 1k2z and 2k1z cases and we see on Fig. 8 that accessing the redshift evolution of the 21 cm signal with data points at different redshifts removes the bias on the morphology parameters obtained with 2k1z.

Looking at these results, it is tempting to think that all the constraining power comes from the three 21 cm data points. However, we have also examined an additional case where we only consider the 21 cm mock data from 3k2z and exclude the pkSZ data point. We find that, while we can recover the morphology parameters with similar, albeit weaker, constraints, the global history parameters are significantly more degenerate. Indeed, the information the pkSZ data point provides on the duration of reionization remains crucial for the accuracy of the forecast. Measuring the pkSZ spectrum at different multipoles would enable tighter constraints on the shape parameter $\kappa$, whether at lower ($\ell = 2000$) or larger ($\ell = 5000$) multipoles. However, note that the former will be more challenging to achieve observationally because of the extremely large amplitude of the primary CMB temperature power spectrum on scales $\ell \lesssim 2000$.

If we have previously kept our results to simplistic cases in order to provide a proof of concept, these further tests show the full potential of our approach. As more and better measurements of the pkSZ and the 21 cm power spectra become available, this method will not only enable constraints on the global history of reionization but also on its morphology.

## 6 CONCLUSIONS

In this work, we examine how the fundamental relation between the patchy kSZ effect and the 21 cm signal can help us constrain the nature of cosmic reionization. In Section 2, we express and relate the 21 cm and kSZ power spectra through the parametrization of the electron power spectrum $P_{ee}(k, z)$ presented in Gorce et al. (2020). The resulting equation (6), while linking both observables, contains non-trivial cross-correlation and second-order terms of the ionized and density fields. We choose to ignore such terms in this work, leading to equation (8). Looking at the semi-numerical RSAGE simulation (Seiler et al. 2019), we find that this assumption leads to overestimating the 21 cm power on all scales and redshifts $z \lesssim 10$.

Aware of this limitation, we use the derived relation in a forecast analysis using an MCMC sampling method: We fit for the reionization mid- and endpoint, $z_{re}$ and $z_{end}$, as well as for two reionization morphology parameters $log_{10}(\alpha_0)$ and $\kappa$. Indeed, for each of set of these parameters, we can derive the electron power spectrum and, in turn, the kSZ and 21 cm power spectra. We generate mock observations of these with the parameter values given in Table 1, assume a ten per cent error bar for the pkSZ data point at $\ell = 3000$, and include thermal noise and sample variance in the 21 cm measurement errors, typically for 1000 h of integration with SKA-Low. We follow a minimalist approach as we do not assume a full detection of the kSZ and 21 cm power spectra over a range of scales, but assess the constraining power of only a few observed data points. Such an approach gives us insight on how much we can learn
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Posterior distributions on reionization mid- and endpoint (upper panel) and the morphology parameters $log_{10}(\alpha_0)$ and $\kappa$ (bottom panel), considering a detection of the pkSZ power spectrum measurement at $\ell = 3000$ and of the 21 cm power spectrum at redshifts of $z = 6.5, 6.5, 7.8$ at $k = 0.1, 0.50, 0.5$ Mpc$^{-1}$ for a 1000 h of integration time with SKA, seen in brown and compared to the 1k2z 1000 h fiducial case shown in blue (see Fig. 5).</figcaption></figure>

a fixed $k$-scale ($k = 0.5$ Mpc$^{-1}$ for $z = 6.5, 7.8, 1k2z$ case), the forecast is sensitive to the redshift evolution of the 21 cm signal and we are able to well constrain the global reionization parameters $z_{re}$ and $z_{end}$. On the other hand, with two data points at fixed redshift but different scales ($z = 6.5$ for $k = 0.1, 0.5$ Mpc$^{-1}, 2k1z$ case), we can constrain the morphology parameters $log_{10}(\alpha_0)$ and $\kappa$. A natural way to extend these results is to venture into the multiscale regime (3k2z case, seen in brown in Fig. 10) by combining both of the aforementioned data sets such that the 21 cm power spectrum is measured at $z = 6, 5, 6.5, 7.8$ and $k = 0.1, 0.5, 0.5$ Mpc$^{-1}$ for 1000 h of integration with the SKA-Low.

We find that the forecast inherits the merits of the 1k2z and 2k1z cases as we access information on both the amplitude and the shape of the 21 cm power spectrum. Compared to our primary 1k2z case (seen in blue in Fig. 10), we improve the accuracy of our constraints of the reionization parameters: We get $z_{re} = 7.37 \pm 0.07$ and $z_{end} = 6.15 \pm 0.03$, significantly reducing the $1\sigma$ standard deviation (especially for $z_{end}$). Consequently, we can tightly constrain the parameter space around the true value and retrieve the correct Thomson optical depth with a precision of $\Delta\tau = \pm 0.001$ (see Appendix D). Hence, our method could provide an independent measurement of the optical depth, separate from the analysis of large-scale CMB polarization anisotropies (Planck Collaboration XLVII

2016), and potentially break well-known degeneracies with other cosmological parameters (Liu et al. 2016). Conversely, large-scale CMB data could be included in the analysis in order to improve constraints on the reionization parameters (Gorce, Douspis & Salvati 2022).

Constraints on the morphology parameters are significantly tighter. We measure $log_{10}(\alpha_0) = 3.12 \pm 0.01$ and $\kappa = 0.164 \pm 0.001$. This improvement is partly explained by the fact that adding data at small scales breaks the degeneracy between the morphology parameters. We can intuitively explain this by referring to the results in Section 5.1.2. The 3k2z case can be understood as the combination of the 1k2z and 2k1z cases and we see on Fig. 8 that accessing the redshift evolution of the 21 cm signal with data points at different redshifts removes the bias on the morphology parameters obtained with 2k1z.

Looking at these results, it is tempting to think that all the constraining power comes from the three 21 cm data points. However, we have also examined an additional case where we only consider the 21 cm mock data from 3k2z and exclude the pkSZ data point. We find that, while we can recover the morphology parameters with similar, albeit weaker, constraints, the global history parameters are significantly more degenerate. Indeed, the information the pkSZ data point provides on the duration of reionization remains crucial for the accuracy of the forecast. Measuring the pkSZ spectrum at different multipoles would enable tighter constraints on the shape parameter $\kappa$, whether at lower ($\ell = 2000$) or larger ($\ell = 5000$) multipoles. However, note that the former will be more challenging to achieve observationally because of the extremely large amplitude of the primary CMB temperature power spectrum on scales $\ell \lesssim 2000$.

If we have previously kept our results to simplistic cases in order to provide a proof of concept, these further tests show the full potential of our approach. As more and better measurements of the pkSZ and the 21 cm power spectra become available, this method will not only enable constraints on the global history of reionization but also on its morphology.

## 6 CONCLUSIONS

In this work, we examine how the fundamental relation between the patchy kSZ effect and the 21 cm signal can help us constrain the nature of cosmic reionization. In Section 2, we express and relate the 21 cm and kSZ power spectra through the parametrization of the electron power spectrum $P_{ee}(k, z)$ presented in Gorce et al. (2020). The resulting equation (6), while linking both observables, contains non-trivial cross-correlation and second-order terms of the ionized and density fields. We choose to ignore such terms in this work, leading to equation (8). Looking at the semi-numerical RSAGE simulation (Seiler et al. 2019), we find that this assumption leads to overestimating the 21 cm power on all scales and redshifts $z \lesssim 10$.

Aware of this limitation, we use the derived relation in a forecast analysis using an MCMC sampling method: We fit for the reionization mid- and endpoint, $z_{re}$ and $z_{end}$, as well as for two reionization morphology parameters $log_{10}(\alpha_0)$ and $\kappa$. Indeed, for each of set of these parameters, we can derive the electron power spectrum and, in turn, the kSZ and 21 cm power spectra. We generate mock observations of these with the parameter values given in Table 1, assume a ten per cent error bar for the pkSZ data point at $\ell = 3000$, and include thermal noise and sample variance in the 21 cm measurement errors, typically for 1000 h of integration with SKA-Low. We follow a minimalist approach as we do not assume a full detection of the kSZ and 21 cm power spectra over a range of scales, but assess the constraining power of only a few observed data points. Such an approach gives us insight on how much we can learn
about reionization with early measurements of both observables. Our findings are as follows.

With as few as 100 h of SKA-Low observations and two 21 cm measurements at $k = 0.5$ Mpc$^{-1}$ and $z = 6.5, 7.8$ (referred to as the 1k2z case, shown in blue in Fig. 5 and throughout this work), we are able to constrain the global reionization parameters $z_{re} = 7.42 \pm 0.13$ and $z_{end} = 6.15 \pm 0.04$ and provide limits on the morphology parameters $\log_{10}(\alpha_0) > 2.64$ and $\kappa < 0.25$ (Table 3). We demonstrate that this constraining power stems from the complementary nature of both data sets: The 21 cm signal is mostly sensitive to the endpoint of reionization, while the pkSZ effect is inherently sensitive to its duration $dz = z_{re} - z_{end}$. Our results show that, even the early stages of SKA-Low operations, its observations will provide valuable insight into the global history of the EoR despite a minimal knowledge of the properties of the first stars and galaxies. The constraints on the reionization global history are further improved when the fitted mock 21 cm measurements are at different redshifts (but a given k-scale), as the data now trace the redshift evolution of the 21 cm signal.

We show that the choice of redshift for the mock data set plays an important role. First, any low-$z$ measurement excludes earlier reionization models, for which the 21 cm signal would be zero. Second, at high redshift, when the global ionization fraction is $x_v \ll 1$, the data will primarily determined by $\log_{10}(\alpha_0)$ and $\kappa$. In our tests, for 21 cm mock measurements at high(er) redshift (hiz case, $z = 7.8$ and $10.4$ at $k = 0.5$ Mpc$^{-1}$), the constraints on the global reionization parameters $z_{re}$ and $z_{end}$ are loosened compared to the 1k2z case, at $z = 6.5$ and $7.8$ (Fig. 8). On the other hand, it is precisely the constraining of the high-$z$ 21 cm power spectrum which allows for firmer constraints on the morphological parameters: $\log_{10}(\alpha_0) = 3.10 \pm 0.17$ and $\kappa = 0.17 \pm 0.02$.

We also observe that the k-scale of the data significantly influences forecast constraints. Indeed, the inside-out nature of cosmic reionization will impact the shape of the 21 cm power spectrum, specifically, its high-k 'tail', which encodes information about reionization morphology on small scales. By choosing two mock 21 cm power spectra at a fixed redshift ($2k_1z, z = 6.5$ at $k = 0.1, 0.5$ Mpc$^{-1}$), we find that the constraints on the global reionization parameters are weakened to upper limits (Fig. 8, in yellow) but constraints on the morphology parameters are improved. However, the recovered values are biased by 6 per cent. We find that major differences between the allowed models (which reach the mid- and endpoint later than the true model) and the true 21 cm power spectrum are most distinguishable for modes $k > 1.0$ Mpc$^{-1}$, for which the model does not possess information. This implies that while multiscale observations of the 21 cm power spectrum at one redshift contain information on the whole process of reionization, it is the highest k modes, which are the most sensitive. However, such measurements are non-trivial because the measurement noise scales with k-scale (see Fig. 1), making them less appealing in the context of early SKA-Low observations.

It is also worth noting that we limit our main results to simplistic cases with two 21 cm power spectrum observations in order to provide a proof-of-concept approach and to better showcase the workings of the forecast. Further tests where we increase the number of data points from two to three (e.g. the 3k2z case seen in Section 5.3) show the full potential of this approach. Highlighted in brown in Fig. 10, we see that as more and better measurements of the pkSZ and the 21 cm power spectrum become available, this method will not only enable constraints on the global history of reionization but also its morphology. We obtain accurate constrains of the reionization parameters $z_{re} = 7.37 \pm 0.07$ and $z_{end} = 6.15 \pm 0.03$ and the morphology parameters $\log_{10}(\alpha_0) = 3.12 \pm 0.01$ and $\kappa =$

$0.164 \pm 0.001$, significantly reducing the $1\sigma$ standard deviation, compared to 1k2z.

Finally, we explore forecast performance on mock data from operating telescopes such as MWA and LOFAR, based current upper limits. (Patil et al. 2017; Mertens et al. 2020; Trott et al. 2020). Naturally, we find the mock data less constraining than in the case of 1k2z (see Fig. 9 in grey and orange) as can be expected from the higher uncertainty on the measurements. However, detections by both telescopes can place firm upper limits on the midpoint of reionization $z_{re} \le 8$. Moreover, the low-redshift observations of MWA (for the same set-up as 1k2z) can also place a lower limit on the end of reionization $z_{end} \ge 6.5$. Meanwhile, data from LOFAR (chosen at $k = 0.1$ Mpc$^{-1}$ for $z = 8.3, 9.1$) is more sensitive to the morphology parameters, notably constraining $\log_{10}(\alpha_0) = 3.15 \pm 0.33$. This implies that even before the first-light of SKA-Low, the ongoing improvements of current 21 cm power spectrum upper limits (see fig. 6 of Raste et al. 2021, for an example) are likely to soon begin constraining the properties of reionization.

The main limitation of our model is the omission of the higher order and cross-correlation terms in our simplified expression of the 21 cm power spectrum in equation (8). Potential future developments would be develop analytic models of the cross-power spectrum and include them in the derivation (McQuinn et al. 2005; Schneider, Schaeffer & Giri 2023). However, this approach would most likely not capture the full complexity present at all k-scales and could be model dependent. Another option would be to model the higher order and the cross-terms as constant biases and account for them as nuisances parameters in the forecast (see the discussion in Appendix C). Lastly, we plan to conduct our analysis for 21 cm data over broader ranges of redshift and k-scales to explore the potential of joint analysis to identify and remove systematics such as residual foregrounds [see Bégin, Liu & Gorce (2022) for such an analysis using the global 21 cm signal].

ACKNOWLEDGEMENTS

The authors thank Adrian Liu, Jack Line, Leon Koopmans, and Cathryn Trott for useful discussions about different parts of this work. We thank the anonymous reviewer for insightful and constructive comments. IG acknowledges the support from the Gustaf och Ellen Kobbs stipendiestiftelse and the Alva and Lennart Dahlmark Research Grant. At the time when most of this work was conducted, AG's work was supported by the McGill Astrophysics Fellowship funded by the Trottier Chair in Astrophysics, as well as the Canadian Institute for Advanced Research (CIFAR) Azrieli Global Scholars program and the Canada 150 Programme. GM's research is supported by the Swedish Research Council project grant 2020-04691_VR.

This research made use of MATPLOTLIB, a PYTHON library for publication-quality graphics (Hunter 2007), of SCIPY, a PYTHON-based ecosystem of open-source software for mathematics, science, and engineering (Jones et al. 2001) - including NUMPY (Oliphant 2006). The MCMC was performed with the EMCEE sampler (Foreman-Mackey et al. 2013) and posterior distributions were drawn thanks to a modified version of the corner PYTHON package (Foreman-Mackey 2016). This research was enabled in part by support provided by Calcul Québec and the Digital Research Alliance of Canada.<fnref n="13" /> Computations were notably performed on the Cedar supercomputer, hosted by the Simon Fraser University.

<footnote n="13">¹³See http://www.alliancecan.ca and https://www.calculquebec.ca.</footnote>
about reionization with early measurements of both observables. Our findings are as follows.

With as few as 100 h of SKA-Low observations and two 21 cm measurements at $k = 0.5$ Mpc$^{-1}$ and $z = 6.5, 7.8$ (referred to as the 1k2z case, shown in blue in Fig. 5 and throughout this work), we are able to constrain the global reionization parameters $z_{re} = 7.42 \pm 0.13$ and $z_{end} = 6.15 \pm 0.04$ and provide limits on the morphology parameters $\log_{10}(\alpha_0) > 2.64$ and $\kappa < 0.25$ (Table 3). We demonstrate that this constraining power stems from the complementary nature of both data sets: The 21 cm signal is mostly sensitive to the endpoint of reionization, while the pkSZ effect is inherently sensitive to its duration $dz = z_{re} - z_{end}$. Our results show that, even the early stages of SKA-Low operations, its observations will provide valuable insight into the global history of the EoR despite a minimal knowledge of the properties of the first stars and galaxies. The constraints on the reionization global history are further improved when the fitted mock 21 cm measurements are at different redshifts (but a given k-scale), as the data now trace the redshift evolution of the 21 cm signal.

We show that the choice of redshift for the mock data set plays an important role. First, any low-$z$ measurement excludes earlier reionization models, for which the 21 cm signal would be zero. Second, at high redshift, when the global ionization fraction is $x_v \ll 1$, the data will primarily determined by $\log_{10}(\alpha_0)$ and $\kappa$. In our tests, for 21 cm mock measurements at high(er) redshift (hiz case, $z = 7.8$ and $10.4$ at $k = 0.5$ Mpc$^{-1}$), the constraints on the global reionization parameters $z_{re}$ and $z_{end}$ are loosened compared to the 1k2z case, at $z = 6.5$ and $7.8$ (Fig. 8). On the other hand, it is precisely the constraining of the high-$z$ 21 cm power spectrum which allows for firmer constraints on the morphological parameters: $\log_{10}(\alpha_0) = 3.10 \pm 0.17$ and $\kappa = 0.17 \pm 0.02$.

We also observe that the k-scale of the data significantly influences forecast constraints. Indeed, the inside-out nature of cosmic reionization will impact the shape of the 21 cm power spectrum, specifically, its high-k 'tail', which encodes information about reionization morphology on small scales. By choosing two mock 21 cm power spectra at a fixed redshift ($2k_{1z}, z = 6.5$ at $k = 0.1, 0.5$ Mpc$^{-1}$), we find that the constraints on the global reionization parameters are weakened to upper limits (Fig. 8, in yellow) but constraints on the morphology parameters are improved. However, the recovered values are biased by 6 per cent. We find that major differences between the allowed models (which reach the mid- and endpoint later than the true model) and the true 21 cm power spectrum are most distinguishable for modes $k > 1.0$ Mpc$^{-1}$, for which the model does not possess information. This implies that while multiscale observations of the 21 cm power spectrum at one redshift contain information on the whole process of reionization, it is the highest k modes, which are the most sensitive. However, such measurements are non-trivial because the measurement noise scales with k-scale (see Fig. 1), making them less appealing in the context of early SKA-Low observations.

It is also worth noting that we limit our main results to simplistic cases with two 21 cm power spectrum observations in order to provide a proof-of-concept approach and to better showcase the workings of the forecast. Further tests where we increase the number of data points from two to three (e.g. the 3k2z case seen in Section 5.3) show the full potential of this approach. Highlighted in brown in Fig. 10, we see that as more and better measurements of the pkSZ and the 21 cm power spectrum become available, this method will not only enable constraints on the global history of reionization but also its morphology. We obtain accurate constrains of the reionization parameters $z_{re} = 7.37 \pm 0.07$ and $z_{end} = 6.15 \pm 0.03$ and the morphology parameters $\log_{10}(\alpha_0) = 3.12 \pm 0.01$ and $\kappa =$

$0.164 \pm 0.001$, significantly reducing the $1\sigma$ standard deviation, compared to 1k2z.

Finally, we explore forecast performance on mock data from operating telescopes such as MWA and LOFAR, based current upper limits. (Patil et al. 2017; Mertens et al. 2020; Trott et al. 2020). Naturally, we find the mock data less constraining than in the case of 1k2z (see Fig. 9 in grey and orange) as can be expected from the higher uncertainty on the measurements. However, detections by both telescopes can place firm upper limits on the midpoint of reionization $z_{re} \le 8$. Moreover, the low-redshift observations of MWA (for the same set-up as 1k2z) can also place a lower limit on the end of reionization $z_{end} \ge 6.5$. Meanwhile, data from LOFAR (chosen at $k = 0.1$ Mpc$^{-1}$ for $z = 8.3, 9.1$) is more sensitive to the morphology parameters, notably constraining $\log_{10}(\alpha_0) = 3.15 \pm 0.33$. This implies that even before the first-light of SKA-Low, the ongoing improvements of current 21 cm power spectrum upper limits (see fig. 6 of Raste et al. 2021, for an example) are likely to soon begin constraining the properties of reionization.

The main limitation of our model is the omission of the higher order and cross-correlation terms in our simplified expression of the 21 cm power spectrum in equation (8). Potential future developments would be develop analytic models of the cross-power spectrum and include them in the derivation (McQuinn et al. 2005; Schneider, Schaeffer & Giri 2023). However, this approach would most likely not capture the full complexity present at all k-scales and could be model dependent. Another option would be to model the higher order and the cross-terms as constant biases and account for them as nuisances parameters in the forecast (see the discussion in Appendix C). Lastly, we plan to conduct our analysis for 21 cm data over broader ranges of redshift and k-scales to explore the potential of joint analysis to identify and remove systematics such as residual foregrounds [see Bégin, Liu & Gorce (2022) for such an analysis using the global 21 cm signal].

## ACKNOWLEDGEMENTS

The authors thank Adrian Liu, Jack Line, Leon Koopmans, and Cathryn Trott for useful discussions about different parts of this work. We thank the anonymous reviewer for insightful and constructive comments. IG acknowledges the support from the Gustaf och Ellen Kobbs stipendiestiftelse and the Alva and Lennart Dahlmark Research Grant. At the time when most of this work was conducted, AG's work was supported by the McGill Astrophysics Fellowship funded by the Trottier Chair in Astrophysics, as well as the Canadian Institute for Advanced Research (CIFAR) Azrieli Global Scholars program and the Canada 150 Programme. GM's research is supported by the Swedish Research Council project grant 2020-04691_VR.

This research made use of MATPLOTLIB, a PYTHON library for publication-quality graphics (Hunter 2007), of SCIPY, a PYTHON-based ecosystem of open-source software for mathematics, science, and engineering (Jones et al. 2001) - including NUMPY (Oliphant 2006). The MCMC was performed with the EMCEE sampler (Foreman-Mackey et al. 2013) and posterior distributions were drawn thanks to a modified version of the corner PYTHON package (Foreman-Mackey 2016). This research was enabled in part by support provided by Calcul Québec and the Digital Research Alliance of Canada.<fnref n="13" /> Computations were notably performed on the Cedar supercomputer, hosted by the Simon Fraser University.

<footnote n="13">¹³See http://www.alliancecan.ca and https://www.calculquebec.ca.</footnote>
## DATA AVAILABILITY

The sampling code used to obtain all the results presented in this article is publicly available at https://github.com/adeliegorce/forecas t.kszx21. The (trained) random forests used to generate the angular patchy kSZ power spectra given a set of cosmological parameters can be found at https://szdb.osups.universite-paris-saclay.fr. Any additional information or data can be requested from the authors.

## REFERENCES

* Alvarez M. A., 2016, ApJ, 824, 118

* Aubert D., Deparis N., Ocvirk P., 2015, MNRAS, 454, 1012

* Battaglia N., Trac H., Cen R., Loeb A., 2013, ApJ, 776, 81

* Bégin J.-M., Liu A., Gorce A., 2022, Phys. Rev. D, 105, 083503

* Bosman S. E. I. et al., 2022, MNRAS, 514, 55

* Chapman E. et al., 2013, MNRAS, 429, 165

* Chen N., Trac H., Mukherjee S., Cen R., 2023, ApJ, 943, 138

* Dixon K. L., Iliev I. T., Mellema G., Ahn K., Shapiro P. R., 2016, MNRAS, 456, 3011

* Douspis M., Aghanim N., Ilić S., Langer M., 2015, A&A, 580, L4

* Douspis M., Salvati L., Gorce A., Aghanim N., 2022, A&A, 659, A99

* Field G. B., 1957, AJ, 62, 15

* Foreman-Mackey D., 2016, J. Open Source Softw., 1, 24

* Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306

* Furlanetto S. R., Zaldarriaga M., Hernquist L., 2004, ApJ, 613, 1

* Furlanetto S. R., McQuinn M., Hernquist L., 2006, MNRAS, 365, 115

* Gaikwad P. et al., 2023, MNRAS, 525, 4093

* Gelman A., Rubin D. B., 1992, Stat. Sci., 7, 457

* George E. M. et al., 2015, ApJ, 799, 177

* Georgiev I., Mellema G., Giri S. K., Mondal R., 2022, MNRAS, 513, 5109

* Gorce A., Ilić S., Douspis M., Aubert D., Langer M., 2020, A&A, 640, A90

* Gorce A., Douspis M., Salvati L., 2022, A&A, 662, A122

* HERA Collaboration, 2023, ApJ, 945, 124

* Howlett C., Lewis A., Hall A., Challinor A., 2012, J. Cosmol. Astropart. Phys., 1204, 027

* Hunter J. D., 2007, Comput. Sci. Eng., 9, 90

* Hutter A., 2018, MNRAS, 477, 1549

* Iliev I. T., Pen U.-L., Bond J. R., Mellema G., Shapiro P. R., 2007, ApJ, 660, 933

* Iliev I. T., Mellema G., Ahn K., Shapiro P. R., Mao Y., Pen U.-L., 2014, MNRAS, 439, 725

* Jones E., Oliphant T., Peterson P. et al., 2001, SciPy: Open source scientific tools for Python. http://www.scipy.org/

* Koopmans L. et al., 2015, in Advancing Astrophysics with the Square Kilometre Array (AASKA14). p. 1, preprint (arXiv:1505.07568)

* Lewis A., 2013, Phys. Rev. D, 87, 103529

* Lewis A., Bridle S., 2002, Phys. Rev. D, 66, 103511

* Lewis A., Challinor A., Lasenby A., 2000, ApJ, 538, 473

* Lidz A., Zahn O., McQuinn M., Zaldarriaga M., Dutta S., Hernquist L., 2007, ApJ, 659, 865

* Liu A., Parsons A. R., Trott C. M., 2014, Phys. Rev. D, 90, 023018

* Liu A., Pritchard J. R., Allison R., Parsons A. R., Seljak U., Sherwin B. D., 2016, Phys. Rev. D, 93, 043013

* Maniyar A., Béthermin M., Lagache G., 2021, A&A, 645, A40

* McQuinn M., Furlanetto S. R., Hernquist L., Zahn O., Zaldarriaga M., 2005, ApJ, 630, 643

* McQuinn M., Zahn O., Zaldarriaga M., Hernquist L., Furlanetto S. R., 2006, ApJ, 653, 815

* Mellema G. et al., 2013, Exp. Astron., 36, 235

* Mertens F. G., Ghosh A., Koopmans L. V. E., 2018, MNRAS, 478, 3640

* Mertens F. G. et al., 2020, MNRAS, 493, 1662

* Mesinger A., McQuinn M., Spergel D. N., 2012, MNRAS, 422, 1403

* Oliphant T., 2006, NumPy: A guide to NumPy. Trelgol Publishing, USA

* Paciga G. et al., 2013, MNRAS, 433, 639

* Park H., Shapiro P. R., Komatsu E., Iliev I. T., Ahn K., Mellema G., 2013, ApJ, 769, 93

* Patil A. H. et al., 2017, ApJ, 838, 65

* Planck CollaborationXLVII, 2016, A&A, 596, A108

* Planck CollaborationI, 2020, A&A, 641, A1

* Pritchard J. R., Loeb A., 2012, Rep. Prog. Phys., 75, 086901

* Raghunathan S., Omori Y., 2023, ApJ, 954, 83

* Raste J., Kulkarni G., Keating L. C., Haehnelt M. G., Chardin J., Aubert D., 2021, MNRAS, 507, 4684

* Reichardt C. L. et al., 2021, ApJ, 908, 199

* , in Ruhl J. Bradford C. M. et al., 2004 et al., eds, Proc. SPIE Conf. Ser. Vol. 5498, Z-Spec: A Broadband Millimeter-wave Grating Spectrometer: Design, Construction, and First Cryogenic Measurements. SPIE, Bellingham, p. 11

* Schneider A., Schaeffer T., Giri S. K., 2023, Phys. Rev. D, 108, 043030

* Seiler J., Hutter A., Sinha M., Croton D., 2018, MNRAS, 480, L33

* Seiler J., Hutter A., Sinha M., Croton D., 2019, MNRAS, 487, 5739

* Shaw L. D., Rudd D. H., Nagai D., 2012, ApJ, 756, 15

* Sunyaev R. A., Zeldovich I. B., 1980, ARA&A, 18, 537

* The HERA Collaboration, 2022, ApJ, 925, 221

* Trott C. M. et al., 2020, MNRAS, 493, 4711

* Wise J. H., 2019, preprint (arXiv:1907.06653)

* Yoshiura S. et al., 2021, MNRAS, 505, 4775

* Zahn O., Zaldarriaga M., Hernquist L., McQuinn M., 2005, ApJ, 630, 657

* Zahn O. et al., 2012, ApJ, 756, 65

## APPENDIX A: RECONSTRUCTING THE PKSZ FROM A MEASUREMENT OF THE 21 CM POWER SPECTRUM

Aware of the limitations of our reconstruction model presented in Section 4.1, we now attempt to use the simplified 21 cm power spectrum model to reconstruct the pkSZ angular power spectrum. That is, we compute numerically the 21 cm power spectrum from our simulation and then use it to reconstruct the electron power spectrum following equation (6). We then plug the reconstructed $P_{ee}(k, z)$ in equation (9) to obtain the patchy kSZ angular power. We compare results for different levels of reconstruction precision in Fig. A1, similarly to Fig. 3. Note that this figure is only intended as a comparison of the potential of the reconstruction, not as a precise estimate of the pkSZ power of the simulation.<fnref n="14" /> Indeed, the power is set to zero on $k$-modes not covered by the simulation. The results are slightly different from the ones presented in Fig. 3. Here, the pkSZ power is systematically under-estimated, as long as all the terms of equation (6) are not included. Despite our approximation largely underestimating the pkSZ power, the reconstructed power still lies within current error bars on pkSZ power spectrum measurements at $\ell = 3000$ (Reichardt et al. 2021).

Instead of reconstructing $P_{21}$ from equation (7) and missing out on crucial power, one could use direct measurements of the 21 cm power spectrum to reconstruct the pkSZ power and compare to CMB observations. To do so, the former must be integrated over a wide range of scales. However, for now, the 21 cm power spectrum has not been measured with sufficient precision on a wide enough range. We must therefore assess which scales contribute the most to the final pkSZ power and are necessary for this reconstruction. The relative contribution of each $k$-mode to the pkSZ $C_\ell$ at $\ell = 3000$ had already been estimated in Gorce et al. (2020), but here we extend this result to the range $1000 \le \ell < 10000$ and show the result in Fig. A2. We

<footnote n="14">¹⁴Remember that all cosmological and reionization parameters are derived from the RSAGE simulation.</footnote>
## DATA AVAILABILITY

The sampling code used to obtain all the results presented in this article is publicly available at https://github.com/adeliegorce/forecas t.kszx21. The (trained) random forests used to generate the angular patchy kSZ power spectra given a set of cosmological parameters can be found at https://szdb.osups.universite-paris-saclay.fr. Any additional information or data can be requested from the authors.

## REFERENCES

* Alvarez M. A., 2016, ApJ, 824, 118

* Aubert D., Deparis N., Ocvirk P., 2015, MNRAS, 454, 1012

* Battaglia N., Trac H., Cen R., Loeb A., 2013, ApJ, 776, 81

* Bégin J.-M., Liu A., Gorce A., 2022, Phys. Rev. D, 105, 083503

* Bosman S. E. I. et al., 2022, MNRAS, 514, 55

* Chapman E. et al., 2013, MNRAS, 429, 165

* Chen N., Trac H., Mukherjee S., Cen R., 2023, ApJ, 943, 138

* Dixon K. L., Iliev I. T., Mellema G., Ahn K., Shapiro P. R., 2016, MNRAS, 456, 3011

* Douspis M., Aghanim N., Ilić S., Langer M., 2015, A&A, 580, L4

* Douspis M., Salvati L., Gorce A., Aghanim N., 2022, A&A, 659, A99

* Field G. B., 1957, AJ, 62, 15

* Foreman-Mackey D., 2016, J. Open Source Softw., 1, 24

* Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306

* Furlanetto S. R., Zaldarriaga M., Hernquist L., 2004, ApJ, 613, 1

* Furlanetto S. R., McQuinn M., Hernquist L., 2006, MNRAS, 365, 115

* Gaikwad P. et al., 2023, MNRAS, 525, 4093

* Gelman A., Rubin D. B., 1992, Stat. Sci., 7, 457

* George E. M. et al., 2015, ApJ, 799, 177

* Georgiev I., Mellema G., Giri S. K., Mondal R., 2022, MNRAS, 513, 5109

* Gorce A., Ilić S., Douspis M., Aubert D., Langer M., 2020, A&A, 640, A90

* Gorce A., Douspis M., Salvati L., 2022, A&A, 662, A122

* HERA Collaboration, 2023, ApJ, 945, 124

* Howlett C., Lewis A., Hall A., Challinor A., 2012, J. Cosmol. Astropart. Phys., 1204, 027

* Hunter J. D., 2007, Comput. Sci. Eng., 9, 90

* Hutter A., 2018, MNRAS, 477, 1549

* Iliev I. T., Pen U.-L., Bond J. R., Mellema G., Shapiro P. R., 2007, ApJ, 660, 933

* Iliev I. T., Mellema G., Ahn K., Shapiro P. R., Mao Y., Pen U.-L., 2014, MNRAS, 439, 725

* Jones E., Oliphant T., Peterson P. et al., 2001, SciPy: Open source scientific tools for Python. http://www.scipy.org/

* Koopmans L. et al., 2015, in Advancing Astrophysics with the Square Kilometre Array (AASKA14). p. 1, preprint (arXiv:1505.07568)

* Lewis A., 2013, Phys. Rev. D, 87, 103529

* Lewis A., Bridle S., 2002, Phys. Rev. D, 66, 103511

* Lewis A., Challinor A., Lasenby A., 2000, ApJ, 538, 473

* Lidz A., Zahn O., McQuinn M., Zaldarriaga M., Dutta S., Hernquist L., 2007, ApJ, 659, 865

* Liu A., Parsons A. R., Trott C. M., 2014, Phys. Rev. D, 90, 023018

* Liu A., Pritchard J. R., Allison R., Parsons A. R., Seljak U., Sherwin B. D., 2016, Phys. Rev. D, 93, 043013

* Maniyar A., Béthermin M., Lagache G., 2021, A&A, 645, A40

* McQuinn M., Furlanetto S. R., Hernquist L., Zahn O., Zaldarriaga M., 2005, ApJ, 630, 643

* McQuinn M., Zahn O., Zaldarriaga M., Hernquist L., Furlanetto S. R., 2006, ApJ, 653, 815

* Mellema G. et al., 2013, Exp. Astron., 36, 235

* Mertens F. G., Ghosh A., Koopmans L. V. E., 2018, MNRAS, 478, 3640

* Mertens F. G. et al., 2020, MNRAS, 493, 1662

* Mesinger A., McQuinn M., Spergel D. N., 2012, MNRAS, 422, 1403

* Oliphant T., 2006, NumPy: A guide to NumPy. Trelgol Publishing, USA

* Paciga G. et al., 2013, MNRAS, 433, 639

* Park H., Shapiro P. R., Komatsu E., Iliev I. T., Ahn K., Mellema G., 2013, ApJ, 769, 93

* Patil A. H. et al., 2017, ApJ, 838, 65

* Planck CollaborationXLVII, 2016, A&A, 596, A108

* Planck CollaborationI, 2020, A&A, 641, A1

* Pritchard J. R., Loeb A., 2012, Rep. Prog. Phys., 75, 086901

* Raghunathan S., Omori Y., 2023, ApJ, 954, 83

* Raste J., Kulkarni G., Keating L. C., Haehnelt M. G., Chardin J., Aubert D., 2021, MNRAS, 507, 4684

* Reichardt C. L. et al., 2021, ApJ, 908, 199

, in Ruhl J. Bradford C. M. et al., 2004 et al., eds, Proc. SPIE Conf. Ser. Vol. 5498, Z-Spec: A Broadband Millimeter-wave Grating Spectrometer: Design, Construction, and First Cryogenic Measurements. SPIE, Bellingham, p. 11

* Schneider A., Schaeffer T., Giri S. K., 2023, Phys. Rev. D, 108, 043030

* Seiler J., Hutter A., Sinha M., Croton D., 2018, MNRAS, 480, L33

* Seiler J., Hutter A., Sinha M., Croton D., 2019, MNRAS, 487, 5739

* Shaw L. D., Rudd D. H., Nagai D., 2012, ApJ, 756, 15

* Sunyaev R. A., Zeldovich I. B., 1980, ARA&A, 18, 537

* The HERA Collaboration, 2022, ApJ, 925, 221

* Trott C. M. et al., 2020, MNRAS, 493, 4711

* Wise J. H., 2019, preprint (arXiv:1907.06653)

* Yoshiura S. et al., 2021, MNRAS, 505, 4775

* Zahn O., Zaldarriaga M., Hernquist L., McQuinn M., 2005, ApJ, 630, 657

* Zahn O. et al., 2012, ApJ, 756, 65

## APPENDIX A: RECONSTRUCTING THE PKSZ FROM A MEASUREMENT OF THE 21 CM POWER SPECTRUM

Aware of the limitations of our reconstruction model presented in Section 4.1, we now attempt to use the simplified 21 cm power spectrum model to reconstruct the pkSZ angular power spectrum. That is, we compute numerically the 21 cm power spectrum from our simulation and then use it to reconstruct the electron power spectrum following equation (6). We then plug the reconstructed $P_{ee}(k, z)$ in equation (9) to obtain the patchy kSZ angular power. We compare results for different levels of reconstruction precision in Fig. A1, similarly to Fig. 3. Note that this figure is only intended as a comparison of the potential of the reconstruction, not as a precise estimate of the pkSZ power of the simulation.<fnref n="14" /> Indeed, the power is set to zero on $k$-modes not covered by the simulation. The results are slightly different from the ones presented in Fig. 3. Here, the pkSZ power is systematically under-estimated, as long as all the terms of equation (6) are not included. Despite our approximation largely underestimating the pkSZ power, the reconstructed power still lies within current error bars on pkSZ power spectrum measurements at $\ell = 3000$ (Reichardt et al. 2021).

Instead of reconstructing $P_{21}$ from equation (7) and missing out on crucial power, one could use direct measurements of the 21 cm power spectrum to reconstruct the pkSZ power and compare to CMB observations. To do so, the former must be integrated over a wide range of scales. However, for now, the 21 cm power spectrum has not been measured with sufficient precision on a wide enough range. We must therefore assess which scales contribute the most to the final pkSZ power and are necessary for this reconstruction. The relative contribution of each $k$-mode to the pkSZ $C_\ell$ at $\ell = 3000$ had already been estimated in Gorce et al. (2020), but here we extend this result to the range $1000 \le \ell < 10000$ and show the result in Fig. A2. We

<footnote n="14">¹⁴Remember that all cosmological and reionization parameters are derived from the RSAGE simulation.</footnote>
<figure><img src="image_12.png" /><figcaption>**Figure A1.** The pkSZ angular power spectrum obtained for the true electron power spectrum measured in the RSAGE const simulation, and for reconstructions based on the 21 cm power spectrum at various precision levels. The colour scheme is similar to Fig. 3. The results are compared to the SPT data point on the pkSZ power (Reichardt et al. 2021).</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure A2.** Physical scales contributing to the final pkSZ power at different angular multipoles $\ell$. The result for $\ell = 3000$ is highlighted in bold.</figcaption></figure>

see that, for all multipoles, most of the power stems from $10^{-3} < k/[\text{Mpc}^{-1}] < 1$ with only the upper limit increasing as $\ell$ increases.

The spherical Fourier modes used in Fig. A2 are effectively the sum of a transversal and a longitudinal mode $k = \sqrt{k_{\perp}^2 + k_{\parallel}^2}$. The Fourier modes probed by an interferometer measuring signal at a frequency $\nu$ corresponding to a redshift $z = \nu_{21}/\nu - 1$ depend on its characteristics. The maximum and minimal $k_{\perp}$ probed will be a function of the bandwidth (of $\nu$) and of the baseline length $b$:

$$k_{\perp} = \frac{2\pi \nu b}{d_c(z) c}, \qquad (A1)$$

where $d_c(z)$ is the comoving distance at redshift $z$. On the other hand, the maximum and minimal $k_{\parallel}$ probed will depend on the bandwidth

<figure><img src="image_14.png" /><figcaption>**Figure A3.** Longitudinal (lower panel) and transversal (upper panel) $k$-modes observable by an interferometer, depending on its characteristics such as baseline length $b$ or frequency resolution. The bandwidth is fixed to 15 MHz. In the lower panel, the frequency resolutions of SKA-Low, HERA, and LOFAR are shown as reference and result in an upper limit on the $k_{\parallel}$ modes that can be observed by these instruments.</figcaption></figure>

B and the frequency resolution $\Delta\nu$ of the instrument:

$$ \left\{ \begin{aligned} k_{\parallel, \text{min}} &= \frac{2\pi}{\alpha(z)} \times \frac{1}{B}, \\ k_{\parallel, \text{max}} &= \frac{2\pi}{\alpha(z)} \times \frac{1}{2\Delta\nu}, \end{aligned} \right. \qquad (A2) $$

where $\alpha(z) \equiv c(1 + z)^2/[\nu_{21}H(z)]$. These limits are presented in Fig. A3 for various instruments, including SKA-Low, HERA, and LOFAR. In particular, for HERA, including the foreground wedge and additional buffer, we find that the accessible Fourier modes are limited to $0.15 \le k/[\text{Mpc}^{-1}] \le 2.05$. Because of the foregrounds, the instrument cannot reach the large-scale modes required to estimate the pkSZ power spectrum.

## APPENDIX B: FULL POSTERIOR DISTRIBUTIONS

In this section, we investigate the possible degeneracies between our model parameters. We choose the 3k2z model ($z = 6.5, 6.5, 7.8$ at $k = 0.1, 0.5, 0.5 \text{ Mpc}^{-1}$) outlined in Section 5.3 for which the Gelman–Rubin convergence criteria in Table 5 is met ($R - 1 \approx 0.001$). The
<figure><img src="image_14.png" /><figcaption>**Figure A1.** The pkSZ angular power spectrum obtained for the true electron power spectrum measured in the RSAGE const simulation, and for reconstructions based on the 21 cm power spectrum at various precision levels. The colour scheme is similar to Fig. 3. The results are compared to the SPT data point on the pkSZ power (Reichardt et al. 2021).</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure A2.** Physical scales contributing to the final pkSZ power at different angular multipoles $\ell$. The result for $\ell = 3000$ is highlighted in bold.</figcaption></figure>

see that, for all multipoles, most of the power stems from $10^{-3} < k/[\text{Mpc}^{-1}] < 1$ with only the upper limit increasing as $\ell$ increases.

The spherical Fourier modes used in Fig. A2 are effectively the sum of a transversal and a longitudinal mode $k = \sqrt{k_{\perp}^2 + k_{\parallel}^2}$. The Fourier modes probed by an interferometer measuring signal at a frequency $\nu$ corresponding to a redshift $z = v_{21}/\nu - 1$ depend on its characteristics. The maximum and minimal $k_{\perp}$ probed will be a function of the bandwidth (of $\nu$) and of the baseline length $b$:

$$k_{\perp} = \frac{2\pi \nu b}{d_c(z) c}, \qquad (A1)$$

where $d_c(z)$ is the comoving distance at redshift $z$. On the other hand, the maximum and minimal $k_{\parallel}$ probed will depend on the bandwidth

<figure><img src="image_16.png" /><figcaption>**Figure A3.** Longitudinal (lower panel) and transversal (upper panel) $k$-modes observable by an interferometer, depending on its characteristics such as baseline length $b$ or frequency resolution. The bandwidth is fixed to 15 MHz. In the lower panel, the frequency resolutions of SKA-Low, HERA, and LOFAR are shown as reference and result in an upper limit on the $k_{\parallel}$ modes that can be observed by these instruments.</figcaption></figure>

B and the frequency resolution $\Delta\nu$ of the instrument:

$$ \left\{ \begin{aligned} k_{\parallel, \text{min}} &= \frac{2\pi}{\alpha(z)} \times \frac{1}{B}, \\ k_{\parallel, \text{max}} &= \frac{2\pi}{\alpha(z)} \times \frac{1}{2\Delta\nu}, \end{aligned} \right. \qquad (A2) $$

where $\alpha(z) \equiv c(1 + z)^2/[v_{21}H(z)]$. These limits are presented in Fig. A3 for various instruments, including SKA-Low, HERA, and LOFAR. In particular, for HERA, including the foreground wedge and additional buffer, we find that the accessible Fourier modes are limited to $0.15 \le k/[\text{Mpc}^{-1}] \le 2.05$. Because of the foregrounds, the instrument cannot reach the large-scale modes required to estimate the pkSZ power spectrum.

## APPENDIX B: FULL POSTERIOR DISTRIBUTIONS

In this section, we investigate the possible degeneracies between our model parameters. We choose the 3k2z model ($z = 6.5, 6.5, 7.8$ at $k = 0.1, 0.5, 0.5 \text{ Mpc}^{-1}$) outlined in Section 5.3 for which the Gelman–Rubin convergence criteria in Table 5 is met ($R - 1 \approx 0.001$). The
<figure><img src="image_15.png" /><figcaption>**Figure B1.** One- and two-dimensional posterior probability distributions for our model parameters and one derived parameter (the optical depth $\tau$) when fitting a detection of the pKSZ power spectrum measurement at $\ell = 3000$ with a 10 per cent uncertainty and of the 21 cm power spectrum at redshifts of $z = 6.5, 6.5, 7.8$ at $k = 0.1, 0.50, 0.5 \text{ Mpc}^{-1}$ for 1000 h of integration with SKA-Low. For each parameter, the thin black vertical line is the true value used to generate the mock data and the dashed blue vertical line is the median of the distribution.</figcaption></figure>

corner plot in Fig. B1 clearly illustrates this, showcasing the relation between the global history and morphology parameters as well as the derived Thomson optical depth $\tau$. The parameter space is well sampled, with clear correlations seen only between either $z_{re}-z_{end}$ and $\log_{10}(\alpha_0)-\kappa$. A correlation is noticeable between the Thomson optical depth and each of the global history parameters. This is more than expected, as $\tau$ is an integral of the ionization history, defined as a function of $z_{re}$ and $z_{end}$ in equation (14). We discuss the constraints on $\tau$ for our different test cases in Appendix D in more detail.

## APPENDIX C: TESTING THE FORECAST WITH RSAGE

We have shown in Section 4.2 that with as few as three data points (one kSZ, two 21 cm), one can recover the global history of reionization and place some lower limits on the reionization morphology parameters. However, in these cases, both the mock data points and the sampled models were generated using our simplified derivation of the 21 cm power spectrum, ignoring third-order terms and cross-correlations, given in equation (7). In reality, of course, the data will include these extra terms. In this section, we will explore how this discrepancy can impact our forecast results.

We generate two new 21 cm data points corresponding to the 1k2z case of Section 4.2 for 100 h of observations with SKA-Low. This time, the points are generated using the full expression of $P_{21}(k, z)$, including the cross- and higher order terms, given in equation (6). As discussed in Section 4.1, these points have an amplitude between 10

and 25 per cent smaller than what was obtained with the simplified expression.

We fit our four reionization parameters to these new points (note that the kSZ data points remain unchanged) through MCMC sampling and show the resulting two-dimensional posterior distributions in Fig. C1. As expected from our discussions in Section 4.1, the general shape of the joint distributions is unchanged, only the results are now biased. As seen in Fig. 3, the simplified model tends to overestimate the 21 cm power for a given redshift, having an amplitude equivalent to a lower ionization level. For this reason, fitting it to unbiased data leads to an underestimate of the endpoint of reionization by $\Delta z_{end} = 0.12$. The mid-point is even more impacted as it is now biased by $\Delta z_{re} = 0.36$, nine times more than with the simplified mock data. Note that these modified reionization histories result in a large optical depth: $\tau = 0.070 \pm 0.002$. Hence, a way to partially mitigate these biases could be to impose a Gaussian prior on $\tau$. Similarly, the posterior distributions of the morphology parameters are shifted compared to the results with the simplified mock data points.

Although these results exhibit a strong bias in all parameters, this bias is nevertheless well understood and easy to model. A potential way of improvement would be to add a nuisance parameter to the fit, as a pre-factor to the model 21 cm power spectrum. This extra parameter would be marginalized over and would account for the extra power produced by the approximations done in equation (7). Another option would be to precisely quantify this bias and systematically subtract it from the models, allowing one to reproduce