Improved constraints on reionisation from CMB observations: A parameterisation of the kSZ effect

We show that, in the context of patchy reionisation, an accurate description of the angular power spectrum of the kinetic Sunyaev-Zel'dovich (kSZ) effect is not possible with simple scaling relations between the amplitude of the spectrum and global parameters such as the reionisation midpoint and its duration. We introduce a new parameterisation of this spectrum, based on a novel description of the power spectrum of the free electrons density contrast Pee(k,z) in terms of reionisation global history and morphology. We directly relate features of the spectrum to the typical ionised bubble size at different stages in the process, and subsequently to the angular scale at which the patchy kSZ power spectrum reaches its maximum. We successfully calibrate our results on a custom set of advanced radiative hydrodynamical simulations and later find our parameterisation to be a valid description of a wide range of other simulations and therefore reionisation physics. In the end, and as long as the global reionisation history is known, two parameters are sufficient to derive the angular power spectrum. Because of their straightforward physical interpretation, related to the morphology of reionisation, these parameters will be constrained by future 21cm intensity mapping experiments. Conversely, a precise measurement of the patchy kSZ power spectrum will constrain these parameters, and consequently the physics of reionisation. This parameterisation is therefore found to be an innovative tool to extract information about reionisation from CMB data.


Introduction
From the launch of the Cosmic Background Explorer (COBE) in 1989 to the publication of the latest results of the Planck satellite in 2018 (Planck Collaboration 2018), the study of the Cosmic Microwave Background (CMB) has triggered a tremendous amount of research. Cosmological parameters have been estimated with exquisite precision and our knowledge of cosmic inflation has been greatly improved. Along the line of sight, the primordial part of the CMB signal is largely modified by the interaction of CMB photons with structures formed later in the Universe. Notably, their interaction with free electrons in the Intergalactic Medium (IGM) modify the shape and amplitude of the measured CMB temperature and polarisation power spectra. The presence of these electrons is the result in particular of cosmic reionisation, an era potentially extending from redshift z ∼ 15 to z ∼ 5 when the first galaxies are thought to have ionised the neutral hydrogen and helium in the surrounding IGM.
CMB photons lose energy from scattering off low-energy electrons. In CMB data analysis, this effect is accounted for when computing the Thomson optical depth. To do so, one needs to assume a global history of reionisation, i.e. a redshiftevolution for the IGM global ionised fraction x e (z). In standard Boltzmann solvers used to compute theoretical predictions in CMB data analysis such as the CAMB code (Lewis et al. 2000, Howlett et al. 2012) 1 , the reionisation scenario used is a steplike transition of x e (z), where the global ionised fraction jumps from 10% to 75% over a (fixed) redshift interval of ∆z = 1.73 (Planck Collaboration 2016a). However, this parameterisation does not match simulations and observations well as we expect the ionisation fraction to slowly rise when the first sources light up, before taking off once about 20% of the IGM is ionised (Robertson et al. 2015, Greig & Mesinger 2016, Gorce et al. 2018). This minimal model can have a huge impact on reionisation constraints: the value of τ inferred from Planck 2016 data varies from 0.066 ± 0.016 for a step-like process to 0.058 ± 0.012 for a more accurate description (Douspis et al. 2015;Planck Collaboration 2016b). It is therefore essential to take the asymmetric evolution of x e (z) into account when trying to accurately constrain reionisation, and global parameters such as the reionisation midpoint z re and duration ∆z are not sufficient.
CMB photons can also gain energy from scattering off electrons with a non-zero bulk velocity relative to the CMB restframe in a process called kinetic Sunyaev-Zel'dovich effect (hereafter kSZ effect, see Zeldovich & Sunyaev 1969;Sunyaev & Zeldovich 1980). This interaction adds power to the CMB temperature spectrum on small angular scales ( 2000, i.e. smaller than about 5 arcminutes), where secondary anisotropies, including kSZ, dominate the signal. The impact of kSZ on the CMB power spectrum is often split between the homogeneous kSZ signal, coming from the Doppler shifting of photons on free electrons homogeneously distributed throughout the IGM once reionisation is over, and the patchy kSZ signal, when CMB photons scatter off isolated ionised bubbles along the otherwise neutral line of sight. Therefore the kSZ power spectrum is sensitive to the duration and morphology of reionisation , Mesinger et al. 2012: for example, the patchy signal is expected to peak around ∼ 2000, corresponding to the typical bubble size during reionisation Iliev et al. 2007).
Secondary anisotropies only dominate the primordial power spectrum on small scales where existing all-sky surveys such as Planck perform poorly. The observational efforts of the groundbased Atacama Cosmology Telescope (ACT) 2 and the South Pole Telescope (SPT) 3 have allowed upper constraints to be put on the amplitude of the kSZ power spectrum at = 3000. Using ACT observations at 148 GHz, Dunkley et al. (2011) find D SZ 3000 ≡ ( + 1) C SZ =3000 /2π = 6.8 ± 2.9 µK 2 at 68% confidence level (C.L.) for the sum of thermal and kinetic SZ. In a first analysis, Reichardt et al. (2012) derive from the three frequency bands used by SPT D kSZ 3000 < 2.8 µK 2 (95% C.L.). This limit is however significantly loosened when anti-correlations between the thermal SZ effect (tSZ) and the Cosmic Infrared Background (CIB) are considered. By combining SPT results with large-scale CMB polarisation measurements, Zahn et al. (2012) are subsequently able to constrain the amplitude of the patchy kSZ by setting an upper limit D patchy 3000 ≤ 2.1 µK 2 (95% C.L.) translated into an upper limit on the duration of reionisation ∆z ≡ z (x e = 0.20) − z (x e = 0.99) ≤ 4.4 (95% C.L.), again largely loosened when CIB×tSZ correlations are considered. Using Planck's large-scale temperature and polarisation (EE) data, combined with ACT and SPT high-measurements, and taking the aforementioned correlations into account, Planck Collaboration (2016b) find a more constraining upper limit on the total kSZ signal D kSZ 3000 < 2.6 µK 2 with a 95% confidence level. Finally, adding new data from SPTpol 4 to their previous results (George et al. 2015), Reichardt et al. (2020) claim the first 3σ detection of the kSZ power spectrum with an amplitude D kS Z 3000 = 3.0±1.0 µK 2 , translated into a confidence interval on the patchy amplitude D pkS Z 3000 = 1.1 +1.0 −0.7 µK 2 using the models of homogeneous signal given in Shaw et al. (2012). These results are further pushed using the scaling relations derived by Battaglia et al. (2013) to obtain an upper limit on the duration of reionisation ∆z < 4.1.
Previous works have focused on relating the amplitude of the kSZ power spectrum at = 3000 to common reionisation parameters such as its duration and its midpoint. Battaglia et al. (2013) use large dark matter simulations (L 2 Gpc/h), post-processed to include reionisation, to construct light-cones of the kSZ signal and estimate its patchy power spectrum. The authors find the scalings D kSZ 3000 ∝z and D kSZ 3000 ∝ ∆z 0.51 wherez is approximately the midpoint of reionisation and here ∆z ≡ z (x e = 0.25)− z (x e = 0.75). Very large box sizes are necessary to capture the large-scale velocity flows contributing to the kSZ power spectrum at high-and results based on insufficiently large simulations will significantly underestimate the power at these scales. Shaw et al. (2012) find that a simulation box of side length 100 Mpc/h would miss 60% of the kSZ power at = 3000. For their own work, Shaw et al. (2012) therefore choose a completely different approach: they use hydrodynamical simulations to map the gas density to the dark matter power spectrum and later include this bias in a purely analytical derivation of the kSZ angular power spectrum. Because the non-linear dark matter power spectrum can be computed using the HALOFIT procedure (Smith et al. 2003) and because the velocity modes can be estimated fully from linear theory under a few assumptions, they avoid the limitations caused by simulation resolution and size mentioned above. With this method, the authors find a power-law dependence on both the reionisation midpoint z re and the optical depth τ for the homogeneous signal. For their most elaborate simulation, dubbed CSF, the cosmology-dependent scaling relations write D kS Z 3000 ∝ τ 0.44 and D kS Z 3000 ∝ z re 0.64 but are independent as one parameter is fixed before varying the other. The authors note that the current uncertainties on cosmological parameters such as σ 8 will wash out any potential constraint on z re and τ obtained from the measurement of the kSZ spectrum.
In this work, we choose to follow a similar approach. We build a comprehensive parameterisation allowing the full derivation of the kSZ angular power spectrum from a known reionisation history and morphology. In Sec. 2, we review the theoretical derivation of the kSZ power spectrum and propose a new parameterisation of the power spectrum of free electrons density contrast, based on the shape of the power spectrum of a bubble field. In Sec. 3, we present the simulations we later use to calibrate this parameterisation. In Sec. 4, we use the resulting expression of P ee (k, z) to compute the patchy kSZ angular power spectrum of our simulations and later apply the same procedure to different types of reionisation simulations. Finally, in Sec. 5, we discuss the physical meaning of our parameters and conclude. All distances are in comoving units and the cosmology used is the best-fit cosmology derived from Planck 2015 CMB data (Planck Collaboration 2016a): h = 0.6774, Ω m = 0.309, Ω b h 2 = 0.02230, Y p = 0.2453, σ 8 = 0.8164 and T CMB = 2.7255 K. Unless stated otherwise, P δδ describes the non-linear total matter power spectrum, x e (z) is the ratio of H ii and He ii ions to protons in the IGM, and the reionisation duration is defined by ∆z = z (x e = 0.25) − z (x e = 0.75).

