Detection of Pairwise Kinetic Sunyaev–Zel’dovich Effect with DESI Galaxy Groups and Planck in Fourier Space

We report a ∼5.2σ detection of the kinetic Sunyaev–Zel’dovich (kSZ) effect in Fourier space, by combining the DESI galaxy groups and the Planck data. We use the density-weighted pairwise kSZ power spectrum as the summary statistic, and the detailed procedure of its measurement is presented in this paper. Meanwhile, we analyze the redshift space group density power spectrum to constrain its bias parameters and photo-z uncertainties. These best-fitted parameters are substituted to a nonlinear kSZ model, and we fit the measured kSZ power spectrum with this model to constrain the group optical depth τ¯ . Selected by a varying lower mass threshold M th, the galaxy group catalogs with different median masses ( M˜ ) are constructed from the Data Release 9 data of the DESI Legacy Imaging Surveys. M˜ spans a wide range of ∼1013–1014 M ⊙/h and the heaviest M˜∼1014 M⊙/h is larger than those of most other kSZ detections. When the aperture photometric filter radius θ AP is set to be 4.′2, the M˜=1.75×1013 M⊙/h group sample at the median redshift z˜=0.64 has the highest kSZ detection signal-to-noise ratio = 5.2. By fitting τ¯s from various samples against their M˜ s, we obtain a linear logτ¯−logM˜ relation: logτ¯=γ(logM˜−14)+logβ , in which γ = 0.55 ± 0.1. We also vary the aperture photometric filter radius and measure the τ¯ profiles of group samples, whose constraints on the baryon distribution within and around dark matter halos will be discussed in a companion paper.


