© ESO 2022

1. Introduction

Gravitational lensing by large-scale structure induces an angular correlation of galaxy shapes on the sky, which in principle contains a wealth of information about the constituents of the Universe. The signal from this ‘weak’ gravitational lensing is determined by both the Universe’s expansion history and the growth of perturbations. Exquisite measurements of hundreds of millions of galaxy shapes that are expected to be made in the coming decade by galaxy surveys, such as the Legacy Survey of Space and Time (LSST) of the Vera C. Rubin Observatory (Ivezić et al. 2019), the Euclid Mission (Laureijs et al. 2011), and the Nancy Grace Roman Space Telescope (Spergel et al. 2015), will allow for the impact of dark energy to be assessed with a precision that has never before been possible. Perhaps most excitingly, the manner of growth of large-scale structure is precisely determined by the underlying gravity theory, and thus general relativity will be robustly tested in a cosmological setting.

The details of how galaxy shapes are distorted is governed by the degree of clustering of large-scale structure between the observed galaxy and Earth. Unfortunately, those structures that contribute most to the correlation of galaxy shapes (over the range of scales where they are best measured) are dense ‘non-linear’ structures whose evolution cannot be calculated via linear (or even quasi-linear) perturbation theory. Structure formation must therefore be simulated, and the clustering statistics extracted from the simulations, in order to compare models to observed data. The simplest clustering statistic is the two-point function, which can be precisely measured from simulations (e.g. Heitmann et al. 2010; although see Smith et al. 2014). Emulators (e.g. Lawrence et al. 2017; Angulo et al. 2021; Euclid Collaboration 2021) can then be constructed to interpolate intelligently between simulation outputs for different cosmological parameters, or physically motivated fitting formulae (e.g. Takahashi et al. 2012; Mead et al. 2021) can be developed.

If gravity were the only important force in the Universe, then the modelling of structure formation at the level of precision required for the measurements from forthcoming surveys would be a solved problem. Unfortunately, non-gravitational processes can have a significant impact on the distribution of matter. Gas cools, contracts, and forms stars, some of which eventually explode as energetic supernovae; black holes form, and complicated feedback mechanisms transfer gravitational energy from accretion events into the gas surrounding galaxies. In cosmology these processes are known collectively as ‘baryonic feedback’, and the complex physics of how they operate is poorly understood in comparison with pure gravitational clustering. While it is possible to run accurate hydrodynamic simulations that incorporate gas physics (e.g. Springel 2010) such as cooling and radiating, the details of processes such as active galactic nucleus (AGN) and supernova feedback must be included schematically, via sub-grid recipes (e.g. Schaye et al. 2010; McCarthy et al. 2017), and it is difficult to assess the numerical convergence and accuracy of such simulations. The unknown impact of feedback on the matter distribution is one of the biggest uncertainties in current lensing analyses (e.g. Chisari et al. 2019b) and if not addressed will limit the ability of forthcoming surveys to test dark energy and gravity theories. Attempts have been made to use these effects of baryon feedback on lensing observables as a probe of galaxy formation (e.g. Harnois-Déraps et al. 2015; Foreman et al. 2016; MacCrann et al. 2017; Huang et al. 2021; Yoon & Jee 2021), but the inferred constraints on baryon feedback have been weak for past and current surveys.

The thermal Sunyaev–Zeldovich (tSZ; Sunyaev & Zeldovich 1972) effect arises when cosmic microwave background (CMB) photons are Compton-scattered off free electrons, predominantly hot electrons in galaxy clusters. This results in a spectral distortion of the CMB black-body spectrum with a characteristic frequency dependence. Maps of this Compton-y signal can be extracted from multi-frequency CMB temperature data (e.g. Planck Collaboration XXII 2016). The tSZ effect therefore traces hot ionised gas in the Universe, which is itself correlated with large-scale structure. This correlation contains useful cosmological information (e.g. Planck Collaboration XXII 2016; Horowitz & Seljak 2017; Bolliet et al. 2018; Makiya et al. 2020; Mitchell et al. 2021) and exhibits a particularly strong scaling with the amplitude of perturbations (Refregier & Teyssier 2002; Komatsu & Seljak 2002). Of particular interest to weak-lensing cosmology is the fact that it is this same hot gas that has been redistributed by baryonic feedback, and therefore the tSZ effect provides an important external constraint on the magnitude of baryonic feedback.

The tSZ–galaxy cross-correlation has been investigated (e.g. Vikram et al. 2017; Makiya et al. 2018; Pandey et al. 2019; Tanimura et al. 2020; Chiang et al. 2020; Koukoufilippas et al. 2020; Schaan et al. 2021; Amodeo et al. 2021; Yan et al. 2021), but constraints on cosmological parameters are dependent on the unknown galaxy bias, which must either be taken from external data or else modelled in the data analysis. Stacking of tSZ observations of large-scale structure filaments (e.g. Tanimura et al. 2019; de Graaff et al. 2019) has shed light on the distribution of diffuse gas. In this paper we utilise measurements of weak lensing and the tSZ effect to place simultaneous constraints on the cosmological model and on the strength of baryonic feedback. The cross-correlation between lensing and the tSZ effect has been measured previously for both CMB lensing (e.g. Hill & Spergel 2014) and for galaxy lensing (e.g. van Waerbeke et al. 2014; Hojjati et al. 2017; Osato et al. 2020; Gatti et al. 2021; Pandey et al. 2021). In such works, model parameters were either fixed or only varied within tight priors, with no coupling between parameters that describe the effect of baryon feedback on the lensing and the tSZ effect. Here we build on these previous works and analyse the tomographic cross-correlation between the tSZ effect and weak lensing from the 1000 deg² release of the Kilo-Degree Survey (KiDS-1000) in a fully Bayesian manner, deriving posteriors on cosmological and baryon feedback parameters while marginalising over redshift uncertainties, the intrinsic alignment (IA) of galaxies, and contamination by the cosmic infrared background (CIB). We then combine the lensing–tSZ cross-correlation with KiDS-1000 cosmic shear in a novel joint analysis. By using the Mead et al. (2020) HMX model, which consistently accounts for gas, stars, and dark matter in the whole joint-data vector, we are able to infer constraints on cosmology, as well as the ‘AGN temperature’ parameter that governs the impact of feedback on both the lensing–tSZ cross-correlation and cosmic shear.

This paper is organised as follows: In Sect. 2 we summarise the lensing and tSZ data, as well as the estimators and summary statistics used; in Sect. 3 we present the modelling of the signal, its covariance, and systematics; Sect. 4 presents the results, including measured data vectors and inferred cosmological constraints; and we summarise our findings in Sect. 5.

2. Data

2.1. Kilo-Degree Survey

We analysed the fourth data release of the Kilo-Degree Survey (Kuijken et al. 2019), KiDS-1000, which spans ∼1000 deg², split equally between an equatorial and southern imaging stripe. KiDS is unique in terms of its complete matched-depth nine-band imaging across the optical to the near-infrared, with relatively uniform and high-resolution seeing in the r-band. This combination of spatial and wavelength resolution provides an exquisite dataset for weak gravitational lensing studies, allowing for sensitive control of data-related systematics (Giblin et al. 2021; Hildebrandt et al. 2021).

Galaxy shear estimates are determined using the model-fitting lensfit method (Miller et al. 2013), validated and calibrated using realistic KiDS image simulations (Kannawadi et al. 2019). The percent-level uncertainty on this calibration correction is accounted for in our covariance calculation (see Sect. 3.4 of Joachimi et al. 2021, for details). Tomographic galaxy redshift distributions are determined using overlapping spectroscopic training data and a self-organising map to account for incompleteness in the training sample (Wright et al. 2020a). The accuracy of the resulting distributions are validated with a cross-correlation clustering analysis and mock survey analysis (Hildebrandt et al. 2021). The measured uncertainty from these tests is used as a prior for a set of correlated nuisance parameters that account for any residual correlated calibration errors in the mean redshift of each tomographic bin (see Sect. 3.3 of Joachimi et al. 2021, and Sect. 3.4.1 for details). The full suite of KiDS-1000 data products, including imaging, catalogues, shear and redshift calibrations, is available to download¹.

2.2. Compton-y maps

The construction of Compton-y maps, which trace the tSZ effect, usually proceeds by taking linear combinations of intensity maps at different frequencies of CMB experiments, potentially with scale-dependent or spatially varying weights. Planck Collaboration XXII (2016) provide two such Compton-y map reconstructions, using two different internal linear combination (ILC) methods: MILCA (modified internal linear combination algorithm; Hurier et al. 2013) and NILC (needlet independent linear combination; Remazeilles et al. 2011). We used the MILCA map as our fiducial y-map in our analysis but confirmed that the results hold for the NILC map as well. Furthermore, we checked our sensitivity to CIB contamination by using Compton-y maps from Yan et al. (2019) where the CIB has been subtracted from the Planck frequency maps. Assessment and mitigation of possible CIB contamination of our y-maps is further discussed in Sect. 3.4.3.

We applied a Galactic mask that covers the 40% of the sky with the strongest dust emissions, as well as a point-source mask, as provided by Planck Collaboration XXII (2016). The extra-galactic lensing signal is assumed not to be correlated with emissions from the Milky Way, such that the Galactic mask primarily serves to reduce noise in the measurement. Furthermore, as the KiDS-1000 footprint avoids the Galactic plane, the Galactic mask has only a small effect on the overlap between the shear and Compton-y maps, such that the unmasked area of the overlap between the KiDS-1000 and Planck data makes up 85% of the footprint of the KiDS-1000 data. The maps are provided in HEALP IX (Górski et al. 2005) format, with nside = 2048, in Galactic coordinates.

The Planck-derived y-maps have a resolution of 10 arcmin, which limits the amount of small-scale information that can be extracted. In contrast, Compton-y maps constructed from the combination of the fourth data release of the Atacama Cosmology Telescope (ACT; Mallaby-Kay et al. 2021) and Planck data, benefit from a much smaller beam of 1.6 arcmin (Madhavacheril et al. 2020). Unlike the Planck-derived maps, the ACT y-maps only partially cover the northern KiDS-1000 patch, while they do not overlap with the southern KiDS-1000 patch at all. After applying the same 40% Galactic and point source masks as for the Planck-derived maps, the overlap between the ACT and KiDS-1000 data is reduced to 30% of the area covered by the KiDS-1000 data, or 35% of the area of the overlap between the Planck-derived y-maps and KiDS-1000.

