Constraining cosmological parameters using the pairwise kinematic Sunyaev-Zel’dovich effect with CMB-S4 and future galaxy cluster surveys

The kSZ effect yields a small change in the CMB blackbody temperate, $\Delta T_{\rm kSZ}$, which is proportional to the product of the free electron number density and line-of-sight velocity. The dependence of the kSZ signal on the bulk velocity of the ionized gas in the cluster offers a unique opportunity to probe the cosmological velocity field [47]. The optical depth of a galaxy cluster to CMB photons is low (typically $\tau_{e}\lesssim 0.01$), so we can use the single-scattering limit. In this limit, the change in the CMB temperature due to the kSZ effect is

where $c$ is the speed of light, $\sigma_{T}$ is the Thomson cross section, $n_{e}$ is the number density of electrons, and $\tau_{e}$ is the Thomson optical depth for CMB photons traversing the galaxy cluster [2]. The unaltered CMB temperature is $T_{\rm CMB}$, while $\bf{v}$ is the cluster velocity and $\hat{\mathbf{r}}$ a unit vector in the line of sight.

On scales smaller than the homogeneity scale ($\sim 350$ Mpc [48]), pairs of galaxy clusters are expected to fall towards one another due to their mutual gravitational pull. This infall will produce a small dipole pattern in the CMB temperature anisotropy at the clusters’ positions due to the kSZ effect (e.g., [49]), which is known as the pairwise kSZ signal. The small pairwise kSZ signal can be detected by stacking many such pairs of clusters. The average pairwise kSZ amplitude $T_{\rm pkSZ}(r)$ for such a stack of clusters at a comoving separation $r$ can be related to the mean pairwise velocity $v_{12}(r)$ of the clusters:

where $\bar{\tau}_{e}$ is the average optical depth of the sample. This equation requires two approximations: (1) the internal motion of the cluster’s gas is negligible; and (2) the optical depth and velocity of the clusters are uncorrelated [50]. We adopt a sign convention so that clusters falling towards one another will have a negative relative velocity $v_{12}(r)$ and negative $T_{\rm pkSZ}$ signal.

The mean pairwise velocity of haloes $v_{12}(r)$ separated by comoving distance $r=|\vec{r}_{2}-\vec{r}_{1}|$ can be analytically modeled in linear theory for a specific cosmology and theory of gravity in terms of the two-point matter correlation function $\xi(r)$ as [51, 52, 53]

where $a$ is the scale factor, $H(a)$ is the Hubble parameter, $f(a)\equiv{\rm d}\ln D/{\rm d}\ln a$ is the growth rate (with $D$ being the linear growth factor), $b$ the mass-averaged halo bias, and $\bar{\xi}$ indicates the average of $\xi(r)$ over a comoving sphere of radius $r$. Fig. 1 compares the theoretical approximation of the pairwise velocity with a measurement obtained from the AGORA simulations, described in Sec. III.

Equations (2) and (3) highlight how measurements of the pairwise kSZ are sensitive to a combination of both cluster astrophysics, through the optical depth $\bar{\tau}_{e}$ and halo bias $b$, and cosmology through the Hubble parameter $H(a)$, the growth rate $f$, and the two-point matter correlation function $\xi(r)$. In particular, the dependence on the growth rate $f$ and matter correlation function $\xi(r)$ makes measurements of the pairwise kSZ signal sensitive to $f\sigma_{8}^{2}$, thus providing complementary information to other cosmological probes such as redshift-space distortions, primarily sensitive to $f\sigma_{8}$ [54]. Therefore, if we break the degeneracy with cluster astrophysics through alternative methods, we can use the pairwise kSZ signal to access information of the velocity field (which on large scales is entirely dependent on the growth of the density field) and probe dark energy and modifications of gravity [53, 55, 56, 57, 58].

As noted in the introduction, the kSZ effect is the fainter of the two SZ effects associated with a galaxy cluster. The larger signal is the tSZ effect, which is due to inverse Compton scattering between CMB photons and the hot electrons of the intracluster medium [59, 4]. The magnitude of the tSZ effect for a galaxy cluster can be quantified by the Compton-$y$ parameter [60], and expressed as an integral over the line of sight:

Here $n_{e}$ is the free electron density, and $T_{e}$ is the electron temperature. The physical constants are the speed of light $c$, Boltzmann constant, $k_{B}$ and electron mass $m_{e}$. Importantly, both the tSZ and kSZ effects depend on the integrated electron number density, albeit with different weightings (temperature versus velocity). Thus the measured tSZ signal can potentially be used to calibrate the optical depth dependence of the kSZ signal (see [61, 50, 62]). We return to this in Sec. V.2.

We use realistic realizations of the millimeter wavelength sky from The Agora: Multi-Component Simulation for Cross-Survey Science [63], which leverages the N-body simulations of the MultiDark Planck 2 Synthetic Skies suite [64]. This simulation is based on Planck’s $\Lambda{\rm CDM}$ model [65] with the following cosmological parameters: $H_{0}=67.77\,{\rm km\,s}^{-1}\,{\rm Mpc}^{-1}$, $\Omega_{m}=0.307115$, $\Omega_{b}=0.048206$, $\sigma_{8}=0.8228$, $n_{s}=0.96$.¹¹1The cosmological parameters listed are the Hubble parameter, matter density, baryonic matter density, current root mean square (rms) of the linear matter fluctuations on scales of $8h^{-1}{\rm Mpc}$, and the spectral index of the primordial scalar fluctuations, respectively. The simulated skies are generated by pasting astrophysical effects onto the dark matter haloes, obtained through a friends-of-friends algorithm applied to the MultiDark Planck 2 $N$-body simulation [64]. The astrophysical modeling in the simulation has been calibrated using observational data and external hydrodynamical simulations.

Outlined below are the main components of the simulated microwave sky.

• •

  The dark matter density field is used to gravitationally lens the CMB sky.

• •

  The tSZ signal from each dark matter halo is added based on the gas profile as described in [66]. The electron thermal pressure profile was calibrated on the hydrodynamical BAHAMAS simulations suite [67].

• •

  The kSZ effect is added in a similar way as the tSZ effect. The same [66] gas profile is used to estimate the electron number density, which is multiplied by the line-of-sight velocity to obtain the kSZ signal from the halo.

• •

  The simulations also include correlated signals from dusty star-forming galaxies, also called the cosmic infrared background (CIB).

• •

  Radio galaxies, which appear as point sources in the simulation maps, have been masked and are not included in the measurements.

We convolve the simulated maps by the expected instrumental beams for each CMB-S4 frequency band and add the expected white and 1/f noise levels to each map. Following [46], the assumed map noise has the form

Here $\Delta_{T}$ is the detector noise level, $\ell_{knee}$ and $\alpha_{knee}$ are used to estimate the atmospheric $1/f$ noise, and $\ell$ is the spherical harmonic index. The assumed values of these parameters for each of the three relevant frequency maps are listed in Table 1.

Mock cluster catalogs are a key ingredient to the pairwise kSZ forecasts of this work. We expect the cluster catalogs for the CMB-S4 measurement will come from an optical survey, such as LSST or 4MOST. Existing optical cluster finders, such as redMaPPer [68], identify galaxy clusters as overdensities of red sequence galaxies. The number of red sequence galaxies in clusters (often referred to as richness $\lambda$) is correlated with halo mass and is commonly used as the mass tracer of galaxy clusters.

This approach to cluster detection can suffer from systematic uncertainties, most significantly projection effects due to the chance alignment of structures along the line of sight [69, 70, 71]. Fully simulating these effects is difficult and beyond the scope of this work, as the magnitude is sensitive to galaxy colors, luminosity, and their environmental dependence [70, 72]. However, [73] demonstrates that an approximate treatment, based on counting simulated galaxies within a cylinder along the line of sight, can capture the effect of projection reasonably well. We follow the same cylinder selection method in this work. However, we use a different projection length of $\pm 30$ Mpc/h since we find this better reproduces the observed number counts $N(\lambda)$ as a function of richness in the DES Year 1 data [69].

