© H. Tanimura et al. 2021

1. Introduction

Galaxy clusters are the largest gravitationally bound structures in the Universe and have been used as probes of cosmology and astrophysics. These massive objects imprint their signature on the cosmic microwave background (CMB) through the Sunyaev-Zel’dovich (SZ) effect (Sunyaev & Zeldovich 1970, 1972, 1980). The SZ effect is caused by the scattering of CMB photons by hot, ionized plasma in the intracluster medium (ICM) giving rise to a change in the CMB temperature. The SZ effect can be classified into two contributions: the thermal Sunyaev-Zel’dovich (tSZ) effect and the kinetic Sunyaev-Zel’dovich (kSZ) effect (e.g., Diaferio et al. 2000).

The tSZ effect is caused by the scattering of CMB photons by free electrons with thermal motions in galaxy clusters, resulting in a characteristic spectral distortion of the CMB blackbody spectrum. This effect allows us to trace gas pressure in these objects, and has been well characterized through measurements of individual massive clusters (e.g., Plagge et al. 2010; Bonamente et al. 2012; Sayers et al. 2013a; Planck Collaboration V 2013; Romero et al. 2015) and statistical measurements of massive to low-mass systems (e.g., Planck Collaboration XI 2013; Greco et al. 2015; Vikram et al. 2017; Hill et al. 2018; Lim et al. 2018; Tanimura et al. 2019, 2020).

The kSZ effect, in turn, is caused by the scattering of CMB photons off the electrons due to the peculiar motion of a galaxy cluster, leading to a Doppler shift of the CMB blackbody spectrum. While kSZ effect is elusive due to its small amplitude and identical spectral shape to the CMB spectrum, kSZ has great potential in constraining both cosmological and astrophysical models. From a cosmological point of view, peculiar velocities of galaxy clusters are found through measurements of the kSZ effect and they allow the amplitude of the growth rate of density fluctuations to be estimated. In coming years, it is expected that kSZ measurements by the advanced Atacama Cosmology Telescope (AdvACT; Henderson et al. 2016) and the third-generation South Pole Telescope (SPT-3G; Benson et al. 2014) will constrain the growth rate of structures to sub-percentage precision at z < 1 (Alonso et al. 2016), which then help to constrain models of dark energy (Bhattacharya & Kosowsky 2008; Ma & Zhao 2014), modified gravity (Mueller et al. 2015a; Bianchini & Silvestri 2016), and massive neutrinos (Mueller et al. 2015b). From an astrophysical point of view, there is debate over whether a significant fraction of diffuse gas is present around halos as a circumgalactic medium, or whether the gas, once expelled because of feedback processes such as star formation, supernovae, and active galactic nuclei (AGNs), is never accreted onto the halos (e.g., Planck Collaboration XI 2013; Anderson et al. 2015; Le Brun et al. 2015). As the kSZ effect is sensitive to the virialized gas and also to the gas surrounding halos, independent of the gas temperature (unlike the tSZ effect), it is well suited to studying these baryons through the distribution of gas around galaxy clusters.

The kSZ signal has so far been detected for a few individual systems (e.g., Sayers et al. 2013b; Adam et al. 2017) or by statistical measurements based on the pairwise method (e.g., Hand et al. 2012; Hernández-Monteagudo et al. 2015; Planck Collaboration Int. XXXVII 2016; Soergel et al. 2016; De Bernardis et al. 2017) or cross-correlation method (Hill et al. 2016). However, the reported significance of these detections is limited to ∼2−4σ. To extract the kSZ signal, most of these statistical measurements rely on a matched filter (e.g., Soergel et al. 2016; Lim et al. 2020) or aperture photometry (e.g., Hernández-Monteagudo et al. 2015; Planck Collaboration Int. XXXVII 2016; Schaan et al. 2016; De Bernardis et al. 2017; Sugiyama et al. 2018). However, due to the large uncertainty in the gas density profile in galaxy clusters, the kSZ signal extracted by a matched filter can be biased if the assumed profile is incorrect (Ferraro & Hensley 2015). Aperture photometry has less dependence on the density profile, however the extracted signal level changes depending on the aperture size because of the varying signal-to-noise ratio (S/N).

In the present analysis, to extract the kSZ signal, we use a stacking method without making any assumptions regarding the spatial distribution of gas around galaxy clusters or setting any aperture. To avoid cancelation of the kSZ signals when all the clusters are simply stacked, we estimate the peculiar velocities of the galaxy clusters along the line of sight (LOS) with which the sign of the kSZ signals is attributed during the stacking. The paper is organized as follows. Section 2 summarizes data sets used in our analyses. Section 3 explains the velocity reconstruction of galaxy clusters. Section 4 explains the stacking method to extract the kSZ signals. Possible systematic errors in our measurements are discussed in Sect. 5. The interpretation of the measurements is presented in Sect. 6. We end this paper with discussions and conclusions in Sects. 7 and 8.

Throughout this work, we adopt the ΛCDM cosmology in Komatsu et al. (2011) with Ω_m = 0.272, Ω_b = 0.046, and H₀ = 70.4 km s⁻¹ Mpc⁻¹. The cosmological parameters are used to estimate peculiar velocities of galaxy clusters and distances to them, and our result may depend on the assumed cosmological parameters. However, we perform our data analysis with Planck cosmology in Planck Collaboration XIII (2016) and obtain consistent results. All masses are quoted in units of solar mass, and M_Δ is mass enclosed within a sphere of radius R_Δ such that the enclosed density is Δ times the critical density at redshift z. Uncertainties are given at the 1σ confidence level.

2. Data

2.1. Galaxy cluster catalog

Wen et al. (2012) and Wen & Han (2015) identified a total of 158 103 Wen-Han-Liu galaxy groups and clusters from the Sloan Digital Sky Survey (SDSS) galaxies in the redshift range between 0.05 and 0.8, of which 89% have spectroscopic redshifts (hereafter WHL galaxy clusters). The masses of the WHL galaxy clusters were estimated from their total luminosity, and these were calibrated by the masses of 1191 clusters using X-ray or tSZ measurements. In our study, we use WHL galaxy clusters that have spectroscopic redshifts at 0.25 <  z <  0.55 and masses of M₅₀₀ >  10^13.5 h⁻¹ M_⊙. In addition, after selections in Sects. 3 and 4.2, the total number of clusters used in our analysis is 30 431.

2.2. Galaxy catalog

The SDSS galaxies in Reid et al. (2016) are composed of 953 193 galaxies in the northern Galactic hemisphere and of 372 542 in the southern Galactic hemisphere, combining the Baryon Oscillation Spectroscopic Survey (BOSS) LOWZ galaxies and constant-mass (CMASS) galaxies galaxies. The completeness of the galaxies is stated to be 99% for CMASS and 97% for LOWZ. Spectroscopic information is available for all the galaxies and their redshifts extend up to z ∼ 0.8. We used the galaxies to compute the galaxy density field, which was then used to estimate peculiar velocities of the WHL clusters. To avoid a bias in the velocity calculation due to the magnitude limit at different redshifts, we limit our analysis to the range of 0.25 <  z <  0.55, in which the number density of the galaxies in the survey volume is fairly flat (see Fig. 11 in Reid et al. 2016).

