© A. Gorce et al. 2020

1. 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 I 2020), 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 that 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 a redshift of 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, that is, a redshift-evolution for the IGM global ionised fraction x_e(z). In standard Boltzmann solvers which are used to compute theoretical predictions in CMB data analysis such as the CAMB code (Lewis et al. 2000; Howlett et al. 2012)¹, the reionisation scenario used is a step-like 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 I 2016). However, this parameterisation does not match simulations and observations well since we expect the ionisation fraction to slowly rise when the first sources light up, before taking off as soon as 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. Int. XLVII 2016). 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 rest-frame in a process called the 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, that is smaller than about 5 arcmin), 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, which come from the Doppler shifting of photons on free electrons that are 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 (McQuinn et al. 2005; Mesinger et al. 2012). For example, the patchy signal is expected to peak around ℓ ∼ 2000, corresponding to the typical bubble size during reionisation (Zahn et al. 2005; 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 ground-based Atacama cosmology telescope (ACT)² and the South Pole telescope (SPT)³ 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 ${\mathcal{D}}_{3000}^{\mathrm{SZ}}\equiv \ell (\ell +1)C_{\ell =3000}^{\mathrm{SZ}}/2\pi =6.8\pm 2.9\mu {\mathrm{K}}^2$ at the 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 ${\mathcal{D}}_{3000}^{\mathrm{kSZ}}<2.8\mu {\mathrm{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 ${\mathcal{D}}_{3000}^{\mathrm{patchy}}\le 2.1\mu {\mathrm{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. Int. XLVII (2016) find a more constraining upper limit on the total kSZ signal ${\mathcal{D}}_{3000}^{\mathrm{kSZ}}<2.6\mu {\mathrm{K}}^2$ with a 95% confidence level. Finally, adding new data from SPTpol⁴ 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 ${\mathcal{D}}_{3000}^{\mathrm{kSZ}}=3.0\pm 1.0\mu {\mathrm{K}}^2$, translated into a confidence interval on the patchy amplitude ${\mathcal{D}}_{3000}^{\mathrm{pkSZ}}=1.1_{-0.7}^{+1.0}\mu {\mathrm{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 ${\mathcal{D}}_{3000}^{\mathrm{kSZ}}\propto \overline{z}$ and ${\mathcal{D}}_{3000}^{\mathrm{kSZ}}\propto \mathrm{\Delta }{z}^{0.51}$ where $\overline{z}$ 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 about 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 ${\mathcal{D}}_{3000}^{\mathrm{kSZ}}\propto {\tau }^{0.44}$ and ${\mathcal{D}}_{3000}^{\mathrm{kSZ}}\propto {z_{\mathrm{re}}}^{0.64}$ but are independent since one parameter is fixed before varying the other. The authors note that the current uncertainties on cosmological parameters such as σ₈ 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 Sect. 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 Sect. 3, we present the simulations we later use to calibrate this parameterisation. In Sect. 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 Sect. 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 I 2016): h = 0.6774, Ω_m = 0.309, Ω_bh² = 0.02230, Y_p = 0.2453, σ₈ = 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). The code used to compute the kSZ power spectrum can be found online⁵.

2. Derivation of the kSZ angular power spectrum

2.1. 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 sight $\widehat{\mathit{n}}$ write

$$\begin{array}{c}\hfill \delta T_{k\mathrm{S}Z}(\widehat{\mathit{n}})=\frac{{\sigma }_T}{c}\int \frac{\mathrm{d}\eta }{\mathrm{d}z}\frac{\mathrm{d}z}{(1+z)}{\mathrm{e}}^{-\tau (z)}{n}_e(z)\mathit{v}·\widehat{\mathit{n}},\end{array}$$(1)

with σ_T being the Thomson cross-section, c the speed of light, η the comoving distance to redshift z and $\mathit{v}·\widehat{\mathit{n}}$ the component of the peculiar velocity of the electrons along the line of sight. As mentioned before, τ is the Thomson optical depth, $\tau (z)=c{\sigma }_{\mathrm{T}}{\int }_0^z{n}_e(z\prime )/H(z\prime ){(1+z\prime )}^2\mathrm{d}z\prime$. n_e is the mean free electrons number density at redshift z from which we derive the density contrast δ_e via $n_e={\overline{n}}_e(1+{\delta }_e)$. We choose the limits of the integral in Eq. (1) depending on the type of signal we are interested in: for homogeneous kSZ, we integrate from 0 to z_end, the redshift when reionisation ends; for patchy kSZ, the main focus of this work, we integrate from z_end to the highest redshift considered in the simulation (here, z_max = 15). The contribution from redshifts larger than the onset of reionisation, when the only free 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 density-weighted 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 as $\stackrel{\mathbf{\sim }}{\mathit{q}}={\stackrel{\mathbf{\sim }}{\mathit{q}}}_E+{\stackrel{\mathbf{\sim }}{\mathit{q}}}_B$. As pointed out by Jaffe & Kamionkowski (1998), when projected along the line of sight, ${\stackrel{\mathbf{\sim }}{\mathit{q}}}_E$ will cancel and only the component of $\stackrel{\mathbf{\sim }}{\mathit{q}}$ perpendicular to k, that is ${\stackrel{\mathbf{\sim }}{\mathit{q}}}_B$, will contribute to the kSZ signal. We want an expression for the kSZ angular power spectrum $C_{\ell }^{k\mathrm{SZ}}\equiv T_{\mathrm{CMB}}^2{|{\stackrel{\sim }{\delta T}}_{\mathrm{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:

$$\begin{array}{c}\hfill C_{\ell }=\frac{8{\pi }^2}{{(2\ell +1)}^3}\frac{{\sigma }_T^2}{c^2}{\displaystyle \int }\frac{{\overline{n}}_e{(z)}^2}{{(1+z)}^2}{\mathrm{\Delta }}_{B,e}^2(\ell /\eta ,z){\mathrm{e}}^{-2\tau (z)}\eta \frac{\mathrm{d}\eta }{\mathrm{d}z}\mathrm{d}z,\end{array}$$(2)

with ${\mathrm{\Delta }}_{B,e}^2(k,z)\equiv k^3{P}_{B,e}(k,z)/(2{\pi }^2)$ and P_B, e the power spectrum of the curl component of the momentum field defined by ${(2\pi )}^3{P}_{B,e}{\delta }_D(\mathit{k}-\mathit{k}\prime )=⟨{\stackrel{\mathbf{\sim }}{\mathit{q}}}_{B,e}(\mathit{k}){\stackrel{\mathbf{\sim }}{\mathit{q}}}_{B,e}^{\ast }(\mathit{k}\prime )⟩$ where δ_D is the Dirac delta function, the tilde denotes a Fourier transform and the asterisk a complex conjugate.

Expanding $⟨{\stackrel{\mathbf{\sim }}{\mathit{q}}}_{B,e}{\stackrel{\mathbf{\sim }}{\mathit{q}}}_{B,e}^{\ast }⟩$, we obtain:

$$\begin{array}{c}\hfill {\stackrel{\mathbf{\sim }}{\mathit{q}}}_{B,e}(\mathit{k})={\displaystyle \int }\frac{{\mathrm{d}}^3\mathit{k}\prime }{{(2\pi )}^3}(\widehat{\mathit{k}}\prime -\mu \widehat{\mathit{k}})\stackrel{\sim }{v}(k\prime ){\stackrel{\sim }{\delta }}_e(|\mathit{k}-\mathit{k}\prime |),\end{array}$$(3)

where $\mu =\widehat{\mathit{k}}·\widehat{\mathit{k}}\prime$, so that

$$\begin{array}{cc}\hfill \frac{⟨{\stackrel{\mathbf{\sim }}{\mathit{q}}}_{B,e}(\mathit{k}){\stackrel{\mathbf{\sim }}{\mathit{q}}}_{B,e}^{\ast }(\mathit{k}\prime )⟩}{{(2\pi )}^3{\delta }_D(|\mathit{k}-\mathit{k}\prime |)}& \equiv \frac{2{\pi }^2}{k^3}{\mathrm{\Delta }}_{B,e}^2(k,z)\hfill \\ \hfill & =\frac{1}{{(2\pi )}^3}{\displaystyle \int }{\mathrm{d}}^{3}k\prime [(1-{\mu }^2)P_{\mathit{ee}}(|\mathit{k}-\mathit{k}\prime |)P_{\mathit{vv}}(k\prime )\hfill \\ \hfill & -\frac{(1-{\mu }^2)k\prime }{|\mathit{k}-\mathit{k}\prime |}{P}_{\mathit{ev}}(|\mathit{k}-\mathit{k}\prime |)P_{\mathit{ev}}(k\prime )],\hfill \end{array}$$(4)

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 $\mathit{v}(\mathit{k})=i\mathit{k}(f\dot{a}/k)\stackrel{\sim }{\delta }(\mathit{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:

$$\begin{array}{c}\hfill P_{\mathit{vv}}(k,z)={\left(\frac{\dot{a}f(z)}{k}\right)}^2{P}_{\delta \delta }^{\mathrm{lin}}(k,z)\end{array}$$(5)

where $P_{\delta \delta }^{\mathrm{lin}}$ is the linear total matter power spectrum. We also assume for the cross-spectrum:

$$\begin{array}{c}\hfill P_{\mathit{ve}}(k,z)\simeq b_{\delta e}(k,z)P_{\delta v}(k,z)=\frac{f\dot{a}(z)}{k}{b}_{\delta e}(k,z)P_{\delta \delta }^{\mathrm{lin}}(k,z),\end{array}$$(6)

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)² = 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² 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

$$\begin{array}{cc}\hfill P_{B,e}(k,z)& =\frac{1}{{(2\pi )}^3}f{(z)}^2\dot{a}{(z)}^2{\displaystyle \int }{\mathrm{d}}^{3}k\prime (1-{\mu }^2)\times \hfill \\ \hfill & [\frac{1}{k^{\prime 2}}{P}_{\mathit{ee}}(|\mathit{k}-\mathit{k}\prime |)P_{\delta \delta }^{\mathrm{lin}}(k\prime ,z)\hfill \\ \hfill & -\frac{b_{\delta e}(k\prime ,z)}{{|\mathit{k}-\mathit{k}\prime |}^2}{b}_{\delta e}(|\mathit{k}-\mathit{k}\prime |,z)P_{\delta \delta }^{\mathrm{lin}}(|\mathit{k}-\mathit{k}\prime |,z)P_{\delta \delta }^{\mathrm{lin}}(k\prime ,z)],\hfill \end{array}$$(7)

which we can plug into Eq. (2) to find the final expression for the kSZ angular power spectrum.

2.2. 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)²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.

Fig. 1. Free electrons density contrast power spectrum for a box filled with enough bubbles of radius R = 15 px = 5.5 Mpc to reach a filling fraction f = 1%. Points are results of a numerical computation of the power spectrum, compared to the theoretical model (solid line). The dotted vertical line corresponds to k = 1/R, the dashed vertical line to 91/4/R, the dashed horizontal line to 4/3πR3/f and the tilted dashed line has slope k−4.

Consider a box of volume V = L³ 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

$$\begin{array}{c}\hfill n_e(\mathit{r})=\frac{{\overline{n}}_e}{f}{\displaystyle \sum _{i=1}^n}\mathrm{\Theta }\left(\frac{|\mathit{r}-{\mathit{a}}_i|}{R}\right),\end{array}$$(8)

where Θ(x) is the Heaviside step function, ${\overline{n}}_e$ is the mean number density of electrons in the box and f the filling fraction of the box (here, f = x_e). ${\overline{n}}_e/f$ is the number of electrons in one bubble divided by its volume and, ignoring overlaps, f = 4/3πR³n/V. Consider the electron density contrast field δ_e on which P_ee(k, z) is built:

$$\begin{array}{c}\hfill {\delta }_e(\mathit{r})=\frac{n_e(\mathit{r})}{{\overline{n}}_e}-1=\frac{1}{f}{\displaystyle \sum _{i=1}^n}\mathrm{\Theta }\left(\frac{|\mathit{r}-{\mathit{a}}_i|}{R}\right)-1,\end{array}$$(9)

represented on Fig. 2 for one of the simulations used in this work. δ_e(r) Fourier–transforms into

$$\begin{array}{c}\hfill {\stackrel{\sim }{\delta }}_e(\mathit{k})=\frac{L^3}{n}W(kR){\displaystyle \sum _{i=1}^n}{\mathrm{e}}^{-i\mathit{k}·{\mathit{a}}_i},\end{array}$$(10)

Fig. 2. Left panel: snapshot of the electron density contrast field for the first of the six simulations, at z = 7.2 and xe = 0.49. Right panel: free electrons power spectrum of the same simulation at fixed redshifts (fixed ionised levels). The shaded area corresponds to scales contributing ${\mathcal{D}}_{3000}^{\mathrm{patchy}}$ the most (see Sect. 3.2) and the solid black line to the field shown in the left panel.

where W is the spherical top-hat window function W(y) = (3/y³)[sin y−y cos y]. Using this expression, and following Bharadwaj & Pandey (2005), the power spectrum of the electron density contrast field writes:

$$\begin{array}{c}\hfill P_{\mathit{ee}}(\mathit{k})=\frac{4}{3}\pi R^3\frac{1}{f}{W}^2(kR),\end{array}$$(11)

which has units Mpc³. Figure 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⁶ to reach a filling fraction f = 1% in a box of 512³ 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:

$$\begin{array}{c}\hfill \begin{array}{ccc}\hfill W(y)& \sim \frac{3}{y^3}\times \frac{y^3}{3}=1\hfill & \hfill \mathrm{as}y\to 0\\ \hfill W(y)& \sim \frac{3}{y^3}\times y=\frac{3}{y^2}\hfill & \hfill \mathrm{as}y\to \infty \end{array}\end{array}$$(12)

so that P_ee(k)∼4/3πR³/f 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⁻⁴ (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. Interestingly, 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 in the free electrons density power spectra of the custom set of simulations used in this work in the early stages of reionisation, as can be seen on the right panel of Fig. 2. 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, α₀ and κ, are defined according to:

$$\begin{array}{c}\hfill P_{\mathit{ee}}(k,z)=\frac{{\alpha }_0{x}_e{(z)}^{-1/5}}{1+{[k/\kappa ]}^3{x}_e(z)}·\end{array}$$(13)

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 α₀ at the start of the reionisation process, when the variance in the free electron field is maximal (see Sect. 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, that is on smaller scales: there is no smaller ionised region than $r_{\text{min}}(z)=2\pi x_e^{1/3}/\kappa$ at this time. For empirical reasons, we choose the power to decrease as k⁻³ and not k⁻⁴ 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⁻¹). 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)² = P_ee(k, z)/P_δδ(k, z) but adapt the parameters to our simulations, which however do not cover redshifts lower than 5.5:

$$\begin{array}{c}\hfill b_{\delta e}{(k,z)}^2=\frac{1}{2}[{\mathrm{e}}^{-k/k_f}+\frac{1}{1+{(gk/k_f)}^{7/2}}].\end{array}$$(14)

We find k_f = 9.4 Mpc⁻¹ 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⁻¹ compared to k ∼ 3 Mpc⁻¹ 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 shown later, the scales mostly contributing to the patchy kSZ signal correspond to modes 10⁻³ <  k Mpc⁻¹ <  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, illustrated in the right panel of Fig. 2, we write the final form for the free electrons density fluctuations power spectrum as

$$\begin{array}{cc}\hfill P_{\mathit{ee}}(k,z)=[f_{\mathrm{H}}-x_e(z)]& \times \frac{{\alpha }_0{x}_e{(z)}^{-1/5}}{1+{[k/\kappa ]}^3{x}_e(z)}\hfill \\ \hfill & +x_e(z)\times b_{\delta e}{(k,z)}^2{P}_{\delta \delta }(k,z),\hfill \end{array}$$(15)

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

3. Calibration on simulations

3.1. 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)³ volume sampled with 1024³ cells at the coarsest level and 1024³ 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 I (2016) cosmology. Simulations were stopped at z ∼ 6, before the full end of reionisation. The dark matter mass resolution is 2.1 × 10⁸ M_⊙ and the stellar mass resolution is 6.1 × 10⁵ 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¹¹ 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) supercomputers, 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. Originally, our simulations 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. Figure 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, shown in the right panel. Figure 3 compares the P_ee(k, z) spectra of the six simulations, 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 translates into different kSZ power spectra.

Fig. 3. Result of the fit of Eq. (15) on the free electrons power spectrum of our six simulations, for three redshift bins (left panels) and three scale bins (right panels). The best-fit is shown as the thick black line with the accompanying 68% confidence interval, and the spectra of the six simulations as thin coloured lines. Error bars on data points are computed from the covariance matrix (see text for details).

3.2. Calibration procedure

We simultaneously fit the power spectra of the six simulations to Eq. (15) on a scale range 0.05 <  k/Mpc⁻¹ <  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)⁸. We sample the parameter space of α₀ and κ on a regular grid (with spacings Δ log α₀ = 0.001 and Δκ = 0.0001) for which we compute the following likelihood:

$$\begin{array}{c}\hfill {\chi }^2={\displaystyle \sum _{\mathrm{n}=1}^6}{\displaystyle \sum _{z_i}\sum _{k_j}}\frac{1}{{\sigma }_e^2}{[P_{\mathit{ee}}^{\mathrm{data}}(k_j,z_i)-P_{\mathit{ee}}^{\mathrm{model}}(k_j,z_i)]}^2,\end{array}$$(16)

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, 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 (α₀, κ) for which the reduced χ² reaches its minimum value of 1.05⁹. The best-fit values, with their 68% confidence intervals are

$$\begin{array}{c}\hfill \begin{array}{cc}& log{\alpha }_0/{\mathrm{Mpc}}^3=3.{93}_{-0.06}^{+0.05}\hfill \\ \hfill & \kappa =0.{084}_{-0.004}^{+0.003}{\mathrm{Mpc}}^{-1}.\hfill \end{array}\end{array}$$(17)

We note a strong correlation between the two parameters due to both physical – see Sect. 5.1 – and analytical reasons. Indeed, 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 originally (∼3500) and the complexity of the evolution of P_ee with k and z, we must 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 left (resp. right) upper 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. The lower panels present the corresponding probability density and cumulative distribution functions. We find that redshifts throughout reionisation contribute homogeneously to the signal, since 50% stems from redshifts z ≤ 7.2, 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⁻³ Mpc⁻¹ <  k <  1 Mpc⁻¹ contribute very marginally to the final signal (about 0.2%), whereas the range 10⁻² <  k/Mpc⁻¹ <  10⁻¹ makes up about 70% of 𝒟₃₀₀₀. Therefore, we choose to only keep data within the redshift range 6.5 <  z <  10.0 (i.e. 7%< x_e <  70%) and the scale range 10⁻³ <  k/Mpc⁻¹ <  1 to constrain our fits. For reference, on Fig. 4, we compare the fit to data points at z = 7.8 (x_e = 0.26) and k = 0.14 Mpc⁻¹ for the first simulation, and find an overall good match.

Fig. 4. Upper panels: Pee(k, z) as fitted on spectra from the fourth simulation, as a function of scales (left panel) and of redshift (right panel). For reference, the fit is compared to data points for z = 7.8 (xe = 0.26) and k = 0.14 Mpc−1 respectively, with corresponding colour. The width of each line represents the contribution of the redshift (resp. scale) of the corresponding scale (resp. redshift) to the final patchy kSZ amplitude at ℓ = 3000. Lower panels: corresponding probability densities (dashed lines) and cumulative distributions (solid lines). Shaded areas correspond to the first 50% of the signal. The dotted vertical line on the lower right panel marks the midpoint of reionisation.

4. Propagation to the kSZ power spectrum

4.1. 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:

$$\begin{array}{c}\hfill {\mathcal{D}}_{3000}^{\mathrm{p}}=0.80\pm 0.06\mu {\mathrm{K}}^2\end{array}$$(18)

and the angular scale at which the patchy angular spectrum reaches its maximum is ${\ell }^{\max }={1800}_{-100}^{+300}$. 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 𝒟₃₀₀₀ = 4.2 μK² 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. Int. XLVII 2016) 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 ${\mathcal{D}}_{3000}^{\mathrm{patchy}}=1.1_{-0.7}^{+1.0}\mu {\mathrm{K}}^2$ (Reichardt et al. 2020), noting that our simulations reionise in a time very close to their constraint $\mathrm{\Delta }z=1.1_{-0.7}^{+1.6}$. The spectrum exhibits the expected bump in amplitude, here around ℓ ∼ 1800, corresponding to larger scales than those found in other works (Iliev et al. 2007; Mesinger et al. 2012), hinting at larger ionised bubbles on average. Figure 5 gives an idea of the variance in the kSZ angular power spectrum for given physics – in 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 𝒟₃₀₀₀ = 0.80 ± 0.06 μK² (at 68% confidence level). Part of this variance can be related to sample variance, since our simulations have a too small side length (L = 128 Mpc h⁻¹) to avoid it (Iliev et al. 2007). 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.

Fig. 5. Results for our six simulations. Upper panel: global reionisation histories, for H II and He II. The dotted horizontal line marks the reionisation midpoint zre. Lower panel: angular kSZ power spectrum after fitting Eq. (15) to the Pee(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.

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, 𝒟₃₀₀₀, 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. 5), hinting at a very different reionisation morphology from ours. When we vary κ and α₀ artificially, by fixing log α₀ = 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 α₀ parameters) and Δz and z_re are not sufficient to derive 𝒟₃₀₀₀. Simulations closer to those used in Battaglia et al. (2013) would likely give larger values for κ and α₀, 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 α₀ = 4.10 and κ = 0.08 Mpc⁻¹. 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 larger scales (ℓ^max = 1400) with a much larger maximum amplitude (𝒟_max = 0.98 μK²) but interestingly the amplitudes at ℓ = 3000 are similar. This suggests that focusing on 𝒟₃₀₀₀ is not sufficient to characterise the kSZ signal.

Fig. 6. Evolution of the amplitude of the patchy power spectrum at ℓ = 3000, 𝒟3000 with the reionisation duration (upper panel) and the reionisation midpoint zre (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.

Fig. 7. Comparison of results for our six initial simulations, corresponding to a late reionisation scenario, and for an additional seventh simulation, corresponding to an early reionisation scenario. Left panel: reionisation histories. Middle panel: patchy kSZ angular power spectra. The data point corresponds to constraints from Reichardt et al. (2020). Right panel: minimal size of ionised regions as a function of global ionised level. Shaded areas correspond to the 68% confidence level on kSZ amplitude propagated from the probability distributions of the fit parameters.

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.

4.2. Tests on other simulations

We now look at 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³ dark matter particles within a 160 Mpc side box, resolving halos of mass ∼4 × 10⁸ 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, shown in the left upper panel of Fig. 8. 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). Applying the fitting procedure to the three simulations, we find that the parameterisation of Eq. (15) is an accurate description of the evolution of their P_ee(k, z) spectra (detailed fit results are given in Appendix B.2). Resulting patchy kSZ angular power spectra are shown in the upper middle panel of Fig. 8. First, we find that rsage fej has the smallest α₀ value, with log α₀ = 2.87 ± 0.04. Because α₀ is the maximum amplitude of the P_ee(k, z) spectrum, built upon the free electrons density contrast field ${\delta }_e(r)=n_e(r)/{\overline{n}}_e-1$, it will scale with the variance of the n_e(r) field. Therefore a smaller α₀ value is equivalent to a smaller field variance at all times. This is consistent with the picture of the different rsage simulations we have: as presented in Seiler et al. (2019), rsage fej exhibits the smallest ionised bubbles on average. For a given filling fraction, a field made of many small bubbles covering the neutral background rather homogeneously will have smaller variance than one made of a few large bubbles. This in turn explains why rsage SFR gives the largest α₀ value (log α₀ = 3.47 ± 0.04), and, later, the largest kSZ amplitude (Fig. 8). Second, the rsage SFR simulation has the smallest value of κ (κ = 0.123 ± 0.004 Mpc⁻¹): the upper right panel of Fig. 8 shows the evolution of $r_{\text{min}}=2\pi /\kappa x_e^{1/3}$ with ionisation level for the three models. Because rsage SFR has the largest ionised bubbles on average (Seiler et al. 2019), this result confirms the interpretation of 1/κ as an estimate of the typical bubble size during reionisation. Additionally, the patchy power spectrum derived from rsage SFR peaks at larger angular scales (ℓ^max ∼ 2400) than for the other simulations, as can be seen in the upper middle panel of the figure. Interestingly, the largest α₀ value leads to the strongest kSZ signal and the smallest κ value to the spectrum whose bump is observed on the largest scales (the smallest ℓ^max). We investigate these potential links in the next section.

Fig. 8. Comparison of results for the three rsage simulations (upper panels) and the three 21CMFAST runs (lower panels) considered. Left panels: reionisation histories. Middle panels: patchy kSZ angular power spectra. The data point corresponds to constraints from Reichardt et al. (2020). Right panels: minimal size of ionised regions as a function of global ionised level. The shaded areas correspond to the 68% confidence interval propagated from the 68% confidence intervals on the fit parameters.

We now turn to three 21CMFAST (Mesinger & Furlanetto 2007; Mesinger et al. 2011) simulations with dimensions L = 160 Mpc for 256³ cells (same box size and resolution as rsage). Between the three runs, we vary the parameter M_turn, the turnover mass, which corresponds to the minimum halo mass before exponential suppression of star formation. For M_turn = 10⁸ M_⊙, the box is fully ionised by z_end = 6.25 and the midpoint of reionisation is reached at z_re = 8.92 for a process lasting Δz = 1.91. For M_turn = 10⁹ M_⊙, we find z_end = 4.69, z_re = 7.11 and Δz = 1.66, which is closest to rsage and our initial six simulations. Finally, M_turn = 10¹⁰ M_⊙ yields z_end = 3.37, z_re = 5.41 and Δz = 1.47. Indeed, the point of these simulations is not only to test the sensitivity of our approach to astrophysical parameters, but also to see the impact of very different reionisation histories on the patchy kSZ power. We find that Eq. (15) again nicely fits the evolution of the P_ee(k, z) spectra of these simulations, as shown in Appendix B.1. The resulting reionisation histories, patchy kSZ spectra and r_min(x_e) are shown in the lower panels of Fig. 8. For M_turn = 10⁸ M_⊙, many small-mass halos are active ionising sources, resulting in an ionising field made of many small bubbles at the start of the process. This translates into this simulation having the largest best-fit κ value of the three (κ = 0.130  ±  0.003 Mpc⁻¹) and so the smallest r_min(x_e). Naturally, the resulting kSZ spectrum peaks at smaller angular scales. For the other extreme case M_turn = 10¹⁰ M_⊙, because the minimal mass required to start ionising is larger, the ionising sources are more efficient and the ionised bubbles larger. Indeed, we find a smaller value of κ = 0.093  ±  0.003 Mpc⁻¹. 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 α₀ = 3.79 ± 0.04 found for this simulation, compared to 3.30 ± 0.03 for the first one. However, this larger value of α₀ this time does not result into the strongest kSZ signal because of the very different reionisation histories of the three simulations. As we have seen in the previous section, the amplitude of the signal is strongly correlated with the duration and midpoint of reionisation, resulting in the first simulation (M_turn = 10⁸ M_⊙), corresponding to the earliest reionisation, having the strongest signal. This again emphasises how essential it is to consider both reionisation morphology and global history to derive the final kSZ spectrum.

These results show that our proposed simple two-parameter expression for P_ee(k, z) can accurately describe different types of simulations, that is different types of physics, further validating the physical interpretation of the parameters α₀ and κ detailed in the next Section.

5. Discussion and conclusions

5.1. Physical interpretation of the parameters

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 (McQuinn et al. 2005; Iliev et al. 2007; Mesinger et al. 2012). To test for this relation, we compute the patchy kSZ power spectrum for a given reionisation history x_e(z) and α₀ but let κ values vary. We find a clear linear relation between κ and ℓ^max as shown in Fig. 9. Despite very different reionisation histories and physics at stake, previous results on the six high-resolution simulations, on 21CMFAST, and on rsage, roughly lie along this line. This means that a detection of the patchy power spectrum in CMB observations would make it possible to directly estimate ℓ^max, giving access to κ without bias from reionisation history, and to the evolution of the typical bubble size. As the growth of ionised regions depends on the physical properties of early galaxies, such as their ionising efficiency or their star formation rate and on the density of the IGM, constraints on κ could, in turn, give constraints on the nature of early light sources and of the early IGM.

Fig. 9. Evolution of the peaking angular scale of the patchy kSZ power spectrum for one given reionisation history but different values of the κ parameter. The red dotted line is the result of a linear regression. Inferences are compared to results for different simulations.

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): ${\alpha }_0{x}_e^{-1/5}\leftrightarrow 4/3\pi R^3/x_e(z)$. Because of the simplicity of the toy model, this 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 α₀ and κ when fitting Eq. (15) to data (recall that R ∝ 1/κ). Finally, since α₀ is independent of redshift, it will be a pre-factor for the left-hand side of Eq. (7), therefore we expect a strong correlation between this parameter and the amplitude of the spectrum at ℓ = 3000 and with the maximum amplitude reached by the spectrum. We confirm this intuition by fixing the reionisation history and κ but varying α₀ on the range 3.0 <  logα₀ <  4.4 and comparing the resulting spectra: there is a clear linear relation between these two parameters and α₀, but in this case results for rsage and 21CMFAST do not follow the correlation. Interestingly, the shape of the different resulting kSZ power spectra is strictly identical (namely, ℓ^max does not change when varying α₀), hinting at the fact that ℓ^max depends only on κ and not α₀ or reionisation history. Therefore it will be possible to make an unbiased estimate of κ from the shape of the measured spectrum. The rsage simulations show that, for a similar reionisation history, a larger value of α₀ will lead to a stronger kSZ signal; but looking at 21CMFAST, we found that an early reionisation scenario can counterbalance this effect and lead to high amplitude despite low α₀ values. This corroborates the results of Mesinger et al. (2012), which find that the amplitude of the spectrum is determined by both the morphology (and so the α₀ value) and the reionisation history. Therefore, fitting CMB data to our parameterisation will likely lead to strongly correlated values of α₀ and parameters such as Δz or z_re. Other methods should be used to constrain the reionisation history and break this degeneracy, such as constraints from the value of the Thomson optical depth, or astrophysical constraints on the IGM ionised level. Conversely, 21 cm intensity mapping should be able to give independent constraints on α₀.

5.2. Conclusions and 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. With as few as these three parameters, we can fully recover the patchy kSZ angular power spectrum, in a way that is quick and easy to forward-model. Our formalism contrasts with current works, which use an arbitrary patchy kSZ power spectrum template enclosing an outdated model of reionisation. Applying it to CMB data will result in obtaining for the first time the actual shape of the patchy kSZ power spectrum, taking consistently into account reionisation history and morphology. In future works, we will apply this framework to CMB observations from SPT and, later, CMB-S4 experiments. Then, the inferred values of α₀ and κ will provide us with detailed information about the physics of reionisation: κ will constrain the growth of ionised bubbles with time and α₀ the evolution of the variance of the ionisation field during EoR, both being related to the ionising properties of early galaxies. The complex derivation of the kSZ signal, based on a series of integrals, leads to correlations between our parameters. For example, a high amplitude of the spectrum can be explained either by a large value of α₀ due to a high ionising efficiency of galaxies, or by an early reionisation. Such degeneracies, however, could be broken by combining CMB data with other observations: astrophysical observations of early galaxies and quasars will help grasp the global history of reionisation and constrain parameters such as Δz and z_re, while 21 cm intensity mapping will help understand reionisation morphology, putting independent constraints on α₀ and κ. 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. Additionally, this result sheds light on the scaling relations observed in previous works by giving them a physical ground. For example, features in the free electron contrast density power spectrum explain the relation between the amplitude at which the patchy kSZ spectrum bumps ℓ^max and the typical bubble size, which was observed empirically in many previous works (McQuinn et al. 2005; Iliev et al. 2007; Mesinger et al. 2012).

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²) for our fiducial simulations. 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² (0.5 to 1.5 μK², respectively) for plausible reionisation scenarios, therefore lying 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⁻³ <  k/Mpc⁻¹ <  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. 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 contrast to previous works, these results are not simulation-dependent as we have tested the robustness of our model by confronting it to different types of simulations, capturing different aspects of the process. However, the analytic formulation of our derivations was calibrated on a relatively small simulation, of side length 128 Mpc h, which could bias our results. To further support our approach, using a larger radiative hydrodynamical simulation would be useful. Additionally, one could derive the kSZ power from lightcones constructed with our simulation, but the limited size of the simulation might lead to a significant underestimate of the kSZ power (Shaw et al. 2012; Alvarez 2016).

Acknowledgments

The authors thank Anne Hutter and Jacob Seiler for kindly providing runs of their simulations, as well as Jonathan Pritchard and Ian Hothi for fruitful discussions at various stages of this analysis. They also thank the referee for useful comments which helped imoprove the quality of these results. AG acknowledges financial support from the European Research Council under ERC grant number 638743-FIRSTDAWN and her work is supported by a PhD studentship from the UK Science and Technology Facilities Council (STFC). This work was initiated during LSS2LSS, a Ψ² thematic programme organised by the Institut d’Astrophysique Spatiale and funded by the Université Paris-Saclay in July 2018 (see https://www.ias.u-psud.fr/LSS2LSS). The authors were granted access to the HPC resources of CINES and IDRIS under the allocation A0070411049 attributed by GENCI (Grand Equipement National de Calcul Intensif) and the Jean-Zay Grand Challenge (CT4) “Émulation de simulations de Réionisation par apprentissage profond”. This work was additionally supported by the Programme National Cosmology et Galaxies (PNCG) of CNRS/INSU with INP and IN2P3, co-funded by CEA and CNES. SI was supported by the European Structural and Investment Fund and the Czech Ministry of Education, Youth and Sports (Project CoGraDS – CZ.02.1.01/0.0/0.0/15_003/0000437). This research made use of astropy, a community-developed core Python package for astronomy (Astropy Collaboration 2013, 2018); matplotlib, a Python library for publication quality graphics (Hunter 2007); scipy, a Python-based ecosystem of open-source software for mathematics, science, and engineering (Jones et al. 2001) – including numpy (Oliphant 2006), and emcee, an implementation of the affine invariant MCMC ensemble sampler (Foreman-Mackey et al. 2013).

References

Appendix A: Variations on the fit

A.1. Six fits for six simulations

Instead of fitting the six simulations simultaneously, we choose to fit each simulation individually to Eq. (15) with the same error bars as the fitting procedure described in Sect. 4. This allows to use the original P_ee(k, z) data points from each simulation, without interpolating them, and the original reionisation history rather than an averaged one. The results are shown in Table A.1, where the maximum likelihood parameters, along with their 68% confidence intervals, and the corresponding values of 𝒟₃₀₀₀ and ℓ^max are given. The six maximum likelihood values of α₀ and κ lie within the 95% confidence interval of the parameter distributions obtained in Sect. 4 and so do the resulting patchy kSZ spectra, as shown in Fig. A.1.

Fig. A.1. Comparison of the patchy kSZ power spectra resulting from one fit on the six simulations (black solid line, with 68% confidence interval as the shaded area) or from six fits (coloured solid lines). The data point corresponds to constraints from Reichardt et al. (2020).

A.2. Attempt at deriving a covariance matrix from a sample of six

Because of the very insufficient number of simulations available to derive a covariance matrix, even when bootstrapping, we choose to average covariance matrices over bins.

Average over z-bins. First, we choose to ignore correlations between scales over redshifts and use a covariance matrix C, average of the 6 × 10 covariance matrices obtained for each simulation and each redshift bin. C has therefore dimensions (20, 20)¹⁰. We fit Eq. (15) to the six simulations, trying to minimise:

$$\begin{array}{c}\hfill {\chi }^2={\displaystyle \sum _{z_i}{X}_i^{\mathrm{T}}{\mathbf{C}}^{-1}{X}_i,}\end{array}$$(A.1)

where $X_i=P_{\mathit{ee}}^{\mathrm{data}}(\{k_j\},x_i)-P_{\mathit{ee}}^{\mathrm{model}}(\{k_j\},x_i)$. We find a minimal reduced χ² of 125, reached for log α₀ = 4.12 and κ = 0.078 Mpc⁻¹ and giving ${\mathcal{D}}_{3000}^{\mathrm{p}}=0.97\mu {\mathrm{K}}^2$ and ℓ^max = 1500. This difference comes from a poor match between the maximum likelihood P_ee(k, z) and the data points on scales 0.1 <  k/Mpc⁻¹ <  0.3. These scales correspond to the power cut-off, so that the value of κ is poorly constrained and, later, the kSZ power spectrum is distorted.

Average over k-bins. Secondly, we choose to ignore correlations between redshifts over scales and use a covariance matrix C, average of the 6 × 20 covariance matrices obtained for each simulation and each scale bin. C has therefore dimensions (10, 10). Comparing the correlation coefficients obtained for the two approaches, we note that the correlations are higher for this approach. We fit Eq. (15) to the six simulations, trying to minimise:

$$\begin{array}{c}\hfill {\chi }^2={\displaystyle \sum _{k_i}{X}_i^{\mathrm{T}}{\mathbf{C}}^{-1}{X}_i,}\end{array}$$(A.2)

where $X_i=P_{\mathit{ee}}^{\mathrm{data}}(k_i,\{x_j\})-P_{\mathit{ee}}^{\mathrm{model}}(k_i,\{x_j\})$. We find a minimal reduced χ² of 4.46, reached for log α₀ = 3.65 and κ = 0.135 Mpc⁻¹ and giving ${\mathcal{D}}_{3000}^{\mathrm{p}}=1.46\mu {\mathrm{K}}^2$ and ℓ^max = 2700. The excess power comes from the fact that the fit systematically overestimate the P_ee power on small scales (k >  0.3 Mpc⁻¹).

Appendix B: Detailed results on rsage and 21CMFAST

B.1. Fits on 21CMFAST

We now fit Eq. (15) to the power spectra of our three 21CMFAST runs. To account for sample variance, we perform 20 realisations of each simulation – the choice of 20 being motivated by Kaur et al. (2020) and computational limitations. From these 20 realisations we derive relative error bars on P_ee(k, z) values, corresponding to the 68% confidence level on the distribution of values for each bin. The results obtained for 21CMFAST and their interpretation are consistent with what is obtained for the other simulations. The upper panel of Fig. B.1 shows the best-fit model for P_ee(k, z), along with snapshot values, for the second simulation.

Fig. B.1. Result of fitting Eq. (13) to the spectra of the 21CMFAST run for Mturn = 109 M⊙ (upper panel) and of rsage fej (lower panel). The error bars correspond to the 68% confidence level on the spectra of 20 realisations of the same 21CMFAST run.

B.2. Fits on rsage

Because we only have one realisation of each rsage simulation, we apply the relative error bars derived from 21CMFAST to the rsageP_ee(k, z) data points. On the scales and redshifts range covered by the fit, the error bars σ(k, z) derived from the 20 realisations of each of the three 21CMFAST simulations follow σ(k, z) = 10^bP_ee(k, z)(k/k₀)^a, where k₀ = 1 Mpc⁻¹, a = −1.12 ± 0.79 and b = −1.74 ± 0.70 have been found by fitting the σ(k, z) values of the 60 simulations simultaneously. We then apply this expression to the spectra of the rsage simulations, a reasonable first approximation of cosmic variance. We fit Eq. (15) to the spectra of the three simulations. The lower panel of Fig. B.1 shows the best-fit model for P_ee(k, z), along with snapshot values, for rsage fej. Note that here, we only show the spectra on the redshift range used for the fit, where the power-law structure is not as striking as for higher redshifts.

All Tables

All Figures

Fig. 1. Free electrons density contrast power spectrum for a box filled with enough bubbles of radius R = 15 px = 5.5 Mpc to reach a filling fraction f = 1%. Points are results of a numerical computation of the power spectrum, compared to the theoretical model (solid line). The dotted vertical line corresponds to k = 1/R, the dashed vertical line to 91/4/R, the dashed horizontal line to 4/3πR3/f and the tilted dashed line has slope k−4.

In the text

Fig. 2. Left panel: snapshot of the electron density contrast field for the first of the six simulations, at z = 7.2 and xe = 0.49. Right panel: free electrons power spectrum of the same simulation at fixed redshifts (fixed ionised levels). The shaded area corresponds to scales contributing ${\mathcal{D}}_{3000}^{\mathrm{patchy}}$ the most (see Sect. 3.2) and the solid black line to the field shown in the left panel.

In the text

Fig. 3. Result of the fit of Eq. (15) on the free electrons power spectrum of our six simulations, for three redshift bins (left panels) and three scale bins (right panels). The best-fit is shown as the thick black line with the accompanying 68% confidence interval, and the spectra of the six simulations as thin coloured lines. Error bars on data points are computed from the covariance matrix (see text for details).

In the text

Fig. 4. Upper panels: Pee(k, z) as fitted on spectra from the fourth simulation, as a function of scales (left panel) and of redshift (right panel). For reference, the fit is compared to data points for z = 7.8 (xe = 0.26) and k = 0.14 Mpc−1 respectively, with corresponding colour. The width of each line represents the contribution of the redshift (resp. scale) of the corresponding scale (resp. redshift) to the final patchy kSZ amplitude at ℓ = 3000. Lower panels: corresponding probability densities (dashed lines) and cumulative distributions (solid lines). Shaded areas correspond to the first 50% of the signal. The dotted vertical line on the lower right panel marks the midpoint of reionisation.

In the text

Fig. 5. Results for our six simulations. Upper panel: global reionisation histories, for H II and He II. The dotted horizontal line marks the reionisation midpoint zre. Lower panel: angular kSZ power spectrum after fitting Eq. (15) to the Pee(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.

In the text

Fig. 6. Evolution of the amplitude of the patchy power spectrum at ℓ = 3000, 𝒟3000 with the reionisation duration (upper panel) and the reionisation midpoint zre (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.

In the text

Fig. 7. Comparison of results for our six initial simulations, corresponding to a late reionisation scenario, and for an additional seventh simulation, corresponding to an early reionisation scenario. Left panel: reionisation histories. Middle panel: patchy kSZ angular power spectra. The data point corresponds to constraints from Reichardt et al. (2020). Right panel: minimal size of ionised regions as a function of global ionised level. Shaded areas correspond to the 68% confidence level on kSZ amplitude propagated from the probability distributions of the fit parameters.

In the text

Fig. 8. Comparison of results for the three rsage simulations (upper panels) and the three 21CMFAST runs (lower panels) considered. Left panels: reionisation histories. Middle panels: patchy kSZ angular power spectra. The data point corresponds to constraints from Reichardt et al. (2020). Right panels: minimal size of ionised regions as a function of global ionised level. The shaded areas correspond to the 68% confidence interval propagated from the 68% confidence intervals on the fit parameters.

In the text

	Fig. 9. Evolution of the peaking angular scale of the patchy kSZ power spectrum for one given reionisation history but different values of the κ parameter. The red dotted line is the result of a linear regression. Inferences are compared to results for different simulations.
In the text

	Fig. A.1. Comparison of the patchy kSZ power spectra resulting from one fit on the six simulations (black solid line, with 68% confidence interval as the shaded area) or from six fits (coloured solid lines). The data point corresponds to constraints from Reichardt et al. (2020).
In the text

	Fig. B.1. Result of fitting Eq. (13) to the spectra of the 21CMFAST run for Mturn = 109 M⊙ (upper panel) and of rsage fej (lower panel). The error bars correspond to the 68% confidence level on the spectra of 20 realisations of the same 21CMFAST run.
In the text