Within the CMB-S4 footprint, the simulated cluster sample contains $N=654,928$ clusters in the redshift range $z\in[0.2,1.2]$, with an optical richness estimate $\lambda>20$. This corresponds to the redshift and richness range for which the future LSST will provide a complete cluster sample [43]. The photometric redshift estimates from LSST reduce the signal-to-noise ratio by half, as shown in previous works [25]. Therefore, we assume 4MOST will provide the spectroscopic redshifts for the cluster catalog.

Following previous works [17, 20, 25], we implement the pairwise kSZ estimator $\hat{T}_{\rm pkSZ}(r)$ introduced by [74], which has the form

This estimator scales the CMB temperature difference at the location of two clusters by a geometrical factor, $c_{ij}=\hat{\bf r}_{ij}\cdot(\hat{\bf r}_{i}+\hat{\bf r}_{j})/2$, to account for the projection of the pair separation $\hat{\bf r}_{ij}=\hat{\bf r}_{i}-\hat{\bf r}_{j}$ onto the line of sight. The temperature difference will depend on the difference between the kSZ signal in the two clusters, which in turn depends on the line of sight projection of the relative velocity between them.

We reconstruct the pairwise kSZ signal in eight equal and linearly spaced bins between comoving pair separation $r$ of 10 and 300 Mpc, with a bin width of $\Delta r=36.25$ Mpc.

To improve signal-to-noise, we use a constrained internal linear combination (cILC) approach [75] to combine the simulated multi-frequency maps for CMB-S4. Combining the different frequency bands in this way will lead to a CMB+kSZ map, while minimizing noise and the influence from frequency-dependent foregrounds, and simultaneously explicitly removing the tSZ signal. We combine data from different frequency channels using the cILC technique in harmonic space as

where $S$ is the signal we want to retrieve from the frequency maps $M^{i}$, $N_{\nu}$ is the total number of frequency channels, and $\omega^{i}$ is the $\ell$-dependent weights for each frequency channel. These weights are tuned to produce a minimum-variance map along with nulling the contribution of the tSZ signal, where the specific frequency response is known, as

where the matrix $\mathcal{C}_{\ell}$ has dimension $N_{\nu}\times N_{\nu}$ and contains the covariance between simulated maps, $\mathcal{F}$ contains the frequency response vector of the desired signal and the foregrounds, and $N$ is used to extract the desired signal as described in [75].

To extract the CMB+kSZ temperature from the cILC combined map for Sec. V.1, we use an aperture photometry filter on the positions of clusters, which is written in map space as

where $\theta_{r}$ is the characteristic filter scale. The aperture photometry acts as a high pass filter for scales that are larger than the filter scale by subtracting the average temperature in the outer ring from the average temperature inside the disc of radius $\theta_{r}$. In contrast to matched filter techniques used in [24, 25], this approach does not assume a specific model for the cluster profile, but it requires that the cluster is contained within the characteristic filter scale $\theta_{r}$ to avoid biases in the temperature estimation $\hat{T}_{0}(\hat{\mathbf{n}}_{i})$.

In this work, we use the radius $\theta_{500c}$ which is defined as the sphere within which the cluster mass is 500 times the critical density of the Universe at the cluster redshift. This aperture could miss some of the gas in the cluster, thus a model of what fraction is covered within that aperture will be necessary for future analyses. The distribution of $\theta_{500c}$ in our cluster catalog is shown in Fig. 2.

Over an extended redshift range, the redshift evolution of foregrounds (like the CIB) can introduce a redshift-dependent bias in the estimated temperatures at the positions of clusters [25]. To mitigate this, we estimate the mean measured temperature as a function of redshift and subtract this mean temperature from the estimated temperatures from the aperture photometry $\hat{T}_{0}(\hat{\mathbf{n}}_{i})$, as

The mean measured temperature at $z_{i}$ is calculated from the weighted sum of contributions of clusters at redshift $z_{j}$ using a Gaussian kernel $G\left(z_{i},z_{j},\Sigma_{z}\right)=\text{exp}[-(z_{i}-z_{j})^{2}/(2\Sigma^{2}_{z})]$. As shown in [17, 24], the choice of $\Sigma_{z}$ does not impact the result significantly, and we set $\Sigma_{z}=0.02$.