2.3. Planck maps

Planck produced all-sky maps in nine frequency bands from 30 to 857 GHz with the angular resolutions ranging from 31′ to 5′. In our study, for the detection of the kSZ signal, we mainly used the 217 GHz frequency map from the Planck 2018 data release (Planck Collaboration III 2020), because it corresponds to the null frequency of the tSZ effect. We also used, for comparison, the Planck CMB maps produced using four different component-separation methods (Planck Collaboration IV 2020): COMMANDER (Optimal Monte-Carlo Markov chAiN Driven EstimatoR; Eriksen et al. 2004, 2008), NILC (Needlet Internal Linear Combination; Delabrouille et al. 2009), SEVEM (Spectral estimation via expectation maximisation; Martínez-González et al. 2003; Leach et al. 2008; Fernández-Cobos et al. 2012), and SMICA (Spectral Matching Independent Component Analysis; Delabrouille et al. 2003; Cardoso et al. 2008). These maps were provided in HEALpix¹ format (Górski et al. 2005) with a pixel resolution of N_side = 2048 (∼1.7 arcmin). To minimize the Galactic and extragalactic contamination, we applied the mask produced by the Planck team for the analysis of the CMB temperature maps, which masks the region around the Galactic plane and the point sources detected at all the Planck frequencies (see Table C.1 in Planck Collaboration IV 2020). In addition, we used a more robust point-source mask, masking radio and infrared sources, which is used for the analysis of the Compton y maps (Planck Collaboration XXII 2016). Combining these two masks excludes ∼50% of the sky.

2.4. Magneticum simulation

The Magneticum simulations are one of the largest cosmological hydrodynamical simulations (Hirschmann et al. 2014; Dolag 2015), and are based on the standard ΛCDM cosmology from Komatsu et al. (2011) with Ω_m = 0.272, Ω_b = 0.046, and H₀ = 70.4 km s⁻¹ Mpc⁻¹. They provide several public datasets² (Ragagnin et al. 2017), of which we use the simulated kSZ light-cone map with an area of ∼1600 deg² (hereafter Magneticum light cone). The corresponding cluster catalog with M₅₀₀ >  3 × 10¹³ h⁻¹ M_⊙ at z ≲ 2 is also provided (Dolag et al. 2016; Soergel et al. 2018). In addition, they also provide the post-processed data of the galaxy catalog and the cluster catalog from the full simulations with a cube box of 640³h⁻¹ Mpc. We use the post-processed data of “Box2b_hr” at z ∼ 0.42 (hereafter Magneticum snapshot), corresponding to the median redshift in our analysis range.

3. Velocity reconstruction

The peculiar velocity of a galaxy cluster at the position, x, can be derived in the linear regime by

$$\begin{array}{c}\hfill \mathit{v}(\mathit{x})=\frac{f(\mathrm{\Omega })aH(a)}{4\pi }{\displaystyle \int }{\mathrm{d}}^3\mathit{y}\delta (\mathit{y})\frac{\mathit{y}-\mathit{x}}{{|\mathit{y}-\mathit{x}|}^3},\end{array}$$(1)

where a is the scale factor, H(a) is the Hubble parameter, f(Ω) is the linear velocity growth rate given by $f(\mathrm{\Omega })\simeq {\mathrm{\Omega }}_{\mathrm{m}}^{0.545}$ (Lahav et al. 1991; Wang & Steinhardt 1998; Linder 2005; Huterer & Linder 2007; Ferreira & Skordis 2010), and δ(y) is the overdensity of matter at the position y.

To compute the matter density field, δ(y), we used the galaxies in Sect. 2.2. Their redshift distances and J2000.0 coordinates were transformed into the comoving Cartesian coordinates:

$$\begin{array}{cc}& X=r(z)cos\alpha cos\delta ,\hfill \\ & Y=r(z)sin\alpha cos\delta ,\hfill \\ & Z=r(z)sin\delta ,\hfill \end{array}$$(2)

where α and δ refer to the J2000.0 right ascension and declination, respectively, and r(z) is the comoving distance at redshift z, with which a number density of galaxies was computed in the SDSS survey field. The galaxy density field can be connected to the matter density field by δ_g = b δ through the galaxy bias, b. The galaxy bias was studied for the LOWZ galaxies (Parejko et al. 2013) and CMASS galaxies (White et al. 2011; Nuza et al. 2013; Rodríguez-Torres et al. 2016), showing b ∼ 2 at scales larger than 10 h⁻¹ Mpc. Since a much larger scale of ∼240 h⁻¹ Mpc, as later described, was considered in our study to calculate the peculiar velocity, we used b = 2. The estimated galaxy bias has a relatively large uncertainty at large scales, which is not included in our analysis. While we did not find a significant bias in the reconstructed peculiar velocities using b = 2 as shown in the following Fig. 1, a more precise estimate of galaxy bias would be needed for a more precise measurement of peculiar velocity.

Fig. 1. x-axis: true LOS velocities of simulated galaxy clusters in the Magneticum simulations. y-axis: LOS velocities of the same galaxy clusters computed using Eq. (1) with the galaxy bias of b = 2 when the RSD effect is not included (left) and included (right).

In practice, we placed a galaxy cluster at the center in a cubic box of 240³h⁻¹ Mpc, in which the box was divided by grid cells of 5³h⁻¹ Mpc. If the cubic box for a galaxy cluster reaches the edge of the SDSS survey, we removed the galaxy cluster from our catalog because it may bias the velocity estimate. We then placed galaxies around the galaxy cluster in the box cells and calculated the overdensities of the galaxies. While finer grids may be better for the velocity estimate, we determined the size so that the grid size is large enough compared to the length expected from redshift-space distortion (RSD): The RSD distortion for a typical velocity of 300 km s⁻¹ is ∼3 h⁻¹ Mpc at the redshift of interest in our analysis. Box cells were then smoothed by a Gaussian kernel of 2 h⁻¹ Mpc to remove sharp grid edges. By applying Eq. (1) to this cubic box, the peculiar velocity of a galaxy cluster was obtained.

We checked our velocity estimates of galaxy clusters using the Magneticum snapshot simulation by comparing the velocities of the simulated galaxy clusters computed using Eq. (1) with the true velocities. To reproduce “real” data in the simulation, we removed the simulated galaxy clusters with M₅₀₀ <  10^13.5 h⁻¹ M_⊙ (minimum mass of the WHL clusters we use), and also removed the simulated galaxies with M_* <  1.8 × 10¹¹ h⁻¹ M_⊙, so that the number density of galaxies is the same as the real data. In this simulation box, we calculated the 3D velocities of the simulated galaxy clusters and compared the velocities in the “LOS” direction, which was defined by placing an observer at the corner of the simulation box. From the comparison, we find that our velocity estimates are well correlated with the true ones as shown in the left panel in Fig. 1 with an uncertainty of ∼180 km s⁻¹ following a Gaussian distribution.