Temperature fluctuations
The CMB temperature anisotropies coming from the scattering of CMB photons off clouds of free electrons with a non-zero bulk velocity v relative to the CMB rest-frame along the line of sightn write with σ T being the Thomson cross-section, c the speed of light, η the comoving distance to redshift z and v ·n the component of the peculiar velocity of the electrons along the line of sight. As mentioned before, τ is the Thomson optical depth, τ(z) = c σ T z electrons in the IGM are leftovers from recombination, is found to be negligible. We define q ≡ v(1 + δ e ) = v + vδ e ≡ v + q e the densityweighted peculiar velocity of the free electrons. It can be decomposed into a divergence-free q B and a curl-free q E components. We write their equivalents in the Fourier domain asq =q E +q B . As pointed out by e.g. Jaffe & Kamionkowski (1998), when projected along the line of sight,q E will cancel and only the component ofq perpendicular to k, i.e.q B , will contribute to the kSZ signal. We want an expression for the kSZ angular power spectrum C kSZ ≡ T 2 CMB |δ T kSZ (k)| 2 where k ≡ /η is the Limber wave-vector and is the multipole moment, which can be related to an angular scale in the sky. In the small angle limit, the kSZ angular power spectrum can be derived from Eq. (1) using the Limber approximation: with ∆ 2 B,e (k, z) ≡ k 3 P B,e (k, z)/(2π 2 ) and P B,e the power spectrum of the curl component of the momentum field defined by (2π) 3 P B,e δ D (k − k ) = q B,e (k)q * B,e (k ) where δ D is the Dirac delta function, the tilde denotes a Fourier transform and the asterisk a complex conjugate.
Expanding q B,eq * B,e , we obtain: where µ =k ·k , so that where 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. In the linear regime, we can write v(k) = ik ( fȧ/k)δ(k), where a is the scale factor and f the linear growth rate defined by f (a) = dlnD/dlna for D the growth function. With this we can compute the velocity power spectrum fully from linear theory and not be limited by the simulation size and resolution: where P lin δδ is the linear total matter power spectrum. We also assume for the cross-spectrum: where the bias b δe is defined by the ratio of the free electrons power spectrum over the non-linear matter power spectrum b δe (k, z) 2 = P ee (k, z)/P δδ (k, z). Although coarse, this approximation only has a minor impact on our results: it implies variations of ∼ 0.05 µK 2 in the power spectrum amplitude (see also Alvarez 2016). The final expression of the power spectrum of the curl component of the momentum field then writes which we can plug into Eq.
(2) to find the final expression for the kSZ angular power spectrum.