INTRODUCTION
The kinetic Sunyaev-Zel'dovich (kSZ) effect is a secondary cosmic microwave background (CMB) effect induced by the scattering of the CMB photon off free electrons with a bulk motion (Sunyaev & Zeldovich 1970, 1972, 1980;Phillips 1995;Birkinshaw 1999).Its detection enables detailed studies of the baryon distribution in the universe and the cosmological model.On one hand, the kSZ signal is linearly proportional to the free electron number density and is independent of the gas temperature.It can be used to constrain the Epoch-of-Reionization (EoR) history (Alvarez 2016;Chen et al. 2023b) or the gas distribution within halos (Schaan et al. 2016;Sugiyama et al. 2018a;Chaves-Montero et al. 2021;Schneider et al. 2022) and filaments (Zheng et al. 2023), and is well suited to search for missing baryons which are believed to mainly reside in low density and low temperature environments.On the other hand, the kSZ effect captures the peculiar velocity field of our universe and records its structure growth information.It can help constrain different aspects of the cosmological model such as Copernican principle (Zhang & Stebbins 2011), primordial non-Gaussianity (Münchmeyer et al. 2019;Kumar et al. 2022), dark energy (Pen & Zhang 2014;Okumura & Taruya 2022), modified gravity (Bianchini & Silvestri 2016;Roncarelli et al. 2018;Zheng 2020;Mitchell et al. 2021), massive neutrinos (Roncarelli et al. 2017) and so on.
Albeit, being rich in cosmological information, the detection of the kSZ effect is difficult in its nature.The amplitude of the kSZ signal is 2 orders of magnitudes lower than the CMB primary fluctuation and several times lower than the thermal Sunyaev-Zel'dovich (tSZ) signal (Zhang et al. 2004).Besides, its spectrum is the same as that of the primary CMB; thus, multifrequency measurements do not help on its detection, in contrast to the tSZ effect.As a result, the current detection signal-to-noise ratio (S/N) of the separated kSZ component from the CMB power spectrum is still low (Dunkley et al. 2013;George et al. 2015;Bleem et al. 2022).Fortunately, the galaxy resides in high electron density environment and its celestial position denotes where the kSZ signal is relatively high.By cross-correlating the CMB data with the galaxy redshift survey data, we up-weight the high kSZ signal regions and enhance the measurement S/N (Hand et al. 2012).Adopting this so-called kSZ tomography technique, we have detected the late-time (post Epoch-of-Reionization) kSZ effect in a significance around 4σ ∼ 7σ, with combinations of different CMB data and galaxy catalogs, e.g., South Pole Telescope (SPT)+Dark Energy Survey (DES) (Soergel et al. 2016;Schiappucci et al. 2023), Atacama Cosmology Telescope (ACT)+SDSS (Schaan et al. 2021;Calafut et al. 2021), Planck/Wilkinson Microwave Anisotropy Probe+unWISE (Kusiak et al. 2021), Planck+DESI(photo-z) (Chen et al. 2022), et al.With the next generation of surveys, S/Ns of such measurements can reach O(100) and beyond (Sugiyama et al. 2018a;Smith et al. 2018).
The kSZ tomography technique can by implemented by various estimators.Most of them are mathematically equivalent to the three-point correlation function (or bispectrum) of the type ⟨ggT ⟩ or ⟨gT T ⟩ (Smith et al. 2018).Here, "g" is the power of the galaxy field, and "T " is the power of the CMB field.The majority of such estimators are ⟨ggT ⟩-like, such as (1) the pairwise estimator, in which we sum the pair difference of the high-pass filtered CMB at galaxy positions within a certain pair separation bin (Hand et al. 2012;Sugiyama et al. 2018a;Gong et al. 2023); (2) the kSZ template method, in which the kSZ template is cross-correlated with the CMB map to extract the kSZ signal, and the core of the kSZ template is the projected galaxy momentum field which is reconstructed from the observed galaxy density field (Ho et al. 2009;Shao et al. 2011); (3) the stacked velocity matched filter method, in which we stack all the high-pass filtered CMB at galaxy locations weighted by the galaxies' velocities to detect the kSZ signal, and the galaxy velocity filed is reconstructed from the galaxy density field (Li et al. 2014;Schaan et al. 2016Schaan et al. , 2021)); (4) the growth reconstructed method, in which we directly compare the velocity field constructed from the galaxy density field with the high-pass filtered CMB data at galaxy locations with the maximum likelihood estimation method and extract the structure growth information (Alonso et al. 2016); (5) the velocity reconstruction method, in which we directly reconstruct the large-scale radial velocity field from the galaxy density-weighted CMB data with a quadratic estimator (Deutsch et al. 2018) or via the machine-learning technique (Wang et al. 2021).The ⟨gT T ⟩-like estimator includes the one by cross-correlating the large-scale structure with the squared high-pass filtered CMB data (DeDeo et al. 2005;Hill et al. 2016;Ferraro et al. 2016;Patki et al. 2023).Besides these three-point ones, higher order estimators for the kSZ tomography with the ⟨gggT ⟩-type has also been proposed (Hernández-Monteagudo et al. 2020;Chaves-Montero et al. 2021;Kuruvilla 2022).
Within all these methods, the pairwise estimator was used to extract the first cosmological kSZ signal (Hand et al. 2012), and is most widely adopted in a real data analysis (Hand et al. 2012;Hernández-Monteagudo et al. 2015;Planck Collaboration et al. 2016;Soergel et al. 2016;De Bernardis et al. 2017;Sugiyama et al. 2018a;Calafut et al. 2021;Chen et al. 2022).We will also implement this estimator in this work.
A dual analysis both in configuration and Fourier space has been a common practice in the galaxy clustering study.Although the correlation function and the power spectrum are Fourier counterparts and contain the same cosmological information in an ideal case, the cosmic variance on large-scales, the shot noise on small scales, and different observational and theoretical systematics restrict the available data within a certain scale range and render different pros and cons to these two statistics (Verde 2007).For example, the correlation function is straightforward to detect observationally while its covariance matrix estimation and redshift space modeling is nontrivial.On the other hand, the power spectrum has a more diagonal covariance matrix and its anisotropic clustering property is easier to model in redshift space, yet its measurement is contaminated by more systematics such as the anisotropic survey window function and the complicated selection function.In this sense, the analysis in both spaces are complementary to each other and are both valuable in constraining cosmological models.Following this spirit, we implement a Fourier space kSZ analysis in this work, while a complementary configuration space kSZ detection from the similar data set has been published in Chen et al. (2022).
The first and only Fourier space kSZ detection until now is from Sugiyama et al. (2018a).We follow the main routine developed in that paper and make several modifications/improvements, including the following: (1) the median group masses ( M ) of our samples range from 10 13 − 10 14 M ⊙ /h, in which the heaviest M ∼ 10 14 M ⊙ /h is larger than those of most other works, and by this we explore a new range of halo mass in the kSZ detection (see also Soergel et al. (2016)); (2) we implement the aperture photometry (AP) filter in the spherical harmonic space (Chen et al. 2022), which can avoid the ambiguity of the pixelized temperature average on the circular filter in configuration space; (3) when applying the kSZ power spectrum estimator including the moving line-of-sight (LOS) effect, we decompose the Legendre polynomial into into a product of spherical harmonics (Hand et al. 2017;Sugiyama et al. 2018b), instead of a Cartesian decomposition (Bianchi et al. 2015;Scoccimarro 2015), and this will reduce the number of needed fast Fouirer transforms (FFTs) if we include the multipole calculation higher than the dipole in the future; ( 4 The paper is organized as follows.In Section 2, we introduce the CMB map and the galaxy group/cluster catalogs we will analyze in this work.In section 3, a kSZ power spectrum model taking into account the photo-z error is developed.In Section 4, we outline the methodology we imply to measure the kSZ power spectrum, and in Section 5, we display the main results of this work.We conclude in Section 6 and more technique details and complementary results are shown in appendices.We will adopt the Planck2018 cosmology (Planck Collaboration et al. 2020) throughout the work, and calculations of cosmological quantities are made via the Nbodykit Python toolkit 1

DATA
The detection of the pairwise kSZ effect requires the synergy between a CMB map and a galaxy (group/cluster) catalog.The angular position of a galaxy or galaxy group indicates a sky location likely with a high kSZ signal.A compensated AP filter is applied to the CMB data on this location, filtering out the contamination from the primordial CMB fluctuation and the detector noise, and extracting the desired kSZ signal.All these individual kSZ signals are then fed into the pairwise kSZ estimator (Equations ( 4) and ( 5)) and the output is the statistical pairwise kSZ effect.This pairwise property of the estimator tackles the directional dependence of the kSZ signal since the result simply vanishes in a trivial average of signals from all locations due to the isotropic peculiar velocity field.In this section we first introduce the data sets that we are going to analyze in this work.

Planck data
1 https://github.com/bccp/nbodykitfrom Hand et al. (2018).The kSZ signal is measured on a full-sky intensity map, High Frequency Instrument (HFI) 217 GHz2 of the public Planck Release 3 data3 .This map has an effective FWHM ≈ 4.87 arcmin and is provided in the HEALPix grid frame (Górski et al. 2005) with N side = 2048.We update the map to N side = 8096 with the zeroth-order interpolation in order to improve the map resolution when applying the AP filter (Section 4.1).
The contamination of the tSZ effect is largely alleviated due to its vanishment in the frequency of 217 GHz.Although its residual is nonzero due to the finite bandwidth, it is not expected to have the directional dependence as the pairwise kSZ signal does, thus will be further eliminated in our estimator and does not bias the signal.The same cleaning procedure happens for different foregrounds on the CMB map when we apply the pairwise estimator.Compared to the CMB foreground cleaned maps such as SMICA, SEVEM, NILC and COMMANDER, the 217 GHz map has higher angular resolution, which is highly desired in the kSZ detection.By making comparison between different maps in Appendix C of Chen et al. (2022), the 217 GHz map is shown to be optimal in our kSZ measurement.

DESI photometric Data Release 9 galaxy groups
We use the galaxy group catalog4 provided by Yang et al. (2021) to identify locations where kSZ signals are large.This catalog is constructed from the Data Release 9 (DR9) of the the DESI Legacy Imaging Survey (Dey et al. 2019), by extending the halo-based group finder developed in Yang et al. (2005).After removing the sky area within |b| < 25. • 0 and around nearby sources or masked pixels, the sky coverage of the catalog remains 9622 deg 2 in the North galactic cap (NGC) and 8601 deg 2 in the South galactic cap (SGC).In total, the group catalog contains 51.44 million galaxy groups in NGC and 45.41 million galaxy groups in SGC with 3D coordinate, richness, halo mass and the total group luminosity.
The halo gas mass is expected to increase with its total mass (e.g., Chen et al. (2022)); therefore, we construct group catalogs by varying a lower group mass threshold M th .From this, we construct group samples with median masses ( M ) ranging from 10 13 − 10 14 M ⊙ /h.In Appendix D, we show how the detection S/N varies with M th , and we recognize the group sample with the highest detection S/N as the baseline sample in this work.
The baseline sample we construct contains ∼ 5.65 × 10 6 heaviest groups in the group catalog, spanning between 0 < z < 0.925 .Their pixelized sky distribution is shown on the left panel of figure 1, and their mass, redshift, angular radius and richness distributions are shown in figure 2. We mark the median values of {M, z, θ halo , richness} for the NGC and SGC subbaseline samples in figure 2. The {M, z} medians will be used in our following calculation of the theoretical kSZ power spectrum.In Section 5, we will present the measured kSZ power spectum and associated results from the baseline sample.
On the lower left panel of figure 2, half of the groups in the baseline sample have an angular radius lower than 1.5 ′ , much smaller than the Planck FWHM, which results in a heavy smoothing effect on the kSZ signal and limits its S/N.In this sense, the Planck data are not ideal for the SZ effect detection, although its all-sky coverage compensates this disadvantage.On the contrary, high resolution CMB experiments such as ACT and SPT are limited by their relatively small sky coverage overlaps with current galaxy surveys (Calafut et al. 2021;Bleem et al. 2022).On the lower right panel of figure 2, half of the groups contain ≤ 3 group members, and most of them reside at high redshifts.The typical photo-z error of the DESI photometric galaxies can be parameterized by σ pho−z = 0.01 + 0.015z (Zhou et al. 2021;Wang et al. 2020;Yang et al. 2021).For groups with richness ≥ 3, Yang et al. (2021) identified that their σ pho−z ranges from 0.005 at z = 0 to 0.015 at z = 1.These priors indicate that our group samples will have typical photo-z errors of σ pho−z = 0.005 ∼ 0.025.Such large redshift uncertainty will heavily dilute group clustering along the LOS.Moreover, regarding that our selected samples span large redshift ranges, the complicated photo-z probability distribution function (PDF), which is redshift dependent and possibly asymmetric, makes the modeling of the galaxy group kSZ and density power spectra in redshift space a challenging task.In Appendix A we test the sensitivity of the optical depth measurement on different photo-z damping models and make sure our results are not significantly biased by the photo-z uncertainty.

THEORY
In this section we present the theoretical model that we will use to fit the measured density-weighted pairwise kSZ power spectra in Section 5, which are generated by the methodology described in Section 4.

kSZ basics
The secondary anisotropy of the CMB temperature caused by the scattering off of CMB photons by free electrons with bulk motions is known as the kSZ effect (Sunyaev & Zeldovich 1972, 1980).The change of the CMB temperature is Here, n e is the physical free electron number density, v • n is the electron peculiar velocity projected along the LOS direction, T CMB = 2.7255K is the average CMB temperature, σ T is the Thomson scattering cross section, and c is the speed of light.The integral dl in this expression is performed along the LOS.
If the CMB photons coming from the n direction are only scattered once by the ith galaxy group, the kSZ effect generated by the cloud of free electrons around this group can be simplified as where v i is the physical bulk motion of the ith group, and τ i = n e σ T dl is its optical depth, depending on its baryon abundance.Assuming further that all groups share approximately the same optical depth, we have in which τ is the average optical depth of the group sample.

Density-weighted pairwise kSZ power spectrum
The trivial average of the kSZ signal in Equation ( 3) vanishes due to its directionl dependence.Fortunately, under the influence of gravity, there is a tendency of two groups moving toward each other on quasi to nonlinear scales.In turn, we can detect the kSZ effect by a pairwise estimator (Sugiyama et al. 2018a): where ξ kSZ is the density-weighted pairwise kSZ correlation function, s ij = s i − s j where s i denotes the redshift space position of the ith group, N represents the number of observed groups, and V is the volume of the survey.
The Fourier counterpart of ξ kSZ is the density-weighted pairwise kSZ power spectrum, Inserting Equation (3) into Equations ( 4) and ( 5), ξ kSZ and P kSZ are proportional to the density-weighted pairwise velocity correlation function ξ pv and power spectrum P pv respectively: where A hidden assumption adopted here is that the group optical depth has no correlation with the group peculiar velocity.
Compared to the traditional pair-weighted statistics, e.g., ξ pair−weight pv = ξ pv /(1 + ξ δ ), the density-weighted statistics encode similar cosmological information, while their modeling is simpler by avoiding the modeling of an additional density correlation function ξ δ in the denominator.A detailed comparison between the density-and pair-weighted statistics can be found in Appendix B of Okumura et al. (2014).
Under the plane-parallel approximation and the condition that the peculiar velocity field is proportional to the linear growth rate f ≡ d ln D/d ln a (which is generally true and D is the linear growth factor), P pv in redshift space can be derived from the redshift space galaxy power spectrum P s (k) by (Sugiyama et al. 2016(Sugiyama et al. , 2017) where a is the cosmic scale factor, H is the Hubble parameter.
In this work, we use the following nonlinear galaxy power spectrum to derive the kSZ power spectrum: In Appendix A, this model is adopted to analyze the redshift space group density power spectrum.Here, the superscript 'spec-z' suggests that we are analyzing a galaxy sample with negligible redshift uncertainties, e.g., from a spectroscopic galaxy survey, b 1 and b 2 are linear and nonlinear galaxy density bias parameters, µ ≡ cos θ denotes the cosine of the angle θ between k and the LOS, P m (k) is the nonlinear dark matter power spectrum at redshift z, a Lorentzian Finger-of-God (FoG) term is adopted and the matter velocity dispersion σ 2 v is calculated by σ 2 v = f 2 H 2 P m dk/6π 26 , assuming a unity volume-weighted galaxy velocity bias (Zheng et al. 2015;Chen et al. 2018).
In order to tackle the nonlinear correlation between density and velocity fields, we include the window function W (k) ≡ P δθ /f P δδ in Equation ( 11) (Zhang et al. 2013).Motivated by the perturbation theory and verified with N-body simulations, W (k) can be described by the following fitting function (Zheng et al. 2013) where ∆ 2 δδ (k, z) = k 3 P δδ /2π 2 is the dimensionless density power spectrum, and the redshift dependence of ∆α can be fitted well by a linear relation (Zheng et al. 2013) ∆α(z) = ∆α(z = 0) + 0.082z = 0.331 + 0.082z . ( By inserting Equation ( 11) into Equations ( 7) and ( 10), we derive the nonlinear density-weighted pairwise kSZ power spectrum as v /H 2 .At linear scales, the above nonlinear model can be simplified to a linear one It is easy to figure out the degeneracy between b 1 and τ , namely that an accurate τ detection requires an accurate b 1 constraint.In the next section, we further notice that the photo-z uncertainty heavily damps the kSZ power spectrum, and this damping is also degenerated with b 1 and τ .This is a unique feature for the photo-z data.In Appendix A, we tackle these issues by applying a joint analysis of kSZ and RSD analysis on the same data set.The tests there show that our τ measurements are minimally biased by the abovementioned degeneacies.

Density-weighted pairwise kSZ power spectrum with photometric redshits
The photometric data we analyze have nonnegligible redshift errors.These photo-z errors randomly distort the LOS positions of galaxy groups and dilute their LOS clustering, which is an analogy of the FoG effect.As discussed in detail in Appendix A, we model the distribution function of the pairwise photo-z errors as an exponential function with zero mean and a dispersion of σ pair pho−z = √ 2σ pho−z , with σ pho−z being the dispersion of the 1-D photo-z error distribution function; then, the galaxy power spectrum of the photo-z data becomes (Davis & Peebles 1983;Ballinger et al. 1996) where D pho−z is the photo-z damping function and σ z ≡ cσ pho−z .Since σ z does not depend on f , the resultant kSZ power spectrum is similarly diluted by the photo-z errors as We will omit the superscript 'pho − z' hereafter. 6In the spec-z data redshift space distortion (RSD) analysis, σ 2 v is usually treated as a free parameter, which effectively absorbs the uncertainty of the Kaiser term modeling.In this photo-z data analysis, the power spectrum damping from the photo-z error is much larger than that induced by σ 2 v .So when we treat the photo-z error dispersion σ 2 pho−z as a free parameter, a direct calculation of σ 2 v here will not trigger severe systematic errors.Solid lines represent multipoles of the nonlinear density-weighted pairwise kSZ power spectrum with photo-z errors (Equation ( 17)) and dashed lines are those without photo-z errors (Equation ( 14)).The red and blue lines are for dipoles (ℓ = 1) and octopoles (ℓ = 3) respectively.The adopted redshift is the median redshift z = 0.64 of the NGC sub-sample, and τ = 3.96×10 −5 is the best-fitted τ value of the baseline sample.
In Appendix A, we measure the galaxy density power spectrum multipoles of group catalogs and fit the measurements with the nonlinear model in Equations ( 11) and ( 16).With f being fixed to the fiducial cosmological value, this fitting determines the posterior distribution of b 1 , b 2 and σ z .We then substitute the best-fitted values of these 3 parameters to Equation ( 14) and fit the kSZ measurement to obtain the constraint on τ .In this way, (1) we model the galaxy density and kSZ power spectra in a self-consistent manner; (2) we self-consistently break the b 1 − τ degeneracy and avoid possible systematic biases coming from a manually selected b 1 − M relation and σ pho−z .; (3) as τ is the only parameter in the kSZ measurement fitting, the S/N of the fitted τ directly represents the S/N of the kSZ detection, which is one of main targets of this work.For a global fitting of the galaxy density and kSZ power spectra simultaneously, which requires calculations of the covariance matrix between two power spectra, we defer it to the future work.
Moreover, we expand P kSZ (k, µ) in terms of Legendre polynomials, Since P kSZ (k, µ) is an odd function on µ, P ℓ kSZ is nonzero only for odd ℓs, in contrast to the galaxy power spectrum P s (k, µ).As shown by solid lines in figure 3, the absolute values of P ℓ=3 kSZ (k, µ) are smaller than those of P ℓ=1 kSZ (k, µ), and in principle, the measurement errors of the former will be larger.Considering the limited S/N of our measurement, we only focus on the dipole P ℓ=1 kSZ in this work.By dashed lines, we also plot the nonlinear kSZ power spectrum multipoles without photo-z errors evaluated by Equation ( 14) in figure 3.These will be signals likely detected by future spectroscopic galaxy surveys, e.g., DESI (DESI Collaboration et al. 2016) or Subaru Prime Focus Spectrograph (PFS) (Takada et al. 2014).Without the smearing effect of photo-z errors, kSZ multipoles exhibit more measurable features at nonlinear scales, and will guarantee a kSZ detection of higher S/N (Sugiyama et al. 2018a).

METHODOLOGY
In this section, we outline the methodology that we follow to measure the kSZ signal and the associated covariance matrix.

Aperture photometry filter
In order to minimize the contamination from the primary CMB, an AP filter is applied to the Planck CMB map (e.g.Sugiyama et al. (2018a); Calafut et al. (2021); Chen et al. (2022)).The AP filter is a 2D compensated top-hat filter: where the averaged CMB temperature within a disk of aperture size θ AP and an annulus of equal area, out to radius √ 2θ AP , are differentiated around each galaxy group to measure the kSZ signal ∆T AP kSZ : in which T map denotes the CMB temperature map.
In this work, we follow Chen et al. ( 2022) and apply the AP filter in the spherical harmonic space, namely where and J 1 is the first Bessel function of the first kind.

Random catalog
We rely on the random catalog to calculate the galaxy group density field in Equation ( 25).We generate a corresponding random catalog for each group sample we construct, by randomly downsizing the random catalog provided along with the DESI Legacy Imaging Survey DR9 data7 .The size of each random catalog is 20 times of the size of the corresponding group sample.Since the group catalog of Yang et al. (2021) is constructed within the continuous region of the DESI sky coverage, we manually discard the discontinuous regions in the random catalog in the first step.The pixelized sky distribution of the random catalog is shown on the right panel of figure 1.
We apply the shuffled method (Ross et al. 2012), which is proven to be robust in the galaxy clustering analysis, to assign redshifts to the random catalog.In this method, each random point is assigned a redshift which is randomly selected from the redshifts of the group sample.

kSZ and galaxy group density fluctuations
In Section 4.1, the AP filtered CMB temperature ∆T AP kSZ of each group is estimated in the spherical harmonic space.To suppress the systematic errors induced by the redshift-dependent foregrounds such as the cosmic infrared background and the residual tSZ signal, we subtract from ∆T AP kSZ its average with a Gaussian weight (Hand et al. 2012): where G(z i , z j , σ z ) = exp(−z 2 ij /2σ 2 z ) with σ z = 0.01, and z ij ≡ z i − z j .The 3D kSZ temperature fluctuation field then reads where N group is the size of the group sample and δ D is the 3D Dirac delta function.
With the help of the random catalog constructed in Section 4.2, we define the 3D galaxy group density field as where N random denotes the number of random points in the random catalog, and the factor α ≡ N group /N random = 0.05 in this work.

Estimator of the pairwise kSZ power spectrum multipoles
We adopt the estimator of the density-weighted pairwise kSZ power spectrum multipoles developed by Sugiyama et al. (2018a), which is an analogy to the estimator of the galaxy power spectrum multipoles (Feldman et al. 1994;Yamamoto et al. 2005), where s 12 = s 1 − s 2 and the unit vector n12 of n 12 = (s 1 + s 2 )/2 denotes the LOS direction to the mid-point of the pair of objects.These multipoles are normalized by (Feldman et al. 1994) where the number density n is directly measured from the random catalog.We do not use the FKP weight in this estimator since it is not derived for an optimal kSZ detection.
The shot noise term disappears in Equation ( 26), as it is canceled out by the subtraction of the two terms in the square bracket.This is also true for the detector noise and the primary CMB anisotropies contaminating the δT (s) measurement.The pairwise estimator guarantees that they will not bias the kSZ power spectrum, although their contributions to the covariance matrix remain.
Under the local plane-parallel approximation, Equation ( 26) can be simplified as (Yamamoto et al. 2005;Sugiyama et al. 2018a) where The evaluation of the Legendre polynomial involved quantity (δn ℓ here) is time consuming.It is proposed in Bianchi et al. (2015); Scoccimarro (2015) that, by decomposing k • ŝ into its Cartesian components, Equation (30) can be expressed as a sum over the Fourier transforms of δn(s) weighted by products of Cartesian vectors.Thus, the application of FFT technique8 can speed up the calculation in several orders of magnitude as compared to a naive summation implementation of Equation ( 30).The number of FFTs required in each Pℓ calculation is (ℓ + 1)(ℓ + 2)/2 in this scenario.
Instead of the Cartesian decomposition, Hand et al. (2017) and Sugiyama et al. (2018b) proposed to expand the Legendre polynomial into a product of spherical harmonics using the spherical harmonic addition theorem, Then, Equation (30) becomes In this case, the implication of only 2ℓ + 1 times of FFTs is needed to evaluate each δn ℓ .Although when ℓ = 1, both methods require 3 FFTs; the improvement will become evident as ℓ increases.Therefore, we adopt the second method to calculate the pairwise kSZ power spectrum dipole, and the calculation of Equation ( 33) heavily refers to the software toolkit nbodykit (Hand et al. 2018).Finally, we select appropriate k bins, k i ∈ [0.05, 0.155] (h/Mpc) with interval ∆k = 0.01h/Mpc, and average the |k| valves within each bin to obtain an effective k value.This tackles the binning discreteness effect, which can be prominent at linear scales.Being dominated by noises and photo-z damping, larger k modes cannot increase the measurement S/N.Furthermore, we average the measured P ℓ kSZ within each k bin from Equation (29) to compute our estimated P ℓ kSZ (k): where N mode is the number of Fourier modes in a k bin (k i ≤ k < k i+1 ).

Measurement details
Referring to Gil-Marín et al. ( 2016); Sugiyama et al. (2018a), the power spectrum multipoles of NGC and SGC are weighted by their effective areas to obtain the combined power spectrum multipoles: where A NGC and A SGC are effective coverage of the contiguous areas of NGC and SGC, namely A NGC = 9622 deg 2 and A SGC = 8601 deg 2 .Equation ( 36) is used analogically to estimate other parameters, such as the median redshift, median mass of the whole group sample, etc.The FFTs are performed on regular grids with a fiducial size of 512 3 .We establish the Cartesian coordinates s = (x, y, z) with the z-axis pointing to the north pole.Then, we locate the group samples into a cuboid of dimensions (L x , L y , L z ), with (L x , L y , L z )[h −1 Mpc] = (2400, 4000, 2600) for NGC, and (2600, 4000, 3500) for SGC.The averaged grid size is ∼ 6h −1 Mpc for both NGC and SGC.The triangular-shaped cloud (TSC) assignment function is used to distribute both the groups and associated δT s onto the grids.This gridding process induces numerical artifacts like the aliasing effect and the window function effect (e.g., Jing (2005)).In Fourier space, we follow the interlacing technique proposed by Hockney & Eastwood (1981); Sefusatti et al. (2016) to alleviate the aliasing effect, and correct the window function (e.g., Equation (18) of Jing ( 2005)) for each k mode on the grids.As illustrated by figure 2 of Hand et al. (2018), these corrections guarantee to reduce the corresponding artifacts down to a subpercent level even at scales near the Nyquist frequency in our case and will not bias the results significantly.We also try the measurement on a 1024 3 grid, and little difference can be observed within the k range we use.More details about the FFT related procedure can be found in Appendix B2 of Sugiyama et al. (2018a).

Covariance matrix estimation
We perform the delete-one jackknife (JK) method (Sugiyama et al. 2018a) to estimate the data covariance through resampling of the date itself.We employ the kmeans algorithm9 on the random catalog to divide the group samples into N JK spatial subregions on the sky (Kwan et al. 2017).We employ N JK = 200, in which N = 95 so that their ratio is close to corresponding sky area ratio.The covariance matrix in the JK scheme is estimated by . Left: Multipoles of the survey window function for the baseline sample, which are used to compute the masked theoretical pairwise kSZ power dipole.Solid and dashed lines indicate the monopole and qradrupole separately.The red and blue lines represent the window functions of NGC and SGC respectively.Right: Masked (solid lines) and unmasked (dashed lines) theoretical predictions of the pairwise kSZ dipole for the baseline sample.The red and blue lines are for NGC and SGC respectively.On the lower panel, the relative difference Diff = P mask /P nomask − 1 induced by the window function effect is presented.The inflation of "Diff" at k ∼ 0.15h/Mpc is due to the zero-crossing of P ℓ=1 kSZ at that scale.
where y a i = P ℓ=1,a kSZ (k i ) is the power spectrum dipole at k i of the ath resampling data configuration, and its expectation value is ȳi = 1 N JK NJK a=1 y a i .
(38) Hartlap et al. (2007) rescales the inverse covariance matrix C −1 JK as to get an unbiased precision matrix Ĉ−1 .Although the original derivation of the Hartlap factor is not based on the JK technique, we apply this factor here to calculate a conservative precision matrix for our parameter inference procedure.

Survey window functions
In this paper, we follow the treatment suggested by Wilson et al. (2017); Beutler et al. (2017) to include the window function effect caused by the masked survey sky in the kSZ power spectrum model.In this methodology, the window function is first calculated in configuration space and multiplies the pairwise kSZ correlation function to give the masked correlation function.This masked correlation function is then Hankel transformed to produce the masked pairwise kSZ power spectrum.
The window function multipoles characterizing the distortion of the survey geometry are calculated by the random catalog and are given by where RR is the number of random pairs in a logarithmic s bin 10 .Here, Q ℓ is normalized to enforce Q 0 (0) = 1.In order to speed up the calculation, for samples with M th ≥ 10 13 M ⊙ /h, we randomly dilute the random catalog to the size of the M th = 10 13 M ⊙ /h group sample, and for samples with M th < 10 13 M ⊙ /h, we dilute the random catalog to 10 Corrfunc.mocks.DDsmu mocks the size of the corresponding group sample size.A convergence test in Appendix B shows that this dilution induces negligible systematic effects on the calculation of P ℓ=1 kSZ in Equation ( 43).The left panel of figure 4 shows the multipole components of the survey window function applied in our analysis.We expect that monopole and quadrupole respectively tend to unit and zero when s → 0, and both of them are equal to zero in dimensions beyond the survey volume.The positive quadrupole values reflect that, for a fixed s bin, there are more random pairs found along the LOS than perpendicular to the LOS, and vice versa.
In turn, the masked kSZ power spectrum including the window function can be theoretically calculated by (Sugiyama et al. 2018a) where ξ ℓ kSZ (s) is given by Here, j ℓ is the spherical Bessel function of order ℓ.For the dipole, the expansion of Equation ( 41) can be linearly truncated as in which the contributions of higher multipoles of the kSZ correlation function (e.g., the octopole ξ ℓ=3 kSZ ) are small and can be ignored.The robustness of this linear truncation is verified in Appendix C. Equation ( 43) is the final expression used to fit the measurement of the pairwise kSZ power dipole in this work.Its detailed derivation is presented in Appendix C of Sugiyama et al. (2018a).
The theoretical predictions of the pairwise kSZ dipoles with and without the survey window functions are shown on the right panel of figure 4. The amplitude of the masked P ℓ=1 kSZ becomes smaller than the true one at large-scales, due to the lack of group pairs beyond the survey volume.The window function effect becomes significant at k < 0.05h/Mpc and is of the order of 10% − 20% around k ∼ 0.01h/Mpc.Since our measurement S/N is dominated by linear scales, this effect is a key systematic error that ought to be corrected in our analysis.

Imaging systematics
It has been well demonstrated that a variety of observation conditions, such as stellar contamination, Galactic extinction, sky brightness, seeing, and airmass, can introduce spurious fluctuations in the measured galaxy density field (e.g., Scranton et al. 2002;Myers et al. 2006;Morrison & Hildebrandt 2015;Kitanidis et al. 2020).These fluctuations, known as imaging systematics, must be addressed in order to obtain reliable results for large-scale clustering studies.The group catalogs used in this work are derived from photometric galaxy catalogs, which may be subject to such biases.
A widely used approach to reduce the imaging systematics is to assign weights to the targets so that the weighted sample does not show any correlation with the various observation conditions (e.g.Ross et al. 2017;Bautista et al. 2018).This can be achieved by minimizing the difference between the predicted and observed target surface density.The predicted target density is simply a function of various imaging maps, which are designed to not contain any cosmological signal.Previous studies assumed a linear functional form (e.g., Ross et al. 2011), which may not be able to capture the nonlinear dependence on imaging systematics in areas with strong contamination, such as close to the Galactic plane (e.g., Ho et al. 2012).Recently, a machine-learning-based approach has been used to capture these complex relationships (e.g., Rezaie et al. 2020;Chaussidon et al. 2022).
We reduce the imaging systematics in our targets using the random forest algorithm, which is available in the GitHub code regressis11 developed by Chaussidon et al. (2022).Xu et al. (2023) adapted the code to remove imaging systematics in autocorrelations and cross-correlations between multiple tomography bins in the galaxy samples from DR9 of the DESI Legacy Imaging Survey.Here, we follow the steps outlined in Xu et al. (2023) to mitigate imaging systematics and refer interested readers to their paper for details.Left: the pairwise kSZ dipole measurement of the baseline sample without (data points with error bars) and with (dashed line) the imaging-systematics correction.Right: The multipoles of redshift space density power spectrum without (data points with error bars) and with (dashed line) the imaging-systematics correction.The lower panel shows the fractional difference between those of "corrected" and "raw" ones.Measurements are made from the baseline sample.
The comparison between the measured P ℓ kSZ with and without correcting imaging systematics is shown on the left panel of figure 5. We see little systematic offsets when we ignore the imaging systematics.On the right panel of figure 5, we also find little differences between redshift space density power spectrum multipoles Pℓ with and without the imaging-systematics correction, in which Pℓ is analyzed in Appendix A to constrain free parameters and determine the theoretical P pv template in the τ fitting.As a result, we consider the imaging systematics as a negligible systematic error source in our analysis, and we will not correct them in this work.

RESULTS
In this section, we presents main results of this paper, including the pairwise kSZ dipoles, the fitted τ values, and the τ − M relationship.

Fitting procedure
We estimate the average optical depth τ of each group sample by fitting the measured kSZ dipole with the template given by Equation ( 43).τ is the only free parameter to be fitted.The χ 2 of the fitting is where Ĉ−1 is the precision matrix, P is the measured kSZ dipole with the vector size N bin = 15, and P is the theoretical with τ as the single variable.
The best-fitted value and the associated statistical error of τ in this linear fitting can be analytically given by (Chen et al. 2022) τbestfit = P T Ĉ−1 P (τ = 1) Finally, we calculate the S/N of the measurement by

Pairwise kSZ dipole and τ
In Appendix D, we vary the lower mass threshold to select the heaviest 1.6 × 10 7 to 1.8 × 10 5 groups.For each sample, we change θ AP and see how the measurement S/N evolves with it (figures 14).We find that the baseline sample, namely the heaviest ∼ 5.65 × 10 6 groups, has the highest S/N when θ AP = 4.2 ′ .We present its dipoles and associated results here.
in the literature, and is the highest S/N kSZ detection in Fourier space currently.As expected, the useful kSZ signal dominates only at linear scales since it is diluted by the photo-z damping and noises (shot noise, CMB detector noise, and primary CMB contamination) at nonlinear scales.The nonlinear kSZ model (solid lines) gives a good fit to the data, and the fitted τ = (3.96± 0.76) × 10 −5 for M = 1.75 × 10 13 M ⊙ /h.On the right panel of figure 6 we present the correlation coefficients of the dipole measurements from baseline sample.The off-diagonal elements are subdominant in the covariance matrix, and the kSZ field on each k scale evolves independently with each other.
By varying θ AP , we display in figure 7 the measured τ profile of the baseline sample.This profile encodes valuable information about the gas distribution within halos, and it can be fitted by some gas density profile model to constrain the gas fraction f gas within halos (Amodeo et al. 2021;Sugiyama et al. 2018a).In a companion paper, we will study this τ profile with adequate profile models, e.g., the baryonification models (Schneider & Teyssier 2015;Giri & Schneider 2021;Aricò et al. 2021), and constrain the halo gas content in our samples.
In addition, we make a null test by targeting N random points within the DESI sky region, where N is the group number of the base line sample.Going through the same pipeline, the measured τ profile of this random sample is consistent with zero, as shown in figure 7.This null test further validates the robustness of our kSZ detections.

τ − M relationship
Next, we explore the relationship between the average optical depth and the median halo mass of a sample.We choose the measured optical depths of different samples in appendix D measured with θ AP = 4.2 ′ and fit the τ − M relationship with a simple power law, log τ = γ(log where β and γ are two free parameters. The results are shown in figure 8.The red data points with error bars are measured τ s associated with the median group sample masses, and the red solid line is the best-fitted Equation (48).The fitted γ = 0.55 +0.1 −0.1 is roughly consistent with the τ ∝ M relation, as expected for massive groups, and this result confirms our motivation of sample selection with the group mass.In Chen et al. (2022), a higher value of γ = 1.25 ± 0.14 is found.Although being deviated from each other, the discrepancy between two results is not as large as it looks here.Since group samples in both works are not independent of each other by definition, the uncertainty of γ is moderately underestimated in both papers.Albeit, with this, the disagreement can come from several factors.(1) A minor one is that two works use different AP filter radii, with θ AP = 3.4 ′ in Chen et al. (2022), and θ AP = 4.2 ′ here.(2) In Chen et al. (2022), most group samples used to fit Equation ( 48) are heavy clusters with masses larger than 10 14 M ⊙ /h, while in our work the mass range of group samples is ranging from 10 13 M ⊙ /h to 10 14 M ⊙ /h.(3) The redshift distribution of the group sample is affected by the mass threshold, so the sample redshift distributions are not identical to each other, and can introduce potential systematic biases in the measured τ − M relationship.(4) Two works treat the density biases (b 1 , b 2 ) and the photo-z uncertainty σ pho−z in different manners.In this work, we apply a joint analysis of kSZ and RSD effects on the same data set.We first constrain b 1 , b 2 and σ pho−z from the measured Pℓ , these determined free parameters are then adopted to determine P pv and fit τ .In Chen et al. (2022), they evaluated b 1 from the Sheth01 b 1 − M model (Sheth et al. 2001).Then, ξ pv was directly measured from the simulation halo catalog constructed to match the calculated b 1 , on which the photo-z uncertainties are directly applied on each halo.
In Chen et al. (2022), a fiducial sample of the richest 1.2 × 10 5 groups with log M = 14.19 was constructed from a previous version of our adopted group catalog, which is based on the Data Release 8 of the the DESI Legacy Imaging Survey.The pairwise kSZ correlation function was measured from the fiducial sample, and the optical depth is fitted to be τ = 1.66 ± 0.35 × 10 −4 , namely a ∼ 5σ detection as well.We confirm that this measurement is consistent with our fitting result regarding Equation ( 48) within errors, and this is a good cross-check of the kSZ data analysis in both configuration and Fourier spaces12 .

CONCLUSION AND DISCUSSION
In this work, we report a ∼ 5.2σ detection of the kSZ effect in Fourier space, by combining the DESI galaxy groups and the Planck data.By implying the pairwise estimator, we measure the density-weighted pairwise kSZ power spectrum of the galaxy group sample, together with the Planck HFI 217 GHz data.The galaxy group catalog is constructed from the DR9 data of the DESI Legacy Imaging Surveys, and is selected by varying the lower limit of the group mass, which is estimated from the group luminosity (Yang et al. 2021).The baseline sample we construct has a lower mass limit of 10 13 M ⊙ /h and a median mass M = 1.75 × 10 13 M ⊙ /h.It gives the highest S/N= 5.2 when the adopted AP filter radius θ AP = 4.2 ′ .By fitting τ s from samples with different median M , we obtain a linear log τ = γ log( M − 14) + log β, in which γ = 0.55 +0.1 −0.1 .The optical depth τ is constrained by implementing a joint analysis of kSZ and RSD effects on the sample group catalog.The RSD analysis determines the group bias parameters b 1 , b 2 and the parameter σ z quantifying the photo-z damping effect.Then, the best-fitted three parameters are substituted to the kSZ analysis and constrain the optial depth τ .The stability of the τ measurement against variations on the photo-z damping functional form and the prior knowledge of b 1 from the b 1 − M model is verified in figure 11.Although the degeneracy exists between τ and b 1 , and between τ and the photo-z damping effect, it does not induce significant bias of the measured τ in this work.Alternatively speaking, in a joint kSZ and RSD analysis, the measured Pℓ determines the shape and amplitude of P ℓ pv and in turn the corresponding τ .The details about the b 1 , b 2 and σ z fitting results, including how they depend on the adopted photo-z damping functional form, are not so important for the purpose of a τ measurement.Indeed, the joint kSZ and RSD analysis in this work is not rigorous, in that we do not apply a global fitting of measured Pℓ and P ℓ kSZ simultaneously.In the future, this joint fitting of kSZ and RSD effects can serve as a standard practice of the large-scale structure analysis, by which we can constrain the cosmic structure growth and the baryonic distribution in and around halos at the same time (Sugiyama et al. 2017;Okumura & Taruya 2022).
The kSZ effect is a good indicator of the baryonic distribution in the universe.By varying the AP filter radius, we also measure the τ profiles of the group samples, which is tightly coupled with the baryon density profiles within halos.In a companion paper, we will study these τ profiles with more adequate gas density profile models, e.g., the baryonification models (Schneider & Teyssier 2015;Giri & Schneider 2021;Aricò et al. 2021), and constrain the associated baryon density profiles (Schneider et al. 2022).In turn, we can constrain the baryonic damping of the dark matter density power spectrum, which is a key systematic error of the weak-lensing cosmology (Jing et al. 2006;Schneider et al. 2022;Chen et al. 2023a;Zheng & Zhang 2024).
With the next generation of CMB experiments and galaxy redshift surveys, the detection S/N of the kSZ effect can reach O(100) and beyond (Sugiyama et al. 2018a;Smith et al. 2018).At that time, the kSZ effect can be robustly used to help measure the cosmic structure growth history and constrain cosmological models (Sugiyama et al. 2017;Zheng 2020;Okumura & Taruya 2022).In this sense, the Fourier space kSZ analysis will has its unique advantage on its easier redshift space modeling and the covariance matrix evaluation.This work can be regarded a necessary preparation of this exciting era. .
Left: window function multipoles evaluated from the random catalog corresponding to the NGC subbaseline sample .Right: theoretical masked and unmasked power spectrum multipoles and their fractional differences.
corresponds to a σ z ∼ O(3000)km/s.It heavily damps the measured monopole P0 even at linear scales.This damping naturally generated a degeneracy between b 1 and σ z , which in turn affects the τ measurement.We will test the potential bias of the τ measurement due to an improper treatment of the photo-z uncertainty in this appendix.Moreover, for a spec-z sample, P2 is positive, and P4 is nearly zero, while in figure 9 P2 becomes negative, and P4 is away from zero in a significant level.The photo-z uncertainty leaves significant imprints on P2 and P4 , from which we can potentially infer its properties.
On the right panel of figure 9, we show the cross-correlation coefficient matrix between Pℓ measurements.The large photo-z uncertainty induces significant mode couplings between k bins.It results in the large cross-correlations between different k bins of P0 .However, this argument does not apply to P2 and P4 which are dominated by the photo-z uncertainty itself.

A.2. Fitting results
During the P ℓ fitting, we include the window function effect in the RSD model, as described in Wilson et al. (2017).As an example, the window function multipoles of the NGC subbaseline sample are shown on the left panel of figure 10, and differences between the masked and unmasked P ℓ models are plotted on the right.As shown, it is necessary to include the window function effect in the model since there are up to 7% differences between models with and without masks at the fitting scales 0.05h/Mpc ≲ k ≲ 0.15h/Mpc.
We exclude multipole data points at k < 0.05h/Mpc during the fitting, shown in the grey regions of figure 9, as these low k data points of P0 and P4 represent a spurious upturn feature when k → 0. As discussed in Section 4.8, accounting for the imaging systematics only mildly alleviates the irrationality.Other systematics such as the angular and radial integral constraints (de Mattia & Ruhlmann-Kleider 2019), or some unknown sources may be responsible for this irrationality.We leave the investigation of it to the future, and in this work, we fit the measured Pℓ only at 0.05h/Mpc ≲ k ≲ 0.15h/Mpc.
The PDF of the photo-z uncertainty is a complicated, possibly asymmetric function in its nature.Modeling its impact on the redshift space galaxy power spectrum is a difficult topic.In particular, improper modeling of the photoz damping can lead to systematic errors of the bias and σ z parameters.Here, we test this possibility by trying three terms of photo-z damping functions, which sequentially correspond to a Gaussian photo-z error PDF, a Lorentzian photo-z error PDF and a Lorentzian pairwise photo-z error PDF (Davis & Peebles 1983;Ballinger et al. 1996): All PDFs are assumed to be symmetric.Different lines on the left panel of figure 9 are the best-fitted RSD models incorporating the above damping functions.Visually, P pho−z s with the Lorentzian D pho−z best fits the data.Therefore, in this work we will mainly adopt the Lorentzian damping term in modeling P pho−z s and P pho−z kSZ .On the left panel of figure 11, we present the fitted b 1 , b 2 and σ z of the baseline sample Pℓ measurements.Contours with different colors represents different D pho−z s.Obviously, the damping functional form imposes a huge impact on the fitting results, implying the importance of a correct photo-z damping effect model in the photo-z galaxy clustering analysis.In particular, not only do we see the b 1 − σ z degeneracy here, it is also noticed that both b 1 and σ z are degenerated with the photo-z damping function form.
We can compare the fitted b 1 and σ z with some prior knowledge we have.As discussed in Section 2.2, we expect that the photo-z uncertainty σ pho−z of our samples to be 0.005 ∼ 0.025.It corresponds to σ z ≡ cσ pho−z = 0.15 × 10 4 ∼ 0.75 × 10 4 km/s.This expectation is closer to the fitted σ z of the Gaussian photo-z damping function, and lower than those of the Lorentzian2 and Lorentzian functions.On the upper right panel of figure 11, we compare the fitted b 1 s with predictions from several b 1 − M models in the literature.Within three D pho−z models, the b 1 − σ z degeneracy indicates a higher b 1 from a higher fitted σ z .Moreover, the fitted b 1 s of low M samples are comparable to model predictions, while, for high M samples, our fitting results are higher than model predictions.Figure 12.Left: Window function multipoles evaluated from the random catalog corresponding to the NGC subbaseline sample.We use the random whose sizes are equal to 2, 1, and 0.5 times of the group sample size Ngroup.The lower panel shows the fractional difference between those of "2 × Ngroup" and "1 × Ngroup" (blue lines) or "0.5 × Ngroup" (green lines).Right: theoretical masked dipoles and their fractional differences.
These discrepancies, although not necessarily meaning that any of the fitted b 1 and σ z are definitely wrong, probably indicate an inadequate model of the photo-z damping effect.Luckily, the main target of this work is to determine an accurate τ , with the accuracy of the fitted b 1 and σ z being less significant.On the lower right panel of figure 11, we compare fitted τ s of different samples when various photo-z damping functions and b 1 − M models are adopted.It is shown that fitted τ s are quite insensitive to these model uncertainties, so we are confident that the main results of this work are immune from our poor understanding of the photo-z uncertainty.

B. CONVERGENCE TEST OF THE Q ℓ ESTIMATION
In order to speed up the window function calculation, for samples with M th ≥ 10 13 M ⊙ /h, we randomly dilute the random catalog to the size of the M th = 10 13 M ⊙ /h group sample, and for samples with M th < 10 13 M ⊙ /h, we dilute the random catalog to the size of the corresponding group sample size.On the left panel of figure 12, we compare the estimated Q 0 and Q 2 of the NGC subbaseline sample by using 1/10, 1/20, and 1/40 of all random points, corresponding to a random sample size of 2, 1, and 0.5 times the group sample size.We observe good convergence of Q ℓ , in particular that for Q 0 , no systematic offsets are presented.On the right panel of figure 12, we compare the masked theoretical dipole evaluated by Equation (43) using different sizes of random samples, and show that the dilution induces little systematic errors in the P ℓ=1 kSZ calculation.Within the k range of our interest, similar tests on other group samples give at most 1% systematic offset of the P ℓ=1 kSZ calculation for the M th > 10 14 M ⊙ /h group sample, which is well within the statistical error of τ estimation and will not affect our conclusions.

C. MASKED KSZ DIPOLE WITH HIGHER ORDER EXPANSION
The next-to-leading order expansion of the pairwise kSZ power dipole with higher multipoles of the kSZ correlation function is shown as P ℓ=1 kSZ (k) = −i4π ds s 2 j ℓ=1 (ks) ξ ℓ=1 kSZ (s) Q 0 (s) +  and in turn, it makes sense to ignore the higher multipoles of the kSZ correlation function in evaluating the theoretical masked kSZ dipole.

D. RESULTS OF GROUP SAMPLES WITH DIFFERENT MASSES
In this work, we select the group sample by its lower mass threshold M th , based on the fact that τ ∝ M , as validated by figure 8.By varying M th and θ AP , the measurement S/Ns of different samples with various AP filter radii are shown on the right panel of figure 14.We choose the sample with the highest S/N as the baseline sample, which includes 5.65 × 10 6 heaviest groups.On the left panel of figure 14, the τ profiles of different group samples are shown, by which the mass dependence of the gas distribution within and around dark matter halos will be studied in a companion paper.

Figure 1 .
Figure 1.The sky distribution of groups (left) and random points (right) of the baseline sample shown in the Galactic coordinate with N side = 256 in the HEALPix grid frame.The color bar represents the number of groups or random points in each pixel.There are ∼ 5.65 × 10 6 groups in the baseline sample, and the random catalog size is 20 times of the corresponding group sample.
) instead of the linear model, we apply a nonlinear kSZ model developed in Zheng (2020); Xiao & Zheng (2023) to fit the measured kSZ power spectrum, and since the galaxy group data are from the photometric observation, we further include the photo-z effect in the model; (5) we measure and analyze the redshift space density power spectrum from the same group catalog to obtain constraints on the galaxy bias parameters b 1 , b 2 and photo-z error dispersion σ pho−z , by which we break the b 1 − τ degeneracy in a self-consistent way and avoid possible systematic errors coming from manually selected b 1 − M relation and σ pho−z .

Figure 2 .
Figure 2. The mass, redshift, angular radius and richness distribution of the baseline sample groups.The median values of {M, z, θ halo , richness} are {1.75 × 10 13 M⊙/h, 0.64, 1.5 ′ , 3} for both NGC (red solid lines) and SGC (blue solid lines) groups.The means of these properties are presented by the dashed lines.

Figure 3 .
Figure 3. Theoretical predictions of the unmasked pairwise kSZ multipoles of the NGC sub-sample of the baseline catalog.Solid lines represent multipoles of the nonlinear density-weighted pairwise kSZ power spectrum with photo-z errors (Equation (17)) and dashed lines are those without photo-z errors (Equation (14)).The red and blue lines are for dipoles (ℓ = 1) and octopoles (ℓ = 3) respectively.The adopted redshift is the median redshift z = 0.64 of the NGC sub-sample, and τ = 3.96×10 −5 is the best-fitted τ value of the baseline sample.
Figure 5.Left: the pairwise kSZ dipole measurement of the baseline sample without (data points with error bars) and with (dashed line) the imaging-systematics correction.Right: The multipoles of redshift space density power spectrum without (data points with error bars) and with (dashed line) the imaging-systematics correction.The lower panel shows the fractional difference between those of "corrected" and "raw" ones.Measurements are made from the baseline sample.

Figure 7 .
Figure 7. Left: Measured τ profiles of the baseline sample by choosing different AP filter radii (θAP).The blue dots with error bars represent the τ profile of the baseline sample, and the black dots with error bars show the results of the null test by targeting a random sample having the same size of the baseline sample.On top of the figure, we label the comoving radius corresponding to θAP.Right: corresponding detection S/Ns.

Figure 8 .
Figure 8.The measured τ − M relationship with θAP = 4.2 ′ .The points with error bars are measurements from various samples with different M th and the solid lines are the best-fitted Equation (48) with β = (0.98+0.15  −0.15 ) × 10 −4 and γ = 0.55 +0.10 −0.10 .We caution that cross-correlations between different τ s are ignored in this fitting, so the fitted error bars for β and γ are underestimated here.

Figure 9 .
Figure 9. Left: Power spectrum multipoles of the baseline sample compared with the theoretical templates with different photo-z damping terms.We exclude data points in the grey regions during the fitting process.Right: Correlation coefficients of the measured baseline sample multipoles.

Figure 11 .
Figure 11.Left: Markov Chain Monte Carlo (MCMC) fitting results of multipoles of power spectrum of our baseline sample with three free parameters.Different colors indicate results of different photo-z damping functions.Upper Right: the b1 − M relationship, where b1s are obtained by MCMC fittings (data points with error bars) or given by theoretical models (solid lines).Lower Right: Measured τ − M relationships when different D pho−z s and b1 − M models are adopted.The black data points with error bars are results when we fit b1, b2, and σz simultaneously.The colored points are results when we first calculate b1 from the corresponding b1 − M model and then fit b2 and σz from Pℓ measurements.The adopted b1 − M models embedded in the colossus Python toolkit (Diemer 2018) are: jing98 (Jing 1998), sheth01 (Sheth et al. 2001), tinker10 (Tinker et al. 2010) and bhattacharya11 (Bhattacharya et al. 2011).

Figure 13 .
Figure13.Left: higher order multipoles of the window function for the NGC sub-sample of the baseline sample.Right: corresponding theoretical kSZ dipoles.The red dashed line represents the true dipole without convolution with the window function, the blue solid line is the masked dipole calculated by Equation (C4), and the green dotted-dashed line is the one given by Equation (43).The relative difference Diff = Equation (43)/Equation (C4) − 1 is of the order of 0.1% within the k range of our interest.

Figure 14 .
Figure 14.Left: the τ profiles of group samples with different lower mass thresholds.N in the legend represents the size of the sample.For a clear presentation of the error bars, data points are shifted horizontally from their real θAPs.Right: corresponding S/Ns of the measurements.