In addition, we checked the effect of RSD on the velocity estimates using the Magneticum snapshot simulation by including the RSD in the simulated galaxies. For simplicity, we did not re-identify galaxy clusters in the redshift space, but used the same galaxy clusters with their positions redshifted. By again comparing the computed velocities of the simulated galaxy clusters with the true velocities, we find that our velocity estimates are well correlated with the true ones as shown in the right panel of Fig. 1, with a larger uncertainty of ∼260 km s⁻¹ following a Gaussian distribution. We also find a negative bias of ∼40 km s⁻¹ relative to the true values, but this bias is not significant given the current uncertainty.

4. Analysis

To detect the kSZ signal, we applied a stacking method for the Planck HFI 217 GHz map or Planck CMB maps. However, because there is equal probability that a cluster will have positive or negative LOS velocity, the associated kSZ signal from clusters cancels out by a simple stacking. Therefore, to avoid the cancelation of equally likely positive and negative kSZ signals, we first performed the stacking after separating galaxy clusters depending on the directions of the LOS velocities estimated in Sect. 3.

In this paper, we define “positive” LOS direction as a radial direction from us: a positive motion is the motion moving away from us and a negative motion is the motion approaching us. It follows that when a galaxy cluster has a positive motion, the CMB is redshifted, resulting in a negative kSZ signal. On the other hand, when a galaxy cluster has a negative motion, the CMB is blueshifted, resulting in a positive kSZ signal.

The stacking was then also performed for each cluster weighted by the LOS velocity. A positive kSZ signal weighted by a negative LOS velocity then has a negative signal and a negative kSZ signal weighted by a positive LOS velocity also has a negative signal. This allows us to set the kSZ signal to a negative value, while other components are canceled out by the positive or negative LOS velocity with equal probability.

4.1. Filtering the cosmic microwave background

The amplitude of the kSZ signal around galaxy clusters is of the order of ∼1 uK, and is dominated by the primordial CMB fluctuations of the order of ∼50 uK. To extract the kSZ signal present at the cluster scales, that is, small angular scales, we applied a spatial filter to the Planck maps. The angular size of the virial radius of the WHL galaxy clusters is 2.2−10.5 arcmin. We therefore filtered out angular scales above 15 arcmin (ℓ ∼ 720) in Fourier space using a smooth function with its response of one below 15 arcmin (ℓ ∼ 720) and zero above 30 arcmin (ℓ ∼ 360). The effect of the filter is shown in Fig. 2 for a simulated power spectrum of the primordial CMB fluctuations. With this filter, the standard deviation of the primordial CMB fluctuations is reduced to ∼40 uK.

Fig. 2. Filter applied to the Planck maps. The filter is a smooth function with its response of one below 15 arcmin (ℓ ∼ 720) and zero above 30 arcmin (ℓ ∼ 360). With this filter, the primordial CMB fluctuations in black are suppressed to the ones in red.

4.2. Stacking

We stacked the filtered Planck maps at the positions of the WHL galaxy clusters and constructed the stacked radial profile. In practice, we placed each galaxy cluster at the center of a two-dimensional grid in “scaled” angular distance in the range of −10 <  θ/θ₅₀₀ <  10, divided into 10 × 10 bins, where θ₅₀₀ is the angular radius of a galaxy cluster calculated with R₅₀₀ provided in the catalog. The Planck maps were scaled accordingly and data were placed on the two-dimensional grids, while data in the masked region were not used. In this process, if more than 20% of the region within 10 × θ₅₀₀ around a galaxy cluster was masked, the galaxy cluster was removed from our catalog, because a large mask may bias our measured profile. We repeated this for the selected 30 431 galaxy clusters and computed the stacked radial profile. In the upper left panel of Fig. 3, the stacked radial profile of the galaxy clusters using the Planck HFI 217 GHz map is shown in black. The stacked profile includes all the Galactic and extragalactic emissions in the frequency in which the cosmic infrared background (CIB) is dominant at small scales as shown in Appendix A. In this step, no velocity weighting is applied and hence positive and negative kSZ signals are canceled out.

Fig. 3. Upper panels, left: stacked radial profile around the 30 431 WHL clusters with the Planck temperature map at 217 GHz (black) with the 1σ uncertainty estimated by a bootstrap resampling (gray). The stacking is performed separately for the 14 093 WHL clusters with negative LOS motions (black dashed line) and for the 16 338 WHL clusters with positive motions (black dash-dotted line). Middle: same as left panel but with an additional LOS velocity cut of 100 km s−1. Right: same as left panel but with an additional LOS velocity cut of 200 km s−1. Lower panels: positive kSZ radial profile (blue) extracted by taking the difference between the black dashed and black lines, and negative kSZ radial profile (red) extracted by taking the difference between the black and black dash-dotted lines.

To extract the kSZ signals of the galaxy clusters, we stacked the galaxy clusters with positive (16 338 clusters) and negative (14 093 clusters) LOS motions separately. As expected, the original stacked profile with all the galaxy clusters (30 431 clusters) was separated into two: one with stacked positive kSZ signal and the other with stacked negative kSZ signal as shown in Fig. 3. As, in these profiles, the same amount of Galactic and extragalactic components are included on average, we can extract the kSZ signals by taking the differences of the two separated profiles relative to the original one, shown in Fig. 3.

To confirm whether or not the separated signals originate from the kSZ, we checked the correlation between the amplitude of the separated signals and the LOS velocity, which is only expected for the kSZ but not for other components such as CMB, CIB, and tSZ. In the upper middle and right panel of Fig. 3, the stacking was performed for the galaxy clusters with the absolute LOS velocity larger than 100 km s⁻¹ or 200 km s⁻¹, and the separated signals are shown in blue and red in the lower panels. The results show that the amplitude of the separated signals increases along with the amplitude of the velocity cut. This trend clearly supports that the separated signals originate from the kSZ.

We assessed the uncertainties of the stacked profile through bootstrap resampling. We drew a random sample of the galaxy clusters with replacement and re-calculated one stacked profile for the new set of galaxy clusters. We repeated this process 1000 times and produced the bootstrapped 1000 stacked profiles, with which the covariance between different radial bins was computed. In Fig. 3, the 1σ statistical uncertainty is represented by the width of the lines, which is a square root of diagonal terms of the covariance matrix.

4.3. Comparison with hydrodynamic simulations