We transformed the ACT y-maps into HEALP IX format (nside = 2048) using PIXELL². Unlike the Planck-derived maps, the ACT maps are provided in equatorial coordinates, the same as the KiDS-1000 data. When measuring the shear–tSZ cross-correlations with the Planck-derived maps, we used a shear catalogue that had been transformed into Galactic coordinates (see e.g. Xia et al. 2020), which avoids interpolation artefacts that would be introduced if the y-map were transformed into a different coordinate system.

2.3. Summary statistics and estimators

Our summary statistic are the angular (cross-)power spectra of the shear field and tSZ effect. We estimated the angular power spectra using the pseudo-C_ℓ method, specifically the MASTER (Hivon et al. 2002) estimator as implemented in the NAMASTER³ (Alonso et al. 2019) software. The method proceeds by first naively estimating the angular power spectra of the masked maps. These estimates are biased and thus give the method its name: pseudo-C_ℓ. The bias introduced by the masks can be calculated and is encoded in the so-called mode-mixing matrices. Given the mode-mixing matrices, an unbiased estimate of the (binned) angular power spectra can be obtained from the biased pseudo-C_ℓ estimates. We refer the interested reader to Alonso et al. (2019), who map these steps out in detail (see their Eqs. (7), (9), (D.2), and (13) for the steps outlined above).

Estimating shear power spectra using pseudo-C_ℓ methods comes with a number of subtleties due to the amount of small-scale power in the masks involved, which can affect the numerical stability of the algorithm. Here we closely followed the methodology of Nicola et al. (2021). We pixelised the shear catalogues into HEALP IX maps $\widehat{\mathit{\gamma }}$ with nside = 2048:

$$\begin{array}{c}\hfill \widehat{\mathit{\gamma }}({\mathit{n}}_i)=\frac{1}{w({\mathit{n}}_i)}{\displaystyle \sum _{j\in {\mathit{n}}_i}{w}_j{\mathit{ϵ}}_j,}\end{array}$$(1)

where n_i denotes the position of the ith pixel and the sum runs over all galaxies in that pixel, which are characterised by their weight w_j and the two ellipticity components, which we denote here with ϵ_j. In our case, the weights w_j are given by the lensfit weights and account for uncertainty in the shape measurement. The weight map w(n_i) is defined as

$$\begin{array}{c}\hfill w({\mathit{n}}_i)={\displaystyle \sum _{j\in {\mathit{n}}_i}{w}_j.}\end{array}$$(2)

The multiplicative shear bias was accounted for by scaling all ellipticity estimates in a tomographic bin by $\mathit{ϵ}\to \frac{1}{1+m_z}\mathit{ϵ}$, where we used the multiplicative bias values m_z in Giblin et al. (2021). The additive shear bias was accounted for by subtracting the weighted ellipticity mean in each tomographic bin, following Giblin et al. (2021). The weight maps of Eq. (2) take the role of the masks for the shear maps of the pseudo-C_ℓ estimator. All cross-power spectra were estimated using their respective masks, that is, the measurements are not restricted to the overlap region.

The auto-power spectra are subject to a noise bias, which we accounted for analytically, using the estimate (Nicola et al. 2021)

$$\begin{array}{c}\hfill N_{\ell }={\mathrm{\Omega }}_{\mathrm{pix}}{\displaystyle {\langle \sum _{j\in {\mathit{n}}_i}{w}_j^2{\sigma }_{\gamma ,j}^2\rangle }_{\mathrm{pix}},}\end{array}$$(3)

where Ω_pix is the pixel area, ${\sigma }_{\gamma ,j}^2=\frac{1}{2}{|{\mathit{ϵ}}_j|}^2$ with $|\mathit{ϵ}{|}^2={ϵ}_1^2+{ϵ}_2^2$, and the average is taken over all pixels. The accuracy of this estimate was confirmed using a large number of Gaussian random field mocks with realistic noise and the real data weight maps of Eq. (2).

The angular power spectra were estimated in 12 logarithmically spaced bins between ℓ = 51 and ℓ = 2952. When discarding the lowest two and highest two ℓ-bins, this binning scheme results in the same eight logarithmically spaced multipole bins between ℓ = 100 and ℓ = 1500 that were used in the KiDS-1000 analyses (Joachimi et al. 2021; Asgari et al. 2021; Heymans et al. 2021; Tröster et al. 2021). Since the KiDS-1000 cosmic shear methodology was validated only on the ℓ ∈ (100, 1500) range, we restricted our analysis to the same ℓ-range. Furthermore, this ℓ-range captures most of the signal, with the signal-to-noise ratio increasing by only ∼10% for the full ℓ ∈ (51, 2952) range.

We note that our pseudo-C_ℓ approach differs from the KiDS-1000 cosmic shear band-power methodology of Asgari et al. (2021) and from the recent pseudo-C_ℓ methodology of Loureiro et al. (2021). Our analysis choice, which does not impact the KiDS-1000 cosmic shear-only constraints (see Appendix B), is motivated by the requirements of the lensing–tSZ cross-correlation analysis and support for the computation of the full spin-2 Gaussian data covariance by NAMASTER (see Sect. 3.3). The angular power spectrum estimates were also checked against those estimated with POLSPICE⁴ (Szapudi et al. 2001; Chon et al. 2004) and found to be in excellent agreement. The measurement scripts, associated tools, and the measurements themselves are being made publicly available⁵.

3. Modelling

3.1. Three-dimensional power spectra

To model our data and extract meaningful cosmological constraints we followed Mead et al. (2020) and used the HMX approach, which is designed to reproduce power spectra measured from the BAHAMAS (McCarthy et al. 2017) suite of cosmological hydrodynamic simulations, which in turn are calibrated to reproduce the observed galaxy stellar mass function and baryon fraction in groups and clusters. HMX works at the three-dimensional level, providing P(k) rather than C_ℓ, because once P(k) is accurate the C_ℓ are automatically also accurate via the projection integrals; modelling and physical understanding are also more straightforward in three dimensions where the fields naturally live. The three-dimensional source of the tSZ effect is the electron pressure, while that of shear is the matter overdensity, and HMX reproduces the electron pressure–matter power spectrum from BAHAMAS at the 15 percent level for z < 2 and k ≲ 7 h Mpc⁻¹ while at the same time providing matter–matter power spectra that are accurate at the 2 percent level⁶. Ultimately we wish to use HMX to extract meaningful constraints from combined shear–shear and shear–tSZ data.

Underpinning HMX is the halo model (e.g. Cooray & Sheth 2002), where halo density profiles are broken down into their constituent dark matter, gas and stellar components. The gas is then given thermal properties such that the halo electron pressure can be simultaneously calculated with the matter. This occurs both at the level of individual haloes and the expelled background gas, which is assumed to be heated by AGN activity. Fitting parameters of the model to data from BAHAMAS provides the link in HMX between gas thermal properties and the density. Any standard halo model cross-correlation calculation (e.g. Hill & Spergel 2014; Ma et al. 2015) provides a degree of cross-talk between the two fields through the common linear spectrum, halo mass function and bias, but most often the halo profiles of the two fields are taken to be independent. For example, to model the shear–tSZ signal, an Arnaud et al. (2010) profile may be adopted for the pressure while a (completely separate) Navarro et al. (1997; NFW) profile may be adopted for the matter. The advantage of HMX is that the halo profiles are also directly coupled to each other, which allows the electron pressure (which probes the hot gas) to directly affect the overall matter profile. The model parameters, consisting of the concentration of dark matter halo density profiles and its redshift evolution, the polytropic index of gas bound in haloes, the halo mass below which haloes have lost half their initial gas content, the ratio of halo temperature to that of the virial equilibrium, and the temperature of the warm-hot intergalactic medium with its redshift evolution, are fit to power spectra from simulations with different feedback strengths. All fitted parameters are then related to a single ‘AGN temperature’ parameter, log₁₀(T_AGN/K), which governs the strength of feedback; increasing the temperature increases the tSZ signal while simultaneously decreasing the matter signal as gas is forced from the halo by the increased AGN activity. This is particularly helpful in gaining an understanding of the impact baryonic feedback may have on the matter power spectrum, currently one of the largest uncertainty in cosmic shear studies at small scales (e.g. Chisari et al. 2019b).

The halo model power spectrum calculation is known to be inaccurate in detail, particularly in the ‘transition region’ between the two- and one-halo terms (Tinker et al. 2005; Valageas & Nishimichi 2011; Fedeli et al. 2014), which is primarily because linear bias is an overly simplistic treatment of halo–halo correlations (Mead & Verde 2021). There are also other errors at larger scales (for example, the lack of perturbation theory) and smaller scales (such as assumed halo sphericity) that arise because of other simplifying assumptions that enter into the standard halo-model calculation. To side-step these problems, rather than predicting the raw power, HMX predicts the ratio of the power spectrum in question to the gravity-only matter–matter spectrum (adopting NFW profiles for all matter) calculated for the same cosmological model. The halo-model prediction for this ratio has been shown to be more robust than the absolute power (e.g. Cataneo et al. 2019) as the underlying cosmology is varied⁷. Once we have this halo-model ratio, we must multiply by an accurate prescription for the matter–matter power, such as provided by an emulator, or a semi-analytic model, such as Mead et al. (2021). Here we used the HMCODE-2016 (Mead et al. 2016) model for this matter–matter power prediction, since the more accurate HMCODE-2020 (Mead et al. 2021) was not yet available when this analysis was performed.

3.2. Angular power spectra

Neither the three-dimensional matter–matter nor the matter–electron pressure power spectra are directly observable. Instead we observe the angular auto and cross-power spectra of the gravitational lensing shear effect induced by the matter distribution and the tSZ effect caused by the electron pressure. The angular power spectrum $C_{\mathcal{l}}^{AB}$ at multipole ℓ between fields A and B can be related to the three-dimensional power spectrum P_AB(k, z) corresponding to these fields through (e.g. Kaiser 1992; LoVerde & Afshordi 2008; Kilbinger et al. 2017)