We estimate the covariance matrix of the binned pairwise kSZ measurement directly from the data using jackknife resampling [76]. The jackknife resampling technique consists of measuring the pairwise kSZ signal by splitting the cluster catalogue into $N_{\rm JK}$ subsamples, removing one of them, and recomputing the pairwise kSZ amplitude from the remaining $N_{\rm JK}-1$ subsamples. This process is repeated until every subsample has been discarded once from the measurement. Then we estimate the covariance matrix as

where $\hat{T}^{\alpha}_{i}$ is the pairwise kSZ signal in separation bin $i$ at jackknife realization $\alpha$, of mean $\bar{T}_{i}$.

The baseline covariance matrix in this work is estimated using the jackknife resampling technique with 1000 subsamples. The off-diagonal elements can be significant, with $15-40\%$ correlations between adjacent bins.

We have tested the robustness of the covariance estimate to the choice of 1000 subsamples, and have found that the estimated covariance is stable. We show the changes in the diagonals of the covariance matrix between 1000 to 10,000 subsamples in Fig. 3. The diagonals agree within 1.5%, while the eigenvalues of the covariance matrices agree within 8%. We conclude that 1000 subsamples are sufficient to provide a robust estimate of the covariance matrix.

We present the expected pairwise kSZ measurement for CMB-S4 in this section. Fig. 4 shows the results for a single redshift bin and Fig. 5 decomposes the contribution from multiple redshift bins. These results have been obtained using the cILC map described in IV.2 and for optically selected clusters with richness $\lambda>20$.

We fit the measured pairwise kSZ signal to a one-parameter model. We then compute the statistical significance of our measurement. The reported signal-to-noise ratio (SNR) is obtained by fixing the cosmological parameters and obtaining the best-fit $\bar{\tau}_{e}$ and its uncertainty by minimizing the $\chi^{2}$ as

where $\hat{T}_{\rm pkSZ}$ is the measured signal, $\tilde{C}$ is the covariance matrix estimated using the jackknife method, and ${T}_{\rm pkSZ}(\bar{\tau}_{e})$ is the model computed using Eq.(2). The SNR is then computed with respect to the null-signal as $\text{SNR}=\sqrt{\Delta\chi^{2}}$, with $\Delta\chi^{2}=\chi^{2}(\bar{\tau}_{e})-\chi^{2}(\bar{\tau}_{e}=0)$.

We forecast that the pairwise kSZ signal will be detected at a significance of $36\,\sigma$ across all bins, and a significance of $(12,20,20,10,8)\,\sigma$ for the multiple redshift bins. These signal-to-noise numbers scale with the number of pairs of clusters, which is a function of the number of clusters in each redshift bin. CMB-S4 in combination with future spectroscopic surveys will detect the pairwise kSZ signal at very high significance.

The ability of the pairwise kSZ signal to constrain the velocity field is weakened by the degeneracy between each cluster’s velocity and optical depth. A possible path to break this degeneracy is to use the observed tSZ signal, which also depends on the product of the number density of free electrons and electron temperature, to independently constrain the clusters’ optical depth. Following [61, 50], we assume that the relationship between the mean integrated Compton-y and mean optical depth can be represented by a power law:

where $A$ and $p$ are parameters to be determined by fitting the simulations. We use noiseless y-maps from the Agora simulations and perform a linear least-squares fit of the relationship between the optical depth and Compton-y parameters measured by applying the aperture photometry filter described in Sec. IV.2 to each halo’s lightcone. To allow for redshift evolution, we perform separate fits for each of the five redshift bins used in the pairwise kSZ forecasts in Fig. 5. The best-fit parameters for each redshift bin are presented in Tab. 2, while Fig. 6 shows the points and model fit for the first and last redshift bins.