We compared our measured kSZ signal with predictions from the Magneticum hydrodynamic simulation. For comparison, we performed the same stacking analysis with the simulated kSZ map and galaxy clusters from the Magneticum light-cone simulation as was performed with the real data. The mass and redshift distributions of the simulated galaxy clusters are matched to the real data. However, because a galaxy catalog is not provided in the Magneticum light-cone simulation, we added uncertainties to the LOS velocities of the simulated galaxy clusters randomly using a Gaussian function with the standard deviation of 260 km s⁻¹ (see Sect. 3). Based on the LOS velocities including the uncertainties, the simulated galaxy clusters were separated according to their positive or negative motions and stacked with the simulated kSZ map separately. (Some cancellation of the kSZ signal due to the uncertainties of the LOS velocities are included in the kSZ profile from the simulation as in the real data.) The result is shown as a green dashed line in Fig. 4 and is compared to the data profile. In the left panel of Fig. 4, we can see that our measured kSZ profiles with positive and negative LOS motions are consistent with the ones from the simulation. We also applied additional LOS velocity cuts of 100 km s⁻¹ and 200 km s⁻¹ to the simulated galaxy clusters in the middle and right panel, respectively. Our measured kSZ profiles show a similar amount of correlation between the amplitude of the kSZ signal and the LOS velocity to that predicted from the simulation.

Fig. 4. Positive kSZ radial profile (blue) around the 14 093 WHL clusters with the Planck temperature map at 217 GHz with the 1σ uncertainty estimated by a bootstrap resampling, and negative kSZ radial profile (red) around the 16 338 WHL clusters, compared to positive kSZ radial profiles (green dashed) and negative kSZ radial profile (green dash-dotted) from the Magneticum hydrodynamic simulations. Middle: same as left panel with an additional LOS velocity cut of 100 km s−1. Right: same as left panel with an additional LOS velocity cut of 200 km s−1.

4.4. Null tests and significance

To estimate the significance of the detected kSZ signal from the 30 431 WHL galaxy clusters, we stacked the filtered Planck maps at the positions of the clusters in the same way, as described above, with each cluster weighted by the LOS velocity as

$$\begin{array}{c}\hfill T(R)=\frac{{\sum }_i{T}_i(R)\times v_{i,\mathrm{LOS}}/{\sigma }_i^2}{{\sum }_i|v_{i,\mathrm{LOS}}|/{\sigma }_i^2},\end{array}$$(3)

where T_i(R) is the temperature value of the ith cluster at the radial distance, R, v_i, LOS is the LOS velocity of the ith cluster, and σ_i is the variance of temperature values within the region we consider (10 × θ₅₀₀) centered on the ith cluster. In this process, the kSZ signals with positive and negative LOS velocities end up with the same sign: a positive kSZ signal weighted by a negative LOS velocity has a negative signal and a negative kSZ signal weighted by a positive LOS velocity also has a negative signal. Therefore, the kSZ signals can be stacked without suffering any cancelation, while other components are canceled out. In addition, a cluster with a low LOS velocity (i.e., a weaker kSZ signal) is underweighted in the stacking. The stacked radial profile with the additional weight is shown in Fig. 5 with the uncertainties estimated by bootstrap. An oscillating angular pattern is seen in the stacked radial profile. This is the result of the convolution of kSZ and CMB with the filter in Fig. 2. The pattern is also seen in the null-test profile in Fig. 6 and the β-model profile in Fig. 7. We refer to this velocity-weighted kSZ profile from now on.

Fig. 5. Velocity-weighted kSZ radial profile around the 30 431 WHL clusters with the Planck temperature map at 217 GHz (black) and with the 1σ uncertainty estimated by a bootstrap resampling, compared to the kSZ radial profiles with the Planck CMB maps produced by different component separations of NILC (green), SMICA (yellow), SMICA-noSZ (purple), SEVEM (cyan), and COMMANDER (orange).

Fig. 6. Velocity-weighted kSZ radial profile around the 30 431 WHL clusters with the Planck temperature map at 217 GHz (black) compared to three null tests. Left panel: clusters are displaced at random positions on the Planck map, and then stacked. This process is repeated 1000 times and the mean of the 1000 random samples is computed (cyan). The 1σ uncertainty is estimated by computing a standard deviation of the 1000 random samples. Middle panel: results from cluster stacking after randomly shuffling the LOS velocities of the clusters; this process is repeated 1000 times and the mean of the 100 velocity-shuffled profiles is computed (yellow). The 1σ uncertainty is estimated by computing a standard deviation of the 1000 velocity-shuffled profiles. Right panel: clusters are stacked with a noise map produced by $(T_{217}^{\mathrm{HM}1}-T_{217}^{\mathrm{HM}2})/2$, where $T_{217}^{\mathrm{HM}1(2)}$ is the half mission 1(2) Planck map at 217 GHz (green). The 1σ uncertainty is estimated by a bootstrap resampling.

Fig. 7. Velocity-weighted kSZ radial profile around the 30 431 WHL clusters with the Planck temperature map at 217 GHz (black), fitted with the β model (red).

First, we estimated the excess of the measured velocity-weighted kSZ profile to the null hypothesis. The signal-to-noise ratio (S/N) can be estimated with

$$\begin{array}{c}\hfill S/N=\sqrt{{\chi }_{\mathrm{data}}^2-{\chi }_{\mathrm{null}}^2},\end{array}$$(4)

where

$$\begin{array}{cc}& {\chi }_{\mathrm{data}}^2={\displaystyle \sum _{i,j}{T}_{\mathrm{data}}{(R_i)}^{\mathrm{T}}(C_{\mathit{ij}}^{-1})T_{\mathrm{data}}(R_j),}\hfill \end{array}$$(5)

$$\begin{array}{cc}& {\chi }_{\mathrm{null}}^2={\displaystyle \sum _{i,j}{T}_{\mathrm{null}}{(R_i)}^{\mathrm{T}}(C_{\mathit{ij}}^{-1})T_{\mathrm{null}}(R_j),}\hfill \end{array}$$(6)

where T_data(R_i) is the temperature value at the R_i bin of the data kSZ profile, T_null(R_i) is the temperature value at the R_i bin under the null hypothesis, which is zero, and C_ij is the covariance matrix of the data profile, estimated by the bootstrap resampling. By measuring the kSZ signal up to 4 × θ₅₀₀, the S/N value was estimated to be 3.5σ.

We also performed a Monte Carlo-based null test to assess the significance of our measurements. In the null test, we displaced the centers of the galaxy clusters at random positions on the sky and then the Planck maps were stacked at these random positions. We repeated this 1000 times to assess the rms fluctuations of the foreground and background signals. The result shows that the average of the null-test profiles is consistent with zero, as shown in cyan in the left panel of Fig. 6 with the rms fluctuations, suggesting that our measurements are unbiased. The mean null-test profile shows the same coherent angular pattern as seen in the data profile in Fig. 5, which is due to the convolution of CMB with our filter in Fig. 2. Using the null-test sample, the significance of the measured kSZ signal can also be estimated with Eq. (4) by replacing T_null(R_i) and C_ij with their corresponding values from the null test. The S/N value was estimated to be 3.4σ and is consistent with the result from bootstrap.