$$\begin{array}{c}\hfill C_{\ell }^{\mathit{AB}}={\displaystyle {\int }_0^{{\chi }_{\mathrm{H}}}\frac{\mathrm{d}\chi }{{\chi }^2}{W}^A(\chi )W^B(\chi )P_{\mathit{AB}}(\frac{\ell +\frac{1}{2}}{\chi },z(\chi )).}\end{array}$$(4)

The integral is taken over the comoving distance χ, from the observer to the distance of the horizon χ_H. Throughout this work we shall assume a flat Λ cold dark matter cosmology with a minimal neutrino mass, as there has been no evidence for extensions to this model (e.g. Planck Collaboration VI 2020; Alam et al. 2021; Tröster et al. 2021). The shear field derives from the three-dimensional matter overdensity field. Hence, for $C_{\mathcal{l}}^{\gamma \gamma }$ the power spectrum P_AB in Eq. (4) is the matter-overdensity power spectrum P_mm. In the case of the angular cross-power spectrum between shear and the tSZ effect, $C_{\mathcal{l}}^{\gamma y}$, the relevant three-dimensional power spectrum is the cross-power spectrum of matter overdensity and electron pressure, P_mP. Unless noted otherwise, the auto- and cross-correlations involving the shear field shall refer to the E-mode of the E/B-mode decomposition of the spin-2 shear field. That is, we shall use $C_{\mathcal{l}}^{\gamma \gamma }$ and $C_{\mathcal{l}}^{\gamma y}$ as a notational shorthand for $C_{\mathcal{l}}^{{\gamma }_E{\gamma }_E}$ and $C_{\mathcal{l}}^{{\gamma }_{E}y}$, respectively. The window function W^γ_i(χ) of the shear field of the ith source sample is given by

$$\begin{array}{c}\hfill W^{{\gamma }_i}(\chi )=\frac{3}{2}{\left(\frac{H_0}{c}\right)}^2{\mathrm{\Omega }}_{\mathrm{m}}\frac{\chi }{a(\chi )}{\displaystyle {\int }_{\chi }^{{\chi }_{\mathrm{H}}}\mathrm{d}{\chi }^{\prime }{n}_i({\chi }^{\prime })\frac{{\chi }^{\prime }-\chi }{{\chi }^{\prime }},}\end{array}$$(5)

where H₀ denotes the present-day Hubble rate, c the speed of light, Ω_m the matter density, a(χ) the scale factor at comoving distance χ, and n_i(χ)dχ = n_i(z)dz the redshift distribution of source sample i. The window function W^y(χ) of the tSZ effect is

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

where σ_T is the Thompson scattering cross-section and m_e is the electron mass. Here the three-dimensional power spectra are given in comoving units, including the electron pressure. That is, the matter–electron pressure cross-power spectrum P_mP has units of [(comoving length)³ × energy/(comoving length)³]=[energy]. The redshift distributions n_i(z) of the source samples, as well as the window functions of Eqs. (5) and (6) are illustrated in Fig. 1.

Fig. 1. Top: redshift distributions ni(z) of the five tomographic source samples (solid lines). The shaded bands indicate the selection window on the photometric redshift point estimate zB. Bottom: line-of-sight window function, W(z), to project the three-dimensional power spectra P(k) (in comoving units) into their angular counterparts, Cℓ. The lensing window functions (solid lines, Eq. (5)) are normalised by the maximum of the window function of the last tomographic bin, while the Compton-y window function (dashed line, Eq. (6)) is normalised by its maximum value.

In order to compare our model of the angular power spectra in Eq. (4) to the measured power spectra described in Sect. 2.3, we need to account for the binning in multipoles, residual effects of the mode-mixing matrices, and E/B-mode mixing. To this end we apply the pseudo-C_ℓ band-power window function ${\mathcal{F}}_{b\ell }^{\mathit{AB}}$, which maps the angular power spectrum $C_{\mathcal{l}}^{AB}$ to the pseudo-C_ℓ band-power power spectrum $C_b^{AB}$ at bin b:

$$\begin{array}{c}\hfill C_b^{\mathit{AB}}={\displaystyle \sum _{\ell }{\mathcal{F}}_{b\ell }^{\mathit{AB}}{C}_{\ell }^{\mathit{AB}}.}\end{array}$$(7)

The pseudo-C_ℓ band-power window function ${\mathcal{F}}_{b\ell }^{\mathit{AB}}$ can be computed analytically and is given by Eq. (16) in Alonso et al. (2019). Since our model does not predict any B-modes, E/B-mode mixing does not affect the model for the E-modes. These window functions are illustrated in Fig. B.2 for the case of E-mode cosmic shear.

3.3. Covariance

To estimate the uncertainty associated with our measured data vector, we largely followed the methodology of Joachimi et al. (2021), with minor modifications to account for the pseudo-C_ℓ estimator and allow the computation of the tSZ parts of the covariance matrix. In contrast to simulation or resampling-based methods, this analytic approach is free of sampling noise and allows us to cleanly separate the different covariance contributions, as well as match the cosmology to that preferred by the data.

We decompose the covariance $\mathsf{C}=\text{Cov}[{\widehat{C}}_{{\ell }_1}^{\mathit{ij}},{\widehat{C}}_{{\ell }_2}^{\mathit{kl}}]$ of the measured angular (cross-)power spectra ${\widehat{C}}_{\ell }^{\mathit{ij}}$ between fields i and j at the effective multipole ℓ into

$$\begin{array}{c}\hfill \mathsf{C}={\mathsf{C}}_{\mathrm{G}}+{\mathsf{C}}_{\mathrm{SSC}}+{\mathsf{C}}_{\mathrm{T}}+{\mathsf{C}}_{\mathrm{m}},\end{array}$$(8)

where C_G is the Gaussian contribution, C_SSC is the super-sample covariance (SSC), C_T is the non-Gaussian contribution due to the non-zero trispectrum, and C_m accounts for the uncertainty in the shear calibration.

We computed the Gaussian covariance term C_G using NAMASTER, accounting for the effect of partially overlapping and weighted sky observations (García-García et al. 2019). While this can be done exactly, the calculation scales as $O({\mathcal{l}}_{\text{max}}^6)$ for the sample variance term. Instead, we used the improved narrow kernel approximation (NKA) introduced in Nicola et al. (2021) for the cosmic shear sample variance and the parts of the covariance that involve tSZ cross-correlations. For cosmic shear, the noise terms can be computed exactly with little extra effort by calculating the different terms in Eq. (2.29) of Nicola et al. (2021) separately. We found that this exact calculation reduces the Gaussian cosmic shear covariance by up to 20% compared to the case where the NKA is used for both sample variance and noise terms. Some of the Compton-y maps have no readily available noise maps necessary for the exact noise computation, and we thus used the NKA for all terms of covariance involving the tSZ effect. As the signal-to-noise ratio of the tSZ cross-correlation is lower than that of the cosmic shear measurement, this is of little concern, however. The sample variance of the cosmic shear and shear–tSZ terms in the Gaussian covariance were computed using the model for the angular power spectra of Sect. 3.2 at the best-fit cosmology of Heymans et al. (2021). By contrast, the tSZ auto-power spectrum term is based on a smoothed version of the measured auto-power spectrum of the y-maps. This choice is motivated by the fact that HMX was not calibrated for accurate predictions of the pressure auto-power spectrum. Since the overlap between the KiDS-1000 footprint and the Compton-y maps constitutes a small fraction of the full y-maps, the measured tSZ auto-spectra are largely independent of the shear spectra.

Osato & Takada (2021) find that the SSC term can be the dominant contribution to the data covariance in analyses of the tSZ auto-power spectrum. In cosmic shear analyses, the SSC term is the dominant non-Gaussian contribution (e.g. Barreira et al. 2018). We thus accounted for SSC in the full cosmic shear, shear–tSZ, and joint covariance. The SSC C_SSC and trispectrum terms C_T were computed using the halo model formalism as implemented in CCL⁸ (Chisari et al. 2019a). The SSC terms of the cosmic shear and shear–tSZ covariances used the power and cross-power spectra of the respective survey masks, while the SSC term of the cross-covariance between the two data vectors used the cross-power spectrum of shear and tSZ masks. The only terms in the full covariance Eq. (8) that did not depend on the shape of the survey footprints are the trispectrum terms, which used the effective survey areas and overlap-areas instead. The halo matter density and pressure profiles were computed with the HMX prescription, using PYHMCODE⁹, such that the modelling of the signal and the covariance is consistent. The version of CCL used in this work only includes the one-halo contribution to the trispectrum term. This is not an issue, however, since for cosmic shear the trispectrum term is subdominant everywhere (Joachimi et al. 2021) and for the tSZ effect it is dominated by the one-halo term (Cooray 2001). The SSC term was found to be sub-dominant for the shear–tSZ cross-correlation but can reach up to ∼20% of the total covariance for off-diagonal terms. This is illustrated in Fig. 2, showing the relative contribution of the SSC to the total covariance. We found that not accounting for the SSC term led to an underestimate of the inferred parameter uncertainties of about 10%.

Fig. 2. Joint cosmic shear and shear–tSZ data covariance. Upper right: contribution of the super sample covariance (SSC) to the joint covariance. Shown is the ratio between the absolute value of the SSC contribution and the absolute value of the total covariance in Eq. (8). The ratio can exceed unity, since different contributions to total covariance can cancel one another out to a certain degree. Lower left: correlation coefficients corresponding to the data covariance matrix.

To account for the uncertainty of the multiplicative shear bias C_m, we followed Joachimi et al. (2021), that is, we assumed that the multiplicative bias corrections are fully correlated between the tomographic bins. Asgari et al. (2021) show that this approach gives very similar constraints to allowing the multiplicative biases to vary, while not requiring sampling over an extra five parameters.

