OPEN ACCESS

# Taking measurements of the kinematic Sunyaev-Zel'dovich effect *forward*: including uncertainties from velocity reconstruction with forward modeling

To cite this article: Nhat-Minh Nguyen et al JCAP12(2020)011

View the [article online](#) for updates and enhancements.

You may also like

- Prospects for kSZ²-Galaxy Cross-correlations during Reionization

Paul La Plante, Jackson Sipple and Adam Lidz

- The 21 cm kSZ-kSZ Bispectrum during the Epoch of Reionization

Paul La Plante, Adam Lidz, James Aguirre et al.

- Detection of Pairwise Kinetic

Sunyaev-Zel'dovich Effect with DESI Galaxy Groups and Planck in Fourier Space

Shaohong Li, Yi Zheng, Ziyang Chen et al.
# Taking measurements of the kinematic Sunyaev-Zel'dovich effect *forward*: including uncertainties from velocity reconstruction with forward modeling

Nhat-Minh Nguyen,ª Jens Jasche,ᵇ Guilhem Lavauxᶜ and Fabian Schmidtª

E-mail: minh@mpa-garching.mpg.de, jens.jasche@fysik.su.se, guilhem.lavaux@iap.fr, fabians@mpa-garching.mpg.de

Received August 19, 2020

Accepted October 1, 2020

Published December 4, 2020

**Abstract.** We measure the kinematic Sunyaev-Zel'dovich (kSZ) effect, imprinted by maxBCG clusters, on the Planck SMICA map of the Cosmic Microwave Background (CMB). Our measurement, for the first time, directly accounts for uncertainties in the velocity reconstruction step through the process of Bayesian forward modeling. We show that this often neglected uncertainty budget typically increases the final uncertainty on the measured kSZ signal amplitude by ≃ 15% at cluster scales. We observe evidence for the kSZ effect, at a significance of ≃ 2σ. Our analysis, when applied to future higher-resolution CMB data, together with minor improvements in map-filtering and signal-modeling methods, should yield both significant and *unbiased* measurements of the kSZ signal, which can then be used to probe and constrain the baryonic content of galaxy clusters and galaxy groups.

**Keywords:** cosmic flows, galaxy clusters, Sunyaev-Zeldovich effect, galaxy clustering

**ArXiv ePrint:** 2007.13721

<figure><img src="image_1.png" /></figure>

<footnote>ªMax-Planck-Institut für Astrophysik, Karl-Schwarzschild-Straße 1, D-85748 Garching, Germany</footnote>

<footnote>ᵇThe Oskar Klein Centre, Department of Physics, Stockholm University, Albanova University Center, SE 106 91 Stockholm, Sweden</footnote>

<footnote>ᶜSorbonne Université, CNRS, UMR 7095, Institut d'Astrophysique de Paris, 98 bis bd Arago, 75014 Paris, France</footnote>
# Contents

<table><tr><td>1 Introduction</td><td>1</td></tr><tr><td>2 Data</td><td>4</td></tr><tr><td>2.1 CMB data</td><td>4</td></tr><tr><td>2.1.1 Planck SMICA CMB map</td><td>4</td></tr><tr><td>2.1.2 Signal extraction</td><td>4</td></tr><tr><td>2.2 Galaxy cluster data</td><td>5</td></tr><tr><td>2.2.1 MaxBCG cluster catalog</td><td>5</td></tr><tr><td>2.2.2 Modeling maxBCG cluster optical depth</td><td>7</td></tr><tr><td>2.3 BORG-SDSS3 reconstructed velocity field</td><td>7</td></tr><tr><td>2.3.1 BORG-SDSS3 reconstruction</td><td>7</td></tr><tr><td>2.3.2 Modeling the large-scale bulk flows of galaxy clusters</td><td>9</td></tr><tr><td>3 Data model and kSZ posterior</td><td>9</td></tr><tr><td>3.1 kSZ likelihood: single angular scale</td><td>10</td></tr><tr><td>3.2 kSZ likelihood: combined signal</td><td>12</td></tr><tr><td>3.3 Modeling photo-z uncertainty</td><td>16</td></tr><tr><td>4 Results</td><td>17</td></tr><tr><td>4.1 Individual-scale signal measurements</td><td>17</td></tr><tr><td>4.2 Combined signal measurement</td><td>18</td></tr><tr><td>5 Null tests for systematics</td><td>20</td></tr><tr><td>5.1 Null tests: individual-scale signal</td><td>21</td></tr><tr><td>5.2 Null tests: combined signal</td><td>21</td></tr><tr><td>6 Discussion and conclusion</td><td>22</td></tr><tr><td>A tSZ contamination and cluster mass cut</td><td>24</td></tr><tr><td>B The kSZ likelihood: mixture weights, MAP mean and variance of the kSZ signal amplitude</td><td>25</td></tr><tr><td>C CMB contribution to covariance matrix of combined kSZ measurement</td><td>27</td></tr><tr><td>D GADGET-2 simulation of the BORG-SDSS3 volume and mock kSZ signal templates</td><td>28</td></tr></table>

## 1 Introduction

Cosmic Microwave Background (CMB) photons, originated from the last-scattering surface, might encounter and interact with moving clouds of free electrons before reaching CMB telescopes or satellites. This phenomenon leaves imprints of cosmic large-scale structures, specifically ionized electron gas inside clusters of galaxies, on the blackbody temperature
anisotropies of the observed CMB, and is often referred to collectively as the Sunyaev-Zel'dovich (SZ) effect [1, 2]. At non-relativistic limit, the coherent part of the electron cloud's motion — following its host cluster and the cosmic flow — can be decoupled from the random thermal part. Since velocity of the former is typically small compared to the speed of light, in the electron gas rest frame $h\nu \ll m_e c^2$ and the photon-electron interaction can be well described by Thomson scattering. The Thomson scattering process essentially introduces a Doppler shift of CMB blackbody temperature in the comoving observer frame. This characteristic shift is specifically described as the kinematic SZ effect.

Consider a single point source, located at position $x$ along direction $\hat{n}$ on the sky. Its kSZ signal is given by [3, 4] — who assumed the case of single, elastic scatterings and, as mentioned above, the regime of low-energy photons (Rayleigh-Jeans limit):

$$ \frac{\Delta T_{\text{kSZ}}(\hat{n})}{T_0} = -\sigma_T \int dl \left( \frac{v_e(x) \cdot \hat{n}}{c} \right) n_e(x) \quad (1.1) $$

where $T_0 = 2.725 \times 10^6 \mu K$ is the CMB blackbody temperature, $n_e(x)$ is the free electron number density at position $x$, and $v_e(x)$ is the peculiar velocity of free electrons, while $\sigma_T$ and $c$ denote the Thomson scattering cross-section and the speed of light in vacuum, respectively. The integral $\int dl$ is performed along the line-of-sight (LOS) $\hat{n}$, from the detector to the last-scattering surface. It is generally assumed that the bulk motion of galaxy clusters follows the large-scale motion of dark matter (DM) [5], i.e. $v_e = v_{DM} = v$, and since the correlation length of the latter is much larger than the physical size of a typical galaxy cluster [6], eq. (1.1) could be further simplified to

$$ \frac{\Delta T_{\text{kSZ}}(\hat{n})}{T_0} = -\tau(x, \hat{n}) (v^{\text{LOS}}(x, \hat{n})/c), \quad (1.2) $$

where $v^{\text{LOS}}(x, \hat{n}) \equiv v(x) \cdot \hat{n}$ denotes the velocity along the LOS $\hat{n}$, and we have defined

$$ \tau(x, \hat{n}) = \sigma_T \int dl \, n_e(x) \quad (1.3) $$

to be the LOS projected optical depth.

The right-hand side (r.h.s) of eq. (1.2) indicates that measurements of the kSZ signal constrain the product $(v(x) \cdot \hat{n})\tau$. Assuming an external constraint on $\tau$ (see, for example, [7]), the kSZ signal then directly measures the peculiar velocity field $v(x)$ and hence allows for constraints not only on modified gravity and Dark Energy models [8–10] but also the sum of neutrino masses [11]. Turning eq. (1.2) the other way around, given $v(x)$, say, reconstructed from galaxy survey data, the kSZ signal directly probes the optical depth $\tau$. The kSZ signal, as can be seen from eq. (1.2), does not depend on the gas temperature and scales linearly with the gas density, similar to the scaling with the gas velocity. Because of that, late-time contribution to the kSZ effect from collapsed, virialized objects appear to be the perfect candidate for probing the otherwise elusive Warm-Hot Intergalactic Medium (WHIM) [12, 13]. This diffuse form of free baryonic gas — whose typical temperature is of order $10^5-10^7 K$ — is too cold to show up in X-ray and thermal SZ measurements. There are mounting evidences suggesting that WHIM might host a large fraction of baryons in our Universe that are still missing compared to the number predicted by our standard cosmological model [13–16]. Further, early contributions from the epoch of reionization [17, 18] can be used to add an additional constraint on the optical depth at re-ionization $\tau_{\text{rei}}$, which would then break
the degeneracy between it and the amplitude of the primordial power spectrum $A_s$ in CMB measurements [19]. Both directions will undoubtedly benefit from upcoming high-resolution CMB experiments and large-volume galaxy redshift surveys. Indeed, the advent of CMB-S4 [20] should allow for the novel application of kSZ tomography, i.e. measurements of kSZ signal at different redshifts. The result of which could then be cross-correlated with datasets from DESI [21] or LSST [22] to constrain either primordial, local non-Gaussianity [23] or the evolution of ionized gas [24].

This wealth of scientific return from kSZ measurements has motivated several attempts to detect this effect using various datasets and estimators, despite the fact that the kSZ signal is deeply buried beneath the primary CMB anisotropies. These efforts have, against the odds, resulted in $\sim 2-4\sigma$ evidence of the kSZ effect using the kSZ pairwise estimator [6, 25-28] and the cross-correlation between CMB maps and a reconstructed velocity field [6, 12, 13, 29, 30]. In this paper, we follow the second approach.

Previously published implementations of this approach [e.g. 6, 12, 13, 29, 30] relied on a reconstruction of the peculiar velocity field around clusters, where $v(\mathbf{x})$ is derived — assuming a certain cosmology with Hubble parameter $H$ and cosmic linear growth rate $f = d \ln \delta/d \ln a$ — by solving the inversion of the linearized continuity equation in either real-space $\mathbf{x}$ [6, 12, 13, 30]

$$ \nabla \cdot \mathbf{v}(\mathbf{x}) = -aHf\delta(\mathbf{x}), \qquad (1.4) $$

or redshift-space $\mathbf{s}$ [29],

$$ \nabla \cdot \mathbf{v}(\mathbf{s}) + f \nabla \cdot [(\mathbf{v}(\mathbf{s}) \cdot \hat{\mathbf{n}}) \hat{\mathbf{n}}] = -aHf\delta(\mathbf{s}), \qquad (1.5) $$

where the DM density field $\delta(\mathbf{x})$ is simply obtained from a smoothed galaxy density field $\delta_g(\mathbf{x})$ by assuming a local, linear bias relation of the form $\delta_g(\mathbf{x}) = b_1\delta(\mathbf{x})$. Specifically, refs. [6, 12] and ref. [13] used galaxy and galaxy group catalogs extracted from SDSS-DR7 [31], while ref. [29] and ref. [30] used CMASS and both CMASS, LOWZ galaxy catalogs obtained from SDSS-DR11 [32] and SDSS-DR12 [33]. This simple method, however, yields only one single realization of the velocity field, regardless of uncertainties in the observed galaxy field. The resulting velocity field thus includes systematics that can potentially bias the kSZ detection and measurement, and the quoted uncertainty on the kSZ signal does not contain the propagated error from the uncertainty in the reconstructed velocity.

In this paper, we derive a posterior for the kSZ signal — as filtered from a temperature anisotropy map at locations of massive clusters — accounting for uncertainties in the velocity reconstruction process by marginalizing over an ensemble of realizations of this field, all of which are compatible with the observed distribution of galaxies [see, e.g. 34, 35]. We apply our method on the Planck SMICA2018 CMB map [36] and the maxBCG cluster catalog [37]. Our ensemble of velocity reconstructions is obtained from the Bayesian forward modeling of galaxy clustering using LOWZ and CMASS galaxy samples presented in [38].

For consistency, except for the estimation of the CMB contribution to the covariance matrix in section 3.2, in this paper we assume the same flat $\Lambda$CDM cosmology assumed by the reconstruction in [38], with $\Omega_r = 0$, $\Omega_K = 0$, $\Omega_m = 0.2889$, $\Omega_b = 0.048597$, $\Omega_\Lambda = 0.7111$, $w = -1$, $n_s = 0.9667$, $\sigma_8 = 0.8159$, $H_0 = 67.74 \text{ km s}^{-1}\text{Mpc}^{-1}$. The paper is structured as follows. In section 2, we describe the datasets used in this work. After formulating the physical models and statistical methods to be applied on the data in section 3, we report our measurements of the kSZ effect from maxBCG clusters and their associated uncertainties in section 4. We then assert the robustness and significance of our measurement by means
of different null tests in section 5. Finally, we summarize our results and discuss relevant systematics as well as future improvements of kSZ measurement in section 6.

# 2 Data

## 2.1 CMB data

We recap here the main features of our CMB data, the Planck SMICA temperature (intensity) map, and subsequently, describe our method of extracting a noisy estimate of the kSZ signal induced by a given galaxy cluster from this map.

### 2.1.1 Planck SMICA CMB map

Our choice of CMB data is the SMICA temperature map from the Planck 2018 release<fnref n="1" /> [36] (SMICA2018 hereafter). SMICA (Spectral Matching Independent Component Analysis) [39] linearly combines Planck frequency channels with multipole-dependent weights, including multipoles up to $ℓ = 4000$ [36], into a single CMB intensity map.

Due to the finite resolution and detector noise associated with any CMB instrument, the observed temperature anisotropy $\Delta T_{\text{obs}}$ is a convolution of the true anisotropy $\Delta T$ — including both primary, i.e. primordial CMB, and secondary anisotropies, e.g. kSZ, thermal SZ (tSZ), integrated Sachs-Wolfe, etc. — with the instrumental beam function $B$,<fnref n="2" /> plus the instrumental noise $\Delta T_{\text{instr}}$, i.e.

$$ \Delta T_{\text{obs}}(\boldsymbol{\theta}_i, \boldsymbol{\theta}) = \int d^2\theta' \Delta T(\boldsymbol{\theta}_i, \boldsymbol{\theta}') B(\boldsymbol{\theta}_i, \boldsymbol{\theta}') + \Delta T_{\text{instr}}(\boldsymbol{\theta}_i, \boldsymbol{\theta}), \quad (2.1) $$

where we have assumed the flat-sky approximation and replaced the three-dimensional vectors $\boldsymbol{x}$ and $\hat{\boldsymbol{n}}$ by the two-dimensional vector $\boldsymbol{\theta}_i$ for a specific cluster $i$. The effective beam function of the SMICA2018 intensity map can be well approximated by a spherically symmetric Gaussian function

$$ B(\boldsymbol{\theta}_i, \boldsymbol{\theta}') = \frac{1}{\sqrt{2\pi\theta_{\text{beam}}^2}} \exp \left[ -\frac{|\boldsymbol{\theta}' - \boldsymbol{\theta}_i|^2}{2\theta_{\text{beam}}^2} \right] \quad (2.2) $$