5. Systematic errors

We performed one null test by displacing the galaxy clusters at random positions on the Planck maps in Sect. 4.4. In this section, we describe two additional null tests that we performed to estimate potential systematic errors in our measurements. Subsequently, we checked the contamination in our measured kSZ signals due to the CMB, tSZ, and CIB, respectively.

5.1. Additional null tests

For the first of the two additional null tests, we randomly shuffled the LOS velocities of the galaxy clusters, and then the clusters were stacked with weights based on the shuffled LOS velocities. This was done to test the effect of correlation between LOS velocities and clusters. One shuffling may not be enough to remove the correlation, and therefore we repeated the random velocity shuffling 1000 times and evaluated the mean and standard deviation of the 1000 stacked profiles with shuffled velocities. The second test consisted in performing the stacking with a noise map produced by $(T_{217}^{\mathrm{HM}1}-T_{217}^{\mathrm{HM}2})/2$, where $T_{217}^{\mathrm{HM}1(2)}$ is the half mission 1(2) Planck map at 217 GHz. The associated results are shown in the middle (yellow) and right (green) panels of Fig. 6 with the uncertainties. As expected, the obtained profiles are both consistent with zero. Including the null test performed in Sect. 4.4, all three null tests suggest that our measurements of the kSZ effect are unbiased.

5.2. Contamination from Galactic and extragalactic emissions

The kSZ signal is only a subdominant component and the measurement may be contaminated by Galactic as well as extragalactic emissions such as CMB, tSZ, and CIB. The CMB, the Galactic emissions, and instrumental noise are uncorrelated with cluster positions and are added as noise in our measurements. The tSZ and CIB are correlated with clusters, but it can be assumed that they are uncorrelated with their LOS velocities and canceled out in our analysis. No bias is expected that can be attributed to other components. For a further check, we compared the kSZ signal extracted with the Planck 217 GHz map to the kSZ signals extracted with the Planck CMB maps produced by different component-separation methods such as NILC, SMICA, SMICA-noSZ, SEVEM, and COMMANDER. We note that the Galactic and extragalactic components are not cleaned in the Planck 217 GHz map, but are cleaned in the Planck CMB maps. The extracted kSZ signals are all consistent, as shown in Fig. 5, regardless of the maps we analyze, with or without the Galactic and extragalactic components. This result suggests that the contamination of our measurements from the Galactic and extragalactic components is minor.

5.3. Contamination of the CMB

While we have shown that the contamination of our measured kSZ signals from the Galactic and extragalactic components is minor, contamination from the CMB may still be present at some level. This is because the null tests that we performed involving random displacements and velocity shuffling of galaxy clusters show no statistical bias as 1000 realizations are considered, but the CMB is one realization and the cancelation by stacking may not be sufficient for the real CMB. Therefore, we estimated a possible level of contamination only from the CMB by simulating 100 CMB maps with different realizations using the PYCAMB³ interface to the CAMB⁴ code (Lewis et al. 2000). We repeated our stacking analysis for the 100 simulated CMB maps, producing 100 CMB profiles. The mean of the 100 CMB profiles is consistent with zero and its standard deviation is ∼23% of the amplitude of our measured kSZ signal. Thus, while the CMB contamination may be present at some level, it should not be significant.

5.4. Contamination of the tSZ

The contamination of the tSZ can also be estimated quantitatively by applying our stacking analysis to the Planck all-sky y maps provided in the Planck 2015 data release⁵ (Planck Collaboration XXII 2016). We used the y map from the modified internal linear combination algorithm (MILCA; Hurier et al. 2013), but the result is consistent using the y map from needlet independent linear combination (NILC; Remazeilles et al. 2013). The Compton y parameter in the Planck HFI 217 GHz map can be calculated with its frequency dependence, given by

$$\begin{array}{c}\hfill \frac{\mathrm{\Delta }{T}_{\mathrm{tsz}}}{T_{\mathrm{CMB}}}=g(\nu )y,\end{array}$$(7)

where g(ν) = xcoth(x/2) − 4 with x = hν/(k_B T_CMB), where h is the Planck constant, k_B is the Boltzmann constant, and T_CMB is the temperature of the CMB. At 217 GHz, T_CMB g(ν) = 0.187 K_CMB is given in Planck Collaboration XXII (2016) for the conversion from the Compton y parameter to CMB temperature. The result of stacking the y map shows that the tSZ contamination is ∼1% of the amplitude of our measured kSZ signal and is thus negligible.

5.5. Contamination of the CIB

Similarly, we estimated the level of contamination from the CIB by stacking the Planck all-sky CIB maps provided in Planck Collaboration Int. XLVIII (2016). However, as the CIB map at 217 GHz was not produced, we scaled the Planck CIB map at 353 GHz to the one at 217 GHz with the power-law spectral index of β ∼ 1.1 (Tucci et al. 2016). The result of stacking the CIB map shows that the CIB contamination is ∼6% of the amplitude of our measured kSZ signal and not significant.

6. Interpretation

The CMB temperature fluctuation caused by the kSZ effect is given by

$$\begin{array}{c}\hfill \frac{\mathrm{\Delta }{T}_{\mathrm{ksz}}}{T_{\mathrm{CMB}}}=-{\sigma }_{\mathrm{T}}{\displaystyle \int }{n}_{\mathrm{e}}\left(\frac{\mathit{v}·\widehat{\mathit{n}}}{c}\right)\mathrm{d}l\simeq -\tau \left(\frac{\mathit{v}·\widehat{\mathit{n}}}{c}\right),\end{array}$$(8)

where σ_T is the Thomson scattering cross section, c is the speed of light, n_e is the electron number density, and $\mathit{v}·\widehat{\mathit{n}}$ represents the peculiar velocity of electrons along the line of sight. In the final transformation, the integral, τ = σ_T∫n_edl, was performed along the line of sight under the approximation that the typical correlation length of LOS velocities (given by $\mathit{v}·\widehat{\mathit{n}}$) is much larger than the density correlation length, and therefore the LOS velocity term can be extracted from the kSZ integral. This is justified by Planck Collaboration Int. XXXVII (2016) who find that the typical correlation length of peculiar velocities is 80−100 h⁻¹ Mpc, well above the typical galaxy correlation length of ∼5 h⁻¹ Mpc.

Physical properties of gas can be estimated by considering a β model (Cavaliere & Fusco-Femiano 1978) for a gas(electron) density profile, given by

$$\begin{array}{c}\hfill n_{\mathrm{e}}(r)=n_{\mathrm{e},0}{(1+\frac{r}{r_{\mathrm{c}}})}^{-3\beta /2},\end{array}$$(9)

