Cross-correlation of the thermal Sunyaev–Zel’dovich and CMB lensing signals in Planck PR4 data with robust CIB decontamination

,


I. INTRODUCTION
The cosmic microwave background (CMB) photons have been travelling through the Universe almost since it began, having been released about 380,000 years after the Big Bang.They carry valuable information not only about the state of the Universe at that time, but also about the structures they have encountered along the way, including information about the growth of structure and astrophysical signals.These late-Universe effects, referred to as CMB "secondary" anisotropies, have long been appreciated as late-Universe large-scale structure (LSS) probes.CMB lensing and the Sunyaev-Zel'dovich effect are particularly well-known examples.
CMB lensing -the anisotropies induced by the weak gravitational lensing deflection of CMB photons as they pass near matter overdensities and underdensities -is an extremely well-understood LSS probe and traces the growth of all of the post-recombination structure in the Universe (along our past lightcone).First detected statistically in cross-correlation with galaxy surveys [3] using WMAP data, it has been detected at high significance by the Atacama Cosmology Telescope (ACT), the South Pole Telescope (SPT), and Planck [4][5][6][7][8].CMB lensing probes matter at nearly all redshifts, but is particularly sensitive to structure in the range z ∼ 0.5 − 5. CMB lensing is usually detected by reconstructing a convergence (κ) map from the observed CMB temperature and polarization anisotropies using optimal quadratic estimators (see, e.g., [9,10]), which exploit the well-understood non-Gaussian statistics induced in a Gaussian primary CMB by lensing.
The Sunyaev-Zel'dovich (SZ) effect [11,12] is sourced by the Compton scattering of CMB photons from free electrons in the late Universe.These are electrons which fill the Universe after "reionization" at z ∼ 7 − 10, when the radiation from the first stars reionized the hydrogen atoms which had been neutral since recombination (the release of the CMB).In particular, the thermal SZ (tSZ) effect is sourced when the electron in the interaction has a high temperature and up-scatters the CMB photon, changing its frequency and sourcing a spectral distortion in the CMB.
The tSZ anisotropies are subdominant with respect to the primary CMB anisotropies, but the unique tSZ frequency dependence can be exploited to separate the primary CMB and the anisotropies induced by the tSZ effect using multi-frequency measurements.While it was first detected by making targeted observations of clusters in millimeter wavelengths [13,14], it is now used to discover clusters at high signal-to-noise [15][16][17] and can be studied at the map (or power spectrum) level by making an all-sky y-map [2,[18][19][20][21][22].As the electrons that source the tSZ effect must have very high temperature, it is mostly a very low-z probe (z < ∼ 1), as the massive clusters with the highest temperatures (which dominate the signal) take nearly a Hubble time to form.
CMB lensing and the tSZ effect provide complementary probes of cosmology.As a very well-theoreticallyunderstood, mostly linear signal, CMB lensing can be used to robustly infer cosmological parameters [7,23,24] via the power spectrum of a lensing convergence map, which is a direct probe of the (redshift-integrated) matter power spectrum.In contrast, the tSZ effect is sensitive to cosmology primarily as a probe of the halo mass function (HMF), which quantifies the abundance of collapsed objects as a function of their mass.In particular, as the tSZ effect probes the highest-mass objects (and thus the rarest density-field peaks) in the Universe, it probes directly the high-mass tail of the HMF, which is especially sensitive to cosmology.This information can be accessed either by making a cluster catalogue from high signal-to-noise ratio tSZ-selected objects or through measuring the power spectrum (or various other N -point statistics, including the bispectrum/skewness [25][26][27] or 1-point PDF [28,29]) of an all-(or partial-)sky tSZ map.Simultaneously, the tSZ effect probes the astrophysics of the intracluster medium (ICM) in which the electrons live, probing the pressure profile of galaxy clusters.Any cosmology constraint derived from the analysis of a single statistic of the tSZ field is thus degenerate with astrophysical uncertainties.
As the tSZ effect is primarily sourced by objects at lower redshifts than those that dominate the CMB lensing redshift kernel, the two signals contain significant independent information, with a correlation coefficient (on perfectly measured signals) predicted to be around 30-40% (for current tSZ models) [1].Their cross-correlation weights the higher-z contribution to the tSZ signal higher than in the tSZ auto-power spectrum, and thus probes higher-z ICM physics (and thus physics of clusters with lower mass) than the tSZ signal alone.In addition, because somewhat lower-mass halos (which are much more common than the highest-mass halos) contribute most of the CMB lensing signal, the tSZ-lensing cross-correlation probes lower-mass halos than the tSZ auto-correlation does.Importantly, the degeneracy between ICM astrophysics and cosmological parameters is in a slightly different direction than for the tSZ auto-spectrum, and so the combination of these two signals (e.g., in a joint analysis of the auto-and cross-spectra) can break these degeneracies and separately constrain astrophysics and cosmology.The tSZ -CMB lensing crosscorrelation signal has been detected only once before, in Ref. [1] (hereafter HS14), using the nominal-mission Planck data.In this work, we measure the signal with the final PR4 (NPIPE) release of the Planck data [30].
To measure the signal, we construct all-sky Compton-y maps (tSZ maps) from the single-frequency NPIPE maps, using a needlet internal linear combination (NILC) method [31].We describe our NILC pipeline and characterize our y-maps in detail in a companion paper [2] (hereafter "Paper I").The needlet ILC method allows us to build a linear combination of frequency maps with weights localized in both pixel and harmonic space.We cross-correlate our ymaps with various publicly available κ maps constructed from the Planck data, each with complementary systematics, including maps built from both the 2018 (PR3) data [7] and the PR4 (NPIPE) data [32].In particular, we use: (1) a PR3 κ map reconstructed with tSZ-deprojected temperature maps, in order to avoid spurious ⟨yyy⟩ three-point correlations in the measured signal; (2) a PR3 κ map with the signal at the location of tSZ clusters restored (it is standard to mask the brightest tSZ sources in the sky in the κ analysis mask), to avoid a systematic low bias to our measurement; and (3) the standard minimum-variance NPIPE κ map [32], which has lower noise than that of the 2018 κ maps, but no mitigation of any of the previously mentioned systematics.Our measurements with all the different κ maps are consistent.
Any measurement of an all-sky mm-wave signal necessarily includes some contamination from the cosmic infrared background (CIB), which is an unresolved, diffuse (at the resolution of CMB experiments) emission sourced by dusty star-forming galaxies at high redshift.As the redshifts that are most efficient for CMB lensing overlap significantly with the redshifts at which the CIB is sourced, the CIB has a very high correlation coefficient with CMB lensing, reaching 80-90% [33,34] -much larger than that of the tSZ signal.Thus, the bias induced by the correlation of residual CIB in a tSZ map with the CMB lensing convergence is a significant systematic in any measurement of the tSZ -CMB lensing signal.We avoid this bias by deprojecting the CIB in our NILC tSZ map, i.e., requiring that the NILC weights have zero response to an assumed CIB spectral energy distribution (SED).However, such a method is sensitive to the modelling of the CIB SED, which encodes its frequency dependence; additionally, it assumes that the CIB can be described as (e.g.) a perfect modified blackbody with no inter-frequency decorrelation.This is not a fully valid assumption, as the CIB does in fact display up to ≈ 10% decorrelation between frequency channels (or higher, if the channels are widely separated), because different frequency channels are sensitive to the emission at different redshift ranges that has been redshifted into the appropriate frequency (with high-frequency channels more sensitive to low-z emission that has experienced less redshift) [35][36][37][38].We avoid these systematics by using a moment-based approach, as suggested in [39].This approach desensitizes our measurement to CIB SED modelling systematics, as we verify using detailed simulations of the microwave sky that include both Galactic foregrounds from PySM3 [40] and extragalactic signals from Websky [41].Using these simulations, we show that we can recover an unbiased measurement of the tSZ -CMB lensing cross-correlation signal after appropriate moment-based CIB deprojection.
We constrain the astrophysics of the halos sourcing the tSZ -CMB lensing cross-correlation by comparing our C yκ ℓ measurements to a theoretical model, in particular by fitting a scale-independent amplitude parameter A to the data, where A = 1 corresponds to the prediction of our fiducial model, based on the pressure profile of Ref. [42].Our tightest constraint (from the NPIPE κ cross-correlation) is A = 0.82 ± 0.21; i.e., our measurement is consistent at ≈ 1σ with the fiducial model.For the other κ maps, our cross-correlation measurements remain broadly consistent, although the measurement with the tSZ-deprojected κ of A = 0.56 ± 0.24 is roughly 2σ lower than the fiducial model.These results are in 1 − 2σ agreement with the measurement of HS14, who found A = 1.10 ± 0.22 for the same fiducial model.Taken at face value, our measurement has similar error bars to those of HS14, which is surprising given the improved data quality in our work, but we emphasize that our signal-to-noise is significantly reduced by our use of the moment deprojection method to clean the CIB.Without this penalty, our signal-to-noise would surpass that of HS14 significantly, but we find that such methods are needed in order to obtain a robust measurement of the signal.We note that most previous such measurements, including, e.g., that of HS14 and [43][44][45], were not validated on detailed, non-Gaussian sky simulations containing tSZ, CIB, lensing, and other fields with realistic correlations, whereas ours is. 1 We thus emphasize that our C yκ ℓ measurement is much more robust to CIB contamination than that of HS14, which employed a simpler CIB-cleaning method that required an (unphysical) assumption that the intrinsic CIB-y cross-correlation vanishes.While our final astrophysical constraints are similar to those of HS14, their robustness and stability are validated with significantly more powerful and thorough methods.
This paper is structured as follows.In Section II we describe the theory of the tSZ signal and CMB lensing signal, and the halo model that we use to model their cross-correlation.In Section III we describe the datasets we use to detect the signal.In Section IV we discuss the foreground mitigation techniques for the ISW and CIB biases, including describing the components we deproject from our NILC y-maps and discussing the frequency coverage in the y-map construction.In Section V we describe our measurement of the tSZ -CMB lensing cross-correlation signal C yκ ℓ and present the results.Section V B in particular explicitly illustrates our use of different combinations of y-maps released in Paper I in order to remove CIB residuals in our final measurement.In Section VI we present the results of our pipeline validation on simulations, demonstrating that our measurement is unbiased to CIB systematics, and comment on the consistency with the previous measurement of this signal.In Section VII we analyze the signal, quantifying the significance of our detection and constraining separately the ICM physics that controls the pressure-mass relation, and cosmological parameters.We discuss our results and conclude in Section VIII.
We assume the cosmology of [48] throughout: {H 0 = 67.32km/s/Mpc; σ 8 = 0.812; n s = 0.96605; Ω b h2 = 0.022383; Ω cdm h 2 = 0.12011}, where H 0 ≡ 100h is the Hubble parameter today; σ 8 is the amplitude of linear fluctuations on a scale of 8h −1 Mpc today (i.e., the linear matter power spectrum integrated over all scales, smoothed with a top-hat window function of 8h −1 Mpc); n s is the spectral index of the primordial fluctuation power spectrum; Ω b h 2 is the physical density of baryons today; and Ω cdm h 2 is the physical density of dark matter today.

II. THEORY
We model the tSZ signal and the CMB lensing signal using the halo model.In this section, we present the theory and model we use to calculate their cross-power spectrum C yκ ℓ .We perform all theory calculations throughout with class sz 2 [49][50][51], which is an extension of the cosmological Boltzmann solver class3 [52].Sample code showing how to use class sz to calculate the observables is available at https://github.com/fmccarthy/tSZ_kappa_halomodel_classsz.

A. The tSZ anisotropies
The tSZ temperature anisotropy at sky position n observed at frequency ν is given by where g ν is the spectral function of the tSZ effect [12]: Here h is Planck's constant; k B is Boltzmann's constant; T CMB = 2.726 K is the mean temperature of the CMB [53][54][55]; and y(n) is the dimensionless (and frequency-independent) Compton y-parameter that quantifies the integral of the electron pressure along the line of sight (LOS): Here σ T is the Thomson scattering cross-section; m e is the electron mass; c is the speed of light; a(χ) is the scale factor at comoving distance χ; and P e (χ, n) is the electron pressure at (χ, n).The integral is done over comoving distance χ from today to the beginning of the epoch of reionization at redshift z ∼ 7 − 10.We model the distribution of P e (χ, n) within the halo model; the details are given below in Section II D.

B. The CMB lensing potential
The CMB lensing convergence at sky position n is given by a line-of-sight integral over δ(χ, n), which denotes the matter overdensity at (χ, n): where the CMB lensing convergence kernel W κ (χ) is given by with χ S the comoving distance to the surface of last scattering, at which the CMB was released (corresponding to z ≃ 1100), and Ω m the matter density parameter.For a review of CMB lensing theory, see [56].

C. The halo model
The halo model (see, e.g., [57] for a review) has long been used to model the distribution of the large-scale structure of our Universe [58,59].Within the halo model, the continuous overdensity field is replaced by a discrete distribution of "halos", which cluster according to the underlying density field but with a bias that is scale-independent on large scales.These halos are modeled as having a spherically symmetric density profile ρ(r), where r is the distance from the center of the halo.In the simplest picture, all halo properties and observables are functions only of the halo mass M and the redshift z.The halos are distributed according to the halo mass function (HMF) dN dM (M, z), which gives the comoving differential number density of halos of mass M at redshift z.Note that the halo model makes several unrealistic simplifying assumptions, such as assuming that all halos are spherical and that all properties of the halo and its constituent galaxy (or galaxies) are functions only of the host halo mass and redshift, neglecting real effects due to scatter and environment.
Within the halo model, all correlation functions are split into various N -halo terms, where N is the number of halos whose profiles appear in that term of the correlation function.For the two-point correlation function, this means that there is a 2-halo term, which describes the correlations between two different halos (due to clustering), and a 1-halo term, which describes the correlations within a single halo.The halo mass function is an important quantity for calculating these correlations, as is the halo bias, which is important on scales where halo clustering is relevant (i.e., the two-halo term).The halo bias quantifies the bias of the halo distribution with respect to the underlying dark matter distribution: where δ is the underlying (continuous) dark matter overdensity; δ h is the halo overdensity; and b h (M, z) is the halo bias.It is assumed that halos are spherically symmetric and that the dark matter within them follows a continuous density profile; in our case we take this to be a Navarro-Frenk-White (NFW) profile [60] ρ(r, M, z) = ρ s r r s where ρ s is a characteristic density and r s is a characteristic NFW scale radius.In our halo model calculations, we use the halo mass function of [61] and the halo bias of [62], as implemented in class sz (for implementation details, see the appendices of [51]).
In general, to define quantities that depend on a halo mass, one requires a definition of halo mass.Several common definitions exist; for example, the spherical overdensity (SO) mass, which is the mass M ∆c,m within the radius R ∆c,m within which the mean density of the halo is ∆ times either the critical density ρ cr (z) (in the case of the subscript c) or the mean background matter density ρ m (z) (in the case of the subscript m).The halo mass function of [61] was defined for M ∆m , at several values of ∆.As the pressure profiles we will use to describe the tSZ signal (see below in Section II D) are defined for M 200c , we use a halo mass function that is defined for this mass, by using M ∆m with a z-dependent ∆ = 200/Ω m (z), which is appropriate as ρ m (z) = Ω m (z)ρ cr (z) (as implemented directly in class sz).
In all our integrals over mass, we use the integration limits 10 10 h −1 M ⊙ < M < 10 15.5 h −1 M ⊙ .In our integrals over z, we use the integration limits 0.005 < z < 10.Note that this neglects the contribution from halos below the lower mass integration limit, for which the HMF is not calibrated.We account for the contribution from these halos using counter-terms [63] (see Appendix B2 of [51] for further details).Additionally, we extrapolate the HMF to higher z than where it was calibrated; following the recommendation of Ref. [61], we deal with this regime by removing the z-evolution of the halo multiplicity function f (σ) and replacing it with its values at z = 2.5 for all z > 2.5.
Finally, we cut all of our radial profiles at r = 2R 200m .This is a choice, and in practice observables depend on this cut-off as the NFW profile does not converge as r → ∞ (although the profiles we use for the tSZ signal do converge).We choose this value as we find that it allows us to reproduce most accurately a measurement of the C yκ ℓ signal from the hydrodynamical simulations from which the tSZ profiles were extracted [64] (see Fig. 16 in Appendix A).Choosing a cut-off for the NFW profile that is not equal to the radius we use to define the mass of our halos leads to some subtleties in the matter power spectrum calculation; we discuss this in Appendix A.

D. tSZ power spectrum within the halo model
In [65,66], a halo model prescription for the calculation of the auto-power spectrum of the tSZ effect was presented.A pressure-mass relation P e (z, M ) is assumed, and (as halos are assumed to be spherically symmetric), a threedimensional pressure profile P 3D (z, M, r) can be written down that describes the pressure at distance r from the center of the halo.The power spectrum, which involves the Fourier transform of this pressure profile, can then be written as where C yy,2h ℓ is the two-halo term that depends on the distribution and clustering of the halos, and C yy,1h ℓ is the one-halo term that probes the distribution of the pressure within each halo.As C yy ℓ is dominated (at ℓ > ∼ 300) by the one-halo term, Ref. [66] neglected the two-halo term; however, for an in-depth discussion and presentation of this (and the one-halo term), see [67] (see also [65]).The final expressions (in the flat-sky and Limber approximations [68]) are where P lin (k, z) is the linear matter power spectrum and both expressions depend on ỹℓ (M, z), the 2-dimensional Limber projection of the 3-dimensional Fourier transform of the pressure profile, which is given by Parameter A0 αm αz P0  [42] for the three-dimensional pressure profile of Equation (13).The parameters A0, αm, and αz are defined in Equation (15).
Here r y is a characteristic scale radius of the three-dimensional pressure profile, for which we use r y = R 200c ; ℓ y = a(χ)χ ry is the characteristic multipole moment of r y ; and y 3D (x, M, z) is the 3-dimensional Compton-y profile of a halo of mass M at redshift z for which a model is required (note the integral in Equation ( 11) is in terms of the dimensionless parameter x = r ry while the model for y 3D (x, M, z) is usually expressed in terms of the dimensionful distance r from the halo center).The Compton-y profile is directly related to the pressure profile P 3D by As our fiducial model for P 3D , we use the generalized NFW-based fitting functions of [42].These are given by where P ∆ , the self-similar amplitude for pressure, provides the physical dimensions of pressure: Here G is Newton's constant and f b ≡ Ω b Ωm is the baryon fraction.Recall that the mean density of the halo within R ∆c is ∆ρ cr (z) (we are working with ∆ = 200).In Equation (13), the amplitude parameter P 0 , the core-radius parameter x c , and the power-law index β are parameters that are fit to hydrodynamical simulations [69], while α and γ are fixed at α = 1 and γ = −0.3[42].P 0 , x c , and β are allowed to have mass and redshift dependence and each are of the form with A 0 , α m , and α z fit separately for each parameter.The best-fit values (to the simulations of [69]) are listed in Table 1 of [42]; for completeness, we repeat them in Table I.This model is implemented directly in class sz.

E. CMB lensing power spectrum within the halo model
We can also model the matter power spectrum within the halo model, which can then be weighted appropriately with the CMB lensing kernel and integrated against redshift to calculate the CMB lensing power spectrum.Again, there are both one-halo and two-halo contributions.The final expressions are similar in form to Equations ( 9) and ( 10), but with ỹℓ replaced by the CMB-lensing-kernel-weighted 2-dimensional mass profile κℓ : with where r s is the characteristic scale radius for κ, which is indeed the NFW scale radius of Equation ( 7), and ℓ s is the characteristic multipole ℓ s = a(χ)χ rs .
Care must be taken when the mass integral is evaluated in the calculation of the matter power spectrum (and thus the CMB lensing power spectrum).In particular, there is a significant contribution to the signal from low-mass halos below the range for which the halo mass function has been calibrated; also, the result depends critically on any low mass cut-off imposed in the integral.This is in contrast to the case of the tSZ power spectrum, where any contribution from low-mass halos is severely downweighted by their low gas pressure amplitude.
We account for the contribution from lower-mass halos as follows.We require that the halo mass function and bias obey the consistency condition This condition ensures that the dark matter is unbiased with respect to itself.Following Ref. [63], we write this as where M cut is the low-mass cut-off of the integral we use in practice when calculating halo model quantities.This allows us to define (and explicitly calculate) a counter-term explicitly as the integral over all the low-mass halos: We can then replace , where ϵ is an infinitesimal mass.This allows us to account for the contribution from the lower-mass halos in the two-halo power spectrum while preserving the consistency conditions and also, importantly, not needing to model any properties of the low-mass halos below the cut-off.Note that this issue is not as serious in the one-halo term, which is dominated by the massive halos.

F. tSZ -CMB lensing power spectrum within the halo model
We can use the halo model formalism to write down an expression for the tSZ -CMB lensing cross-correlation C yκ ℓ .As in HS14, we use one factor of ỹℓ (M, z) and one factor of κℓ (M, z): Again, we must replace the contribution from the low-mass halos to κℓ in the two-halo term with an appropriate counter-term as described above.However, while the counter-terms for y are also well-defined according to the formalism of [63], we choose not to add them.The counter-terms account for the contribution of the halos below the low-mass limit, which we take to be M = 10 10 h −1 M ⊙ ; this is already extrapolated far below the mass range in which the y profiles were fit, and indeed we expect very little contribution from such low-mass objects.Thus, the counter-term addition in C yκ ℓ probes no gas in such low-mass halos, and only the clustering of the low-z lensing halos with the gas in the larger-mass halos through the 2-halo term.

Comparison to simulations
In Appendix A, we explicitly compare our halo model prediction to a measurement of the tSZ -CMB lensing cross-power spectrum in the hydrodynamical simulations from which the pressure profile model was calibrated [64].While there is good agreement at high ℓ, at ℓ < ∼ 2000 (which is the regime where we measure the signal), the halo model under-predicts the total signal.
As the simulations were performed with different cosmological parameters than those we use for our theory model, we cannot use the measurement from simulations as a template with which to model the signal.However, we can use it to calibrate a ratio has a significant contribution even from z < 0.1, while C yκ ℓ is more evenly distributed in the redshift range 0 < z < 1, and also receives non-negligible contributions (≈ 25% of the total signal) from z > 1.
simulation), which we can assume to be cosmology-independent at first order.We can then use this ratio to rescale our halo model calculation to account for this mismatch, which is likely due to unbound gas blown out of halos into the IGM.
While we choose to retain the un-rescaled halo model theory prediction in all plots and when we quote our fiducial constraints, we will also quote constraints on an amplitude of this rescaled halo model prediction when we constrain the model in Section VII B. Of course, the best-fit amplitude for this model will necessarily be lower than that of the unrescaled model, as the rescaled model predicts a higher signal by ∼ 25% over our multipole range of interest.