Another approach to the estimation of the data covariance is through large ensembles of simulated mock data. This approach has some considerable downsides for the present analysis, however. Obtaining mock data that accurately models both the tSZ effect, as well as weak lensing, over large sky areas and with a sufficient number of realisations to reduce the sampling noise, is a challenging proposal. While Tröster et al. (2019) demonstrate a deep-learning method to augment weak-lensing-optimised N-body simulations with gas physics that accurately predict the tSZ effect, this approach is not easily amenable to different cosmological parameters. Furthermore, accounting for the SSC is difficult, unless one employs costly ‘separate universe’ simulations (Li et al. 2014). Finally, the covariance is dominated by the Gaussian term (Joachimi et al. 2021), which can be accurately computed analytically, even for large survey areas.

Our data covariance model was checked against jackknife estimates based on various patch sizes and jackknife estimators. We found agreement at the 10% level on average but with large scatter between different jackknife estimates. The jackknife approaches tended to slightly overestimates (underestimates) the covariance at small (large) scales. However, since the jackknife patches are heteroskedastic due to inhomogeneous masks and weight maps, as well as breaking the assumption of independence between the jackknife patches at large scales, these jackknife estimates need to be interpreted with caution.

3.4. Systematics

3.4.1. Photometric redshift uncertainties

The redshift distributions of the KiDS-1000 source samples have been estimated and validated to a high degree of accuracy (Hildebrandt et al. 2021). To account for the residual uncertainty in the redshift distributions of the source samples, we followed Asgari et al. (2021) and allowed the redshift distribution for each tomographic bin to shift. We imposed a multivariate Gaussian prior on these shifts, where the mean and covariance were derived from the calibration of the redshift estimation method to realistic mock data. We refer the interested reader to Wright et al. (2020b), Hildebrandt et al. (2021), and Asgari et al. (2021) for further details on the estimation of the redshift distributions and their uncertainties.

3.4.2. Intrinsic alignments

Since the tSZ effect is correlated with large-scale structure, source galaxies aligning with the same large-scale structure can mimic a gravitational lensing signal, and thus the correlation between the tSZ effect and these IA effects could bias our shear–tSZ cross-correlation analysis. We accounted for this IA effect using the same model as in the KiDS-1000 analyses (see Joachimi et al. 2021; Asgari et al. 2021; Heymans et al. 2021). Specifically, we employed the non-linear alignment (NLA; Hirata & Seljak 2004; Bridle & King 2007) model and allowed the amplitude parameter A_IA to vary, while keeping the redshift-evolution parameter η_IA fixed. Asgari et al. (2021) show that letting η_IA vary serves to slightly loosen the constraints on S₈ but does not affect its central value. For a detailed discussion of the choice of IA model for the KiDS-1000 data, we refer the interested reader to Appendix D of Heymans et al. (2021). We found that IA has a small but significant impact on the shear–tSZ cross-correlation signal, in agreement with findings of Gatti et al. (2021) and Pandey et al. (2021). Its relative contribution to the total signal, as well as its fraction of the measurement error, is comparable to the IA contribution in the case of cosmic shear.

3.4.3. Cosmic infrared background

The tSZ maps are contaminated by CIB (Hurier 2015; Planck Collaboration XXII 2016; Planck Collaboration XXIII 2016). To assess and correct for this contamination in our shear–tSZ cross-correlation, we considered two complementary approaches. The first approach models the effect of the CIB contamination based on the observed correlation between the shear field and maps of the CIB, while the second approach uses Compton-y maps where the CIB has been subtracted during the reconstruction process.

For the first, model-based approach, we followed the simple but commonly used method of Hill & Spergel (2014) and Vikram et al. (2017), where the contamination of the observed Compton-y map y_obs is modelled as

$$\begin{array}{c}\hfill y_{\mathrm{obs}}=y_{\mathrm{true}}+{\alpha }_{\mathrm{CIB}}{I}_{\mathrm{CIB}}+N,\end{array}$$(9)

where y_true is the map of the uncontaminated tSZ effect, I_CIB is a map tracing the CIB, α_CIB is a scalar parameter, and N denotes the collection of noise and other contaminants that are not expected to correlate with the shear field, unlike the other two terms.

This model is quite simplistic; the use of a single frequency CIB map assumes that the CIB is not varying with frequency, and the single, scalar parameter α_CIB assumes that the contamination is not varying spatially. Neither of these assumptions strictly hold true: the CIB is the result of complex astrophysical processes, whose detailed frequency dependence is not constant, and the reconstruction algorithms used for the y-maps in this work employ scale-dependent filters, resulting in scale-dependent ILC coefficients, and would thus require a scale-dependent α_CIB as well.

In practice, these assumptions appear to be sufficient, however, at least for the applications considered in the present work, since we found CIB contamination to be negligible (see Sect. 4.2). Compton-y maps using scale-independent ILC weights (e.g. van Waerbeke et al. 2014; Hill & Spergel 2014) are consistent with those using MILCA and NILC reconstructions used here and the CIB is strongly correlated between different frequency bands (e.g. Planck Collaboration XXX 2014; Lenz et al. 2019; Stein et al. 2020). Nonetheless, future work will have to revisit these approximations. Potential improvements include the use of a multi-band analysis, such as in Chiang et al. (2020), where the signal in each frequency band is modelled separately, or the use of other CIB maps (e.g. Lenz et al. 2019).

Under the model of Eq. (9), the cross-power spectrum between the shear field and the observed Compton-y map is then given by

$$\begin{array}{c}\hfill C_{\ell }^{\gamma y_{\mathrm{obs}}}=C_{\ell }^{\gamma y_{\mathrm{true}}}+{\alpha }_{\mathrm{CIB}}{C}_{\ell }^{\gamma I_{\mathrm{CIB}}},\end{array}$$(10)

where $C_{\mathcal{l}}^{\gamma y_{\text{true}}}$ corresponds to the shear–tSZ cross-correlation model in the absence of CIB contamination and $C_{\mathcal{l}}^{\gamma I_{\text{CIB}}}$ is a model for the shear-CIB cross-correlation. The scaling parameter α_CIB can be found by comparing the measured cross-correlation between the Compton-y map and the CIB map with their auto-correlations. Here we used the 545 GHz CIB map from Planck Collaboration XLVIII (2016) as a proxy for I_CIB in Eq. (9) together with α_CIB = (2.3 ± 6.6)×10⁻⁷ MJy⁻¹ sr, as derived in Alonso et al. (2018).

To model $C_{\mathcal{l}}^{\gamma I_{\text{CIB}}}$ we chose a phenomenological approach, since modelling the CIB within the HMX framework is not possible at this point. Instead, we measured the tomographic cross-power spectra between the KiDS-1000 source samples and the 545 GHz CIB map (smoothed to a beam size of 10′ and masked with the y-map mask) and then fitted a Gaussian process (GP) to these observed cross-power spectra. The shear–CIB cross-spectra, together with the GP model and its uncertainty are shown in Fig. 3. The Gaussian data covariance was computed in an iterative process. First, the covariance was computed by NAMASTER using the cosmic shear model of Sect. 3.2 and the measured and smoothed auto-power spectrum of the CIB map, with zero cross-correlation. After fitting the GP model using this initial covariance, the covariance was recomputed, now including the GP model for the cross-correlation term. Finally, the GP model was refitted using the new covariance. The GP model was implemented using GPYTORCH¹⁰ (Gardner et al. 2018). The model is a multi-task GP (Bonilla et al. 2006), with a fixed data covariance. The fitted parameters consist of the length-scale of the squared exponential kernel, which accounts for the correlation between ℓ-modes, and the task covariance, which accounts for the correlation between the tomographic bins and is parameterised as the sum of a rank-one matrix and a diagonal component, for a total of 11 free parameters. Being a multi-task GP, the model thus accounts for the correlation between the redshift bins, as well as between ℓ-bins. Using a fixed data covariance ensures that the GP regression does not fit to spurious or excess (co)-variances. This formalism can be extended to account for the correlation between the shear–tSZ and shear–CIB cross-correlations, such that the GP model is conditioned on the shear–tSZ model. Due to the small effect of the CIB on the present analysis, we leave this sophistication to future work.

Fig. 3. Measured angular cross-power spectra between KiDS-1000 shear and the 545 GHz CIB map (blue points). The best-fit Gaussian process (GP) model is shown as the orange line, with the orange shaded band indicating the 68th percentiles of the GP model uncertainty.

To account for residual biases in our treatment of the CIB contamination, we marginalised over the α_CIB parameter in our fiducial analysis, using a Gaussian prior based on the value derived in Alonso et al. (2018) but with a variance that is four times larger. Alternatively, being a linear parameter, the marginalisation over α_CIB could have been accounted for analytically by an appropriate term in the data covariance (e.g. Bridle et al. 2002), analogous to the treatment of the multiplicative shear bias.

In our second approach to assess the sensitivity of our cross-correlation measurements to CIB contamination, we computed the shear–tSZ cross-correlation using CIB-subtracted y-maps from Yan et al. (2019). These maps subtract the CIB contributions in the 353, 545, and 845 GHz bands, based on the CIB maps constructed in Planck Collaboration XLVIII (2016). The CIB contributions at lower frequencies are extrapolated from these maps, using a grey body spectrum with the best-fit parameters from Planck Collaboration XXX (2014).

A different approach was used in the construction of the ACT Compton-y map (Madhavacheril et al. 2020) and the y-maps of the recent Dark Energy Survey Year 3 (DES Y3) cross-correlation analyses (Gatti et al. 2021; Pandey et al. 2021). During the ILC procedure, rather than projecting out the primary CMB, they deproject a grey body spectrum. The grey body spectrum is parameterised by the dust temperature T_d and the emissivity index β of the grey body spectrum, which modifies the slope of the black-body spectrum at low frequencies. Gatti et al. (2021) and Pandey et al. (2021) varied the emissivity index β between 1 and 1.4. This is significantly lower than the β = 1.75 ± 0.06 derived from the CIB power spectrum in Planck Collaboration XXX (2014) and which we used in our CIB-subtracted y-maps. Other works, such as the WebSky mocks (Stein et al. 2020) use the model parameters derived in Viero et al. (2013) based on Herschel CIB data, which prefer β = 1.6 ± 0.1. The effective dust temperature is T_d = 24 K in Gatti et al. (2021), while Yan et al. (2019) use an effective dust temperature at z = 1.2 of T_d = 32.4 K. However, at Planck frequencies the effect of the dust temperature and its redshift evolution is subdominant to the large changes in the emissivity index considered here. To check our CIB-subtracted Compton-y map against changes in β, we also generated a map that extrapolates the CIB to the lower frequency bands assuming a value of β = 1.2. In Sect. 4.2, we show that our results are robust with respect to all of these choices.