The power spectrum of free electrons density contrast
In Shaw et al. (2012), the authors choose to describe the behaviour of the free electrons power spectrum in terms of a biased matter power spectrum: they take P ee (k, z) ≡ b δe (k, z) 2 P δδ (k, z) and calibrate b δe (k, z) on their simulations, either extrapolating or assuming a reasonable behaviour for the scales and redshifts not covered by the simulations. However, because P ee describes the free electrons density fluctuations, it has a relatively simple structure, close to the power spectrum of a field made of ionised spheres on a neutral background, shown in Fig. 1, and using a bias is not necessary.
Consider a box of volume V = L 3 filled with n fully ionised bubbles of radius R, randomly distributed throughout the box so that their centres are located at a i for i ∈ {1, n}. The density of free electrons in the box follows where Θ (x) is the Heaviside step function,n e is the mean number density of electrons in the box and f the filling fraction of the box (here, f = x e ).n e / f is the number of electrons in one bubble divided by its volume and, ignoring overlaps, f = 4/3πR 3 n/V. Consider the electron density contrast field δ e on which P ee (k, z) is built: represented on Fig. 2 for one of the simulations used in this work. δ e Fourier-transforms intõ where W is the spherical top hat window function W(y) = (3/y 3 ) sin y − y cos y . Using this expression, and following Bharadwaj & Pandey (2005), the power spectrum of the electron density contrast field writes: which has units Mpc 3 . Fig. 1 gives an example of such a power spectrum. We have generated an ionisation field made of enough bubbles of radius R = 15 px = 5.5 Mpc 5 to reach a filling fraction f = 1% in a box of 512 3 pixels and side length L = 128/h Mpc. We compare the expression in Eq. (11) with power spectrum values computed directly from the 3D field and find a good match. On very small or very large scales, the window function behaves as: so that P ee (k) ∼ 4/3πR 3 / f i.e. is constant (see dashed horizontal line on the figure) on very large scales and has higher amplitude for smaller filling fractions. On small scales, the toy model power spectrum decreases as k −4 (see tilted dashed line on the figure). The intersection point of the horizontal and tilted dashed lines on the figure corresponds to k = 9 1/4 /R (dashed vertical line), hinting at a relation between the cut-off frequency and the bubble size. Note that Xu et al. (2019) find a similar feature, also related to the typical bubble size, in the bias between the H i and matter fields. This behaviour is close to what we observe when we compute the free electrons density power spectra of our six simulations, that can be seen at different stages in the reionisation process in Fig. 3. Therefore, we choose in this work to use a direct parameterisation of the scale and redshift evolution of P ee (k, z) during reionisation and calibrate it on our simulations. The parameters, α 0 and κ, are defined according to: .
In log-space, on large scales, P ee has a constant amplitude which, as mentioned above, depends on the filling fraction and therefore reaches its maximum α 0 at the start of the reionisation process, when the variance in the free electron field is maximal (see Sec. 5 The bubble radii actually follow a Gaussian distribution centred on 15 px with standard deviation 2 px. 5.1). It then slowly decreases as x e (z) −1/5 . Before the onset of reionisation, despite the few free electrons left over after recombination, the amplitude of P ee is negligible. This constant power decreases above a cut-off frequency that increases with time, following the growth of ionised bubbles, according to κx e (z) −1/3 . There is no power above this frequency, i.e. on smaller scales: there is no smaller ionised region than r min (z) = 2πx 1/3 e /κ at this time. For empirical reasons, we choose the power to decrease as k −3 and not k −4 as seen in the theoretical power spectrum on small scales. This difference can be explained by the fact that in our simulations, small ionised regions will keep appearing as new sources light up, maintaining power on scales smaller than the typical bubble size. Additionally, the density resolution will allow correlations between regions within a given bubble, whereas in the toy models ionised bubbles are only filled with ones. The complexity of the electron density contrast field is illustrated for one of the six simulations used in this work on Fig.  2: the underlying matter field is visible within the ionised regions.
Once reionisation is over and all IGM atoms are ionised, the fluctuations in free electrons density follow those of dark matter on large scales (k < 1 Mpc −1 ). On smaller scales, gas thermal pressure induces a drop in P ee (k, z) compared to the dark matter. To describe this evolution at low redshifts, we choose the same parameterisation as Shaw et al. (2012), given in Eq. (14), to describe the gas bias b δe (k, z) 2 = P ee (k, z)/P δδ (k, z) but adapt the parameters to our simulations, which however do not cover redshifts lower than 5.5: We find k f = 9.4 Mpc −1 and g = 0.5, constant with redshift.
Our values for k f and g are quite different from those obtained by Shaw et al. (2012), as in their work power starts dropping between 0.05 and 0.5 Mpc −1 compared to k ∼ 3 Mpc −1 for our simulations. This can be explained by our simulations making use of adaptive mesh refinement, therefore resolving very well the densest regions, so that our spectra are more sensitive to the thermal behaviour of gas. This model, where k f and g are constant parameters, is a very basic one. It will however be sufficient for this work since we focus on the patchy component of the kSZ effect, at z ≥ 5.5. Additionally, as we will show later, the scales mostly contributing to the patchy kSZ signal correspond to modes 10 −3 < k/Mpc −1 < 1 where P ee follows the matter power spectrum, so that a precise knowledge of b δe (k, z) is not required. In the future, if we want to apply our results to constrain reionisation with the measured CMB temperature power spectrum, we will need a better model as the observed signal will be the sum of homogeneous and patchy kSZ, with the former dominating on all scales. To account for the smooth transition of P ee from a power-law to a biased matter power spectrum, we write the final form for the free electrons density fluctuations power spectrum as for f H = 1 + Y p /4X p 1.08, with Y p and X p the primordial mass fraction of helium and hydrogen respectively. The total matter power spectrum P δδ is computed using the Boltzmann integrator CAMB (Lewis et al. 2000;Howlett et al. 2012) for the linear terms and the HALOFIT procedure for the non-linear contributions (Smith et al. 2003).