G. Halo mass and redshift dependence of the signal
To source the tSZ effect, the electrons involved in the Compton scattering process must have high temperature.They inherit this high energy from the gravitational collapse of massive structures, and thus the highest-temperature electrons are in the most massive clusters (T e ∼ M 2/3 in a self-similar ICM model [70]).These are primarily located at low z, as the formation of such large structures requires nearly a Hubble time.Thus, the tSZ effect is primarily a low-z probe, with most (> 90% of the contribution to C yy ℓ in the range ℓ < 2000) of the signal sourced at z < 1 (and > 95% of the signal for ℓ < 1000).The cross-correlation with CMB lensing, while also being a low-z probe (with closer to ∼ 70% of the contribution to C yκ ℓ being sourced at z < 1 for ℓ < 2000), is weighted more towards high redshifts than C yy ℓ , as the objects with high efficiency for CMB lensing are located at higher z due to the lensing kernel.We illustrate this in Figure 1, where we plot C yy ℓ (z > z min ) and C yκ ℓ (z > z min ), i.e., the contributions to the respective signals from z > z min .
In Figure 2 we show the dependence of the signals on the highest mass included in the halo model integral.These plots illustrate the significantly lower mass range probed by C yκ ℓ compared to C yy ℓ .Roughly 50% of the C yκ ℓ signal comes from halos with M > 10 14 h −1 M ⊙ ; the equivalent number for the auto-power spectrum is ≈ 90%.This is due to the combination of two reasons: first, the CMB lensing signal is sourced at higher z, where there are fewer high-mass halos, and so is preferentially probing a redshift regime with lower-mass halos than the tSZ signal.Second, at every z, the highest mass halos that source the largest tSZ signal are very rare, and contribute fractionally less of the total mass at a given redshift than the lower-mass halos.Thus most of the CMB lensing signal (which is weighted by overall mass) is sourced by the lower-mass halos that host most of the mass.

