© ESO 2021

1. Introduction

The study of large-scale structure (LSS) is a major topic in modern cosmology. The theoretical framework of LSS in the linear regime is well established and has been constrained by multiple observations (see, for example, Dodelson & Schmidt 2020 for a detailed description). On small scales (∼1 Mpc), the growth of structure is driven by the combination of the non-linear gravitational collapse and baryonic processes in the intergalactic gas (van Daalen et al. 2011; Semboloni et al. 2011; Fedelia 2014; Mead et al. 2015). Although the latter is challenging to model, an increasing number of multi-wavelength sky surveys reach high redshifts and high angular resolution (Catinella et al. 2010; Heymans et al. 2012; de Jong et al. 2013; Abbott et al. 2016; Aihara et al. 2018, for example), which extend our understanding of the late-time history of the Universe and make us sensitive to subtle and complicated small-scale physics. In addition, observations of different ‘tracers’ make it possible to study different aspects of LSS. In summary, it is a golden age dominated by surveys that shed light on the role of small-scale physics in LSS formation and evolution.

For many years, cross-correlations between different LSS tracers have been an important tool for helping us understand relations between underlying physics (Hill & Spergel 2014; Van Waerbeke et al. 2014; Kirk et al. 2016; Hojjati et al. 2017; Singh et al. 2017; Ammazzalorso et al. 2020). Compared to other LSS tracers, the galaxy distribution is easier to measure with high precision. Cross-correlations between galaxy positions and LSS has proved to be a powerful tool for studying different properties of LSS. For example, the cross-correlation between galaxy positions and cosmic infrared background (CIB) has been used to probe the properties of dust in star-forming galaxies (Serra et al. 2014); the cross-correlation between galaxy position and 21 cm emission is useful for studying the cosmic reionisation history (Lidz et al. 2008); Kuntz (2015) probes the cross-correlation between galaxy positions and cosmic microwave background (CMB) lensing to study the galaxy bias and lensing amplitude. In this study, we focus on the cross-correlation between the galaxy distribution and the thermal Sunyaev-Zeldovich (tSZ) effect (denoted as the ‘gy’ cross-correlation hereafter) to probe the ‘gas pressure bias’, defined as the multiplication of the mean electron pressure and gas bias, the ratio of the gas overdensity to the mass overdensity, as a proxy of intergalactic gas properties.

The tSZ effect (Zeldovich & Sunyaev 1969; Sunyaev & Zeldovich 1972) is the distortion of the CMB energy spectrum due to inverse Compton scattering by high-energy electrons; tSZ effect is therefore a tracer of the projected intergalactic gas pressure. Since warm to hot intergalactic gas is typically concentrated in galaxy clusters, one can use the tSZ effect to detect galaxy clusters (Planck Collaboration VIII 2011; Hincks et al. 2010). In addition, cross-correlations between tSZ, galaxy clustering, and weak lensing are useful for studying the properties of the diffuse gas as well as the mass distribution of galaxy clusters (Hojjati et al. 2017; Makiya et al. 2018; Koukoufilippas et al. 2020). In order to probe intergalactic hot gas, the tSZ effect has several advantages over the X-ray emission originating from Bremsstrahlung. Firstly, its amplitude, characterized by the Compton y parameter, does not depend on the cluster redshift, while X-ray surface brightness scales with (1 + z)⁻⁴, which makes tSZ sensitive to higher redshifts. Secondly, the y parameter depends linearly on the density of gas particles, while X-ray brightness has a quadratic dependence. The X-ray emission is thus more affected by the clumpiness of gas. In addition, the characteristic frequency dependence of tSZ makes it possible to be fully extracted against other sources of radiation such as CMB, Galactic dust thermal emission, and synchrotron emission (Remazeilles et al. 2011a), while X-ray spectra highly depend on the temperature and composition of sources.

Makiya et al. (2018), Pandey et al. (2019), and Koukoufilippas et al. (2020) report measurements of gy with galaxy data from the 2MASS photometric redshift survey and WISE × SuperCOSMOS, the Dark Energy Survey redMaGiC sample, and the 2MASS redshift survey, respectively. In our study, we use the galaxy sample from the fourth data release of the Kilo-Degree Survey (KiDS; Kuijken et al. 2019) and the Compton y map from the 2015 data release of the Planck mission (Planck Collaboration XXII 2016). The previous studies have larger sky coverage but lower survey depth, while KiDS covers only about 2% of the sky but goes as deep as z ∼ 1. Chiang et al. (2020) has also measured the tSZ signal up to z ∼ 1, but uses a quasar catalogue at high redshift. Quasars might have a strong feedback effect in the cross-correlation, which is difficult to model. In contrast, our measurement uses a pure galaxy sample that goes to the highest redshift to date.

The galaxy distribution is a biased tracer of the mass distribution, and when used in cross-correlation studies, the galaxy bias is degenerate with the bias of the other tracer. For example, in the gy measurements, the galaxy bias and the gas pressure bias are degenerate. In Koukoufilippas et al. (2020), Pandey et al. (2019), and Makiya et al. (2018), this problem was addressed by measuring the galaxy bias from the galaxy auto-correlation function, which requires careful modelling of the auto-correlation noise. For KiDS, the field-to-field depth variation is large, which makes galaxy auto-correlations challenging to model with precision. We therefore adopt an alternative approach, which consists of measuring the galaxy bias using the CMB lensing convergence as the mass proxy, via the ‘gκ’ cross-correlation, which is also adopted in studies such as Ferraro et al. (2015). The gravitational lensing effect of CMB photons (Lewis & Challinor 2006) is an unbiased mass tracer of LSS that has been used to cross-correlate with other tracers to study mass clustering and galaxy bias (Bianchini et al. 2015; Singh et al. 2017; Hurier et al. 2019). In this work, we measure the cross-correlation between galaxy positions and the Planck CMB lensing map (Planck Collaboration VIII 2020) to independently constrain the galaxy bias and to eliminate the need for modelling the galaxy auto-correlation function. It should be noted that the noise modelling of the galaxy distribution also affects the cross-correlations, but only at the covariance level. In this study, we focus on the linear scale properties of the gas and galaxy position while modelling the cross-correlations on non-linear scales with simple one-parameter models. We do not attempt to extract any cosmological information from them. However, with future improvements in data quality, this approach could in principle be generalised to probe non-linear scales.

We note that the galaxy bias could also be constrained from the cross-correlation between foreground galaxy positions and background galaxy shear, known as galaxy-galaxy lensing. However, at high precision, the interpretation of galaxy lensing requires the modelling of non-lensing effects such as the source-lens clustering and the intrinsic alignments (Hamana et al. 2002; Hall & Taylor 2014; Valageas 2013) and shape measurement residual systematics. These are extensively studied in their own right, but in this work we intend to highlight the feasibility of using the galaxy distribution as a proxy for the mass distribution with CMB lensing as the galaxy-mass calibration tool. Although CMB lensing has a much lower signal-to-noise than galaxy lensing for a given sky area, it extends to much higher redshift and is immune to most of the non-lensing effects that can contaminate galaxy lensing.

This paper is structured as follows: In Sect. 2 we describe the theoretical model we use for the cross-correlations; Sect. 3 introduces the dataset that we are using; Sect. 4 presents the method to measure cross-correlations, as well as our estimation of covariance matrix, likelihood, and systematics; Sect. 5 presents the results; Sect. 6 discusses the results and summarises our conclusion. Throughout this study, we assume a flat ΛCDM cosmology with fixed cosmological parameters from Planck Collaboration VI (2020) as our fiducial cosmology: (h, Ω_ch², Ω_bh², σ₈, n_s) = (0.676, 0.119, 0.022, 0.81, 0.967). The impacts of fixing cosmological parameters are discussed in Sect. 6.

2. Models

We measure the angular cross-correlation in harmonic space. In general, the angular cross-correlation between two projected tracers, u and v, at scales ℓ ≳ 10 are well computed by the Limber approximation (Limber 1953; Kaiser 1992):

$$\begin{array}{c}\hfill C_{\ell }^{\mathit{uv}}={\displaystyle {\int }_0^{{\chi }_{\mathrm{H}}}\frac{\mathrm{d}\chi }{{\chi }^2}{W}^u(\chi )W^v(\chi )P_{\mathit{UV}}(k=\frac{\ell +1/2}{\chi },z(\chi )),}\end{array}$$(1)

where χ denotes the comoving distance; χ_H is the comoving distance to the horizon; W^u(χ) is the radial kernel of tracer u; P_UV(k, z) is the three-dimensional (3D) cross-power spectrum of associated 3D fluctuating tracers U and V:

$$\begin{array}{c}\hfill ⟨{\delta }_U(\mathit{k}){\delta }_V(\mathit{k}\prime )⟩={(2\pi )}^3\delta (\mathit{k}+\mathit{k}\prime )P_{\mathit{UV}}(k).\end{array}$$(2)

In this work, the fluctuating physical quantity for tSZ is the 3D electron pressure fluctuation ΔP_e; for galaxy number counts it is the 3D galaxy overdensity δ_G; for CMB lensing it is the 3D mass overdensity δ_M. We note that, throughout this paper, the two-dimensional (2D) projected tracers (projected galaxy number, Compton y and lensing convergence κ) are labelled as lower case letters g, y, and κ; while the corresponding 3D tracer (galaxy number distribution, electron pressure, and mass distribution) are labelled as capital letters G, P, and M.

In this work we measure $C_{\ell }^{\mathrm{g}\mathit{y}}$ and $C_{\ell }^{\mathrm{g}\kappa }$. The angular fluctuation of galaxy number density Δ_g is the 2D projection of 3D number-density fluctuations:

$$\begin{array}{c}\hfill {\mathrm{\Delta }}_{\mathrm{g}}(\widehat{\mathit{\theta }})={\displaystyle {\int }_0^{{\chi }_{\mathrm{H}}}\mathrm{d}\chi \frac{H(z)}{c}{n}_{\mathrm{g}}(z){\delta }_{\mathrm{G}}(\chi (z)\widehat{\mathit{\theta }},z),}\end{array}$$(3)

where c is the speed of light; H(z) is the Hubble constant at redshift z; ${\delta }_{\mathrm{G}}(\chi (z)\widehat{\mathit{\theta }},z)$ is the 3D galaxy number density fluctuation, and n_g(z) is the normalised redshift distribution of galaxies, which depends on the sky survey. At large scales, we model the number density fluctuations so that it is proportional to the underlying mass fluctuation: δ_G = ⟨b_g⟩δ_M where ⟨b_g⟩ is the mean galaxy bias of the galaxy population and δ_M is the total mass overdensity. The galaxy kernel is given as

$$\begin{array}{c}\hfill W^{\mathrm{g}}(\chi )=\frac{H(\chi )}{c}{n}_{\mathrm{g}}[z(\chi )].\end{array}$$(4)