where the 5-arcmin full-width-half-maximum (FWHM) resolution translates into $\theta_{\text{beam}} \approx 5.0 \text{ arcmin}/\sqrt{8 \ln(2)} \approx 2.1 \text{ arcmin}$.

### 2.1.2 Signal extraction

To extract the kSZ signal from a CMB map, an aperture photometry (AP) filter of radius $\theta_f$ is applied at the location of all clusters. The extracted flux $\Delta T^{\theta_f}$ can then be expressed as a convolution of the observed flux $\Delta T_{\text{obs}}$ with a radial weight function $W^{\theta_f}$ associated with that AP filter, i.e.

$$ \Delta T^{\theta_f}(\boldsymbol{\theta}_i) = \int d^2\theta \, W^{\theta_f}(\boldsymbol{\theta} - \boldsymbol{\theta}_i) \Delta T_{\text{obs}}(\boldsymbol{\theta}_i, \boldsymbol{\theta}) \quad (2.3) $$

Specifically, the spherically symmetric weight function $W^{\theta_f}$ is given by

$$ W^{\theta_f}(\boldsymbol{\theta} - \boldsymbol{\theta}_i) = \begin{cases} 1 & 0 \le |\boldsymbol{\theta} - \boldsymbol{\theta}_i| < \theta_f \\ -1 & \theta_f \le |\boldsymbol{\theta} - \boldsymbol{\theta}_i| < \sqrt{2}\theta_f \\ 0 & \text{otherwise.} \end{cases} \quad (2.4) $$

<footnote n="1">¹https://pla.esac.esa.int/</footnote>

<footnote n="2">²Strictly speaking, there is also a convolution with the pixel window function of the HEALPix map, which we include in our analysis but omit in the equations here and below, for readability.</footnote>
This compensated filter is designed to reduce contributions from primary CMB anisotropies, as well as other low-redshift sources of contamination,<fnref n="3" /> which vary on scales larger than $\theta_f$ in the extracted flux. Thus, as $\theta_f$ increases, so does contamination from these sources in the filtered flux. Consequently, the signal is underestimated at small filter sizes, where parts of the signal fall outside the inner disks and are thus subtracted out. Once the whole cluster is encompassed by the filter, the signal should asymptote to a limiting value<fnref n="4" /> while the uncertainty should increase due to increasing contamination.

In our analysis, we measure the kSZ signal while varying the filter size $\theta_f$. Specifically for the individual-scale measurements of the kSZ signal in section 3.1, section 4.1, and section 5.1, we adopt an adaptive aperture photometry (AAP) filter whose radius $\theta_{f,i}$ scales with the effective apparent size of cluster $i$. This ensures that the filter always probes the same fraction of baryonic gas for each cluster — assuming a universal gas profile. In practice, the application of an AP or AAP filter amounts to taking the difference between the pixel-averaged temperature anisotropies within the inner disk and that within the outer ring. For this estimate of the kSZ signal, the primary noise source for large filter sizes is still the primary CMB, while for small filter sizes — where only a very limited number of pixels are encompassed by the inner disk or outer ring — the instrumental noise dominates.

Our method of extracting the kSZ flux in this paper is similar to that in [6, 12, 28, 29, 40]. It is worth mentioning here that the typical apparent size of maxBCG clusters selected for our analysis is very close to the Planck beam (see section 2.1 and section 2.2.1). We defer a more optimal filtering method (which would require more specific assumptions on the form of the gas profile), such as the matched filter [see, e.g. 13, 26, 41], to applications on CMB data with higher resolutions.

## 2.2 Galaxy cluster data

In this section, we first review relevant properties of our galaxy cluster data, the maxBCG cluster catalog. We then detail how we model the cluster optical depth.

### 2.2.1 MaxBCG cluster catalog

The public version of maxBCG catalog, a volume-limited, red-sequence galaxy cluster sample, includes clusters identified from the SDSS data. These clusters span a scaled richness $N_{200} = 10 - 188$ and a redshift range of $z = 0.1 - 0.3$ [37]. We use the mean richness-mass relation given by eq. (A15) in [42] to convert the scaled richness $N_{200}$ of maxBCG clusters into $M_{200}$, the cluster total mass within the $R_{200}$ radius, accounting for the difference in mass definitions and cosmology [see 43, appendix and references therein for details]. If we assume that the projected gas distribution in cluster $i$, with physical size $R_{200,i}$ and angular diameter distance $D_{A,i}$, can be approximated by a Gaussian profile (cf. eq. (2.6)), then the width of its profile — in the flat-sky approximation — is given by

$$ \theta_{\text{eff},i} = \sqrt{\theta_{200,i}^2 + \theta_{\text{beam}}^2} \qquad (2.5) $$

where $\theta_{200,i} = R_{200,i}/D_{A,i}$. As mentioned earlier in section 2.1.2, for the case of AAP filter, the radius of each filter is adapted to each specific cluster as $\theta_f = \phi_f \theta_{200,i}$. It is worth

<footnote n="3">³Contributions from structures below the redshift range of maxBCG clusters and the BORG-SDSS3 volume, as well as CMB foregrounds, generally manifest themselves as large-scale anisotropies in the observed CMB.</footnote>

<footnote n="4">⁴The exact value of this asymptotic limit depends on various factors, it however should be proportional to the free baryonic fraction within the clusters [29].</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** HEALPix projected map of 3512 clusters selected for our analysis, with 908 clusters with BCG spectroscopic redshift in red circles and 2604 clusters with photometric redshift in blue circles. The size of the circles are scaled with $log_{10}(M_{200})$.</figcaption></figure>

pointing out that an overall shift in the $N_{200} - M_{200}$ relation would affect only the signal amplitude but neither the significance nor the signal-to-noise (S/N) of the kSZ detection. The catalog mean cluster mass and redshift are $M_{200} = 1.288 \times 10^{14} M_{\odot}$ and $z = 0.23$ respectively. Below, we describe various selection cuts that we apply on the original maxBCG catalog, and the resulting sub-sample of maxBCG clusters used in our analysis.

Firstly, to avoid any possible tSZ contamination, we exclude clusters whose $M_{200} > 0.85 \times 10^{14} M_{\odot}$ from our analysis (see appendix A). Next, we select only clusters within regions where the BORG-SDSS3 reconstruction are well-constrained by data, i.e. sky regions where LOWZ and CMASS galaxies are actually observed. In addition, we remove clusters outside of the Planck 2018 common confidence mask recommended for temperature analysis [36].

This leaves us with a final sub-sample consisting of 3512 clusters from the original maxBCG catalog. We show in figure 1 a HEALPix map<fnref n="5" /> [44, 45] of these clusters in Galactic coordinates. We further divide our cluster sample into two datasets:

spec-z: this set includes 908 clusters whose brightest member galaxies (BCG) have spectroscopically determined redshifts, denoted by $z_{spec}$;

photo-z: this set includes 2604 clusters for which only photometric redshifts as inferred by the maxBCG algorithm, denoted by $z_{photo}$, are available.

The mean and median mass of our cluster sample (including both datasets) are $M_{200}^{mean} = 7.18 \times 10^{13} M_{\odot}$ and $M_{200}^{median} = 7.31 \times 10^{13} M_{\odot}$, while the mean and median redshift are $z^{mean} = z^{median} = 0.257$. This results in a mean and a median apparent size of $\theta_{eff}^{mean} = \theta_{eff}^{median} = 3.9'$. The full histograms of redshift and apparent angular size for both sets of clusters are shown in the left and right panels of figure 2, respectively. The apparent size distributions shown in right panel of figure 2 are rather concentrated; their standard deviations are ~ 0.25 arcmin. This actually implies that the AAP filter would not yield significant improvement, as compared to the traditional AP filter, for both sets of cluster considered here. We nevertheless adopt the former for the individual-scale measurements, as already noted in section 2.1.2, since it offers a more physical interpretation of the variation

<footnote n="5">⁵https://healpix.sourceforge.io</footnote>
<figure><img src="image_3.png" /><figcaption>Figure 2. Redshift (left) and apparent size (right) distributions of 3512 maxBCG clusters selected for our analysis. The vertical lines present the mean of redshift distributions for spec-z and photo-z set, respectively.</figcaption></figure>

of the signal with scale. For the combined measurements in section 3.2, section 4.2, and section 5.2, we adopt the latter as it simplifies the amount of computation involved.

### 2.2.2 Modeling maxBCG cluster optical depth

Eq. (1.3) is written for a point source. If a cluster is resolved, we need to integrate eq. (1.3) over the angular extent $\theta_{\text{eff},i}$ of the cluster. For example, assuming a spherically symmetric Gaussian profile for the electron gas, eq. (1.3) becomes

$$ \tau(\theta_i, \theta) = \frac{\tau_{0,i}}{\sqrt{2\pi\theta_{\text{eff},i}^2}} \int_0^\theta d^2\theta \exp\left(-\frac{(\theta - \theta_i)^2}{2\theta_{\text{eff},i}^2}\right), \quad (2.6) $$

where $\tau_{0,i}$ is the integrated optical depth specific for cluster $i$. We will return to this point below in the discussion following eq. (3.1).

As the fraction of electrons held by the neutral gas is presumably small, we will assume that all baryons in all the clusters being considered here are fully ionized, i.e. setting $f_{\text{free}} = 1$ in $N_e = f_{\text{free}}f_{\text{gas}}M_{200,i}/\mu_e m_p$. We further adopt a universal gas-mass fraction $f_{\text{gas}} = f_b \equiv \Omega_b/\Omega_m = 0.16$ following the cosmological baryon abundance and a mean particle weight per electron $\mu_e = 1.17$. Our expression for the integrated cluster optical depth defined in eq. (2.6) then becomes

$$ \tau_{0,i} = \frac{\sigma_T}{D_{A,i}^2} \frac{f_b M_{200,i}}{\mu_e m_p}. \quad (2.7) $$

## 2.3 BORG-SDSS3 reconstructed velocity field

Below, we briefly summarize how the velocity field employed in our analysis is obtained through the process of Bayesian forward modeling, highlighting key features of the BORG algorithm and its output. We then detail how we model the large-scale bulk flow of maxBCG clusters using the reconstructed velocity field based on BORG-SDSS3 outputs.

### 2.3.1 BORG-SDSS3 reconstruction

We employ the non-linear velocity field set by initial conditions within the SDSS3-BOSS survey volume. The latter is constrained in [38] using the Bayesian forward-modeling algorithm BORG [34]. Taking the LOWZ and CMASS galaxies from the SDSS3-BOSS DR13
release [46, 47] as input, BORG systematically explores a 256³-dimensional parameter space of the three-dimensional initial conditions — augmented by the galaxy bias parameters. The algorithm does this essentially by following these three steps repeatedly:

1. Forward-evolving each realization of a *Gaussian* initial matter field from $z_{\text{ini}} \sim 1000$ to $z=0$. The amplitude of the initial fluctuations is set by the fixed, pre-chosen set of cosmological parameters. In [38], the gravitational forward model used to evolve the initial fluctuations is the first-order Lagrangian perturbation theory [48, 49].

2. Apply a galaxy bias transformation, together with other observational effects, including redshift-space distortion, light-cone effect, survey selection function, foreground contamination, etc. (see section 2.2-2.6 of [38] for further details), on the *evolved* matter density field to achieve a *predicted* galaxy field.

3. Compare the *predicted* and *observed* galaxy field (here LOWZ and CMASS) to simultaneously constrain the initial and evolved matter density fields.

The result is a fully probabilistic inference of the cosmic matter and velocity fields, taking into account known and unknown systematic effects (within some pre-defined limits [38]). The setup of the inference is given as follows. The initial conditions are generated on a comoving grid consisting of 256³ cells and covering a comoving volume of 4000³ h⁻³ Mpc³, which corresponds to a grid resolution of $L_{\text{grid}} = 15.624 h^{-1} \text{ Mpc}$. The SDSS3-BOSS data are projected in a sub-volume with the observer located at $\mathbf{x} = \{200, 0, -1700\} h^{-1} \text{ Mpc}$ with respect to the center. A total number of 10360 MCMC samples was collected [38]. Initial power-spectrum analysis in [38] showed that the MCMC chain converged after ~2000 samples, i.e. the chain has passed its burn-in phase and reached a stationary distribution after the first ~2000 samples. Here, we consistently remove all samples whose identifier is less than $s = 2000$. To facilitate the storage and processing of these samples, the chain is thinned by a factor of 10, i.e. we only include 1-in-every-10 samples in our analysis; more details can be found in [38].

Taking initial conditions constrained by the BORG-SDSS3 reconstruction [38] as input, we run DM-only simulations with the same cosmological parameters adopted by the reconstruction using the *temporal* comoving Lagrangian Acceleration algorithm [50] (tCOLA hereafter), and the cloud-in-cell (CIC) projection of particles, to obtain the large-scale velocity field at the maxBCG catalog mean redshift $z = 0.23$. This includes in total 837 tCOLA simulations, one for each of our constrained initial conditions, at the resolution of $N_{\text{part}} = 1024^3$. These are used for the estimation of cluster LOS velocities as well as their uncertainties.

We additionally generate a GADGET-2 [51] simulation at resolution of $N_{\text{part}} = 2048^3$ from initial conditions specified by the sample $s = 9000$ of BORG-SDSS3 reconstruction (see appendix D for details). We use this full N-body, high resolution simulation to specifically:

1. estimate the small-scale motion of clusters unresolved by the BORG-SDSS3 reconstruction and tCOLA re-simulation (see section 2.3.2),

2. verify that our kSZ estimators are unbiased (see section 3.1),

3. measure the cluster signal profile (see section 3.2).
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Histogram of $(v_{\text{NBODY},i}^{\text{LOS}} - v_{\text{tCOLA},i}^{\text{LOS}})$ measured from DM halos within LOWZ-CMASS volume in our GADGET-2 and tCOLA simulations of the BORG-SDSS3 sample $s = 9000$. We apply the same redshift and mass cuts that are applied for maxBCG clusters. The vertical line represents the sample mean. The sample variance is measured at $\sigma^2 = 4.67 \times 10^4 \text{ km}^2\text{s}^{-2}$.</figcaption></figure>

### 2.3.2 Modeling the large-scale bulk flows of galaxy clusters

We model the large-scale bulk flow of galaxy clusters in eq. (1.2) as a sum of two components:

$$v^{\text{LOS}}(\mathbf{x}, \hat{\mathbf{n}}) = v_{L}^{\text{LOS}}(\mathbf{x}, \hat{\mathbf{n}}) + \epsilon_{S}^{\text{LOS}}(\mathbf{x}, \hat{\mathbf{n}}), \qquad (2.8)$$

where $v_{L}^{\text{LOS}}$ is the large-scale LOS bulk-flow estimated from the BORG-SDSS3 reconstruction posterior while $\epsilon_{S}^{\text{LOS}}$ is the unresolved small-scale LOS velocity. We further assume that, for all clusters:

$$\epsilon_S^{\text{LOS}} \sim \mathcal{N}(0, \sigma_{\epsilon_S^{\text{LOS}}}), \qquad (2.9)$$

where $\sigma_{\epsilon_S^{\text{LOS}}} \simeq 206.12 \text{ kms}^{-1}$ is estimated from the standard deviation of $(v_{\text{NBODY},i}^{\text{LOS}} - v_{\text{tCOLA},i}^{\text{LOS}})$ distribution (see figure 3) in which $v_{\text{NBODY},i}^{\text{LOS}}$ and $v_{\text{tCOLA},i}^{\text{LOS}}$ refer to the LOS velocity of halo $i$ as respectively measured from the previously mentioned GADGET-2 simulation and from tCOLA simulation of the same BORG-SDSS3 sample, $s = 9000$.

## 3 Data model and kSZ posterior

Given an inferred ensemble of the large-scale LOS velocity of each galaxy cluster $i$, $v_{L,i}^{\text{LOS}}$, our goal is to construct a likelihood for the extracted kSZ signals at locations of all galaxy clusters in our sample. Below, we will derive this likelihood in two cases of input data:

1. The single-cluster signal is extracted at individual physical scales. This yields multiple measurements of the signal, each using information from a specific scale.

2. The single-cluster signal is simultaneously extracted at multiple scales. This yields one single measurement of the signal combining information from all scales.
While the former can be applied for an analysis focusing on the study of cluster gas profile, we expect the latter to be a more sensitive measurement, as information from all scales is combined and the impact of CMB noise on large scales can be reduced. Our derivations assume that the kSZ measurements at individual cluster locations are independent, i.e. there is no significant overlap between the AP/AAP filters. We verified this is indeed the case for our cluster sample selected from the maxBCG catalog (see section 2.2.1).

## 3.1 kSZ likelihood: single angular scale

Let us express our data model as

$$ \Delta T_{\text{kSZ},i}/T_0 = -\alpha^{\theta_f} \tau_i v_{L,i}^{\text{LOS}}/c - \tau_i \epsilon_{S,i}^{\text{LOS}}/c + \epsilon_{0,i}^{\theta_f}, \quad (3.1) $$

where we have introduced $\alpha^{\theta_f}$ as the amplitude of the kSZ signal from the large-scale bulk flow of cluster i. Here and below, we simply adopt $\tau_i = \tau_{0,i}$. For cases where the filter size is smaller than that of the cluster, this means letting $\alpha^{\theta_f}$ absorb all specific details about the cluster gas profile.

We expect to measure a value of $\alpha^{\theta_f}$ consistent with zero in the case of no detection, whereas a value of order of unity at filter sizes that are large enough to encompass the whole cluster, i.e. $\theta_f \ge \sqrt{\theta_{\text{vir}}^2 + \theta_{\text{beam}}^2}$, corresponds to the expectation from the simple model of optical depth discussed in section 2.2.2. The exact value of $\alpha$ depends on that of $f_{\text{free}}$ (cf. eq. (2.7)) and on the amount of ionized gas outside but associated with the cluster. Further, it is also sensitive to systematics in the amplitude of the reconstructed velocity field and in the weak lensing mass calibration of galaxy clusters [42]. Hence a direct interpretation of $\alpha$ as $f_{\text{gas}}$ would require careful modeling of these uncertainties. In this paper, we focus on the significance of the kSZ signal detection, which is fortunately not significantly affected by any of those.

In eq. (3.1), the first noise term on the r.h.s. is induced by unresolved small-scale LOS velocity $\epsilon_{S,i}^{\text{LOS},s}$ (see section 2.3.2). Note that this term scales with cluster optical depth $\tau_i$. It is thus typically negligible for the clusters considered in our analysis. The second noise term $\epsilon_{0,i}^{\theta_f}$ denotes the residual of primary CMB anisotropies plus inhomogeneous instrumental noise, which we assume to be a Gaussian random noise with zero mean and variance $(\sigma_{0,i}^{\theta_f})^2$.

Eq. (3.1) holds so long as the tSZ plus other foreground contaminations cancel out. A significant degree of cancelation is expected since they are, to first order, uncorrelated with the LOS large-scale velocity $v_{L,i}^{\text{LOS}}$. Given the cluster mass cut introduced in section 2.2.1, we confirmed that this condition holds (see appendix A).

Both the signal amplitude and noise are functions of the AP filter size $\theta_f$ (or AAP filter scale $\phi_f$) — as indicated by the superscript $\theta_f$. However, for the sake of readability, we will omit the superscript $\theta_f$ in all following equations.

The corresponding likelihood distribution for a single velocity realization — as inferred by a BORG-SDSS3 sample s — is given by

$$ \mathcal{P}\left(\left\{\Delta T_{\text{kSZ},i}/T_0\right\} \middle| \alpha, \left\{\tau_i v_{L,i}^{\text{LOS},s}/c\right\}\right) = \prod_i \frac{1}{\sqrt{2\pi\sigma_i^2}} \exp \left\{ -\frac{1}{2\sigma_i^2} \left( \frac{\Delta T_{\text{kSZ},i}}{T_0} + \alpha \tau_i \frac{v_{L,i}^{\text{LOS},s}}{c} \right)^2 \right\} \quad (3.2) $$
where

$$ \sigma_i^2 = \sigma_{0,i}^2 + (\tau_i/c)^2 \sigma_{\epsilon_S^{\text{LOS}}}^2. \quad (3.3) $$

We now seek to construct a posterior distribution for $\alpha$, marginalizing over $N$ velocity realizations, i.e.

$$ \begin{aligned} \mathcal{P}(\alpha | \{\Delta T_{\text{kSZ},i}/T_0\}) &= \int d\{x_i\} \mathcal{P}(\alpha, \{x_i\} | \{\Delta T_{\text{kSZ},i}/T_0\}) \\ &\propto \mathcal{P}(\alpha) \int d\{x_i\} \mathcal{P}(\{x_i\}) \mathcal{P}(\{\Delta T_{\text{kSZ},i}/T_0\} | \alpha, \{x_i\}), \end{aligned} \quad (3.4) $$

where we have used $\{x_i\} = \{\tau_i v_{\text{LOS},i}/c\}$ and introduced the prior on $\alpha$ explicitly as $\mathcal{P}(\alpha)$. The BORG-SDSS3 reconstruction provides a sampled approximation

$$ \mathcal{P}(\{x_i\}) \approx \frac{1}{N} \sum_{s=1}^{N} \delta^D(\{x_i\} - \{x_i\}^s), \quad (3.5) $$

where $\delta^D(\{x_i\})$ denotes the Dirac delta distribution. We thus can rewrite eq. (3.4) as

$$ \mathcal{P}(\alpha | \{\Delta T_{\text{kSZ},i}/T_0\}) \propto \mathcal{P}(\alpha) \frac{1}{N} \sum_{s=1}^{N} \prod_{i} \frac{1}{\sqrt{2\pi \sigma_i^2}} \exp \left\{ -\frac{1}{2\sigma_i^2} \left( \frac{\Delta T_{\text{kSZ},i}}{T_0} + \alpha x_i^s \right)^2 \right\}, \quad (3.6) $$

$$ \propto \mathcal{P}(\alpha) \sum_{s}^{N} \lambda_s \frac{\exp\left(-\frac{(\alpha - \mu_s)^2}{2(\sigma_s)^2}\right)}{\sqrt{2\pi (\sigma_s)^2}}, \quad (3.7) $$

which is a mixture of Gaussian distributions. Each component of the mixture consists of an individual velocity realization — indexed by $s$ — associated with a mixture weight $\lambda_s$, which is given by (see appendix B for a detailed derivation)

$$ \lambda_s = \frac{\exp\left[\omega_s + \frac{1}{2}\ln(2\pi (\sigma_s)^2)\right]}{\sum_s^N \exp\left[\omega_s + \frac{1}{2}\ln(2\pi (\sigma_s)^2)\right]}, \quad (3.8) $$

in which

$$ \sigma_s^2 = \left[ \sum_i \left( \frac{x_i^s}{\sigma_i} \right)^2 \right]^{-1}, \quad (3.9) $$

and

$$ \omega_s = \frac{\mu_s^2}{(2\sigma_s^2)}, \quad \mu_s = \frac{\sum_i [(\Delta T_{\text{kSZ},i}/T_0) x_i^s] / \sigma_i^2}{\sum_i (x_i^s / \sigma_i)^2}. \quad (3.10) $$

Note that the weights $\lambda_s$ give preference to better-fitting realizations of the peculiar velocity field, and $\sum_s^N \lambda_s = 1$.

For simplicity, in what follows, we assume a *uniform* prior on $\alpha$ such that $\mathcal{P}(\alpha) = 1$ for all BORG-SDSS3 velocity realizations. The posterior mean estimate of $\alpha$ is then given by (see appendix B for details)