III. DATA A. Intensity data
We analyze Compton-y maps constructed by applying a NILC pipeline to the single-frequency maps from the Planck NPIPE data release (PR4) [30].Our baseline data set includes all maps from 30 to 545 GHz.We compare various choices of component deprojection in these maps in order to minimize residual CIB contamination.We discuss the y-maps and our needlet ILC pipeline in detail in Paper I.

B. CMB lensing data
Maps of the CMB lensing potential ϕ (or convergence κ) can be reconstructed from CMB temperature and polarization anisotropy maps by taking advantage of the well-theoretically-understood non-Gaussianities and statistical anisotropies induced in an underlying (assumed Gaussian and statistically isotropic) primary CMB map by the intervening gravitational potential.Optical quadratic estimators (QEs) [9,10] have been derived for such an operation.The lowest-noise estimation of κ from Planck data used the globally optimal estimator ("generalized minimum variance (GMV) QE") of Ref. [10] applied to the PR4 Planck NPIPE maps [32]; previous measurements used the minimum-variance (MV) QE of Ref. [9] (the Hu-Okamoto QE).
This lowest-noise NPIPE measurement is not necessarily the optimal one for our analysis, as the lensing analysis as part of the 2018 Planck data release (PR3) [7] provides κ estimations with various analysis choices that are more suitable for a robust measurement of C yκ ℓ .We consider the following κ maps for our C yκ ℓ measurement: 1.The κ reconstruction from the 2018 Planck release which applied the Hu-Okamoto QE to tSZ-deprojected CMB maps to reconstruct a κ map free of any bias due to residual tSZ signal 4 .As any C yκ ℓ measurement using a QE is in fact a ⟨T T T ⟩ three-point function (with one T leg coming from the y estimate and two coming from the quadratic-in-T κ estimate), the deprojection of y in the temperature map used to estimate κ avoids any potential ⟨yyy⟩ contribution from the bispectrum of the Compton-y field in the C yκ ℓ measurement (or mixed bispectra of the form ⟨yyT CIB ⟩ or ⟨yyT radio ⟩). 5.The κ reconstruction from the 2018 Planck release without the tSZ-selected galaxy cluster mask applied. 6It is standard to provide κ maps with analysis masks that mask out these clusters, as the lensing reconstruction can become biased at these locations.Measurement of the C yκ ℓ signal using such a mask, which is highly correlated with the tSZ signal, would bias the signal low (although most of the tSZ -CMB lensing cross-correlation comes from much lower halo masses than those found in such a catalog, as shown in Figure 2).
3. The lowest-noise estimation of κ from the NPIPE Planck data. 7e measure C yκ ℓ separately using all three of the above κ reconstructions from Planck (NPIPE; tSZ-deprojected; tSZ-clusters-unmasked) and compare the signal and the signal-to-noise ratio (SNR) in every case.We find broadly consistent signals with all maps, with the NPIPE measurement having the highest SNR.However, we find below that the tSZ-deprojected measurement is slightly lower, by ≈ 1σ. 8