The tSZ signal is parametrised by the Compton-y parameter, given by:

$$\begin{array}{c}\hfill y(\widehat{\mathit{\theta }})=\frac{{\sigma }_{\mathrm{T}}}{m_{\mathrm{e}}{c}^2}{\displaystyle {\int }_0^{{\chi }_{\mathrm{H}}}\frac{\mathrm{d}\chi }{1+z}{P}_{\mathrm{e}}(\chi \widehat{\mathit{\theta }}),}\end{array}$$(5)

where σ_T is the Thomson scattering cross-section; and m_e is the electron mass; P_e is the gas electron pressure. Like the galaxy number density fluctuation, we also model the intergalactic gas overdensity as a linearly biased tracer of the underlying mass fluctuation at large scales (Goldberg & Spergel 1999; Van Waerbeke et al. 2014), so the pressure fluctuation ΔP_e ≡ P_e − ⟨P_e⟩=⟨P_e⟩δ_gas = ⟨b_yP_e⟩δ_m, where δ_gas denotes the gas overdensity; b_y is the gas bias; and ⟨P_e⟩≡⟨n_ek_BT_e⟩ is the mean electron pressure in gas halos. The combination ‘gas pressure bias’ ⟨b_yP_e⟩ is the mean pressure weighted by gas bias, which is related to the thermodynamics of gas inside halos. The tSZ kernel is given by

$$\begin{array}{c}\hfill W^y(\chi )=\frac{{\sigma }_{\mathrm{T}}}{m_{\mathrm{e}}{c}^2}\frac{1}{1+z(\chi )}·\end{array}$$(6)

The CMB lensing kernel is the mass fluctuation convolved with the lensing kernel:

$$\begin{array}{c}\hfill W^{\kappa }(\chi )=\frac{3H_0^2{\mathrm{\Omega }}_{\mathrm{m}}}{2a{c}^2}\chi \frac{{\chi }_{\mathrm{CMB}}-\chi }{{\chi }_{\mathrm{CMB}}},\end{array}$$(7)

where χ_CMB is the comoving distance to the last-scattering surface at z ∼ 1100; a denotes the scale factor. Therefore, in the linear regime, all three tracers are modelled as linearly biased mass fluctuation convolved with respective kernels, so the linear cross-power spectrum P_UV(k) is the linearly biased matter power spectrum:

$$\begin{array}{c}\hfill \begin{array}{cc}& P_{\mathrm{GP}}^{\mathrm{lin}}(k)=⟨b_{\mathrm{g}}⟩⟨b_y{P}_{\mathrm{e}}⟩P^{\mathrm{lin}}(k)\hfill \\ \hfill & P_{\mathrm{GM}}^{\mathrm{lin}}(k)=⟨b_{\mathrm{g}}⟩P^{\mathrm{lin}}(k),\hfill \end{array}\end{array}$$(8)

where P^lin(k) is the linear matter power spectrum. For simplicity, hereafter we omit the redshift dependence in the notation of ⟨b_g⟩ and ⟨b_yP_e⟩.

Following the method given in Hang et al. (2021), to account for non-linear effects at small scales, we model the non-linear portion of our power spectra as an unknown amplitude multiplied by a physical model template:

$$\begin{array}{c}\hfill \begin{array}{cc}& P_{\mathrm{GP}}^{\mathrm{nl}}(k)=⟨c_{\mathrm{g}y}⟩T_{\mathrm{GP}}^{\mathrm{nl}}(k),\hfill \\ \hfill & P_{\mathrm{GM}}^{\mathrm{nl}}(k)=⟨c_{\mathrm{g}\kappa }⟩T_{\mathrm{GM}}^{\mathrm{nl}}(k),\hfill \end{array}\end{array}$$(9)

where T^nl(k) is the additional non-linear templates. ⟨c_gy⟩ and ⟨c_gκ⟩ are two re-scaling parameters that account for differences between the amplitudes of non-linear gy, gκ, and dark matter cross-correlations. The total power spectrum is then modelled as:

$$\begin{array}{c}\hfill \begin{array}{cc}& P_{\mathrm{GP}}(k)=⟨b_{\mathrm{g}}⟩⟨b_y{P}_{\mathrm{e}}⟩P^{\mathrm{lin}}(k)+⟨c_{\mathrm{g}y}⟩T_{\mathrm{GP}}^{\mathrm{nl}}(k),\hfill \\ \hfill & P_{\mathrm{GM}}(k)=⟨b_{\mathrm{g}}⟩P^{\mathrm{lin}}(k)+⟨c_{\mathrm{g}\kappa }⟩T_{\mathrm{GM}}^{\mathrm{nl}}(k).\hfill \end{array}\end{array}$$(10)

In the standard halo model for non-linear matter power spectrum, the two-halo (linear part) and one-halo (non-linear part) terms both depend on profile amplitudes, so their amplitudes should be correlated. However, there might be some uncertainties in the modelling that might break this correlation, for example, the uncertainty in halo mass function. In addition, Koukoufilippas et al. (2020) points out that the profiles of tracers in Fourier space may not be independent, so the authors introduce a free parameter ρ_yg for the 1-halo term to account for it. In our study, c_gy and c_gκ account for combinations of such uncertainties to the first order.

To ensure that our constraints on linear bias are robust to the exact non-linear models, we try three well-used models for the non-linear power spectrum templates as well as trying a purely linear model:

(1) The first is the halo model: the non-linear power spectrum is given by the one-halo term of the halo model (Cooray & Sheth 2002; Seljak 2000):

$$\begin{array}{c}\hfill T_{\mathit{UV}}^{\mathrm{nl}}(k)=P_{\mathit{UV}}^{1\mathrm{h}}(k)\equiv {\displaystyle {\int }_0^{\infty }\mathrm{d}M\frac{\mathrm{d}n}{\mathrm{d}M}{p}_U(k\mid M)p_V(k\mid M),}\end{array}$$(11)

where dn/dM is the halo mass function and p_U(k ∣ M) is the profile of the tracer U with mass M in Fourier space:

$$\begin{array}{c}\hfill p_U(k\mid M)\equiv 4\pi {\displaystyle {\int }_0^{\infty }\mathrm{d}r{r}^2\frac{sin(kr)}{\mathit{kr}}{p}_U(r\mid M).}\end{array}$$(12)

For CMB lensing, we take the profile that enters the halo model as the dark matter halo profile, which is typically modelled via the Navarro-Frenk-White profile (NFW; Navarro et al. 1996):

$$\begin{array}{c}\hfill p_{\mathrm{M}}(r\mid M)={\rho }_{\mathrm{NFW}}(r\mid M)\propto \frac{1}{r/r_{\mathrm{s}}{(1+r/r_{\mathrm{s}})}^2},\end{array}$$(13)

where r_s is the characteristic radius of a dark matter halo, which relates to the halo mass that we take from the mass-concentration relation.

The galaxy population in a halo is divided into centrals and satellites, the abundances of which we relate to halo mass via a Halo Occupation Distribution (HOD) model (Zheng et al. 2005; Peacock & Smith 2000):

$$\begin{array}{c}\hfill \begin{array}{cc}& N_{\mathrm{c}}(M)=\frac{1}{2}[1+\mathrm{erf}\left(\frac{log(M/M_{min})}{{\sigma }_{\ln \mathrm{M}}}\right)]\hfill \\ \hfill & N_{\mathrm{s}}(M)=N_{\mathrm{c}}(M)\mathrm{\Theta }(M-M_0){\left(\frac{M-M_0}{M_1}\right)}^{{\alpha }_{\mathrm{s}}},\hfill \end{array}\end{array}$$(14)

where N_c(M) and N_s(M) are the mean number of central and satellite galaxies respectively. M₁, M₀, M_min, σ_M, α_s are free parameters in principle. The galaxy density profile is then:

$$\begin{array}{c}\hfill p_{\mathrm{g}}(k\mid M)={\overline{n}}_{\mathrm{g}}^{-1}[N_{\mathrm{c}}(M)+N_{\mathrm{s}}(M)p_{\mathrm{s}}(k\mid M)],\end{array}$$(15)

where ${\overline{n}}_{\mathrm{g}}$ is the mean galaxy number density. We assume that central galaxies exist at the halo centre while satellites follow the underlying matter distribution, so the satellite profile p_s is the NFW profile. In this work, we fix HOD parameters to {σ_M, α_s, log₁₀M₁, log₁₀M₀, log₁₀M_min}={0.15, 1, 13, 11.86, 11.68} as constrained from Zheng et al. (2005). Here masses are in the unit of h⁻¹ M_⊙. While this may not be a correct description of our galaxy population, we see later that our conclusions are unaffected by the details of our non-linear model.

The y signal derives via the electron pressure profile p_e(r, M, z), which we take from Arnaud et al. (2010):

$$\begin{array}{cc}\hfill p_{\mathrm{e}}(r,M,z)=& 1.65{(h/0.7)}^2\mathrm{eV}{\mathrm{cm}}^{-3}\hfill \\ \hfill & \times E^{8/3}(z){\left[\frac{M}{3\times {10}^{14}(0.7/h)M_{\odot }}\right]}^{2/3+{\alpha }_{\mathrm{p}}}p(x),\hfill \end{array}$$(16)

where x ≡ r/r₅₀₀ and E(z)≡H(z)/H₀. r₅₀₀ is the radius that encloses a region with average density equal to 500 times the critical density of the Universe. The parameter α_p = 0.12 as given by Arnaud et al. (2010). The self-similar part of the pressure profile p(x) is given by (Nagai et al. 2007):

$$\begin{array}{c}\hfill p(x)\equiv \frac{P_0{(0.7/h)}^{3/2}}{{\left(c_{500}x\right)}^{\gamma }{[1+{\left(c_{500}x\right)}^{\alpha }]}^{(\beta -\gamma )/\alpha }}·\end{array}$$(17)

The parameters in p(x) are taken as the best-fitted values from Planck Collaboration V (2013): {P₀, c₅₀₀, α, β, γ}={6.41, 1.81, 1.33, 4.13, 0.31}.

We note that, in halo model, the gas pressure bias can be expressed as:

$$\begin{array}{c}\hfill ⟨b_y{P}_{\mathrm{e}}⟩(z)={\displaystyle {\int }_0^{\infty }\mathrm{d}M\frac{\mathrm{d}n}{\mathrm{d}M}(z)b_{\mathrm{h}}(M,z){\int }_0^{\infty }\mathrm{d}r4\pi r^2{p}_{\mathrm{e}}(r,M,z),}\end{array}$$(18)

where

$$\begin{array}{c}\hfill {\displaystyle {\int }_0^{\infty }\mathrm{d}r4\pi r^2{p}_{\mathrm{e}}(r,M,z)\equiv E_{\mathrm{T}}(M,z),}\end{array}$$(19)

which means that