An important caveat to this approach is that we have calibrated the Compton-y parameter to the integrated number of free electrons, but not the velocity-weighted integral of the number density that shows up in the pairwise kSZ signal. We examine the magnitude of this difference by comparing the mean optical depth of the simulated clusters to the velocity-weighted mean optical depth of the same clusters in each redshift bin in Fig. 7. We find significant offsets between the two quantities, with the mean optical depth differing in each redshift bin by $(6,9,4,11,0)\sigma$, respectively, relative to the velocity-weighted values (the $\sigma$ here is only sample variance and does not include noise). When considering the full sample over all redshifts, the difference is $4\,\sigma$. The discrepancy is likely due to factors like the optically selected catalog containing line-of-sight effects, which the friends-of-friends catalog does not consider. Future work should analyse why this happens and consider approaches to deal with these differences and the magnitude of the biases that would be incurred by using the mean optical depth.

The pairwise kSZ is a direct measure of the infall velocity between pairs of clusters, and thus the gravitational attraction felt by these clusters. We consider two extensions to $\Lambda{\rm CDM}$ that should lead to observable differences in this infall. The first is the $\Lambda{\rm CDM}+\gamma$ model, where the growth index $\gamma$ is defined from the growth rate $f\approx{\Omega_{m}}^{\gamma}$. For general relativity, $\gamma\approx 0.55$, but alternate theories of gravity predict other values of $\gamma$ [77]. The second extension is $w{\rm CDM}$, where the dark energy equation of state $w$ is allowed to vary from the $-1$ value of a cosmological constant.

We use a Fisher information formalism to predict expected constraints on cosmological parameters in each model. The covariance between two parameters $p_{\mu}$ and $p_{\nu}$, from [57], is given by the Fisher matrix, which is calculated as

where $\tilde{C}^{-1}$ is the covariance between two pairwise kSZ bins estimated from the data in Sec. IV.4.

We determine which parameters have a higher effect on the pairwise kSZ signal by computing their log-derivatives with respect to the parameter $p_{\mu}$ being

The log-derivatives give an estimate of how much each parameter contributes to the pairwise kSZ signal, i.e., the higher the log-derivative, the higher the plausible constraints on that parameter. We show these derivatives in Tab. 3, where we can see that the main cosmological parameters that contribute to the signal are $\sigma_{8}$, $\gamma$ and $H_{0}$.

We show the combined constraints of the main parameters in Fig. 8(a) for the $\Lambda{\rm CDM}+\gamma$ and $w{\rm CDM}$ cosmologies. The plotted constraints include the pairwise kSZ measurement in all redshift bins, with a prior on the cluster optical depths. These constraints involve adding different redshift bins, using the mean optical depth $\bar{\tau}_{e}$ values and uncertainties from the $y$-map as a prior, and adding Planck’s 2018 covariance from the TTTEEE low$\ell$+low$E$ $\Lambda{\rm CDM}$ and $w{\rm CDM}$ MCMC chains [65] as a prior to the cosmological parameters. We present the $\Lambda{\rm CDM}+\gamma$ and $w{\rm CDM}$ cosmological parameters in Tab. 4.

In the $\Lambda{\rm CDM}+\gamma$ model, the growth index $\gamma$ covers the entire prior range of $0<\gamma\leq 1.0$. Adding the pairwise kSZ measurement tightly constrains $\gamma$, with an uncertainty of $\sigma(\gamma)=0.02$. The pairwise kSZ data do not appreciably improve constraints on the other $\Lambda{\rm CDM}$ parameters, with the largest fractional improvement being only 4% for $\sigma_{8}$. There are mild degeneracies between $\gamma$ and both $\sigma_{8}$ and $H_{0}$, suggesting that some improvement might be obtained by combining the data with external probes of large-scale structure or the expansion rate. On the other hand, since Planck leaves the posterior of $\gamma$ the same as the prior, we obtain a 96% improvement. For the $w{\rm CDM}$ cosmology, we can see how the pairwise kSZ improves on all the fiducial values from Planck. $H_{0}$, $\sigma_{8}$ and $w_{0}$ are improved by 38%, 46%, and 42%, respectively.