3.4.4. Radio sources

Unresolved radio sources are another potential contamination in the Compton-y map. Shirasaki (2019) suggests that unresolved radio sources can bias shear–tSZ cross-correlations by up to 30%. Their study however differs in key assumptions from ours. To calculate the expected contamination from unmasked radio sources, they assumed a flux cut of 400 mJy. In our analysis, we used the point-source and 40% Galactic mask provided in Planck Collaboration XXII (2016). The mask covers sources with fluxes > 300 mJy of the Second Planck Catalogue of Compact Sources (PCCS2; Planck Collaboration XXVI 2016), with none of the unmasked sources lying within the KiDS-1000 footprint. The flux cut assumed in Shirasaki (2019) therefore overestimates the number and flux from unmasked sources. Considering the sizeable bias estimated in Shirasaki (2019), even a lower flux cut might still give rise to a significant amount of radio contamination in our Compton-y maps. We therefore followed a similar approach to our treatment of the CIB contamination, namely estimate the effect on the shear–tSZ cross-correlation based on the cross-correlation of the shear field with the contaminant. To this end, we measured the angular cross-power spectrum between the KiDS-1000 source samples and the Planck 100 GHz map as a tracer of radio emissions. The cross-correlation with Gaussian errors is shown in Fig. 4. Unlike the cross-correlation with the CIB, we detected no significant correlation between the shear field and the 100 GHz map. We hence concluded that our estimates of the shear–tSZ cross-correlation are not significantly affected by contamination of our y-maps due to unresolved radio sources. This finding is also in agreement with Planck Collaboration XXII (2016), who find that radio sources are a subdominant contribution to the Compton-y auto-power spectrum, as well as with Gatti et al. (2021), who find that masking radio sources does not affect their shear–tSZ cross-correlation.

Fig. 4. Measured angular cross-power spectra between KiDS-1000 shear and the Planck 100 GHz intensity map.

3.5. Inference pipeline

Our inference pipeline is based on a modified version¹¹ of COSMOSIS¹² (Zuntz et al. 2015). The linear matter power spectrum is calculated using CAMB¹³ (Lewis et al. 2000). The non-linear matter and electron pressure power spectra are computed using HMX, after which the modelling proceeds as described in Sects. 3.2 and 3.4. The computational cost of the model prediction is dominated by the calculation of the non-linear power spectra, specifically the calculation of the halo profiles, since the matter and electron-pressure profiles have no analytic Fourier transforms, necessitating numerical transforms. To mitigate this, we evaluated the three-dimensional power spectra at eight redshifts between z = 0 and z = 1.5, which reduces the evaluation time of a single likelihood to ∼30 s. The smoothness of the power spectra and kernels in the integral of Eq. (4) limits the numerical error incurred by this coarse redshift sampling. We found that this redshift sampling is accurate to 1% of the measurement uncertainty for the shear–tSZ cross-correlation and to 0.5% for cosmic shear. We assumed a Gaussian likelihood, given by the model described above; a data vector consisting of either the cosmic shear or shear–tSZ data, as described in Sect. 2.3, or the concatenation of the two in the case of the joint analysis of cosmic shear and shear–tSZ; and the data covariance described in Sect. 3.3.

In our fiducial setup, we sampled 13 parameters, listed in Table 1, using the nested sampling implementation of MULTINEST (Feroz & Hobson 2008; Feroz et al. 2009, 2019), using 500 or 1000 live points, and an efficiency parameter of 0.3. The maximum of the posterior was found by running multiple Nelder-Mead (Nelder & Mead 1965) optimisations, starting from samples with the highest posterior values.

4. Results

Figure 5 presents the measurements of our fiducial data vector, consisting of the pseudo-C_ℓ estimates of the KiDS-1000 cosmic shear angular power spectra and the shear–tSZ angular cross-power spectra, using the Planck MILCA Compton-y map. Overplotted are the 68 and 95 percentiles of the posterior predictive distribution (PPD) and samples of the translated posterior distribution (TPD; Köhlinger et al. 2019) for the joint data vector, coloured by their value of log₁₀(T_AGN/K), which characterises the strength of AGN feedback in our model. The PPD describes fictional observations, conditioned on the data and model, and is a useful diagnostic to assess consistency within and between data vectors (e.g. Gelman et al. 2013; Köhlinger et al. 2019; Doux et al. 2021). The TPD translates the model parameter posterior into data-space, thus capturing the distribution of model predictions, conditioned on the data.

Fig. 5. Measured angular power spectra of KiDS-1000 cosmic shear (bottom left) and the angular cross-power spectra between KiDS-1000 shear and the tSZ effect (top right). Shown here is the cross-correlation using the MILCAPlanck Compton-y map. The grey bands indicate the 68th and 95th percentile of the posterior predictive distribution (PPD) of the joint-analysis posterior. Coloured lines are samples of the translated posterior distribution (TPD), with the colour corresponding to the baryon feedback strength. Since the TPD is derived from the full joint-posterior, the baryon feedback strength only accounts for part of the variance of the shown TPD samples.

The cosmic shear signal is detected at 20.9 σ signficance, while the shear–tSZ cross-correlation is detected at 9.4 σ significance. The significance was estimated as $\sqrt{{\mathit{d}}^{\mathrm{T}}{\mathsf{C}}^{-1}\mathit{d}}$, where d is the pseudo-C_ℓ data vector and C is the data covariance, given by Eq. (8). The B-modes of our measurements are consistent with zero, with details discussed in Appendix A.

Parameter constraints on Ω_m and σ₈ derived from the cosmic shear, shear–tSZ, and joint data are shown in Fig. 6. Our cosmic shear constraints exhibit the usual ‘cosmic banana’ degeneracy and are in good agreement with the KiDS-1000 band-power analysis of Asgari et al. (2021, cf. Appendix B). The cosmic shear signal scales approximately as $\propto {\sigma }_8^2{\mathrm{\Omega }}_{\text{m}}$ and is thus usually summarised with the parameter $S_8={\sigma }_8\sqrt{{\mathrm{\Omega }}_{\text{m}}/0.3}$. Depending on the specific estimator, the exact scaling can deviate slightly from S₈, such that the parameter ${\mathrm{\Sigma }}_8^{\alpha }={\sigma }_8{\left({\mathrm{\Omega }}_{\text{m}}/0.3\right)}^{\alpha }$, with an appropriately chosen α, is a better description of the direction orthogonal to the degeneracy. For example, the σ₈–Ω_m degeneracy for the band-power estimator in Asgari et al. (2021), as well as the pseudo-C_ℓ estimator in this work, is better described by ${\mathrm{\Sigma }}_8^{0.58}$ than S₈.

Fig. 6. Constraints from KiDS-1000 cosmic shear (blue), the cross-correlation of shear and the tSZ effect from Planck (orange), and their joint analysis (green) on Ωm and σ8.

The shear–tSZ cross-correlation is more sensitive to σ₈ (e.g. Refregier & Teyssier 2002; Komatsu & Seljak 2002; Mead et al. 2020) than cosmic shear, resulting in a degeneracy direction that is shallower in σ₈. We find that the cross-correlation signal scales approximately as $\propto {\sigma }_8^5{\mathrm{\Omega }}_{\text{m}}$. We use the parameter ${\mathrm{\Sigma }}_8^{0.2}$ when comparing different shear–tSZ cross-correlation constraints in the remainder of this analysis. In our fiducial analysis, which uses the Planck MILCA Compton-y map, we infer

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_8^{0.2}=0.{723}_{-0.032}^{+0.046}.\end{array}$$(11)

The different degeneracy directions of the cosmic shear and shear–tSZ cross-correlations allow the joint analysis of the two to break these degeneracies and provide tight cosmological constraints. We find

$$\begin{array}{cc}\hfill {\mathrm{\Omega }}_{\mathrm{m}}& =0.{342}_{-0.037}^{+0.042}\hfill \\ \hfill {\sigma }_8& =0.{692}_{-0.039}^{+0.05}\hfill \\ \hfill S_8& =0.{751}_{-0.017}^{+0.02}\hfill \end{array}$$(12)

for the maximum of the marginal posterior and marginal highest posterior density credible intervals (CIs). The S₈ constraints are improved by a factor of 40% over cosmic shear alone and are in excellent agreement with the KiDS-1000 3 × 2pt analysis of Heymans et al. (2021). This improvement is driven by the breaking of the residual degeneracy between S₈ and Ω_m, since S₈ does not optimally capture the direction orthogonal to the σ₈–Ω_m degeneracy for the estimator considered here. Using instead ${\mathrm{\Sigma }}_8^{0.58}$, as in Asgari et al. (2021), which is better constrained than S₈ by the pseudo-C_ℓ estimator, we find a 15% improvement of the constraining power of the shear–tSZ + cosmic shear joint analysis over cosmic shear alone.

The joint analysis also strongly tightens the constraints on Ω_m and σ₈. This is to be expected, however, since the cosmic shear constraints on these parameters are prior-dominated, and any breaking of the cosmic banana will yield comparatively tight constraints compared to cosmic shear alone (see e.g. Tröster et al. 2021, for examples of this effect in the joint analyses of KiDS-1000 cosmic shear with CMB lensing or super novae).

The shear–tSZ cross-correlation prefers lower values of σ₈ than those inferred from the CMB by Planck; the discrepancy between our inferred value of ${\mathrm{\Sigma }}_8^{0.2}$ and that of Planck Collaboration VI (2020; using the TTTEEE+lowE likelihood) is 2.2 σ. This trend is also seen in the recent shear–tSZ cross-correlation analysis with DES Y3 data (Gatti et al. 2021), although they impose Planck priors on σ₈ and Ω_m, making a direct comparison with the results in this work difficult. The constraints on S₈ and ${\mathrm{\Sigma }}_8^{0.58}$ inferred from the joint analysis are in 3.4 σ and 2.9 σ tension with those of Planck, similar to previous joint analyses of KiDS-1000 with other datasets (Heymans et al. 2021; Tröster et al. 2021), as well as KiDS-1000 cosmic shear by itself (Asgari et al. 2021). The full parameter-space constraints are shown in Fig. C.1.