Description of the simulations
The simulations we use in this work were produced with the EMMA simulation code (Aubert et al. 2015) and previously used in Chardin et al. (2019). The code tracks the collisionless dynamics of dark matter, the hydrodynamics of baryons, star formation and feedback, and the radiative transfer using a moment-based method (see Aubert et al. 2018;Deparis et al. 2019). This code adheres to an Eulerian description, with fields described on grids, and enables adaptive mesh refinement techniques to increase the resolution in collapsing regions. Six simulations with identical numerical and physical parameters were produced in order to make up for the limited physical size of the box and the associated sample variance. They only differ in the random seeds used to generate the initial displacement phases, resulting in 6 different configurations of structures within the simulated volumes. Each run has a (128 Mpc/h) 3 volume sampled with 1024 3 cells at the coarsest level and 1024 3 dark matter particles. Refinement is triggered when the number of dark matter particles exceeds 8, up to 6 refinement levels. Initial conditions were produced using MUSIC (Hahn & Abel 2013) with a starting redshift of z = 150, assuming Planck Collaboration (2016a) cosmology. Simulations were stopped at z ∼ 6, before the full end of reionisation. The dark matter mass resolution is 2.1 × 10 8 M and the stellar mass resolution is 6.1 × 10 5 M . Star formation proceeds according to standard recipes described in Rasera & Teyssier (2006), with an overdensity threshold equal to 20 to trigger the gas-to-stellar particle conversion with a 0.1 efficiency: such values allow the first stellar particles to appear at z ∼ 17. Star particles produce ionising radiation for 3 Myr, with an emissivity provided by the Starburst99 model for a Top-Heavy initial mass function and a Z = 0.001 metallicity (Leitherer et al. 1999). Supernova feedback follows the prescription used in Aubert et al. (2018): as they reach an age of 15 million years, stellar particles dump 9.8 × 10 11 J per stellar kg in the surrounding gas, 1/3 in the form of thermal energy, 2/3 in the form of kinetic energy. Using these parameters, we obtain a cosmic star formation history consistent with constraints by Bouwens et al. (2015) and end up with 20 millions stellar particles at z = 6. The simulations were produced on the Occigen (CINES) and Jean-Zay (IDRIS) super- computers, using CPU architectures : a reduced speed of light of 0.1c has been used to reduce the cost of radiative transfer. Table 1 gives the midpoint z re and end of reionisation z end for each simulation, as well as the duration of the process, defined as the time elapsed between global ionisation fractions of 25% and of 75%. The upper panel of Fig. 5 shows the interpolated reionisation histories, where data points correspond to the snapshots available for each simulation. Note that the simulations stop before the reionisation process is fully over, therefore we need to extrapolate x e (z) to find the z end value given in Table 1. Additionally, the simulations originally do not include the first reionisation of helium. We correct for this by multiplying the IGM ionised fraction of hydrogen x H ii measured in the simulations by f H = 1 + Y p /4X p 1.08. Because we limit our study to redshifts z > 5.5, the second reionisation of helium is ignored. Fig. 2 shows the electron density contrast field for the first of our six simulations, close to the midpoint of reionisation. The complexity of the structure of this field is summarised in its power spectrum, visible in Fig. 3. The figure presents P ee (k, z) for each simulation, taken either at fixed redshift (first column) or fixed scale (right column). Despite identical numerical and physical parameters and very similar reionisation histories, the six simulations have different free electrons density power spectra, which translate into different kSZ power spectra.