$$\begin{array}{c}\hfill ⟨b_y{P}_{\mathrm{e}}⟩(z)={\displaystyle {\int }_0^{\infty }\mathrm{d}M\frac{\mathrm{d}n}{\mathrm{d}M}(z)b_{\mathrm{h}}(M,z)E_{\mathrm{T}}(M,z).}\end{array}$$(20)

Therefore, the gas pressure bias directly links to the thermal energy of dark matter halos.

We took the halo mass function and halo bias needed in the halo model from the fitting formulae of Tinker et al. (2008, 2010), respectively. We chose this halo model as our fiducial model.

(2) HALOFIT non-linear model: we isolate the purely non-linear part of the HALOFIT model (Smith et al. 2003; Takahashi et al. 2012) by taking the full model and subtracting linear theory:

$$\begin{array}{c}\hfill T^{\mathrm{nl}}(k)=P^{\mathrm{HF}}(k)-P^{\mathrm{lin}}(k),\end{array}$$(21)

where P^HF(k) is the HALOFIT matter power spectrum. While the HALOFIT model was calibrated to the matter power spectrum only, we hope that the non-linear shape is general enough to capture the correct shape of the galaxy–Pressure and galaxy–matter power spectra that are relevant to our cross-correlations. Any amplitude differences will be absorbed by our multiplicative non-linear coefficients. This non-linear template is also used in Hang et al. (2021).

(3) Constant non-linear model: The non-linear power spectra are constants:

$$\begin{array}{c}\hfill T_{\mathit{UV}}^{\mathrm{nl}}(k)=P_{\mathit{UV}}^{1\mathrm{h}}(0),\end{array}$$(22)

where $P_{\mathit{UV}}^{1\mathrm{h}}(0)$ is the one-halo term of the halo-model power spectrum at k = 0 Mpc⁻¹. This model should only work on large scales where the one-halo region can be treated as a point source.

To make sure that our non-linear models are not sensitive to the precise shape of the non-linear power spectrum, we only fit C_ℓ’s in each tomographic bin within an angular scale corresponding to k < 0.7 Mpc⁻¹ via the Limber approximation. In addition, Mead & Verde (2021) points out that the halo model is not accurate in the transition between the one- and two-halo regions. To attempt to correct for this, we follow Koukoufilippas et al. (2020) and multiply the power spectra given by (1) and (3) with a scale-dependent quantity:

$$\begin{array}{c}\hfill R(k)\equiv \frac{P^{\mathrm{HF}}(k)}{P_{\mathrm{MM}}^{\mathrm{hm}}(k)},\end{array}$$(23)

where P^HF(k) is from HALOFIT and $P_{\mathrm{MM}}^{\mathrm{hm}}(k)$ is the matter–matter power evaluated via the halo model (with NFW profiles). Although Mead & Verde (2021) showed that the correction required to the halo model in the transition region is not universal, and instead depends on the tracers being modelled, it was shown that attempting some correction is better than no correction at all.

Previous studies on galaxy or gas cross-correlations such as Van Waerbeke et al. (2014), Bianchini et al. (2015), Kuntz (2015) treat the galaxy or gas distribution to be proportional to mass distribution on all scales, which means they only have two free parameters: ⟨b_g⟩ or ⟨b_yP_e⟩. However, it should be noted that this model is not physically accurate because the galaxy and gas distributions have significantly different non-linear details compared to the matter distribution. Pandey et al. (2019) use a linear bias model for the galaxy distribution over the scales R ≳ 10 Mpc, to which our measurements are not sensitive. However, as is indicated in Sugiyama et al. (2020), to apply a fully linear model, one needs a scale cut of at least 12 Mpc. The free prefactors for our three models account for different non-linear amplitudes, but our models are probably still not accurate for the details of the non-linear shape. To model the shapes more accurately, one would need to constrain HOD and pressure profile parameters as well as to significantly improve the treatment of the transition region, which is not feasible in this analysis given the noise level of current data. We leave this to future study. In summary, with fixed cosmological parameters, by measuring $C_{\ell }^{\mathrm{g}\mathit{y}}$ and $C_{\ell }^{\mathrm{g}\kappa }$ we can independently constrain and compare ⟨b_g⟩, ⟨b_yP_e⟩, ⟨c_gy⟩ and ⟨c_gκ⟩ with the three non-linear models introduced above.

3. Data

3.1. KiDS Data

We use the lensing catalogue provided by the fourth data release of the KiDS (Kuijken et al. 2019) as our galaxy sample. KiDS is a sky-survey project, which measures the positions and shapes of galaxies using the VLT Survey Telescope (VST) at the European Southern Observatory (ESO). It is primarily designed for weak-lensing applications. The footprint of KiDS DR4 (also called KiDS-1000) is divided into a northern and southern patch, with total coverage of 1006 deg² of the sky (corresponding to a fraction of f_sky = 2.2%). The footprint is shown in the upper panel of Fig. 1. High-quality images are produced with VST-OmegaCAM. Combining with the VISTA Kilo-degree INfrared Galaxy survey (VIKING; Edge et al. 2013), the observed galaxies are photometrically measured in nine optical and near-infrared bands ugriZYJHK_s. The KiDS survey covers redshifts z ≲ 1.5, which makes it a useful dataset to trace the history of different components of the LSS into the early Universe. For each galaxy in the lensing catalogue, the ellipticities are measured with the LENSFIT algorithm (Miller et al. 2013). We only use the ‘gold subsample’ (Wright et al. 2020) of the lensing catalogue since the redshift distribution is more accurately calibrated in this subsample. We present the information of the galaxy sample that we use in Table 1. Note, however, that we do not use the shape information in our fiducial analysis. In Appendix A we use the shape information to replace CMB lensing as an alternative measurement and sanity check.

Fig. 1. KiDS-1000 footprint and masks for the Plancky map and CMB lensing map. Regions in purple are masked out.

We perform a tomographic measurement of cross-correlations by dividing the galaxy catalogue into five redshift bins according to the best-fit photometric redshift z_B of each galaxy. These are the same redshift bins used in the KiDS-1000 cosmology papers (Asgari et al. 2021; Heymans et al. 2021; Tröster et al. 2021). The redshift distribution of each bin is calibrated using Self-Organising Maps (SOM) as described in Wright et al. (2020), Hildebrandt et al. (2021). We note that the SOM-calibrated redshift distributions in this study are not exactly the same as Hildebrandt et al. (2021) in which the redshift distributions are calibrated with a galaxy sample weighted by the lensfit weight, while in this work the redshift distributions are calibrated with the raw, unweighted sample. The redshift distributions of the 5 tomographic bins are shown in Fig. 2. Galaxy overdensity maps are produced for each tomographic bin in the HEALPIX (Gorski et al. 2005) format with NSIDE = 1024, corresponding to a pixel size of 3.4 arcmin. For each tomographic bin, the galaxy overdensity in each pixel is given as

$$\begin{array}{c}\hfill {\mathrm{\Delta }}_{\mathrm{g},i}=\frac{n_i-\overline{n}}{\overline{n}},\end{array}$$(24)

Fig. 2. Redshift distributions of the five tomographic bins of the KiDS gold galaxy sample.

where i denotes the pixel index, n_i is the number of galaxies in the ith pixel and $\overline{n}$ is the average galaxy number of all the pixels in the footprint and the given redshift bin. The galaxy mask for the cross-correlation measurement is just the KiDS footprint, which is presented in the upper panel of Fig. 1.

3.2. tSZ data

We use the all-sky Compton-y-pagination map presented in Planck Collaboration XXII (2016) for the tSZ signal. The y map that we used is constructed with the Modified Internal Linear Combination Algorithm (MILCA; Hurier et al. 2013), which properly suppresses scale-dependent contamination, and projects out the CMB signal. Planck has also released another y map constructed with the Needlet Internal Linear Combination (NILC; Basak & Delabrouille 2012) method. Although the MILCA and NILC y maps agree with each other in most of the relevant studies, the NILC map turns out to be noisier (Planck Collaboration XXII 2016). Therefore, we apply the Planck MILCA map in this study and take the NILC y map as a consistency check. Both y maps have a beam full width at half Maximum (FWHM) of 10 arcmin.

Before calculating the gy cross-correlation, we mask out the Milky Way and point sources with a joint mask of the Planck 60% Galactic mask and point source mask. The combined mask is shown in the middle panel of Fig. 1. The mask and y map are originally provided in the HEALPIX format with NSIDE = 2048 and we degrade them into NSIDE = 1024 to match the resolution of the KiDS galaxy overdensity map.

To evaluate the CIB contamination in the galaxy-tSZ cross-correlation, we also introduced the Planck 545 GHz CIB intensity map as a CIB template (Planck Collaboration XLVIII 2016). The CIB intensity map was generated with the generalised-NILC (Remazeilles et al. 2011b, GNILC) method and has an angular resolution of 5 arcmin. We first convolved the CIB map with a $\sqrt{{10}^2-5^2}=8.66$ arcmin Gaussian filter to match its resolution to the 10 arcmin of the y map before degrading the map to NSIDE = 1024.

3.3. CMB lensing data

We used the tSZ-deprojected Planck CMB lensing map from the 2018 Data Release (Planck Collaboration VIII 2020) to measure the galaxy-CMB Lensing cross-correlation. The map is provided in the format of the spherical harmonic transformation of the lensing convergence κ_ℓm, which is related to the lensing potential ϕ via

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

There might be reconstruction bias in the CMB lensing map, which is typically descriped by a lensing amplitude parameter, A_L. Planck Collaboration VIII (2020) reports a value of A_L that is slightly greater than 1, with a significance of ∼2σ. Efstathiou & Gratton (2020) points out that this might be due to the fluctuations in the temperature power spectrum at high ℓ. In gκ cross-correlation, A_L degenerates with galaxy bias. This analysis cannot break the degeneracy, so we take the prior information, A_L = 1.

We first transformed κ_ℓm back into a HEALPIXκ map with NSIDE = 1024. The corresponding mask is provided along with the CMB lensing data. It is shown in the lower panel of Fig. 1.

4. Measurements

4.1. Cross-correlation measurements

The cross-correlation between two sky maps, that are smoothed with the beam window function b_beam(ℓ), is related to the real C_ℓ with

$$\begin{array}{c}\hfill {\widehat{C}}_{\ell }^{\mathit{uv}}=C_{\ell }^{\mathit{uv}}{b}_{\mathrm{beam}}^u(\ell )b_{\mathrm{beam}}^v(\ell )b_{\mathrm{pix}}^u(\ell )b_{\mathrm{pix}}^v(\ell ),\end{array}$$(26)

where ${\widehat{C}}_{\ell }^{u\mathit{v}}$ denotes the smoothed C_ℓ between sky map u and v; and b_pix(ℓ) is the pixelisation window function. In our analysis we take the Gaussian window function which is given by