We note here that the considerable challenges of the accurate and consistent modelling of the cross-correlation signal warrant revisiting these results in the future, when the modelling has matured and more data have become available. In particular, accounting for the dependences of, and interactions between, cosmology, baryon feedback, and IA requires both models that consistently include these effects and simulations that sufficiently cover this parameter space to validate said models. These points are further discussed in Sect. 4.3.

4.1. Goodness of fit

We find that the data are well described by our model. The χ² of the pseudo-C_ℓ cosmic shear data with respect to the model of the maximum a posteriori (MAP) estimate is 150.7, with 117.1effective degrees of freedom. The effective model dimensionality was estimated following Raveri & Hu (2019), and accounts for the fact that many of the varied parameters in the model are prior-dominated and as such do not contribute extra degrees of freedom to the model. Our fiducial shear–tSZ cross-correlation has a χ² of 28.0, with 36.2effective degrees of freedom, while the joint analysis has a χ² of 179.3, with 155.4effective degrees of freedom. The goodness-of-fit values of our measurements are summarised in Table 2.

4.2. Alternative Compton-y maps and effect of the CIB

To assess the degree to which y-map reconstruction and CIB contamination treatment affect our results, we derived constraints using a range of Compton-y maps and compared the resulting Ω_m − σ₈ and ${\mathrm{\Sigma }}_8^{0.2}$ constraints to those derived based on our fiducial setup, consisting of the Planck MILCAy-map and CIB marginalisation using the GP model described in Sect. 3.4.3. The Ω_m − σ₈ constraints based on the other Compton–y maps are shown in Fig. 7. The constraints on ${\mathrm{\Sigma }}_8^{0.2}$ are compared in Fig. 8.

Fig. 7. Effect of different Compton-y map reconstruction methods on the shear–tSZ cross-correlation constraints on Ωm and σ8. Top: Compton-y maps derived from Planck data. Bottom: Compton-y maps from ACT.

Fig. 8. Marginal posterior densities of ${\mathrm{\Sigma }}_8^{0.2}={\sigma }_8{\left({\mathrm{\Omega }}_{\text{m}}/0.3\right)}^{0.2}$ from the cross-correlation of KiDS and different Compton-y maps. The orange band and dashed line indicate the marginal 68th percentile CI and mode of the marginal posterior of the fiducial setup.

The Planck-derived maps yield tighter constraints than those from ACT, owing to the larger overlap between the KiDS-1000 and Planck data. Among the ACT maps, we find that the y-map where no other component has been deprojected leads to weaker constraints, with a tail to low values of ${\mathrm{\Sigma }}_8^{0.2}$. We explain this effect with a slightly larger data covariance due to a higher auto-power spectrum, as well as anomalously low data points in the shear–tSZ cross-correlation of the Compton-y map with no deprojected component compared to the maps with either the CMB or CIB projected out. We discuss these points in more detail in Appendix D. Both the CMB-deprojected and CIB-deprojected y-maps give very similar constraints, with no systematic shift between them. If the cross-correlation was significantly affected by CIB contamination, we would expect to find different constraints in these two cases.

We find good agreement between all maps and no evidence that the CIB causes significant biases in our constraints. This is in contrast to the finding in Gatti et al. (2021), who find that deprojecting a grey body spectrum during the construction of their Planck-derived Compton-y map leads to significantly different cross-correlation signals compared to the case where no such deprojection was applied. The origin of this discrepancy is currently unknown, however, and will require a detailed comparison of the different Compton-y maps. We leave this for future work.

The mask of our Compton-y maps is potentially correlated with the lensing signal, as the point source mask includes some clusters. Since these clusters are sources of both tSZ and weak lensing signals, this could in principle cause the estimated signal to be biased low. As the number of masked point sources in the KiDS-1000 footprint is low, this effect is likely not relevant in our case.

4.3. Baryon feedback

A key goal of this analysis has been to constrain baryonic feedback through the gas physics probed by the tSZ effect. In our model, the effect of baryonic feedback, specifically AGN feedback, is parameterised by the parameter log₁₀(T_AGN/K). The fiducial BAHAMAS model assumes log₁₀(T_AGN/K) = 7.8, since this value reproduces the observed present-day galaxy stellar mass function and the hot gas mass fractions of groups and clusters in the simulations (McCarthy et al. 2017). The constraints on this parameter, as inferred from cosmic shear alone, the shear–tSZ cross-correlation, and the joint analysis are shown in Fig. 9. While the KiDS-1000 cosmic shear band-power and pseudo-C_ℓ data by themselves have weak preference for low feedback strength (Asgari et al. 2021; Tröster et al. 2021), the shear–tSZ cross-correlation data prefer strong feedback. For the joint analysis we find

$$\begin{array}{c}\hfill {log}_{10}\left(\frac{T_{\mathrm{AGN}}}{\mathrm{K}}\right)=7.{96}_{-0.48}^{+0.2},\end{array}$$(13)

Fig. 9. Marginal posterior densities of the baryon feedback strength, parameterised by log10(TAGN/K). The green band and dashed line indicate the 68th percentile marginal CI and MAP, respectively, of the joint analysis of cosmic shear and the shear–tSZ cross-correlation. The grey lines indicate the range of the feedback strength implemented in the BAHAMAS suite of hydrodynamical simulations, which our HMX model is calibrated on.

in good agreement with the BAHAMAS model. We note that this constraint does not strictly pass the criterion of Asgari et al. (2021) for a parameter to count as constrained by the data: the marginal posterior density at the (lower) prior edge is not less than 13% of the posterior density at the marginal mode.

This result should be interpreted with caution, however. The model for the matter and electron pressure power spectra underlying our analysis, HMX, was calibrated against the BAHAMAS suite of hydrodynamical simulations. The range of feedback strength considered in BAHAMAS ranges from 7.6 to 8.0, indicated by the grey band in Fig. 9, and was chosen such that the simulations reproduce the observed galaxy stellar mass function and gas fraction in groups and clusters. Outside of this range the model extrapolates from the simulations it was calibrated on. While the physical nature of the model makes this extrapolation viable to a certain degree, there is significant modelling uncertainty at the edges of the log₁₀(T_AGN/K) ∼ U(7.1, 8.5) prior we use in our inference pipeline. Furthermore, the fitting of the HMX model parameters and their relation to log₁₀(T_AGN/K) is done at a single, fixed cosmology. It has been shown that the effects of cosmology and galaxy formation are largely independent (e.g. van Daalen et al. 2011; Mummery et al. 2017; van Daalen et al. 2020; Pfeifer et al. 2020). Residual interplay between feedback parameters and cosmology, for example through the baryon fraction, can affect the inferred values of log₁₀(T_AGN/K) to a small degree, however (Mead et al. 2021). To address these shortcomings, improvements in the modelling and its calibration will require suites of simulations that cover a wider range of cosmologies and feedback parameters. These simulations also need to cover larger volumes in order to suppress the sampling variance of the electron pressure power spectra at large scale (cf. Sect. 3.1 and Mead et al. 2020).

The shear–tSZ cross-correlations using the ACT maps do not constrain log₁₀(T_AGN/K). This can be attributed to the generally weaker constraints of the ACT maps due to the much smaller overlap with the KiDS-1000 data and the associated statistical scatter.

While the shear–tSZ cross-correlation by itself prefers strong feedback, the joint constraints on the feedback parameter agree with the BAHAMAS model. This shift in the preferred feedback strength is due to the interaction of IA and baryon feedback: both have the effect of a broadband suppression of the shear–tSZ cross-power spectra, such that high IA amplitudes can be compensated to a certain degree by weak baryon feedback. The joint analysis with cosmic shear restricts the range of allowed IA amplitudes, excluding the parts of parameter space with strong feedback. The IA model used in this analysis is quite simplistic, and the HMX model is getting extrapolated to feedback strengths outside the range it is calibrated on. To better understand this IA–baryon feedback interaction, especially away from their fiducial parameters, a more holistic model is required, for example through including the halo-model IA formalism of Fortuna et al. (2021) in HMX, as well as expanding the calibration set of hydrodynamical simulations.

One might have expected the shear–tSZ cross-correlation to yield tighter constraints on baryon feedback, since it is sensitive to both the effect of AGN feedback on gas physics and its suppression of the matter power spectrum. At large scales, k ≲ 0.3 h Mpc⁻¹, stronger feedback increases the matter-electron pressure cross-power spectrum due to heating of the warm-hot intergalactic medium. These scales are largely outside of our scale cuts of ℓ ∈ (100, 1500), however. At small scales, k ≳ 0.3 h Mpc⁻¹, increasing the feedback strength suppresses the matter-electron pressure cross-power spectrum (see Fig. 9 in Mead et al. 2020). For the scales considered in this analysis, the effect of AGN feedback is thus largely restricted to a suppression of the angular cross-power spectrum over a large range of multipoles, similar to its effect on cosmic shear. This is illustrated in Fig. 5, where the colouring of the samples of the TPD indicates the effect of log₁₀(T_AGN/K) on the model prediction: except at scales of ℓ ∼ 100, the effect of log₁₀(T_AGN/K) is limited to a suppression of the angular power spectra. Conversely, ongoing and future large-area surveys, such as DES Y3, LSST, or Euclid, have access to larger angular scales and thus could make use of the different behaviour at small and large scales to derive tighter constraints on baryon feedback.

5. Conclusion

We conducted a pseudo-C_ℓ cosmic shear analysis of the KiDS-1000 weak gravitational lensing data and find excellent agreement with previous KiDS-1000 angular power spectrum analyses (Asgari et al. 2021; Loureiro et al. 2021). We then measured the cross-correlation of KiDS-1000 weak lensing and the tSZ effect, detecting the cross-correlation at a 9.4 σ significance.

