© The Authors 2023

1. Introduction

The standard cosmological model successfully explains current observations with a spatially flat universe composed of 68% dark energy and 32% matter, the latter dominated by nonbaryonic cold dark matter (80%) and the rest as ordinary baryonic matter (e.g., Planck Collaboration VI 2020). Gravity arranges dark matter into filaments that the baryons, in the form of diffuse gas and galaxies, fill out, the whole making the cosmic web (Bond et al. 1996). Nongravitational mechanisms involving cooling and feedback drive the galaxy evolution and determine the relative distributions of gas, galaxies, and dark matter within the web. Measurements of their joint distribution tell us about these mechanisms. Notably, only about 10–15% of the baryons reside as stars and gas in galaxies today, implying that feedback prevented the majority from cooling into stars and dispersed it as a diffuse gas into galaxy halos (the circumgalactic medium, CGM) and the broadly distributed intergalactic medium (IGM; Fukugita & Peebles 2006; Tumlinson et al. 2017). Understanding these processes is central to galaxy formation theory.

Large astronomical surveys spanning the electromagnetic spectrum now enable a cosmic census of the dark matter, galaxy, and gas distributions. Spectroscopic galaxy surveys provide detailed information on the galaxy distribution as a function of galaxy properties, such as stellar mass. Over large scales, galaxies trace the dark matter field with a linear bias (for a review, see Desjacques et al. 2018). This bias is crucial in galaxy clustering studies since it is degenerate with σ₈, the amplitude of the matter fluctuations. On small scales, the relation between galaxies and the underlying dark matter distribution becomes more complex and nonlinear.

Gravitational lensing directly probes the dark matter distribution. Lensing distorts the images of distant galaxies and permits reconstruction of the foreground mass distribution (e.g., Kaiser & Squires 1993; Bartelmann & Schneider 2001; Chang et al. 2018). Lensing also leaves a characteristic imprint on the anisotropies of the cosmic microwave background (CMB), which presents a back-light probing the matter distribution across the observable Universe (Lewis & Challinor 2006). The effect was first detected through a cross-correlation of Wilkinson Microwave Anisotropy Probe (WMAP) sky maps and tracers of large-scale structure (Smith et al. 2007; Hirata et al. 2008). The most recent work from ); 2008PhRvD..78d3520H. The most recent work from Planck (Planck Collaboration VIII 2020) measured the lensing signal in the CMB at 40σ significance. The South Pole Telescope (SPT) detected the lensing amplitude at 39σ significance using temperature and polarization measurements, and at 10.1σ significance when using only polarization; and the Atacama Cosmology Telescope (ACT) has constructed a lensing map covering 2100 square degrees (Wu et al. 2019; Das et al. 2014; Darwish et al. 2021).

The distribution of diffuse ionized gas can be probed through the thermal and kinetic Sunyaev-Zeldovich (tSZ and kSZ) effects. Energy transferred to CMB photons when scattered by free hot electrons creates the unique spectral signature of the tSZ effect relative to the pure thermal CMB spectrum. The amplitude of the tSZ measures the gas pressure along the line of sight. Any bulk peculiar motion of the gas relative to the comoving frame modifies the CMB brightness via the Doppler effect to generate the kSZ signal, which is proportional to the product of the bulk velocity and the Thomson optical depth along the line of sight (Carlstrom et al. 2002; Sunyaev & Zeldovich 1972).

With galaxy and CMB surveys now covering large overlapping areas of extragalactic sky, these three probes can be combined to study the joint distributions of dark matter, galaxies, and gas. Cross-correlations between CMB lensing and galaxies have been performed in a number of studies (e.g., Kuntz 2015; Bianchini et al. 2015; Pullen et al. 2016; Singh et al. 2017; Doux et al. 2018; Aguilar Faúndez et al. 2019; Hang et al. 2021; Shekhar Saraf et al. 2022; Kitanidis & White 2021; Darwish et al. 2021). Most studies aim either to measure the cross-correlation amplitude (which should be equal to 1 in the standard cosmological model) or to constrain cosmological parameters (primarily σ₈ and Ω_m). Galaxies have also been cross-correlated with the tSZ effect to perform a (sometimes tomographic) measurement of the gas pressure (Vikram et al. 2017; Koukoufilippas et al. 2020), or to measure the mass bias (Makiya et al. 2018). Planck Collaboration XI (2013), and Greco et al. (2015) stacked Planck tSZ measurements of the locally brightest galaxies to measure the relation between the galaxy mass and the thermal energy of ionized gas in and around galaxy halos. Amodeo et al. (2021) probed the gas thermodynamics of CMASS galaxies using cross-correlations with the measurements of tSZ and kSZ effects from the ACT, whereas Vavagiakis et al. (2021) also used these effects to probe the baryon content of SDSS DR15 galaxies. Other studies have cross-correlated the tSZ effect with CMB lensing (Hill & Spergel 2014; Pandey et al. 2022), and Yan et al. (2021) cross-correlated the three probes to constrain the redshift evolution of the galaxy and gas pressure bias factors.

Most cross-correlation studies aim at constraining cosmological parameters or at testing fundamental physics; for instance, cross-correlations break the degeneracy between the galaxy bias and cosmology. As mentioned, some studies binned the galaxies in redshift, for example, to obtain tomographic information on dark matter or gas pressure distributions. A number of studies have examined the clustering of galaxies binned according to their properties (color, stellar mass, and star formation rate; Swanson et al. 2008a; Marulli et al. 2013; Durkalec et al. 2018; Guo et al. 2019; Ishikawa et al. 2020; Simon & Hilbert 2021). By binning on galaxy properties, it is possible not only to study the galaxy distribution, but also to probe how the dark matter and gas (or any other cross-correlated field) distributions depend on galaxy properties.

In this work, we perform a cosmic census of the joint distribution of different matter components with the aim of better understanding their relations. We study the dependence of galaxy, gas, and mass distributions on stellar mass, using a halo model formalism to interpret our measurements of galaxy-galaxy, galaxy-convergence, and galaxy-tSZ cross-spectra. The halo model formalism builds upon Cooray & Sheth (2002), Cacciato et al. (2013), van den Bosch et al. (2013) and on More et al. (2015), who also performed an analysis of the clustering of CMASS galaxies binned by stellar mass to constrain the cosmology.

The paper is organized as follows: Sect. 2 details the data sets we used, and Sect. 3 presents our method for measuring the cross-power spectra. Section 4 describes our modeling approach. We present and discuss our results in Sect. 5 and conclude in Sect. 7. The cosmology used in this paper is fixed and follows Planck Collaboration VI (2020) (TT, TE, EE+lowE+lensing+BAO, where T and E denote the CMB temperature and polarization, and BAO stands for baryon acoustic oscillation).

2. Data

2.1. Galaxy catalog

We used the Sloan Digital Sky Survey (SDSS) Baryon Oscillation Spectroscopic Survey (BOSS) Data Release 12 (DR12) CMASS galaxy sample¹ (Alam et al. 2015). The BOSS survey observed galaxies up to redshift z = 0.7 and was divided into two samples (LOWZ and CMASS). The LOWZ sample contains low-redshift galaxies, and the CMASS sample consists of galaxies in the redshift range 0.43 < z < 0.7.

2.1.1. Subsample definition

We restricted our work to galaxies in the redshift range 0.47 < z < 0.59, following More et al. (2015), in order to work with a sample of nearly constant (as a function of redshift) number density. We also followed More et al. (2015) and used the Portsmouth stellar population synthesis code (Maraston et al. 2013) to retrieve galaxy stellar masses. Maraston et al. (2013) obtained the stellar masses by fitting model spectral energy distributions to u, g, r, i, z magnitudes. They developed two different models, depending on whether a galaxy was assumed to be star forming or passive. We used their results for a Kroupa (Kroupa 2001) initial mass function and chose between the star-forming or passive catalogs according to the simple color criteria they advocate: a star-forming catalog for galaxies with g − i < 2.35, and the passive model otherwise. We then defined four subsamples by stellar mass, M_*, such that log M_*/M_⊙ > 10.8, log M_*/M_⊙ > 11.1, log M_*/M_⊙ > 11.25, and log M_*/M_⊙ > 11.4. The corresponding number of galaxies in each subsample is 473596, 396298, 250964, and 124493, respectively. Figure 1 shows the distribution of CMASS galaxies (in the restricted reshift range 0.47 < z < 0.59) as a function of their color g − i and the stellar mass given by Maraston et al. (2013). Most galaxies are passive, especially at high stellar mass. In the first (second, third, and fourth) stellar mass bin, only 19.51% (13.77%, 11.27%, and 8.90%) of the galaxies are star forming.

Fig. 1. Number of galaxies in a color (g − i) and stellar mass pixel. The horizontal orange lines show the stellar mass thresholds we used for the bins. The vertical red line indicates the color criteria we used to determine whether a galaxy should be considered as passive or star forming.

2.1.2. Galaxy maps and masks

We constructed maps and corresponding masks of the galaxy distribution from the catalogs following an approach similar to that of Reid et al. (2016) and Loureiro et al. (2019). The SDSS collaboration released acceptance and veto masks in MANGLE² format (Swanson et al. 2008b). The acceptance mask represents the completeness of the observations in a given region of the sky; it represents the fraction of galaxies for which a spectrum could be obtained compared to the total number of targets. The veto masks exclude regions of the sky in which no galaxies were observed, for example, because of point spread function (PSF) modeling failures, reduction pipeline timeout due to too many blended objects, regions around the center posts of the plates, or where higher priority objects were observed. We refer to Anderson et al. (2014) for details. In total, about 5% of the covered sky is excluded by the veto masks.

We transformed these MANGLE masks into a single high-resolution binary mask in HEALPix³ (Górski et al. 2005; Zonca et al. 2019) format using its Python implementation healpy with N_SIDE = 8192 (corresponding to a pixel size of 0.43 arcmin), such that any pixel from a veto mask or with a completeness (in the acceptance mask) lower than 0.7 was rejected. We then degraded this binary mask to a lower-resolution mask (N_SIDE = 4096, corresponding to a pixel size of 0.86 arcmin) where the value of each pixel is the mean of its subpixels in the high-resolution binary mask. This yielded our final completeness mask, whose values are noted C_pix following Loureiro et al. (2019).

As in their study (even though the resolution of the mask we used is different), when we computed the galaxy maps, we only accepted pixels with completeness C_pix > 0.8. This condition is different from the previous cut that was performed on completenesses lower than 0.7 because that cut was made before to degrade the mask resolution to N_SIDE = 4096. As a result, some regions that were not excluded from the high-resolution mask can eventually be rejected in the lower-resolution mask if the neighboring regions are not complete enough. Figure 2 shows the mask.