$$\begin{array}{c}\hfill b_{\mathrm{beam}}(\ell )=exp(-\ell (\ell +1){\sigma }^2/2),\end{array}$$(27)

where $\sigma =\mathit{FWHM}/\sqrt{8ln2}$. The pixelisation window function corresponding to NSIDE = 1024 is provided by the HEALPIX package. We note that for the KiDS galaxy map, FWHM = 0.

We use POLSPICE to estimate the angular cross-power spectra. Mode-coupling due to mask and beam smoothing are corrected during this process. Fourier ringings are reduced by setting the internal parameters of POLSPICE apodizesigma = 60 and thetamax = 60 deg. The measured angular power spectra are binned into 10 linear bins from ℓ = 100 to ℓ = 1100. The high limit of ℓ corresponds to the Planck beam, which has a size of 10 arcmin.

4.2. Covariance matrix

We combined two methods to estimate the covariance matrix in our C_ℓ measurement: One is jackknife resampling, the other is an analytical Gaussian covariance matrix. For jackknife resampling, we generated 415 jackknife samples by masking out pixels corresponding to NSIDE = 32 (which has a size of 1.83°) in turn from the KiDS galaxy overdensity map. Smaller jackknife pixels would fail to estimate the variance at large scales, while with larger jackknife pixels we would not have enough realisations. So we chose this intermediate jackknife pixel size to balance the pixel size with the sample size. However, as shown in Fig. 3, the jackknife pixel with NSIDE = 16 (corresponding to an angular size of 3.66°) gives a consistent standard deviation. The cross-correlations $C_{\ell }^{\mathrm{g}\mathit{y}}$ and $C_{\ell }^{\mathrm{g}\kappa }$ are measured with each of these jackknife samples and the covariance matrix is calculated with

$$\begin{array}{c}\hfill {\mathrm{Cov}}^{\mathrm{JK}}(C_{\ell }^{\mathit{uv}},C_{{\ell }^{\prime }}^{\mathit{wz}})=\frac{N_{\mathrm{JK}}-1}{N_{\mathrm{JK}}}{\displaystyle \sum _{n=1}^{N_{\mathrm{JK}}}\mathrm{\Delta }{C}_{\ell }^{uv,(n)}\mathrm{\Delta }{C}_{{\ell }^{\prime }}^{wz,(n)},}\end{array}$$(28)

Fig. 3. Standard deviation calculated from the diagonal of ${\mathrm{Cov}}^{\mathrm{JK}}(C_{\ell }^{\mathrm{g}\mathit{y}},C_{\ell }^{\mathrm{g}\mathit{y}})$ and ${\mathrm{Cov}}^{\mathrm{ANA}}(C_{\ell }^{\mathrm{g}\mathit{y}},C_{\ell }^{\mathrm{g}\mathit{y}})$ in the third tomographic bin. We also plot the standard deviation from the jackknife covariance matrix with Nside = 16 as a consistency check.

where N_JK = 415 denotes the number of jackknife samples; uv, wz ∈ {gy, gκ}; $\mathrm{\Delta }{C}_{\ell }^{u\mathit{v},(n)}$ is the difference between the cross-correlation of the nth jackknife sample and the mean cross-correlation over all samples.

Since the KiDS footprint is only ∼2% of the sky, it is hard to generate enough jackknife samples to fully recover the true covariance matrix; as such, the off-diagonal components of the jackknife covariance matrix are noisy. In addition, since different jackknife regions have slightly different shapes, we could not recover the mode-coupling in the covariance matrix associated with the whole map geometry. To better estimate the off-diagonal components and account for mode-coupling accurately, we also estimated the covariance matrix using an analytical method.

The main contribution to the covariance matrix is from a Gaussian random field:

$$\begin{array}{c}\hfill {\mathrm{Cov}}^{\mathrm{G}}(C_{\ell }^{\mathit{uv}},C_{{\ell }^{\prime }}^{\mathit{wz}})={\delta }_{\ell {\ell }^{\prime }}\frac{C_{\ell }^{\mathit{uw}}{C}_{{\ell }^{\prime }}^{\mathit{vz}}+C_{\ell }^{\mathit{uz}}{C}_{{\ell }^{\prime }}^{\mathit{vw}}}{f_{\mathrm{sky}}(2\ell +1)},\end{array}$$(29)

where f_sky = 2.2% is the sky fraction. Sky masks introduce non-zero coupling between different ℓ. To account for this, we used the method given by Efstathiou (2004) and García-García et al. (2019) and implemented in the NAMASTER package (Alonso et al. 2019)¹. The auto-power spectra in Eq. (29) are directly measured from maps so that noise auto-spectra can also be included; the cross-power spectra are instead calculated from the theoretical model described in Sect. 2, since their measurements are significantly noisier.

The non-Gaussian term includes a connected contribution resulting from the small-scale non-linear clustering of the tracers, related to the trispectrum of the tracers. According to Koukoufilippas et al. (2020), Barreira et al. (2018), and Nicola et al. (2020), this contribution is only significant for low redshifts z ≲ 0.2; therefore, we neglect it in our covariance matrix. Another non-Gaussian contribution is the super-sample covariance (Takada & Hu 2013, SSC) resulting from mode mixing between observed in-survey and the unobserved out-of-survey modes, we also ignore this in this work.

To calculate the analytical covariance matrices, we need to use model parameters that we do not know a priori. So we follow Koukoufilippas et al. (2020) and take a two-step fitting: we first take an assumed value for the model parameters to calculate these covariance matrices and then use these to find the best-fit parameters. We then update the covariance matrix, using these best-fit parameters, and fit for the parameters again. The best-fit parameters from this second round of fitting are taken to be our fiducial results.

The diagonal components of the jackknife and analytical covariance matrices generally agree with each other (see Fig. 3 as an example), this justifies that we can ignore the non-Gaussian contribution in the covariance matrix. To ensure that we recover realistic error bar sizes, we combine the variance estimated from the jackknife covariance matrix with the analytical covariance matrix, as in Koukoufilippas et al. (2020), to account for the coupling between different modes caused by masks and non-Gaussianities while avoiding the statistical noise in the jackknife estimator. Therefore our final covariance matrix is:

$$\begin{array}{c}\hfill {\mathrm{Cov}}_{\mathit{ij}}={\mathrm{Cov}}_{\mathit{ij}}^{\mathrm{ANA}}\sqrt{\frac{{\mathrm{Cov}}_{\mathit{ii}}^{\mathrm{JK}}{\mathrm{Cov}}_{\mathit{jj}}^{\mathrm{JK}}}{{\mathrm{Cov}}_{\mathit{ii}}^{\mathrm{ANA}}{\mathrm{Cov}}_{\mathit{jj}}^{\mathrm{ANA}}}}·\end{array}$$(30)

The correlation coefficient matrices of all the tomographic bins are shown in Fig. 4.

Fig. 4. Full correlation coefficient matrices of covariance matrices (defined in Eq. (30)) of the gy and the gκ cross-correlations in each tomographic bin. Each covariance matrix consists of four sub-matrices, corresponding to the covariances of $C_{\ell }^{\mathrm{g}\mathit{y}}$ and $C_{\ell }^{\mathrm{g}\kappa }$ (block diagonals), as well as their cross-covariance. For each sub-matrix, the pixels show the cross coefficient between binned Cℓ’s, where ℓ bins are defined in Sect. 4.1.

4.3. Likelihood

Since we are working with a wide ℓ range, there are many degrees of freedom in each ℓ bin. According to the central limit theorem, the bin-averaged C_ℓ’s obey a Gaussian distribution around their true values. Thus we assume that the measured power spectra follow a Gaussian likelihood:

$$\begin{array}{c}\hfill -2lnL(\mathit{D}\mid \mathit{q})={\chi }^2\equiv {(\mathit{D}-\mathit{M}(\mathit{q}))}^{\mathrm{T}}{\mathrm{Cov}}^{-1}(\mathit{D}-\mathit{M}(\mathit{q})),\end{array}$$(31)

where q ≡ {⟨b_g⟩ × ⟨b_yP_e⟩,⟨b_g⟩, ⟨c_gy⟩, ⟨c_gκ⟩} stands for our model parameters (for the ‘linear model’ we only have two parameters q ≡ {⟨b_g⟩ × ⟨b_yP_e⟩,⟨b_g⟩}); the data vector $\mathit{D}\equiv (C_{\ell }^{\mathrm{g}\mathit{y}},C_{\ell }^{\mathrm{g}\kappa })$ is a concatenation of measured gy and gκ cross-correlations; M(q) is the cross-correlation predicted by model described in Sect. 2 with parameter q. The likelihoods are calculated separately in each tomographic bin. We note that the gas pressure bias ⟨b_yP_e⟩, which we are primarily interested in, is the ratio between the first two model parameters.

We sample the posterior distribution of model parameters using the Markov chain Monte Carlo method (MCMC) using the emcee package (Foreman-Mackey et al. 2013). We take flat priors for all four model parameters:

$$\begin{array}{c}\hfill \begin{array}{cc}& 0\le ⟨b_{\mathrm{g}}⟩\le 3,\hfill \\ \hfill & 0\le ⟨b_{\mathrm{g}}⟩\times ⟨b_y{P}_{\mathrm{e}}⟩\le 9,\hfill \\ \hfill & 0\le ⟨c_{\mathrm{g}\kappa }⟩\le 10,\hfill \\ \hfill & 0\le ⟨c_{\mathrm{g}y}⟩\le 10.\hfill \\ \hfill \end{array}\end{array}$$(32)

The lower boundaries of the prior is a physical constrain, indicating that the parameters cannot be negative; the upper boundaries are set so that at least 5σ of the marginalised posterior distribution falls in these ranges.

We are not doing a joint analysis of all the tomographic bins, so the correlations between different redshift bins are not relevant to our analysis. Following Crocce et al. (2015), Koukoufilippas et al. (2020), and Hang et al. (2021), our model parameters are fit independently for each of the tomographic bins. The theoretical model is calculated using the Core Cosmology Library package (Chisari et al. 2019).

4.4. Systematics

4.4.1. CIB contamination

The CIB radiation is the accumulated emission from early galaxy populations spanning a large range of redshifts, mostly generated from dust thermal radiation around extragalactic star-formation regions (Hauser & Dwek 2001). The tSZ map is contaminated by residual CIB (Hurier 2015; Yan et al. 2019), which dominates extragalactic signals at high frequency and high redshifts. This residual might contaminate our galaxy-tSZ cross-correlation. We follow the method in Koukoufilippas et al. (2020) to model the CIB contamination in the y map as a factor α_CIB times a CIB template map, which is taken to be the Planck CIB intensity map at 545 GHz (Planck Collaboration XLVIII 2016):