IV. FOREGROUND MITIGATION
In our Paper I we present a set of Compton-y maps with various deprojection choices to remove other components from the final y-map.There are, in particular, (at least) two components that must be removed in order to make an unbiased detection of the ⟨yκ⟩ signal: the integrated Sachs-Wolfe (ISW) signal [71] signal and the cosmic infrared background (CIB) signal.These two components are present in the mm-wave intensity maps that we use to construct our y-maps and are correlated with large-scale structure (and thus with κ).Thus, they can induce biases in our C yκ ℓ measurement if they are not sufficiently removed from the y-map prior to the cross-correlation measurement.We describe these signals and the associated biases below, and our techniques to mitigate them.

A. Mitigating the ISW biases
The ISW signal [71] is sourced when primary CMB photons travel through a time-varying gravitational potential.In a constant potential, a photon will not gain or lose energy as it enters and leaves the potential.However, if the potential evolves while the photon passes through it, it can result in the photon requiring less or more energy to leave the potential than it gained entering the potential.This effect generates anisotropies in the CMB during the periods in the history of the Universe when potentials were time-varying.During matter domination, potentials are constant.However, during the more recent dark energy domination, potentials have been decaying, and so there is an ISW signal imprinted on the primary CMB, which is most significant on large angular scales as these correspond to the projection of late-time linear modes.Note that the ISW effect preserves the blackbody spectrum of the primary CMB photons.
The ISW-κ cross-correlation on large scales is dominant (or comparable) to the tSZ-κ signal up to ℓ ∼ 100−200 (see, e.g., Fig. 1 of [72]) at typical tSZ-sensitive frequencies (e.g., ≈ 100 GHz), so it is important to deproject the primary CMB on these scales when building the y-map for use in our tSZ -CMB lensing cross-correlation measurement.Fortunately, such a deprojection is simple due to the well-known blackbody SED of the primary CMB anisotropies.Ideally, we would deproject the CMB everywhere in our NILC y-map, but we will find that we do not have enough frequency coverage to deproject the CMB and the CIB along with its first moments on small scales.As such, we only deproject the CMB in the first five needlet scales in our y-map; we will refer to this as CMB 5 -deprojection.

CIB mitigation: constrained ILC
The CIB is highly correlated with the CMB lensing signal (e.g., [34]), and it is imperative to remove this signal from the y-map before measuring C yκ ℓ .Thus, we need to deproject the CIB using an estimate of its frequency dependence.However, unlike the tSZ and the CMB signals, which display no frequency decorrelation and which have SEDs that are well-understood from first principles and can thus be calculated theoretically, the CIB is not described perfectly by one SED.It is sourced by the line-of-sight-integrated thermal emission of different objects; in particular, different frequency channels are sensitive to slightly different objects, as source emission at different redshifts will be redshifted into different frequency bands.This leads to frequency decorrelation between the CIB channels.However, as the correlation coefficients are ≳ 90% at the frequencies of interest [35][36][37][38], it is still possible to clean the CIB using multifrequency measurements.Its SED does not need to be known or modeled in order for the CIB to be (partially) cleaned in a standard, unconstrained ILC.However, if we wish to explicitly deproject it, we must model its SED.We model the CIB SED as a modified blackbody: where B ν (T ) is the Planck function The SED Θ CIB ν depends on two parameters: the spectral index β and the dust temperature T eff CIB .In Paper I, we fit this SED to estimates of the CIB monopole, and find best-fit parameters of β = 1.77,T eff CIB = 10.14 K.However, there is significant uncertainty on these parameters; additionally, it is not a perfect approximation to the CIB emission to describe it as an exact modified blackbody.Thus, in addition to deprojecting the CIB, we also deproject the first moment(s) of the CIB with respect to these parameters.We describe this approach in detail in Paper I; in short, it involves additionally deprojecting components whose SEDs are exactly defined by the derivative with respect to the appropriate parameter of the modified blackbody SED.We refer to these components as δβ (for the β-moment) and δT eff CIB (for the T eff CIB moment).

CIB mitigation: decreased frequency coverage
We note that, while using more frequency channels in the NILC allows for better characterization of the spatial structures of foregrounds and lower variance in the final map, it also requires better characterization of the CIB SED in order to deproject it.Thus, we explore whether using fewer frequencies than standard can lead to an improved CIB deprojection, in particular by dropping the 545 GHz information from the NILC.We additionally always leave out the 857 GHz channel, both because it is not calibrated as precisely as the other channels, and because it requires extrapolation of the CIB modified blackbody SED to higher frequencies than we have been using to constrain it.

C. Analysis masks
After performing the NILC analysis to build the y-map, non-negligible residual foreground power remains in the dustiest areas of the sky (e.g., close to the Galactic plane), especially in the δβ-deprojected tSZ map.This can add significant variance to our measurement.To avoid this, we mask the dustiest areas of the sky, which we define by thresholding the highest (1 − f sky ) × 100% of pixels in the Planck 857 GHz single-frequency map, where f sky is the fraction of sky on which we wish to make our measurement.We use f sky = 0.8 as our fiducial choice; we explore effects of varying the masked sky area in Appendix B.
Additionally, the reconstructed CMB lensing maps from Planck are provided with analysis masks applied.We multiply our f sky -thresholded mask by this analysis mask.Finally, we also multiply this mask with the union of the Planck HFI and LFI point source masks, which were used in the preprocessing of the single-frequency maps before performing the NILC [2](as a result, the resulting maps contain no information at the locations of these point sources).In principle, this slightly suppresses the Compton-y signal, due to the tSZ-source correlation, but this effect is negligible at Planck signal-to-noise [73,74].We apodize the resulting mask with an apodization scale of 10 FIG.3: The various masks used in our analysis; the top two rows include the unapodized masks that we multiply to make the "combined" masks in the bottom row.Every "combined" mask is made of the product of the f sky = 80% threshold mask with the point source mask and one of the lensing masks; the final product has a 10 ′ apodization applied as described in the text.
arcminutes, using the "C1" apodization procedure of NaMaster [75]9 .Most of the sky area that is covered by the f sky = 80% threshold mask is indeed covered by the κ analysis masks, such that our overall masks are essentially limited by the size of these κ masks.Note that, as is standard, the κ analysis masks remove regions containing massive tSZ-selected clusters [32], and so the use of these masks has the risk of biasing our tSZ -CMB lensing cross-correlation signal low.To avoid this, it is desirable to use a κ map which does not mask these clusters, and so we perform a measurement with a Planck 2018 κ reconstruction in which these clusters are not masked.However, as discussed in Section III B, we also perform a measurement using the Planck 2018 κ reconstruction performed on tSZ-deprojected temperature maps, which avoids a potential ⟨yyy⟩ intrinsic bispectrum bias in our measurement, as well as on the κ reconstruction from the NPIPE maps, which has the lowest noise.These final two κ reconstructions are provided with analysis masks that remove tSZ clusters, meaning that there is potential for the signal to be biased low.However, we find that the measurements with the three κ reconstructions are all statistically consistent, indicating that any low bias due to the masking of tSZ clusters is negligible.
The relevant masks are shown in Figure 3.The combination of the f sky = 0.8 threshold mask and the lensing analysis mask does not remove a large sky area, with the remaining sky fraction decreasing from 80% to 61.57% for the tSZ-deprojected κ analysis mask; 62.74% for the tSZ-clusters-restored mask; and 62.24% for the NPIPE mask.