Fig. 2. CMASS, CMB lensing, and Compton-y masks in Galactic coordinates using the Mollweide projection.

Each galaxy was weighted to account for a number of observational effects, such as fibre collisions or redshift failures. Following Anderson et al. (2014), we assigned a weight to each galaxy,

$$\begin{array}{c}\hfill w_{\mathrm{tot}}=(w_{\mathrm{cp}}+w_{\mathrm{zf}}-1)w_{\mathrm{star}}{w}_{\mathrm{see},}\end{array}$$(1)

where w_cp and w_zf account for missing redshifts due to close pairs and redshift failures by upweighting the nearest neighbor galaxy, respectively. The w_star and w_see weights correct for the fact that fewer galaxies are observed in regions with high stellar densities and due to seeing conditions, respectively (for more details, see Ross et al. 2012). We then constructed a galaxy count map as

$$\begin{array}{c}\hfill n_{\mathrm{p}}=\{\begin{array}{cc}\frac{1}{C_{\mathrm{pix}}^p}{\sum }_{i\in p}{w}_{\mathrm{tot}}^i\hfill & \text{if}{C}_{\mathrm{pix}}^p>0.8\hfill \\ 0\hfill & \text{otherwise},\hfill \end{array}\end{array}$$(2)

where n_p is the weighted number of galaxies in pixel p. As mentioned previously, we dropped pixels with a completeness lower than 0.8. Finally, we built overdensity maps defined as

$$\begin{array}{c}\hfill {\delta }_{\mathrm{p}}=(\frac{n_{\mathrm{p}}}{\overline{n}}-1),\end{array}$$(3)

where $\overline{n}$ is the mean of the weighted galaxy count map.

2.2. (CMB lensing) convergence map

We used the CMB lensing map provided by the Planck⁴ Collaboration in their 2018 release and described in Planck Collaboration VIII (2020). The map was computed using a minimum-variance estimate from temperature and polarization data and covers 67% of the sky. It provides spherical harmonic coefficients for the lensing convergence, κ, with N_SIDE = 4096, together with the corresponding mask given at N_SIDE = 2048. The convergence is related to the lensing potential, ϕ, through

$$\begin{array}{c}\hfill {\kappa }_{\ell m}=\frac{\ell (\ell +1)}{2}{\varphi }_{\ell m}.\end{array}$$(4)

We used healpy to transform these spherical harmonics into a real-space sky map.

2.3. tSZ map

For the tSZ effect, we used the all-sky Compton-y parameter map given in Planck Collaboration XXII (2016) and released as part of the 2015 Planck data release. The map and associated mask are provided in HEALPix format with N_SIDE = 2048. We applied the Planck 60% Galactic mask combined with the point-source mask to remove Galactic emission and radio and infrared sources contaminating the y map. We also corrected for cosmic infrared background (CIB) contamination as described in Sect. 3.4.

3. Measurement method

3.1. Estimating the power spectra

We decomposed each map, a (galaxy density map, convergence map, and Compton-y map), a scalar function of the sky position $\widehat{\theta }$, into spherical harmonics as

$$\begin{array}{c}\hfill a(\widehat{\theta })={\displaystyle \sum _{\ell m}{a}_{\ell m}{Y}_{\ell m}(\widehat{\theta }).}\end{array}$$(5)

The angular power (cross)-spectrum between two maps (a and a′), $C_{\mathcal{l}}^{aa\prime }$, is defined by

$$\begin{array}{c}\hfill ⟨a_{\ell m}{a}_{\ell \prime m\prime }^{\prime \ast }⟩=C_{\ell }^{aa\prime }{\delta }_{\ell \ell \prime }{\delta }_{mm\prime },\end{array}$$(6)

where angle brackets indicate an average over the theoretical ensemble of possible skies, and δ represents the Kronecker delta function. An unbiased estimator of the power spectrum, given full-sky coverage, is

$$\begin{array}{c}\hfill C_{\ell }^{aa\prime }=\frac{1}{2\ell +1}{\displaystyle \sum _{m=-\ell }^{\ell }}{a}_{\ell m}{a}_{\ell m}^{\prime \ast }.\end{array}$$(7)

Incomplete sky coverage and masking, however, limit our maps to less than full-sky coverage, as shown in Fig. 2, coupling the different multipoles and reducing the resolution in harmonic space. We therefore binned the multipoles and calculated the bin-bin coupling matrix using the NaMaster⁵ code from Alonso et al. (2019).

The code also corrects for beam smoothing of the Compton-y map. Following Planck Collaboration XXII (2016), we adopted a circular Gaussian beam of 10 arcmin full width at half maximum (FWHM), resulting in beam harmonic coefficients

$$\begin{array}{c}\hfill b_l=exp-\ell (\ell +1){\sigma }^2/2,\end{array}$$(8)

with $\sigma =\mathrm{FWHM}/\sqrt{8ln2}$. Additionally, we used the NaMaster apodization function with a smoothing scale of 10 arcmin to apodize the tSZ and lensing convergence masks.

The multipole binning followed a linear scheme from ℓ = 30 to 190 and logarithmic from 190 to 3800. We used the same scheme for every power spectrum, except for the cross-spectrum between galaxies and Compton-y parameter, for which we only considered multipoles below ℓ = 1250 because, even though we corrected for beam smoothing, the Gaussian approximation might bias the high multipoles since the beam smoothing mainly affects the small scales (i.e., high multipoles). Moreover, there may remain some contamination from point sources and CIB (see Sect. 3.4).

3.2. Computing the noise bias

The only noise spectrum we computed is the galaxy shot noise because we only considered auto-power spectra for galaxies, where the shot noise contributes a noise bias term. The noise being uncorrelated between the different maps, however, means that there is no noise bias term in the cross-power spectra between the galaxy map and the maps of CMB lensing and the tSZ effect. The galaxy shot noise is given by

$$\begin{array}{c}\hfill N_l^{\mathit{gg}}=\frac{4\pi f_{\mathrm{sky}}}{N},\end{array}$$(9)

where N is the weighted number of galaxies in each stellar mass bin, and f_sky is the observed sky fraction.

3.3. Covariance estimate

In order to estimate the covariance of our measurements, we used the formula derived by Shekhar Saraf et al. (2022), which is a generalization of the usual Gaussian covariance matrix that allows us to use the full sky area observed by each probe instead of restricting the three probes to their common overlapping sky region. The expression for the covariance matrix is given by

$$\begin{array}{cc}\hfill {\mathrm{Cov}}_{LL\prime }^{AB,CD}& =\frac{{\delta }_{LL\prime }}{(2{\ell }_L+1)\mathrm{\Delta }\ell f_{\mathrm{sky}}^{\mathit{AB}}{f}_{\mathrm{sky}}^{\mathit{CD}}}[f_{\mathrm{sky}}^{AC,BD}(C_L^{\mathit{AC}}+N_L^{\mathit{AC}})\hfill \\ \hfill & \times (C_L^{\mathit{BD}}+N_L^{\mathit{BD}})+f_{\mathrm{sky}}^{AD,BC}(C_L^{\mathit{AD}}+N_L^{\mathit{AD}})(C_L^{\mathit{BC}}+N_L^{\mathit{BC}})],\hfill \end{array}$$(10)

where A, B, C, and D are indices for galaxies, g, CMB convergence, κ, or Compton-y parameter, y; Δℓ is the multipole bin-width; and $f_{\mathrm{sky}}^{\mathit{AB}}$ denotes the common sky fraction covered by fields A and B.

3.4. CIB decontamination

Redshifted thermal dust emission from star-forming galaxies generates the cosmic infrared background (CIB) that contaminates measurements of the tSZ signal. Since the CIB comes from galaxies, it traces large-scale structures and likely correlates with our other observations. We used a method proposed in a number of studies (e.g., Hill & Spergel 2014; Vikram et al. 2017; Alonso et al. 2018; Koukoufilippas et al. 2020; Yan et al. 2021) to mitigate this CIB contamination.

The observed y map was modeled as

$$\begin{array}{c}\hfill y_{\mathrm{obs}}(\widehat{\theta })=y_{\mathrm{true}}(\widehat{\theta })+{ϵ}_{\mathrm{CIB}}c(\widehat{\theta }),\end{array}$$(11)

where ${\mathit{y}}_{\mathrm{obs}}(\widehat{\theta })$ and ${\mathit{y}}_{\mathrm{true}}(\widehat{\theta })$ are the observed and decontaminated y maps, respectively, and $c(\widehat{\theta })$ is the CIB contamination weighted by some coefficient ϵ_CIB. The cross-power spectrum then becomes

$$\begin{array}{c}\hfill C_{\ell }^{gy,\mathrm{obs}}=C_{\ell }^{gy,\mathrm{true}}+{ϵ}_{\mathrm{CIB}}{C}_{\ell }^{\mathit{gc}}.\end{array}$$(12)

To estimate $C_{\mathcal{l}}^{gc}$, we cross-correlated our galaxy maps with the Planck 545 GHz map that was used as a template for the CIB contamination. Alonso et al. (2018) then estimated the cross-correlation between this map and the Compton-y map such that

$$\begin{array}{c}\hfill C_{\ell }^{cy,\mathrm{obs}}=C_{\ell }^{cy,\mathrm{true}}+{ϵ}_{\mathrm{CIB}}{C}_{\ell }^{\mathit{cc}}.\end{array}$$(13)

Using models (Planck Collaboration XXX 2014; Planck Collaboration XXIII 2016) for $C_{\mathcal{l}}^{cy,\text{true}}$ and $C_{\mathcal{l}}^{cc}$, Alonso et al. (2018) determined ϵ_CIB = (2.3 ± 6.6)×10⁻⁷ (MJy/sr)⁻¹.

We used this value to remove the CIB bias term. For our galaxy sample, this resulted in only a small correction (2% at most) at high multipoles. Given the high uncertainty in ϵ_CIB, we also subtracted the bias term using the 2σ upper value of ϵ_CIB, and this again remained a small correction (up to 13.2%) compared to the uncertainties on our measurements (see Fig. 3). At higher multipoles or with the increasing accuracy of future CMB experiments, the CIB contamination could become an issue for this type of study.

Fig. 3. Impact of CIB decontamination on yg cross-power spectrum. Top: Measured yg cross-power spectrum (in blue) along with the associated uncertainties. The orange (green) curve is the CIB-decontaminated curve using the best-fit (the 2σ upper) value for ϵCIB. Bottom: Relative difference between the CIB-decontaminated curves and our raw measured yg cross-power spectrum. On the scales we considered, the effect of the CIB decontamination is 2.0% at most with the best-fit value of ϵCIB, and 13.2% at most with the 2σ upper value.

4. Modeling approach

4.1. Power spectrum line-of-sight projection

To model our angular power spectra, we computed the three-dimensional theoretical power spectra (see Sect. 4.2) for each probe and projected them along the line of sight using the Limber approximation (Limber 1953), which is applicable except for the lowest multipoles, with kernels W^A(z) adapted to each probe,

$$\begin{array}{c}\hfill C_{\ell }^{\mathit{AB}}={\displaystyle \int }\frac{\mathrm{d}z}{c}\frac{H(z)}{\chi {(z)}^2}{W}_A(z)W_B(z)P^{\mathit{AB}}(k=\frac{\ell }{\chi (z)},z),\end{array}$$(14)

where P^AB is the three-dimensional cross-power spectrum between tracers A and B. In this expression, χ(z) is the comoving distance to redshift z, H(z) is the Hubble parameter at z, c is the speed of light, and k is the comoving wavenumber.

The kernels for the different probes (galaxies, CMB lensing convergence, and Compton-y parameter) are given by

$$\begin{array}{cc}\hfill W_g(z)& =\frac{1}{n_{\mathrm{tot}}}\frac{\mathrm{d}n}{\mathrm{d}z}\hfill \end{array}$$(15)

$$\begin{array}{cc}\hfill W_{\kappa }(z)& =\frac{3}{2}{\mathrm{\Omega }}_m{H}_0^2\frac{(1+z)}{H(z)}\frac{\chi (z)}{c}\left(\frac{\chi (z_{\ast })-\chi (z)}{\chi (z_{\ast })}\right)\hfill \end{array}$$(16)

$$\begin{array}{cc}\hfill W_y(z)& =\frac{1}{1+z},\hfill \end{array}$$(17)

where $\frac{\mathrm{d}n}{\mathrm{d}z}$ is the galaxy redshift distribution, n_tot is the total number of galaxies in the sample, and z_* is the redshift of the last scattering. We neglected the magnification bias contribution to the galaxy kernel, which is more important for high-redshift galaxy samples. We show our derivation for the Compton-y kernel in Appendix C, and more details on the CMB gravitational lensing kernel can be found in Bartelmann & Schneider (2001) and Lewis & Challinor (2006).

4.2. Halo model

In the halo model formalism (for a review, see Cooray & Sheth 2002), all dark matter is assumed to reside in halos. We defined a halo as the spherical region in which the average density is Δ = 200 times the cosmic mean density. This region has a radius r_200m, and the associated mass is M_200m. If not specified otherwise, we simply note M_200m = M. In the literature, halos are also defined as overdensities relative to the critical mass density, and several overdensity values are considered, mainly Δ = 200 or Δ = 500. Conversions between the different definitions require the adoption of a specific halo mass profile, as well as a concentration value.

We worked with M_500c for the gas profile and converted between M_200m and M_500c. In order to make this conversion, we followed Hu & Kravtsov (2003) and assumed a Navarro-Frenk-White (NFW) profile (Navarro et al. 1997) and a concentration following Dolag et al. (2004), to remain consistent with the halo mass function derived by Tinker et al. (2008), who made use of this concentration-mass relation. The conversion was performed using

$$\begin{array}{c}\hfill M_{500\mathrm{c}}=M_{200\mathrm{m}}\frac{{\mathrm{\Delta }}_{500\mathrm{c}}(z)}{{\mathrm{\Delta }}_{200\mathrm{m}}}{c}_{200\mathrm{m}}^{-3}x{\left(\frac{{\mathrm{\Delta }}_{500\mathrm{c}}(z)}{{\mathrm{\Delta }}_{200\mathrm{m}}}f(c_{200\mathrm{m}}^{-1})\right)}^{-3}.\end{array}$$(18)

The quantities in this equation are defined in Appendix B. We probed the NFW profile, and therefore finding a non-NFW profile would require use of the observed profile when performing the mass conversion.

The matter distribution is then defined by the large-scale distribution of halos and the small-scale distribution inside individual halos. For the galaxy distribution, the approach is similar: we specified the spatial distribution of galaxies within halos and how galaxies populate these halos. The latter is given by the halo occupation distribution (HOD).

4.2.1. HOD implementation

Galaxies were split into centrals and satellites. Central galaxies reside at the center of their host halos, and satellites orbit around the centrals within the halo. We adopted an HOD parameterization following Zheng et al. (2005). The mean numbers of central and satellite galaxies in a halo of mass M are

$$\begin{array}{cc}\hfill ⟨N_{\text{c}}|M⟩& =\frac{1}{2}[1+\mathrm{erf}\left(\frac{logM-log{M}_{\min }}{{\sigma }_{logM}}\right)]f_{\mathrm{inc}}(M)\hfill \end{array}$$(19)

$$\begin{array}{cc}\hfill ⟨N_{\text{s}}|M⟩& =⟨N_c|M⟩{\left(\frac{M-M_0}{M_1}\right)}^{{\alpha }_g}\mathcal{H}(M>M_0),\hfill \end{array}$$(20)

where ℋ is the Heaviside function, and erf is the error function, defined as

$$\begin{array}{c}\hfill \mathrm{erf}(x)=\frac{2}{\sqrt{\pi }}{\displaystyle {\int }_0^x}{e}^{-t^2}\mathrm{d}t.\end{array}$$(21)

The mean is taken over all halos of mass M.

The f_inc function was introduced by More et al. (2015) and enables us to take into account the fact that the CMASS sample is incomplete at low stellar masses. Leauthaud et al. (2016) showed that CMASS is about 80% complete at redshift 0.55 at stellar mass log(M/M_⊙) = 11.4, while the completeness is very low at a low stellar mass. The f_inc function is defined as

$$\begin{array}{c}\hfill f_{\mathrm{inc}}(M)=\max (0,\min (1+{\alpha }_{\mathrm{inc}}(logM-log{M}_{\mathrm{inc}}))),\end{array}$$(22)

so that it is a function of the halo mass M and takes values between 0 and 1. This function is defined such that the completeness is 1 for masses M > M_inc. The completeness therefore increases when M_inc decreases or when α_inc decreases.

However, Leauthaud et al. (2016) noted that modeling the galaxy incompleteness with this function is only approximate because this approach assumes that the missing galaxies are a random selection at fixed halo mass and stellar mass, which is not the case because CMASS galaxies are color selected, and there is a correlation between galaxy colors and stellar masses. We leave evaluation of this effect for future work and refer to Saito et al. (2016) for an analysis focused on dealing with CMASS stellar mass incompleteness.

In order to reduce degeneracies, we restricted the number of free parameters by setting M₀ = M_min. Physically, this meant that we imposed that satellite galaxies started forming in halos with a mass high enough to have 0.5 central galaxies on average. Nevertheless, this is not a strong condition because M₀ is difficult to constrain, but has little impact. We also fixed α_g = 1, which is found by a number of studies (e.g., Zheng et al. 2005; Zehavi et al. 2011; More et al. 2015). Since the HOD model we used is phenomenological, this means that we restricted our analysis to a subclass of HOD parameterization with less freedom. This, however, is a more important assumption that we discuss in Sect. 5.2.1.

Additionally, when running our Monte Carlo Markov chain (MCMC) pipeline (see Sect. 5) we added a Gaussian prior on log M_inc centered around log M_min with a width σ = 0.6, which is about twice the 1σ found by More et al. (2015). Adding this prior explicitely assumes that log M_inc should be on the same order of magnitude as log M_min, which is motivated by the idea that low stellar mass galaxies (which are missing in CMASS) probably reside in low-mass halos whose limit we estimate roughly at M_min (the value at which halos contain 0.5 central galaxies on average). Finally, for clarity, we recall the free parameters in the HOD that were varied, which are M_min, σ_{log M}, M₁, α_inc, and M_inc.

4.2.2. Implementing the halo model

The matter distribution is described as a sum of the one-halo and two-halo terms that describe the distribution of matter inside an individual halo and the distribution of halos in space, respectively. In harmonic space, we write

$$\begin{array}{c}\hfill P(k,z)=P^{1\mathrm{h}}(k,z)+P^{2\mathrm{h}}(k,z),\end{array}$$(23)

where P denotes the power spectrum. We used a theoretical framework closely following the approach described in van den Bosch et al. (2013) and Cacciato et al. (2013) to compute each term.

Specifically, van den Bosch et al. (2013) showed that these terms can be computed as

$$\begin{array}{cc}\hfill P_{\mathit{xy}}^{1\mathrm{h}}(k,z)& ={\displaystyle \int }\mathrm{d}M{H}_x(k|M,z)H_y^{\ast }(k|M,z)n(M,z)\hfill \end{array}$$(24)

$$\begin{array}{cc}\hfill P_{\mathit{xy}}^{2\mathrm{h}}(k,z)& ={\displaystyle \int }\mathrm{d}{M}_1{H}_x(k|M_1,z)n(M_1,z)\hfill \\ \hfill & \times {\displaystyle \int }\mathrm{d}{M}_2{H}_{\mathrm{y}}^{\ast }(k|M_2,z)n(M_2,z)Q(k|M_1,M_2,z),\hfill \end{array}$$(25)

where x and y take values in {m, c, s, y}, standing for matter, central galaxies, satellite galaxies, and Compton-y parameter, respectively. For galaxies, these expressions, with the H_x functions described below, are valid assuming that the occupation numbers of centrals and satellites are independent and that the number of satellites follows a Poisson distribution (for more details, see van den Bosch et al. 2013).

The halo mass function, n, is described in the next section, whereas in the simplest case, the function Q reduces to b(M₁)b(M₂)P^lin(M), where b is the halo bias, and P^lin refers to the linear-theory power spectrum. However, in order to account for the radial dependence of the halo bias and the halo exclusion effect, we employed the more detailed approach described in Appendix A.

van den Bosch et al. (2013) gave the functions H_x for the components

$$\begin{array}{cc}\hfill H_{\mathrm{m}}(k|M,z)& =\frac{M}{{\overline{\rho }}_{\mathrm{m}}}{\stackrel{\sim }{u}}_{\mathrm{h}},\hfill \end{array}$$(26)

$$\begin{array}{cc}\hfill H_{\mathrm{c}}(k|M,z)& =\frac{⟨N_{\mathrm{c}}|M⟩}{{\overline{n}}_{\mathrm{g}}(z)},\hfill \end{array}$$(27)

$$\begin{array}{cc}\hfill H_{\mathrm{s}}(k|M,z)& =\frac{⟨N_{\mathrm{s}}|M⟩}{{\overline{n}}_{\mathrm{g}}(z)}{\stackrel{\sim }{u}}_{\mathrm{s}}(k|M,z),\hfill \end{array}$$(28)

with $\stackrel{\sim }{u}$ the Fourier transform of the normalized (meaning that its integral over real space is equal to 1) radial profile of halo mass (subscript h) and satellite distribution (subscript s) (see Sect. 4.2.4), and ${\overline{n}}_{\mathrm{g}}(z)$ the comoving number density of galaxies at redshift z, computed through

$$\begin{array}{c}\hfill n_{\mathrm{g}}(z)={\displaystyle \int }[⟨N_{\mathrm{c}}|M⟩+⟨N_{\mathrm{s}}|M⟩]n(M,z)\mathrm{d}M.\end{array}$$(29)

For the Compton-y parameter, we added

$$\begin{array}{c}\hfill H_y(k|M,z)=\frac{c}{H}\stackrel{\sim }{y}(k|M,z),\end{array}$$(30)

where $\stackrel{\sim }{y}$ is the Fourier transform of the radial distribution $\mathit{y}(r|M)=\frac{{\sigma }_{\mathrm{T}}}{m_{\mathrm{e}}{c}^2}{P}_{\mathrm{e}}(r|M)$ such that

$$\begin{array}{c}\hfill \stackrel{\sim }{y}(k|M)={\displaystyle {\int }_0^{\infty }}4\pi r^{2}y(r|M)\frac{sin(kr)}{\mathit{kr}}\mathrm{d}r,\end{array}$$(31)

with m_e the electron mass and σ_T the Thomson scattering cross-section. In practice, the integral was performed from r_min = 10⁻⁶r_500c to r_max = 50r_500c. We used the pressure profile, P_e(r|M), from Arnaud et al. (2010) with the parameters taken from Planck Collaboration V (2013),

$$\begin{array}{cc}\hfill P_{\mathrm{e}}(r|M)& =1.65{(h/0.7)}^2\mathrm{eV}{\mathrm{cm}}^{-3}\hfill \\ \hfill & \times E^{8/3}(z){\left(\frac{(1-b_{\mathrm{h}})M_{500c}}{3\times {10}^{14}(0.7/h)M_{\odot }}\right)}^{2/3+{\alpha }_{\mathrm{p}}}p(r|(1-b_{\mathrm{h}})M_{500c}),\hfill \end{array}$$(32)

where p is the generalized Navarro-Frenk-White (GNFW) profile (see Sect. 4.2.4 for more details), α_p = 0.12 and E(z) = H(z)/H₀. This profile was fit on clusters assuming that they are in hydrostatic equilibrium (HSE). Since this may not be the case, we converted the true halo mass into the effective hydrostatic mass with the b_h term, also known as the hydrostatic bias. This profile was calibrated on the M_500c mass definition, so that we converted from the M_200m used mostly in this work into M_500c through Eq. (18).

The cross-correlation power spectra involving galaxies have contributions from central and satellite galaxies (which are not individually observable),

$$\begin{array}{cc}\hfill P_{\mathrm{gg}}^{1\mathrm{h}}(k,z)& =2P_{\mathrm{cs}}^{1\mathrm{h}}(k,z)+P_{\mathrm{ss}}^{1\mathrm{h}}(k,z)\hfill \end{array}$$(33)

$$\begin{array}{cc}\hfill P_{\mathrm{gg}}^{2\mathrm{h}}(k,z)& =P_{\mathrm{cc}}^{2\mathrm{h}}(k,z)+2P_{\mathrm{cs}}^{2\mathrm{h}}(k,z)+P_{\mathrm{ss}}^{2\mathrm{h}}(k,z)\hfill \end{array}$$(34)

$$\begin{array}{cc}\hfill P_{\mathrm{gm}}(k,z)& =P_{\mathrm{cm}}(k,z)+P_{\mathrm{sm}}(k,z)\hfill \end{array}$$(35)

$$\begin{array}{cc}\hfill P_{\mathrm{g}y}(k,z)& =P_{\mathrm{c}y}(k,z)+P_{\mathrm{s}y}(k,z),\hfill \end{array}$$(36)

where the one-halo and two-halo terms in the last two equations are identical, so that we do not specify them.

4.2.3. Halo mass function and halo bias

The halo mass function gives the comoving number of halos of mass M per unit volume and logarithmic mass. In the halo model, we assumed that the halo mass function depends only on halo mass and redshift. We used the halo mass function that Tinker et al. (2008) defined and calibrated on numerical simulations,

$$\begin{array}{c}\hfill n(M,z)=\frac{{\overline{\rho }}_{\mathrm{m}}}{M^2}\nu f(\nu )\frac{\mathrm{d}ln\nu }{\mathrm{d}lnM},\end{array}$$(37)

with

$$\begin{array}{c}\hfill f(\nu )={\alpha }_t[1+{({\beta }_t\nu )}^{-2{\mathrm{\Phi }}_t}]{\nu }^{2{\eta }_t}{e}^{-{\gamma }_t{\nu }^2/2},\end{array}$$(38)

where Tinker et al. (2008) found β_t = β_t, 0(1 + z)^0.20, ϕ_t = ϕ_t, 0(1 + z)^−0.08, η_t = η_t, 0(1 + z)^0.27, and γ_t = γ_t, 0(1 + z)^−0.01 with β_t, 0 = 0.589, γ_t, 0 = 0.864, ϕ_t, 0 = −0.729, η_t, 0 = −0.243. The value of α_t was chosen so that matter is not biased with respect to itself, which is imposed by

$$\begin{array}{c}\hfill {\displaystyle \int }b(\nu )f(\nu )\mathrm{d}\nu =1,\end{array}$$(39)

with b the halo bias function and ν the peak height, defined as

$$\begin{array}{c}\hfill \nu (z)=\frac{{\delta }_{\mathrm{sc}}(z)}{\sigma (M)}.\end{array}$$(40)

As pointed out by van den Bosch et al. (2013), the critical threshold for spherical collapse derived by Navarro et al. (1997) is well approximated in the case of a flat universe by

$$\begin{array}{c}\hfill {\delta }_{\mathrm{sc}}(z)=0.15{(12\pi )}^{2/3}\frac{{\mathrm{\Omega }}_{\mathrm{m}}{(z)}^{0.0055}}{D(z)}.\end{array}$$(41)

The remaining term required to compute the peak height, the linear matter variance on the Lagrangian scale of the halo, is given by

$$\begin{array}{c}\hfill {\sigma }^2(M)=\frac{1}{2{\pi }^2}{\displaystyle \int }{P}_{\mathrm{mm}}^{\mathrm{lin}}(k,0){\stackrel{\sim }{W}}^2(kR)k^2\mathrm{d}k,\end{array}$$(42)

where $P_{\mathrm{mm}}^{\mathrm{lin}}$ is the linear matter power spectrum that we calculated using CLASS⁶ (Blas et al. 2011), and $\stackrel{\sim }{W}$ is the Fourier transform of the real space top-hat filter, equal to

$$\begin{array}{c}\hfill \stackrel{\sim }{W}(kR)=\frac{3(sinkR-kRcoskR)}{{(kR)}^3},\end{array}$$(43)

with $R={(3M/4\pi {\overline{\rho }}_{\mathrm{m}})}^{1/3}$.

The halo distribution depends on halo mass, which is accounted for through the mass-dependent halo bias function, b. Numerical simulations show that massive halos tend to be more clustered than smaller halos, but with the constraint Eq. (39). We made use of the simulation-fitted halo bias provided by Tinker et al. (2010),

$$\begin{array}{c}\hfill b(\nu )=1-\frac{A_t{\nu }^{a_t}}{{\nu }^{a_t}+{\delta }_{\mathrm{sc}}{(z=0)}^{a_t}}+B_t{\nu }^{b_t}+C_t{\nu }^{c_t},\end{array}$$(44)

with

$$\begin{array}{ccc}& A_t\hfill & \hfill =1.0+0.24y_t{e}^{-{(4/y_t)}^4},\\ \hfill & a_t\hfill & \hfill =0.44y_t-0.88,\\ \hfill & B_t\hfill & \hfill =0.183,\\ \hfill & b_t\hfill & \hfill =1.5,\\ \hfill & C_t\hfill & \hfill =0.019+0.107y_t+0.19e^{-{(4/y_t)}^4},\\ \hfill & c_t\hfill & \hfill =2.4,\\ \hfill & y_t\hfill & \hfill =log\mathrm{\Delta }.\end{array}$$

The halo mass function provided by Tinker et al. (2008) was normalized at z = 0 to account for all the mass in the Universe: ∫f(ν)dν = 1. At higher redshifts, Eq. (39) approximately imposes this property. However, with our framework, integrals were performed over mass rather than peak height, ν. In practice, when integrating over mass, we may find $\int Mn(M)\mathrm{d}M\ne {\overline{\rho }}_{\mathrm{m}}$, violating this condition. This comes from the fact that even though σ(M) tends to +∞ as M approaches 0, the divergence is only logarithmic with k so that Eq. (42) varies very slowly with mass at the low-mass end. As a result, we chose to slightly modify the value of f at low ν; more precisely, Tinker et al. (2008) stated that the mass function is calibrated for σ⁻¹ > 0.25, which corresponds approximately to M = 10¹⁰ M_⊙, whereas the behavior of the fitting function at lower masses is arbitrary. We therefore computed the mass integrals between M₈ = 10⁸ M_⊙ and M₁₆ = 10¹⁶ M_⊙ and modified the mass function below M₁₀ = 10¹⁰ M_⊙ so that ${\int }_0^{\nu (M_{10})}f(\nu )\text{d}\nu ={\int }_{\nu (M_8)}^{\nu (M_{10})}f(\nu )\text{d}\nu$. We also adapted the halo bias over this mass range so that it equaled its f-weighted mean. We checked that varying the mass boundaries has no significant impact on the resulting power spectrum.

4.2.4. Generalized Navarro-Frenk-White profile

We use the GNFW profile to model the small-scale distribution of dark matter mass and satellite galaxies. This profile, proposed by Nagai et al. (2007), is a generalization of the NFW profile (Navarro et al. 1997),

$$\begin{array}{c}\hfill \rho (r|M)\propto \frac{1}{{(r/r_{\mathrm{s}})}^{\gamma }{(1+{(r/r_{\mathrm{s}})}^{\alpha })}^{(\beta -\gamma )/\alpha }},\end{array}$$(45)

where r_s is a characteristic radius related to r_200m through the concentration parameter c by r_200m = cr_s, and r_200m is defined by $M_{200\mathrm{m}}=(\frac{4\pi }{3})200{\rho }_{\mathrm{m}}{r}_{200\mathrm{m}}^3$. The GNFW profile reduces to the standard NFW profile when α = γ = 1 and β = 3. In order to compute the different functions $\stackrel{\sim }{u}$, we took the Fourier transform of Eq. (45), normalized to unity as required by Eqs. (26)–(28),

$$\begin{array}{c}\hfill \stackrel{\sim }{u}(k|M)={\displaystyle {\int }_0^{r_{200\mathrm{m}}}}4\pi r^2\frac{sin(kr)}{\mathit{kr}}\frac{\rho (r|M)}{M}\mathrm{d}r.\end{array}$$(46)

For the dark matter and galaxy distributions, we fixed α and γ to unity and varied only β in Eq. (45) to test whether the matter and satellite galaxies follow NFW profiles. This parameter is denoted β_m ( β_s) when it relates to the matter (satellite) profile.

In the same way as for the mass conversion, we used the concentration parameter given by Dolag et al. (2004). Even though our profiles are different from NFW (when β ≠ 3), changing β affects the profile so much that the value of the concentration has very little impact. The concentration we used was then

$$\begin{array}{c}\hfill c(M,z)=\frac{c_0}{1+z}{\left(\frac{M}{{10}^{14}{h}^{-1}{M}_{\odot }}\right)}^{{\alpha }_c},\end{array}$$(47)

where c₀ = 9.59 and α_c = −0.102.

For the gas distribution, Planck Collaboration V (2013) derived a GNFW profile for the gas pressure in Eq. (32),

$$\begin{array}{c}\hfill p(r|M)=\frac{P_0}{{(r/r_{\mathrm{s}})}^{\gamma }{(1+{(r/r_{\mathrm{s}})}^{\alpha })}^{(\beta -\gamma )/\alpha }},\end{array}$$(48)

with α = 1.33, β = 4.13, γ = 0.31, P₀ = 6.41, and using a concentration parameter c = 1.81.

We sought to probe the gas pressure in the vicinity of galaxies to determine whether it depends on stellar mass. We had the choice to either assume a hydrostatic mass bias (in Eq. (32)) and vary the pressure normalization P₀, or to fix P₀ and measure the hydrostatic bias. The two approaches are not equivalent because varying the hydrostatic bias has an effect on the profile shape by changing the value of r_s through r_s(1 − b_h)^1/3. We chose the second option because this also enabled us to search for any stellar mass dependence of the hydrostatic bias (assuming that the pressure profile from Arnaud et al. (2010) is indeed universal). As for the HOD, we recall the parameters introduced in this section that were varied for the profiles: β_m and β_s in the matter and satellite profiles, and b_h, the hydrostatic bias.

4.3. Galaxy-matter cross-correlation amplitude

We added another parameter for the cross-correlation between CMB lensing convergence and galaxy maps, the cross-correlation amplitude A. This modifies Eq. (14) to

$$\begin{array}{c}\hfill C_{\ell }^{\kappa \mathrm{g}}=A{\displaystyle \int }\frac{\mathrm{d}z}{c}\frac{H(z)}{\chi {(z)}^2}{W}_{\kappa }(z)W_{\mathrm{g}}(z)P^{\kappa \mathrm{g}}(k=\frac{\ell }{\chi (z)},z).\end{array}$$(49)

In theory, this parameter should be equal to unity, but many studies (Giannantonio et al. 2016; Singh et al. 2017; Hang et al. 2021; Shekhar Saraf et al. 2022; Kitanidis & White 2021) have found different values. This parameter can be used to check for consistency, because a deviation from unity could indicate a possible failure either in the data processing or in the theoretical modeling.

It might also indicate the presence of stochasticity (e.g., Dekel & Lahav 1999; Giannantonio et al. 2016), since stochasticity would change the relation between galaxy and matter overdensities, which is assumed to be statistically deterministic and parameterized through the galaxy bias. If a stochastic component were to break this deterministic relation, it would change the galaxy auto-power spectrum, but not the matter-galaxy cross correlation. The galaxy bias measured from the galaxy auto-power spectrum would differ from the one measured with the cross correlation with CMB lensing, and as a result, A would capture this deviation.

Additionally, A could also measure deviations from General Relativity, for instance, if the two Newtonian potentials entering the perturbed Friedmann-Lemaître-Robertson-Walker line element were to be different. This is the goal of the μ, γ parameterization (Zhao et al. 2009) or of the E_G statistics (Zhang et al. 2007). We did not introduce this parameter for the cross correlation between the Compton-y and galaxy maps because it would be degenerate with the hydrostatic mass bias, which already probes the amplitude of $C_{\mathcal{l}}^{yg}$.

5. Method and results

5.1. Method

With our theoretical model described in the previous section, we used the MCMC sampler emcee⁷ (Foreman-Mackey et al. 2013) to fit our observations. As a summary, our considered observations are the galaxy auto power spectrum, $C_{\mathcal{l}}^{gg}$ (16 observational points after multipole binning), the galaxy-CMB lensing convergence cross correlation, $C_{\mathcal{l}}^{\kappa g}$ (16 observational points), and the galaxy-Compton-y parameter cross-correlation, $C_{\mathcal{l}}^{yg}$ (12 observational points). We also added the galaxy abundance as an observational constraint, which is approximately constant with redshift because we only selected galaxies with redshift z ∈ [0.47, 0.59] where the observations are the most complete. In order to account for the slight variation in the abundance with redshift, we adopted a 20% uncertainty for the abundance constraint in the likelihood.

We varied the parameters M_min, σ_{log M}, M₁, α_inc, and M_inc in the HOD model, adding a flat prior between 0 and 1 on α_inc and a Gaussian prior with σ = 0.6 on log M_inc (see Sect. 4.2.1 for the discussion of this prior). We also varied β_m and β_s in the matter and satellite profiles. However, as we show in Sect. 5.2.2, β_m cannot be constrained with the current data and was therefore fixed to 3 (the value for which the GNFW is simply an NFW) to avoid adding some small amount of noise in the fitting procedure.

Finally, the hydrostatic bias, b_h, and the CMB lensing convergence-galaxy cross-correlation amplitude, A, were also varied. We recall that we fixed the cosmological parameters, for which we take the values from Planck Collaboration VI (2020) (TT, TE, EE+lowE+lensing+BAO): H₀ = 67.66 km s⁻¹ Mpc⁻¹, Ω_bh² = 0.2242, Ω_ch² = 0.11933, τ = 0.0561, n_s = 0.9665, and σ₈ = 0.8102. Even though some of the relations (mass-concentration relation, halo mass function, and bias, etc.) were derived with a slightly different cosmology, we assumed that they are still valid in the Planck cosmology and that a slight difference would be absorbed in our phenomenological parameterization (in the HOD model, e.g.).

We fit each stellar mass threshold bin separately to reduce computation time, meaning that we did not use cross correlations between the different mass bins. As a result, our observational data vectors consist of 45 observational points, and we have eight free parameters. Assuming a Gaussian likelihood,

$$\begin{array}{c}\hfill ln\mathcal{L}=-\frac{1}{2}[{(X(\theta )-X^{\mathrm{obs}})}^T{C}^{-1}(X(\theta )-X^{\mathrm{obs}})+\frac{{(n_{\mathrm{g}}(\theta )-n_{\mathrm{g}}^{\mathrm{obs}})}^2}{{\sigma }_{n_{\mathrm{g}}}^2}].\end{array}$$(50)

where X^T denotes the transpose of a vector X, θ is the parameter vector, X the concatenation of the three angular (cross-) power spectra, C the covariance matrix, n_g the galaxy abundance, and σ_ng the galaxy abundance uncertainty. This likelihood does not show the Gaussian prior added on M_inc. Using the MCMC sampler emcee, we searched for the maximum of this likelihood.

5.2. Results

Figures 4–6 show the theoretical (cross)-power spectra fitted to the observations for each stellar mass bin. The agreement is good overall, and we succeeded in reproducing the spectral amplitudes and shapes at both small and large scales. The goodness of the fits can be better evaluated with the reduced χ² values given in Table 1. The table also gives the parameter medians together with the 68% credible intervals, which are estimated via the 16th and 84th percentiles. We used the python package GetDist⁸ (Lewis 2019) to show in Fig. 7 the degeneracies between the parameters for the less restrictive stellar mass threshold bin (log(M_*/M_⊙) > 10.8), except for the β_m parameter, which we were not able to constrain (see Sect. 5.2.2).

Fig. 4. Observed and best-fit galaxy angular auto-power spectra, $C_{\mathcal{l}}^{\text{gg}}$, for the four different stellar mass threshold bins. The same colors are used for the observed error bars and the corresponding theoretical best fit in each stellar mass bin.

Fig. 5. Observed and best-fit angular convergence-galaxy cross-power spectra, $C_{\mathcal{l}}^{\kappa \text{g}}$, for the four different stellar mass threshold bins.

Fig. 6. Observed and best-fit angular tSZ-galaxy cross-power spectra, $C_{\mathcal{l}}^{y\text{g}}$, for the four different stellar mass threshold bins.

5.2.1. Galaxy distribution

Table 1 shows that M_min increases with the stellar mass, suggesting that galaxies with more stellar mass prefer to occupy more massive dark matter halos (see Eq. (19)), as could be naively expected. At the same time, σ_{log M} also increases and is highly correlated with M_min (see Fig. 7). This is due to the fact that as minimum halo mass for a central galaxy, M_min, increases, σ_{log M} can also increase to compensate for this, offsetting the impact on the mean number of central galaxies in lower-mass halos. This trend and correlation agree with the results found by More et al. (2015).

Fig. 7. Posterior distributions for the stellar mass threshold bin with log(M*/M⊙) > 10.8. The βm parameter is not shown because it could not be constrained.

The parameters α_inc and M_inc, which were included to account for galaxy incompleteness at the low stellar mass end (see Eqs. (19) and (22)), are not well constrained and remain prior dominated, as shown in Fig. 7. As a result, these parameters mainly act as a source of noise representing our lack of knowledge about the galaxies that are missing in the sample. In particular, we see degeneracies between σ_{log M}, α_inc and M_inc; this comes from the fact that their effect on the power spectra is very similar. On the one hand, M_inc and α_inc model the fact that some galaxies are missing in low-mass halos (meaning that these galaxies exist, but are not part of the sample), whereas on the other hand, σ_{log M} changes the number of galaxies in low-mass halos.

Figure 4 shows that clustering increases more rapidly at high multipoles for galaxies with high stellar mass (e.g., compare the slope of the red curve to the blue curve at high ℓ). Since small-scale (high multipole) clustering is primarily driven through the one-halo term, this implies that massive galaxies tend to be more concentrated inside individual halos. In our model, this leads to an increase in β_s (see Eq. (45)), placing satellite galaxies closer to their central galaxy. We indeed find that β_s increases with stellar mass, and we conclude that higher-mass satellite galaxies have steeper radial profiles. This conclusion might change when α_g (see Eq. (20)) were allowed to vary because this also increases the clustering at small scales by raising the number of satellite galaxies. The two parameters are not completely degenerate, however, because α_g changes the total number of galaxies, in contrast to β_s. Because our HOD model is mainly constrained by the galaxy auto-power spectrum, we chose not to vary α_g to reduce the degeneracies between the parameters.

The parameter M₁, which controls the number of galaxies at high halo mass (see Eq. (20)), is insensitive to stellar mass. This can be understood from the top and middle panels of Fig. 8, which show the number of galaxies inside individual halos as a function of halo mass. More galaxies of low stellar mass lie inside less massive halos (see the previous discussion on M_min), whereas the number of galaxies inside high-mass halos is approximately the same for each stellar mass bin. Similarly, the bottom panel of Fig. 8 shows the proportion of galaxies lying in halos of different masses. There is a higher proportion of low stellar mass galaxies in low-mass halos and of high stellar mass galaxies in the most massive halos. The large uncertainties at low halo mass come from the poor constraints on α_inc and M_inc and show that no galaxy in our sample may be residing in a halo with this mass.

Fig. 8. Galaxy distribution as a function of the host halo mass. Top: Total (central + satellite) number of galaxies contained in halos as a function of halo mass, M, evaluated at the sample median redshift. The four curves correspond to the different stellar mass threshold bins. The shaded bands represent the 1σ uncertainty. Middle: Total number of galaxies per comoving volume and halo mass range (in units of Mpc−3 $M_{\odot }^{-1}$) as a function of halo mass. Bottom: Total number of galaxies per halo mass range normalized by the number of galaxies in each stellar mass bin. This fraction of galaxies lies in a halo of mass M. The quantity V appearing in the label stands for the comoving volume.

After constraining the HOD, we computed the satellite fraction in each subsample as

$$\begin{array}{c}\hfill f_{\mathrm{s}}(z)=\frac{\int \mathrm{d}M⟨N_{\mathrm{s}}|M⟩n(M,z)}{\int \mathrm{d}M[⟨N_{\mathrm{c}}|M⟩+⟨N_{\mathrm{s}}|M⟩]n(M,z)}.\end{array}$$(51)

The evolution of the satellite fraction with galaxy stellar mass is shown in Fig. 9. The satellite fraction decreases with stellar mass, so that the most massive galaxies tend to be centrals. This trend is expected, and is similar to what More et al. (2015) found, even though the values are slightly different. These differences may arise because the observations are not exactly the same: These authors also included galaxy weak lensing. Moreover, the satellite fraction is not an observable, but a derived parameter, so that modeling the differences can lead to different results.

Fig. 9. Derived satellite fraction as a function of galaxy stellar mass. The error bars depict the 1σ credible intervals.

Making use of this HOD, we can also derive the galaxy linear bias describing the statistical relation between the galaxy distribution and the underlying dark matter field. The galaxy bias, b_g, is formally defined as δ_g = b_gδ, where δ_g and δ are the galaxy and matter overdensity fields, respectively. In the power spectrum formalism, this yields $P_{\text{gg}}(k)=b_{\text{g}}{}^2{P}_{\text{mm}}(k)$, so that using Eq. (25) in the large-scale linear limit leads to

$$\begin{array}{c}\hfill b_{\mathrm{g}}={\displaystyle \int }\mathrm{d}M[⟨N_{\mathrm{c}}|M⟩+⟨N_{\mathrm{s}}|M⟩]\frac{n(M,z)}{{\overline{n}}_{\mathrm{g}}(z)}b(M,z).\end{array}$$(52)

Figure 10 shows the evolution of galaxy bias with stellar mass, together with the 1σ interval. The increase in linear galaxy bias with stellar mass is clear, which is obvious from the galaxy power spectra in Fig. 4 (where the amplitude of the power spectra increases with galaxy stellar mass). This follows from our result that more massive galaxies populate higher-mass halos (e.g., Fig. 8), which are intrinsically more biased, that is, the halo bias increases with halo mass (this can be seen in the halo bias formula from Tinker et al. (2010), where the bias increases with ν and therefore with the halo mass). The important implication is that galaxy clustering depends on the properties of the galaxies under consideration; not all galaxies are interchangeable tracers of the dark matter field.

Fig. 10. Derived galaxy bias as a function of galaxy stellar mass. The error bars represent the 1σ credible intervals.

5.2.2. Dark matter distribution

The galaxy-CMB lensing convergence cross-correlation amplitude, A (introduced in Eq. (49)), is not sensitive to stellar mass, and its value is consistent with unity for each stellar mass bin. This supports our theoretical modeling as well as the data processing.

Singh et al. (2017) found A = 0.78 ± 0.13 using the same data (Planck CMB lensing and the CMASS survey, but without division into stellar mass bins) and working in configuration space. The difference between their result and ours may arise from our use of a restricted redshift range (0.47 < z < 0.59), from the use of different scale-binning schemes, and/or the fact that we have included much smaller scales in our study. This might indicate at an effect on large scales (either on the data processing or theoretical side) because Fig. 5 shows that our best-fit model falls in the upper part of the 1σ interval on large scales (small multipoles).

We therefore made another run for the first stellar mass bin, using only multipoles up to ℓ = 492 for the CMB lensing-galaxy cross-correlation. We obtained $A=0.{86}_{-0.11}^{+0.13}$ (compared to $0.{981}_{-0.092}^{+0.091}$ when using all the multipoles available in our analysis), so that we can confirm that A indeed shifts toward lower values when only the largest scales are considered. This value is still consistent with unity. As mentioned previously, the remaining difference between our result and the result from Singh et al. (2017) might be due to the restricted redshift range and to differences in the scale binning schemes. The difference in our results is slight, however, when the uncertainties are taken into account.

The other parameter (apart from HOD parameters) used to fit the galaxy-CMB lensing cross-correlation power spectrum is β_m, which enters into the matter distribution profile parameterization (Eq. (45)). Unfortunately, the uncertainties are too large to quote a constraint on β_m, meaning that we cannot constrain the matter density profile even when A is fixed to 1. However, we examine in Sect. 6 how well β_m may be constrained by the future Simons Observatory CMB experiment.

5.2.3. Gas distribution

The Compton-y cross-correlation probes the gas pressure in the vicinity of CMASS galaxies, and with our approach, this is quantified through the hydrostatic mass bias, b_h (see Eq. (32). As shown in Table 1 or in Fig. 11, this bias increases slightly with stellar mass, indicating that the gas pressure is lower around high stellar mass galaxies. We can also express the trend with the halo-bias weighted thermal pressure, defined as

$$\begin{array}{c}\hfill ⟨b{P}_{\mathrm{e}}⟩(z)={\displaystyle \int \mathrm{d}Mn(M,z)b(M,z){\int }_0^{\infty }\mathrm{d}r4\pi r^2{P}_{\mathrm{e}}(r|M,z),}\end{array}$$(53)

Fig. 11. Hydrostatic mass bias (blue) and halo-bias weighted thermal pressure (orange) as a function of galaxy median stellar mass. The halo-bias weighted thermal pressure is consistent with the results from Pandey et al. (2019) when compared with the similar redshift bin. The two sets of error bars are shifted along the x-axis for visualization.

where we recall that b is the halo bias. This quantity is especially interesting because it is directly measured in the power spectrum on linear scales. This bias is also directly linked to the pressure profile integral, and because the halo bias increases around high stellar mass galaxies, finding a decreasing ⟨bP_e⟩ means that the pressure profile integral must decrease with stellar mass. We note, however, that the data remain consistent with a constant value.

The dependence of ⟨bP_e⟩ on stellar mass is also shown in Fig. 11, and the values found (approximately 0.18 to 0.21 meV cm⁻³) agree quite well with Pandey et al. (2022), who performed the cross correlation between tSZ and galaxy shears. The hydrostatic bias we infer takes values between 0.55 and 0.62, which also agrees very well with Pandey et al. (2022); 1 − b_h = 0.56 ± 0.03) and is slightly lower than Makiya et al. (2018); 1 − b_h = 0.649 ± 0.041) and Ibitoye et al. (2022); 1 − b_h = 0.67 ± 0.03).