$$\begin{array}{c}\hfill \widehat{y}(\mathit{\theta })=y(\mathit{\theta })+{\alpha }_{\mathrm{CIB}}{I}_{\mathrm{CIB}}(\mathit{\theta }),\end{array}$$(33)

where $\widehat{\mathit{y}}$ denotes the contaminated y map and I_CIB denotes the CIB intensity in 545 GHz. So the measured galaxy-tSZ cross-correlation is given by:

$$\begin{array}{c}\hfill C_{\ell }^{\mathrm{g}\widehat{\mathrm{y}}}=C_{\ell }^{\mathrm{g}y}+{\alpha }_{\mathrm{CIB}}{C}_{\ell }^{\mathrm{g}{I}_{\mathrm{CIB}}},\end{array}$$(34)

where the galaxy-CIB cross-correlation $C_{\mathcal{l}}^{\text{g}{I}_{\text{CIB}}}$ can be directly measured with the Planck CIB map. For α_CIB, we take the value α_CIB = (2.3 ± 6.6)×10⁻⁷ (MJy sr⁻¹)⁻¹ reported in Alonso et al. (2018).

4.4.2. Cosmic magnification

The measured galaxy overdensity depends not only on the real galaxy distribution but also on lensing magnification induced by the line-of-sight mass distribution (Schneider 1989; Narayan 1989). This magnification, or so-called cosmic magnification, has two effects on the measured galaxy overdensity: (i) overdensities along the line-of-sight cause the local angular separation between source galaxies to increase, so the galaxy spatial distributions are diluted and cross-correlation is suppressed; (ii) lenses along line-of-sight magnify the flux of source galaxies such that fainter galaxies enter the observed sample, so the overdensity increases. These effects bias galaxy-related cross-correlations, especially for high-redshift galaxies (Hui et al. 2007; Ziour & Hui 2008; Hilbert et al. 2009). To take these two effects into account, we modify the expression for the galaxy overdensity (Hui et al. 2007):

$$\begin{array}{c}\hfill {\widehat{\delta }}_{\mathrm{g}}(z)=⟨b_{\mathrm{g}}⟩{\delta }_{\mathrm{m}}(z)+2(2.5s-1)\kappa (z),\end{array}$$(35)

where the second term on the right-hand side of the equation is the cosmic magnification contribution. Here κ(z) is the line-of-sight integral of the lensing convergence to the galaxy redshift z; s is the slope of the logarithmic cumulative number counts of our galaxy sample at magnitude limit m_lim-pagination-14pt

$$\begin{array}{c}\hfill {s\equiv \frac{\partial {log}_{10}N(<m)}{\partial m}|}_{m=m_{lim}}.\end{array}$$(36)

The KiDS gold galaxy sample has an r-band magnitude limit of 25, but the completeness limit is around 24 (Wright et al. 2019). To properly estimate s, one needs a galaxy sample from a deeper survey that has a completeness magnitude of at least 25, as well as the same redshift distribution in each tomographic bin as the KiDS gold galaxy sample. In addition, the slopes for different tomographic bins should be different, but we do not have a good way to estimate them, so we take a simple method to estimate s. Given that the logarithmic cumulative number counts of the galaxy sample is nearly linear with respect to m (see Fig. 6 of Wright et al. 2019)², we estimate s by extrapolating the slope at the completeness magnitude (which is 24) to the magnitude limit, which is 25 (Hildebrandt 2016). The resulting slope is 0.29. We also estimate the error of s by assigning Poisson error in the count bins and measure an error on the slope. This yields an error of 0.001, which is subdominant. We also try other magnitudes to extrapolate from and find that our best-fitting parameters only change marginally. We note that the 2 × 2.5sκ(z) term accounts for flux magnification and 2 × κ(z) term accounts for angular dilution, so when s = 0.4 both effects cancel out.

4.4.3. Uncertainty of the redshift distribution

Uncertainties in the galaxy redshift distributions could affect galaxy cross-correlations. We estimate the uncertainties of the SOM redshift distributions using the same method as described in Hildebrandt et al. (2021), which gives uncertainties of the mean redshift at a level of ∼0.02 in all 5 tomographic bins. To evaluate the impact of this uncertainty, we shift the fiducial redshift distribution n_g(z) in (4) by δ_z  =  {−0.02, −0.01, 0, 0.01, 0.02}. We fit the model parameters with these shifted redshift distributions to see if this changes our results.

5. Results

We estimate the galaxy-tSZ and the galaxy-CMB lensing cross-correlations of KiDS galaxies as well as corresponding covariance matrices in each of the 5 tomographic bins with the methods described in Sect. 4. Figure 5 shows our measurements of $C_{\ell }^{\mathrm{g}\kappa }$ (left column) and $C_{\ell }^{\mathrm{g}\mathit{y}}$ (right column) bandpowers with red dots. The bandpowers are calculated as the mean C_ℓ in each ℓ bin and the error bars are given by the square root of the diagonals of covariance matrices. Each row corresponds to one of the tomographic bins. The measured cross-correlations are over-plotted with the best-fit fiducial model (with the halo model non-linear matter power spectrum) with green lines. Shaded regions are angular scales corresponding to k_cut > 0.7 Mpc⁻¹, which are not included in the model-fit³. We set this threshold because we do not think our simple models will capture the details of the non-linear regime. We note that the red dots show our fiducial results with both CIB contamination and cosmic magnification corrected. In all the tomographic bins, CIB contributions are at a level of ∼1% even with the most conservative level of α_CIB (with α_CIB = (2.3 + 6.6)×10⁻⁷ (MJy sr⁻¹)⁻¹), and can be neglected. We validate this claim by also fitting the raw, CIB contaminated gy cross-correlation with our model and find no significant difference between the best-fit ⟨b_yP_e⟩ values and the fiducial fitting. In order to evaluate the impact of cosmic magnification, we also fit a model without cosmic magnification correction given by Eq. (35). Comparing with our fiducial results, we find that if the cosmic magnification is neglected, ⟨b_g⟩ in all the redshift bins will be slightly under-estimated at a level of ∼1%. We also show the best-fit parameters from shifted redshift distributions in Fig. 6. Points with error bars in different colours correspond to different shifts of redshift distributions δ_z. Our fiducial results have δ_z = 0. From the plot, we conclude that a redshift bias of δ_z ∼ 0.02 would only have a marginal effect on our results. It should be noted that the constraint in the first redshift bin gets mostly affected. This might be due to the fact that the redshift distribution of this bin is the narrowest, which makes it more sensitive to a redshift error.

Fig. 5. Measurements of $C_{\ell }^{\mathrm{g}\kappa }$ and $C_{\ell }^{\mathrm{g}\mathit{y}}$ (red dots with error bars) over-plotted with best-fit models (green lines). The non-linear model shown here is the halo model non-linear template. We also plot the linear part with dashed lines. Each row shows results in each tomographic bin. Shaded regions are scales corresponding to ℓ scales within the threshold kcut > 0.7 Mpc−1, which are not included the model-fit.

Fig. 6. Best fit ⟨byPe⟩ with varying redshift distribution shifts. The non-linear model is the halo model non-linear template. The best-fit parameter values and errors are calculated as the modes and standard deviations of the Gaussian kernel density estimation (KDE) fits of marginalised posterior distribution.

An example of the fiducial MCMC posterior (that corresponds to the halo model non-linear power spectrum) distribution of {⟨b_yP_e⟩, ⟨b_g⟩, ⟨c_gy⟩, ⟨c_gκ⟩} is shown in Fig. 7. We note that the linear biases that we are fitting are the normalisation of the two linear cross-correlations, namely {⟨b_g⟩ × ⟨b_yP_e⟩, ⟨b_g⟩}. Their posteriors are Gaussian in the linear region because they both linearly depend on the linear cross-correlations. The gas pressure bias ⟨b_yP_e⟩ is then the ratio of two Gaussian parameters, so its posterior distribution is asymmetric, as shown in Fig. 7. A summary of the fiducial fitting results of ⟨b_g⟩ and ⟨b_yP_e⟩ is given in Table 2. The best-fit parameter values and errors are calculated as the modes and standard deviations of the Gaussian kernel density estimation (KDE) fittings of the marginalised posterior distributions. We evaluate the constraining power of both linear bias parameters with the method given by Asgari et al. (2021). That is, we calculate the values of the marginalised posterior at both extremes of the prior distribution, and compare them with 0.135, the ratio between the peak of a Gaussian distribution and the height of the 2σ confidence level. If the posterior at the extreme is higher than 0.135, then the parameter boundary is not well constrained. We find that the lower bound of ⟨b_yP_e⟩ of the fifth redshift bin does not meet this criterion (this can also be seen from the fact that the 2σ lower bound of ⟨b_yP_e⟩ in the fifth redshift bin is below zero). However, it should be noted that the lower extremes of the bias parameters are physical limits, which could not be extended. For the parameters that are not bounded, KDE might give inaccurate results. To test that, we calculate the best-fit ⟨b_yP_e⟩ in the fifth redshift bin without smoothing the marginalised posterior with a KDE kernel. This changes the value from 0.16 to 0.134. The difference is well below the constraining error.

Fig. 7. 68% and 95% contours of the posterior distribution of {⟨byPe⟩, ⟨bg⟩, ⟨cgy⟩, ⟨cgκ⟩} in the third redshift bin with the halo model non-linear template in the third redshift bin. The contours and 1D posteriors have been smoothed for the purpose of presentation.

In order to evaluate the goodness-of-fit, we calculate the χ² in each tomographic bin, and then calculate the corresponding probability-to-exceed (PTE) given the degree-of-freedom. Heymans et al. (2021) adopts the criterion PTE > 0.001 (corresponding to a ∼3σ deviation) to be acceptable. We find our fittings of all the tomographic bins with all the three non-linear models meet this criterion. The pure linear model does not fit well, especially in low redshift bins. Thus we conclude that all of our non-linear models fit well with our data.