V. THERMAL SZ -CMB LENSING CROSS-CORRELATION MEASUREMENT
In this section, we present our measurement of the C yκ ℓ signal.In particular, we describe in Section V A our estimate of the cross-correlation data points and associated covariance.In Section V B we discuss systematics induced by the CIB, in particular by comparing the data points as measured using different assumptions for the parameters of the CIB SED in the deprojection scheme used in the y-map construction (this analysis provides an explicit example of how to use the CIB-cleaned y-maps of Paper I to remove CIB residual bias).In Section V C we show our final results for the different κ maps we use.
Following from the standard expression for the covariance in the measurement of a cross-power spectrum of Gaussian fields, (where δ ℓℓ ′ is the Kronecker delta function), the error bars on the power spectrum Ĉyκ ℓ measured on a fraction of the sky f sky are given by where C AB ℓ are the true underlying power spectra and N AB ℓ are the noise power spectra (including instrumental noise and residual foregrounds).In our case, we use measured auto-power spectra Ĉyy ℓ term, and we assume N yκ ℓ = 0, i.e., that all sources of noise bias and foreground power on this measurement are zero (this assumes that any Galactic foregrounds that remain in the tSZ map do not appear in the κ reconstruction and that the Planck instrument noise has vanishing skewness).These assumptions are robust due to the fact that the noise power spectra of the y-map and CMB lensing map highly dominate the error bars on C yκ ℓ .Again, we use NaMaster to decouple the mask in the final calculation of our error bars (using the function gaussian covariance()); note that this requires us to interpolate our measured values of Ĉyy ℓ and Ĉκκ ℓ , as this function requires a theoretical estimate of C XY ℓ at all multipoles.Note that the decoupling of the mask with Namaster induces off-diagonal (ℓ ̸ = ℓ ′ ) terms in the covariance matrix. 10e measure C yκ ℓ using y-maps with various deprojection choices and frequency coverage.In Figure 4, we plot the data points for various deprojection choices in the NILC.We plot separately the points for the cases when we have different frequency coverage in our NILC (the "standard frequency coverage" case corresponds to when we include all frequencies except for 857 GHz).In general, we see an extremely strong signal in the case when we do not deproject the CIB (labeled "CMB deprojection") in the figures.When we deproject the CIB, the signal is significantly lowered, indicating that the non-CIB-deprojected measurement is biased.Deprojecting further moments results in points with larger error bars but much less scatter (with respect to the chosen parameters of the CIB SED), eventually converging on a stable, robust measurement.We discuss the different deprojections below.

B. Systematics due to residual CIB contamination
Our first step to remove the CIB contribution to our measurement is a simple deprojection in the NILC of a component with the appropriate frequency scaling for the CIB.We discuss this deprojection in detail in Paper I, in The tSZ -CMB lensing cross-power spectrum, C yκ ℓ , as measured for different frequency coverage (left panel including all frequencies from 30 -545 GHz and right panel with 545 GHz removed) and different deprojection combinations (as labeled in the legend).Note that, for the case when we exclude 545 GHz, we do not have enough frequency channels to simultaneously deproject the CIB and both δβCIB and δT eff CIB on the smallest scales; it would be possible to do so on large scales, but we do not consider such a measurement.The solid black curve shows our fiducial theoretical model, with the one-halo (two-halo) term shown in dashed (dash-dotted).Note that, in the legend, CMB 5 refers to the CMB deprojected on the largest five needlet scales.Also note that some points are slightly offset on the x axis to allow for better visualization.Both plots were made with the SZ-deprojected reconstruction of the κ map.
which we directly fit a modified blackbody to the CIB monopole calculated in Ref. [36], and find best-fit values for the spectral index β and the dust temperature T eff CIB of (T eff CIB , β) = (10.14K, 1.77).These are the values that we use for our standard CIB deprojection, as presented in Figure 4.
However, there is significant uncertainty on the CIB SED, including significant parameter degeneracy between β and T eff CIB , as well as uncertainties on the CIB monopole estimate to which the SED is fit, along with the fact that the CIB is not perfectly described by a modifed blackbody SED (and that the CIB monopole is not necessarily the appropriate observable to consider in this case, as it is sourced at different redshifts and thus has different frequency dependence to the CIB anisotropies).Thus we also make our measurement for y-maps on which various values of (β,T eff CIB ), which are drawn from the posterior for these parameters [2] (with values as indicated in the figure legend), have been used in the deprojection.The resulting changes in the data points are shown in Figure 5, where we note that we have included the primary CMB deprojection in the first five needlet scales in all cases.It is clear that there is a significant systematic uncertainty associated with the exact CIB SED used in the deprojection, as the spread of these data points is much larger than the statistical error bars.However, we note that on small scales (ℓ > ∼ 1000), when we drop 545 GHz from the NILC, this systematic is significantly reduced (see the right panel of Figure 5).Thus, in order to reduce the systematic uncertainty associated with the CIB SED deprojection, we proceed to measure the tSZ -CMB lensing cross-correlation using y-maps from which the first moment of the CIB with respect to β has been deprojected.We show the resulting measurements in Figure 6.We note that the scatter in the "standard-frequency" case is still significantly larger than the size of the statistical error bars for some multipole bins.However, in the "no-545 GHz" case, this behavior is no longer seen.Thus, for the "no-545 GHz" case, we use this deprojection configuration to make our measurement, as we have achieved stability.For the final data points that we use in our analysis, we use the deprojection done with the best-fit SED (β = 1.77,T eff CIB = 10.14 K.) However, for the standard-frequency case it still remains to find a stable configuration for which we can proceed.Thus, we additionally deproject the first moment of the CIB SED with respect to T eff CIB .The resulting points are shown in Figure 7.These points are stable, i.e., the variation with the CIB SED used in the deprojection is much less than the statistical error bars for all multipole bins.For the final analysis we use the best-fit CIB SED from Paper I  FIG.5: The tSZ -CMB lensing cross-power spectrum, C yκ ℓ , measured on maps with the CIB deprojected, for various choices of the CIB SED parameters.The SED parameters used here are obtained by taking ten samples from the posterior of (β, T eff CIB ) which we calculate in Paper I. On the left, we show the standard-frequency-coverage case (where we include 545 GHz but no 857 GHz), and on the right the case where we remove 545 GHz from the NILC.In all cases, the primary CMB is deprojected in the first five needlet scales of the NILC.The points are systematically offset along the x axis for better visibility of overlapping points.Both plots were made with the SZ-deprojected reconstruction of the κ map.
to perform the deprojections.

Summary of the CIB removal and deprojection choices
To summarize, our method for removing the CIB residual bias is as follows: 1. Make a standard, non-deprojected (or CMB-deprojected) NILC y-map with a given frequency coverage and measure the observable (C yκ ℓ ). 2. Make a CIB-deprojected (or CMB+CIB-deprojected) NILC map for various choices of the CIB SED and measure the observable.If the measurement is stable to the choice of CIB SED, use this measurement.If it is not, move on to Step 3.
3. Make a CIB+δβ-deprojected (or CMB+CIB+δβ-deprojected) NILC y-map for various choices of the CIB SED and measure the observable.If the measurement is stable to the choice of CIB SED, use this measurement.If it is not, move on to Step 4.
4. Make a CIB+δβ + δT eff CIB -deprojected (or CMB+CIB+δβ + δT eff CIB -deprojected) NILC y-map for various choices of the CIB SED and measure the observable.If the measurement is stable to the choice of CIB SED, use this measurement.If not, continue deprojecting higher moments of the CIB SED until stability is found.
For the standard-frequency-coverage case, we obtain stability in our measurement after Step 4 (see Figure 7).For the no-545-GHz case, we obtain stability after Step 3 (see Figure 6, right-hand side).

C. Results for the different κ maps
We also explore the changes in our measurement when we use different κ maps.Using the CMB 5 +CIB+δβdeprojected no-545 GHz y-map and the CMB 5 +CIB+δβ+δT eff CIB -deprojected standard-frequency y-map described above, we show in Figure 8 the signal as measured with different κ maps.These include the κ map released with an analysis mask that restores the signal at the location of tSZ clusters; the κ map reconstructed from tSZ-deprojected  FIG.6: The tSZ -CMB lensing cross-power spectrum, C yκ ℓ , measured on y-maps with the CIB and the first moment of the CIB SED with respect to β deprojected, for various choices of the CIB SED parameters.The SED parameters are the same as those in Figure 5. On the left, we show the standard-frequency-coverage case (where we include 545 GHz but no 857 GHz), and on the right the case where we remove 545 GHz from the NILC.In all cases, the primary CMB is deprojected in the first five needlet scales of the NILC.The points are systematically offset along the x axis for better visibility of overlapping points.Both plots were made with the SZ-deprojected reconstruction of the κ map.
temperature maps; and the κ map reconstructed from the Planck PR4 (NPIPE) map.It is encouraging that our measurements with different κ reconstructions look broadly consistent, as they have different sensitivity to systematics (in particular, the tSZ-deprojected map should not suffer from a potential ⟨yyy⟩ bias from spurious tSZ signal in the κ map).
While it is encouraging that our measurements are broadly consistent, we continue to use all three κ maps to present our results, although for ease of presentation in Figures 4, 5, 6, and 7 we have chosen the tSZ-deprojected κ map for the plots.The corresponding results with the different κ maps are similar.