Using the HMX model, which allows for the consistent modelling of both cosmic shear and the shear–tSZ cross-correlation – including the effect of baryon feedback – we derived constraints on cosmology and the baryon feedback parameter from the shear–tSZ cross-correlation. While the inferred value of the amplitude of the shear–tSZ signal, as parameterised by ${\mathrm{\Sigma }}_8^{0.2}$, is fully consistent with other KiDS analyses, it is low by 2.2 σ compared to CMB constraints from Planck. Finally, we conducted a joint analysis of our cosmic shear and shear–tSZ cross-correlation data, taking the full cross-covariance between these data vectors into account. By breaking degeneracies, the inferred cosmological parameter constraints of the joint analysis are tighter than those of either probe by itself. For example, the constraints on S₈ of the joint analysis improve by 40% over cosmic shear by itself. We find that our cross-correlation data prefer strong baryon feedback, while we infer log₁₀(T_AGN/K) = ${7.96}_{-0.48}^{+0.2}$ from the joint analysis, consistent with the BAHAMAS hydrodynamical simulations. Unlike recent work by Gatti et al. (2021) on shear–tSZ cross-correlation using DES Y3 data, we find no evidence of significant CIB contamination in our Planck-derived Compton-y maps.

While the analysis presented here demonstrates the power of shear–tSZ cross-correlations and their joint analyses with cosmic shear to constrain both cosmology and baryonic processes of galaxy formation, there are still significant challenges in the detailed interpretation of these results. The modelling of the non-linear matter spectrum to the accuracies required by current and future weak-lensing surveys is already a challenging proposition. Accounting for the electron-pressure field required to predict the tSZ effect adds a further level of complexity to the modelling framework. While the consistent modelling of various probes in HMX is a significant improvement on earlier works that treat the profiles of different components as independent, its current accuracy for the matter–electron pressure power spectrum is limited at the 15 percent level and will need to be improved for future high-precision measurements, for example through the inclusion of non-linear halo biasing in the core halo model (Mead & Verde 2021). Improvements in the modelling must also go hand in hand with improvements of the simulation suites used for calibration, both in their coverage of the cosmology and baryon feedback parameter spaces and in their volumes to suppress sampling variance.

If these challenges can be overcome, the joint analysis of weak lensing from forthcoming galaxy surveys, such as LSST and Euclid, with low-noise and high-resolution measurements of the tSZ effect from future CMB experiments, such as the Simons Observatory (Ade et al. 2019), can be used to extract more information on galaxy formation processes and cosmology from these datasets than when analysed independently. Another avenue would be to combine the measurements considered here with other datasets. The combination of weak lensing, the clustering of galaxies, and their cross-correlations into ‘3 × 2pt’ analyses has provided some of the tightest constraints on cosmology from the low-redshift Universe (e.g. DES Collaboration 2018, 2022; Heymans et al. 2021). A lensing–tSZ 3 × 2pt analysis, through inclusion of the tSZ auto-correlation in the measurements considered in this work, would thus be a logical next step and would build on previous work of joint analyses of the tSZ auto-power spectrum and lensing (Osato et al. 2020; Makiya et al. 2020). Beyond the shear–tSZ cross-correlation considered in this work, cross-correlations of galaxies with both the tSZ effect and the kinetic Sunyaev–Zeldovich effect have been shown to be a potent probe of the thermodynamics of haloes (e.g. Battaglia et al. 2017; Schaan et al. 2021; Amodeo et al. 2021). Adding these various Sunyaev–Zeldovich-effect cross-correlations into the framework of 3 × 2pt analyses could therefore improve the constraining power of the data that are already available.

Acknowledgments

The figures in this work were created with MATPLOTLIB (Hunter 2007) and GETDIST (Lewis 2019), making use of the NUMPY (Oliphant 2006) and SCIPY (Jones et al. 2001) software packages. This project has received significant funding from the European Union’s Horizon 2020 research and innovation programme. We thank and acknowledge support from: the European Research Council under grant agreement No. 647112 (T.T., A.M., C.H., and M.A.) and No. 770935 (H.Hi., A.D., and A.H.W.) in addition to the Marie Skłodowska-Curie grant agreements No. 797794 (T.T.) and No. 702971 (A.M.). We also acknowledge support from the Leverhulme Trust (T.T.); 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 (C.H.); the University of British Columbia, Canada’s NSERC, and CIFAR (Z.Y. and L.V.W.); the Beecroft Trust and the Science and Technology Facilities Council through an Ernest Rutherford Fellowship, grant reference ST/P004474 (D.A.); the Polish Ministry of Science and Higher Education through grant DIR/WK/2018/12, and the Polish National Science Center through grants no. 2020/38/E/ST9/00395, 2018/30/E/ST9/00698, and 2018/31/G/ST9/03388 (M.B.); the Deutsche Forschungsgemeinschaft Heisenberg grant Hi 1495/5-1 (H.Hi.); the Royal Society and Imperial College (K.K.); the NSFC of China under grant 11973070, the Shanghai Committee of Science and Technology grant No.19ZR1466600, the Key Research Program of Frontier Sciences, CAS, Grant No. ZDBS-LY-7013, and CMS-CSST-2021-A01 (HYS). Based on data products from observations made with ESO Telescopes at the La Silla Paranal Observatory under programme IDs 177.A-3016, 177.A-3017 and 177.A-3018. Author contributions: All authors contributed to the development and writing of this paper. The authorship list is given in three groups: the lead authors (T.T., A.M., C.H.) 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: B-modes

Figure A.1 presents the B-modes of our measurements, which are expected to be consistent with zero in the absence of systematics. To be more specific, for cosmic shear the γ_Bγ_B component of the angular power spectra is shown, while for the shear–tSZ cross-correlation, the γ_By component is depicted. Both pseudo-C_ℓ measurements are consistent with zero, with ${\chi }_{{\gamma }_B{\gamma }_B}^2=128.3$ and ${\chi }_{{\gamma }_{B}y}^2=36.5$ for 120 and 40 degrees of freedom, respectively. This translates to probabilities to exceed (PTEs) of 0.29and 0.63, which makes us confident in the robustness of our measurements. The same is true for the other y-maps considered in this work. Furthermore, the B-modes remain consistent with zero at the same level of significance when considering the whole ℓ ∈ (51, 2952) range they were measured on, rather than the ℓ ∈ (100, 1500) range used for parameter inference.

Fig. A.1. B-modes of the measurements considered in this work. The cosmic shear BB component is shown on the left, for which we find ${\chi }_{{\gamma }_B{\gamma }_B}^2=128.3$, corresponding to a PTE of 0.29. The B-mode of the shear–tSZ cross-correlation is shown in the upper right. It has ${\chi }_{{\gamma }_{B}y}^2=36.5$, for a PTE of 0.63.

Appendix B: Comparison with previous KiDS-1000 cosmic shear results

The cosmic shear analysis presented here closely followed that of Asgari et al. (2021) but differs in two aspects: the choice of angular power spectrum estimator and the modelling of baryonic feedback.

B.1. Comparison of the band-power and pseudo-C_ℓ estimators

Different two-point statistics estimators can give slightly different results (e.g. Hikage et al. 2019; Hamana et al. 2020; Asgari et al. 2021), as they are sensitive to different scales and are thus differently affected by statistical fluctuations in the underlying data. Here we compared our pseudo-C_ℓ estimates and derived cosmological constraints to those from the band-power estimator used in Asgari et al. (2021). That band-power estimator is based on the integral transform of the two-point correlation function (e.g. Schneider et al. 2002; Becker & Rozo 2016; Tröster et al. 2017; van Uitert et al. 2018; Joachimi et al. 2021) and has attractive properties, such as being insensitive to the footprint geometry nor requiring subtraction of the noise bias.

We find excellent agreement between the measured data vectors, illustrated in Fig. B.1. Since both are estimators for the angular power spectrum, this is to be expected to a certain degree but it is reassuring the see this level of agreement. The uncertainties differ somewhat between the estimators, likely due to differences in the details of the data covariance calculation, for example, in the treatment of the mixed term of the Gaussian contribution. The window functions for the EE-component that translate the model prediction into data-space, accounting for example for residual effects of the mode-mixing matrices and E/B-mode mixing, are shown in Fig. B.2 for the two estimators. The largest differences are at large scales, where mask effects and the finite angular range on which the correlation function were measured become important.

Fig. B.1. Comparison of KiDS-1000 angular power spectrum estimates. The pseudo-Cℓ estimates used in this work are shown in orange, and the band-power estimates of Asgari et al. (2021) are shown in blue.

Fig. B.2. Window functions of Eq. (7), relating the smooth angular power spectrum model to the binned prediction compatible with the measured data vector. The differences between the multipole bin (grey bands) and the window function are due to masking and weighting of the shear maps in the case of the pseudo-Cℓ estimator (orange) and the finite angular range used for the band-power estimator (blue).

We derived cosmological constraints from both data vectors, taking care to homogenise the inference pipeline as much as possible. We use the same HMCODE-2016 model (Mead et al. 2016) for the matter power spectrum, as well as the multiplicative shear bias values from Giblin et al. (2021) used in Asgari et al. (2021). The marginal posteriors in the Ω_m-S₈ plane and for S₈ are shown in Figs. B.3 and B.4, demonstrating good agreement between the two estimators. Due to differences in the way the data covariances for the two estimators were computed, small differences in the parameter constraints are to be expected.

Fig. B.3. Comparison of the constraints on Ωm and S8 between the band-power analysis of Asgari et al. (2021) and the pseudo-Cℓ analysis of this work (top) and effect of different non-linear matter power spectrum modelling approaches on the cosmology constraints of the pseudo-Cℓ cosmic shear measurements (bottom).

Fig. B.4. Marginal posterior densities of S8, for the same data and models as in Fig. B.3. The blue band and dashed line indicate the 68th percentile marginal CI and MAP of the Asgari et al. (2021) band-power analysis, respectively.