We show the redshift dependence of ⟨b_g⟩ and ⟨b_yP_e⟩ in Fig. 8. The dots with different colours are the constraints from our results with three non-linear power spectrum models, marginalised over the non-linear parameters ⟨c_gy⟩ and ⟨c_gκ⟩. The plots show that the constraints on ⟨b_g⟩ and ⟨b_yP_e⟩ are consistent with different non-linear power spectrum models, indicating that our measurements are insensitive to the details of the non-linear cross-correlations. To further verify this argument, we repeat the model fitting with different scale cuts k_cut (modes of C_ℓ with a scale smaller than k_cut are removed from the model fitting procedure) and plot the best-fit values of ⟨b_g⟩ and ⟨b_yP_e⟩ as a function of k_cut in Fig. 9. The plot shows that the constraints do not change significantly as different scales are removed, so we conclude that our constraints are robust to non-linear details. To highlight the importance of transition region error, we also fit the halo model and the constant non-linear model without correction with R(k) defined in Eq. (23). We find that this correction changes the best-fit parameter value by a few percent, worst in the lower redshift bin (about 10%), but the differences are below the 1σ level. This is because the Planck beam size ensures that the data are noisy at these scales. However, future studies with higher resolution should be sensitive to these systematics. We also acknowledge that different non-linear models affect ⟨b_yP_e⟩ in the lowest redshift at a level of 0.5 − 1σ because Planck beam does not smooth out the non-linear details as completely as high redshift bins. In Fig. 9, we also plot the fitting of ⟨b_g⟩ and ⟨b_yP_e⟩ with a pure linear model in purple (i.e. ⟨c_gκ⟩ and ⟨c_gy⟩ are both fixed to be zero). We find that the pure linear model gives ⟨b_g⟩ values that are higher than non-linear models on all scales, and have the tendency to merge with the fiducial fitting with low k_cut. The gas pressure bias ⟨b_yP_e⟩ is the ratio between linear amplitudes of gκ and gy cross-correlations, so it could be close to the real value even if the linear model gives biased amplitudes. This result highlights the necessity to include some form of non-linear model in the fitting.

Fig. 8. Constraints of ⟨bg⟩ and ⟨byPe⟩ in each tomographic bin. The best-fitting parameter values and errors are calculated as the modes and standard deviations of the Gaussian KDE fittings of marginalised posterior distributions. Dots with different colours are correspond to the different non-linear power spectrum models. The grey line shows the best-fit model from Chiang et al. (2020).

Fig. 9. Constraints of ⟨bg⟩ and ⟨byPe⟩ for different non-linear matter power spectrum models with different scale cuts kcut. We also plot the fitting of ⟨bg⟩ and ⟨byPe⟩ with a pure linear model in purple (i.e. ⟨cgκ⟩ and ⟨cgy⟩ are both fixed to be zero). We note that for low redshifts and low kcut, we do not have enough degrees-of-freedom so ⟨bg⟩ and ⟨byPe⟩ are not presented. In addition, for high redshift bins, large-scale cuts are beyond the high limit of ℓ, so with these kcut values, the constraints do not change.

In our analysis, we use the whole KiDS gold lensing galaxy sample, in which there are many blue galaxies that are distributed out to a large distance from cluster centres (Croton et al. 2007). The best-fit ⟨b_g⟩ values in all the redshift bins are consistent with one, which suggests that our galaxy sample is a good tracer of the dark matter distribution. This should be contrasted with, for example, luminous red galaxies (LRG), which are strongly biased tracers of mass (Zehavi et al. 2005) because LRGs are known to be clustered around halo centres.

We test the robustness of our model fitting to different Compton y reconstruction method by replacing the Planck MILCA y map with the NILC y map. The constraints on ⟨b_yP_e⟩ is shown in Fig. 10, which indicates a consistency between two y maps. Thus we conclude that our measurement is not sensitive to different y reconstruction methods.

Fig. 10. Constraints of ⟨byPe⟩ with Planck MILCA and NILC y map.

Figure 11 compares our constraints on ⟨b_yP_e⟩ to previous studies. These studies relied on the cross-correlation of data from different surveys with Planck tSZ data. Van Waerbeke et al. (2014; orange dot) uses the lensing data from the RCSLenS sample; Koukoufilippas et al. (2020; purple dots) uses the 2MPZ and WISE × SuperCosmos samples; Pandey et al. (2019; red dots) uses the DES sample; Chiang et al. (2020) uses the galaxy samples from SDSS, BOSS for low redshifts, and QSO samples from SDSS, BOSS and eBOSS for high redshifts. All these studies used the galaxy auto-correlation to measure ⟨b_g⟩, except Van Waerbeke et al. (2014) who used the mass distribution measured from weak gravitational lensing. In our approach, we use the galaxy-CMB lensing cross-correlation to constrain ⟨b_g⟩. It is remarkable that these measurements, from very different surveys and with very different estimators give consistent results for ⟨b_yP_e⟩. However, it should be noted that although our measurement goes to higher redshift, our constraining power is significantly weaker than previous studies in the same redshift range. This is because (1) the limitation of the KiDS footprint makes it less sensitive to linear scales; (2) the CMB lensing map is noisy. With future sky surveys having wider sky coverage and CMB surveys having lower noise levels, these two drawbacks can be improved. The grey line in Fig. 11 is the best-fit redshift dependence of the tSZ halo model given by Chiang et al. (2020). Although this work does not constrain that model, it is introduced and discussed in Appendix B for forecasting future studies. With an agreement with these previous results as well as the halo model prediction, our tomographic measurement provides insights into the thermal history of the LSS.

Fig. 11. Constraints of ⟨byPe⟩ in each tomographic bin. Our results with the halo model non-linear matter power spectrum is presented as black dots with error bars. The best-fit parameter values and errors are calculated as the modes and standard deviations of the Gaussian KDE fittings of marginalised posterior distributions. Results from previous studies are also plotted as well as the best-fit ⟨byPe⟩ model given by Chiang et al. (2020).

The linear bias assumption might break down on small scales where baryonic effects become significant. However, in our analysis, the Planck beam makes our measurements insensitive to these effects. We leave as a future work the generalisation of our results to the small scales when data with higher resolution are available. Such a situation can be handled with a more sophisticated model, for example that from Mead et al. (2020).

In our study, our approach consists of using CMB lensing as a way to constrain ⟨b_g⟩ of a galaxy sample. One can independently measure ⟨b_g⟩ using weak gravitational lensing as a replacement for CMB lensing. Appendix A shows the results when CMB lensing is replaced by the KiDS weak lensing signal. We find that the constraints on ⟨b_yP_e⟩ are consistent with CMB lensing. This result validates our approach and highlights the fact the CMB lensing and galaxy lensing can be used independently to calibrate the mass of a galaxy distribution using a very different source redshift screen. Figure A.2 also shows that the highest redshift bin is significantly noisier when galaxy lensing is used compared to CMB lensing, while galaxy lensing provides higher signal-to-noise for lower redshifts. This is a direct consequence of the very different source redshift between CMB and galaxy lensing, and it illustrates the fact that lensing signal-to-noise decreases dramatically when the lenses are close to the sources, as expected.

6. Discussion and conclusion

In this work we use the galaxy sample from the fourth KiDS Data Release, the Plancky map and Planck CMB lensing map to probe the redshift dependence of galaxy bias of KiDS galaxies and gas pressure bias from the galaxy × tSZ and the galaxy × κ_CMB cross-correlations. We assume that, in the linear region, both tSZ y parameter and galaxy overdensity are proportional to the underlying mass fluctuation, with the proportionality parametrised by galaxy bias ⟨b_g⟩ and gas pressure bias ⟨b_yP_e⟩, which is consistent with our measurement being restricted to large angular scales. To account for the non-linear effects, we also model the non-linear power spectra of gy and gκ cross-correlations as rescaled non-linear templates. We tried three kinds of non-linear templates: halo model, HALOFIT, and constant, all of which yield consistent constraints of ⟨b_g⟩ and ⟨b_yP_e⟩, indicating that our measurements are not yet sensitive to the non-linear details. However, with an additional inconsistent constraint with a purely linear model, we emphasise the necessity to consider non-linear cross-correlations. ⟨b_g⟩ and ⟨b_yP_e⟩ are constrained for galaxies from five tomographic bins within z ≲ 1, which counts amongst the furthest distance probed from this kind of analysis. The reduced χ² of the best-fit parameter values indicate that our model fits the data well.

The best-fitting galaxy bias is close to 1, which indicates that the KiDS galaxy sample is an unbiased tracer of the underlying mass distribution. In previous works (for example Koukoufilippas et al. 2020; Chiang et al. 2020; Pandey et al. 2019) the authors used galaxy auto spectra to constrain galaxy bias, which is subject to modelling uncertainties and auto-correlated noise. Our approach avoids this problem by using the CMB lensing to calibrate the mass from the galaxy distribution. In Appendix A, we will show that our results are unchanged when we replace CMB lensing by galaxy lensing from KiDS.

Figure 11 shows our constraints on ⟨b_yP_e⟩ in each of the tomographic bins as well as the measurements from previous studies (Van Waerbeke et al. 2014; Pandey et al. 2019; Koukoufilippas et al. 2020; Chiang et al. 2020). Our result agrees well with them. We also compare our result with predictions of the halo model (Chiang et al. 2020) and find good agreement. Our tomographic measurement of ⟨b_yP_e⟩ confirms the evolution of biased thermal energy in halos into the high redshift regime. In addition, the gas bias b_y, estimated to be ∼3.5 (Chiang et al. 2020), parametrises the link between gas and dark matter halo; the mean electron pressure $⟨P_{\mathrm{e}}⟩=⟨n_{\mathrm{e}}⟩k_{\mathrm{B}}{\overline{T}}_{\mathrm{e}}$ is associated with the thermal dynamic property of electrons. Based on CMB constraints, the average electron number density is ⟨n_e⟩∼0.25 m⁻³ (Hinshaw et al. 2013). Taking these values into account, the mean electron temperature ${\overline{T}}_{\mathrm{e}}$ is at a level of ${\overline{T}}_{\mathrm{e}}\sim {10}^6$ K, which is consistent with the estimated temperature of ‘missing baryons’ (Cen & Ostriker 1999). This means that if the tSZ signal were entirely from intergalactic gas, it could account for all the missing baryons within the temperature range 10⁵ − 10⁷ K (Bregman 2007). To confirm this, we need a halo model for diffuse baryons that can properly describe the spatial distribution of gas within dark matter halos, which we leave to future work. Our study consolidates our understanding of intergalactic gas at high redshift, which, combined with future tomographic measurements on ⟨b_yP_e⟩, will improve our understanding on the thermal history of the Universe, as well as the evolution of links between gas and dark matter halos.

The uncertainty we find for ⟨b_yP_e⟩ is larger than previous studies because KiDS has a smaller sky coverage compared with those surveys, which makes it less sensitive to the linear regime where the majority of constraining power lies. Uncertainties of ⟨b_g⟩ and ⟨b_yP_e⟩ are both dominated by sample variance on the linear scales. Future sky surveys such as the Rubin Observatory Legacy Survey of Space and Time (LSST; Abell et al. 2009) and the Euclid survey (Laureijs et al. 2010) will cover a larger fraction of the sky, making it possible to yield tighter constrain on linear biases. In addition, future CMB-S4 and Simons Observatory-like experiments will provide CMB lensing and y maps with lower noise levels (Hadzhiyska et al. 2019) and with higher angular resolution, which will improve both the constraining capacity of galaxy bias and sensitivity to small-scale physics. With these improvements, one could model the non-linear cross-correlation with a more sophisticated model, such as the full halo model with a GNFW profile (Arnaud et al. 2010) for the tSZ and the HOD model for the galaxy distribution. We make a forecast on gκ and gy cross-correlations with such a model with hypothetical LSST, Euclid, and CMB-S4 surveys in Appendix B, which yields a tight constraint on ⟨b_yP_e⟩. Our forecast highlights the validity of multi-tracer analysis for future sky surveys.

