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

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 $\simeq15\%$ at cluster scale. We observe evidence for the kSZ effect, at a significance of $\simeq2\sigma$. 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 baryonic content of galaxy clusters and galaxy groups.


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ν 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 directionn 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): where T 0 = 2.725 × 10 6 µ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 σ T and c denote the Thomson scattering cross-section and the speed of light in vacuum, respectively. The integral dl is performed along the line-of-sight (LOS)n, from the detector to the lastscattering 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 where v LOS (x,n) ≡ v(x) ·n denotes the velocity along the LOSn, and we have defined τ (x,n) = σ T dl n e (x) (1.3) to be the LOS projected optical depth. The r.h.s of Eq. (1.2) indicates that measurements of the kSZ signal constrain the product (v(x) ·n)τ . Assuming an external constraint on τ (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][9][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 τ . 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, latetime 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][14][15][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 τ 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 2 − 4σ evidence of the kSZ effect using the kSZ pairwise estimator [6,[25][26][27][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(x) is derived -assuming a certain cosmology with Hubble parameter H and cosmic linear growth rate f = d ln δ/d ln a -by solving the inversion of the linearized continuity equation in either real-space x [6,12,13,30] ∇ · v(x) = −aHf δ(x), (1.4) or redshift-space s [29], (1.5) where the DM density field δ(x) is simply obtained from a smoothed galaxy density field δ g (x) by assuming a local, linear bias relation of the form δ g (x) = b 1 δ(x). Specifically, Ref. [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 reconstruction 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 Sec. 3.2, in this paper we assume the same flat ΛCDM cosmology assumed by the reconstruction in [38], with Ω r = 0, Ω K = 0, Ω m = 0.2889, Ω b = 0.048597, Ω Λ = 0.7111, w = −1, n s = 0.9667, σ 8 = 0.8159, H 0 = 67.74 kms −1 Mpc −1 . The paper is structured as follows. In section Sec. 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 Sec. 3, we report our measurements of the kSZ effect from maxBCG clusters and their associated uncertainties in section Sec. 4. We then assert the robustness and significance of our measurement by means of different null tests in Sec. 5. Finally, we summarize our results and discuss relevant systematics as well as future improvements of kSZ measurement in section Sec. 6.

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.

Planck SMICA CMB map
Our choice of CMB data is the SMICA temperature map from the Planck 2018 release 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 ∆T obs is a convolution of the true anisotropy ∆Tincluding 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 2 , plus the instrumental noise ∆T instr , i.e.
where we have assumed the flat-sky approximation and replaced the three-dimensional vectors x andn by the two-dimensional vector θ 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 where the 5-arcmin full-width-half-maximum (FWHM) resolution translates into θ beam ≈ 5.0 arcmin/ 8 ln(2) ≈ 2.1 arcmin.

Signal extraction
To extract the kSZ signal from a CMB map, an aperture photometry (AP) filter of radius θ f is applied at the location of all clusters. The extracted flux ∆T θ f can then be expressed as a convolution of the observed flux ∆T obs with a radial weight function W θ f associated with that AP filter, i.e.
Specifically, the spherically symmetric weight function W θ f is given by This compensated filter is designed to reduce contributions from primary CMB anisotropies, as well as other low-redshift sources of contamination 3 , which vary on scales larger than θ f in the extracted flux. Thus, as θ f increases, so does contamination from these sources in the filtered flux. Consequently, the signal would be 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 4 while the uncertainty should increase due to increasing contamination. In our analysis, we measure the kSZ signal while varying the filter size θ f . Specifically for the individual-scale measurements of the kSZ signal in Sec. 3.1, Sec. 4.1, and Sec. 5.1, we adopt an adaptive aperture photometry (AAP) filter whose radius θ 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 pixelaveraged 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 Sec. 2.1 and Sec. 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.

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.

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 Note 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 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 × 10 14 M from our analysis (see App. 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 Fig. 1 a HEALPix map 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.

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 θ eff,i of the cluster. For example, assuming a spherically symmetric Gaussian profile for the electron gas, Eq. (1.3) becomes where τ 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 free = 1 in N e = f free f gas M 200,i /µ e m p . We further adopt a universal gas-mass fraction f gas = f b ≡ Ω b /Ω m = 0.16 following the cosmological baryon abundance and a mean particle weight per electron µ e = 1.17. Our expression for the integrated cluster optical depth defined in Eq. (2.6) then becomes (2.7)

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.

borg-SDSS3 reconstruction
We employ the non-linear velocity fields reconstructed using the forward-modeling algorithm borg [34,38] based on the observed large-scale structure as traced by LOWZ and CMASS galaxies from the SDSS3-BOSS DR13 release [46,47]. This reconstruction delivers a predicted cosmic velocity field within the SDSS volume that fully accounts for many physical and observational effects, including galaxy bias, light-cone effect, survey geometry, plus other selection and multiplicative systematic effects (see Section 2.2-2.6 of [38] for further details). The borg algorithm systematically explores the high-dimensional parameter space N dim = 256 3 consisting of the three-dimensional initial conditions at z ∼ 1000, augmented by the galaxy bias parameters. This is only possible thanks to the introduction and development of Hamiltonian Monte Carlo (HMC) sampling technique [48,49] for LSS inference [34,50] (see also [38] and references therein). The final density fields at z = 0 is linked to the initial conditions by a first-order Lagrangian perturbation theory (LPT) forward model. borg thus searches for three-dimensional matter distributions that are physically compatible with the constraints from the observed galaxy distribution. 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. The setup of the inference is given as follows. The initial conditions are generated on a comoving grid consisting of 256 3 cells and covering a comoving volume of 4000 3 h −3 Mpc 3 , which corresponds to a grid resolution of L grid = 15.624 h −1 Mpc. The SDSS3-BOSS data are projected in a sub-volume with the observer located at x = {200, 0, −1700} h −1 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 ∼ 1000 samples. Here, to guarantee that we only include MCMC samples after convergence, we further 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 [51] (tcola hereafter), and the cloudin-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 part = 1024 3 . These are used for the estimation of cluster LOS velocities as well as their uncertainties.
We additionally generate a GADGET-2 [52] simulation at resolution of N part = 2048 3 from initial conditions specified by the sample s = 9000 of borg-SDSS3 reconstruction (see App. D for details). We use this full N-body, high resolution simulation to speficically:

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: where v LOS L is the large-scale LOS bulk-flow estimated from the borg-SDSS3 reconstruction posterior while LOS S is the unresolved small-scale LOS velocity. We further assume that, for all clusters:  Fig. 3) in which v LOS nbody,i and v LOS tcola,i 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.

Data model and kSZ posterior
Given an inferred ensemble of the large-scale LOS velocity of each galaxy cluster i, v LOS L,i , 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 Sec. 2.2.1).

kSZ likelihood: single angular scale
Let us express our data model as where we have introduced α θ f as the amplitude of the kSZ signal from the large-scale bulk flow of cluster i. Here and below, we simply adopt τ i = τ 0,i . For cases where the filter size is smaller than that of the cluster, this means letting α θ f absorb all specific details about the cluster gas profile. We expect to measure a value of α θ 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. θ f ≥ θ 2 vir + θ 2 beam , corresponds to the expectation from the simple model of optical depth discussed in Sec. 2.2.2. The exact value of α depends on that of f 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 α as f 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. Both the signal amplitude and noise are functions of the AP filter size θ f (or AAP filter scale ϕ f ) -as indicated by the superscript θ f . However, for the sake of readability, we will omit the superscript θ 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 We now seek to construct a posterior distribution for α, marginalized over N velocity realizations, i.e.
where we have used {x i } ≡ {τ i v LOS,i /c} and introduced the prior on α explicitly as P(α).
The borg-SDSS3 reconstruction provides a sampled approximation where δ D ({x i }) denotes the Dirac delta distribution. We thus can rewrite Eq. (3.4) as 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 λ s , which is given by (see App. B for a detailed derivation) , (3.8) in which and (3.10) Note that the weights λ s give preference to better-fitting realizations of the peculiar velocity field, and N s λ s = 1. For simplicity, in what follows, we assume a uniform prior on α such that P(α) = 1 for all borg-SDSS3 velocity realizations. The posterior mean estimate of α is then given by (see App. B for details) while its variance is given by 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(α > 0) and no-kSZ hypothesis P(α ≤ 0): 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 σ 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 normal (x ≤ −sign.) = P(α ≤ 0) , (3.14) and sign. would be our significance measured in the familiar unit of σ. Note that Eq. (3.11) is identical to Eq. (9) in [29] if one takes all λ 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 Fig. 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 -10 − 20%. It is worth pointing out also that different realizations yield different estimates of α s itself. A combination of, say, high α s and low σ α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 Tab. 1 of [53]. 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 App. 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 λ s=9000 → 1 while the other weights go to zero, correctly singling out the original borg-SDSS3 sample used to generate the mock signal template.

kSZ likelihood: combined signal
To combine measurements at different filter sizes θ f , it is necessary to modify Eqs. (3.1)-(3.7) to include the cluster gas profile f (θ f ) as 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 App. D). We show below the measurement of f θ f at locations of corresponding DM halos identified by Rockstar halo finder 6 [54,55], 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 θ f does not depend on the cluster effective apparent size.
where we have again omitted the superscript θ f and instead expressed quantities that depend on θ f as vectors to stress the fact that all measurements at different filter scales are now combined in Eq. (3.16). Note also that α is now simply a scalar instead of a function of θ f . Eq. (3.2) can then be rewritten as where the covariance matrix C i for cluster i is given by (3.18) 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 τ 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 instrumentspecific mocks to estimate [56,57]. 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 [56,57] (see details in App. 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 [58]. In the flat-sky limit, where W (lθ f ) is the Fourier transform of the AP filter The CMB covariance matrix in Eq. (3.18) is then explicitly given by (see App. C) where C CMB denotes the Planck 2018 ΛCDM best fit angular power spectrum. We show in Fig. 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.
The expressions of the posterior mean α n and its variance σ 2 α are similar to those in Eq. (3.11) and Eq. (3.12) with modifications to µ s and ω s as described in App. B.  [10,11,12] bins. As indicated in the right panel of Fig. 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 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 µ −0.005 in z photo (see Fig. 7). Note that when drawing from the Gaussian distribution in Eq. (3.25), we limit the range of z r photo,i in all redshift realizations to within [0.05, 0.5]. So our Gaussian distributions of photometric redshift error are actually truncated.

Modeling photo-z uncertainty
We then introduce an additional sum over all N r redshift realizations in Eq. (3.7), Note that, in the case of the AAP filter, for all redshift realizations, we keep the filter size in Eq.

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 α ϕ f nover velocity and redshift realizations (cf. Eq. (3.11)) -and the associated significance (cf. Eq. (3.14)).

Individual-scale signal measurements
We show in Fig. 8 the inferred value of α ϕ f n as a function of the AAP filter scale ϕ f , using spec-z set alone (top panel) or photo-z set and both sets (bottom panel). We emphasize again that the 1σ uncertainty shown here as shaded band and reported in Tab. 1 includes also uncertainties in the reconstructed large-scale velocity field.
As can be noted from Fig. 8, we obtain consistent results between the two datasets specz and photo-z, which can be combined as shown in red in the bottom panel of Fig. 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 Fig. 8 show that most of the information come 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 loose its constraining power.
In addition, we provide the significance at each filter scale for each dataset in Tab. 1. All sets show peaks of significance at ϕ = 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 Fig. 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 λ s in the three panels on the right of Fig. 9, which show λ s for the combined set, are consistent with those of λ s in the leftmost and middle columns, which respectively show λ 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. The posterior mean of the kSZ signal amplitude α n ± 1σ α 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.

Combined signal measurement
We show in Fig. 10 our measurements of combined signal amplitude α θ f n (cf. Eq. (3.11)), as a cumulative function of the AP filter radius θ 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 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  Figure 9: Distributions of the mixture weight λ 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 ϕ = 0.9. Note that we do not sample the redshifts of clusters in the spec-z set, which is why the λ n are evenly distributed among redshift realizations (y-axes) in the leftmost column. 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 Tab. 2. As in Tab. 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 [59,60], 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, 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.

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 spec-z photo-z spec-z + photo-z α ϕ f n sign. α ϕ f n sign. α ϕ f n sign.  Table 2: The posterior mean of the kSZ combined signal amplitude α n ± 1σ α 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 θ f = [3,7] arcmin.
approximated by the simple ratio S/N= α ϕ f n /σ α . Below, we employ this approximation to ease the numerical computation of the significance in these null tests.

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 Fig. 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 Fig. 12 the histograms of S/N for increasing individual filter scale from ϕ = 0.9 to ϕ = 1.1. The histogram of ϕ = 0.9 (top row, left panel) appears consistent with our reported significances of 1.7 for the spec-z dataset (cf. Tab. 1).

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 Fig. 11. Our null test using the Planck simulations also recovers S/N ratios consistent with zero-mean Gaussian distributions. We show in Fig. 12   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, in Sec. 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 α can be easily introduced for cases wherein the cluster gas profile is known from complementary measurements, e.g. X-ray, tSZ, etc. 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 Sec. 4, moderate evidence of the signal at 2σ -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 Fig. 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 [61].
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 Sec. 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.

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 . We have found that, specifically for our sub-sample 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×10 13 M . This is demonstrated in Fig. 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 × 10 14 M . 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 LOS = 300 kms −1 for all clusters. Both panels of Fig. 13 show signs of a significant tSZ contamination when including clusters more massive than 8.5 × 10 13 M .

B The kSZ likelihood: mixture weights, MAP mean and variance of the kSZ signal amplitude
Here, we first derive the individual mixture weight λ s for each component of the Gaussian mixture distribution of large-scale bulk-flow amplitude α s in Sec. 3. Let us denote and Then, the likelihood on the r.h.s. of Eq. (3.7) can be re-written as We further denote  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 α includes also the uncertainty from the velocity reconstruction. We next compute ln P ({A i }|α, {B s i }) for measurements at all θ f scales in a similar fashion to how we arrived at Eq. (B.11):

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 Sec. 3 (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 ∆T at 1000 random points in each of the 1000 SMICA2018-like CMB mocks in Fig. 14. The analytical calculatrion using the Planck 2018 best-fit ΛCDM power spectrum is in extremely good agreement with the numerical estimate using SMICA2018-like CMB realizations.

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 [52] 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 8 [54,55] 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 speficically taken from sample 9000. The high resolution allows us to achieve a correct halo mass function (HMF) down to M h = 2 × 10 13 h −1 M at redshift z = 0.23, as shown in Fig. 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]:  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.