Due to higher redshift uncertainties and/or higher CMB map noise levels, current Stage-3 measurements [20, 25] only detect the pairwise kSZ measurement at $\sim 5\sigma$, compared to $36\sigma$ for CMB-S4. As a result, the Stage-3 measurements are largely unable to constrain cosmological models. For instance, using the Stage-3 uncertainties we estimate a constraint of $\gamma$ at $\sigma(\gamma)\simeq 0.5$, which is 25 times larger than the forecast for CMB-S4. Thus, while Stage-3 measurements can provide initial insights, the advanced capabilities of CMB-S4 are crucial for achieving the precision needed to significantly constrain non-standard cosmological models.

We also observe a correlation between $\bar{\tau}_{e}$ and all the cosmological parameters, which conveys the importance of properly constraining $\bar{\tau}_{e}$ through techniques like the Compton $y$ values of the clusters. We explore the full range of cosmological parameters of the $w{\rm CDM}$ $+\gamma+w_{a}$ that are not as well constrained by the pairwise kSZ in Appendix A.

In this section, we assess the impact of systematic effects for the pairwise kSZ measurement. We consider the following sources.

• •

  Mass scatter: In order to match the mass range from simulations, where the masses of clusters are known, to the optical cluster catalog that is selected in richness, we need a good understanding of how to obtain an accurate representation of the mass range under analysis. This is of particular interest because the analytical model in Eq. 3, which is used to infer the optical depth of clusters, depends on the mass range of interest and changing the typical cluster mass could significantly bias this results. In this work, we have selected the simulation sample using the relation given in [78]. However, we need to take into account that these are estimated, and therefore a scatter in the cluster mass of the optical data can occur. To model this scatter, we draw mass errors from a normal distribution with width $\sigma_{{\rm ln}(M_{500c})}=0.7$, which is an estimation of the scatter for a Gaussian error model in our optically selected catalog. We then use these errors to compute the pairwise kSZ signal from the simulations, obtaining an average decrease of the signal detection of $\sim 1\sigma$, which corresponds to about a $3\%$ reduction on the signal strength.

• •

  Mis-centering: The measured pairwise kSZ signal can be diluted due to the fact that the position of clusters estimated from the optical survey catalog might not coincide with the location of the cluster kSZ signal. This mis-centering has a larger impact on clusters that are not fully relaxed or are merging, where the potential minimum is not located on the brightest cluster galaxy. To model the mis-centering, we draw RA and DEC errors from a normal distribution with width $\sigma_{\rm RA,\;DEC}=0.5\,{\rm arcmin}$. This value corresponds to an uncertainty of $0.2$ Mpc at the mean redshift of the cluster catalog $\bar{z}=0.7$. This uncertainty is $3\sigma$ higher than the mean mis-centering of the optical catalog in comparison with high-resolution X-ray measurements [79]. This error can dilute the signal up to $5\sigma$, which corresponds to a $14\%$ decrease in the significance of the measurement.

To account for these effects, we introduce an amplitude parameter $a_{sel}$ to the pairwise kSZ measurement as

This amplitude will constrain the amount of dilution that occurs to the pairwise kSZ measurement due to mass scatter and mis-centering of the clusters. We estimate this amplitude by calculating the ratio of the pairwise kSZ bins with respect to ones obtained after the selection effects are added to the simulations. Following [80], we estimate the uncertainty as a normal approximation to the distribution of the ratio of two random variables, obtaining $a_{sel}=1.23\pm 0.04$.

If we include $a_{sel}$ in the Fisher formalism, using the previously obtained value as the central value, we obtain the constraints as shown in Fig. 8(a) for $\Lambda{\rm CDM}+\gamma$, and Fig. 8(b) for $w{\rm CDM}$. We observe that $a_{sel}$ is highly correlated with all the parameters, but particularly with $\gamma$ and $w_{0}$. This decreases the improvements coming from the pairwise kSZ, as shown in the numbers inside the parenthesis in Tab. 4. This decrease is worse for $w{\rm CDM}$ cosmology, seeing a reduction of around $15\%$ for all the parameters. Meanwhile, for $\Lambda{\rm CDM}+\gamma$, there is only a 3% reduction on the constraining power of $\gamma$. Thus, the selection effects can decrease the effectiveness of the pairwise kSZ measurement and have to be taken into account.