where n_e, 0 is the central electron density, r is the cluster radial extension, and r_c is the core radius of the electron distribution. In our study, we used β = 0.86 and r_c = 0.2 × R₅₀₀ from the measurements of the South Pole Telescope clusters (Plagge et al. 2010). We can express the data profile as a geometrical projection of the density profile with n_e(r),

$$\begin{array}{c}\hfill \tau (R)={\sigma }_{\mathrm{T}}{\displaystyle \int }\frac{2r{n}_{\mathrm{e}}(r)}{\sqrt{r^2-R^2}}\mathrm{d}r,\end{array}$$(10)

where R is the tangential distance from a galaxy cluster. (We represent the 3D distance with the lowercase letter r, and the 2D distance on a map with the uppercase letter R.)

We fit this model to the measured kSZ profile using the average of the LOS velocities (in absolute value) of the WHL galaxy clusters estimated in Sect. 3. However, the average velocity is overestimated because of the uncertainties of the LOS velocities. (The Gaussian distribution of the estimated LOS velocities have a larger standard deviation than the distribution of the true LOS velocities because of the uncertainties.) The average velocity can be corrected with the uncertainty of the LOS velocities investigated in Sect. 3. The uncertainty on the LOS velocities also induces the decrease in the amplitude of the measured kSZ signal (Nguyen et al. 2020) because it causes some cancelation of the kSZ signal. This can be corrected analytically with the uncertainty on the LOS velocities estimated in Sect. 3. Including these corrections in the model, we fit this β model to the data. The result of the model fitting is shown in Fig. 7. The reduced χ² value is 1.3. We note that we can see a coherent angular pattern in the model profile similar to the data profile, which is due to the convolution of the β profile with our filter in Fig. 2. Defining the optical depth of a galaxy cluster within R₅₀₀ as

$$\begin{array}{c}\hfill {\tau }_{\mathrm{e},500}={\displaystyle {\int }_0^{R_{500}}{\sigma }_{\mathrm{T}}{n}_{\mathrm{e}}(r)\mathrm{d}V,}\end{array}$$(11)

the fitting result provides the average optical depth of the WHL clusters:

$$\begin{array}{c}\hfill {\overline{\tau }}_{\mathrm{e},500}=(3.0\pm 0.9)\times {10}^{-3}.\end{array}$$(12)

Offsets between cluster centers from optical data and centers of gas distribution in clusters may have an impact on the estimated optical depth. Rozo & Rykoff (2014) studied the offsets between X-ray cluster centers (proxy of gas center) and optical cluster centers and showed that ∼70% of the WHL clusters have offsets of less than ∼0.1 Mpc from X-ray clusters. This latter distance corresponds to ∼0.2 arcmin at the median redshift of our sample (z ∼ 0.44), which is smaller than the angular resolution of the Planck maps (∼5 arcmin). To check the contribution of this offset to the measurement of the optical depth, we randomly displaced the centers of the WHL clusters by distances drawn from a Gaussian distribution with the standard deviation of ∼0.1 Mpc and repeated the stacking. The result shows that the amplitude of our measured kSZ profile decreases only by ∼3% and the effect is minor.

We can also estimate a total gas mass in a galaxy cluster defined as

$$\begin{array}{c}\hfill M_{\mathrm{gas},500}={\displaystyle {\int }_0^{R_{500}}{n}_{\mathrm{e}}(r){\mu }_{\mathrm{e}}{m}_{\mathrm{p}}\mathrm{d}V,}\end{array}$$(13)

where μ_e = 1.148 is the mean molecular weight of electrons (Arnaud et al. 2010), and m_p is the mass of proton. Through our measurement of the optical depth, the average gas mass in the WHL clusters can be estimated to be ${\overline{M}}_{\mathrm{gas},500}\sim 1.2\times {10}^{13}$h⁻¹ M_⊙. This provides a gas mass fraction of f_gas, 500 = M_gas, 500/M₅₀₀ = 0.12 ± 0.04 for our sample with an average mass of M₅₀₀ ∼ 1.0 × 10¹⁴ h⁻¹ M_⊙.

The Magneticum simulation also provides gas masses of the simulated clusters within R₅₀₀, and reveals a gas mass of M_gas, 500 ∼ 1.3 × 10¹³ h⁻¹ M_⊙ for the cluster with M₅₀₀ ∼ 1.0 × 10¹⁴ h⁻¹ M_⊙ (average mass of our cluster sample). This corresponds to a gas mass fraction of f_gas, 500 ∼ 0.13 and is consistent with our result.

7. Discussion

The kSZ effect is sensitive to cosmological parameters such as the growth rate of density perturbations. However, we find that our result is consistent using the WMAP and Planck cosmology and the current level of our kSZ detection does not allow us to constrain cosmological parameters. We therefore used our measurement of the kSZ signal to constrain the average optical depth (or gas mass) of the cluster sample. Our result is compared to other measurements as follows.

Gas masses in groups and clusters of galaxies were studied in X-rays by Gonzalez et al. (2013). The authors combined measurements with XMM-Newton and Chandra observations from Vikhlinin et al. (2006), Sun et al. (2009), and Sanderson et al. (2013), and showed the X-ray gas fraction within R₅₀₀ as a function of M₅₀₀. The relation shows f_gas, 500 ∼ 0.1 at M₅₀₀ ∼ 1.3 × 10¹⁴ M_⊙ (average mass of our sample), which is slightly lower than our result of f_gas, 500 = 0.12 ± 0.04, but is consistent within our measurement uncertainties.

Gas in halos was studied with the kSZ in Planck Collaboration Int. XXXVII (2016). The authors detected the pairwise kSZ signal at the positions of central galaxies from the SDSS DR7 data using the Planck and WMAP maps and measured an average optical depth of $\overline{\tau }=(1.4\pm 0.5)\times {10}^{-4}$. This value is one order of magnitude lower then our value of ${\overline{\tau }}_{\mathrm{e},500}=(3.0\pm 0.9)\times {10}^{-3}$, but this discrepancy can be attributed to the difference in halo mass: the halo mass of their sample is M₅₀₀ ∼ 0.3 × 10¹⁴ M_⊙, which is significantly lower than the average mass of our sample, M₅₀₀ ∼ 1.3 × 10¹⁴ M_⊙.