As pointed out by Salvati et al. (2018), a mild difference is known on this parameter because tSZ number cluster counts combined with the observation of CMB primary anisotropies lead to low hydrostatic bias values of about 1 − b_h = 0.6, whereas estimates from weak-lensing and X-ray observations favor higher values of about 1 − b_h = 0.8 − 0.85. The hydrostatic mass bias we find is therefore consistent with CMB and tSZ cluster count observations (Planck Collaboration XXIV 2016). However, we emphasize that we kept the cosmology fixed, and allowing the cosmology to vary would broaden the error bars.

Assuming that the profile from Arnaud et al. (2010) is universal and well calibrated, we find a hint that the hydrostatic bias increases (or 1 − b_h decreases) with galaxy stellar mass. This is not statistically significant, but it is a tantalizing glimpse of the potential of this type of study to establish relations between galaxy properties and the gas surrounding them in the circumgalactic and intergalactic media. This may prove a powerful probe of galaxy feedback regulating the baryon cycle, galaxy evolution, and star formation. We explore this potential further in Sect. 6. Interestingly, some recent studies using cluster gas fractions (Eckert et al. 2019; Wicker et al. 2023, and references therein) also find indications that the hydrostatic mass bias increases with cluster mass.

6. Future potential: Forecast

6.1. Forecast specifications

