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

Ivelin Georgiev1 ID , Adélie Gorce2,3 ID , and Garrelt Mellema1 ID
1The Oskar Klein Centre, Department of Astronomy, Stockholm University,
AlbaNova, SE-10691 Stockholm, Sweden
2Department of Physics and Trottier Space Institute, McGill University, 3600
University Street, Montreal, QC H3A 2T8, Canada
3Institut d’Astrophysique Spatiale, CNRS, Université Paris-Saclay, 91405
Orsay, France E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

During the Epoch of Reionization (EoR), the ultraviolet radiation from the
first stars and galaxies ionised the neutral hydrogen of the intergalactic
medium, which can emit radiation through its 21 cm hyperfine transition. This
21 cm signal is a direct probe of the first light sources in the early
Universe. 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 effect (pkSZ), observed as a foreground to the
primary cosmic microwave background temperature anisotropies on small scales.
Despite recent promising measurements by ground-based telescopes, placing
constraints on reionization from pkSZ observations is a non-trivial task,
subject to strong model dependence.

In this work, 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 are able to jointly fit mock 21 cm and pkSZ data points. We confirm
that these two 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
hours 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 global 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.

###### keywords:

observations – reionization – cosmological parameters

††pubyear: 2023††pagerange: Constraining cosmic reionization by combining the
kinetic Sunyaev-Zel’dovich and the 21 cm power spectra–D