A similar study to that of Planck Collaboration Int. XXXVII (2016) was performed with the kSZ by Soergel et al. (2016) using clusters with photometric redshifts. The authors detected the pairwise kSZ signal by combining a cluster catalog from the Dark Energy Survey (DES) with the CMB map from the South Pole Telescope Sunyaev-Zel’dovich Survery (SPT-SZ) and measured an average optical depth of $\overline{\tau }=(3.75\pm 0.89)\times {10}^{-3}$ for the clusters with M₅₀₀ ∼ (1 − 3)×10¹⁴ M_⊙. This value is higher than our value of ${\overline{\tau }}_{\mathrm{e},500}=(3.0\pm 0.9)\times {10}^{-3}$, but this discrepancy may again be related to the difference in mass: the average mass of our cluster sample, M₅₀₀ ∼ 1.3 × 10¹⁴ M_⊙, corresponds to the lower end of their cluster mass range. Soergel et al. (2016) also estimated the gas mass fraction within R₅₀₀ to be f_gas, 500 = 0.08 ± 0.02 assuming a β model for the density distribution. Our value of f_gas, 500 = 0.12 ± 0.04 is slightly higher than their result, but consistent within ∼1σ when we consider the uncertainties on both measurements.

Lim et al. (2020) also detected the kSZ signals from galaxy groups and clusters in the mass range of 2 × 10¹² M_⊙ < M₅₀₀ < 2 × 10¹⁴ M_⊙ by combining the cluster catalog from Yang et al. (2007) and the Planck frequency maps at 100, 143, and 217 GHz. Surprisingly, their results show that the gas fraction in halos is consistent with the universal baryon fraction in their cluster mass range (f_gas ∼ 0.17). To compare our results with theirs, we converted our measurements to their estimator of ${\stackrel{\sim }{K}}_{500}$, which is the intrinsic kSZ signal scaled to redshift z = 0 and to a fixed angular diameter distance. Our value translates to ${\stackrel{\sim }{K}}_{500}=(1.4\pm 0.4)\times {10}^{-2}$ at M₅₀₀ ∼ 1.3 × 10¹⁴ M_⊙, and the value is consistent with their value within ∼1σ.

Lim et al. (2020) claim that the gas fraction in halos is consistent with the universal baryon fraction down to the mass of 2 × 10¹² M_⊙, but this result is not consistent with those of Soergel et al. (2016) or the X-ray measurements by Gonzalez et al. (2013), showing a lower gas fraction than the universal value. One possible explanation for this discrepancy may be because of the filter used in Lim et al. (2020) to extract the signals. Indeed, Lim et al. (2020) used a matched filter to extract the kSZ signal, which may induce a bias, if the assumed profile is incorrect (Ferraro & Hensley 2015). However, given the large Planck beam, they did not find a significant bias in their results by testing with an incorrect profile. In addition, while they assumed the β profile from Plagge et al. (2010), we used the same β profile and find that it closely fits our measured kSZ profiles with the reduced χ² of 1.3. Therefore, the assumed profile does not explain the difference. On the other hand, our measured kSZ profiles show extended kSZ signal beyond R₅₀₀ and this signal does not seem to be consistent with the result of Lim et al. (2020) that the gas fraction within R₅₀₀ of halos is approximately equal to the universal baryon fraction. However, the beam of the Planck maps (∼5 arcmin) is equivalent to the average angular size of our cluster sample (∼4 arcmin) and this coarse beam prevents a definitive conclusion. Another possible explanation for this discrepancy may be due to the difference in cluster redshifts. Lim et al. (2020) studied a cluster sample at z <  0.12, while the sample in Soergel et al. (2016) is at z ∼ 0.5. Thus, the evolution of the gas in halos may explain the difference in the kSZ measurements. However, the X-ray measurement in Gonzalez et al. (2013) was also applied to local clusters at z ∼ 0.1 and the evolution does not seem to explain the difference. So far, the reason for the discrepancy is unknown. Our result is statistically consistent with Lim et al. (2020) as well as Soergel et al. (2016) and Gonzalez et al. (2013), and more precise measurements of gas are needed in order to come to firm conclusions.

8. Conclusion

In this paper, we present the first direct detection of the kSZ signal with a significance of 3.5σ. The measurement was performed by stacking the Planck temperature map at 217 GHz from the Planck 2018 data release at the positions of Wen-Han-Liu (WHL) galaxy clusters constructed from the SDSS galaxies. If all the clusters are simply stacked, the kSZ signals are canceled out because of the equal probability of clusters showing positive or negative kSZ signals. To avoid this cancelation, we estimated the peculiar velocities of the galaxy clusters along the LOS through the galaxy density field computed from the SDSS galaxies, which is related to the matter density field with the galaxy bias of b = 2 estimated in other studies. Using the LOS velocities as a weight in the stacking, the positive and negative kSZ signals are turned to the same sign and are added together, while other components are canceled out by the positive or negative LOS velocity with equal probability. The measured kSZ signals show a clear correlation with the amplitude of the LOS velocities. As a result of the stacking, we obtain an average kSZ profile for the galaxy clusters with a mass of M₅₀₀ ∼ 1.0 × 10¹⁴ h⁻¹ M_⊙, showing an extended distribution of gas around the galaxy clusters beyond 3 × R₅₀₀.

The kSZ signal is a subdominant component and our measurement may be contaminated by Galactic and extragalactic emission. The CMB, the Galactic emission, and instrumental noise are uncorrelated with cluster positions and are added as noise in our measurements. The tSZ and CIB are correlated with clusters, but uncorrelated with their LOS velocities. We performed three different null tests, the results of which are consistent with zero, indicating that our measurements are unbiased. Possible contamination from the CMB, tSZ, and CIB was also investigated but we did not find any significant contamination from these components.

Based on our kSZ measurement, we estimated the average optical depth and found ${\overline{\tau }}_{\mathrm{e},500}=(3.0\pm 0.9)\times {10}^{-3}$ for clusters with the mass of M₅₀₀ ∼ 1.0 × 10¹⁴ h⁻¹ M_⊙ assuming a β model. This provides an average gas mass of the galaxy clusters of ${\overline{M}}_{\mathrm{gas},500}\sim 1.2\times {10}^{13}$ h⁻¹ M_⊙, leading to a gas fraction of f_gas, 500 = 0.12 ± 0.04 within R₅₀₀. We compared our measurement of gas fraction with that calculated from the Magneticum hydrodynamic simulations as well as the measurements from X-rays (Gonzalez et al. 2013) and kSZ (Soergel et al. 2016), which show that the gas mass fraction is lower than the cosmic baryon fraction. We also compared our result with that of Lim et al. (2020), who, on the other hand, claim that the gas fraction in halos is approximately equal to the universal baryon fraction down to the mass of 2 × 10¹² M_⊙, which is in disagreement with the measurements by Gonzalez et al. (2013) and Soergel et al. (2016). Our result is statistically consistent with the results of these authors, and more precise measurements of gas are necessary in order to make firm conclusions.

Acknowledgments