We forecast how well future CMB experiments will be able to constrain the value of β_m, A, and b_h taking the Simons Observatory (SO) specifications as an example (Ade et al. 2019). To do this, we used our theoretical model to create fiducial galaxy auto-power spectrum, $C_{\mathcal{l}}^{\text{gg}}$, cross-correlation galaxy–CMB convergence spectrum, $C_{\mathcal{l}}^{\kappa \text{g}}$, and galaxy–tSZ cross-power spectrum, $C_{\mathcal{l}}^{y\text{g}}$, all assuming an NFW profile for the matter distribution (β_m = 3) and ideal galaxy-matter coefficient (A = 1).

We adopted the first stellar mass bin (of the CMASS sample), meaning that we kept the same $\frac{\mathrm{d}n}{\mathrm{d}z}$ (number of galaxies, dn, per redshift slice, dz) and that we used our rounded best-fit HOD parameters as well as the best-fit hydrostatic bias. In order to also form an idea of potential improvement by using a future galaxy survey experiment, we considered very approximate specifications of the Dark Energy Spectroscopic Instrument (DESI; DESI Collaboration 2016). Considering a sample at similar redshift, the specifications of the luminous red galaxy (LRG) sample can be approximated as a sky coverage of 14 000 square degrees, and a $\frac{\mathrm{d}n}{\mathrm{d}z\mathrm{d}{\mathrm{deg}}^2}\approx 1000$. Using the same redshift range as previously (0.47 < z < 0.59), we obtain approximately 1.7 × 10⁶ galaxies, even though the real LRG survey of DESI will be at slightly higher redshift. This is about four times the number of galaxies we have in our first CMASS stellar mass bin. Because this is only a very rough approximation of the specifications of DESI, we call this case DESI-like, and refer to the cross-correlation between CMASS and SO as the CMASS × SO case, and between the DESI-like survey and SO as the DESI-like × SO case.