Calibration procedure
We simultaneously fit the power spectra of the six simulations to Eq. (15) on a scale range 0.05 < k/Mpc −1 < 1.00 (20 bins), corresponding to the scales which contribute the most to the signal at = 3000 (see next paragraph), and a redshift range of 6.5 ≤ z ≤ 10.0 (10 bins), corresponding to the core of the reionisation process (0.07 < x e < 0.70). 6 We sample the parameter space of α 0 and κ on a regular grid (with spacings ∆ log α 0 = 0.001 and ∆κ = 0.0001) for which we compute the following likelihood: where {z i } and {k j } are the redshift and scale bins and the first sum is over the six simulations. Because our sample of six simulations is not sufficient to derive a meaningful covariance matrix, A&A proofs: manuscript no. main we choose to ignore correlations between scales across redshifts and use the diagonal of the covariance matrix to derive error bars σ e for each data point. We refer the interested reader to a discussion of this choice in Appendix A. We choose the best-fit as the duplet (α 0 , κ) for which the reduced χ 2 reaches its minimum value of 1.05 7 . The best-fit values, and their 68% confidence intervals are the following: We note a strong correlation between the two parameters due to both physical -see Sec. 5.1 -and analytical reasons: the value of κ impacts the low frequency amplitude of the P ee (k, z) model. The best-fit model, compared to the P ee (k, z) spectra of the six simulations Eq. (15) is fitted on, can be seen in Fig. 3 for three different redshift bins (left-hand column) and three different scale bins (right-hand column). Overall, we see a good agreement between the fit and the data points on the scales of interest, despite the simplicity of our model. Given the large number of P ee (k, z) data points (∼ 3500) and the complexity of the evolution of P ee with k and z, we must 7 The raw value is χ 2 ∼ 2500. limit our fits to given ranges. In order to assess what scales and redshifts contribute the most to the final kSZ signal, we look at the evolution of the integrand on z in Eq. (2) with time and at the evolution of the integral on k in Eq. (7) with scales. The results are shown in Fig. 4. The upper (resp. lower) panel presents the evolution of P ee (k, z) with scales (resp. redshift) after applying the fitting procedure described above. The width of each line represents the contribution of the redshift (resp. scale) of the corresponding colour to the final patchy kSZ amplitude at = 3000. We find that redshifts throughout reionisation contribute homogeneously to the signal, as 50% of the signal comes from redshifts z ≤ 7.2, i.e. slightly before the midpoint z re = 7.0. Redshifts on the range 6.5 < z < 8.5 contribute the most as they represent about 75% of the final kSZ power. Conversely, redshifts z > 10 contribute to only 0.4% of the total signal. On the lower panel, we see that scales outside the range 10 −3 Mpc −1 < k < 1 Mpc −1 contribute very marginally to the final signal (about 0.2%), whereas the range 10 −2 < k/Mpc −1 < 10 −1 makes up about 70% of D 3000 . Therefore, we choose to only keep data within the redshift range 6.5 < z < 9.0 (i.e. 7% < x e < 70%) and the scale range 10 −3 < k/Mpc −1 < 1 to constrain our fits. For reference, we compare the fit to data points at z = 7.8 (x e = 0.26) and k = 0.14 Mpc −1 for the first simulation, which show a good match overall.