VI. PIPELINE VALIDATION A. Application on Websky simulations
In the previous section, we performed various null tests and stability tests on the data that inspire confidence in our results.The stability with respect to different CIB SED parameter choices encourages us that we have truly deprojected the CIB, and the stability of our data points with respect to different sky fractions (see Appendix B) reassures us that we are not seeing signal from residual Galactic foregrounds in the CMB lensing reconstruction correlating with those in the y-map.However, we wish to also validate our pipeline on simulations, to ensure that performing the δβ and δT eff CIB deprojections satisfactorily removes the CIB contamination and preserves our signal of interest.
As such, we build simulations of the microwave sky at the frequencies we use in our NILC, in particular {30, 44, 70, 100, 143, 217, 353, 545} GHz.We include simulated cosmological signals as well as the Galactic foregrounds, using the Websky [41] simulations for the cosmological signals and the Python Sky Model (PySM311 ) [40] for the Galactic foregrounds.From Websky, we include the lensed primary CMB (blackbody in its frequency dependence), the CIB, and the tSZ field. 12The Websky CIB model, which is based on the halo model of Refs.[35,76] FIG.7: The -lensing cross-power spectrum, C yκ ℓ , measured on y-maps with the CIB deprojected as well as its first moments with respect to both β and T eff CIB , for various choices of the CIB SED parameters.The SED parameters are the same as those in Figure 5.In all cases, the primary CMB is deprojected in the first five needlet scales of the NILC.The points are systematically offset along the x axis for better visibility of overlapping points.This plot was made with the SZ-deprojected reconstruction of the κ map.CIB deprojection applied to the standard frequency coverage (CMB 5 +CIB+δβ deprojection applied to the no-545-GHz frequency coverage).The measurements are broadly consistent for all three κ maps.Some points are slightly offset along the x axis for better visibility of overlapping points.
appropriate light-cone evolution and inter-frequency decorrelation at the frequencies at which it is simulated (ν > 93 GHz).At lower frequencies than this, we simply scale the lowest-frequency CIB map13 by the SED of Equation ( 25) (note that the CIB is extremely subdominant to other sky components at these frequencies).
The Galactic foregrounds that we include are those corresponding to the [d1,s1,a1,f1] model of PySM3.Thus we include components due to Galactic dust (d1); synchrotron radiation (s1); anomalous microwave emission (a1); Frequency (GHz) Beam FWHM (arcmin) Noise power spectrum amplitude (µK arcmin) Noise (µK  II: Characterizing features of the Planck experiment, in particular the beam FWHM and white noise levels for each frequency.We quote the noise both in µK arcmin and in µK 2 ; the conversion between the two is given by Equation ( 29).See Ref. [77] for the LFI characterization and Ref. [78] for the HFI characterization.Note that we have applied a Jy/sr-to-µK conversion factor to the values quoted for 545 GHz, as we simulate the maps in µK.
and free-free emission (f1).The dust is the dominant foreground at high frequencies, and the synchrotron at low frequencies.We refer the reader to the PySM3 documentation for details of these models, which are built from analyses of Planck, WMAP, and other data sets.We find that these simulations reproduce appropriately well (for the purposes of this pipeline validation, which we stress is not aiming to assess the level of Galactic foreground removal but instead the recovery of the true cosmological signal for our NILC deprojection choices) the measured power spectra of the Planck data both on the full sky and on a patch of sky masked by the smallest-area Galactic mask released by the Planck collaboration, which has a sky area of 20%.
As we use the true, simulated Websky κ map as our proxy for the CMB lensing map, it would be futile to build an overly-complicated Galactic foreground map, as our pipeline validation is not sensitive to the true source of systematics from Galactic foregrounds, i.e., the correlations between residual foregrounds in the κ reconstruction and the NILC y-map.A true pipeline validation would propagate the Galactic biases through the lensing reconstruction pipeline by performing lensing reconstruction on the simulated single-frequency maps (or, ideally, on an ILC combination or for various independent simulations of the noise).Such a thorough procedure is outside the scope of this work, as we do not perform lensing reconstruction ourselves but instead use publicly-available κ maps.As we find results that are stable to large changes in the sky area used in our measurements (see Appendix B), we do not expect these Galactic foreground biases to be a significant issue.However, for future measurements of this signal (e.g., with forthcoming high-resolution CMB data), such a simulation validation would be interesting and useful.
After adding the Galactic foregrounds to the cosmological signals, we convolve our maps with Gaussian beam window functions with full-width-at-half maxima (FWHMs) corresponding to those of Planck.These FWHM values are listed in Table II.Finally, we add isotropic Gaussian white noise corresponding to the Planck noise levels, which are also listed in Table II.In particular, we simulate fields with white-noise power spectra corresponding to where T n is the noise level in µK arcmin and N ℓ is the noise power spectrum in µK 2 .Thus our final sky model is where T CMB (n) is the Websky lensed CMB temperature map, y(n) is the Websky Compton-y parameter map; T CIB (ν, n) is the Websky CIB temperature map at frequency ν; T Galaxy (ν, n) is the simulated Galactic temperature map from PySM3; g ν is the spectral response function of the tSZ effect; W (•) is the beam-convolving operation in real space (which corresponds to a filtering by the beam window function in harmonic space); and T n (ν, n) is the simulated instrumental noise.Note that all of our simulations are in units of µK CMB , i.e., units in which the CMB has unit SED; thus we need to the CIB maps, which are supplied in Jy/sr, to µK CMB .FIG. 9: Validation on the simulations.Note that the "true signal" that we compare to is estimated by using healpy.anafast()on the full-sky maps.In both cases, we have applied our NILC pipeline with various deprojections (as labeled) to the simulated microwave sky maps to build y-maps, and measured C yκ ℓ using the true κ map from Websky, with appropriate Gaussian noise added (left) to approximate the scatter in our own measurement on data, and also with no lensing noise added (right) to illustrate the biases associated with the different deprojections in the zero-lensing-noise case.Some points are slightly offset along the x axis for better visibility of overlapping points.
After applying our NILC pipeline (see Paper I for details on the pipeline) to construct y-maps from the simulated intensity maps, we measure C yκ ℓ using the simulated κ map provided by Websky.While we do not include any systematics due to the fact that the true κ is reconstructed from temperature data, we do include a realization of Gaussian noise with power spectrum matching the noise power spectrum of the 2018 Planck lensing reconstruction.We also measure the cross-correlation on a realization of the lensing map with no noise included, to demonstrate the effectiveness of the CIB moment deprojections in the limit of zero CMB lensing noise.
The resulting data points, measured with the CIB deprojections done using the default best-fit SED, are shown in Figure 9. Comparing the measurement on the left of Figure 9, which includes appropriate Gaussian lensing noise, with the measurement on the right, which is performed on the noiseless κ map, it is clear that (especially for the cases when > 1 components are deprojected) there is significant scatter due to the lensing noise.Note that the deprojection that we use for our true data points corresponds to the red points, where both δβ and δT eff CIB are deprojected.It is also clear, in both cases, that in the case when no CIB components are deprojected, there is a significant bias on the measurement.The CIB deprojection removes some of this bias, but (as is especially clear on the right) a significant amount remains.The moment deprojections remove significantly more of this bias.We investigate this further in Figure 10 where we vary the SED used for the CIB and moment deprojections to illustrate the reduction in scatter they yield.
In Figure 10, we show the signal measured on the simulations for the different moment deprojections as we vary the CIB SED.Again, we show the points measured with lensing noise on the left, and without lensing noise on the right; the biases and scatter are easier to see when there is no lensing noise.We see that deprojecting δβ indeed reduces some of the bias induced when the CIB alone is deprojected with different SEDs; the bias is removed by jointly deprojecting moments with respect to both parameters δβ and δT eff CIB , in the bottom plots, when we see scatter no greater than ≈ 1σ and an unbiased (by eye) measurement.
Our overall conclusion from these simulation tests is that a CIB+δβ + δT eff CIB -deprojected map can be used to make an unbiased estimate of C yκ ℓ , for a wide range of SED parameters assumed for the CIB modified blackbody spectrum.This reinforces our findings from the actual data in the previous section, and substantiates the overall robustness of our measurement.FIG.10: Deprojecting the CIB and various first moments thereof, measured on simulations with appropriate lensing noise (left) and no lensing noise (right).We use the same values for the CIB SED parameters as in the samples we use when we vary the SED parameters for our deprojections on the data (e.g., in Figure 5).The points are systematically offset along the x axis for better visibility of overlapping points.

VII. ANALYSIS
Our final data points, measured with all three κ maps, with the maximally-deprojected y maps (CMB 5 +CIB+δβdeprojection for the no-545 GHz case and CMB 5 +CIB+δβ+δT eff CIB -deprojection for the standard-frequency case) are shown in Figure 11.Also included on this plot, for comparison, are the data points from HS14; as the ℓ-binning schemes are similar (in particular, the widths of the ℓ-bins are identical), the error bars can be compared directly by eye.It is clear that on large scales, our measurement has lower SNR than HS14, although on small scales our error bars are slightly smaller.The data points measured without 545 GHz and with the standard frequency coverage are broadly consistent.We proceed to fit theoretical models to the points in Figure 11.
The tSZ effect is a powerful probe both of cosmology and of the physics of the ICM.Most of the signal is sourced in massive clusters in the late Universe (z < ∼ 1-1.5), which are rare objects; thus, the tSZ signal probes the tail of the HMF.This signal can be accessed through various tSZ statistics, for example through the power spectrum [18,19,49,66,79], cluster counts [80][81][82], or the 1-point PDF [28,29].In every case, while there is high sensitivity to the cosmological parameters that govern the distribution of the halos (σ 8 , which governs the amplitude of matter clustering on scales of 8 h −1 Mpc and Ω m , which parametrizes the matter density of the Universe), there is also a degeneracy with the parameters that describe the tSZ flux-mass relation (i.e., pressure-mass relation), which is required to translate from Our final data points for the tSZ -CMB lensing cross-power spectrum, measured using the CMB 5 +CIB+δβdeprojected map with 545 GHz removed (orange), and the CMB 5 +CIB+δβ+δT eff CIB -deprojected map with standard frequency coverage (blue).We show the data points for all three κ maps as indicated in the titles, and deproject the primary CMB in the first five needlet scales.We also include the data points from HS14 (green) for comparison.We note that, on large scales, our error bars are significantly larger than those of HS14, as we pay a large SNR penalty due to our many deprojections needed to robustly clean CIB contamination (the width of our multipole bins is the same as in HS14).The black curve is not a fit to the data, but rather shows our fiducial theoretical model (see Section II), with the one-halo (two-halo) term shown in dashed (dash-dotted).Some points are slightly offset along the x axis for better visibility of overlapping points.
the HMF to the flux distribution that is measured.Thus, these parameters, which are poorly constrained in regimes where ICM physics is not well understood, must be jointly constrained with the cosmological parameters.In our case, where we consider just one statistic of the tSZ effect, it is not possible to fully break the degeneracy between ICM physics and cosmology.However, as different tSZ statistics depend on various combinations of these parameters in different ways, this degeneracy can in principle be broken by combining our data with different tSZ probes, such as the auto-power spectrum.We leave such an approach to future work and instead only constrain ICM physics while holding the cosmology fixed, or alternatively constrain cosmology while holding the pressure-mass relation fixed.
The model that we constrain is described in Section II.When we fix the cosmology, and unless otherwise stated, we use the cosmological parameters from the best-fit ΛCDM model to the Planck TT+TE+EE+lowE+lensing Plik likelihood [48]: {Ω b h 2 =0.022383;Ω c h 2 = 0.12011; H 0 = 67.32km/s/Mpc; A s = 2.101 × 10 −9 ; n s = 0.96605}, where Ω b h 2 is the physical density of baryons today; Ω c h 2 is the physical density of cold dark matter today; H 0 ≡ 100h km/s/Mpc is the Hubble parameter; and A s and n s are respectively the amplitude and spectral index of scalar fluctuations with a pivot scale of 0.05 Mpc −1 .For this cosmology, σ 8 = 0.8120.

A. Likelihood for C yκ ℓ
To constrain the models for C yκ ℓ , we construct a simple Gaussian likelihood L: where with C yκ ℓ (p) the theoretical model for the signal evaluated for parameters p (calculated with class sz), Ĉyκ ℓ the cross-correlation data points, and C ℓℓ ′ the full covariance matrix.The theoretical model for the signal is a function of cosmological parameters (e.g., σ 8 and Ω m ) and pressure profile parameters (e.g., A P0 0 and A β 0 , cf.Equations ( 13) and ( 15)).
We explore the parameter posteriors with Markov chain Monte Carlo (MCMC) sampling, using the Metropolis-Hastings algorithm implementation of Refs.[83,84] in Cobaya14 [85,86].Priors on the sampled parameters are discussed below.We run all chains until they are converged with a Gelman-Rubin convergence criterion [87] of |R − 1| < 0.001 (for the P 0 fits) or |R − 1| < 0.01 (for the cosmology fits).We quote each parameter's best-fit value as the point that maximizes the posterior, which is determined using the BOBYQA [88][89][90] minimizer implemented in Cobaya.We use GetDist15 [91] to analyze our chains and plot our parameter posteriors.

B. Constraints on pressure profile parameters
We constrain the normalization of the amplitude of the pressure-mass relation, A P0 0 (see Equations ( 13) and ( 15)).For ease of notation, we redefine this amplitude parameter as simply P 0 , i.e., P 0 ≡ A P0 0 .This procedure is equivalent to simply fitting a one-parameter amplitude that rescales the entire theoretical model.In our MCMC sampling, we use a wide, uninformative, flat prior on P 0 ; we use no other priors, so the posterior is given by the log-likelihood in Equation (32).For convenience, we further rescale P 0 by its fiducial value given in Table I, i.e., we plot P 0 /P fid , where P fid = 18.1 is the fiducial amplitude of the pressure-mass relation.
The constraints on P0 P fid for different data combinations are summarized in Table III, and the posteriors for P 0 /P fid are shown in Figure 12.For comparison, we also constrain the model using the measurements from HS14.Using the HS14 data, we find a very similar result to that quoted in HS14, with In both cases we use the maximally-deprojected y maps: CMB 5 +CIB+δβ-deprojection for the no-545 GHz case, and CMB 5 +CIB+δβ+δT eff CIBdeprojection for the standard-frequency case.For all of the analyses shown here, the signal is consistent with the fiducial prediction.We find P0/P fid consistent with unity at the 1-1.5σ level in all cases except for the tSZ-deprojected κ map analysis in the right panel, which is consistent with unity at ≈ 2σ.
in this work with the counter-term prescription, which is different to the extrapolation of the HMF used in HS14, which integrated directly the contribution of all halos with 10 5 h −1 M ⊙ < M 200c < 5 × 10 15 h −1 M ⊙ .Third, we cut our halo profiles at 2r 200m , while HS14 do so at 1.5r vir .Nevertheless, these differences all lead to relatively small changes in the template signal, and hence the central value of P0 P fid in our reanalysis of HS14 is quite close to theirs.All of our constraints are consistent with each other to within 1σ.We get much better fits from the standard frequency case than the no-545 GHz case, as well as tighter constraints.The tightest constraint comes from the standard-frequency-coverage map in combination with the NPIPE κ map, resulting in P0 P fid = 0.82 ± 0.21.The SZdeprojected results are slightly (≈ 1σ) slower, resulting in P0 P fid = 0.56 ± 0.24, perhaps hinting at a ⟨yyy⟩ bias; the SZ-clusters-restored measurement is very slightly higher than the NPIPE measurement.It would be interesting to perform the measurement on a tSZ-deprojected + cluster-signal-restored NPIPE map, to see if these conclusions hold.
It is notable that our final constraints are not tighter than those from HS14, with all of our error bars comparable to or larger than theirs, despite the higher-quality data we use.However, our results are significantly more robust to potential CIB contamination due to the battery of deprojection operations implemented in our analysis, which are also responsible for inflating our error bars significantly.We show posteriors in Figure 12 for all combinations of our fully-deprojected y-maps (i.e., CIB+δβ+CMB 5 -deprojected for the no-545 GHz case and CIB+δβ+δT eff CIB +CMB 5deprojected for the standard-frequency case) with the three κ maps we consider.We also list the best-fit values in Table III.It is generally true that the tSZ-deprojected κ map yields a lower measurement (by ≈ 1σ) of P 0 than the other two maps.However, they are all statistically consistent.It is also generally true that including 545 GHz and deprojecting δT eff CIB improves the quality of the fit (and also decreases the error bars), with the reduced χ 2 values decreasing in every case.

Including a correction from hydrodynamical simulations in the theory
In Appendix A, we compare our halo model to a measurement of C yκ ℓ from the hydrodynamical simulations from which our P − M relations were calibrated (see Figure 16).It is clear that some portion of the signal is unaccounted for by the halo model in the regime where we are measuring C yκ ℓ .We can incorporate this by measuring the ratio of the halo model prediction to the simulation measurement (as calculated for cosmological parameters corresponding to those in the simulation, as we do in Appendix A).Under the assumption that this ratio is cosmology-independent, we can then rescale the halo model prediction at the Planck cosmology we are using in this work, in order to include this contribution.When we do this, and constrain an amplitude for the model, the posteriors are shown in Figure 13.As expected, the posterior means are shifted to lower amplitudes, since the fiducial model prediction is now higher.with the exact values of the power-law indices depending on the multipole ℓ.Importantly, this is a different cosmological dependence to that of the auto-power spectrum, which displays the following dependence16 : As the relative dependences between cosmology and astrophysics (represented by P 0 in Equations ( 34) and ( 35)) are different, this means that while both probes separately cannot separate cosmology and astrophysics, a combined analysis can break this degeneracy (see, e.g., HS14).However, in this work we focus on constraining cosmology only from the cross-power spectrum, leaving a joint analysis to future work.Due to the degeneracy between σ 8 and Ω m , we can only constrain simultaneously the product of these parameters.We find the tightest constraints on the parameter combination σ 8 (Ω m ) ≈0. 25 , which is very close to the best-determined parameter combination σ 8 (Ω m ) 0.26 found in HS14.When we vary the cosmological parameters, we only vary Ω cdm and σ 8 , with H 0 held fixed to H 0 = 67.32km/s/Mpc and Ω b h 2 fixed to 0.02383.We impose a flat prior on σ 8 (0 < σ 8 < 2) and a flat prior on Ω cdm with an upper limit that corresponds to 0 < Ω m < 1.We thus obtain Ω m as a derived parameter in the chains.

VIII. DISCUSSION
In this work, we have used our CIB-cleaned y-maps from Ref. [2] to measure and constrain the cross-correlation of the tSZ effect and CMB lensing in our Universe, with Planck PR4 data.The tSZ effect has been cross-correlated with low-z tracers of LSS (in particular galaxy catalogs and weak lensing maps) many times (e.g.[43,46,47,93,94]); however, the cross-correlation with CMB lensing has only been measured once before, in HS14.The difficulty in this measurement lies in the fact that at the high redshifts probed by CMB lensing, the tSZ signal is subdominant in our mm-wave sky maps to the CIB.Much of our work has thus focused on removing the CIB from this measurement in a model-independent way.The result is a less statistically significant, but much more robust, measurement than that of HS14.
We have chosen to remove the CIB by making constrained ILC maps from the Planck single-frequency maps.It has been standard for previous y-LSS cross-correlations to neglect the CIB (for low-z tracers), or to remove only a component corresponding to the CIB SED, by modelling the CIB as a modified blackbody with fixed parameters β and T eff CIB describing its SED (e.g., in [46,47]).However, due to the uncertainty on the values of these parameters for the CIB, we have extended this approach by measuring our cross-correlation on moment-deprojected maps [39].This approach allows us to account for errors incurred due to choosing incorrect parameter values for the SED, as well as for the spatial variation and line-of-sight variation of the CIB SED.
We have verified extensively that such an approach should measure an unbiased signal, by applying our NILC pipeline to detailed simulations of the microwave sky, incorporating appropriate inter-frequency decorrelation of the CIB.We have checked that this approach works on the simulations regardless of the SED parameters chosen to deproject the CIB.We have also tested that our measurement on the actual data is stable to the choice of parameters used for the CIB SED.We report our measurements with two different y-maps: one with all frequency channels up to 545 GHz used in the NILC and with the CIB and its two first-order moments δβ and δT eff CIB deprojected, and one without 545 GHz used, which has the CIB and only δβ deprojected.Both maps have the CMB deprojected in the first five needlet scales to avoid ISW bias in the measurement.
As such, we have made the most robust CIB-free, unbiased measurement of the tSZ-κ cross-correlation to date.We have made the measurement with three different choices for the κ map used in the cross-correlation, each chosen either to minimize a certain systematic (a κ map made with tSZ-deprojected T maps to minimize residual ⟨yyy⟩ bias; a κ map which restores signal in the regions of the sky where the tSZ clusters are originally masked in the reconstruction; and the PR4 NPIPE κ map, which has the largest signal-to-noise).Our measurements are broadly consistent, both on a point-by-point basis and in the overall constraint we make on the amplitude of the theory model, although the tSZ-deprojected amplitude is lower than the other amplitudes by ≈ 1σ.Overall, our highest-significance detection of the cross-correlation is at the level of ≈ 3.6σ, with the PR4 κ map and the no-545-GHz CMB 5 +CIB+δβ-deprojected y-map.Note that if the data points perfectly matched the fiducial theory template, the SNR in this case would reach ≈ 5σ.
While our detection is not currently at extremely high significance, we note that the noise in the κ map is a limiting factor that will be lowered very soon.Indeed, it will be extremely interesting to make this measurement with the much lower-noise-per-mode ACT DR6 CMB lensing map [8,24].Additionally, improvements in the small-scale measurement of the y-map [22] using the ACT DR6 data will allow us to probe this signal on smaller scales than we have in this work.The methods to remove the CIB developed here will be of great use for such a measurement, as well as other future LSS-y cross-correlation measurements, which will need to remove CIB bias (see, e.g., [95]).Further exploration of the trade-off between ILC deprojections (to mitigate CIB biases) and the measurement SNR will also be of interest, e.g., using the methods of [96].Such a measurement will allow us to learn more about currently unprobed astrophysics in the ICM of high-z groups and clusters.
Finally, we note that it would be of great interest to combine an analysis of C yκ ℓ with a measurement of the autopower spectrum (or indeed of any other statistic of the y-map), to break the degeneracies between cosmology and astrophysics and to access the large potential that the y signal contains to constrain cosmology.x cut = 1 x cut = 2 x cut = 2.5 measurement from Ref. [64] FIG.16: A comparison between our halo-model calculation of the signal and the direct measurement from simulations in Ref. [64].The halo model calculation for our fiducial cutoff choice (xcut = 2) is shown in orange.It is clear that we miss some power on large scales, which arises from diffuse, unbound gas outside of halos and/or contributions from low-mass halos.
the total amount of matter in the Universe (i.e., the counter-terms required to account for lower-mass halos become negative to compensate).17: The measured data points on different sky areas.In all cases, we use the maximally-deprojected y-maps (CMB 5 +CIB+δβ+δT eff CIB for standard frequencies and CMB 5 +CIB+δβ-deprojected for the no-545 GHz case), and the tSZdeprojected κ map.Overall, we do not see a systematic variation in the data points as we decrease the sky area (and thus measure on regions of the sky that are less impacted by Galactic foregrounds).On the left is the standard frequency coverage case, and on the right is the no-545-GHz case.

CFIG. 1 :
FIG. 1: The contribution to C yy ℓ and C yκ ℓ from different redshifts.C yy ℓ FIG. 2: C yy ℓ and C yκ ℓ as calculated with different maximum masses in the halo model integral.The C yy ℓ signal is highly concentrated in the highest mass range with M > 10 14.5 h −1 M⊙.The C yκ ℓ signal probes lower-mass halos, with significant contributions down to M > 10 13.5 h −1 M⊙.The units of M in the legend are h −1 M⊙.

fsky = 80% threshold mask 0 1 Planck 1 SZ 1 SZ 1 NPIPE 1 SZ 1 SZ 1 NPIPE
HFI+LFI point source masks: fsky = 93.86%0 -deprojected κ analysis mask: fsky = 66.28% 0 -clusters-restored κ analysis mask: fsky = 67.37%0 κ analysis mask: fsky = 67.02%0 -deprojected combined mask (10 apodized): fsky = 61.57%0 -clusters-restored combined mask (10 apodized): fsky = 62.74% 0 combined mask (10 apodized): fsky = 62.24% 0 1 use our fiducial theoretical model (see Section II) to calculate the (sub-dominant) C yκ FIG.4:The tSZ -CMB lensing cross-power spectrum, C yκ ℓ , as measured for different frequency coverage (left panel including all frequencies from 30 -545 GHz and right panel with 545 GHz removed) and different deprojection combinations (as labeled in the legend).Note that, for the case when we exclude 545 GHz, we do not have enough frequency channels to simultaneously deproject the CIB and both δβCIB and δT eff CIB on the smallest scales; it would be possible to do so on large scales, but we do not consider such a measurement.The solid black curve shows our fiducial theoretical model, with the one-halo (two-halo) term shown in dashed (dash-dotted).Note that, in the legend, CMB 5 refers to the CMB deprojected on the largest five needlet scales.Also note that some points are slightly offset on the x axis to allow for better visualization.Both plots were made with the SZ-deprojected reconstruction of the κ map.

FIG. 8 :
FIG.8:The tSZ -CMB lensing cross-correlation measured on different κ maps, as labeled in the legend.The left (right) panel uses a NILC y-map built with the CMB 5 +CIB+δβ+δT eff CIB deprojection applied to the standard frequency coverage (CMB 5 +CIB+δβ deprojection applied to the no-545-GHz frequency coverage).The measurements are broadly consistent for all three κ maps.Some points are slightly offset along the x axis for better visibility of overlapping points.
for varying x cut FIG.17:The measured data points on different sky areas.In all cases, we use the maximally-deprojected y-maps (CMB 5 +CIB+δβ+δT eff CIB for standard frequencies and CMB 5 +CIB+δβ-deprojected for the no-545 GHz case), and the tSZdeprojected κ map.Overall, we do not see a systematic variation in the data points as we decrease the sky area (and thus measure on regions of the sky that are less impacted by Galactic foregrounds).On the left is the standard frequency coverage case, and on the right is the no-545-GHz case.

TABLE I :
The best-fit values from Ref.