We adopted the forecast CMB lensing and tSZ noise curves that are publicly provided by the SO Collaboration (Ade et al. 2019). For the lensing, we considered the noise curve of the minimum variance estimator combining TT, TE, EE, EB, and TB measurements assuming no deprojection and the baseline sensitivity. For the tSZ noise curve, we used the baseline sensitivity curve without deprojection obtained via internal linear combination (ILC) component separation.

These noise curves give the covariance matrix through Eq. (10). We assumed that the common sky fraction is the same as between CMASS and Planck for the CMASS × SO case, even though in reality, SO will not cover the full CMASS (or DESI) area. Similarly, for the DESI-like × SO case, we adopted a sky coverage of 14 000 square degrees. Thus, our forecasts illustrate the potential of possible future combinations of CMB and galaxy surveys with ideal overlap, rather than specific predictions for the actual combinations of the planned SO and DESI surveys.

We adapted the binning scheme to the noise curves, which are provided up to high multipoles (ℓ_max = 5000 for the lensing and ℓ_max = 7980 for the tSZ). We therefore used multipoles between ℓ = 30 and ℓ = 7980 for the galaxy auto-power spectrum (making 21 logarithmically spaced bins), as well as for the cross correlation with the tSZ (using these high multipoles will require a thorough beam modeling and point source and CIB decontamination), and between ℓ = 30 and ℓ = 5000 for the cross correlation with the CMB lensing. The noise curves for multipoles lower than ℓ = 80 are not provided for the tSZ, and in the range ℓ ∈ [30, 80], we used the noise from the Planck experiment.