Results on our six simulations
Now that we have a fitted P ee (k, z), we can compute the kSZ angular power spectrum using Eq. (2). We find: D p 3000 = 0.80 ± 0.06 µK 2 (18) and the angular scale at which the patchy angular spectrum reaches its maximum is max = 1800 +300 −100 . The angular patchy power spectrum is shown on the lower panel of Fig. 5. The error bars correspond to the propagated 68% confidence interval on the fit parameters. The amplitude of the homogeneous signal largely dominates that of the patchy signal, being about 4 times larger. The total kSZ amplitude reaches D 3000 = 4.2 µK 2 and so slightly exceeds the upper limits on the total kSZ amplitude given by SPT and Planck when SZxCIB correlations are allowed (resp. Reichardt et al. 2020;Planck Collaboration 2016b) but is however within the error bars of the ACT results (Sievers et al. 2013). With respect to the patchy signal, the amplitude is in perfect agreement with the claimed detection by the SPT at D patchy 3000 = 1.1 +1.0 −0.7 µK 2 (Reichardt et al. 2020), noting that our simulations reionise in a time very close to their constraint ∆z = 1.1 +1.6 −0.7 . The spectrum exhibits the expected bump in amplitude, here around ∼ 1800, i.e. larger scales than those found in other works (Iliev et al. 2007;Mesinger et al. 2012), hinting at larger ionised bubbles on average. Fig. 5 gives an idea of the variance in the kSZ angular power spectrum for given physics -in  15) to the P ee (k, z) data points from our six simulations (thick solid line) compared to the spectra obtained when interpolating the data points for each simulation (thin solid lines). Error bars correspond to the propagation of the 68% confidence interval on the fit parameters. The data point corresponds to constraints from Reichardt et al. (2020) at = 3000. particular a given matter distribution, and very similar reionisation histories: the distribution among simulations gives a reionisation midpoint defined at z re = 7.10 ± 0.06, corresponding to a range of kSZ power spectrum amplitude D 3000 = 0.80±0.06 µK 2 (at 68% confidence level), partly due to the limited size of our simulations (L = 128 Mpc/h). We compare in Fig. 5 the kSZ power spectrum resulting from fitting Eq. (15) on our six simulations simultaneously to the six spectra obtained when interpolating the P ee (k, z) data points available for each simulation: the six interpolated spectra lie withing the confidence limits of our best-fit.
Fixing the fit parameters to their most likely value for the fourth simulation, we artificially vary the reionisation history and compute the corresponding power spectrum. We successively fix the reionisation redshift but increase its duration ∆z or fix the duration but shift the midpoint z re . This corresponds to a scenario where the reionisation morphology is exactly the same, but happens later or earlier in time. We find clear scaling relations between the amplitude of the signal at = 3000, D 3000 , and both the reionisation duration ∆z and its midpoint z re . However, they are sensibly different from the results of Battaglia et al. (2013) as can be seen in Fig. 6. Even after rescaling to their z re = 8 and cosmology, we get a much lower amplitude. Note also that their patchy spectra bump around = 3000, whereas in our simulations the power has already dropped by = 3000 (Fig.  Fig. 6. Evolution of the amplitude of the patchy power spectrum at = 3000, D 3000 with the reionisation duration (upper panel) and the reionisation midpoint z re (lower panel), for different values of our parameters. Error bars correspond to the dispersion of kSZ amplitude at = 3000 (68% confidence level) propagated from errors on the fit parameters. The diamond data point corresponds to a seventh simulation, with reionisation happening earlier. In both panels, results are compared to those of Battaglia et al. (2013), rescaled to our cosmology. 5). When we vary κ and α 0 artificially, by fixing logα 0 = 3.54 instead of 3.70 as before, there is still a scaling relation, but both the slope and the intercept change. All of this demonstrates that the amplitude of the patchy signal largely depends on the physics of reionisation (here via the κ and α 0 parameters) and ∆z and z re are not sufficient to derive D 3000 . Simulations closer to those used in Battaglia et al. (2013) would likely give larger values for κ and α 0 , therefore increasing the amplitude to values closer to the authors' results. To confirm this, we generate a new simulation, with same resolution and box size but with twice as much star formation as in the six initial simulations, therefore reionising earlier (z re = 7.94) but on a similar redshift interval (∆z = 1.20). Applying the fitting procedure described above, we find logα 0 = 4.10 and κ = 0.08 Mpc −1 . The resulting patchy kSZ power spectrum can be seen in Fig. 7, along with the reionisation histories and the evolution of the typical bubble size r min = 2π/κx e (z) 1/3 . Results for this simulation are compared with what was obtained for our six simulations. The kSZ spectrum corresponding to an early reionisation scenario bumps at smaller scales ( max = 1400) with a much larger maximum amplitude (D max = 0.98 µK 2 ) but interestingly the amplitudes at = 3000 are similar. This suggests that focusing on D 3000 is not sufficient to characterise the kSZ signal.
These results corroborate the work of Park et al. (2013), who found that the scalings derived in Battaglia et al. (2013) are largely dependent on the simulations they were calibrated on, and therefore cannot be used as a universal formula to constrain reionisation. Notably, an asymmetric reionisation history x e (z) naturally deviates from this relation. Global parameters such as ∆z and z re are not sufficient to accurately describe the patchy kSZ signal, and one needs to take the physics of reionisation into account to get an accurate estimation of not only the shape, but also the amplitude of the power spectrum. Additionally, limiting ourselves to the amplitude at = 3000 to constrain reionisation can be misleading.