We carefully evaluate the systematics in our data that could cause bias in our model fitting. The main systematics considered in this study are the cosmic magnification in galaxy overdensity measurements, CIB contamination in tSZ map and uncertainties in the redshift distributions. Though all of these systematics are not significant in our measurement due to low signal-to-noise in our data, they will become significant for future surveys with large sky coverage. In principle, the first two systematics affect high redshifts more significantly, so future studies with a deeper redshift reach must carefully take these into account. It should also be noted that the parameter constraints in this study are quoted under the assumption of fixed cosmological parameters from Planck (Planck Collaboration VIII 2020). The amplitude of angular cross-correlations is closely related to galaxy and gas biases as well as σ₈ and possible reconstruction bias in the CMB lensing map. So these parameters are strongly degenerate, and this is the fundamental limit of this analysis. We could robustly constrain cosmological parameters as well as galaxy and gas biases by combining more cross-correlation measurements, like galaxy clustering and cosmic shear. Once again, we leave this to future studies.

This work shows the potential to study LSS by combining different cross-correlations measurements. cross-correlation is known to be immune to auto-correlated noise. A combination of different cross-correlations can break the degeneracies between model parameters. Specifically, in our fiducial measurements, we do not use cosmic shear, which is affected by intrinsic alignments and shape miscalibration. Instead, we use CMB lensing as a non-biased tracer of LSS to independently constrain ⟨b_g⟩. We provide a sanity check in Appendix A by replacing the CMB lensing map with the KiDS shear map and perform the same analysis, which gives consistent results for all tomographic redshift bins; this validates our fiducial method. However, the results from galaxy lensing at high redshift are noisier than our fiducial results, which indicates the advantage of using CMB lensing as a proxy for the mass distribution, especially at high redshift. Future work could combine CMB lensing and galaxy lensing as independent mass tracers, which could yield tighter constraints on LSS properties. Future surveys will also provide denser galaxy samples in wider ranges and deeper reaches in the sky, as well as cleaner CMB lensing and y maps, which will make this method more promising for multi-tracer cosmology.

Acknowledgments

We thank Prof. Peter Schneider, Prof. Cheng Li, Prof. Houjun Mo, and the anonymous referee for fruitful discussions. ZY and LVW acknowledge support by the University of British Columbia, Canada’s NSERC, and CIFAR. We acknowledge support from the European Research Council under grant numbers 647112 (TT, MA, CH, AJM) and 770935 (AHW, HHi). DA acknowledges support from the Beecroft Trust, and from the Science and Technology Facilities Council through an Ernest Rutherford Fellowship, grant reference ST/P004474. MB is supported by the Polish National Science Center through grants no. 2020/38/E/ST9/00395, 2018/30/E/ST9/00698 and 2018/31/G/ST9/03388, and by the Polish Ministry of Science and Higher Education through grant DIR/WK/2018/12. CH acknowledges support from the Max Planck Society and the Alexander von Humboldt Foundation in the framework of the Max Planck-Humboldt Research Award endowed by the Federal Ministry of Education and Research. HHi is supported by a Heisenberg grant of the Deutsche Forschungsgemeinschaft (Hi 1495/5-1). NK is funded by the Science and Technology Facilities Council (STFC). KK acknowledges support from the Royal Society and Imperial College HYS acknowledges the support from NSFC of China under grant 11973070, the Shanghai Committee of Science and Technology grant No.19ZR1466600 and Key Research Program of Frontier Sciences, CAS, Grant No. ZDBS-LY-7013. Author contributions: All authors contributed to the development and writing of this paper. The authorship list is given in three groups: the lead authors (ZY & LW) followed by two alphabetical groups. The first alphabetical group includes those who are key contributors to both the scientific analysis and the data products. The second group covers those who have either made a significant contribution to the data products, or to the scientific analysis.

References

Appendix A: Alternative method: Galaxy-galaxy lensing cross-correlations

In this appendix, we measure the cross-correlation between galaxy overdensity and lensing-induced galaxy shear as an alternate to CMB lensing. This measurement serves as a sanity check of our fiducial result. We replace the CMB lensing map with the galaxy shear map generated with all the galaxies in the fifth tomographic bin of the gold KiDS lensing shear catalogue and measure the cross-correlation between galaxy overdensity in each tomographic bin and galaxy lensing convergence κ_gl. The kernel of κ_gl is given by:

$$\begin{array}{c}\hfill W^{{\kappa }_{\mathrm{gl}}}(\chi )=\frac{3{\mathrm{\Omega }}_{\mathrm{m}}{H}_0^2}{2a{c}^2}g(\chi ),\end{array}$$(A.1)

with

$$\begin{array}{c}\hfill g(\chi )={\displaystyle {\int }_{\chi }^{{\chi }_{\mathrm{H}}}}\mathrm{d}{\chi }^{\prime }{n}_{\mathrm{g}}(z({\chi }^{\prime }))\chi \frac{{\chi }^{\prime }-\chi }{{\chi }^{\prime }},\end{array}$$(A.2)

where χ_H is the comoving distance to the horizon. For the non-linear region, we take the rescaled halo model described in Sect. 2.

The KiDS lensing shear catalogue (Giblin et al. 2021) provides ellipticities (e₁, e₂) of each galaxy. From them we construct a map triplet of ellipticities (0, −e₁, e₂) (Harnois-Déraps et al. 2017) as an analogue to CMB temperature-polarisation map triplet (T, Q, U). The minus sign on e₁ is due to the different convention of positive x direction between lensing experiments and CMB polarisation experiments. We measure the cross-correlation between the KiDS galaxy overdensity maps and shear map triplet with POLSPICE. The output ‘Δ_gE’ mode (similar to the ‘TE’ mode of CMB) is the cross-correlation between galaxy overdensity and galaxy-galaxy lensing convergence gκ_gl cross-correlation that we want. The ‘Δ_gB’ mode (similar to the ‘TB’ mode of CMB) should be zero, and we take it as a null test.

The covariance matrices are calculated the same way as our fiducial measurement. The cosmic magnification is also corrected. For scales with k > 0.3 h Mpc⁻¹, non-linear galaxy bias becomes significant (Heymans et al. 2021). Besides, we want to match the physical scale of the cross-correlations between gκ_gl and gκ_CMB. Intrinsic alignment is negligible in galaxy-galaxy lensing measurements (Blazek et al. 2012). Take these factors into account; we only use angular scales 100 < ℓ < 200 for the lowest redshift bin and 100 < ℓ < 300 for the rest. The measurements and best-fit models are plotted in Fig. A.1. The ‘Δ_gB’ mode is consistent with zero in all redshift bins. The gκ_gl cross-correlation measured at high redshift are very noisy because for those bins the shear maps contain many galaxies that are in front of lenses, which only contributes to random noise from their intrinsic alignments as well as their cosmic shear from the foreground.

Fig. A.1. Galaxy-shear cross-correlation. Left column: ‘ΔgE’ mode, which corresponds to gκgl cross-correlation and right column: ‘ΔgB’ mode as a null test. Green lines are the best-fit models. We find that the ‘ΔgB’ mode is consistent with zero.

We present the constraints of our model parameters with the galaxy-galaxy lensing cross-correlation in Fig. A.2. The constraints on both ⟨b_g⟩ and ⟨b_yP_e⟩ are consistent. The errors of CMB lensing and galaxy lensing are comparable in the low redshift bins, which is expected if both mass proxies are cut at a comparable scale where the ellipticity noise is subdominant, and the error is mainly driven by sampling variance. This is a proof that CMB lensing and galaxy lensing are not only consistent (which is impressive given ∼12 billion years of separation between the two source populations), but also have comparable signal-to-noise in the sampling variance dominated regime. The uncertainty of the ⟨b_g⟩ in the last two redshift bins is very high due to low signal-to-noise in these bins because in these bins most of the source galaxies are actually in front of lens galaxies. Thus they only contribute noise. This indicates that CMB lensing out-performs galaxy lensing when cross-correlated with high-redshift galaxies. This sanity check not only proves the reliability of our fiducial results but also shows the advantages of CMB lensing experiments over galaxy lensing with the KiDS data. For future studies, a combination of these measurements will yield even better constraint on models of interest.

Fig. A.2. Constraints of model parameters with the galaxy-galaxy lensing and the galaxy-tSZ cross-correlations. Upper panel: constraints on ⟨bg⟩ and lower panel: ⟨byPe⟩ overplotted with halo model results from Chiang et al. (2020).

Appendix B: Forecasting the constraining power of full halo model for future sky surveys

In this section, we forecast the constraining power on the full halo model (Cooray & Sheth 2002; Seljak 2000) parameters from gκ and gy cross-correlations measured from future sky surveys. The galaxy catalogue is assumed to be taken from a LSST/Euclid-like survey; the CMB lensing and tSZ data are from the CMB-S4 (Abazajian et al. 2016) experiment. LSST and Euclid will cover a much wider area of the sky than KiDS, which makes them sensitive to the linear region as well as improves the signal-to-noise. CMB-S4 will achieve a higher angular resolution than Planck, which makes it possible to reveal the details of non-linear cross-correlations.

With such improvements, it is possible to use data from these future surveys to probe the full halo model for galaxies and the tSZ effect. The general halo model divides the power spectrum into the two-halo term, which accounts for the correlation between different halos, and the one-halo term, which accounts for correlations within the same halo, so that

$$\begin{array}{c}\hfill P_{\mathit{UV}}(k)=P_{UV,1\mathrm{h}}(k)+P_{UV,2\mathrm{h}}(k).\end{array}$$(B.1)

Both terms are related to the profiles of U and V in Fourier space:

$$\begin{array}{c}\hfill \begin{array}{cc}& P_{UV,1\mathrm{h}}(k)={\displaystyle {\int }_0^{\infty }}\mathrm{d}M\frac{\mathrm{d}n}{\mathrm{d}M}⟨p_U(k\mid M)p_V(k\mid M)⟩\hfill \\ \hfill & P_{UV,2\mathrm{h}}(k)=⟨b_U⟩(k)⟨b_V⟩(k)P^{\mathrm{lin}}(k)\hfill \\ \hfill & ⟨b_U⟩(k)\equiv {\displaystyle {\int }_0^{\infty }}\mathrm{d}M\frac{\mathrm{d}n}{\mathrm{d}M}{b}_{\mathrm{h}}(M)⟨p_U(k\mid M)⟩,\hfill \end{array}\end{array}$$(B.2)

where P^lin(k) is the linear power spectrum; dn/dM is the halo mass function; b_h is the halo bias and p_U(k ∣ M) is the profile of the tracer U with mass M in Fourier space:

$$\begin{array}{c}\hfill p_U(k\mid M)\equiv 4\pi {\displaystyle {\int }_0^{\infty }}\mathrm{d}r{r}^2\frac{sin(kr)}{\mathit{kr}}{p}_U(r\mid M).\end{array}$$(B.3)

For the one-halo term, one needs to calculate the correlation between different profiles. We take the one-parameter model from Koukoufilippas et al. (2020) to account for the cross-correlation between abundances of u and v:

$$\begin{array}{c}\hfill ⟨p_U(k\mid M)p_V(k\mid M)⟩=(1+{\rho }_{\mathit{UV}})⟨p_U(k\mid M)⟩⟨p_V(k\mid M)⟩.\end{array}$$(B.4)

The details of profiles have been introduced in Sect. 2. For the HOD profile we still fix σ_M = 0.15 and α_s = 1, and let {M₀, M₁, M_min} vary.

For the tSZ profile, we must note that the mass term in this formula is calibrated with X-ray observations and are possibly biased (Planck Collaboration XX 2014). A ‘hydrostatic bias’ b_H is introduced, so that the mass term in (16) is replaced by (1 − b_H)M. In the modelling, we also multiply the power spectra with correction factor R(k, z) defined in (23). Now we can see that the non-linear bias ⟨c_gκ⟩ and ⟨c_gy⟩ are directly related to ρ_GM, ρ_GP, and b_H.

We construct the Fisher matrix:

$$\begin{array}{c}\hfill F_{\alpha \beta }=\frac{\partial {\mathit{M}}^{\mathrm{T}}}{\partial q_{\alpha }}{\mathrm{Cov}}^{-1}\frac{\partial \mathit{M}}{\partial q_{\beta }},\end{array}$$(B.5)

where M(q) is the cross-correlations given by full halo model described above. We also assume that the covariance matrix does not depend on parameters. The free parameters are q ∈ {log₁₀M₁, log₁₀M₀, log₁₀M_min, b_H, ρ_GP, ρ_GM}. The Fisher matrix is calculated at the best-fit parameter values, which are assumed to be {13, 11.68, 11.86, 0.16, −0.5, 0}.

The hypothetical measurement should be extended to non-linear region, so we set the ℓ region to be 50 < ℓ < 3000. The Gaussian covariance matrix is calculated the same way as Eq. (29), with $C_{\ell }^{u\mathit{v}}$ the real angular power spectra signal calculated from best-fit model plus noise spectrum:

$$\begin{array}{c}\hfill C_{\ell }^{\mathit{uv}}=C_{\ell }^{uv,\mathrm{signal}}+N_{\ell }^{\mathit{uv}}.\end{array}$$(B.6)

For cross-correlations, $N_{\ell }^{u\mathit{v}}=0$. The galaxy count noise is taken as the shot noise:

$$\begin{array}{c}\hfill N_{\ell }^{\mathrm{gg}}=1/\overline{N},\end{array}$$(B.7)

where $\overline{N}$ is the mean number of galaxies per steradian on the survey. In this section, we present the forecast for the tomographic bin of a LSST/Euclid-like survey corresponding to the last tomographic bin of KiDS (0.9 < Z_B ≤ 1.2). The redshift distribution and galaxy numbers are from LSST Science Collaboration (2009). The information of LSST and Euclid surveys are summarised in Table B.1. Since both has similar survey coverage and galaxy number density, we only show the forecast for LSST for clarity. We take the noise power spectra of CMB lensing and the tSZ effect for CMB-S4 presented in Shirasaki et al. (2019) and Schaan et al. (2017), respectively.

Since the hypothetical measurement extends into non-linear region, we include the non-Gaussian covariance matrix, which is given by:

$$\begin{array}{cc}\hfill {\mathrm{Cov}}^{\mathrm{NG}}(C_{\ell }^{\mathit{uv}},C_{{\ell }^{\prime }}^{\mathit{wz}})=& {\displaystyle {\int }_0^{\infty }}\mathrm{d}\chi \frac{W^u(\chi )W^v(\chi )W^w(\chi )W^z(\chi )}{4\pi f_{\mathrm{sky}}{\chi }^6}\hfill \\ \hfill & \times T_{\mathit{UVWZ}}(k=\frac{\ell +1/2}{\chi },k\prime =\frac{\ell \prime +1/2}{\chi }),\hfill \end{array}$$(B.8)

where T_UVWZ(k) is the trispectrum. Using the halo model, the trispectrum is decomposed into one- to four- halo terms. Here we only take the one-halo term into account since it dominates the scales we are interested in (Pielorz et al. 2010):

$$\begin{array}{cc}\hfill T_{\mathit{UVWZ}}^{1\mathrm{h}}(k,k\prime )\equiv & {\displaystyle {\int }_0^{\infty }}\mathrm{d}M\frac{\mathrm{d}n}{\mathrm{d}M}\hfill \\ \hfill & \times ⟨p_U(k\mid M)p_V(k\mid M)p_W(k\prime \mid M)p_Z(k\prime \mid M)⟩.\hfill \end{array}$$(B.9)

With these pieces of information, we calculate the covariance matrix as Cov = Cov^G + Cov^NG. The confidence contours from the Fisher matrix is calculated and presented in Fig. B.1.

Fig. B.1. The 95% and 68% confidence contours of model parameters of the full halo model. The confidence contours are calculated from the Fisher matrix for LSST/Euclid-like galaxy surveys and CMB-S4.

The parameter that we are interested in is b_H, which is constrained as 0.16 ± 0.04 for such LSST/Euclid-like surveys. The gas pressure bias ⟨b_yP_e⟩ can then be calculated according to (20), yielding a constraint (0.37 ± 0.02) [meV cm⁻³]. Although the constraining power for galaxy parameters {M₁, M₀, M_min} is weak, the constraining power for gas pressure bias is relatively strong, indicating that galaxy-CMB lensing cross-correlation is a valid method to break the degeneracy between parameters of galaxy distribution and other large-scale tracers when measuring galaxy cross-correlations. Our forecast indicates that future sky surveys like LSST and Euclid as well as CMB-S4 will obtain data that can be used to conduct such measurements.

All Tables

All Figures

	Fig. 1. KiDS-1000 footprint and masks for the Plancky map and CMB lensing map. Regions in purple are masked out.
In the text

	Fig. 2. Redshift distributions of the five tomographic bins of the KiDS gold galaxy sample.
In the text

Fig. 3. Standard deviation calculated from the diagonal of ${\mathrm{Cov}}^{\mathrm{JK}}(C_{\ell }^{\mathrm{g}\mathit{y}},C_{\ell }^{\mathrm{g}\mathit{y}})$ and ${\mathrm{Cov}}^{\mathrm{ANA}}(C_{\ell }^{\mathrm{g}\mathit{y}},C_{\ell }^{\mathrm{g}\mathit{y}})$ in the third tomographic bin. We also plot the standard deviation from the jackknife covariance matrix with Nside = 16 as a consistency check.

In the text

Fig. 4. Full correlation coefficient matrices of covariance matrices (defined in Eq. (30)) of the gy and the gκ cross-correlations in each tomographic bin. Each covariance matrix consists of four sub-matrices, corresponding to the covariances of $C_{\ell }^{\mathrm{g}\mathit{y}}$ and $C_{\ell }^{\mathrm{g}\kappa }$ (block diagonals), as well as their cross-covariance. For each sub-matrix, the pixels show the cross coefficient between binned Cℓ’s, where ℓ bins are defined in Sect. 4.1.

In the text

Fig. 5. Measurements of $C_{\ell }^{\mathrm{g}\kappa }$ and $C_{\ell }^{\mathrm{g}\mathit{y}}$ (red dots with error bars) over-plotted with best-fit models (green lines). The non-linear model shown here is the halo model non-linear template. We also plot the linear part with dashed lines. Each row shows results in each tomographic bin. Shaded regions are scales corresponding to ℓ scales within the threshold kcut > 0.7 Mpc−1, which are not included the model-fit.

In the text

	Fig. 6. Best fit ⟨byPe⟩ with varying redshift distribution shifts. The non-linear model is the halo model non-linear template. The best-fit parameter values and errors are calculated as the modes and standard deviations of the Gaussian kernel density estimation (KDE) fits of marginalised posterior distribution.
In the text

	Fig. 7. 68% and 95% contours of the posterior distribution of {⟨byPe⟩, ⟨bg⟩, ⟨cgy⟩, ⟨cgκ⟩} in the third redshift bin with the halo model non-linear template in the third redshift bin. The contours and 1D posteriors have been smoothed for the purpose of presentation.
In the text

	Fig. 8. Constraints of ⟨bg⟩ and ⟨byPe⟩ in each tomographic bin. The best-fitting parameter values and errors are calculated as the modes and standard deviations of the Gaussian KDE fittings of marginalised posterior distributions. Dots with different colours are correspond to the different non-linear power spectrum models. The grey line shows the best-fit model from Chiang et al. (2020).
In the text

Fig. 9. Constraints of ⟨bg⟩ and ⟨byPe⟩ for different non-linear matter power spectrum models with different scale cuts kcut. We also plot the fitting of ⟨bg⟩ and ⟨byPe⟩ with a pure linear model in purple (i.e. ⟨cgκ⟩ and ⟨cgy⟩ are both fixed to be zero). We note that for low redshifts and low kcut, we do not have enough degrees-of-freedom so ⟨bg⟩ and ⟨byPe⟩ are not presented. In addition, for high redshift bins, large-scale cuts are beyond the high limit of ℓ, so with these kcut values, the constraints do not change.

In the text

	Fig. 10. Constraints of ⟨byPe⟩ with Planck MILCA and NILC y map.
In the text

Fig. 11. Constraints of ⟨byPe⟩ in each tomographic bin. Our results with the halo model non-linear matter power spectrum is presented as black dots with error bars. The best-fit parameter values and errors are calculated as the modes and standard deviations of the Gaussian KDE fittings of marginalised posterior distributions. Results from previous studies are also plotted as well as the best-fit ⟨byPe⟩ model given by Chiang et al. (2020).

In the text

	Fig. A.1. Galaxy-shear cross-correlation. Left column: ‘ΔgE’ mode, which corresponds to gκgl cross-correlation and right column: ‘ΔgB’ mode as a null test. Green lines are the best-fit models. We find that the ‘ΔgB’ mode is consistent with zero.
In the text

	Fig. A.2. Constraints of model parameters with the galaxy-galaxy lensing and the galaxy-tSZ cross-correlations. Upper panel: constraints on ⟨bg⟩ and lower panel: ⟨byPe⟩ overplotted with halo model results from Chiang et al. (2020).
In the text

	Fig. B.1. The 95% and 68% confidence contours of model parameters of the full halo model. The confidence contours are calculated from the Fisher matrix for LSST/Euclid-like galaxy surveys and CMB-S4.
In the text