We used the same pipeline as in our main study to fit our set of nine parameters. Here, we are primarily interested in β_m, A, and 1 − b_h in the case CMASS × SO, to see how well the constraints will be improved with SO, and we do not expect much improvement on the HOD parameters because in this case, we used the same galaxy survey as in the main study of this paper. In the DESI-like × SO case, however, we are interested in potential constraints on the nine parameters and do expect some improvement on the HOD parameters due to the lower shot noise and larger sky coverage.

Table 2 and Fig. 12 show the constraints that can be achieved when we cross correlate the first stellar mass bin from CMASS and the DESI-like survey with SO over the aforementioned multipole ranges. Figure 12 shows that the estimated parameters recover their fiducial values within 1σ, showing that the results are unbiased.

Fig. 12. Forecast constraints with future experiments. Red: Forecast constraints when the first stellar mass bin from a CMASS survey is cross correlated with full overlap with the future SO CMB experiment, using the noise curves provided by the SO collaboration. Blue: Same, but using the galaxy shot noise and sky coverage of a DESI-like experiment with full overlap with SO (see text for details). The dashed lines show the fiducial parameter values we used for the forecast.

6.2. Results in the CMASS × SO case

In the case of the cross correlation between CMASS and SO, we find that the SO CMB lensing measurement will probe the matter distribution profile and constrain β_m at the level of σ = 0.91. Depending on how strongly β_m is correlated with galaxy stellar mass, it may be possible to detect a trend. Furthermore, when A was fixed to unity, we find σ = 0.66, implying that the constraint will be better if we can assume that galaxies are perfect tracers of dark matter (e.g., no stochasticity).