The authors thank an anonymous referee for the useful comments and suggestions. This research has been supported by the funding for the ByoPiC project from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme grant agreement ERC-2015-AdG 695561. The authors acknowledge fruitful discussions with the members of the ByoPiC project (https://byopic.eu/team). Furthermore, S.Z. acknowledge support by the Israel Science Foundation (grant no. 255/18). We thank LSS2LSS science program of the Ψ2 initiative of the Institut Pascal, funded by the Paris-Saclay University, in which this research has been initiated. This publication used observations obtained with Planck (http://www.esa.int/Planck), an ESA science mission with instruments and contributions directly funded by ESA Member States, NASA, and Canada. The authors thank Klaus Dolag and Antonio Ragagnin for providing the Magneticum simulations, and Adi Nusser for fruitfull discussions.

References

Appendix A: Components at HFI 217 GHz

The stacking of the Planck temperature map at 217 GHz shows an excess at the positions of the WHL clusters, as seen for example in the black line in the upper left panel of Fig. 3. We investigated the contribution to this excess from the Galactic and extragalactic components.

In the Galactic emissions, we considered the dust and CO emission. The other Galactic emissions such as synchrotron are minor at 217 GHz according to Planck Collaboration X (2016). The all-sky Galactic dust map was provided in Planck Collaboration Int. XLVIII (2016). However, as the dust map at 217 GHz was not produced, we scaled the Planck dust map at 353 GHz to the one at 217 GHz with the power-law spectral index of β ∼ 1.5 (Planck Collaboration X 2016). The CO emission map at 230 GHz, because of the CO rotational transition of J = 2 − 1, was produced in Planck Collaboration X (2016). In the extragalactic emissions, we consider the tSZ and CIB component and use the tSZ and CIB maps scaled at 217 GHz, as described in Sects. 5.4 and 5.5.

We stacked these maps at the positions of the WHL clusters and compare the stacked signals in Fig. A.1. We note that all these maps were convolved by an additional beam, giving rise to a common beam size of 10 arcmin in FWHM, which is the largest beam of the Planck tSZ map. The dominant contribution at small scales is from the CIB and that from the rest of the components is minor. The signals at large scales cannot be well explained by the components we consider. This may be due to an unknown spectral index of CIB or relatively uncertain CO emission. However, the signals at large scales do not contribute to our kSZ signals at small scales (< 4 × R₅₀₀) and we do not take the difference into account in this paper.

Fig. A.1. Radial profile around the 30 431 WHL clusters with the Planck temperature map at 217 GHz (black) with contributions from CMB (blue), CIB (magenta), the Galactic dust (cyan), tSZ (green), and CO emissions (yellow) as well as the total (red).

Appendix B: Contamination to the kSZ signal

In order to evaluate the contamination of the measured kSZ signal attributable to the Galactic and extragalactic emissions, we used the maps of CMB, CIB, tSZ, and Galactic dust convolved to a common beam size of 10 arcmin in FWHM, as described in Appendix A. We repeated the same stacking analysis in Sect. 4 for these maps as for the HFI map at 217 GHz. The results are shown in Fig. B.1 and the contribution from each component is summarized in Sect. 5.

Fig. B.1. Contamination in our measured kSZ profile (black) due to contributions from CMB (blue), CIB (magenta), the Galactic dust (cyan), and tSZ (green).

All Figures

	Fig. 1. x-axis: true LOS velocities of simulated galaxy clusters in the Magneticum simulations. y-axis: LOS velocities of the same galaxy clusters computed using Eq. (1) with the galaxy bias of b = 2 when the RSD effect is not included (left) and included (right).
In the text

	Fig. 2. Filter applied to the Planck maps. The filter is a smooth function with its response of one below 15 arcmin (ℓ ∼ 720) and zero above 30 arcmin (ℓ ∼ 360). With this filter, the primordial CMB fluctuations in black are suppressed to the ones in red.
In the text

Fig. 3. Upper panels, left: stacked radial profile around the 30 431 WHL clusters with the Planck temperature map at 217 GHz (black) with the 1σ uncertainty estimated by a bootstrap resampling (gray). The stacking is performed separately for the 14 093 WHL clusters with negative LOS motions (black dashed line) and for the 16 338 WHL clusters with positive motions (black dash-dotted line). Middle: same as left panel but with an additional LOS velocity cut of 100 km s−1. Right: same as left panel but with an additional LOS velocity cut of 200 km s−1. Lower panels: positive kSZ radial profile (blue) extracted by taking the difference between the black dashed and black lines, and negative kSZ radial profile (red) extracted by taking the difference between the black and black dash-dotted lines.

In the text

Fig. 4. Positive kSZ radial profile (blue) around the 14 093 WHL clusters with the Planck temperature map at 217 GHz with the 1σ uncertainty estimated by a bootstrap resampling, and negative kSZ radial profile (red) around the 16 338 WHL clusters, compared to positive kSZ radial profiles (green dashed) and negative kSZ radial profile (green dash-dotted) from the Magneticum hydrodynamic simulations. Middle: same as left panel with an additional LOS velocity cut of 100 km s−1. Right: same as left panel with an additional LOS velocity cut of 200 km s−1.

In the text

	Fig. 5. Velocity-weighted kSZ radial profile around the 30 431 WHL clusters with the Planck temperature map at 217 GHz (black) and with the 1σ uncertainty estimated by a bootstrap resampling, compared to the kSZ radial profiles with the Planck CMB maps produced by different component separations of NILC (green), SMICA (yellow), SMICA-noSZ (purple), SEVEM (cyan), and COMMANDER (orange).
In the text

Fig. 6. Velocity-weighted kSZ radial profile around the 30 431 WHL clusters with the Planck temperature map at 217 GHz (black) compared to three null tests. Left panel: clusters are displaced at random positions on the Planck map, and then stacked. This process is repeated 1000 times and the mean of the 1000 random samples is computed (cyan). The 1σ uncertainty is estimated by computing a standard deviation of the 1000 random samples. Middle panel: results from cluster stacking after randomly shuffling the LOS velocities of the clusters; this process is repeated 1000 times and the mean of the 100 velocity-shuffled profiles is computed (yellow). The 1σ uncertainty is estimated by computing a standard deviation of the 1000 velocity-shuffled profiles. Right panel: clusters are stacked with a noise map produced by $(T_{217}^{\mathrm{HM}1}-T_{217}^{\mathrm{HM}2})/2$, where $T_{217}^{\mathrm{HM}1(2)}$ is the half mission 1(2) Planck map at 217 GHz (green). The 1σ uncertainty is estimated by a bootstrap resampling.

In the text

	Fig. 7. Velocity-weighted kSZ radial profile around the 30 431 WHL clusters with the Planck temperature map at 217 GHz (black), fitted with the β model (red).
In the text

	Fig. A.1. Radial profile around the 30 431 WHL clusters with the Planck temperature map at 217 GHz (black) with contributions from CMB (blue), CIB (magenta), the Galactic dust (cyan), tSZ (green), and CO emissions (yellow) as well as the total (red).
In the text

	Fig. B.1. Contamination in our measured kSZ profile (black) due to contributions from CMB (blue), CIB (magenta), the Galactic dust (cyan), and tSZ (green).
In the text