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

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 ℓ = 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


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  from the Cosmic Microwave Background (CMB).As CMB photons emitted during the recombination epoch travel through the intergalactic medium, they scatter off the free electrons produced during reionization.The Thomson scattering generates a polarization signal on scales larger than the horizon scale during reionization while it suppresses the temperature anisotropies on angular scales lower than the horizon size during reionization.Assuming an instantaneous reionization history, mea- * E-mail: ivelin.georgiev@astro.su.se surements of  = 0.051 ± 0.006 by the Planck Collaboration et al. (2020) indicate that the midpoint of EoR lies around a redshift of  re ∼ 8. On the other hand, observations of the fluctuations of the Lyman- optical depth caused by the Gunn-Peterson effect in high- quasar spectra (Bosman et al. 2022) and the inferred low mean free path of ionising photons ( MFP , Gaikwad et al. 2023) hint that reionization may extend past redshift six and complete by  end ∼ 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 (Mc-Quinn 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 Telescope1 (Ruhl et al. 2004) and the Planck satellite2 (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  pkSZ 3000 = 3.0 ± 1.0 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 upper limit on duration of reionization from 25 to 75 % of Δ < 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 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(Furlanetto et al. , 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  MFP through an empirical formula  trans ≈ 2/ 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)3 (Mertens et al. 2020), the Hydrogen Epoch of Reionization Array (HERA)4 (The HERA Collaboration et al. 2022;HERA Collaboration et al. 2023), the Murchison Widefield Array (MWA)5 (Trott et al. 2020;Yoshiura et al. 2021), the Giant Metrewave Radio Telescope (GMRT) 6 (Paciga et al. 2013), as well as the forthcoming Square Kilometre Array (SKA) 7 (Koopmans et al. 2015).However, no detection has been made and the derived upper limits set only weak constraints on astrophysical quantities (The HERA Collaboration 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- 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.In Bégin et al. (2022), the authors already demonstrate this complementarity at the level of the global 21 cm signal, showing that a combined analysis makes it possible to reconstruct a model-independent reionisation history, with no assumed parameterisation of the redshift-evolution.Here, we push this analysis a step further by including second-order statistics in the assessment and considering the 21 cm power spectrum.We place constraints on cosmic reionization by 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 constraining 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.Unless stated otherwise, we address comoving megaparsec as Mpc.

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.

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 with (Pritchard & Loeb 2012) Here,   =   / ρ − 1 is the baryonic density fluctuation and the bar denotes a spatial average and  H i is the fluctuation of the neutral fraction  H i field, i.e. the fraction of the IGM baryons that are neutral.Naturally,  H i = 1 −  H ii , where  H ii is the ionisation field.
The pre-factor  0 includes most of the spatially averaged information, including cosmological parameters, and the spin temperature information  s .We assume the limit where during cosmic reionization the gas in the IGM is sufficiently heated by the first ionising sources so  s ≫  CMB and the spin-temperature dependence in  0 is dropped.
Let us consider the electron fraction field where  H is the fractional quantity of electrons per Hydrogen atom.
If the first reionization of Helium is considered,  H ≃ 1.08.We distinguish the mass-weighted   and volume-weighted   average ionised fractions of H i atoms.By definition, ⟨  ⟩ ≡  H   , and, if we take the fluctuations of each element in equation ( 3), we have 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  H i = 1 −  H ii , we have where the spatial-and time-dependence have been omitted for simplicity.Taking the Fourier transform and averaging, we find where  21 is the power spectrum of the 21 cm signal fluctuations,   of the electron density fluctuations,   of the ionisation fluctuations,   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 ( 6) as Some further simplifications of equation ( 7) are assumed: We identify the mass-weighted   and the volume-weighted   ionised fractions and approximate the baryon power spectrum as a biased dark matter power spectrum   (, ) =  2  ()   (, ), such that the fully simplified expression used in this work, unless stated otherwise, is 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   .We will now look for a similar relation for the kSZ angular power spectrum.

Derivation of the patchy kSZ angular power spectrum
The angular power spectrum of the pkSZ effect at multipole ℓ can be derived from the electron density power spectrum   under some assumptions 8 ., such that (Mesinger et al. 2012;Gorce et al. 2020): with   being the Thomson scattering cross-section, n () and () are the mean electron density and the comoving distance for redshift , respectively.Moreover,  , is the power spectrum of the curl component of the momentum field q , such that (2 where   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  and  at a certain wave number , as where where  = k • k′ and the -dependencies have been omitted for simplicity.  (, ) is the power spectrum of the free electrons density fluctuations and   is the free electrons density -velocity cross-spectrum.The latter is computed as where  is the scale factor,  the linear growth rate and the bias is defined by the ratio   (, ) 2 ≡   (, )/   (, ).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.

Electron power spectrum
We assume the   (, ) parameterisation introduced in Gorce et al. (2020), that is 8 These assumptions include the Limber approximation, the assumption that the velocity power spectrum is a biased linear matter power spectrum, and the omission of third and fourth order correlation terms (Gorce et al. 2020).The latter implies variations of the order of 10% (∼ 0.05 2 ) in the patchy kSZ amplitude (see Alvarez 2016, for details) where  is the electron drop-off frequency, which can be related to the typical size of ionised bubbles during the EoR, log 10 ( 0 ) is the large-scale   amplitude, related to the variance of the electron field during reionization.The redshift-independent baryon-dark matter bias   () is given by the adapted Shaw et al. (2012) parameterisation, that is where   = 9.4 Mpc −1 and  = 0.5 are constant with redshift and calibrated on the EMMA simulation (Aubert et al. 2015), which includes coupled radiative transfer and hydrodynamics and is therefore sensitive to the thermal and reionization history.The global reionization history   () in equation ( 12) is defined according to the following parameterisation (Douspis et al. 2015;Gorce et al. 2022): where  early = 20 is the redshift for which   ( early ) ≈ 10 −4 , the ionisation leftover from recombination.Random forests are used to speed up computations 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 CAMB 9 (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.

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 , we fit for the reionization mid-and endpoint,  re and  end , as well as for the reionization morphology parameters log 10 ( 0 ) and , 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   .With the latter, we obtain the dimensionless 21 cm power spectrum Δ 2 21 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: where uncertainty in the measurements Δ 2 noise , Δ 2 var , and  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  parameter used in this method represents the variance of the chain means compared to the mean of chain variances.In the cases where  ≲ 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 %.

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 where  21 () is the redshifted 21 cm wavelength,  is the bandwidth centred on redshift  and Δ  is its length in comoving Mpc.We write the system noise  sys = 100 + 300( obs / 0 ) −2.55 K, for  0 = 150 MHz and  obs () = / 21 () is the observed 21 cm frequency.The total integration time is  int ,   () is the comoving distance to redshift , and  base the number of baselines, roughly equal to the number of antennas squared  2  .The total collecting area of the telescope is  coll , such that  coll ≡    2  where   is the radius of the antenna.We have  core the core area of the array and  eff the effective collecting area of each antenna.We use the values presented in Table 2 to estimate the noise for four different experiments: MWA, HERA, LOFAR, and SKA-Low.Note that for the latter and LOFAR, we consider each station as a single antenna.
We also include the contribution of sample variance to the uncertainty of the measurement following the relation (see eqs. 9 &10 of Mellema et al. 2013): where  is the logarithmic binning of Δ =   and  is the survey  size, which can be expressed as and  is the field of view.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  ⪅ 0.3 Mpc −1 because of its strong -dependency (Δ 2 var ∝  −3/2 ).The opposite is true of the noise estimate for MWA, which has a larger field of view and as Δ 2 var ∝ 1/ 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  pkSZ 3000 = 3.0 ± 1.0 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).

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 -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  esc = 0.2.In Fig. 2, we show simulation slices of the differential surface brightness temperature of the 21 cm signal.The cosmology used is consistent with Planck Collaboration et al. ( 2020) and summarised in Table 1.The reionisation mid-and endpoint values given in Table 1 are obtained by fitting the asymmetric parameterisation from equation ( 14) to the reionisation history of the simulation.Similarly, the values of the morphology parameters log 10 ( 0 ) and  from Table 1 are derived by fitting the RSAGE electron density power spectrum with equation ( 12) (see Gorce et al. 2020, for details of the fit).We have checked that the baryon power spectrum from RSAGE is well described using equation ( 13) for the given values of   and .Overall, we find the model of the electron density power spectrum in equation ( 12) to be a good fit for the redshift and -scales explored in this work.In the following section, we therefore limit our inspections to the accuracy of the 21 cm power spectrum reconstruction.

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

Reconstructing the power spectra
We want to understand the role of each contributing term in equation (6), especially the ones we have discarded in our simplified  21 reconstruction (equation 8), in order to assess its accuracy.We reconstruct and examine the 21 cm power spectrum from the RSAGE simulation described in section 3. The top panel of Fig. 3 showcases the true  21 plotted against different reconstruction precision levels: -The recons no cross case corresponds to the highest simplification level, corresponding to equation ( 8), and the approximation we will use in the remaining of this work.Here, it is computed with the true   (, ) spectrum obtained from the sim ulation, and not equation ( 12), therefore, we only assess the accuracy of the 21 cm reconstruction.The contributions from both third-order terms and cross-terms are discarded.We also include an additional case, recons no cross (using   ), where we lift the assumption of equating the mass-weighted and volume-weighted ionization fractions.-For recons with cross, we add the cross power spectra contributions: −2  ()  (, ) 10 , -For recons with third order, we further add the three-point power spectrum contribution: −2  () , (, ) 10 .This corresponds to the full form of the 21 cm power spectrum, without simplification.
We look at the relative contribution of each of the terms in equation (6) to the 21 cm power spectrum in the RSAGE simulation.This is represented in the middle panel of Fig. 3 for  = 0.141 Mpc −1 .We see that for  > 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 (  ≲ 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  21 look identical in the upper panel.Beyond   ≈ 10%, when the first luminous sources have reionised their local over-densities and 10 Note that this contribution can be positive or negative, as illustrated in Fig. 3.The coloured lines correspond to different reconstruction levels (see text for details).recons no cross corresponds to equation ( 8), which is used for the forecast in section 4.2.Additionally, recons no cross (with   ) corresponds to equation ( 8) but for   ≠   .In recons with cross, the cross power spectrum is included: −2  ()   (, ).Lastly, recons with third order adds three-point power spectrum contribution to recons with cross: −2  ()  , (, ), which corresponds to the full form of the 21 cm power spectrum in equation ( 6).Middle Panel: Ratio of each component in equation ( 6) against the total 21 cm power spectrum.In this panel, solid (dashed) lines represent a positive (negative) contribution.Bottom Panel: Reionization history of the RSAGE simulation.
the ionisation field is correlated with the density field,  21 is primarily determined by the difference between the electron density power spectrum and the   cross-spectrum.The third-order term   , 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 − 2  ()] factor present in equation ( 7).
We can form a clearer understanding of the impact of each term by examining the top panel of Fig. 3, where the power spectrum of the 21 cm signal is calculated according to equation (8), in blue, while the full model presented in equation ( 6), in green.We see that the simplified expression in the recons no cross model overestimates the true power throughout reionization.Early on, recons no cross follows the true recons with third order model.With the onset of reionization, a positive bias of an order of two appears between the true spectrum and its model, reaching a global maximum at the midpoint of reionization.Additionally, in the simplified model, reionization is delayed by four per cent compared to the true model.Hence, our simplified model will tend to overestimate  21 and will be biased towards late reionization scenarios, the implica-tions of which are discussed in section 4.2.Additionally, we note that by removing the   ≡   assumption and replacing  2    by  2    in equation ( 7) results in a higher amplitude (seen in purple) compared to the recons no cross model (light blue) using the volume-weighted fraction.The increased bias of the model is due to the mass-weighted ionisation fraction being larger than the volume-weighted ionisation fraction   >   at all redshifts in inside-out models of reionisation (see, e.g., Dixon et al. 2016).We have compared this across the -scales accessible within the RSAGE simulation volume, and we find the   ≡   simplification leads to a  21 model closer to the truth when excluding the higher-order contributions.
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 -scales and redshifts.Generally, the negative cross-correlation between the density and ionisation fields weakens at higher -scales as well as with the progression of reionization.We see that by excluding the higherorder power spectrum and varying the -scale, the recons with cross model converges to the recons with no cross model on large -scales.In Appendix A, we look at the possibility of using equation ( 7) to reconstruct the pkSZ angular power spectrum from an interferometric measurement of the 21 cm power spectrum.
In summary, the different terms comprising our expression of the 21 cm power spectrum, described in equation ( 6), are non-trivial.The importance of each term presented in Fig. 3 is shown to evolve with redshift and -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 scales, resulting in a positive bias on the end of reionization, which we can easily quantify.In contrast, including the cross-power spectrum   contribution to the model introduces an evolving amplitude bias with -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 mock data points and each sampled model.
For the sake of clarity, we list below the different layers of assumptions for our model in equation ( 8) and how we address them.
(i) The baryon power spectrum follows equation 13 when using the CAMB-derived theory non-linear matter power spectrum.Comparing to the spectrum directly measured from RSAGE, we find a good match with the model across all redshifts and -scales covered in this work.(ii) The global reionisation history follows equation 14.We have compared the evolution obtained with equation 14 and the parameters of Table 1 to the RSAGE history seen at the bottom panel Fig. 3 and found to be a good match for the redshifts considered in this paper.(iii) The mass-and volume-weighted ionisation fractions   and   are identical.As seen in the upper panel of Fig. 3, the positive bias induced by excluding the higher order terms is slightly alleviated for the recons no cross model for the inside-out reionization models considered in this work, for which   <   .(iv) The true electron power spectrum in RSAGE follows equation ( 12).
This model has been compared against the RSAGE simulation within this work while the best-fit morphology parameters are reported The 21 cm power spectrum RSAGE (seen in black) compared against our model (seen in fuchsia) for redshifts  = 6.5, 7.8 (  = 0.87, 0.37) as thick solid and dashed lines, respectively.Our model 21 cm power spectrum has been generated using equation ( 8) for the reference parameters in Table 1 and the CAMB non-linear matter power spectrum.
in Table 1 (see appendix B.2 of Gorce et al. 2020, for a detailed discussion).
We conclude this assessment of our model uncertainties by comparing, in Fig. 4, the true 21 power spectrum from RSAGE against the one generated by our fully simplified parameterisation, given in equation ( 8).Overall, the parameterisation performs remarkably well on the -scales considered in this work ( ∼ 0.1 Mpc −1 ).At high -scales ( > 1.0 Mpc −1 ), the model in equation ( 12) scales as  −3 and quickly approaches the baryon power spectrum.The result is a negative 21 cm power spectrum as can be observed in Fig. 4 at z = 6.5.Therefore, the model should not be used as is at very small scales.At low-, differences between the model and the RSAGE power spectrum will be sensitive to sample variance due to the simulation box size and hence have not been studied within this work.
We discuss the effect of the aforementioned bias on our results in section C, where we derive the mock data using the RSAGE simulation and no approximation.In Sec. 6, we discuss future improvements and possible avenues in modelling the higher-order terms in our model.

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 ℓ = 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.Posterior distributions on reionization mid-and endpoint (upper panel) and the morphology parameters log 10 (  0 ) and  (bottom panel), considering a detection of the pKSZ power spectrum at ℓ = 3000 and of the 21 cm power spectrum at redshifts  = 6.5, 7.8 and  = 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.

Forecast with SKA-Low
Our primary scenario is based on detection of the 21 cm power spectrum at redshifts  = 6.5 and 7.8, both at  = 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 ℓ = 3000 with a ten per cent error bar such that   =3000 = 0.86 ± 0.09  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 -scale.
In Fig. 5, we show the joint probability distributions for our four sampled parameters (see appendix B for a discussion on parameter degeneracies).The best-fit parameters, their 1 error, and Gelman-Rubin convergence parameters are presented in Table 3 11 .
11 The Gelman-Rubin parameters for morphology parameters of the 1000 h case, given in Table 3, are of order  − 1 ≈ 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 Our results indicate that we can recover the true values of  re and  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,  = 0.066 ± 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 -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. 5, we can place loose lower limits on log 10 ( 0 ), although the 1000 hrs case is more constraining on both distributions: Adding integration time additionally provides a lower limit on .Note that the ranges of the axes in the figure are smaller than those of the priors from Table 1 (while the limit on  is extended), confirming that these limits are not only informed by the priors.We further investigate whether the choice of -scale has an effect on these results in section 5.1.2.
To get more insight into the results above, we look in more detail at the dependence of our two observables on the model parameters.Fig. 6 presents the reionization history, the pkSZ power spectrum, and the 21 cm power spectrum at  = 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  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  end decreases the duration of the EoR, suppressing the pkSZ amplitude.This explains the strong anti-correlation observed between the two parameters in Fig. 5.Note that for the model where  end = 7, naturally, the 21 cm power spectrum is zero at  = 6.5.
On the other hand, log 10 ( 0 ) and  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 ( 0 ) or , as seen in the third and fourth columns of Fig. 6, results in a similar impact on both observables.The log 10 ( 0 ) parameter is effectively a measurement of the large-scale amplitude of   at high redshift (before the reionization midpoint), whilst the   power breaks as  −3 for scales  > , such that  can be interpreted as the minimal size of ionised regions during reionization.Hence, varying log 10 ( 0 ) directly impacts the amplitude of both observables, whilst varying  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  increases.This dependence on the observables on these two parameters explains the degeneracy between the morphology parameters seen in Fig. 5.For the 1k2z case, where both data points are the posteriors of relevant parameters do not change significantly as iterations increase, confirming that the results still hold good qualitative significance.at a fixed -scale, a lower value log 10 ( 0 ) can be compensated by a higher  value, within the measurement error of the data.However, this degeneracy can be reduced in the presence of 21 cm data points at different -scales or by kSZ data points at different multipoles, since  impacts the shape of both power spectra.We explore the -scale dependency of our chosen data in section 5.1.2.

Constraining power of each dataset
In this subsection, we examine each dataset's ability to individually impose constraints on reionization.In Fig. 7, we present the joint posterior distribution of the sampled parameters for three different cases assuming 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. 5.To better understand the role of each parameter in each of the aforementioned cases, we firstly group our analysis by examining the global reionization parameters  re and  end , then the astrophysical log 10 ( 0 ) and , in Fig 7.
Examining the two-dimensional  re and  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 d ≡  re −  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. 6).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  = 6.5 naturally excludes models where reionization finishes earlier than that redshift, hence, the forecast gives a lower limit on  end .Lastly, models at the lower left of the two-dimensional posterior ( re ≲ 7 and  end ≲ 6) are excluded by the 21 cm data at  = 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 ( 0 ) and , 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. 6: Increasing either of the morphology parameters results in a boosting of the pkSZ and 21 cm power spectra, resulting in a degeneracy between them.The pkSZ data seems to favour a low , which in order to match the measurement is then compensated by a low log 10 ( 0 ).Consequently, while the true pkSZ power spectra peaks at ℓ ≈ 3000, here there are multiple models for which the pkSZ reaches its maximum amplitude at ℓ ≈ 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  re and  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  end , it cannot constrain  re .Conversely, we find that the pkSZ measurement mainly constrains the duration of reionization (d =  re −  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 -scales and we discuss how we can further constrain them in the following section.

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 -ranges at which the data points were selected in section 4.2.1, as well as the noise models and study the impact on our results.We also explore how the choice of input parameters affects the forecast.Note that for all cases in this discussion, we limit ourselves to the fixed number of data points as in section 4.2.1 to quantify the role of each data point and its driving physical processes.The interested reader can refer to section 5.3 for results on how an additional 21 cm data point can significantly improve the forecast.

Choice of redshift
We explore the constraining ability of the forecast based on the choice of redshifts for the 21 cm observation.To do so, we examine the SKA-Low 1000h case (hiz) and vary the redshift of the observed 21 cm power spectrum from  = 6.5 and 7.8 to  = 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  re −  end and log  0 −  are presented in Fig. 8 in purple.Compared to the fiducial case (in blue) two main differences are apparent.First, the data are now compatible with models whose mid-and endpoint are generally at higher redshifts (early and rapid reionization scenarios).
For comparison, for the 1k2z cases in section 4.2.1, the fact that the 21 cm signal was non-zero at  = 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  = 10.4 has limited constraining power on the global history parameters but provides better constraints on morphology parameters, as illustrated in the bottom panel of Fig. 8, in orange.The reason for this is that at that redshift,   ≪ 1, such that  21 ∼   +   2   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).

Choice of 𝑘-scales
Next, we investigate the constraining power of the forecast using two detections of the 21 cm power spectrum at  = 6.5 and  = 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- 'tail' which encodes information about the IGM on small scales.Fig. 8 illustrates the constraining power of this approach and results for 2k1z are shown in yellow, in contrast to the fiducial 1k2z case, shown in blue.While the data have clear constraining power on the  re and  end values, the two-dimensional posterior distribution exhibits a tail and is biased towards late reionization models ( end = 5.29 ± 0.55, whilst the true value is at redshift 6.15).Observing a non-zero 21 cm data point close to the true  end naturally excludes early reionization models.Despite the increased measurement error, the low- (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- 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 -scales results in a tighter constraining power on log 10 ( 0 ) and , although the two-dimensional posterior is biased towards lower values of the parameters, corresponding to later reionization models.Examining and comparing the 21 cm power spectra in Fig. 6 reveals that the differences in -scale are most prominent for  > 1.0 Mpc −1 scales.Since the high -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 -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  re - end posterior as well as biased morphology parameters.One way to improve constraints on  end would be to consider 21 cm data for  > 1.0 Mpc −1 scales, or to limit the prior on  end to values allowed by measurements of Ly  absorption in quasar spectra (Bosman et al. 2022).

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 considered12 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 -scale of  = 0.1 Mpc −1 for redshifts  = 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. 9 and Table. 4.
Results for the MWA case are overall consistent with section 4.2.1, with key differences emerging due to the higher uncertainty on the 21 cm power spectrum, which is approximately 1000 times larger than that of the SKA (see Fig. 1).As shown in orange in the upper panel of Fig. 9, we are still able to constrain  re = 7.02 ± 0.49.However, the spread in the distribution has increased compared to the SKA-Low case in blue (for which  re = 7.39 ± 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  end < 7.5.A lower constraining power is also observed for the morphology parameters.We can primarily obtain upper limits for log 10 ( 0 ) < 3.75 and  > 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  = 6.5 has a higher amplitude and a lower value of , and the pkSZ power spectrum peaks at ℓ ≈ 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 ℓ = 2000 or more data on the 21 cm power spectrum at lower- values.
Regarding the LOFAR case, results are roughly similar to what was observed in Sec.5.1.2where the 21 cm mock data is taken at a fixed -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- 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. 9, we can constrain  re < 8.0 and can only exclude only late reionization models for which  end > 5.5.Compared to the MWA case, having higher- 21 cm data leads to higher constraining power on the morphology parameters.While we can only place an upper limit on  > 0.1, the constraint on log 10 ( 0 ) = 3.15 ± 0.33 is fairly comparable to that of 1k2z 1000 hours (where log 10 ( 0 ) = 3.04 ± 0.32).While this is mostly related to the choice of data, for example, the hiz case can constrain 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.

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 -scale ( = 0.5 Mpc −1 for  = 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  re and  end .On the other hand, with two data points at fixed redshift but different scales ( = 6.5 for  = 0.1, 0.5 Mpc −1 , 2k1z case), we can constrain the morphology parameters log 10 ( 0 ) and .A natural way to extend these results is to venture into the multi-scale regime (3k2z case, seen in brown in Fig. 10) by combining both of the aforementioned data sets such that the 21 cm power spectrum is measured at  = 6, 5, 6.5, 7.8 and  = 0.1, 0.5, 0.5 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. 10), we improve the accuracy of our constraints of the reionization parameters: We get  re = 7.37 ± 0.07 and  end = 6.15 ± 0.03, significantly reducing the 1 standard deviation (especially for  end ).Consequently, we can tightly constrain the parameter space around the true value and retrieve the correct Thomson optical depth with a precision of Δ = ±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).10.Posterior distributions on reionization mid-and endpoint (upper panel) and the morphology parameters log 10 (  0 ) and  (bottom panel), considering a detection of the pKSZ power spectrum measurement at ℓ = 3000 and of the 21 cm power spectrum at redshifts of  = 6.5, 6.5, 7.8 at  = 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, 5).
Constraints on the morphology parameters are significantly tighter: We measure log 10 ( 0 ) = 3.12 ± 0.01 and  = 0.164 ± 0.001.This improvement is partly explained by the fact that adding data at small scales breaks the degeneracy between the morphology parameters.We can intuitively explain this by referring to the results in section 5.1.2:The 3k2z case can be understood as the combination of the 1k2z and 2k1z cases and we see on Fig. 8 that accessing the redshift-evolution of the 21 cm signal with data points at different redshifts removes the bias on the morphology parameters obtained with 2k1z.
Looking at these results, it is tempting to think that all the constraining power comes from the three 21 cm data points.However, we have also examined an additional case where we only consider the 21 cm mock data from 3k2z and exclude the pkSZ data point.We find that, while we can recover the morphology parameters with similar, albeit weaker, constraints, the global history parameters are significantly more degenerate.Indeed, the information the pkSZ data point provides on the duration of reionization remains crucial for the accuracy of the forecast.Measuring the pkSZ spectrum at different multipoles would enable tighter constraints on the shape parameter , whether at lower (ℓ = 2000) or larger (ℓ = 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 ℓ ≲ 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.

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   (, ) presented in Gorce et al. (2020).The resulting equation ( 6), while linking both observables, contains non-trivial cross-correlation and second-order terms of the ionised and density fields.We choose to ignore such terms in this work, leading to equation ( 8).Looking at the semi-numerical RSAGE simulation (Seiler et al. 2019), we find that this assumption leads to overestimating the 21 cm power on all scales and redshifts  ≲ 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,  re and  end , as well as for two reionization morphology parameters log 10 ( 0 ) and .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 ℓ = 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  = 0.5 Mpc −1 and  = 6.5, 7.8 (referred to as the 1k2z case, shown in blue in Fig. 5 and throughout this work), we are able to constrain the global reionization parameters  re = 7.42 ± 0.13 and  end = 6.15 ± 0.04 and provide limits on the morphology parameters log 10 ( 0 ) > 2.64 and  < 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 =  re −  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 -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- measurement excludes earlier reionization models, for which the 21 cm signal would be zero.Second, at high redshift, when the global ionisation fraction is   ≪ 1, the data will primarily determined by log 10 ( 0 ) and .In our tests, for 21 cm mock measurements at high(er) redshift (hiz case,  = 7.8 and 10.4 at  = 0.5 Mpc −1 ), the constraints on the global reionization parameters  re and  end are loosened compared to the 1k2z case, at  = 6.5 and 7.8 (Fig. 8).On the other hand, it is precisely the constraining of the high- 21 cm power spectrum which allows for firmer constraints on the morphological parameters: log 10 ( 0 ) = 3.10 ± 0.17 and  = 0.17 ± 0.02.
We also observe that the -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- 'tail', which encodes information about reionization morphology on small scales.By choosing two mock 21 cm power spectra at a fixed redshift (2k1z,  = 6.5 at  = 0.1, 0.5 Mpc −1 ), we find that the constraints on the global reionization parameters are weakened to upper limits (Fig. 8, in yellow) but constraints on the morphology parameters are improved.However, the recovered values are biased by 6%.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  > 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- modes, which are the most sensitive.However, such measurements are non-trivial because the measurement noise scales with -scale (see Fig 1), making them less appealing in the context of early SKA-Low observations.
It is also worth noting that we limit our main results to simplistic cases with two 21 cm power spectrum observations in order to provide a proof-of-concept approach and to better showcase the workings of the forecast.Further tests where we increase the number of data points from two to three (e.g., the 3k2z case seen in section 5.3) show the full potential of this approach: Highlighted in brown in Fig. 10, we see that as more and better measurements of the pkSZ and the 21 cm power spectrum become available, this method will not only enable constraints on the global history of reionization but also its morphology.We obtain accurate constrains of the reionization parameters  re = 7.37 ± 0.07 &  end = 6.15 ± 0.03 and the morphology parameters log 10 ( 0 ) = 3.12 ± 0.01 &  = 0.164 ± 0.001, significantly reducing the 1 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. 9 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  re ≤ 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  end ≥ 6.5.Meanwhile, data from LOFAR ( chosen at  = 0.1 Mpc −1 for  = 8.3, 9.1) is more sensitive to the morphology parameters, notably constraining log 10 ( 0 ) = 3.15 ± 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 higherorder and cross-correlation terms in our simplified expression of the 21 cm power spectrum in equation ( 8).Potential future developments would be develop analytic models of the cross-power spectrum and include them in the derivation (McQuinn et al. 2005;Schneider et al. 2023).However, this approach would most likely not capture the full complexity present at all -scales and could be model-dependent.Another option would be to model the higher-order and the cross terms as constant biases and account for them as nuisances parameters in the forecast (see the discussion in appendix C).Lastly, we plan to conduct our analysis for 21 cm data over broader ranges of redshift and -scales to explore the potential of joint analysis to identify and remove systematics such as residual foregrounds (see Bégin et al. (2022) for such an analysis using the global 21 cm signal).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  ∥ modes that can be observed by these instruments. = 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 ( − 1 ≈ 0.001).The corner plot in Fig. B1 clearly illustrates this, showcasing the relation between the global history and morphology parameters as well as the derived Thomson optical depth .The parameter space is well sampled, with clear correlations seen only between either  re −  end and log 10 ( 0 ) − .A correlation is noticeable between the Thomson optical depth and each of the global history parameters.This is more than expected, as  is an integral of the ionisation history, defined as a function of  re and  end in equation ( 14).We discuss the constraints on  for our different test cases in appendix.D in more detail.

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 crosscorrelations, 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  21 (, ), including the cross-and higher-order terms, given in equation ( 6).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. C1.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 Δ end = 0.12.The midpoint is even more impacted as it is now biased by Δ re = 0.36, nine times more than with the simplified mock data.Note that these modified reionization histories result in a large optical depth:  = 0.070 ± 0.002.Hence, a way to partially mitigate these biases could be to impose a Gaussian prior on .Similarly, the posterior distributions of the morphology parameters are shifted compared to the results with the simplified mock data points.
Although these results exhibit a strong bias in all parameters, this bias is nevertheless well understood and easy to model.A potential way of improvement would be to add a nuisance parameter to the fit, as a pre-factor to the model 21 cm power spectrum.This extra parameter would be 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   () using the global history parameters  end and  re (see equation 14) and can, therefore, derive the Thomson optical depth .With this in mind, we can investigate the constraining power of the forecast on  from the combined 21 cm and pkSZ power spectra analysis.Figure D1 presents the constraint on  for each of the cases discussed throughout this work.The dashed line represents a Gaussian centred on the true model value  = 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 -scale for the 21 cm mock data points.Because reionization is a time-dependent .One-and two-dimensional posterior probability distributions for our model parameters and one derived parameter (the optical depth ) when fitting a detection of the pKSZ power spectrum measurement at ℓ = 3000 with a ten per cent uncertainty and of the 21 cm power spectrum at redshifts of  = 6.5, 6.5, 7.8 at  = 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.
and inhomogeneous process, the 21 cm spectrum from the EoR is expected to evolve both with redshift and with -scale.There is a characteristic -scale below which the scale-dependence of  21 is governed by the matter power spectrum.This characteristic scale evolves as reionization progresses (Furlanetto et al. 2006).Hence, -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 ( = 6.5, 7.8) exclude larger values of , data limited to  = 0.5 Mpc −1 slightly favour earlier reionization scenarios.The hiz case (in purple,  = 0.5 Mpc −1 for  = 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  synonymous with a later reionization.It is only when the case where both 21 cm mock data is at a fixed  = 6.5 but at both high and low scales  = 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 -scale of the 21 cm power spectrum (see Sec. 5.1.2),resulting in a tighter constraint on .
Naturally, the 3k2z case (seen in yellow), which combines the 1k2z and 2k1z mock data, can further constrain  without bias and with a lower uncertainty.3 and 5).The dashed vertical line represents the 'true' value  = 0.065.

Figure 1 .
Figure1.Measurement errors on the 21 cm power spectrum from RSAGE (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 -scales are presented for redshifts  = 6.5, 7.8 as thick solid and dashed lines, respectively.

Figure 2 .
Figure 2. Slices of the differential surface brightness temperature  21 for the RSAGE simulation.The upper panel is at a redshift of  = 7.83 and has a volume-averaged ionisation fraction of   = 0.37, while for the bottom panel, the same values are  = 6.52 and   = 0.87.Both redshifts correspond to the mock data points used in Section 4.2.

Figure 3 .
Figure 3. Upper Panel: The dimensionless 21 cm power spectrum at  = 0.141 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 (8), which is used for the forecast in section 4.2.Additionally, recons no cross (with   ) corresponds to equation (8) but for   ≠   .In recons with cross, the cross power spectrum is included: −2  ()   (, ).Lastly, recons with third order adds three-point power spectrum contribution to recons with cross: −2  ()  , (, ), which corresponds to the full form of the 21 cm power spectrum in equation (6).Middle Panel: Ratio of each component in equation (6) against the total 21 cm power spectrum.In this panel, solid (dashed) lines represent a positive (negative) contribution.Bottom Panel: Reionization history of the RSAGE simulation.

Figure 6 .
Figure 6.Evolution of the reionization history (first column), the pkSZ angular power spectrum (second column), and the dimensionless 21 cm power spectrum at redshifts  = 6.5, 7.8, considered in section 4.2 (last two columns, respectively) as a function of our four model parameters.The values of  are shown as vertical dotted lines in the final row with the corresponding colour.Note the difference in the shape of the 21 cm power spectrum at high- compared to Fig. 1 is discussed in Sec.4.1.

Figure 7 .
Figure7.Posterior distributions on reionization mid-and endpoint (upper panel) and the morphology parameters log 10 (  0 ) and  (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.

Figure 8 .
Figure8.Posterior distributions on reionization mid-and endpoint (upper panel) and the morphology parameters log 10 (  0 ) and  (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 ℓ = 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.

Figure 9 .
Figure 9. Posterior distributions on reionization mid-and endpoint (upper panel) and the morphology parameters log 10 (  0 ) and  (bottom panel) when fitting 21 cm and pkSZ power spectra jointly, for 1000 hours of MWA, LO-FAR, & SKA observations, seen in grey, orange, and blue.The vertical and horizontal lines show the 'true' values used to generate the mock data.

Figure A3 .
Figure A3.Longitudinal (lower panel) and transversal (upper panel) -modes observable by an interferometer, depending on its characteristics such as baseline length  or frequency resolution.The bandwidth is fixed to 15 MHz.In the lower panel, the frequency resolutions of SKA-Low, HERA and LOFAR are shown as reference and result in an upper limit on the  ∥ modes that can be observed by these instruments.

Figure C1 .Figure D1 .
Figure C1.Posterior distributions on reionization global (upper panel) and morphology (lower panel) parameters when assuming two measurements of the 21 cm power spectrum at  = 6.5, 7.8 and  = 0.5 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 (6), in orange.Vertical and horizontal black lines correspond to the 'true' values of the parameters, used to generate the mock data.

Table 1 .
Reference cosmological and reionization parameters used to generate the mock data, based on the tuned to best fit the RSAGE simulation.All parameters apart from  re ,  end , log 10 (  0 ), and  are fixed.

Table 2 .
Parameters describing the mock observations by different telescopes used to derive our 21 cm power spectrum errors (equation 17).

Table 3 .
Best fit values of the model parameters for the 1k2z cases, their corresponding 1 uncertainty and Gelman-Rubin parameter.Best-fit values and errors are also given for derived parameters such as .

Table 4 .
Best fit values of the distributions in Fig.9and their corresponding 1 uncertainty as well as the Gelman-Rubin convergence diagnostic for each parameter.

Table 5 .
Best fit values of the distributions in and their corresponding 1 uncertainty as well as the Gelman-Rubin convergence diagnostic for each parameter.