The expected constraint on 1 − b_h is nearly twice as strong as the constraint obtained with Planck; the improvement arises from the lower tSZ noise and access to much higher multipoles. Measuring this parameter to this precision would test the gas pressure dependence on stellar mass indicated in our study, and hence provide important information on feedback processes. Additionally, these higher multipoles would place constraints not only on the hydrostatic bias (i.e., the pressure profile normalization), but also on the pressure profile. However, as mentioned previously, using these high multipoles will require mitigation of CIB and point source contamination.

As expected, there is no significant improvement on the HOD parameter constraints. The reason is that the HOD parameters are mainly constrained by the galaxy auto-power spectrum because the cross-spectra already constrain other parameters.

6.3. Results in the DESI-like × SO case

In this case, the improvement on all parameters (except α_inc and M_inc) is significant. The parameter β_m entering the matter profile is constrained to σ = 0.30, and even σ = 0.25 when A is fixed to unity. This case study shows the exciting potential of probing the small-scale matter distribution and its dependence on galaxy properties.

The hydrostatic bias constraint gains another factor of two compared to the CMASS × SO case. At this level, it would be possible to study feedback processes with high accuracy; for instance, relating gas thermal energy and mass to host active galactic nuclei (AGN) activity.

The constraints on the galaxy distribution parameters improve by a factor of two when compared to the CMASS × SO case. This comes from the much lower shot noise that enables us to benefit from the smaller scales, as well as from the larger sky coverage.

This forecast shows that using cross correlations from SO and a DESI-like survey with full overlap would lead to improved understanding of the relations between dark matter, gas, and galaxies. It appears that correlating these data will enable researchers to study both the large- and small-scale distributions of all these components, a systematic cosmic census.

7. Conclusion

The advent of multiple wide-area surveys sensitive to galaxy properties, total mass, and gas mass enables a comprehensive study of all components of the cosmic web though cross-correlation studies. We examined the joint distribution of galaxies, dark matter, and gas as a function of galaxy stellar mass. We divided CMASS galaxies into stellar mass bins and cross correlated them with the Planck CMB lensing convergence and Compton-y maps. Making use of the halo model, we fit an HOD and probed the large-scale structure of the Universe as well as the small-scale galaxy distribution through measurement of the galaxy and linear gas bias parameters and the galaxy halo distribution profile. We find in Sect. 5.2.1 that the galaxy bias increases with stellar mass (see Fig. 10) because more massive galaxies preferably reside in higher-mass halos that are most strongly biased. We also find that most massive galaxies are central galaxies (see Fig. 9). Satellite galaxies tend to be less massive, and their radial distribution with their host halo seems to be steeper as their stellar mass increases.

We measured the gas pressure around galaxies and find a hint that the pressure decreases slightly in the vicinity of the most massive galaxies, as shown in Fig. 11. The hydrostatic bias agrees with previous studies, which also reported a slight increase with halo mass (Eckert et al. 2019; Wicker et al. 2023). Finally, cross correlating the Planck CMB lensing convergence map with galaxy maps, we measured a cross-correlation amplitude (lensing parameter A) consistent with unity, implying that galaxies trace the dark matter field well; as a corollary, we found no evidence for stochasticity. However, the noise in the Planck CMB lensing convergence map does not permit any constraint on the matter profile inside halos. In the near future, we forecast in Sect. 6 that next generation CMB experiments will be able to do this by quantifying in particular expectations for SO as well as a DESI-like survey.

Cross-correlation studies such as the one presented here enable a comprehensive study of all components of the cosmic web, their interactions, and their coevolution. This is a powerful new technique with exciting potential to reveal details of large-scale cosmic structure and its emergence.

Acknowledgments

We thank Jean-Baptiste Melin for useful discussions about mass definitions conversions and gas pressure profiles. We acknowledge the use of the python libraries matplotlib (Hunter 2007), numpy (Harris et al. 2020) and scipy (Virtanen et al. 2020). Some of the results in this paper have been derived using the healpy and HEALPix packages (Górski et al. 2005).

References

Appendix A: Radius-dependent halo bias and halo exclusion

For our theoretical model to be accurate on intermediate and small scales, we must account for the radial dependence of halo bias and the effects of halo exclusion. These two effects are dealt with through Q(k|M₁, M₂, z) in Eq. (25). We used the recipe given by van den Bosch et al. (2013) and Cacciato et al. (2013), which we summarize here for completeness.

We write the halo 2 pt correlation function, ξ_hh, in terms of the matter 2 pt function, ξ_mm, as