$$ \langle \alpha \rangle_s = \sum_s^N \lambda_s \mu_s = \sum_s^N \lambda_s \frac{\sum_i [(\Delta T_{\text{kSZ},i}/T_0) x_i^s] / \sigma_i^2}{\sum_i (x_i^s / \sigma_i)^2}, \quad (3.11) $$
while its variance is given by

$$ \sigma_{\alpha}^{2} = \sum_{s}^{N} \lambda_{s} \sigma_{s}^{2} + \sum_{s}^{N} \lambda_{s} (\mu_{s} - \langle \alpha \rangle_{s})^{2}. \qquad (3.12) $$

In the Bayesian language, the significance of each kSZ measurement in supporting the positive-kSZ hypothesis (against the no-kSZ hypothesis) would be quantified through the Bayes factor. In our particular case, this factor can be directly computed as the ratio between two cumulative distribution functions (CDF) of positive-kSZ $P(\alpha > 0)$ and no-kSZ hypothesis $P(\alpha \le 0)$:

$$ BF = \frac{P(\alpha > 0 | \{\Delta T_{\text{kSZ},i}/T_0\})}{P(\alpha \le 0 | \{\Delta T_{\text{kSZ},i}/T_0\})}. \qquad (3.13) $$

Previous kSZ measurements in literature, however, usually reported their significance in terms of distance between the measurement and the value of the null hypothesis, measured in multiples of standard deviation $\sigma$ of a normal distribution [see, e.g. 6, 12, 13, 29, 30]. To facilitate a direct comparison, we adopt the same unit by matching the CDF of the normal distribution to that of our posterior distribution, i.e.

$$ P_{\text{normal}}(x \le -\text{sign.}) = P(\alpha \le 0), \qquad (3.14) $$

and sign. would be our significance measured in the familiar unit of $\sigma$.

Note that eq. (3.11) is identical to eq. (9) in [29] if one takes all $\lambda_s = 1$ and neglects the uncertainty in the velocity reconstruction process. As can be seen in eq. (3.12), our estimator properly includes systematic uncertainties in this step. To better illustrate this point, in figure 4, we compare the uncertainty on the inferred amplitude when one only has a single velocity realization $s$, and when one has an ensemble of realizations $\{s\}$. The single realization case consistently *underestimates* the uncertainty by — from small to large scales — $\simeq 10-20\%$. It is worth pointing out also that different realizations yield different estimates of $\alpha_s$ itself. A combination of, say, high $\alpha_s$ and low $\sigma_{\alpha_s}$, can significantly bias the significance of any kSZ detection that does not account for uncertainties in velocity reconstruction. We note that this trend might explain the large variation in detection significance of previous kSZ measurements, as, for example, summarized in table 1 of [52]. In our case, the significance can be biased by as high as a factor of 2 if one cherry-picks only one single best-fitting BORG-SDSS3 velocity realization.

We further test our estimator in eq. (3.11) on mock input data where a kSZ signal template — generated by DM halos in the GADGET-2 simulation of the BORG-SDSS3 sample $s = 9000$ (see appendix D), assuming a Gaussian gas profile — is injected into a SMICA2018-like CMB map (including both primary CMB and instrumental noise). By artificially varying the noise level, we verified that our estimator is indeed *unbiased*. In the limit of vanishing noise, we obtain $\lambda_{s=9000} \to 1$ while the other weights go to zero, correctly singling out the original BORG-SDSS3 sample used to generate the mock signal template.

## 3.2 kSZ likelihood: combined signal

To combine measurements at different filter sizes $\theta_f$, it is necessary to modify eqs. (3.1)–(3.7) to include the cluster gas profile $f(\theta_f)$ as

$$ \tau_i = \tau_{i,0} f^{\theta_f}, \qquad (3.15) $$
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Ratio between the uncertainty on $\langle \alpha \rangle_s$ when having only one single realization of the peculiar velocity field, as in the cases reported so far in the literature, and that when having instead an ensemble of realizations, as in our case. Each color point represents a realization of the reconstructed velocity field, corresponding to a BORG-SDSS3 MCMC sample.</figcaption></figure>

where we have assumed that this profile is universal. In this work, we obtain an estimate this profile by applying our pipeline on a pure-kSZ signal template generated from individual DM particles found in the GADGET-2 high-resolution simulation (see details in appendix D). We show below, in figure 5, the measurement of $f^{\theta_f}$ at locations of corresponding DM halos identified by Rockstar halo finder<fnref n="6" /> [53, 54], an adaptive hierarchical friends-of-friends (FoF) algorithm in six-dimensional phase-space. Note that we use the same velocity field to assign LOS velocity to the DM particles and halos. For simplicity, in case of the combined signal, we restrict ourselves to the AP filter whose $\theta_f$ does not depend on the cluster effective apparent size.

Eq. (3.1) now becomes

$$ \overrightarrow{\Delta T}_{\text{kSZ},i}/T_0 = -\alpha \tau_{i,0} \vec{f} (v_{L,i}^{\text{LOS}}/c) - \tau_{i,0} \vec{f} (\epsilon_{S,i}^{\text{LOS},s}/c) + \vec{\epsilon}_{0,i}, \quad (3.16) $$

where we have again omitted the superscript $\theta_f$ and instead expressed quantities that depend on $\theta_f$ as vectors to stress the fact that all measurements at different filter scales are now combined in eq. (3.16). Note also that $\alpha$ is now simply a scalar instead of a function of $\theta_f$. Eq. (3.2) can then be rewritten as

$$ \mathcal{P} \left( \left\{ \overrightarrow{\Delta T}_{\text{kSZ},i}/T_0 \right\} | \alpha, \tau_{i,0} \vec{f} \left( v_{L,i}^{\text{LOS},s}/c \right) \right) \propto \prod_i |\mathbf{C}_i|^{-1/2} \exp \left\{ -\frac{1}{2} \left[ \overrightarrow{\Delta T}_{\text{kSZ},i}/T_0 + \alpha \tau_{i,0} \vec{f} \left( v_{L,i}^{\text{LOS},s}/c \right) \right]^T (\mathbf{C}_i)^{-1} \times \left[ \overrightarrow{\Delta T}_{\text{kSZ},i}/T_0 + \alpha \tau_{i,0} \vec{f} \left( v_{L,i}^{\text{LOS},s}/c \right) \right] \right\} \quad (3.17) $$

where the covariance matrix $\mathbf{C}_i$ for cluster $i$ is given by