Tests on other simulations
We apply the same fitting procedure to the rsage simulation, described in Seiler et al. (2019), to test the robustness of our parameterisation. This simulation starts off as an N-body simulation (Seiler et al. 2018), containing 2400 3 dark matter particles within a 160 Mpc side box, resolving halos of mass ∼ 4 × 10 8 M with 32 particles. Galaxies are evolved over cosmic time following the Semi-Analytic Galaxy Evolution (SAGE) model of Croton et al. (2016), modified to include an improved modelling of galaxy evolution during the Epoch of Reionisation, including the feedback of ionisation on galaxy evolution. The semi-numerical code cifog (Hutter 2018a,b) is used to generate an inhomogeneous ultraviolet background (UVB) and follow the evolution of ionised hydrogen during the EoR. Three versions of the rsage simulation are used, each corresponding to a different way of modelling the escape fraction f esc of ionising photons from their host galaxy into the IGM. The first, dubbed rsage const, takes f esc constant and equal to 20%. The second, rsage fej, considers a positive scaling of f esc with f ej , the fraction of baryons that have been ejected from the galaxy compared to the number remaining as hot and cold gas. In the last one, rsage SFR, f esc scales with the star formation rate and thus roughly with the halo mass. Because they are based on the same dark matter distribution, the three simulations start reionising at similar times (z ∼ 13), but different source properties lead to different reionisation histories. In rsage SFR, the ionised bubbles are statistically larger than the other two simulations at a given redshift: this results into rsage SFR reaching 50% of ionisation at z re = 7.56 vs. z re = 7.45 and z re = 7.37 for rsage const and rsage fej respectively, and the full ionisation being achieved in a shorter time. For more details, we refer the interested reader to Seiler et al. (2019). We find that from the three different patchy kSZ power spectra presented in the upper panel of Fig. 8, we can clearly distinguish between the SFR model and the other two. The lower panel of the figure shows the evolution of r min (z) = 2π/κx e (z) 1/3 with redshift for the three models. We confirm the findings of Seiler et al. (2019), i.e. that an escape fraction scaling with the star formation history leads to larger ionised bubbles on average: the patchy power spectrum peaks for larger angular scales (smaller max ) for rsage SFR than for the other simulations and has larger r min values on average. Because it has larger bubbles, it also has a larger amplitude on average.
We further apply the fitting procedure to two 21CMFAST (Mesinger & Furlanetto 2007;Mesinger et al. 2011) simulations with dimensions L = 160 Mpc for 256 3 cells but different minimal halo mass. For M min = 10 8 M , the box is fully ionised by z end = 7.93 and the midpoint of reionisation is reached at z re = 9.45 for a process lasting ∆z = 2.17. For M min = 10 10 M , we find z end = 4.97, z re = 5.91 and ∆z = 1.34. We find that the parameterisation of Eq. (15) is again an accurate description of the evolution of the P ee (k, z) spectra of the two simulations. Because of the very noisy aspect of the ionisation field generated by 21CMFAST, with many small sources creating many small ionised regions at the start of the process, the best-fit κ value found for M min = 10 8 M (κ = 0.100 ± 0.008 Mpc −1 ) is larger than for our six high resolution simulations, reducing the minimal size of bubbles r min (z) at all times. Naturally, the resulting kSZ spectrum, seen in the top panel of Fig. 8, peaks at smaller angular scales. On the contrary, for M min = 10 10 M , because the minimal mass required to start ionising is larger, the ionising sources will be naturally more efficient and the ionised bubbles should be larger. Indeed, fitting Eq. (15) to the second 21CMFAST simulation, we find a smaller value of κ = 0.076 ± 0.014 Mpc −1 , corresponding to larger bubbles, as can be seen in the bottom panel of Fig. 8. With larger bubbles, we also expect the variance in the ionisation field at the start of the process to be higher than if many small ionised regions cover the neutral background. This corresponds to the larger value of log α 0 = 4.07 ± 0.25 found for the second simulation, compared to 3.63 ± 0.01 for the first one.
These results show that our proposed simple two-parameter expression for P ee (k, z) can accurately describe different types of simulations, i.e. different types of physics, further validating the physical interpretation of the parameters α 0 and κ detailed in the next Section.