$$\begin{array}{c}\hfill {\xi }_{\mathit{hh}}(r,z|M_1,M_2)=\{\begin{array}{cc}b(M_1,z)b(M_2,z)\zeta (r,z){\xi }_{\mathit{mm}}(r,z)\hfill & \text{if}r>r_{\min }\hfill \\ -1\hfill & \text{otherwise,}\hfill \end{array}\end{array}$$(A.1)

where ζ(r, z) describes the radial dependence of the halo bias and r_min = max(r_vir(M₁),r_vir(M₂)) ensures the halo exclusion (i.e., that the center of a halo does not reside within another halo). The radial bias function is given by

$$\begin{array}{c}\hfill \zeta (r,z)=\{\begin{array}{cc}{\zeta }_0(r,z)\hfill & \text{if}r\ge r_{\psi }\hfill \\ {\zeta }_0(r_{\psi },z)\hfill & \text{if}r\le r_{\psi },\hfill \end{array}\end{array}$$(A.2)

where r_ψ is defined through

$$\begin{array}{c}\hfill log[{\zeta }_0(r_{\psi },z)/{\xi }_{\mathit{mm}}(r_{\psi },z)]=\psi ,\end{array}$$(A.3)

with ψ = 0.9. The function ζ₀ is fit to N-body simulations by Tinker et al. (2005),

$$\begin{array}{c}\hfill {\zeta }_0(r,z)=\frac{{(1+1.17{\xi }_{\mathit{mm}}(r,z))}^{1.49}}{{(1+0.69{\xi }_{\mathit{mm}}(r,z))}^{2.09}}.\end{array}$$(A.4)

Appendix B: Mass conversion

In this appendix, we derive Eq. (18) by applying the recipe described in Hu & Kravtsov (2003). The NFW profile has the form

$$\begin{array}{c}\hfill \rho (r)=\frac{{\rho }_{\mathrm{s}}}{\left(\frac{r}{r_{\mathrm{s}}}\right){(1+\frac{r}{r_{\mathrm{s}}})}^2}.\end{array}$$(B.1)

Integrating the mass profile until radii r_200m and r_500c gives

$$\begin{array}{cc}\hfill M_{200\mathrm{m}}& =4\pi {\rho }_{\mathrm{s}}{r}_{200\mathrm{m}}^{3}f(r_{\mathrm{s}}/r_{200\mathrm{m}}),\hfill \end{array}$$(B.2)

$$\begin{array}{cc}\hfill M_{500\mathrm{c}}& =4\pi {\rho }_{\mathrm{s}}{r}_{500\mathrm{c}}^{3}f(r_{\mathrm{s}}/r_{500\mathrm{c}}),\hfill \end{array}$$(B.3)

where f(x) = x³(ln(1+x⁻¹)−(1+x)⁻¹). Using the mass definitions, we also have

$$\begin{array}{cc}\hfill M_{200\mathrm{m}}& =\frac{4\pi r_{200\mathrm{m}}^3}{3}{\mathrm{\Delta }}_{200\mathrm{m}}{\rho }_{\mathrm{m}}(z),\hfill \end{array}$$(B.4)

$$\begin{array}{cc}\hfill M_{500\mathrm{c}}& =\frac{4\pi r_{500\mathrm{c}}^3}{3}500{\rho }_{\mathrm{c}}(z)\hfill \end{array}$$(B.5)

$$\begin{array}{cc}& =\frac{4\pi r_{500\mathrm{c}}^3}{3}{\mathrm{\Delta }}_{500}(z){\rho }_{\mathrm{m}}(z),\hfill \end{array}$$(B.6)

with Δ_200m = 200 and ${\mathrm{\Delta }}_{500}(z)=500\frac{{\rho }_{\mathrm{c}}(z)}{{\rho }_{\mathrm{m}}(z)}=\frac{500H{(z)}^2}{{\mathrm{\Omega }}_{\mathrm{m}}{H}_0^2{(1+z)}^3}$.

Using Eqs. (B.2) and (B.3), we have

$$\begin{array}{c}\hfill {\rho }_{\mathrm{s}}=\frac{M_{200\mathrm{m}}}{4\pi r_{200\mathrm{m}}^{3}f(r_{\mathrm{s}}/r_{200\mathrm{m}})}=\frac{M_{500\mathrm{c}}}{4\pi r_{500\mathrm{m}}^{3}f(r_{\mathrm{s}}/r_{500\mathrm{c}})}.\end{array}$$(B.7)

Plugging in the mass definitions, we can relate the two radii

$$\begin{array}{cc}& \frac{{\mathrm{\Delta }}_{200\mathrm{m}}{\rho }_{\mathrm{m}}(z)}{3f(r_{\mathrm{s}}/r_{200\mathrm{m}})}=\frac{{\mathrm{\Delta }}_{500}(z){\rho }_{\mathrm{m}}(z)}{3f(r_{\mathrm{s}}/r_{500\mathrm{c}})},\hfill \end{array}$$(B.8)

$$\begin{array}{cc}& f(r_{\mathrm{s}}/r_{500\mathrm{c}})=\frac{{\mathrm{\Delta }}_{500}(z)}{{\mathrm{\Delta }}_{200\mathrm{m}}}f(r_{\mathrm{s}}/r_{200\mathrm{m}}),\hfill \end{array}$$(B.9)

$$\begin{array}{cc}& \frac{r_{\mathrm{s}}}{r_{500\mathrm{c}}}=x\left(\frac{{\mathrm{\Delta }}_{500}(z)}{{\mathrm{\Delta }}_{200\mathrm{m}}}f(r_s/r_{200\mathrm{m}})\right)=x\left(\frac{{\mathrm{\Delta }}_{500}(z)}{{\mathrm{\Delta }}_{200\mathrm{m}}}f(c_{200\mathrm{m}}^{-1})\right),\hfill \end{array}$$(B.10)

where x is the reciprocal function of f. For completeness, we copy the fitting formula derived by Hu & Kravtsov (2003),

$$\begin{array}{c}\hfill x(f)={(a_1{f}^{2p}+{(3/4)}^2)}^{-1/2}+2f,\end{array}$$(B.11)

with p = a₂ + a₃lnf + a₄(lnf)² and (a₁, a₂, a₃, a₄) = (0.5116, −0.4283, −3.13 × 10⁻³, −3.52 × 10⁻⁵).

Plugging Eq. (B.10) into Eq. (B.6), we obtain

$$\begin{array}{c}\hfill M_{500\mathrm{c}}=\frac{4\pi }{3}{\mathrm{\Delta }}_{500}(z){\rho }_{\mathrm{m}}(z)r_s^{3}x{\left(\frac{{\mathrm{\Delta }}_{500}(z)}{{\mathrm{\Delta }}_{200\mathrm{m}}}f(c_{200\mathrm{m}}^{-1})\right)}^{-3}.\end{array}$$(B.12)

Finally, using r_s = r_200m/c_200m and Eq. (B.4) to express $r_{200\mathrm{m}}={\left(\frac{3M_{200\mathrm{m}}}{4\pi {\mathrm{\Delta }}_{200\mathrm{m}}{\rho }_{\mathrm{m}}(z)}\right)}^{1/3}$, we find

$$\begin{array}{c}\hfill M_{500\mathrm{c}}=M_{200\mathrm{m}}\frac{{\mathrm{\Delta }}_{500}(z)}{{\mathrm{\Delta }}_{200\mathrm{m}}}{c}_{200\mathrm{m}}^{-3}x{\left(\frac{{\mathrm{\Delta }}_{500}(z)}{{\mathrm{\Delta }}_{200\mathrm{m}}}f(c_{200\mathrm{m}}^{-1})\right)}^{-3}.\end{array}$$(B.13)

Appendix C: Compton-y parameter angular power spectrum

In this appendix, we derive the expression for the Compton-y angular power spectrum and its kernel. We first establish a result on the halo distribution that will prove useful in the derivation.

We model the halo number density field, n(M,r), as a Poisson realization, with mean $\overline{n}(M,\mathit{r})$, of the underlying continuous matter overdensity field, δ(r). Taking expectations over just the Poisson ensemble, ⟨ ⋅ ⟩_Poisson, for fixed δ(r), we have

$$\begin{array}{cc}\hfill {⟨n(M,\mathit{r})⟩}_{\mathrm{Poisson}}& =\overline{n}(M,\mathit{r})=\overline{n}(M,z)(1+b(M)\delta (\mathit{r})),\hfill \end{array}$$(C.1)

$$\begin{array}{cc}\hfill {⟨n(M_1,{\mathit{r}}_1)n(M_2,{\mathit{r}}_2)⟩}_{\mathrm{Poisson}}& =\overline{n}(M_1,{\mathit{r}}_1)\overline{n}(M_2,{\mathit{r}}_2)+\overline{n}(M_1,{\mathit{r}}_1){\delta }_D({\mathit{r}}_1-{\mathit{r}}_2){\delta }_D(M_1-M_2),\hfill \end{array}$$(C.2)

where b(M) is the halo bias and $\overline{n}(M,z)$ is the average number density of halos of mass M at redshift z. In this appendix, $\overline{n}(M,z)$ is the halo mass function, which was denoted simply n in the main body of the paper. In the above, δ_D is the Dirac delta function. Taking the full expectation, ⟨⟨⋅⟩⟩ = ⟨⟨⋅⟩_Poisson⟩, over both the Poisson ensemble, using the above equations, and the δ(r)-field ensemble, ⟨ ⋅ ⟩, we find

$$\begin{array}{cc}\hfill ⟨⟨n(M_1,{\mathit{r}}_1)n(M_2,{\mathit{r}}_2)⟩⟩& =⟨\overline{n}(M_1,{\mathit{r}}_1)\overline{n}(M_2,{\mathit{r}}_2)+\overline{n}(M_1,{\mathit{r}}_1){\delta }_D({\mathit{r}}_1-{\mathit{r}}_2){\delta }_D(M_1-M_2)⟩\hfill \end{array}$$(C.3)

$$\begin{array}{cc}\hfill & =⟨\overline{n}(M_1,z_1)[1+b(M_1)\delta ({\mathit{r}}_1)]\overline{n}(M_2,z_2)[1+b(M_2)\delta ({\mathit{r}}_2)]+\overline{n}(M_1,z_1)[1+b(M_1)\delta ({\mathit{r}}_1)]{\delta }_D({\mathit{r}}_1-{\mathit{r}}_2){\delta }_D(M_1-M_2)⟩\hfill \end{array}$$(C.4)

$$\begin{array}{cc}& =\overline{n}(M_1,z_1)\overline{n}(M_2,z_2)+b(M_1)b(M_2)\overline{n}(M_1,z_1)\overline{n}(M_2,z_2)⟨\delta ({\mathit{r}}_1)\delta ({\mathit{r}}_2)⟩+\overline{n}(M_1,z_1){\delta }_D({\mathit{r}}_1-{\mathit{r}}_2){\delta }_D(M_1-M_2).\hfill \end{array}$$(C.5)

The following derivation of the Compton-y angular power spectrum is inspired by Taburet et al. (2011), who adopted the approach introduced by Hernández-Monteagudo et al. (2006). The Compton-y parameter is equal to the line-of-sight integral of the electron pressure,

$$\begin{array}{cc}\hfill y& =\frac{k_{\mathrm{B}}{\sigma }_{\mathrm{T}}}{m_{\mathrm{e}}{c}^2}{\displaystyle \int }dl{n}_{\mathrm{e}}{T}_{\mathrm{e}}\hfill \end{array}$$(C.6)

$$\begin{array}{cc}& ={\displaystyle \int }adr{\displaystyle \sum _{i=1}^N}{y}_{3\mathrm{D},\mathrm{c}}({\mathit{w}}_i)u({\mathit{w}}_i-\mathit{r}),\hfill \end{array}$$(C.7)