$$ C_i^{\theta_f \theta_f'} = \left\langle \frac{\Delta T_{\text{CMB}}^{\theta_f}(\theta_i)}{T_0} \frac{\Delta T_{\text{CMB}}^{\theta_f'}(\theta_i)}{T_0} \right\rangle + \left\langle \frac{\Delta T_{\text{instr}}^{\theta_f}(\theta_i)}{T_0} \frac{\Delta T_{\text{instr}}^{\theta_f'}(\theta_i)}{T_0} \right\rangle + \tau_{i,0}^2 f^{\theta_f} f^{\theta_f'} \sigma_{\epsilon_S^{\text{LOS}}}^2, \quad (3.18) $$

<footnote n="6">⁶https://bitbucket.org/gfcstanford/rockstar/</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 5.** The profile $f^{\theta_f}$ at locations of DM halos measured from our mock kSZ template generated by member particles in the GADGET-2 simulation detailed in appendix D. As a consistency check, we divide our DM halo sample of $M_{200} \le 1.5 \times 10^{14} M_{\odot}$ into three mass bins as shown in the plot. Below $M_{200} \le 10^{14} M_{\odot}$, the profile is not very sensitive to $M_{200}$. Only the result for $M_{200} \le 0.85 \times 10^{14} M_{\odot}$ is used in our analysis.</figcaption></figure>

in which we have separated the primary CMB anisotropy and the instrumental noise into the first and second term, respectively. Although, for completeness, we do include the last term on the r.h.s. of eq. (3.18) in our covariance estimate, we note that it is negligible compared to the first two terms — as it scales with $\tau^2$ — and the exclusion of this term would affect neither the signal amplitude nor significance. The second noise contribution is highly inhomogeneous due to the scanning strategy of the Planck satellite and requires instrument-specific mocks to estimate [55, 56]. To this end, it is worth mentioning that Planck does provide a limited set of 300 noise and residual systematics simulations for SMICA2018 [36]. For this work, we instead choose to generate 2500 instrumental noise maps in which the noise value of each pixel is drawn from a zero-mean Gaussian whose variance is given by the corresponding temperature intensity variance in the Planck 2018 HFI Sky Map at frequency 143GHz [55, 56] (see details in appendix C). While this estimate is likely conservative since SMICA is a weighted linear combination of multiple frequency channels, we expect it to be robust and stable, especially at small filter sizes where this term dominates. The first term could be estimated analytically using the Planck 2018 best-fit ΛCDM power spectrum [57]. In the flat-sky limit,

$$ \Delta T_{\text{CMB}}^{\theta_f}(\boldsymbol{\theta}_i) = \int \frac{d\ell}{(2\pi)^2} \exp(i\ell \cdot \boldsymbol{\theta}_i) (\pi\theta_f^2) W(\ell\theta_f) \Delta T_{\text{CMB}}^{\text{obs}}(\ell), \quad (3.19) $$

where $W(l\theta_f)$ is the Fourier transform of the AP filter

$$ W(l\theta_f) = 2 \left[ W_{\text{TH}}(l\theta_f) - W_{\text{TH}}(\sqrt{2}l\theta_f) \right], \quad W_{\text{TH}}(l\theta_f) = 2 \frac{J_1(l\theta_f)}{l\theta_f}, \quad (3.20) $$

while $\ell$ is the two-dimensional wavevector perpendicular to the LOS, $\ell = |\ell|$, and

$$ \Delta T_{\text{CMB}}^{\text{obs}}(\ell) = \Delta T_{\text{CMB}}(\ell) B(\ell). \quad (3.21) $$

The CMB covariance matrix in eq. (3.18) is then explicitly given by (see appendix C)

$$ C_{\text{CMB},i}^{\theta_f \theta_f'} = \frac{\pi \theta_f^2 (\theta_f')^2}{2T_0^2} \int_0^\infty d\ell \, \ell W(\ell\theta_f) W(\ell\theta_f') C_\ell^{\text{CMB}}, \quad (3.22) $$
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Top: correlation matrices of the primary CMB (left panel) and instrumental noise (right panel), as described by the first and second term on the r.h.s. of eq. (3.18). Bottom: the total noise correlation matrix, evaluated using the full r.h.s. of eq. (3.18).</figcaption></figure>

where $C_{\ell}^{\text{CMB}}$ denotes the Planck 2018 ΛCDM best fit angular power spectrum. We show in figure 6 the CMB correlation matrix evaluated from eq. (3.22), the instrumental noise correlation matrix estimated, for one single cluster, from 2500 instrumental noise mocks, and the corresponding total correlation matrix of the AP measurements (including both primary CMB and instrumental noise contributions) at different filter sizes estimated by eq. (3.18). Note that the last two vary between cluster locations due to the inhomogeneity of the Planck instrumental noise.

We use the covariance matrix estimated by eq. (3.18) for both the individual-scale and combined measurements of the signal. Specifically, for the individual-scale case which employs the AAP filter, we interpolate $\sigma_{\text{CMB},i}^2$ (cf. eq. (3.3)) from a 1024x1024 $C_{\text{CMB}}^{\theta_f \theta_f'}$ matrix.

The posterior of $\alpha$ given measurements at all filter sizes can be constructed similarly to eq. (3.7)

$$ \mathcal{P}(\alpha | \{\overrightarrow{\Delta T}_{\text{kSZ},i}/T_0\}) \propto \mathcal{P}(\alpha) \frac{1}{N} \sum_{s=1}^{N} \prod_{i} |\vec{C}_i|^{-1/2} \times \exp \left\{ -\frac{1}{2} \left[ \overrightarrow{\Delta T}_{\text{kSZ},i}/T_0 + \alpha \vec{x}_i^s \right]^T \left( \vec{C}_i \right)^{-1} \left[ \overrightarrow{\Delta T}_{\text{kSZ},i}/T_0 + \alpha \vec{x}_i^s \right] \right\} \quad (3.23) $$

where we have similarly used $\vec{x}_i^s \equiv \{\tau_{i,0} \vec{f}(v_{L,i}^{\text{LOS},s}/c)\}$.
<figure><img src="image_8.png" /><figcaption>Figure 7. Histograms of differences between BCG spectroscopic redshift and cluster photometric redshift for 2129 maxBCG clusters below $M_{200} = 0.85 \times 10^{14} M_{\odot}$. We separate the clusters into three bins of richness $N_{200} = 10$ (top, left panel), $N_{200} = 11$ (top, right panel), $N_{200} = 12$ (bottom panel) but detect no dependence of photo-z error on cluster richness. For clarity, the mean $\mu$ and standard deviation $\sigma$ of each histogram are additionally shown.</figcaption></figure>

The expressions of the posterior mean $\langle \alpha \rangle_n$ and its variance $\sigma_\alpha^2$ are similar to those in eq. (3.11) and eq. (3.12) with modifications to $\mu_s$ and $\omega_s$ as described in appendix B.

## 3.3 Modeling photo-z uncertainty

To account for uncertainties in $v_{L,i}^{\text{LOS},s}$ induced by photometric redshift error (see left panel of figure 7), we introduce an additional sampling step. Specifically, we generate a sample of $N_r$ realizations of maxBCG cluster positions in redshift space, in which we keep the redshifts of spec-z clusters fixed while sampling those of photo-z clusters as

$$z_{\text{photo},i}^r = z_{\text{photo},i}^0 + \delta z_i^r, \qquad (3.24)$$

wherein $z_{\text{photo},i}^0$ is the fiducial photometric redshift of cluster i and

$$\delta z_i^r \sim \mathcal{N}(0, \sigma_z). \qquad (3.25)$$

Here, $\mathcal{N}(0, \sigma_z)$ is a zero-mean Gaussian with a standard deviation of $\sigma_z = 0.015$. The latter value follows the scatter between spectroscopic and photometric redshifts of reference clusters $\sigma_z = \langle \sigma[z_{\text{spec}}(N_{200}) - z_{\text{photo}}(N_{200})] \rangle$, averaged over three $N_{200} = [10, 11, 12]$ bins. As
indicated in the right panel of figure 7, this scatter shows no strong dependence on richness for the clusters considered in our analysis. Note that we have assumed that there is no bias in the fiducial $z_{\text{photo}}$, as the public version of maxBCG catalog is reportedly corrected for this effect [37]. Additionally, we verified that our results are not significantly affected if we are to correct for a possible small bias of the order $\mu \simeq -0.005$ in $z_{\text{photo}}$ (see figure 7).

Note that when drawing from the Gaussian distribution in eq. (3.25), we limit the range of $z_{\text{photo},i}^r$ in all redshift realizations to within [0.05, 0.5]. So our Gaussian distributions of photometric redshift error are actually truncated.

We then introduce an additional sum over all $N_r$ redshift realizations in eq. (3.7),

$$ \mathcal{P}(\alpha | \{\Delta T_{\text{kSZ},i} / T_0\}) = \mathcal{P}(\alpha) \frac{1}{N_r} \sum_{r=1}^{N_r} \frac{1}{N_s} \sum_{s=1}^{N_s} \prod_i \frac{1}{\sqrt{2\pi\sigma_i^2}} \exp \left\{ -\frac{1}{2\sigma_i^2} \left( \frac{\Delta T_{\text{kSZ},i}}{T_0} + \alpha \tau_i v_{L,i}^{\text{LOS},sr} / c \right)^2 \right\}. \quad (3.26) $$

Note that, in the case of the AAP filter, for all redshift realizations, we keep the filter size in eq. (3.26) fixed as $\theta_f = \varphi_f \theta_{200,i}^0$ where $\theta_{200,i}^0$ is computed from the fiducial values $z_{\text{photo},i}^0$. Given the range in which the apparent size of selected maxBCG clusters varies (see right panel of figure 2), this approximation does not affect our estimator in any significant way.

Since the two indices $s$ and $r$ in eq. (3.26) are mathematically equivalent — in the sense that they only appear in $v_{L,i}^{\text{LOS},sr}$ — we can rewrite them as

$$ n \equiv \{s,r\} \quad (3.27) $$

and

$$ v_{L,i}^{\text{LOS},n} \equiv v_{L,i}^{\text{LOS},s} (z_{\text{photo},i}^r) \quad (3.28) $$

so that

$$ \mathcal{P}(\alpha | \{\Delta T_{\text{kSZ},i} / T_0\}) = \mathcal{P}(\alpha) \frac{1}{N} \sum_{n=1}^{N} \prod_i \frac{1}{\sqrt{2\pi\sigma_i^2}} \exp \left\{ -\frac{1}{2\sigma_i^2} \left( \frac{\Delta T_{\text{kSZ},i}}{T_0} + \alpha x_i^n \right)^2 \right\}, \quad (3.29) $$

where $N = N_r N_s$ and $x_i^n \equiv \tau_i v_{L,i}^{\text{LOS},n} / c$.

The derivation of the estimator $\langle \alpha \rangle_n$ and its uncertainty $\sigma_\alpha^2$ is then almost identical to that in section 3 and appendix B, only with $s$ replaced by $n$.

# 4 Results

Below we report the results of our measurement of the kSZ effect, imprinted by selected maxBCG clusters on the SMICA2018 CMB map, in terms of the posterior mean $\langle \alpha^{\varphi_f} \rangle_n$ — over velocity and redshift realizations (cf. eq. (3.11)) — and the associated significance (cf. eq. (3.14)).

## 4.1 Individual-scale signal measurements

We show in figure 8 the inferred value of $\langle \alpha^{\varphi_f} \rangle_n$ as a function of the AAP filter scale $\varphi_f$, using spec-z set alone (top panel) or photo-z set and both sets (bottom panel). We emphasize again that the $1\sigma$ uncertainty shown here as shaded band and reported in table 1 includes also uncertainties in the reconstructed large-scale velocity field.
<table><caption>**Table 1** The posterior mean of the kSZ signal amplitude $⟨α⟩_n \pm 1\sigma_α$ measured at individual AAP filter scales, and the detection significance (cf. eq. (3.14)) are shown for, from left to right, spec-z, photo-z and both datasets. The maximum significance in each case is highlighted.</caption><thead><tr><th rowspan="2">φ<sub>f</sub></th><th colspan="2">spec-z</th><th colspan="2">photo-z</th><th colspan="2">spec-z + photo-z</th></tr><tr><th>⟨αφ<sub>f</sub>⟩<sub>n</sub></th><th>sign.</th><th>⟨αφ<sub>f</sub>⟩<sub>n</sub></th><th>sign.</th><th>⟨αφ<sub>f</sub>⟩<sub>n</sub></th><th>sign.</th></tr></thead><tbody><tr><td>0.7</td><td>0.38 ± 0.52</td><td>0.76</td><td>-0.39 ± 0.36</td><td>-1.06</td><td>-0.01 ± 0.29</td><td>-0.09</td></tr><tr><td>0.8</td><td>0.76 ± 0.68</td><td>1.13</td><td>0.14 ± 0.54</td><td>0.27</td><td>0.44 ± 0.37</td><td>1.16</td></tr><tr><td>0.9</td><td>1.45 ± 0.84</td><td><b>1.71</b></td><td>0.75 ± 0.61</td><td><b>1.19</b></td><td>0.99 ± 0.39</td><td><b>2.21</b></td></tr><tr><td>1.0</td><td>0.65 ± 1.04</td><td>0.63</td><td>0.80 ± 0.79</td><td>0.98</td><td>0.82 ± 0.62</td><td>1.27</td></tr><tr><td>1.1</td><td>1.02 ± 1.32</td><td>0.78</td><td>0.29 ± 1.09</td><td>0.30</td><td>0.63 ± 0.74</td><td>0.85</td></tr><tr><td>1.2</td><td>0.80 ± 1.67</td><td>0.48</td><td>0.84 ± 1.41</td><td>0.61</td><td>0.85 ± 0.99</td><td>0.84</td></tr><tr><td>1.3</td><td>0.30 ± 2.05</td><td>0.15</td><td>1.84 ± 1.60</td><td>1.09</td><td>0.98 ± 1.19</td><td>0.81</td></tr><tr><td>1.4</td><td>-0.35 ± 2.60</td><td>-0.14</td><td>1.83 ± 2.22</td><td>0.85</td><td>0.69 ± 1.53</td><td>0.46</td></tr><tr><td>1.5</td><td>-0.67 ± 3.18</td><td>-0.21</td><td>2.58 ± 2.51</td><td>0.99</td><td>0.99 ± 1.92</td><td>0.51</td></tr><tr><td>1.6</td><td>-0.18 ± 3.76</td><td>-0.05</td><td>4.44 ± 3.11</td><td>1.32</td><td>1.80 ± 2.36</td><td>0.76</td></tr><tr><td>1.7</td><td>0.65 ± 4.58</td><td>0.14</td><td>3.84 ± 3.39</td><td>1.08</td><td>2.12 ± 2.52</td><td>0.83</td></tr></tbody></table>

As can be noted from figure 8, we obtain consistent results between the two datasets spec-z and photo-z, which can be combined as shown in red in the bottom panel of figure 8. Note that the slightly larger uncertainty region of $⟨α⟩_n$ measurement using spec-z set — compared to that of $⟨α⟩_n$ measurement using photo-z set — is mostly caused by their limited quantity, as suggested by the ratio between the two uncertainties being approximately constant across all filter scales. Both panels of figure 8 show that most of the information comes from small scales. As the filter scale increases, the AAP estimate picks up more and more contribution from primary CMB anisotropies, as well as other large-scale sources of contamination, and quickly loses its constraining power.

In addition, we provide the significance at each filter scale for each dataset in table 1. All sets show peaks of significance at $\varphi = 0.9$, close to the angular cluster scale as one should expect. Further, for the photo-z set, the photo-z uncertainty shows a bigger impact at small filter sizes, which is also an expected behavior, since on large scales, the primary CMB is still the dominant noise source.

For illustration, we additionally show the mixture weights (cf. eq. (3.8)) in figure 9, since it can provide some insights on how the information in spec-z and photo-z sets are combined. One can see that the distributions of $\lambda_s$ in the three panels on the right column of figure 9, which show $\lambda_s$ for the combined set, are consistent with those of $\lambda_s$ in the leftmost and middle columns, which respectively show $\lambda_s$ for spec-z and photo-z set. This implies that eq. (3.29) consistently combines information from the spec-z and photo-z to simultaneously and correctly pick out the better-fitting BORG-SDSS3 samples and redshift realizations.

## 4.2 Combined signal measurement

We show in figure 10 our measurements of combined signal amplitude $⟨α^{\theta_f}⟩_n$ (cf. eq. (3.11)), as a cumulative function of the AP filter radius $\theta_f$, for spec-z or photo-z set separately and both sets combined.

As can be seen from both cases of individual and combined signal, the uncertainty in photometric redshift, when accounted for by double sampling, reduces the constraining power of the photo-z set, i.e. adding the clusters from the photo-z set does not substantially improve
<figure><img src="image_9.png" /><figcaption>**Figure 8.** The posterior mean of the kSZ signal amplitude $\langle\alpha\rangle_n$ measured at different AAP filter scales $\phi_f = [0.7, 1.6]$ using spec-z (top panel, red) or photo-z (bottom panel, blue) and both sets (bottom panel, red). The shaded regions denote the corresponding $1\sigma$ uncertainties, which, in our case, including both uncertainties in CMB anisotropies and the reconstructed velocity field.</figcaption></figure>

our significance. The most physically important factor is whether the photo-z error can be well-approximated by a *symmetric* distribution, for example, the Gaussian, centered on the fiducial redshift. As clusters are virialized, collapsed objects, they are more likely to form in overdense regions. It is then reasonable to suspect that the photo-z fluctuations are not symmetrically distributed around the inferred value by the maxBCG algorithm, but rather biased towards overdense regions along the respective line of sight. Further detailed investigation is thus required to identify the optimal way to combine information from both spectroscopic and photometric data, subject to current constraints on computational resources. We defer such an investigation to future work.

The significance of the cumulative combined signal for each dataset is summarized in table 2. As in table 1, most of the information is limited to filter sizes below and around the apparent size of maxBCG clusters in our samples. Above that scale, CMB noise severely limits our significance. This suggests that data from current CMB experiments such as ACTpol and SPT-3G with their much higher resolutions, ~1 arcmin [58, 59], can improve our significance significantly. On the other hand, even at the modest resolution of ~2–3 arcmin, we expect a future experiment such as CMB-S4, with much lower instrumental noise,
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Distributions of the mixture weight $\lambda_n$ over 837 BORG-SDSS3 velocity realizations with mcmc identifier 2000-10360 (x-axes) and 400 redshift realizations (y-axes) for the cases of spec-z (left), photo-z (middle) and both datasets combined (right) at $\phi = 0.9$. Note that we do not sample the redshifts of clusters in the spec-z set, which is why the $\lambda_n$ are evenly distributed among redshift realizations (y-axes) in the leftmost column.</figcaption></figure>

to also yield a significantly improved measurement. In both cases, the details of the gas profile would become important and one would almost certainly need to go beyond the Gaussian profile and the AP filter case considered here.

# 5 Null tests for systematics

In this section, we further assert the significance of our measurement by performing two null tests on mock data in which we

1. shuffle the position of the clusters in our analysis, and

2. replace the SMICA2018 by a set of 300 SMICA-like mock maps, also taken from Planck 2018 data release [36].

Since the photo-z set adds very little information, in the tests below, we will only focus on the spec-z set. We have found that the significance computed from eq. (3.14) can be rather well approximated by the simple ratio $S/N=\langle\alpha^{\phi_f}\rangle_n/\sigma_\alpha$. Below, we employ this approximation to ease the numerical computation of the significance in these null tests.
<figure><img src="image_11.png" /><figcaption>**Figure 10.** The posterior mean of the kSZ combined signal amplitude $\langle\alpha\rangle_n$ combining measurements at progressively larger AP filter sizes, $\theta_f = [3.0, 11.0]$ arcmin, using spec-z (top panel, red) or photo-z (bottom panel, blue) and both sets (bottom panel, red). Going from left to right, each data point combines information from all previous points. The shaded regions denote the corresponding $1\sigma$ uncertainties, including both uncertainties in CMB anisotropies and the reconstructed velocity field.</figcaption></figure>

## 5.1 Null tests: individual-scale signal

For the case of individual-scale measurements using the AAP filter, we measure a signal amplitude consistent with zero for the spec-z set with cluster positions being shuffled, as can be seen in the left panel of figure 11. When applying our pipeline on the set of 300 SMICA2018 simulations (including CMB and instrumental noise) provided by Planck [36], we also recover S/N ratios consistent with zero-mean Gaussian distributions. We show in the top panels of figure 12 the histograms of S/N for increasing individual filter scale from $\varphi = 0.9$ to $\varphi = 1.1$. The histogram of $\varphi = 0.9$ (top row, left panel) appears consistent with our reported significances of $\simeq 1.7$ for the spec-z dataset (cf. table 1).

## 5.2 Null tests: combined signal

For the case of cumulative multi-scale measurements using the AP filter, we also measure a cumulative signal amplitude consistent with zero for spec-z set with cluster positions being shuffled, as shown in figure 11. Our null test using the Planck simulations also recovers
<table><caption>**Table 2.** The posterior mean of the kSZ combined signal amplitude $⟨α⟩_n \pm 1\sigma_α$ measured at increasing maximum AP filter size, and the corresponding detection significance (cf. eq. (3.14)) are shown for spec-z, photo-z and both datasets (left to right). The maximum significance in each case is highlighted. Note that, for readability we only list here results for $\theta_f = [3, 7]$ arcmin.</caption><thead><tr><th rowspan="2">θ<sub>f</sub> [arcmin]</th><th colspan="2">spec-z</th><th colspan="2">photo-z</th><th colspan="2">spec-z + photo-z</th></tr><tr><th>⟨αϕ<sub>f</sub>⟩<sub>n</sub></th><th>sign.</th><th>⟨αϕ<sub>f</sub>⟩<sub>n</sub></th><th>sign.</th><th>⟨αϕ<sub>f</sub>⟩<sub>n</sub></th><th>sign.</th></tr></thead><tbody><tr><td>3.0</td><td>1.23 ± 0.81</td><td>1.51</td><td>0.87 ± 1.04</td><td><b>0.81</b></td><td>1.44 ± 0.57</td><td><b>2.09</b></td></tr><tr><td>3.5</td><td>1.12 ± 0.66</td><td>1.67</td><td>-0.02 ± 0.80</td><td>-0.01</td><td>0.72 ± 0.49</td><td>1.39</td></tr><tr><td>4.0</td><td>0.70 ± 0.65</td><td>1.07</td><td>-0.21 ± 0.76</td><td>-0.26</td><td>0.41 ± 0.51</td><td>0.79</td></tr><tr><td>4.5</td><td>0.67 ± 0.66</td><td>1.03</td><td>-0.32 ± 0.78</td><td>-0.40</td><td>0.34 ± 0.51</td><td>0.67</td></tr><tr><td>5.0</td><td>0.70 ± 0.66</td><td>1.07</td><td>-0.53 ± 0.78</td><td>-0.66</td><td>0.32 ± 0.54</td><td>0.60</td></tr><tr><td>5.5</td><td>1.16 ± 0.60</td><td><b>1.88</b></td><td>-0.86 ± 0.68</td><td>-1.15</td><td>0.59 ± 0.51</td><td>1.12</td></tr><tr><td>6.0</td><td>0.96 ± 0.61</td><td>1.56</td><td>-0.66 ± 0.75</td><td>-0.90</td><td>0.54 ± 0.53</td><td>0.98</td></tr><tr><td>6.5</td><td>0.58 ± 0.67</td><td>0.91</td><td>-1.07 ± 0.55</td><td>-1.63</td><td>0.08 ± 0.62</td><td>0.11</td></tr><tr><td>7.0</td><td>0.55 ± 0.56</td><td>0.98</td><td>-0.88 ± 0.55</td><td>-1.41</td><td>0.06 ± 0.52</td><td>0.12</td></tr></tbody></table>

<figure><img src="image_12.png" /><figcaption>**Figure 11.** Left panel: individual-scale signal amplitude measured using spec-z sample but with sky positions of the clusters shuffled, plotted as a function of AAP filter scale. Right panel: combined signal amplitude measured using spec-z set but with sky positions of the clusters shuffled, plotted as a cumulative function of AP filter sizes.</figcaption></figure>

S/N ratios consistent with zero-mean Gaussian distributions. We show in figure 12 the histograms of S/N for cumulative multi-scale from $\theta_f = 3.0 - 3.5$ arcmin to $\theta_f = 3.0 - 5.5$ arcmin. The histogram of $\theta_f = 3.0 - 5.5$ arcmin (bottom row, right panel) again shows that our reported significance of $\simeq 1.9$ for the spec-z set (cf. table 2) is consistent with the 95% confidence interval.

# 6 Discussion and conclusion

kSZ measurements have shown promising potential to become a new probe for both ionized baryons and Dark Energy with the next generation of CMB experiments. It is thus important to reduce systematic uncertainties in cosmology inference from future kSZ measurements. So far, one of the often neglected sources of systematics is that in the reconstructed velocity. Using a systematic-free ensemble of inferred velocity fields within the SDSS3-BOSS volume,
<figure><img src="image_13.png" /><figcaption>**Figure 12.** Top: histograms of the S/N of individual-scale signal for spec-z set, measured at $\phi_f = [0.9, 1.0, 1.1]$, from 300 simulations of SMICA2018 (including CMB and SMICA-like instrumental noise). The vertical dashed line presents the mean in each histogram. For clarity, the mean $\mu$ and standard deviation $\sigma$ of each histogram are also shown in each plot. Bottom: histograms of the S/N of multi-scale signal for spec-z set, measured cumulatively between $\theta_f = 3.0$ and $\theta_f = 5.5$ arcmin, from 300 simulations of SMICA2018 (including CMB and SMICA-like instrumental noise). The vertical dashed line presents the mean in each histogram. For clarity, the mean $\mu$ and standard deviation $\sigma$ of each histogram are also reported. The y-axis of each histogram displays the BORG-SDSS3 velocity realization count.</figcaption></figure>

in section 3, we have developed a robust likelihood for the kSZ effect that, for the first time, accounts for uncertainties in the velocity reconstruction process into the final uncertainty on the measured signal. As such, the significance of our kSZ measurement can be thought of as a rigorously conservative estimate. It is worth pointing out that extending or optimizing our framework for specific studies is relatively straightforward. For example, one can assume specific cluster gas profiles and adopt more sophisticated filtering techniques. Another advantage of our approach is that a prior on $\alpha$ can be easily introduced for cases wherein the cluster gas profile is known from complementary measurements, e.g. X-ray, tSZ, etc.

In addition, we would like to remark on one of the main concerns regarding our approach: the numerical expense of including the uncertainty from velocity reconstruction through forward modeling. Here, the major bottleneck is, however, the initial conditions reconstruction process in [38] of which this work is built on. A particular question is, to which resolution can this method achieve with larger volume surveys? To this end, it is worth pointing out that there are several performance improvements have been developing and implementing within the BORG framework, including the flexibility of having different resolutions for the initial and evolved matter density field. Technically, one can then achieve the same resolution with, for example, a survey of two-to-four times effective volume size, while having to sample only the same parameter space (of the initial conditions and bias parameters). All things considered, the advantage of including the uncertainty in velocity reconstruction (so as to not bias any inferred quantity, for example, the optical depth measured from the kSZ effect) far outweighs the disadvantage of numerical complexity it introduces. This is especially important as recent and upcoming kSZ measurements keep improving their S/N ratios [60, 61].
We apply our method to measure the kSZ signal, imprinted by large-scale bulk flow of selected maxBCG clusters, in the SMICA2018 CMB map. We observe, as reported in section 4, moderate evidence of the signal at $\sim 2\sigma$ — including velocity uncertainty. In this simple demonstration, the sensitivity of our kSZ measurement is limited by several factors, including the simple AP filter approach, the typical cluster scale being close to the resolution of the SMICA2018 map, and the current level of instrumental noise. Yet, we find that ignoring systematic uncertainty in the reconstructed velocity could already lead to a significant bias in the kSZ measurement (cf. eq. (3.12) and figure 4). It would be interesting to study how this would affect kSZ measurements using more sophisticated filtering methods, for example, the matched-filter approach, and data from higher resolution and sensitivity CMB experiments, e.g. CMB-S4 [20], SPT-3G [62].

For the sole purpose of kSZ detection, another important factor is the sheer number of clusters. In fact, our main dataset, the maxBCG spec-z sample, is quite limited in number at only 908 clusters. To overcome this issue, we have explored in section 3.3 the possibility of including clusters with only photometric redshifts, the maxBCG photo-z sample, by sampling each clusters’ redshift from a symmetric Gaussian distribution centered on the fiducial value. However, we only saw a modest increase in detection significance when including photometric clusters. We leave a detailed investigation on whether the true underlying distribution of redshift fluctuations is Gaussian, and whether more information can be extracted from photometric clusters, for future work.

In this work, we have ignored the thermal SZ (tSZ) signal from clusters, and opted to remove the most massive clusters from the analysis. One could also think of simultaneously modeling and measuring both kSZ and tSZ signals for each cluster, so that there is no need for the removal of massive clusters. Nothing, in principle, prevents this simple extension of the data model in eq. (3.1). In fact, we plan to pursue this approach in a follow-up study.

**Data and code availability.** For the BOSS-SDSS3 reconstruction data, please contact [https://www.aquila-consortium.org/contact/](https://www.aquila-consortium.org/contact/). The numerical implementation of the kSZ likelihood presented in this paper is available upon request to the corresponding author. Public releases of both will be available in the future, accompanying follow-up works.

## Acknowledgments

We would like to thank Eiichiro Komatsu, Jochen Weller, Franz Elsner, Alex Barreira, Giovani Cabass, Seunghwan Lim, Emmanuel Schaan, Naonori Sugiyama, Daisuke Nagai, Shun Saito, as well as the members of the Aquila Consortium for inspiring and helpful discussions. MN and FS acknowledge support from the Starting Grant (ERC-2015-STG 678652) “GrInflaGal” of the European Research Council. GL acknowledges financial support from the ILP LABEX (under reference ANR-10-LABX-63) which is financed by French state funds managed by the ANR within the Investissements d’Avenir programme under reference ANR-11-IDEX-0004-02. This work was supported by the ANR BIG4 project, grant ANR-16-CE23-0002 of the French Agence Nationale de la Recherche. This work is done within the Aquila Consortium.<fnref n="7" />

## A tSZ contamination and cluster mass cut

Ref. [29] pointed out that the tSZ contamination becomes important for rare, massive clusters with total mass larger than $10^{14} M_{\odot}$. We have found that, specifically for our sub-sample

<footnote n="7">⁷https://aquila-consortium.org</footnote>
<figure><img src="image_14.png" /><figcaption>Figure 13. The average AAP filter output from SMICA2018 CMB map at locations of maxBCG clusters below $1.1 \times 10^{14} M_{\odot}$, binned in $M_{200}$ (left panel) and cumulative (right panel). The dot-dashed lines represent the typical kSZ signal amplitude of clusters in the $M_{200}$ bin of the same color, assuming a universal $v^{\text{LOS}} = 300 \text{ km s}^{-1}$.</figcaption></figure>

of maxBCG clusters after redshift and survey mask cuts and the SMICA2018 CMB map, it is necessary to remove all clusters with mass greater than $8.5 \times 10^{13} M_{\odot}$. This is demonstrated in figure 13, where we show the average AAP filter output, as a function of filter scale, at locations of clusters divided in two equi-log $M_{200}$ bins, up to $M_{200} = 1.1 \times 10^{14} M_{\odot}$. We especially check for the cancellation of tSZ signal (and other possible foreground contaminations) by comparing the average AAP filter output to the typical amplitude of the kSZ signal of clusters in each corresponding mass bin. For simplification, we assume a LOS velocity $v^{\text{LOS}} = 300 \text{ km s}^{-1}$ for all clusters. Both panels of figure 13 show signs of a significant tSZ contamination when including clusters more massive than $8.5 \times 10^{13} M_{\odot}$.

## B The kSZ likelihood: mixture weights, MAP mean and variance of the kSZ signal amplitude

Here, we first derive the individual mixture weight $\lambda_s$ for each component of the Gaussian mixture distribution of large-scale bulk-flow amplitude $\alpha_s$ in section 3. Let us denote

$$A_i = \Delta T_{\text{kSZ},i} / T_0, \qquad (B.1)$$

and

$$B_i^s = -\tau_i v_{L,i}^{\text{LOS},s} / c. \qquad (B.2)$$

Then, the likelihood on the r.h.s. of eq. (3.7) can be re-written as

$$\ln \mathcal{P}(\{A_i\} | \alpha, \{B_i^s\}) = -\frac{1}{2} \sum_i \left[ \frac{A_i}{\sigma_i} \right]^2 - \frac{1}{2} \left[ \sum_i \left( \frac{B_i^s}{\sigma_i} \right)^2 \right] \left[ \alpha^2 - 2\alpha \frac{\sum_i \frac{A_i B_i^s}{\sigma_i^2}}{\sum_i \left( \frac{B_i^s}{\sigma_i} \right)^2} \right] - \frac{1}{2} \sum_i \ln[\sigma_i^2]. \qquad (B.3)$$

We further denote

$$\mu_s = \frac{\sum_i \left[ (A_i B_i^s) / \sigma_i^2 \right]}{\sum_i (B_i^s / \sigma_i)^2}, \qquad \sigma_s^2 = 1 / \left[ \sum_i \left( \frac{B_i^s}{\sigma_i} \right)^2 \right], \qquad \gamma = \sum_i \left( \frac{A_i}{\sigma_i} \right)^2, \qquad (B.4)$$
and

$$ \delta = \sum_i \ln[\sigma_i^2], \qquad (B.5) $$

so that we can shorten eq. (B.3) to

$$ \ln \mathcal{P} = -\frac{1}{2}\gamma + \frac{\mu_s^2}{2\sigma_s^2} - \frac{1}{2}\delta - \frac{1}{2}\sigma_s^2 (\alpha - \mu_s)^2. \qquad (B.6) $$

The first and third terms on the r.h.s. of eq. (B.6) neither vary between BORG-SDSS3 samples nor depend on $\alpha$, thus they can be dropped from eq. (B.6), such that

$$ \ln \mathcal{P} \propto \frac{1}{2\sigma_s^2} [\mu_s^2 - (\alpha - \mu_s)^2]. \qquad (B.7) $$

We can now simply write

$$ \mathcal{P} \left( \alpha | \left\{ \Delta T_{\text{kSZ},i} / T_0 \right\}, \left\{ \tau_i v_{L,i}^{\text{LOS},s} / c \right\} \right) \propto \mathcal{P}(\alpha) \frac{1}{N} \sum_s^N e^{\omega_s} \sqrt{2\pi \sigma_s^2} \frac{\exp\left(-\frac{(\alpha-\mu_s)^2}{2\sigma_s^2}\right)}{\sqrt{2\pi \sigma_s^2}} \qquad (B.8) $$

where $\omega_s \equiv \frac{\mu_s^2}{2\sigma_s^2}$. The normalized version will then be

$$ \begin{aligned} \mathcal{P} \left( \alpha | \left\{ \Delta T_{\text{kSZ},i} / T_0 \right\}, \left\{ \tau_i v_{L,i}^{\text{LOS},s} / c \right\} \right) &= \mathcal{P}(\alpha) \frac{\frac{1}{N} \sum_s^N e^{\omega_s} \sqrt{2\pi (\sigma_s)^2} \frac{\exp\left(-\frac{(\alpha-\mu_s)^2}{2(\sigma_s)^2}\right)}{\sqrt{2\pi (\sigma_s)^2}}}{\frac{1}{N} \sum_s^N e^{\omega_s} \sqrt{2\pi (\sigma_s)^2}} \\ &= \mathcal{P}(\alpha) \sum_s^N \lambda_s \frac{\exp\left(-\frac{(\alpha-\mu_s)^2}{2(\sigma_s)^2}\right)}{\sqrt{2\pi (\sigma_s)^2}}, \end{aligned} \qquad (B.9) $$

where

$$ \lambda_s = \frac{e^{\omega_s + \frac{1}{2} \ln[2\pi(\sigma_s)^2]}}{\sum_s^N e^{\omega_s + \frac{1}{2} \ln[2\pi(\sigma_s)^2]}}. \qquad (B.10) $$

The mean of this distribution can be derived as follows, assuming a uniform prior, i.e. $\mathcal{P}(\alpha) = 1$,

$$ \begin{aligned} \langle \alpha \rangle_s &= \int d\alpha \, \alpha \sum_s^N \lambda_s \frac{\exp\left(-\frac{(\alpha-\mu_s)^2}{2(\sigma_s)^2}\right)}{\sqrt{2\pi (\sigma_s)^2}} \\ &= \sum_s^N \lambda_s \int d\alpha \, \alpha \frac{\exp\left(-\frac{(\alpha-\mu_s)^2}{2(\sigma_s)^2}\right)}{\sqrt{2\pi (\sigma_s)^2}} = \sum_s^N \lambda_s \mu_s, \end{aligned} \qquad (B.11) $$

which is precisely eq. (3.11).

Similarly, we can explicitly work out the variance of the Gaussian mixture distribution:

$$ \begin{aligned} \sigma_{\alpha}^{2} &= \langle (\alpha - \langle \alpha \rangle_{s})^{2} \rangle \\ &= \int d\alpha \mathcal{P} \left( \alpha | \left\{ \Delta T_{\text{kSZ},i} / T_0 \right\}, \left\{ \tau_i v_{L,i}^{\text{LOS},s} / c \right\} \right) (\alpha - \langle \alpha \rangle_{s})^{2} \end{aligned} $$
$$ 
\begin{aligned}
&= \sum_s^N \lambda_s \int d\alpha \frac{\exp\left(-\frac{(\alpha-\mu_s)^2}{2(\sigma_s)^2}\right)}{\sqrt{2\pi (\sigma_s)^2}} (\alpha - \langle\alpha\rangle_s)^2 \\
&= \sum_s^N \lambda_s (\sigma_s)^2 + \sum_s^N \lambda_s [(\mu_s - \langle\alpha\rangle_s)^2].
\end{aligned} 
\quad (B.12) 
$$

So the result is the sum of the average noise variances and the variance of the mean estimate (cf. eq. (3.12)). The second term clearly shows that our uncertainty on $\alpha$ includes also the uncertainty from the velocity reconstruction.

We next compute $\ln \mathcal{P}(\{\mathbf{A}_i\}|\alpha, \{\mathbf{B}_i^s\})$ for measurements at all $\theta_f$ scales in a similar fashion to how we arrived at eq. (B.9):

$$ 
\begin{aligned}
\ln \mathcal{P}(\{\mathbf{A}_i\}|\alpha, \{\mathbf{B}_i^s\}) &= -\frac{1}{2} \sum_i [\mathbf{A}_i - \alpha \mathbf{B}_i^s]^{\top} (\mathbf{C}_i)^{-1} [\mathbf{A}_i - \alpha \mathbf{B}_i^s] - \frac{1}{2} \sum_i \ln |\mathbf{C}_i| \\
&\propto -\frac{1}{2} \sum_i (\mathbf{B}_i^s)^{\top} (\mathbf{C}_i)^{-1} \mathbf{B}_i^s \left[ \alpha^2 - 2\alpha \frac{\sum_i \mathbf{A}_i^{\top} (\mathbf{C}_i)^{-1} \mathbf{B}_i^s}{\sum_i (\mathbf{B}_i^s)^{\top} (\mathbf{C}_i)^{-1} \mathbf{B}_i^s} \right] 
\end{aligned} 
\quad (B.13) 
$$

where we have omitted terms that do not vary between BORG-SDSS3 samples. Eq. (B.13) is similar to eq. (B.7) with

$$ \mu_s = \frac{\sum_i \mathbf{A}_i^{\top} (\mathbf{C}_i)^{-1} \mathbf{B}_i^s}{\sum_i (\mathbf{B}_i^s)^{\top} (\mathbf{C}_i)^{-1} \mathbf{B}_i^s} \quad (B.14) $$

and

$$ \sigma_s^2 = 1 / \left[ \sum_i (\mathbf{B}_i^s)^{\top} (\mathbf{C}_i)^{-1} \mathbf{B}_i^s \right]. \quad (B.15) $$

## C CMB contribution to covariance matrix of combined kSZ measurement

We provide here a detailed derivation of the CMB covariance matrix term (cf. eq. (3.22)) that contributes to the covariance matrix of the combined kSZ measurement described in section 3.2. Let us plug eq. (3.19) into the first term on the r.h.s. of eq. (3.18), to obtain

$$ 
\begin{aligned}
\mathbf{C}_{\text{CMB},i}^{\theta_f \theta_f'} &= \left( \frac{\pi \theta_f \theta_f'}{T_0} \right)^2 \left\langle \int \frac{d\ell}{(2\pi)^2} \int \frac{d\ell'}{(2\pi)^2} \right. \\
&\qquad \left. \exp \left[ i (\ell - \ell') \cdot \boldsymbol{\theta}_i \right] \hat{W}(\ell \theta_f) \Delta T_{\text{CMB}}^{\text{obs}}(\ell) \hat{W}(\ell' \theta_f') \Delta T_{\text{CMB}}^{\text{obs}}(\ell') \right\rangle \\
&= \frac{\pi \theta_f^2 (\theta_f')^2}{2T_0^2} \int_0^\infty d\ell \, \ell \, \hat{W}(\ell \theta_f) \hat{W}(\ell \theta_f') C_{\ell}^{\text{CMB}}.
\end{aligned} 
\quad (C.1) 
$$

We further validate the analytical computation of the CMB covariance matrix in eq. (C.1) by comparing its diagonal elements with the sample variance of $\Delta T_{\text{CMB}}^{\theta_f}(\boldsymbol{\theta}_i)/T_0$ computed at 1000 random points in each of the 1000 SMICA2018-like CMB mocks in figure 14. The analytical calculation using the Planck 2018 best-fit ΛCDM power spectrum is in extremely good agreement with the numerical estimate using SMICA2018-like CMB realizations.
<figure><img src="image_15.png" /><figcaption>Figure 14. Comparison between analytical (red) and numerical (blue) estimates of CMB contribution to the diagonal of combined kSZ signal covariance matrix.</figcaption></figure>

## D GADGET-2 simulation of the BORG-SDSS3 volume and mock kSZ signal templates

In order to generate the mock template of kSZ signal within the BORG-SDSS3 volume, we use DM particles and halos from a GADGET-2 [51] simulation with DM-only at a very high resolution of $N_{part} = 2048^3$. The initial conditions for the simulation is taken from the BORG-SDSS3 sample $s = 9000$. The halos are identified as main halos by the Rockstar halo finder algorithm<fnref n="8" /> [53, 54] with a minimum number of 20 particles per halo. The cosmology and box size of this simulation agree with those of our BORG-SDSS3 reconstruction; the initial conditions are specifically taken from sample 9000. The high resolution allows us to achieve a correct halo mass function (HMF) down to $M_h = 2 \times 10^{13} h^{-1} M_{\odot}$ at redshift $z = 0.23$, as shown in figure 15.

For the validity test of our estimator, we model the gas profile of each halo $i$ with a Gaussian profile [see, e.g. 28, 29]:

$$n_{e,i}(\theta) = \frac{N_{e,i}}{\sqrt{2\pi\theta_i^2}} \exp\left(-\frac{\theta^2}{2\theta_i^2}\right) \qquad (D.1)$$

where $\theta_i^2 = \theta_{200,i}^2 + \theta_{\text{beam}}^2$ and $N_{e,i} = (f_b M_{200,i})/(\mu_e m_p)$. For the measurement of the signal profile in section 3.2, we directly generate the kSZ template using all individual DM particles within the same volume analyzed in this work. The LOS velocity of each particle (or halo) is interpolated from the tCOLA simulation of BORG-SDSS3 sample 9000.

<footnote n="8">⁸https://bitbucket.org/gfcstanford/rockstar</footnote>
<figure><img src="image_16.png" /><figcaption>**Figure 15.** The number density of DM halos identified by Rockstar in our high-resolution GADGET-2 simulation of BORG-SDSS3 sample 9000, as compared to the Tinker 2008 halo mass function [63], at redshift $z = 0.23$.</figcaption></figure>

References

[1] R.A. Sunyaev and Y.B. Zeldovich, *The observations of relic radiation as a test of the nature of X-ray radiation from the clusters of galaxies*, Comm. Astrophys. Space Phys. **4** (1972) 173.

[2] J.P. Ostriker and E.T. Vishniac, *Generation of microwave background fluctuations from nonlinear perturbations at the ERA of galaxy formation*, Astrophys. J. Lett. **306** (1986) L51 [INSPIRE].

[3] P.R. Phillips, *Calculation of the kinetic Sunyaev-Zeldovich effect from the Boltzmann equation*, Astrophys. J. **455** (1995) 419.

[4] M. Birkinshaw, *The Sunyaev-Zel'dovich effect*, Phys. Rept. **310** (1999) 97 [astro-ph/9808050] [INSPIRE].

[5] V. Desjacques, D. Jeong and F. Schmidt, *Large-scale galaxy bias*, Phys. Rept. **733** (2018) 1 [arXiv:1611.09787] [INSPIRE].

[6] PLANCK collaboration, *Planck 2015 results. XXIV. Cosmology from Sunyaev-Zeldovich cluster counts*, Astron. Astrophys. **594** (2016) A24 [arXiv:1502.01597] [INSPIRE].

[7] M.S. Madhavacheril, N. Battaglia, K.M. Smith and J.L. Sievers, *Cosmology with the kinematic Sunyaev-Zeldovich effect: breaking the optical depth degeneracy with fast radio bursts*, Phys. Rev. D **100** (2019) 103532 [arXiv:1901.02418] [INSPIRE].

[8] R. Keisler and F. Schmidt, *Prospects for measuring the relative velocities of galaxy clusters in photometric surveys using the kinetic Sunyaev-Zel'dovich effect*, Astrophys. J. Lett. **765** (2013) L32 [arXiv:1211.0668] [INSPIRE].

[9] Y.-Z. Ma and G.-B. Zhao, *Dark energy imprints on the kinematic Sunyaev-Zel'dovich signal*, Phys. Lett. B **735** (2014) 402 [arXiv:1309.1163] [INSPIRE].

[10] E.-M. Mueller, F. de Bernardis, R. Bean and M.D. Niemack, *Constraints on gravity and dark energy from the pairwise kinematic Sunyaev-Zel'dovich effect*, Astrophys. J. **808** (2015) 47 [arXiv:1408.6248] [INSPIRE].
[11] E.-M. Mueller, F. de Bernardis, R. Bean and M.D. Niemack, *Constraints on massive neutrinos from the pairwise kinematic Sunyaev-Zel'dovich effect*, Phys. Rev. D 92 (2015) 063501, [arXiv:1412.0592] [INSPIRE].

[12] C. Hernández-Monteagudo et al., *Evidence of the missing baryons from the kinematic Sunyaev-Zel'dovich effect in Planck Data*, Phys. Rev. Lett. 115 (2015) 191301, [arXiv:1504.04011] [INSPIRE].

[13] S. Lim, H. Mo, H. Wang and X. Yang, *The detection of missing baryons in galaxy halos with kinetic Sunyaev-Zel'dovich effect*, arXiv:1712.08619.

[14] M. Fukugita, C.J. Hogan and P.J.E. Peebles, *The cosmic baryon budget*, Astrophys. J. 503 (1998) 518 [astro-ph/9712020] [INSPIRE].

[15] J. Shull, B.D. Smith and C.W. Danforth, *The baryon census in a multiphase intergalactic medium: 30% of the baryons may still be missing*, Astrophys. J. 759 (2012) 23 [arXiv:1112.2706] [INSPIRE].

[16] F. Nicastro et al., *Observations of the missing baryons in the warm-hot intergalactic medium*, Nature 558 (2018) 406 [arXiv:1806.08395] [INSPIRE].

[17] W. Hu, *Reionization revisited: secondary CMB anisotropies and polarization*, Astrophys. J. 529 (2000) 12 [astro-ph/9907103] [INSPIRE].

[18] C. Dvorkin and K.M. Smith, *Reconstructing patchy reionization from the cosmic microwave background*, Phys. Rev. D 79 (2009) 043003 [arXiv:0812.1566] [INSPIRE].

[19] M.A. Alvarez, S. Ferraro, J.C. Hill, R. Hložek and M. Ikape, *Mitigating the optical depth degeneracy using the kinematic Sunyaev-Zel'dovich effect with CMB-S4*, arXiv:2006.06594 [INSPIRE].

[20] CMB-S4 collaboration, *CMB-S4 science book, first edition*, arXiv:1610.02743 [INSPIRE].

[21] DESI collaboration, *The DESI experiment, a whitepaper for Snowmass 2013*, arXiv:1308.0847 [INSPIRE].

[22] LSST collaboration, *LSST: from science drivers to reference design and anticipated data products*, Astrophys. J. 873 (2019) 111 [arXiv:0805.2366] [INSPIRE].

[23] M. Münchmeyer, M.S. Madhavacheril, S. Ferraro, M.C. Johnson and K.M. Smith, *Constraining local non-Gaussianities with kinetic Sunyaev-Zel'dovich tomography*, Phys. Rev. D 100 (2019) 083508 [arXiv:1810.13424] [INSPIRE].

[24] Y.-Z. Ma, *Constraining the ionized gas evolution with CMB-spectroscopic survey cross-correlation*, Nucl. Phys. B 920 (2017) 402 [arXiv:1707.03348] [INSPIRE].

[25] N. Hand et al., *Evidence of galaxy cluster motions with the kinematic Sunyaev-Zel'dovich effect*, Phys. Rev. Lett. 109 (2012) 041101 [arXiv:1203.4219] [INSPIRE].

[26] DES, SPT collaboration, *Detection of the kinematic Sunyaev-Zel'dovich effect with DES year 1 and SPT*, Mon. Not. Roy. Astron. Soc. 461 (2016) 3172 [arXiv:1603.03904] [INSPIRE].

[27] F. De Bernardis et al., *Detection of the pairwise kinematic Sunyaev-Zel'dovich effect with BOSS DR11 and the Atacama Cosmology Telescope*, JCAP 03 (2017) 008 [arXiv:1607.02139] [INSPIRE].

[28] N.S. Sugiyama, T. Okumura and D.N. Spergel, *A direct measure of free electron gas via the kinematic Sunyaev-Zel'dovich effect in Fourier-space analysis*, Mon. Not. Roy. Astron. Soc. 475 (2018) 3764 [arXiv:1705.07449] [INSPIRE].

[29] ACTPOL collaboration, *Evidence for the kinematic Sunyaev-Zel'dovich effect with the Atacama Cosmology Telescope and velocity reconstruction from the Baryon Oscillation Spectroscopic Survey*, Phys. Rev. D 93 (2016) 082002 [arXiv:1510.06442] [INSPIRE].
[30] H. Tanimura, S. Zaroubi and N. Aghanim, *Direct detection of the kinetic Sunyaev-Zel'dovich effect in galaxy clusters*, arXiv:2007.02952 [inSPIRE].

[31] SDSS collaboration, *The seventh data release of the Sloan Digital Sky Survey*, Astrophys. J. Suppl. 182 (2009) 543 [arXiv:0812.0649] [inSPIRE].

[32] BOSS collaboration, *The clustering of galaxies in the SDSS-III Baryon Oscillation Spectroscopic Survey: baryon acoustic oscillations in the Data Releases 10 and 11 Galaxy samples*, Mon. Not. Roy. Astron. Soc. 441 (2014) 24 [arXiv:1312.4877] [inSPIRE].

[33] B. Reid et al., *SDSS-III Baryon Oscillation Spectroscopic Survey data release 12: galaxy target selection and large scale structure catalogues*, Mon. Not. Roy. Astron. Soc. 455 (2016) 1553 [arXiv:1509.06529] [inSPIRE].

[34] J. Jasche and B.D. Wandelt, *Bayesian physical reconstruction of initial conditions from large scale structure surveys*, Mon. Not. Roy. Astron. Soc. 432 (2013) 894 [arXiv:1203.3639] [inSPIRE].

[35] J. Jasche and G. Lavaux, *Physical bayesian modelling of the non-linear matter distribution: new insights into the Nearby Universe*, Astron. Astrophys. 625 (2019) A64 [arXiv:1806.11117] [inSPIRE].

[36] PLANCK collaboration, *Planck 2018 results. IV. Diffuse component separation*, Astron. Astrophys. 641 (2020) A4 [arXiv:1807.06208] [inSPIRE].

[37] SDSS collaboration, *A MaxBCG catalog of 13,823 galaxy clusters from the Sloan Digital Sky Survey*, Astrophys. J. 660 (2007) 239 [astro-ph/0701265] [inSPIRE].

[38] G. Lavaux, J. Jasche and F. Leclercq, *Systematic-free inference of the cosmic matter density field from SDSS3-BOSS data*, arXiv:1909.06396 [inSPIRE].

[39] J.-F. Cardoso, M. Martin, J. Delabrouille, M. Betoule and G. Patanchon, *Component separation with flexible models. Application to the separation of astrophysical emissions*, arXiv:0803.1814 [inSPIRE].

[40] B. Soergel, A. Saro, T. Giannantonio, G. Efstathiou and K. Dolag, *Cosmology with the pairwise kinematic SZ effect: calibration and validation using hydrodynamical simulations*, Mon. Not. Roy. Astron. Soc. 478 (2018) 5320 [arXiv:1712.05714] [inSPIRE].

[41] M. Li, R.E. Angulo, S.D.M. White and J. Jasche, *Matched filter optimization of kSZ measurements with a reconstructed cosmological flow field*, Mon. Not. Roy. Astron. Soc. 443 (2014) 2311 [arXiv:1404.0007] [inSPIRE].

[42] E. Rozo et al., *Constraining the scatter in the mass-richness relation of maxBCG clusters with weak lensing and X-ray data*, Astrophys. J. 699 (2009) 768 [arXiv:0809.2794] [inSPIRE].

[43] SDSS collaboration, *Cross-correlation weak lensing of SDSS galaxy clusters II: cluster density profiles and the mass-richness relation*, arXiv:0709.1159 [inSPIRE].

[44] K.M. Gorski et al., *HEALPix — A Framework for high resolution discretization, and fast analysis of data distributed on the sphere*, Astrophys. J. 622 (2005) 759 [astro-ph/0409513] [inSPIRE].

[45] A. Zonca et al., *healpy: equal area pixelization and spherical harmonics transforms for data on the sphere in python*, J. Open Source Softw. 4 (2019) 1298.

[46] K.S. Dawson et al., *The Baryon Oscillation Spectroscopic Survey of SDSS-III*, Astron. J. 145 (2012) 10.

[47] SDSS-III collaboration, *The eleventh and twelfth data releases of the Sloan Digital Sky Survey: final data from SDSS-III*, Astrophys. J. Suppl. 219 (2015) 12 [arXiv:1501.00963] [inSPIRE].

[48] Y.B. Zel'Dovich, *Gravitational instability: an approximate theory for large density perturbations*, Astron. Astrophys. 500 (1970) 13.
[49] T. Buchert, Lagrangian theory of gravitational instability of Friedman-Lemaitre cosmologies and the 'Zel'dovich approximation', Mon. Not. Roy. Astron. Soc. 254 (1992) 729.

[50] S. Tassev, M. Zaldarriaga and D. Eisenstein, Solving large scale structure in ten easy steps with COLA, JCAP 06 (2013) 036 [arXiv:1301.0322] [INSPIRE].

[51] V. Springel, The cosmological simulation code GADGET-2, Mon. Not. Roy. Astron. Soc. 364 (2005) 1105 [astro-ph/0505010] [INSPIRE].

[52] PLANCK collaboration, Planck intermediate results. LIII. Detection of velocity dispersion from the kinetic Sunyaev-Zeldovich effect, Astron. Astrophys. 617 (2018) A48 [arXiv:1707.00132] [INSPIRE].

[53] P.S. Behroozi, R.H. Wechsler and H.-Y. Wu, The rockstar phase-space temporal halo finder and the velocity offsets of cluster cores, Astrophys. J. 762 (2013) 109 [arXiv:1110.4372] [INSPIRE].

[54] A. Knebe et al., Galaxies going MAD: the galaxy-finder comparison project, Mon. Not. Roy. Astron. Soc. 428 (2013) 2039 [arXiv:1210.2578] [INSPIRE].

[55] PLANCK collaboration, Planck 2015 results. VIII. High frequency Instrument data processing: calibration and maps, Astron. Astrophys. 594 (2016) A8 [arXiv:1502.01587] [INSPIRE].

[56] PLANCK collaboration, Planck 2018 results. III. High frequency instrument data processing and frequency maps, Astron. Astrophys. 641 (2020) A3 [arXiv:1807.06207] [INSPIRE].

[57] PLANCK collaboration, Planck 2018 results. VI. Cosmological parameters, Astron. Astrophys. 641 (2020) A6 [arXiv:1807.06209] [INSPIRE].

[58] R.J. Thornton et al., The Atacama Cosmology Telescope: the polarization-sensitive ACTPol instrument, Astrophys. J. Suppl. 227 (2016) 21 [arXiv:1605.06569] [INSPIRE].

[59] J.E. Austermann et al., SPTpol: an instrument for CMB polarization measurements with the South Pole Telescope, Proc. SPIE Int. Soc. Opt. Eng. 8452 (2012) 84521E [arXiv:1210.4970] [INSPIRE].

[60] E. Schaan et al., The Atacama Cosmology Telescope: combined kinematic and thermal Sunyaev-Zel'dovich measurements from BOSS CMASS and LOWZ halos, arXiv:2009.05557 [INSPIRE].

[61] S. Amodeo et al., The Atacama Cosmology Telescope: modelling the gas thermodynamics in BOSS CMASS galaxies from kinematic and thermal Sunyaev-Zel'dovich measurements, arXiv:2009.05558 [INSPIRE].

[62] SPT-3G collaboration, SPT-3G: a next-generation cosmic microwave background polarization experiment on the South Pole Telescope, Proc. SPIE Int. Soc. Opt. Eng. 9153 (2014) 91531P [arXiv:1407.2973] [INSPIRE].

[63] J.L. Tinker et al., Toward a halo mass function for precision cosmology: the limits of universality, Astrophys. J. 688 (2008) 709 [arXiv:0803.2706] [INSPIRE].