##  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 ionised 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 onto dark matter (DM) haloes. The radiation produced by these
young stars and galaxies ionised the neutral hydrogen in the intergalactic
medium (IGM), forming H ii ‘bubbles’ which progressively grew and overlapped
as new ionising 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 recombination travel through the intergalactic
medium, they scatter off the free electrons produced by reionization,
resulting in a boost of the CMB polarisation and a suppression of temperature
anisotropies on large scales. Assuming an instantaneous reionization history,
measurements of $\tau=0.051\pm 0.006$ by the Planck Collaboration et al.
(2020) indicate that the midpoint of EoR lies around a redshift of
$z_{\mathrm{re}}\sim 8$. On the other hand, observations of the fluctuations
of the Lyman-$\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 ionising photons ($\lambda_{\mathrm{MFP}}$, Gaikwad et al., 2023) hint
that reionization may extend past redshift six and complete by
$z_{\mathrm{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 ionised IGM of the post-EoR Universe (Shaw et al., 2012), while
the ‘patchy’ kSZ (hereafter, pkSZ) effect has been shown to depend strongly on
the morphology of H ii bubbles during reionization (McQuinn et al., 2005;
Mesinger et al., 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 bubbles, respectively (Zahn et al., 2005;
Iliev et al., 2007; Gorce et al., 2020). By combining data from the South Pole
Telescope111https://pole.uchicago.edu (Ruhl et al., 2004) and the Planck
satellite222www.esa.int/planck (Planck Collaboration et al., 2020), George et
al. (2015) and, later, Reichardt et al. (2021) have constrained the amplitude
of the pkSZ angular power spectrum to
$D^{\mathrm{pkSZ}}_{\mathrm{3000}}=3.0\pm 1.0\,\mu\mathrm{K}^{2}$ using the
post-reionization models described in Shaw et al. (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 % 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\mathrm{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 et al., 2012). One way to get around this
computational challenge is using a parameterised 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 ionising
bubbles during reionization (Furlanetto et al., 2004, 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 ionising photons $\lambda_{\mathrm{MFP}}$ through an empirical
formula $k_{\mathrm{trans}}\approx 2/\lambda_{\mathrm{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)333www.lofar.org (Mertens et al.,
2020), the Hydrogen Epoch of Reionization Array
(HERA)444https://reionization.org (The HERA Collaboration et al., 2022; HERA
Collaboration et al., 2023), the Murchison Widefield Array
(MWA)555www.mwatelescope.org (Trott et al., 2020; Yoshiura et al., 2021), the
Giant Metrewave Radio Telescope (GMRT)666www.gmrt.ncra.tifr.res.in (Paciga et
al., 2013), as well as the forthcoming Square Kilometre Array
(SKA)777www.skatelescope.org (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 et al., 2022; HERA
Collaboration et al., 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., Liu et al., 2014; Chapman et al.,
2013; Mertens et al., 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. Indeed, the
kSZ is sensitive to slowly varying reionization histories, a section of
parameter space difficult to access for 21 cm observations because of the
brightness of foregrounds along such modes (Bégin et al., 2022). We place
constraints on cosmic reionization by utilising the fundamental relationship
between the power spectra of the 21 cm and the pkSZ signal from the EoR, which
we formalise through a simple yet effective parametric method.

This paper is organised 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 Monte-Carlo Markov Chain 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 construing it is on
the parameters of reionization. In section 5, we outline the caveats of the
forecast and explore special cases. We summarise our conclusion and elaborate
on future improvements in section 6.

##  2 Methods

In this section, we study the connection between the neutral and ionised
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 auto-correlation power spectrum. Using a parameterisation 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 H
i atom results in the emission or absorption of photons with a 21 cm
wavelength (Field, 1957). The abundance of H i 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 Epoch of Reionization 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

$\displaystyle\delta
T_{21}(\mathbf{r},z)=T_{0}(\mathbf{r},z)\,x_{\text{H\,{i}}~{}}(\mathbf{r},z)\left[1+\delta_{b}(\mathbf{r},z)\right]\
,$ (1)

with (Pritchard & Loeb, 2012)

$\displaystyle T_{0}(\mathbf{r},z)$ $\displaystyle\approx 27~{}\mathrm{mK}\
\bigg{(}\frac{1+z}{10}\bigg{)}^{1/2}\bigg{(}\frac{T_{\mathrm{s}}(\mathbf{r},z)-T_{\mathrm{CMB}}(z)}{T_{\mathrm{s}}(\mathbf{r},z)}\bigg{)}$
$\displaystyle\times\bigg{(}\frac{\Omega_{\mathrm{b}}}{0.044}\frac{h}{0.7}\bigg{)}\times\frac{\Omega_{\mathrm{b}}h\sqrt{0.15(1+z)}}{0.023\sqrt{10\Omega_{\mathrm{m}}}}.$
(2)

Here, $\delta_{b}=\rho_{b}/\bar{\rho}_{b}-1$ is the baryonic density
fluctuation and the bar denotes a spatial average and
$\delta_{\text{H\,{i}}~{}}$ is the fluctuation of the neutral fraction
$x_{\text{H\,{i}}~{}}$ field, i.e. the fraction of the IGM baryons that are
neutral. Naturally, $x_{\text{H\,{i}}~{}}=1-x_{\text{H\,{ii}}~{}}$, where
$x_{\text{H\,{ii}}~{}}$ is the ionisation field. The pre-factor $T_{0}$
includes most of the spatially averaged information, including cosmological
parameters, and the spin temperature information $T_{\mathrm{s}}$. We assume
the limit where during cosmic reionization the gas in the IGM is sufficiently
heated by the first ionising sources so $T_{\mathrm{s}}\gg T_{\mathrm{CMB}}$
and the spin-temperature dependence in $T_{0}$ is dropped.

Let us consider the electron fraction field

$x_{e}\equiv f_{\mathrm{H}}\,x_{\text{H\,{ii}}~{}}(1+\delta_{b}),$ (3)

where $f_{\mathrm{H}}$ is the fractional quantity of electrons per Hydrogen
atom. If the first reionization of Helium is considered, $f_{\mathrm{H}}\simeq
1.08$. We distinguish the mass-weighted $x_{m}$ and volume-weighted $x_{v}$
average ionised fractions of H i atoms. By definition, $\langle
x_{e}\rangle\equiv f_{\mathrm{H}}\,x_{m}$, and, if we take the fluctuations of
each element in equation (3), we have

$\displaystyle x_{m}(1+\delta_{e})$
$\displaystyle=x_{v}(1+\delta_{i})(1+\delta_{b}),$ (4)

where we refer to the density and ionisation field perturbations as ‘b’ and
‘i’, respectively. Including these definitions in equation (1) and using the
fact that $x_{\text{H\,{i}}~{}}=1-x_{\text{H\,{ii}}~{}}$, we have

$\displaystyle\frac{\delta T_{21}}{T_{0}(z)}$
$\displaystyle=(1+\delta_{b})-x_{m}(1+\delta_{e})$ (5)
$\displaystyle=(1+\delta_{b})\left[1-x_{v}(1+\delta_{i})\right],$

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

$\displaystyle\frac{P_{21}(k,z)}{T_{0}(z)^{2}}$
$\displaystyle=P_{bb}(k,z)+x_{m}(z)^{2}P_{ee}(k,z)$
$\displaystyle-2x_{v}(z)\Big{(}P_{bb}(k,z)+P_{bi}(k,z)+P_{bi,b}(k,z)\Big{)},$
(6)

where $P_{21}$ is the power spectrum of the 21 cm signal fluctuations,
$P_{ee}$ of the electron density fluctuations, $P_{ii}$ of the ionisation
fluctuations, $P_{bb}$ of the baryon over-density, 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 (2.1.1) as

$\displaystyle\frac{P_{21}(k,z)}{T_{0}(z)^{2}}=\left[1-2x_{v}(z)\right]P_{bb}(k,z)+x_{m}(z)^{2}P_{ee}(k,z).$
(7)

Some further simplifications of equation (7) are assumed: We identify the
mass-weighted $x_{m}$ and the volume-weighted $x_{v}$ ionised fractions and
approximate the baryon power spectrum as a biased dark matter power spectrum
$P_{bb}(k,z)=b_{\delta b}(k)P_{\delta\delta}(k,z)$, such that the fully
simplified expression used in this work, unless stated otherwise, is

$\displaystyle\frac{P_{21}(k,z)}{T_{0}(z)^{2}}=[1-2x_{v}(z)]\,b_{\delta
b}(k)P_{\delta\delta}(k,z)+x_{v}(z)^{2}P_{ee}(k,z).$ (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 the Limber
approximation and some assumptions on the velocity power spectrum, such that
(Mesinger et al., 2012; Gorce et al., 2020):

$\displaystyle C_{\ell}=$
$\displaystyle\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)\,{\rm
e}^{-2\tau(z)}\,\eta(z)\,\frac{\mathrm{d}\eta}{\mathrm{d}z}\,\mathrm{d}z,$ (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 $\mathbf{q}_{B,e}$ such that
$(2\pi)^{3}P_{B,e}\,\delta_{D}({\bf k}-{\bf
k}^{\prime})=\langle{\bf\tilde{q}}_{B,e}({\bf k})\
{\bf\tilde{q}}_{B,e}^{*}({\bf k}^{\prime})\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^{2}_{\mathrm{a,b}}(k)=k^{3}P_{\mathrm{a,b}}(k)/(2\pi^{2})$. We have

$\displaystyle\frac{\langle{\bf\tilde{q}}_{B,e}({\bf k})\
{\bf\tilde{q}}_{B,e}^{*}({\bf
k}^{\prime})\rangle}{(2\pi)^{3}{\delta_{D}}(|{\bf k}-{\bf
k}^{\prime}|)}\equiv$
$\displaystyle\frac{2\pi^{2}}{k^{3}}\Delta^{2}_{B,e}(k,z)$ (10)
$\displaystyle=$
$\displaystyle\frac{1}{(2\pi)^{3}}\int\mathrm{d}^{3}k^{\prime}\,\left[(1-\mu^{2})\,P_{ee}(|{\bf
k}-{\bf k}^{\prime}|)\,P_{vv}(k^{\prime})\right.$
$\displaystyle\left.-\frac{(1-\mu^{2})\,k^{\prime}}{|{\bf k}-{\bf
k}^{\prime}|}P_{ev}(|{\bf k}-{\bf k}^{\prime}|)\,P_{ev}(k^{\prime})\,\right],$

where where $\mu={\bf\hat{k}}\cdot{\bf\hat{k}}^{\prime}$ 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_{ve}(k,z)=\frac{f\dot{a}(z)}{k}b_{\delta
e}(k,z)P_{\delta\delta}^{\operatorname{lin}}(k,z)$ (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)$ parameterisation introduced in Gorce et al.
(2020), that is

$\displaystyle P_{ee}(k,z)=$ $\displaystyle
f_{\mathrm{H}}\left[1-x_{v}(z)\right]\times\frac{\alpha_{0}x_{v}(z)^{-1/5}}{1+[k/\kappa]^{3}x_{v}(z)}$
(12) $\displaystyle+x_{v}(z)b_{\delta b}(k,z)^{2}P_{\delta\delta}(k,z),$

where $\kappa$ is the electron drop-off frequency, which can be related to the
typical size of ionised 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-dark matter bias
$b_{\delta b}(k)$ is given by the adapted Shaw et al. (2012) parameterisation,
that is

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

where $k_{f}=9.4\,\mathrm{Mpc}^{-1}$ and $g=0.5$ are constant with
redshift.888In this work we address comoving megaparsec as Mpc. The global
reionization history $x_{v}(z)$ in equation (12) is defined according to the
following parameterisation (Douspis et al., 2015; Gorce et al., 2018):

$\displaystyle z$ $\displaystyle<z_{\mathrm{end}},\quad x_{v}(z)=1,$ (14)
$\displaystyle z$ $\displaystyle\geqslant z_{\mathrm{end}},\quad
x_{v}(z)=\bigg{(}\frac{z_{\mathrm{early}}-z}{z_{\mathrm{early}}-z_{\mathrm{end}}}\bigg{)}^{\alpha},$
(15)

where $z_{\mathrm{early}}=20$ is the redshift for which
$x_{v}(z_{\mathrm{early}})\approx 10^{-4}$, the ionisation leftover from
recombination. Random forests are used to speed up computing times and predict
the pkSZ power spectrum given the parameter set; see Gorce et al. (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 CAMB999Available at
https://github.com/cmbant/CAMB. (Lewis et al., 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 Monte-Carlo Markov chain 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
Monte-Carlo Markov Chain sampler (MCMC, Lewis & Bridle, 2002; Lewis, 2013),
modified as described in Gorce et al. (2022). Rather than sampling the Thomson
optical depth $\tau$, we fit for the reionization mid- and endpoint,
$z_{\mathrm{re}}$ and $z_{\mathrm{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 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}_{\mathrm{tot}}=-\frac{(\Delta_{21,\mathrm{true}}^{2}-\Delta_{21,\mathrm{model}}^{2})^{2}}{2(\Delta^{2}_{\mathrm{noise}}+\Delta^{2}_{\mathrm{var}})}-\frac{(C_{\ell,\mathrm{true}}^{\mathrm{pkSZ}}-C_{\ell,\mathrm{model}}^{\mathrm{pkSZ}})^{2}}{2\sigma^{2}_{\mathrm{pkSZ}}},$
(16)

where uncertainty in the measurements $\Delta^{2}_{\mathrm{noise}}$,
$\Delta^{2}_{\mathrm{var}}$, and $\sigma_{\mathrm{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 asses 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 marginalised posterior
probability’s maximum, which is more suitable for skewed distributions.
Confidence intervals correspond to intervals with the highest probability
density at 68 %.

Parameter | Ref. value | Prior
---|---|---
$h$ | 0.681 | N.A.
$\Omega_{b}$ | 0.050 | N.A.
$\Omega_{m}$ | 0.302 | N.A.
$n_{s}$ | 0.96 | N.A.
$A_{s}$ | $2.10\times 10^{-9}$ | N.A.
$z_{\mathrm{re}}$ | 7.37 | [5.0, 10.0]
$z_{\mathrm{end}}$ | 6.15 | [4.5, 9.]
$dz=z_{\mathrm{re}}-z_{\mathrm{end}}$ | 1.22 | [0.5, 5.5]
$\log[\alpha_{0}/\mathrm{Mpc}^{3}]$ | 3.12 | [2.5, 4.5]
$\kappa/\mathrm{Mpc}^{-1}$ | 0.145 | [0.05, 0.25]
Table 1: Reference cosmological and reionization parameters used to generate
the mock data, based on the RSAGE simulation. All parameters apart from
$z_{\mathrm{re}}$, $z_{\mathrm{end}}$, log${}_{10}(\alpha_{0})$, and $\kappa$
are fixed.

####  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^{2}_{\mathrm{noise}}(k,z)\equiv
k^{3/2}\lambda_{21}(z)\times\frac{2}{\pi}\sqrt{\frac{D^{2}_{c}(z)\,\Delta
D_{c}}{A_{\mathrm{eff}}}}\frac{T_{\mathrm{sys}}(z)^{2}}{B\,t_{\mathrm{int}}}\frac{A_{\mathrm{core}}A_{\mathrm{eff}}}{A_{\mathrm{coll}}^{2}},$
(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_{\mathrm{sys}}=100+300(\nu_{\mathrm{obs}}/\nu_{0})^{-2.55}$ K, for
$\nu_{0}=150$ MHz and $\nu_{\mathrm{obs}}(z)=c/\lambda_{21}(z)$ is the
observed 21 cm frequency. The total integration time is $t_{\mathrm{int}}$,
$D_{c}(z)$ is the comoving distance to redshift $z$, and $n_{\mathrm{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_{\mathrm{coll}}$, such that $A_{\mathrm{coll}}\equiv N_{a}\pi R_{a}^{2}$
where $R_{a}$ is the radius of the antenna. We have $A_{\mathrm{core}}$ the
core area of the array and $A_{\mathrm{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.

Parameters | SKA-Low | MWA | LOFAR | HERA
---|---|---|---|---
$B$, MHz | 10 | 10 | 10 | 10
$N_{\mathrm{ant}}$ | 224 | 128 | 12 | 350
$A_{\mathrm{eff}}$, m2 | 600 | 21.5 | 804 | 150
$R_{\mathrm{core}}$, m | 500 | 150 | 150 | 150
$\theta$, deg2 | 32 | 24.72 | 3.72 | 82
Table 2: Parameters describing the mock observations by different telescopes
used to derive our 21 cm power spectrum errors (equation 17).

We also include the contribution of sample variance to the uncertainty of the
measurement following the relation (see eqs. 9 &10 of Mellema et al., 2013):

$\Delta^{2}_{\mathrm{var}}(k,z)=0.01\,\bigg{(}\frac{k}{0.1\,\mathrm{Mpc}^{-1}}\bigg{)}^{-3/2}\bigg{(}{\frac{V}{1\,\mathrm{cGpc}^{3}}}\bigg{)}^{-1/2}\bigg{(}{\frac{\epsilon}{0.5}}\bigg{)}^{-1/2}\Delta^{2}_{21},$
(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\,\mathrm{Gpc}^{3}\times\bigg{(}\frac{\theta}{5^{\circ}}\bigg{)}^{2}\bigg{(}\frac{B}{12\,\mathrm{MHz}}\bigg{)}\,\big{[}(1+z)^{1/2}-2\big{]},$
(19)

and $\theta$ is the field of view.

Figure 1: Measurement errors on the 21 cm power spectrum (seen in black) for
1000 hours 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.

An example of measurement errors for an integration time of 1000 hours 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\lessapprox 0.3$ Mpc-1
because of its strong $k$-dependency ($\Delta^{2}_{\mathrm{var}}\propto
k^{-3/2}$). The opposite is true of the noise estimate for MWA, which has a
larger field of view and as $\Delta^{2}_{\mathrm{var}}\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 ten per cent uncertainty on the data point throughout this work. We base
this on the current constraint of the total kSZ amplitude
$D^{\mathrm{pkSZ}}_{\mathrm{3000}}=3.0\pm 1.0\,\mu\mathrm{K}^{2}$ (Reichardt
et al., 2021), assuming that in future observations longer observation times
will allow to reduce the noise, whilst numerous frequency channels and
improved modelling will help characterising and removing other CMB foregrounds
(Maniyar et al., 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 parameterisation
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}$ dark matter particles and a volume of (160 Mpc)3.
The radiative transfer is conducted using the semi-numerical 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 ionising photons $f_{\mathrm{esc}}=0.2$. The cosmology used is consistent
with Planck Collaboration et al. (2020) and summarised in Table 1. The mid-
and endpoint for the EoR are $z=7.37,6.15$, respectively. In Fig. 2, we show
simulation slices of the differential surface brightness temperature of the 21
cm signal.

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 ionisation 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.

##  4 Results

In this section, we first look at the uncertainties associated with our
simplified derivation of the 21 cm power spectrum given an electron power
spectrum by comparing it to the full spectrum computed from the RSAGE
simulation. In Sec. 4.2, we then apply this derivation to jointly fit mock 21
cm and kSZ data points to a reionization model.

###  4.1 Reconstructing the power spectra

We want to understand the role of each contributing term in equation (2.1.1),
especially the ones we have discarded in our simplified $P_{21}$
reconstruction (equation 7), 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:

-

The recons no cross case corresponds to the highest simplification level,
corresponding to equation (7), and the approximation we will use in the
remaining of this work. The contributions from both third-order terms and
cross-terms are discarded.

-

For recons with cross, we add the cross power spectra contributions:
$-2x_{v}(z)P_{bi}(k,z)$101010Note that this contribution can be positive or
negative, as illustrated in Fig. 3.,

-

For recons with third order, we further add the three-point power spectrum
contribution: $-2x_{v}(z)P_{bi,b}(k,z)$10. This corresponds to the full form
of the 21 cm power spectrum, without simplification.

Figure 3:  Upper Panel: The dimensionless 21 cm power spectrum at
$k=0.141~{}\mathrm{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
text for details). recons no cross corresponds to equation (7), which is used
for the forecast in section 4.2. 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 (2.1.1). Middle Panel: Ratio of each component in
equation (2.1.1) 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.

We look at the relative contribution of each of the terms in equation (2.1.1)
to the 21 cm power spectrum in the RSAGE simulation. This is represented in
the middle panel of Fig. 3 for $k=0.141~{}\mathrm{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
ionisation level ($x_{v}\lesssim 10^{-4}$), the ionisation 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\%$, when the first luminous sources have
reionised their local over-densities and the ionisation field is correlated
with the density field, $P_{21}$ is primarily determined by the difference
between the electron density power spectrum and the $P_{ib}$ 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 ionisation 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 (7), in blue, while the full model presented
in equation (2.1.1), 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. 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
ionisation 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 (2.1.1), 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, whilst keeping its
limitations in mind. We use this model to derive our true data points and each
sampled model. 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 Sec. 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 midpoint, 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 et al., 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, whilst 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 hrs
of SKA-Low observations before turning to mock observations by MWA and LOFAR.

####  4.2.1 Forecast with SKA-Low

Figure 4: 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$ Mpc-1 for either 100 and 1000 hours 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.

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 et al., 2014), considering noise levels associated with
100 and 1000 hours 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^{pkSZ}_{l=3000}=0.86\pm 0.09\,\mu 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.

In Fig. 4, we show the joint probability distributions for our four sampled
parameters (see appendix B for a discussion on parameter degeneracies). The
best-fit parameters, their 1$\sigma$ error, and Gelman-Rubin convergence
parameters are presented in Table 3111111The 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 marginalise 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..

Our results indicate that we can recover the true values of $z_{\mathrm{re}}$
and $z_{\mathrm{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 Sec. 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 ionising sources (marginalising 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 midpoint of the EoR, the 21 cm power spectrum
decreases in amplitude.

For both cases presented in Fig. 4, we can place loose lower limits on
log${}_{10}(\alpha_{0})$, although the 1000 hrs 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.

Models |  |  | $z_{\mathrm{re}}$ | $z_{\mathrm{end}}$ | log${}_{10}(\alpha_{0})$ | $\kappa$ | $\tau$ | d$z$
---|---|---|---|---|---|---|---|---
Label | Data | True | 7.37 | 6.15 | 3.12 | 0.16 | 0.0649 | 1.22
1k2z | z = 6.5, 7.8 | 1000h | 7.39 $\pm$ 0.14 | 6.15 $\pm$ 0.04 | 3.04 $\pm$ 0.32 | 0.18 $\pm$ 0.04 | 0.0651 $\pm$ 0.0018 | 1.24 $\pm$ 0.17
| k = 0.5 Mpc-1 | R - 1 | 0.3 | 1.5 | 1.6 | 0.5 | N/A | N/A
1k2z | z = 6.5, 7.8 | 100h | 7.42 $\pm$ 0.13 | 6.15 $\pm$ 0.04 | 2.95 $\pm$ 0.31 | 0.19 $\pm$ 0.04 | 0.0655 $\pm$ 0.0017 | 1.27 $\pm$ 0.16
| k = 0.5 Mpc-1 | R - 1 | 0.1 | 0.1 | 0.3 | 0.3 | N/A | N/A
Table 3: Best fit values of the model parameters for the 1k2z cases, their
corresponding 1$\sigma$ uncertainty and Gelman–Rubin parameter. Best-fit
values and errors are also given for derived parameters such as $\tau$.

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. 5 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_{\mathrm{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_{\mathrm{end}}$ decreases the duration of the EoR, suppressing the pkSZ
amplitude. This explains the strong anti-correlation observed between the two
parameters in Fig. 4. Note that for the model where $z_{\mathrm{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. 5, 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 midpoint), whilst the $P_{ee}$ power breaks
as $k^{-3}$ for scales $k>\kappa$, such that $\kappa$ can be interpreted as
the minimal size of ionised regions during reionization. Hence, varying
log${}_{10}(\alpha_{0})$ directly impacts the amplitude of both observables,
whilst 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. 4. 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.

Figure 5: 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
correspoding colour.

####  4.2.2 Constraining power of each dataset

Figure 6:  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 hours of SKA-Low observations. The
vertical and horizontal lines show the ‘true’ values used to generate the mock
data.

In this subsection, we examine each dataset’s ability to individually impose
constraints on reionization. In Fig. 6, we present the joint posterior
distribution of the sampled parameters for three different cases assuming a
SKA-Low detection after integrating over 1000 hours. We compare a case where
the forecast is run only with the 21 cm data points (o21), a case with the
pkSZ measurement only (opkSZ), and, lastly, the joint forecast previously
illustrated in Fig. 4. 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_{\mathrm{re}}$ and $z_{\mathrm{end}}$, then
the astrophysical log${}_{10}(\alpha_{0})$ and $\kappa$, in Fig 6.

Examining the two-dimensional $z_{\mathrm{re}}$ and $z_{\mathrm{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 $\mathrm{d}z\equiv
z_{\mathrm{re}}-z_{\mathrm{end}}$ as the rate of the growth of ionising
bubbles impacts the amplitude of the signal (see also the first and second
column of Fig. 5). Models with an extended reionization period result in a
larger amplitude of the pkSZ power spectrum as more free electrons interact
with the ionised medium, and reciprocally (McQuinn et al., 2005; Mesinger et
al., 2012; Battaglia et al., 2013). Note that the posterior distribution for
the opkSZ case covers a region of parameter space where the Universe has
reionised 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 o21 case, the data constraints tightly constrain the end of
reionization but only place an upper limit on the midpoint. Compared to the
pkSZ-only case, the o21 case is not sensitive to the duration of reionization
and favours extended models of the EoR, as previously noted by Bégin et al.
(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_{\mathrm{end}}$.
Lastly, models at the lower left of the two-dimensional posterior
($z_{\mathrm{re}}\lesssim 7$ and $z_{\mathrm{end}}\lesssim 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 anti-correlated and the pkSZ data alone only provides a biased
measurement. This can be better understood by examining the two lower rows of
Fig. 5: 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 ionisation 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 o21 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).

To summarise, we find that with as few as 100 hours of integration with SKA-
Low, we can successfully constrain the global reionization parameters
$z_{\mathrm{re}}$ and $z_{\mathrm{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_{\mathrm{end}}$, it cannot
constrain $z_{\mathrm{re}}$. Conversely, we find that the pkSZ measurement
mainly constrains the duration of reionization
(d$z=z_{\mathrm{re}}-z_{\mathrm{end}}$). Naturally, the combined analysis
benefits from the complimentarity of probes. While 100 hours of integration
are sufficient to obtain upper limits on the reionization morphology
parameters, it is only with 1000 hours 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 & 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.

###  5.1 Importance of the choice of data points

####  5.1.1 Choice of redshift

Figure 7: 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 hours of SKA-Low observations compared to the fiducial case from Sec. 4.2
(in blue). The vertical and horizontal lines show the ‘true’ values used to
generate the mock data.

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 et al., 2022). The resulting two-dimensional posterior
distributions for $z_{\mathrm{re}}-z_{\mathrm{end}}$ and
$\log\alpha_{0}-\kappa$ are presented in Fig. 7 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 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. 7, 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

Figure 8:  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 hours 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. Models |  |  | $z_{\mathrm{re}}$ | $z_{\mathrm{end}}$ | log${}_{10}(\alpha_{0})$ | $\kappa$ | $\tau$ | dz
---|---|---|---|---|---|---|---|---
Label | Data | True | 7.37 | 6.15 | 3.12 | 0.16 | 0.0649 | 1.22
MWA | z = 6.5, 7.8 | 1000h | 7.02 $\pm$ 0.49 | 5.39 $\pm$ 0.54 | 3.25 $\pm$ 0.45 | 0.13 $\pm$ 0.04 | 0.0623 $\pm$ 0.0051 | 1.53 $\pm$ 0.49
| k = 0.1 Mpc-1 | R - 1 | 0.01 | 0.008 | 0.03 | 0.03 | N/A | N/A
LOFAR | z = 8.3, 9.1 | 1000h | 7.30 $\pm$ 0.87 | 5.74 $\pm$ 0.80 | 3.15 $\pm$ 0.33 | 0.14 $\pm$ 0.03 | 0.0643 $\pm$ 0.0094 | 1.38 $\pm$ 0.83
| k = 0.1 Mpc-1 | R - 1 | 0.06 | 0.03 | 0.1 | 0.04 | N/A | N/A
Table 4: Best fit values of the distributions in Fig. 8 and their
corresponding 1$\sigma$ uncertainty as well as the Gelman–Rubin convergence
diagnostic for each parameter.

Figure 9: 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 hours of integration time with SKA, seen in brown and compared to the 1k2z 1000h fiducial case shown in blue (see Fig, 4). Models |  |  | $z_{\mathrm{re}}$ | $z_{\mathrm{end}}$ | log${}_{10}(\alpha_{0})$ | $\kappa$ | $\tau$ | dz
---|---|---|---|---|---|---|---|---
Label | Data | True | 7.37 | 6.15 | 3.12 | 0.16 | 0.0649 | 1.22
2k1z | z = 6.5, 6.5 | 1000h | 7.22 $\pm$ 0.21 | 5.29 $\pm$ 0.55 | 2.97 $\pm$ 0.08 | 0.15 $\pm$ 0.01 | 0.0649 $\pm$ 0.0016 | 1.90 $\pm$ 0.40
hiz | z = 7.8, 10.4 | 1000h | 7.28 $\pm$ 0.33 | 6.02 $\pm$ 0.52 | 3.09 $\pm$ 0.15 | 0.17 $\pm$ 0.02 | 0.0642 $\pm$ 0.0025 | 1.27 $\pm$ 0.19
| k = 0.1, 0.1 Mpc-1 | R - 1 | 0.7 | 0.7 | 0.5 | 0.5 | N/A | N/A
3k2z | z = 6.5, 6.5, 7.8 | 1000h | 7.37 $\pm$ 0.07 | 6.15 $\pm$ 0.03 | 3.12 $\pm$ 0.01 | 0.164 $\pm$ 0.001 | 0.0649 $\pm$ 0.0010 | 1.22 $\pm$ 0.10
| k = 0.1, 0.5, 0.1 Mpc-1 | R - 1 | 0.02 | 0.02 | 0.005 | 0.004 | N/A | N/A
Table 5: Best fit values of the distributions in and their corresponding
1$\sigma$ uncertainty as well as the Gelman–Rubin convergence diagnostic for
each parameter.

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 hours 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 utilised the fact that
reionization is an evolving process, whilst here, we make use of its multi-
scale 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-
ionised 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’ which encodes information about the IGM on small scales.

Fig. 7 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_{\mathrm{re}}$
and $z_{\mathrm{end}}$ values, the two-dimensional posterior distribution
exhibits a tail and is biased towards late reionization models
($z_{\mathrm{end}}=5.29\pm 0.55$, whilst the true value is at redshift 6.15).
Observing a non-zero 21 cm data point close to the true $z_{\mathrm{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 midpoint of reionization would either re-ionise 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. 5 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_{\mathrm{re}}$-$z_{\mathrm{end}}$ posterior as well as biased morphology
parameters. One way to improve constraints on $z_{\mathrm{end}}$ would be to
consider 21 cm data for $k>1.0$ Mpc-1 scales, or to limit the prior on
$z_{\mathrm{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 considered121212Note that we do
not consider HERA in this section, since its ideal noise levels are lower than
the ones for the SKA (Fig. 1). 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$ 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 hours. The resulting two-
dimensional posterior probability distributions of the cosmological and
morphology parameters are presented in Fig. 8.

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. 8, we are still able to
constrain $z_{\mathrm{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_{\mathrm{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_{\mathrm{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
Sec. 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 data point further reduces the ability of the forecast
to place upper limits on the global reionization parameters (see Sec. 5.1).
For example, as seen in grey in the upper panel of Fig. 8, we can constrain
$z_{\mathrm{re}}<8.0$ and can only exclude only late reionization models for
which $z_{\mathrm{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 hours (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 et al. (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 could obtain from early,
partial 21 cm data, as a detection will only be achieved incrementally. On the
one hand, with with two data points at 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_{\mathrm{re}}$ and $z_{\mathrm{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 multi-scale regime (3k2z case, seen in brown in Fig. 9)
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.1$ Mpc-1 for 1000
hours 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. 9), we
improve the accuracy of our constraints of the reionization parameters: We get
$z_{\mathrm{re}}=7.37\pm 0.07$ and $z_{\mathrm{end}}=6.15\pm 0.03$,
significantly reducing the 1$\sigma$ standard deviation (especially for
$z_{\mathrm{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 polarisation anisotropies (Planck
Collaboration et al., 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 et al., 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.
7 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 parameterisation of the electron power spectrum
$P_{ee}(k,z)$ presented in Gorce et al. (2020). The resulting equation
(2.1.1), while linking both observables, contains non-trivial cross-
correlation and second-order terms of the ionised and density fields. We
choose to ignore such terms in this work, leading to equation (7). 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_{\mathrm{re}}$ and $z_{\mathrm{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 hours 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 hours 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. 4 and throughout this work), we are able to constrain the global
reionization parameters $z_{\mathrm{re}}=7.42\pm 0.13$ and
$z_{\mathrm{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
d$z=z_{\mathrm{re}}-z_{\mathrm{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 ionisation 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_{\mathrm{re}}$ and $z_{\mathrm{end}}$ are loosened compared to the 1k2z
case, at $z=6.5$ and 7.8 (Fig. 7). 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 (2k1z,
$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. 7, in yellow) but
constraints on the morphology parameters are improved. However, the recovered
values are biased by 6%. 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
multi-scale 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. 9, 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_{\mathrm{re}}=7.37\pm 0.07$ & $z_{\mathrm{end}}=6.15\pm 0.03$ and the
morphology parameters log${}_{10}(\alpha_{0})=3.12\pm 0.01$ & $\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. (Trott et al.,
2020; Patil et al., 2017; Mertens et al., 2020). Naturally, we find the mock
data less constraining than in the case of 1k2z (see Fig. 8 in gray 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_{\mathrm{re}}\leq 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_{\mathrm{end}}\geq 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 (7). Potential future developments would be develop
analytic models of the cross-power spectrum and include them in the derivation
(McQuinn et al., 2005; Schneider et al., 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).

## Acknowledgements

The authors thank Adrian Liu, Jack Line, Leon Koopmans, and Cathryn Trott for
useful discussions about different parts of this work. IG acknowledges the
support form the Gustaf och Ellen Kobbs stipendiestiftelse and the Alva and
Lennart Dahlmark Research Grant. AG’s work is partly 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 Canada131313See http://www.alliancecan.ca
and https://www.calculquebec.ca.. Computations were notably performed on the
Cedar supercomputer, hosted by the Simon Fraser University.

## Data Availability

The sampling code used to obtain all the results presented in this article is
publicly available at https://github.com/adeliegorce/forecast_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 (2016) Alvarez M. A., 2016, ApJ, 824, 118
  * Battaglia et al. (2013) Battaglia N., Trac H., Cen R., Loeb A., 2013, ApJ, 776, 81
  * Bégin et al. (2022) Bégin J.-M., Liu A., Gorce A., 2022, Phys. Rev. D, 105, 083503
  * Bosman et al. (2022) Bosman S. E. I., et al., 2022, MNRAS, 514, 55
  * Chapman et al. (2013) Chapman E., et al., 2013, MNRAS, 429, 165
  * Chen et al. (2023) Chen N., Trac H., Mukherjee S., Cen R., 2023, ApJ, 943, 138
  * Douspis et al. (2015) Douspis M., Aghanim N., Ilić S., Langer M., 2015, A&A, 580, L4
  * Douspis et al. (2022) Douspis M., Salvati L., Gorce A., Aghanim N., 2022, A&A, 659, A99
  * Field (1957) Field G. B., 1957, AJ, 62, 15
  * Foreman-Mackey (2016) Foreman-Mackey D., 2016, The Journal of Open Source Software, 1, 24
  * Foreman-Mackey et al. (2013) Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306
  * Furlanetto et al. (2004) Furlanetto S. R., Zaldarriaga M., Hernquist L., 2004, ApJ, 613, 1
  * Furlanetto et al. (2006) Furlanetto S. R., McQuinn M., Hernquist L., 2006, MNRAS, 365, 115
  * Gaikwad et al. (2023) Gaikwad P., et al., 2023, MNRAS, 525, 4093
  * Gelman & Rubin (1992) Gelman A., Rubin D. B., 1992, Statistical Science, 7, 457
  * George et al. (2015) George E. M., et al., 2015, ApJ, 799, 177
  * Georgiev et al. (2022) Georgiev I., Mellema G., Giri S. K., Mondal R., 2022, MNRAS, 513, 5109
  * Gorce et al. (2018) Gorce A., Douspis M., Aghanim N., Langer M., 2018, A&A, 616, A113
  * Gorce et al. (2020) Gorce A., Ilić S., Douspis M., Aubert D., Langer M., 2020, A&A, 640, A90
  * Gorce et al. (2022) Gorce A., Douspis M., Salvati L., 2022, A&A, 662, A122
  * HERA Collaboration et al. (2023) HERA Collaboration et al., 2023, ApJ, 945, 124
  * Howlett et al. (2012) Howlett C., Lewis A., Hall A., Challinor A., 2012, J. Cosmology Astropart. Phys., 1204, 027
  * Hunter (2007) Hunter J. D., 2007, Computing In Science & Engineering, 9, 90
  * Hutter (2018) Hutter A., 2018, MNRAS, 477, 1549
  * Iliev et al. (2007) Iliev I. T., Pen U.-L., Bond J. R., Mellema G., Shapiro P. R., 2007, ApJ, 660, 933
  * Iliev et al. (2014) Iliev I. T., Mellema G., Ahn K., Shapiro P. R., Mao Y., Pen U.-L., 2014, MNRAS, 439, 725
  * Jones et al. (2001) Jones E., Oliphant T., Peterson P., et al., 2001, SciPy: Open source scientific tools for Python, http://www.scipy.org/
  * Koopmans et al. (2015) Koopmans L., et al., 2015, in Advancing Astrophysics with the Square Kilometre Array (AASKA14). p. 1 (arXiv:1505.07568), doi:10.22323/1.215.0001
  * Lewis (2013) Lewis A., 2013, Phys. Rev. D, 87, 103529
  * Lewis & Bridle (2002) Lewis A., Bridle S., 2002, Phys. Rev. D, 66, 103511
  * Lewis et al. (2000) Lewis A., Challinor A., Lasenby A., 2000, ApJ, 538, 473
  * Lidz et al. (2007) Lidz A., Zahn O., McQuinn M., Zaldarriaga M., Dutta S., Hernquist L., 2007, ApJ, 659, 865
  * Liu et al. (2014) Liu A., Parsons A. R., Trott C. M., 2014, Phys. Rev. D, 90, 023018
  * Liu et al. (2016) Liu A., Pritchard J. R., Allison R., Parsons A. R., Seljak U., Sherwin B. D., 2016, Phys. Rev. D, 93, 043013
  * Maniyar et al. (2021) Maniyar A., Béthermin M., Lagache G., 2021, A&A, 645, A40
  * McQuinn et al. (2005) McQuinn M., Furlanetto S. R., Hernquist L., Zahn O., Zaldarriaga M., 2005, ApJ, 630, 643
  * McQuinn et al. (2006) McQuinn M., Zahn O., Zaldarriaga M., Hernquist L., Furlanetto S. R., 2006, ApJ, 653, 815
  * Mellema et al. (2013) Mellema G., et al., 2013, Experimental Astronomy, 36, 235
  * Mertens et al. (2018) Mertens F. G., Ghosh A., Koopmans L. V. E., 2018, MNRAS, 478, 3640
  * Mertens et al. (2020) Mertens F. G., et al., 2020, MNRAS, 493, 1662
  * Mesinger et al. (2012) Mesinger A., McQuinn M., Spergel D. N., 2012, MNRAS, 422, 1403
  * Oliphant (2006) Oliphant T., 2006, NumPy: A guide to NumPy, USA: Trelgol Publishing, http://www.numpy.org/
  * Paciga et al. (2013) Paciga G., et al., 2013, MNRAS, 433, 639
  * Park et al. (2013) Park H., Shapiro P. R., Komatsu E., Iliev I. T., Ahn K., Mellema G., 2013, ApJ, 769, 93
  * Patil et al. (2017) Patil A. H., et al., 2017, ApJ, 838, 65
  * Planck Collaboration et al. (2016) Planck Collaboration et al., 2016, A&A, 596, A108
  * Planck Collaboration et al. (2020) Planck Collaboration et al., 2020, A&A, 641, A1
  * Pritchard & Loeb (2012) Pritchard J. R., Loeb A., 2012, Reports on Progress in Physics, 75, 086901
  * Raghunathan & Omori (2023) Raghunathan S., Omori Y., 2023, ApJ, 954, 83
  * Raste et al. (2021) Raste J., Kulkarni G., Keating L. C., Haehnelt M. G., Chardin J., Aubert D., 2021, MNRAS, 507, 4684
  * Reichardt et al. (2021) Reichardt C. L., et al., 2021, ApJ, 908, 199
  * Ruhl et al. (2004) Ruhl J., et al., 2004, in Bradford C. M., et al., eds, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series Vol. 5498, Z-Spec: a broadband millimeter-wave grating spectrometer: design, construction, and first cryogenic measurements. pp 11–29 (arXiv:astro-ph/0411122), doi:10.1117/12.552473
  * Schneider et al. (2023) Schneider A., Schaeffer T., Giri S. K., 2023, Phys. Rev. D, 108, 043030
  * Seiler et al. (2018) Seiler J., Hutter A., Sinha M., Croton D., 2018, MNRAS, 480, L33
  * Seiler et al. (2019) Seiler J., Hutter A., Sinha M., Croton D., 2019, MNRAS, 487, 5739
  * Shaw et al. (2012) Shaw L. D., Rudd D. H., Nagai D., 2012, ApJ, 756, 15
  * Sunyaev & Zeldovich (1980) Sunyaev R. A., Zeldovich I. B., 1980, ARA&A, 18, 537
  * The HERA Collaboration et al. (2022) The HERA Collaboration et al., 2022, ApJ, 925, 221
  * Trott et al. (2020) Trott C. M., et al., 2020, MNRAS, 493, 4711
  * Wise (2019) Wise J. H., 2019, arXiv e-prints, p. arXiv:1907.06653
  * Yoshiura et al. (2021) Yoshiura S., et al., 2021, MNRAS, 505, 4775
  * Zahn et al. (2005) Zahn O., Zaldarriaga M., Hernquist L., McQuinn M., 2005, ApJ, 630, 657
  * Zahn et al. (2012) 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 (2.1.1). 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. 10,
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 simulation141414Remember that all cosmological and reionization
parameters are derived from the RSAGE simulation.. 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 (2.1.1)
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\leq\ell<10\,000$ and show the result in Fig. 11. We see
that, for all multipoles, most of the power stems from
$10^{-3}<k/[\mathrm{Mpc}^{-1}]<1$ with only the upper limit increasing as
$\ell$ increases.

Figure 10: 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).  Figure 11: Physical scales
contributing to the final pkSZ power at different angular multipoles $\ell$.
The result for $\ell=3000$ is highlighted in bold.

The spherical Fourier modes used in Fig. 11 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}{d_{c}(z)}\frac{\nu b}{c},$ (20)

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
$B$ and the frequency resolution $\Delta\nu$ of the instrument:

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

where $\alpha(z)\equiv c(1+z)^{2}/[\nu_{21}H(z)]$. These limits are presented
in Fig. 12 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\leq
k/[\mathrm{Mpc}^{-1}]\leq 2.05$. Because of the foregrounds, the instrument
cannot reach the large-scale modes required to estimate the pkSZ power
spectrum.

Figure 12: 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 and bandwidth. 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.

##  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$
Mpc-1) outlined in Sec. 5.3 for which the Gelman-Rubin convergence criteria in
Table 5 is met ($R-1\approx 0.001$). The corner plot in Fig. 13 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_{\mathrm{re}}-z_{\mathrm{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 ionisation history, defined as a function of $z_{\mathrm{re}}$
and $z_{\mathrm{end}}$ in equation (14). We discuss the constraints on $\tau$
for our different test cases in appendix. D in more detail.

Figure 13: 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 ten 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$ Mpc-1 for 1000 hours 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.

##  Appendix C Testing the forecast with RSAGE

We have shown in Sec. 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 Sec.
4.2, for 100 hours 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 (2.1.1). As discussed in Sec. 4.1, these
points have an amplitude between 10 to 25% 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. 14. As expected from our
discussions in Sec. 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 ionisation level. For this reason, fitting it
to unbiased data leads to an underestimate of the endpoint of reionization by
$\Delta z_{\mathrm{end}}=0.12$. The midpoint is even more impacted as it is
now biased by $\Delta z_{\mathrm{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.

Figure 14: Posterior distributions on reionization global (upper panel) and
morphology (lower panel) parameters when assuming two measurements of the 21
cm power spectrum at $z=6.5$, 7.8 and $k=0.5\,\mathrm{Mpc}^{-1}$ for 100 hours
of observation with SKA-Low. The mock 21 cm data points used for the fit are
generated either with the simplified equation (7), in purple, or with the full
equation (2.1.1), in orange. Vertical and horizontal black lines correspond to
the ‘true’ values of the parameters, used to generate the mock data.

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 marginalised 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 the
results of Sec. 4.2. However, this idea would be less successful when applied
to real data, as this bias is model-dependent: From one simulation to another,
the amplitude and shape of the cross- and higher-order terms missing in
equation (7) will vary (Lidz et al., 2007; Georgiev et al., 2022).

##  Appendix D Forecast constraints on the Thomson Optical Depth

In section. 2.1.3, we derived the reionization history $x_{v}(z)$ using the
global history parameters $z_{\mathrm{end}}$ and $z_{\mathrm{re}}$ (see
equation 14) and can, therefore, derive the Thomson optical depth $\tau$. With
this in mind, we can investigate the constraining power of the forecast on
$\tau$ from the combined 21 cm and pkSZ power spectra analysis. Figure 15
presents the constraint on $\tau$ for each of the cases discussed throughout
this work. The dashed line represents a Gaussian centred on the true model
value $\tau=0.0649$ (Table 1) and with standard deviation the uncertainty
reported in Planck Collaboration et al. (2020). First, we note that, for each
mock data set, we are able to recover the true value of the optical depth
(values are biased by less than 1%) while achieving a tighter constraining
power than with current large-scale CMB data. Despite this, small variations
in constraining power are visible between cases.

Specifically, for our fiducial 1k2z case (see Sec. 4.2.1), there is a visible
positive bias, deviating from the true value by a fraction of a per cent. The
bias arises from the choice of a fixed $k$-scale for the 21 cm mock data
points. Because reionization is a time-dependent and inhomogeneous process,
the 21 cm spectrum from the EoR is expected to evolve both with redshift and
with $k$-scale. There is a characteristic $k$-scale below which the scale-
dependence of $P_{21}$ is governed by the matter power spectrum. This
characteristic scale evolves as reionization progresses (Furlanetto et al.,
2006). Hence, $k$-scales lower than this characteristic scale will be less
sensitive to reionization, and vice versa (see fig. 5 and 6 Georgiev et al.,
2022). While mock data in the latter half of reionization ($z=6.5,7.8$)
exclude larger values of $\tau$, data limited to $k=0.5$ Mpc-1 slightly favour
earlier reionization scenarios. The hiz case (in purple, $k=0.5$ Mpc-1 for
$z=7.8,10.4$) is biased towards lower values of the Thomson optical depth for
similar reasons. The forecast is conducted on data at higher redshift, in the
very early stages of reionization. The result is then a per cent bias of the
derived optical depth and increased uncertainty, preferring a lower value of
$\tau$ synonymous with a later reionization. It is only when the case where
both 21 cm mock data is at a fixed $z=6.5$ but at both high and low scales
$k=0.1,0.5$ Mpc-1 (2k1z in yellow) are considered, that we can recover an
unbiased estimate. This is because for the 2k1z mock data, the progress of
reionization is encoded differently in each of the $k$-scale of the 21 cm
power spectrum (see Sec. 5.1.2), resulting in a tighter constraint on $\tau$.
Naturally, the 3k2z case (seen in yellow), which combines the 1k2z and 2k1z
mock data, can further constrain $\tau$ without bias and with a lower
uncertainty.

Figure 15: Posterior distributions of the Thomson optical depth for the range
of models discussed within this work (see Tables 3 and 5). The dashed vertical
line represents the ‘true’ value $\tau=0.065$.