Our methodology differs in a subtler manner from the pseudo-C_ℓ methodology of Loureiro et al. (2021), however, even though both approaches make use of the pseudo-C_ℓ formalism. In the present work the data vector was deconvolved with the mode-mixing matrices, the mode-mixing matrices were calculated using the weight map, and we used analytic estimates for the noise bias and data covariance. On the other hand, Loureiro et al. (2021) forward-modelled the effect of the mode mixing, used a binary mask to compute the mode-mixing matrices, and estimated their noise bias and data covariance from simulated mock data.

B.2. Impact of baryon modelling

The second significant difference between the analysis choices of the present work and Asgari et al. (2021) is in the non-linear model for the matter power spectrum. Asgari et al. (2021) used HMCODE-2016 (Mead et al. 2016), which was calibrated on the OWLS suite of hydrodynamical simulations. The HMX model used here has been calibrated against the BAHAMAS suite of simulations instead. Importantly, it models gas, stars, and dark matter separately. The model therefore does not have a dark-matter-only limit, in contrast to HMCODE-2016.

Figures B.3 and B.4 show the effects of the different matter power spectrum model choices, comparing the HMCODE-2016 and HMX models. The comparison also includes the HMCODE-2020 (Mead et al. 2021) model, which has been optimised for the modelling of the matter power spectrum. Like HMX, the baryonic feedback has been calibrated on BAHAMAS but unlike HMX, it only predicts the matter power spectrum and thus cannot model the cross-correlations considered in the present analysis. As discussed in appendix A of Tröster et al. (2021), HMCODE-2020 allows for slightly higher values of S₈ because the finite amount of AGN feedback in the model always causes a degree of suppression of power. The same effect can be observed for HMX, for the same reasons: the physical model for gas in the model prevents it from behaving as if it were dark matter.

Appendix C: Other parameter constraints

Figure C.1 shows the marginal posteriors for all 13 sampled parameters in our fiducial cosmic shear, shear–tSZ, and joint analyses.

Fig. C.1. Constraints on all sampled parameters from KiDS-1000 cosmic shear, the cross-correlation of shear and the tSZ effect, and their joint analysis.

Appendix D: ACT measurements

Figure D.1 shows the measured shear–tSZ angular cross-power spectra for different ACT Compton-y maps provided by Madhavacheril et al. (2020):

• tilec_single_tile_BN_comptony_map_v1.2.0_joint, corresponding to the case where no extra component has been deprojected;

• tilec_single_tile_BN_comptony_deprojects_cmb_map_v1.2.0_joint, where the primary CMB has been deprojected;

• tilec_single_tile_BN_comptony_deprojects_cib_map_v1.2.0_joint, where a CIB grey body spectrum has been deprojected. (See Sects. 3.4.3 for further discussion of the CIB treatment in the ACT maps.)

Fig. D.1. Measured angular cross-power spectra between KiDS-1000 shear and the ACT DR4 Compton-y maps.

The auto-power spectrum of the ACT y map with no deprojection is higher than that of the CMB and CIB deprojected maps on the ℓ-range we used in our analysis. This leads to 5-15 % larger errors on the shear–tSZ cross-correlation compared to the maps with deprojected components. The increase in the data covariance by itself does not explain the significantly weaker parameter constraints seen in Fig. 7, however. Furthermore, the best-fit χ² is higher than those of the maps with deprojected components, even though the data covariance is larger. Inspection of the measured pseudo-C_ℓ, shown in Fig. D.1, reveals that the highest redshift bin exhibits two data points that appear low, considering the strong correlation between the maps. To asses the impact of these data points, we analysed the shear–tSZ pseudo-C_ℓ for the three ACT maps with the range ℓ ∈ (387.3, 1069.3) excised for the two highest redshift bins. This removes the tail to low values in the ${\mathrm{\Sigma }}_8^{0.2}$ posterior and brings all constraints into good agreement. Since we have no a priori reason to reject these data points, we refrain from further discussion of the constraints from these restricted data vectors. These findings warrant a more in-depth analysis of the ACT Compton-y maps, which we leave for future work.

Appendix E: PYHMCODE

We provide an easy to use and unified Python interface to HMCODE and HMX, called PYHMCODE¹⁴. A usage example is shown below. It assumes a CCL cosmology object ccl_cosmology and arrays of the linear matter power spectrum pofk_lin, redshift z, and wave number k (in units of  h Mpc⁻¹) and computes the auto- and cross-power spectra between the matter and electron-pressure fields.

import pyhmcode
from pyhmcode.halo_profile_utils import
    ccl2hmcode_cosmo

# Create the pyhmcode cosmology object.
hmcode_cosmology = ccl2hmcode_cosmo(
    ccl_cosmo=ccl_cosmology,
    pofk_lin_k_h=k,
    pofk_lin_z=z,
    pofk_lin=pofk_lin,
    log10_T_heat=7.8)

# Create the halo model object, which holds
# information on the specific halo model
# to use. E.g., the HMCode or HMx version.
hmcode_model = pyhmcode.Halomodel(
    pyhmcode.HMx2020_matter_pressure_w_temp_scaling)

# Now we can compute the non-linear power spectrum,
# given the cosmology, halo model,
# and a list of fields.
hmcode_pofk = pyhmcode.calculate_nonlinear_power_spectrum(
    cosmology=hmcode_cosmology,
    halomodel=hmcode_model,
    fields=[pyhmcode.field_matter,
            pyhmcode.field_electron_pressure])

All Tables

All Figures

Fig. 1. Top: redshift distributions ni(z) of the five tomographic source samples (solid lines). The shaded bands indicate the selection window on the photometric redshift point estimate zB. Bottom: line-of-sight window function, W(z), to project the three-dimensional power spectra P(k) (in comoving units) into their angular counterparts, Cℓ. The lensing window functions (solid lines, Eq. (5)) are normalised by the maximum of the window function of the last tomographic bin, while the Compton-y window function (dashed line, Eq. (6)) is normalised by its maximum value.

In the text

Fig. 2. Joint cosmic shear and shear–tSZ data covariance. Upper right: contribution of the super sample covariance (SSC) to the joint covariance. Shown is the ratio between the absolute value of the SSC contribution and the absolute value of the total covariance in Eq. (8). The ratio can exceed unity, since different contributions to total covariance can cancel one another out to a certain degree. Lower left: correlation coefficients corresponding to the data covariance matrix.

In the text

	Fig. 3. Measured angular cross-power spectra between KiDS-1000 shear and the 545 GHz CIB map (blue points). The best-fit Gaussian process (GP) model is shown as the orange line, with the orange shaded band indicating the 68th percentiles of the GP model uncertainty.
In the text

	Fig. 4. Measured angular cross-power spectra between KiDS-1000 shear and the Planck 100 GHz intensity map.
In the text

Fig. 5. Measured angular power spectra of KiDS-1000 cosmic shear (bottom left) and the angular cross-power spectra between KiDS-1000 shear and the tSZ effect (top right). Shown here is the cross-correlation using the MILCAPlanck Compton-y map. The grey bands indicate the 68th and 95th percentile of the posterior predictive distribution (PPD) of the joint-analysis posterior. Coloured lines are samples of the translated posterior distribution (TPD), with the colour corresponding to the baryon feedback strength. Since the TPD is derived from the full joint-posterior, the baryon feedback strength only accounts for part of the variance of the shown TPD samples.

In the text

	Fig. 6. Constraints from KiDS-1000 cosmic shear (blue), the cross-correlation of shear and the tSZ effect from Planck (orange), and their joint analysis (green) on Ωm and σ8.
In the text

	Fig. 7. Effect of different Compton-y map reconstruction methods on the shear–tSZ cross-correlation constraints on Ωm and σ8. Top: Compton-y maps derived from Planck data. Bottom: Compton-y maps from ACT.
In the text

	Fig. 8. Marginal posterior densities of ${\mathrm{\Sigma }}_8^{0.2}={\sigma }_8{\left({\mathrm{\Omega }}_{\text{m}}/0.3\right)}^{0.2}$ from the cross-correlation of KiDS and different Compton-y maps. The orange band and dashed line indicate the marginal 68th percentile CI and mode of the marginal posterior of the fiducial setup.
In the text

Fig. 9. Marginal posterior densities of the baryon feedback strength, parameterised by log10(TAGN/K). The green band and dashed line indicate the 68th percentile marginal CI and MAP, respectively, of the joint analysis of cosmic shear and the shear–tSZ cross-correlation. The grey lines indicate the range of the feedback strength implemented in the BAHAMAS suite of hydrodynamical simulations, which our HMX model is calibrated on.

In the text

	Fig. A.1. B-modes of the measurements considered in this work. The cosmic shear BB component is shown on the left, for which we find ${\chi }_{{\gamma }_B{\gamma }_B}^2=128.3$, corresponding to a PTE of 0.29. The B-mode of the shear–tSZ cross-correlation is shown in the upper right. It has ${\chi }_{{\gamma }_{B}y}^2=36.5$, for a PTE of 0.63.
In the text

	Fig. B.1. Comparison of KiDS-1000 angular power spectrum estimates. The pseudo-Cℓ estimates used in this work are shown in orange, and the band-power estimates of Asgari et al. (2021) are shown in blue.
In the text

Fig. B.2. Window functions of Eq. (7), relating the smooth angular power spectrum model to the binned prediction compatible with the measured data vector. The differences between the multipole bin (grey bands) and the window function are due to masking and weighting of the shear maps in the case of the pseudo-Cℓ estimator (orange) and the finite angular range used for the band-power estimator (blue).

In the text

	Fig. B.3. Comparison of the constraints on Ωm and S8 between the band-power analysis of Asgari et al. (2021) and the pseudo-Cℓ analysis of this work (top) and effect of different non-linear matter power spectrum modelling approaches on the cosmology constraints of the pseudo-Cℓ cosmic shear measurements (bottom).
In the text

	Fig. B.4. Marginal posterior densities of S8, for the same data and models as in Fig. B.3. The blue band and dashed line indicate the 68th percentile marginal CI and MAP of the Asgari et al. (2021) band-power analysis, respectively.
In the text

	Fig. C.1. Constraints on all sampled parameters from KiDS-1000 cosmic shear, the cross-correlation of shear and the tSZ effect, and their joint analysis.
In the text

	Fig. D.1. Measured angular cross-power spectra between KiDS-1000 shear and the ACT DR4 Compton-y maps.
In the text