Physical interpretation of the parameters: theoretical framework
Many previous works have empirically related the angular scale at which the patchy kSZ power spectrum reaches its maximum max to the typical size of bubbles during reionisation (Mc-Quinn et al. 2005;Iliev et al. 2007;Mesinger et al. 2012). To test this relation, we compute the patchy kSZ power spectrum for a given reionisation history x e (z) and α 0 but let κ values vary. We find a clear linear relation between κ and max as shown in Fig. 9. We see that previous results on the six high-resolution simulations and on 21CMFAST lie along this line, whereas the values found for the rsage simulations are slightly off. Such a relationship means that a detection of the patchy power spectrum in CMB observations would make it possible to directly estimate max , giving access to κ and consequently to the evolution of the typical bubble size with reasonable error bars. From Fig. 8, it looks like smaller ionised bubbles will lead to a higher amplitude for the kSZ spectrum, in contradiction with the suggestion made in Mesinger et al. (2012), an additional argument in favour of the necessity to consider the physics and morphology of reionisation to infer the amplitude of the patchy kSZ power spectrum. Additionally, we can link the theoretical expression of the large-scale amplitude of the bubble power spectrum in Eq. (11) with our parameterisation of P ee (k, z) in Eq. (13): α 0 x −1/5 e ↔ 4/3πR 3 /x e (z). Because of the simplicity of the toy model, this Article number, page 9 of 12 A&A proofs: manuscript no. main relation is not an equivalence. For example, contrary to the toy model, in our simulations, the locations of the different ionised bubbles are correlated, following the underlying dark matter distribution and this correlation will add power to the spectrum on large scales. This analogy can however explain the correlation observed between α 0 and κ when fitting Eq. (15) to data (recall that R ∝ 1/κ). Finally, since α 0 is independent of redshift, it will be a prefactor for the left-hand side of Eq. (7), therefore we expect a strong correlation between this parameter and the amplitude of the spectrum D 3000 , confirmed when we fix the reionisation history and κ but vary α 0 on the range 3.0 < log α 0 < 4.4. The shape of the kSZ spectrum, however, does not change in this case. This further confirms that the morphology of reionisation, through the variance of the ionisation field, has a strong impact on the amplitude of the patchy kSZ power spectrum.

Conclusions & prospects
In this work, we have used state-of-the-art reionisation simulations (Aubert et al. 2015) to calibrate an analytical expression of the angular power spectrum of the kSZ effect stemming from patchy reionisation. We have shown that describing the shape, but also amplitude of the signal only in terms of global parameters such as the reionisation duration ∆z and its midpoint z re is not sufficient: it is essential to take the physics of the process into account. In our new proposed expression, the parameters can be directly related to both the global reionisation history x e (z) and to the morphology of the process, through the typical size of ionised regions and the variance of the ionisation field. With as few as these three parameters, we can fully recover the patchy kSZ angular power spectrum. Conversely, a precise measurement of the spectrum will constrain these parameters and provide us with detailed information about the physics of reionisation. This result sheds light on the scaling relations observed in previous works, for example between the amplitude of the spectrum at = 3000 and the reionisation midpoint (Battaglia et al. 2013) or between the scale at which the spectrum peaks and the typical size of ionised bubbles (Iliev et al. 2007). We have tested the robustness of our model by confronting it to different types of simulations, capturing different aspects of the process, namely the rsage simulation (Seiler et al. 2018), and by comparing our results to ways of deriving the kSZ power spectrum which are closer to data, such as interpolating the different spectra used in the computation.
On average, our results are in good agreement with previous works, despite a low amplitude of the patchy kSZ angular power spectrum at = 3000 (∼ 0.80 µK 2 ). There is undoubtedly a bump around scales ∼ 2000 that can be related to the typical bubble size and the amplitude of the total (patchy) kinetic SZ spectrum ranges from 4 to 5 µK 2 (0.5 to 1.5 µK 2 , respectively) for plausible reionisation scenarios, i.e. lies within the error bars of the latest observational results of ACT (Sievers et al. 2013) and SPT (Reichardt et al. 2020). We have found that the majority of the patchy kSZ signal stems from scales 10 −3 < k/Mpc −1 < 1 and from the core of the reionisation process (10% < x e < 80%), ranges on which we must focus our efforts to obtain an accurate description. However, this analysis does not consider third-and fourth-order components of the kSZ signal, which can represent as much as 10% of the total signal (Alvarez 2016), and uses a coarse approximation for the electrons density -velocity cross spectra.
In future work, we will use the parameterisation proposed in this work to constrain reionisation with both current and future CMB observations from SPT and CMB-S4 experiments. In order to constrain α 0 and κ, these measurements could be combined with ionisation maps from 21cm experiments such as the Square Kilometre Array (SKA, Mellema et al. 2015). The main challenge remains to separate first the kSZ signal from other foregrounds, and then the patchy kSZ signal from the homogeneous one. To solve the first part of this problem, Calabrese et al. (2014) suggest to subtract the theoretical primary power spectrum (derived from independent cosmological parameter constraints obtained from polarisation measurements) from the observed one so that the signal left is the kSZ power spectrum alone. Secondly, one would need a good description of the homogeneous spectrum, similar to the results of Shaw et al. (2012) but updated with more recent simulations, in order to estimate how accurately one can recover the patchy signal.