Cross Correlation between the Thermal Sunyaev–Zeldovich Effect and the Integrated Sachs–Wolfe Effect

We present a joint cosmological analysis of the power spectra measurement of the Planck Compton parameter and the integrated Sachs–Wolfe (ISW) maps. We detect the statistical correlation between the Planck thermal Sunyaev–Zeldovich (tSZ) map and ISW data with a significance of a 3.6σ confidence level (CL), with the autocorrelation of the Planck tSZ data being measured at a 25σ CL. The joint auto- and cross-power spectra constrain the matter density to be Ωm=0.317−0.031+0.040 , the Hubble constant to be H0=66.5−1.9+2.0kms−1Mpc−1 , and the rms matter density fluctuations to be σ8=0.730−0.037+0.040 at the 68% CL. The derived large-scale structure S 8 parameter is S8≡σ8(Ωm/0.3)0.5=0.755±0.060 . If using only the diagonal blocks of covariance matrices, the Hubble constant becomes H0=69.7−1.5+2.0kms−1Mpc−1 . In addition, we obtain the constraint of the product of the gas bias, gas temperature, and density as bgasTe/(0.1keV)n¯e/1m−3=3.09−0.380+0.320 . We find that this constraint leads to an estimate on the electron temperature today as Te=(2.40−0.300+0.250)×106K , consistent with the expected temperature of the warm–hot intergalactic medium. Our studies show that the ISW–tSZ cross correlation is capable of probing the properties of the large-scale diffuse gas.


INTRODUCTION
Understanding the composition and distribution of baryonic matter in the Universe is a crucial step toward unraveling the mysteries of its formation and evolution.While baryons make up a small fraction of the total matter-energy content (Aiola et al. 2020;Planck Collaboration et al. 2020), they play a vital role in the processes that shape the cosmic structures that we observe today.Early estimates and numerical simulations show that most baryons are "missing," whereas the baryons that are already made into stars and galaxies constitute a small portion of the total baryon budget (Fukugita et al. 1998;Fukugita & Peebles 2004;Cen & Ostriker 2006;Shull et al. 2012).Accurate constraints on the missing baryons are valuable to improve our understanding of the fundamental parameters that govern the Universe.
To this aim, several techniques have been adopted.For instance, Nicastro et al. (2018) and Nevalainen et al. (2019) used quasars to search for absorption features that correspond to baryons in the form of neutral gas and the intracluster medium.The dispersion measure of fast radio bursts (FRBs) provides an integrated measurement of baryon density (Ω b h 2 ) out to localized host galaxies at high redshifts (Macquart et al. 2020;Yang et al. 2022).Quantified deuterium abundance in the high-redshift metal-poor damped Lyα system (DLAs) can reveal the primordial deuterium abundance, which leads to precise measurements of baryon density (Cooke et al. 2018).
The tSZ effect (Sunyaev & Zeldovich 1972) is produced by the inverse Compton scattering of the cosmic microwave background (CMB) photons by high-energy electrons.It is sensitive to the properties of the hot gas in galaxy clusters, including its density and temperature.Therefore, by measuring the tSZ effect, we can study the distribution and physical properties of the ionized baryons, which are important for understanding the large-scale structure (LSS) of the Universe and the processes of cluster formation and evolution.The effect imprints a distortion in the CMB blackbody spectrum, which can be expressed as where S SZ (x) = x coth(x/2)−4 is the spectral distortion function, x ≡ hν/k B T 0 , and T 0 = 2.725 K is the background CMB temperature.The Compton-y parameter in Equation ( 1) is the integral of the electron pressure along the line of sight: where m e , T e , and n e are the mass, temperature, and number density of electrons, respectively.σ T , k B , and c are the Thomson scattering cross section, the Boltzmann constant, and the speed of light, respectively.Unlike X-ray luminosity, the tSZ signal is linearly proportional to the baryon density and the electron tem-perature, making the detection of gas with low densities and medium temperatures possible.To achieve this aim, stacking analyses have been used at the locations of clusters and superstructures to achieve the goal of detecting such gas.Recently, Tanimura et al. (2019) stacked ∼ 10 5 pairs of luminous red galaxies on the Planck tSZ map and obtained the first detection of warm gas along the filamentary structures.A complementary study of the gas density and temperature was performed in Tanimura et al. (2020), by stacking 24,544 filaments (across 18-30 Mpc scales) identified by the DisPerSE method (Sousbie 2011).de Graaff et al. ( 2019) conducted a similar study with CMASS galaxies at higher redshifts and obtained results consistent with the predicted WHIM properties.Other relevant studies can be found in Mittaz et al. (1998), Eckert et al. (2015), Vavagiakis et al. (2021), Kusiak et al. (2021) and Bonamente et al. (2022).
Similarly, the tSZ auto-angular power spectrum and cross-correlation with other LSS tracers have been extensively studied to trace gas distribution, calibrate cluster masses, and characterize cosmic structures on large scales.This is due to the particular advantage of the tSZ effect of not being very sensitive to redshift evolution, which leaves the possibility of cross-correlating with other LSS tracers to pull out the signal of interest.For example, Van Waerbeke et al. (2014) crosscorrelated the tSZ map with the Canada-France-Hawaii weak lensing data and obtained a ∼ 6σ confidence level (CL) detection.Ma et al. (2015) and Hojjati et al. (2017) found that the detected signal in Van Waerbeke et al. (2014) corresponds to ∼ 50% of the baryons lying beyond the virial radius of the halos with a temperature of 7 × 10 5 K ≤ T ≤ 3 × 10 8 K. Hurier et al. (2019) combined the tSZ effect, X-rays, and weak-lensing auto-and cross-correlation angular power spectra to estimate the cluster mass bias and obtained a result consistent with the Planck measurement at 2σ CL.Ma et al. (2021) used the tSZ-lensing cross-correlation to constrain the pressure profile of galaxy clusters.A series of other studies based on cross-correlating the tSZ results with other tracers to measure the cluster mass bias can be found in Ma (2017), Bolliet et al. (2018), Li et al. (2018), Makiya et al. (2018), Salvati et al. (2019), Koukoufilippas et al. (2020) and Ibitoye et al. (2022).
The cross-correlation technique equally allows us to study the statistical relationship between LSS tracers and other cosmological probes, providing valuable insights into the underlying cosmological parameters such as the Hubble parameter (h), matter density (Ω m ), and the amplitude of matter density fluctuations (σ 8 ).For example, Osato et al. (2019) cross-correlated the tSZ maps from Planck with the Subaru Hyper Suprime-Cam (HSC) year one data and achieved the constraint of Ω m = 0.3149±0.008and σ 8 = 0.8304±0.014(68% CL.).Tröster et al. (2021)  Our study pursues a similar scientific goal.We constrain the cosmological parameters including the Hubble parameter, matter density, and the amplitude of matter density fluctuations, together with the astrophysical parameter W SZ (the product of the mean electron density ne , electron temperature T e , gas bias b gas at redshift z = 0).To achieve this, we cross-correlate the tSZ effect with the integrated Sachs-Wolfe (ISW) effect (Sachs & Wolfe 1967), which characterizes the largest-scale perturbations in the Universe, and use this to provide valuable constraints on WHIM properties and the underlying cosmology.The ISW effect is a secondary anisotropy of the CMB, which arises when a CMB photon from the last scattering surface enters and leaves the timeevolving gravitational potential, which can change its net energy.Therefore, it causes an additional temperature anisotropy on large scales and is sensitive to the growth rate of cosmic structures.Similarly, because the tSZ effect is also sensitive to the growth rate of cosmic structures, the cross-correlation of the ISW with the tSZ effect can provide useful insight into the gas properties on very large scales (Creque-Sarbinowski et al. 2016).
We will employ the ISW and tSZ maps obtained with the Planck satellite for this analysis.The Planck Collaboration detected the ISW signal, with significances from ∼ 2.5σ to ∼ 4σ CL, by cross-correlating the Planck CMB map with radio sources from the NVSS catalog, galaxies from the optical Sloan Digital Sky Survey (SDSS) and the Wide-field Infrared Survey Explorer (WISE) surveys, and the Planck 2015 convergence lensing map (Planck Collaboration et al. 2014a, 2016a).In addition, the Planck Collaboration also provided different maps of the ISW fluctuations, from different combinations of the abovementioned LSS tracers.We will consider the potential systematics from the cosmic infrared background (CIB) in the cross-correlation, which can contaminate both ISW and tSZ maps at z < 1.We aim to conduct a joint analysis of both the ISW and tSZ effects and constrain the aforementioned cosmological parameters and gas quantities on very large scales.
This paper is arranged as follows.In Sec. 2, we detail the observables, including the theoretical models and  data sets employed.In Sec.3, we present the detailed cross-correlation analysis and demonstrate the significance of the cross-correlated signal.In Sec. 4, we carry out the power spectrum analysis.In Sec. 5, we describe the parameter constraints and their cosmological inference.We present our conclusions in Sec. 6.Throughout the paper, we assume a spatially flat ΛCDM model in which, apart from the parameters being varied in Table 2, we fix the other cosmological parameters to be n s = 0.965, τ = 0.0540, and ln(10 10 A s ) = 3.043 (Planck Collaboration et al. 2020).

ISW Map
The ISW effect (Sachs & Wolfe 1967) is a secondary anisotropy of the CMB caused by the time-varying gravitational potentials of the LSS, which dominates the CMB on large scales.At low redshifts, it shows a direct signature of dark energy in a ΛCDM Universe because of the potential decay effect.ISW is also sensitive to the evolution of the growth factor, so it can test alter- native theories of gravity along with other cosmological probes.However, because the ISW signal is mostly on large scales that are dominated by cosmic variance, the number of modes extractable is limited (Maniyar et al. 2019).Thus, the ISW effect is only detectable by crosscorrelating with other LSS observations (Crittenden & Turok 1996;Planck Collaboration et al. 2016a).Figure 1 is the Planck's ISW map delivered in HEALPix 1 format (Górski et al. 2005) with an original pixel resolution N side = 64.We notice that the resolution of the ISW map is 160 ′ in full width at half maximum (FWHM).map has a negligible contribution from foreground residuals.

Compton parameter map
In this analysis, we employ the full-sky Comptony map provided by the Planck satellite mission (Figure 2;Planck Collaboration et al. 2016b).
The Compton-y map is the result of the combination of Planck individual frequency maps convolved to a resolution of 10 ′ using either the Modified Internal Linear Combination Algorithm (MILCA) method (Hurier et al. 2013) or the Needlet Independent Linear Combination (NILC) algorithm (Remazeilles et al. 2011).Because these two maps are consistent in the limit of uncertainties, for the rest of this work we only present the results for the NILC map.To block out spurious contamination and Galactic foregrounds that can affect our results, we combine the 40% Galactic mask with the point-source To match the resolution of the ISW map, we deconvolve the Compton-y map with a 10 ′ Gaussian beam and then convolve with a 160 ′ Gaussian beam to make it the same resolution with the ISW. 3 We then downgrade the Compton-y map to a coarser pixelization of N side = 64, the same as the ISW map.We will use this reprocessed y-map in all subsequent analyses.

Measurement
We now describe how we obtained our auto-and crosscorrelation measurements.The Planck ISW and Compton y-parameter maps described in Sections.2.1 and Sec.2.2 are used to measure the projected autocorrelation angular power spectrum of each observable and the cross-correlation angular power spectrum of the two.This analysis is conducted using the publicly available pseudo-C ℓ MASTER algorithm (Hivon et al. 2002) implemented in the NaMaster package4 (Alonso et al. 2019).The code permits dealing with spin 0 and 2 signals, e.g., the case of the CMB temperature and polarization.However, the case considered in this work always deals with scalar quantities.The software is efficient in handling complex masks and high-resolution pixelization.It also executes E-or B-mode purification in both flat-sky approximation and in full-sky mode.We then made a sanity check with the anafast subroutine included in the HEALPix software package (Górski et al. 1998;Zonca et al. 2019).The power spectra measured using both software packages are in good agreement.
By construction, NaMaster accepts all sky maps in the form of HEALPix maps exclusively with RING ordering and produces the C T1T2 ℓ power spectrum if two spin−0 fields are given, C T1E2 ℓ , C T1B2 ℓ for one spin−0 field and one spin > 0 field, and power spectra for two given spin > 0 fields (the subscripts "1" and "2" represent the field or the observable of interest).Hence, since both the ISW and  7).The diagonal blocks show the standard correlations for the TT, yy, and Ty spectra; the offdiagonal blocks show the additional cross-correlations between different spectra.The correlation values are shown for pairs of the effective band power multipoles defined in Table 1.
Compton parameter maps are given as spin−0 fields, the measured auto-and cross-angular power spectra5 are C TT ℓ , C yy ℓ and C Ty ℓ on full sky respectively.To compare them with the theoretical prediction, we show in Figure 3 the resulting y-auto (yy), ISW-auto (TT), and ISW-y cross-power spectra (Ty).
The sensitive ℓ-range in this study is ℓ ∈ [16,190].The lower cut is due to the ISW map not having modes for ℓ < 8 (Sachs & Wolfe 1967) and the significant cosmic variance between ℓ ∈ [8, 16], which can bias our result.The high-ℓ cut is due to the limited resolution of the ISW map (N side = 64).The data points we show in Figure 3 are the binned-ℓ ones, to reduce the empirical scatter between adjacent multipole amplitudes ℓ.We adopt the same scheme for binning as in Planck Collaboration et al. (2016c).In addition, this makes the power spectrum coupling matrix invertible to some extent.One may notice that the Planck Collaboration et al. (2016c) used 19 multipole bins in total, but we used 10 bins for the range of ℓ eff = 18.0 to ℓ eff = 181.0 in total.The binning is done by evaluating each bin's average power spectrum values.The resulting values are quoted in Table 1.
We remind the reader that the y-auto power spectrum presented in the upper left of Figure 3 was measured in full sky while the one presented on the left of Figure 11 (2016b) shows that the C yy ℓ measurement is greatly affected by Galactic emissions below ℓ < 20, but not so much for ℓ ≫ 20.Therefore, it is not surprising to see the consistency between our study and the Planck results for 18.0 < ℓ < 152.5 (see Appendix C for more details).

Covariance matrix calculation
An accurate estimate of the Gaussian covariance can be calculated using NaMaster.We then generated 1000 realizations of the Compton-y map from the fiducial power spectrum.Each of these maps is supplied as a spin-0 field to the code with the same weighting/masking scheme, to output binned power spectra C i ℓ (i.e.Ĉyy,i ℓ ); the resultant covariance can then be calculated from the binned power spectra.
This covariance matrix estimated from NaMaster is Gaussian and does not account for correlation between different multipoles.However, hydrodynamical simulations showed that SZ fluctuations can indeed be non-Gaussian (Seljak et al. 2001;Zhang et al. 2002).Therefore, to obtain accurate constraints on cosmological parameters, it is important to include the non-Gaussian nature of the SZ fluctuations.To this aim, we calculated the connected non-Gaussian covariance matrix using the pyccl6 code (Chisari et al. 2019) to capture the coupling between the measurements at multipoles ℓ and ℓ ′ .This makes the offdiagonal elements nonzero, by an amount that depends on the parallel configurations of the connected trispectrum.The resultant error bar plotted in the upper left panel of Figure 3 is then calculated from the sum of the Gaussian and non-Gaussian covariance matrix in the SZ solution.
On the other hand, to accurately account for the statistical information encoded in the ISW sky map, and its cross-correlation with the Compton map, we used Monte Carlo (MC) simulations to generate a set of N sim ISW sky maps (where N sim = 900), as described in Section 3.2.From each of these simulated maps, we output the binned power spectra as ĈTT,i ℓ and their respective cross-correlation power spectra with the tSZ map as ĈTy,i ℓ .We then compute the covariance matrices using ). (3) The average of the resampled power spectra is given by where A stands for either TT, yy, or Ty, and N band = 10 is the number of bins in ℓ-space.
To maximize and optimize the information encoded in the measurements, we concatenate all three measured power spectra into a vector: Then the total covariance matrix is a 3×3 block matrix, and each block is an N band × N band matrix, as Cov TT,Ty T Cov Ty,Ty Cov Ty,yy Cov TT,yy T Cov Ty,yy T Cov yy,yy where, for instance, the Cov TT,Ty is the cross-covariance obtained using ĈTT,i ℓ and ĈTy,i ℓ in Equation (3).To evaluate the correlation between different blocks of the covariance matrix, we calculate the correlation coefficient matrix with Equation (6) using the expression We show the correlation matrix in Figure 4.One can see that the offdiagonal matrix blocks, i.e. (TT,yy), (Ty,yy), and (TT,Ty) have much smaller values than the diagonal blocks, indicating a small correlation between the power spectra C yT ℓ , C yy ℓ , and C TT ℓ .We can also calculate the significance of our observed spectra against the null detection as S/N ≡ χ 2 , where and the signal-to-noise ratios (S/Ns) for Ty and yy are 3.6σ and 25σ, respectively, with no significance on the TT correlation (See Section 3.2 for more details).In Sec.5.4, we see that using C TT ℓ only is very weak for constraining cosmological parameters.

CROSS-CORRELATION DETECTION
In this section, we summarize the approach adopted to characterize the tSZ and the ISW by means of a null hypothesis test, by comparing the observed crosscorrelation for the data with that obtained from the actual tSZ and a set of simulated ISWs constrained in the same way as the ISW data.This detection approach is similar to the one used in the literature to detect the ISW effect (e.g., Boughn & Crittenden 2004;Vielva et al. 2006).

Simulations
The ISW simulations used to characterize the null hypothesis are the ones described in Section 5 of Planck Collaboration et al. (2016a).These simulations contain the statistical properties of the Planck ISW map and are built by applying the linear covariance-based (LCB) filter introduced in Barreiro et al. (2008) and extended in Manzotti & Dodelson (2014) and Bonamente et al. (2022).They start from coherent galaxy density number simulations of different LSS tracers (the same used to derive the Planck ISW map), which are afterward combined with the LCB to map out the expected ISW anisotropies.Among the different survey combinations described in Planck Collaboration et al. (2016a), we adopt the one containing information from the CMB, and the NVSS, SDSS, WISE, and Planck lensing LSS tracers (see Section 5 in Planck Collaboration et al. 2016a, for details).We used 1000 simulations generated in this way to characterize the mean value and the full covariance of the tSZ-ISW cross-correlation, as explained in the next subsection.

Detection
To quantify the statistical significance of our measurement, we perform a robust null hypothesis test from the MC simulations described in 3.1.We remark that with our simulation package, we produced 1000 simulated ISW sky maps.Among these, 900 were used to compute the covariance matrix, as detailed in Section 2.4, while the remaining 100 (name it, N sig ) were used to quantify the significance as follows.We first measured the spectra from these maps and grouped them into a matrix with size N sig × N band dimensions.We then computed the χ 2 associated with each of these power spectra and their respective S/Ns7 using Equation (8) (see Figure 5), from which we estimated their mean and dispersion (rms value).We remark that when such a procedure is followed using the ISW simulations, we are getting a distribution of this statistic (S/R) for the null hypothesis (e.g., there is no correlation).For the crosspower spectrum measurement, the true S/N is 5.63σ (using the ISW sky map obtained from Planck), shown as a black solid line in the lower left panel of Figure 5 and as a black dashed line in the lower right panel of Figure 5.The realizations from 100-simulations are distributed with S/N = 3.12 ± 0.7.Therefore, under the no-correlation hypothesis, the true data deviate from the simulated ones with (5.63-3.12)/0.7 ≃ 3.6σ.This implies that the no-correlation hypothesis is rejected by 3.6σ (≃ 99.97% CL.), indicating the significance of true ISW-tSZ cross-correlation is achieved at 3.6σ CL.

POWER SPECTRUM ANALYSIS
One can see from Figure 3 that the measurements of the ISW auto-power spectrum and the ISW-tSZ crosspower spectrum are mainly on large scales.This is because the ISW effect primarily exists at low multipoles, due to the late-time decay of the gravitational potential.For this reason, the ISW map from Planck has a low resolution (θ FWHM = 160 ′ ) that smears out the structures beyond ℓ ≳ 100.Therefore, in the following, we will calculate the theoretical power spectrum based on large-scale gas bias models.

Compton-y parameter
We adopt a large-scale gas bias model for the tSZ effect.Following Goldberg & Spergel (1999) and Van Waerbeke et al. (2014), the gas density contrast is given by δ gas (x) = b gas (z)δ m (x = χn, χ(z)), where δ m is the matter density contrast and b gas (z) ∝ (1+z) −1 is the gas bias.Since we correlate the fluctuation of gas density, we have: where ne0 and b gas,0 are the electron number density and the gas bias at the present day, respectively.In this model of large-scale gas bias (Goldberg & Spergel 1999;Van Waerbeke et al. 2014), the spatial fluctuation of the gas temperature is ignored and the average temperature is proportional to a, i.e., T e ∝ a.Therefore, for the electron temperature, we have T e (z) = T e (0)a.Substituting Equation (9) into Equation (2), we have where we have substituted the relation between differential proper distance and the comoving distance dl = adχ.W SZ becomes a constant factor, where In the following, we will omit "0" in the subscripts for brevity.

Temperature fluctuation of ISW effect
The temperature fluctuation caused by this effect is (Cooray 2002;Taburet et al. 2011): The Poisson equation relates the gravitational potential with the density fluctuation, which can be inverted to calculate the potential value (see, e.g.Creque-Sarbinowski et al. 2016): where ∇ x is the gradient in comoving coordinates.

Auto-and cross-power spectra
The TT angular power spectrum or the temperature correlation represents the mathematical breakdown of the anisotropic temperature map using spherical harmonics analysis.The tSZ auto-power spectrum is calculated by carrying out the spherical harmonics decomposition on the sky (Planck Collaboration et al. 2014b;Van Waerbeke et al. 2014;Ma et al. 2015;Creque-Sarbinowski et al. 2016): where we have defined, in the second equality, ∆ SZ (z) ≡ W SZ /χ.P m ((ℓ + 1/2)/χ, z) is the linear matter power spectrum, which we compute using the camb code (Lewis et al. 2000).
Similarly, we use Limber approximation to derive the ISW auto-power spectrum weighted by the ISW kernel (Creque-Sarbinowski et al. 2016): where is the redshift-dependent kernel of the ISW. 8The crosspower spectrum of the ISW and tSZ becomes where we have used the definition below Equation ( 15).We note that both Equations ( 18) and ( 15) cannot be compared directly with the measurements described in Sections 2.1 and 2.2, because they do not represent the exact information from maps until we factor in the foreground contributions in each data set.

Modeling of the contaminations
The power spectra of the yy and Ty both suffer from contamination.Due to the imperfect cleaning of the Compton-y map by the component separation algorithms, the resulting map contains residual foregrounds from the clustered CIB, emission due to the infrared (IR) point sources, radio point-source emission, and the instrumental correlated noise. 9These foreground components must therefore be accounted for in the y power spectrum calculation.We assume that the residual foregrounds in the map keep the same shapes as 8 Our definition of the kernel has an additional 1/D(z) factor to Equation (6) in Creque-Sarbinowski et al. ( 2016), because we use the power spectrum at redshift z. 9 The instrumental noise is called "correlated" because certain ℓdependent features are introduced in the noise due to the NILC algorithm procedure, which makes the noise no longer white.their original power spectra but with reduced amplitudes, due to the cleaning process.Therefore, we model the residual foregrounds by scaling them with a set of amplitudes referred to as nuisance parameters, i.e., (A CIB , A IR , A RS , A CN ).Then the measured C yy ℓ power spectrum can be written as (see also, e.g.Makiya et al. 2018;Ibitoye et al. 2022): where C tSZ ℓ is the theoretical tSZ effect power spectrum.
are the CIB, IR sources, radio sources, and the correlated noise spectra that may remain in the Compton-y map, which are tabulated in Bolliet et al. (2018).These foregrounds mainly dominate on small scales.For example, the correlated noise power spectrum may only match the yy spectrum am-plitude at the multipoles ℓ ≃ 2742 (Ibitoye et al. 2022).In our analysis, because of the resolution of the y-map N side = 64, the sensible ℓ-range is up to ℓ max ≃ 192 (see Figure 3).Thus, the contribution from the instrumental correlated noise is negligible.Similarly, we neglected the IR and radio point-source components because they are also subdominant on large angular scales and would only contribute significantly on scales beyond ℓ ∼ 200 (see Table 1 of this paper and also Table 3 in Bolliet et al. 2018 and Table 3 and Figure 6 in Makiya et al. 2018).Therefore, to recover the observed yy power spectrum, the CIB power spectrum weighted by the A CIB coefficient would be sufficient.Hence, Equation ( 19) is reduced to One can see from the upper left panel of Figure 3 that for the low-resolution y-map we are using, the CIB contributes the most at ℓ ≥ 50, where the other foregrounds are subdominant.We have also verified that including other foregrounds does not necessarily improve our results (see Figure 10).For these foregrounds, we use the same power spectrum templates used in the original Planck analyses, Planck Collaboration et al. (2014b) and Planck Collaboration et al. (2016a) which is also shown in Table 3 of Bolliet et al. (2018).
We also calculate the redshift dependence of the three cross-correlation power spectra and see what range of redshifts they are most sensitive to.We take the logarithmic derivative of the three power spectra with respect to the redshift: where we have substituted ∆ XY ℓ (z) = ∆ X ℓ (z)∆ Y ℓ (z).X, Y can be either the ISW (auto) or SZ (auto) or both (cross).In Equation ( 21), we fix the model parameters at the best-fitting values in Table 2.We show the redshift dependence in Figure 6.One can see that although the redshift distributions for the SZ and ISW auto-power spectra are comparable, that of ISW auto-power spectrum falls off at z ∼ 3. On the other hand, the tSZ-ISW cross-power spectrum has a broad distribution, mainly covering 0.1 < z ≤ 10.
Because of the dependence of the ISW on the redshift regime z ≥ 1, it is necessary to consider the effect of CIB contamination (extragalactic dust) in the yT cross-correlation (Makiya et al. 2018).We estimate the level of CIB contamination in the yT crosscorrelation by using the three individual frequency CMB maps at 353, 545, and 857 GHz from the Planck legacy archive.We degraded the resolution of each CIB map from N side = 2048 to N side = 64 to match the pixelization of the ISW.We remind the reader that the CIB pixelized maps are at a resolution of 10 ′ .Therefore, we also deconvolve the CIB maps with a 10 ′ Gaussian beam and then convolve them with a 160 ′ Gaussian beam to match the ISW map.We then measure the cross-correlation between all three CIB maps and the ISW map and show them in Figure 7.We find that the cross-correlation of each of the CIB frequency maps with the ISW map is positive, as expected.This is because at z > 1, where ISW is relevant, the CIB is non-negligible.Similarly, the CMB photons that travel through the time-evolving gravitational potentials that underlie the large-scale cosmic structures are the sources of the CIB (Ilić et al. 2011).As shown in Figure 7, the amplitude of the cross-correlation is sensitive to the frequency, which may be because the CIB probes galaxies on low-frequency bands at higher redshifts, while ISW is sensitive to galaxies at low redshifts (Maniyar et al. 2019).Therefore, since the power spectra of the CIB-ISW cross-correlations are very similar in shape but with different amplitudes, we compute the average of the three power spectra, C CIBavg−ISW ℓ , as the final CIB contamination model.
The Ty cross-correlation data are then modeled as where C tSZ−ISW ℓ is the theoretical prediction of the ISW and tSZ cross-power spectrum (Equation ( 18)).C CIBavg−ISW ℓ is the average of the three cross-power spectra of the ISW and CIB.We then introduce a dimensional parameter B CIB (with unit sr • Myr −1 ).This parameter estimates the CIB contamination level in the yT cross-correlation analyses.We remind the reader that B CIB is different from but related to the dimensionless A CIB .The latter represents the amplitude of the CIB contamination in the yy autocorrelation.The two parameters are related because they model the same foreground contamination in different maps.
We do not expect to have a significant contribution of galactic foregrounds, because, as indicated in Section 2.1, the ISW map is almost free of foreground contamination.Also, for our multipole regime, the NILC tSZ maps seem to be quite robust against different masking scenarios (see Figure 11 in Planck Collaboration et al. 2016b).In addition, a significant fraction of the ISW signal comes from LSS tracers, which are not correlated with the galactic foregrounds.

CONSTRAINTS ON PARAMETERS
The theoretical prediction modeled in Sec. 2 allows us to connect the astrophysical, cosmological, and foreground parameters with the observed autoand cross-power spectra.The cosmological parameters that the power spectra are sensitive to include the fractional matter density Ω m , Hubble parameter h ≡ H 0 /100 km s −1 Mpc −1 , and the rms matter density fluctuations at 8 h −1 Mpc, i.e., σ 8 .The astrophysical parameter involved is mainly the one related to W SZ , i.e. the product of b gas (gas bias), ne (the electron number density today), and T e (the electron temperature at present).This joint parameter quantifies the amplitude of the warm baryonic dark matter on the largest scales.The foreground parameter A CIB gauges the amplitude of the foreground spectrum to the measured C yy ℓ .B CIB captures the amount of the CIB contamination at the level of the ISW-tSZ cross-correlation.In the following subsection, we describe the details of our methodology to constrain these parameters 10 .

Likelihood method
Here we compare the theoretical power spectra with the observed spectra.Let the theoretical spectra be C th ℓ (Θ) and the observed spectra be C obs ℓ .The predictions are controlled by the parameter set Θ ≡ (Ω m , h, σ 8 , W SZ , B CIB , A CIB ), while the observed power spectra C obs ℓ ≡ (C TT ℓ , C Ty ℓ , C yy ℓ ).With the full covariance matrix Cov tot calculated in Sec.2.4, we formulate the χ 2 function as 10 Appendix B and Figure 14 show the comparison between results from using emcee and ultranest packages.
Table 2. Results of Parameter Estimation.Note.For each parameter, we report the range over which it was fitted with a flat prior (second column) and best-fit values with 68% CL (third column).All parameters are dimensionless, apart from BCIB in [sr MJy −1 ] (We have converted the dimensional constraint on W SZ to the joint constraint on bgasTe(0)ne).The first three rows on the upper partition are the constraints on the cosmological parameters.The fourth row is the derived constraint on S8.The middle partition is the constraint on the gas parameter and the lower one is for foreground parameters.Notice that for the h value, we quote the diagonal block's constraint here (Equation ( 33)) for the reason given in Section 5.The likelihood function L(Θ) ∼ e −χ 2 /2 and the final posterior probability distribution for the parameters P (Θ|d) are related by where P (Θ) is the prior probability function, for which we assume flatness for all parameters within the ranges listed in Table 2. L(d|Θ) is the likelihood function of the data, which includes both contributions from the cosmic variance and noise; it also considers the effect of masks on the data.We then employ the Markov Chain Monte Carlo (MCMC) technique11 to explore the parameter space.We then scale the constraints on σ 8 by Ω m to derive S 8 ≡ σ 8 (Ω m /0.3) 0.5 .The final posterior distributions on our parameter sets are plotted in Figure 8.

Constraints on LSS parameters
Starting with the fundamental parameters, we estimated the matter density parameter, Hubble constant, and the level at which matter clusters (the amplitude of  Additionally, we used the measured correlation functions to derive a constraint on the degenerate combination: S 8 ≡ σ 8 (Ω m /0.3) 0.5 = 0.755 ± 0.06 (68% CL). ( 26) Our constraint on Ω m broadly agrees with previous literature, as shown in Figure 9.However, we see that the matter density of the Universe, Ω m is slightly correlated with the tSZ gas parameter.This is indeed expected, as the amplitude of the tSZ effect (set by the gas parameter) is proportional to the line-of-sight integral of the electron pressure times the electron density along the line of sight.The electron pressure is related to the temperature of the gas, which is determined by the gravitational potential of the cluster and the energy input from various sources, such as accretion shocks and feedback from active galactic nuclei.On the other hand, the electron density is related to the matter density of the gas.Hence, the tSZ effect is a sensitive probe of the underlying matter density of the Universe and should be correlated.In addition, we see that the tSZ gas parameter is correlated with the Hubble parameter h but anticorrelated with σ 8 .Although the tSZ effect and the Hubble parameter are not directly related to each other, both can be correlated, because they are both sensitive to the properties of galaxy clusters.The tSZ effect in a sample of galaxy clusters can be used to estimate their masses, and distances, which in turn can be used to constrain the Hubble parameter.Hence, W SZ and h can be  correlated.However, there is no direct explanation for the anticorrelation between the tSZ gas parameter and the amplitude of the matter power spectrum on scales of 8 h −1 Mpc.This is because they are sensitive to different aspects of the LSS of the Universe.While σ 8 characterizes the strength of the clustering of matter on large scales in the Universe, W SZ relates to the properties of the gas in clusters, so an anticorrelation might be possible.To further prove whether they are both correlated would depend on several factors, including the redshift range and mass range of the specific galaxy clusters contributing to the tSZ effect, as well as the assumed cosmological model, which is beyond the scope of this current work.
We now discuss our constraints on the S 8 parameter in comparison with those obtained from other measurements.First, we would like to compare our estimate with Planck Collaboration et al. (2016b), considering the fact that we utilize the same tSZ data.We remind the reader that Planck Collaboration et al. (2016b) parameterized the amplitude of matter fluctuation as S 8 = σ 8 (Ω m /0.28) 3/8 to obtain a value of S 8 = 0.80 +0.01 −0.03 and S 8 = 0.90 +0.01 −0.03 for a mass bias of 0.2 and 0.4 respectively.Using the same parameterization, our results yield S 8 = σ 8 (Ω m /0.28) 3/8 = 0.768 ± 0.05, a value that is consistent in 0.6σ (for mass bias = 0.2), and 2.3σ (for mass bias = 0.4), respectively.In Figure 9, we further show the comparison of our results with measurements obtained from other analyses, such as the measurements of the E-mode polarization and temperature-E correlation of the CMB from the South Pole Telescope-3G (SPT-3G) 2018 data Dutcher et al. 2021), the Atacama Cosmology Telescope (ACT) Data Release 4 (DR4) data (Aiola et al. 2020), the HSC lensing and cross-correlation with tSZ (Osato et al. 2019), and the measurements from the Planck "TT, TE, EE+lowE" and " TT, TE, EE+lowE+lensing" data (Planck Collaboration et al. 2020).Our estimate is consistent with the KiDS lensing cross-correlation with tSZ (Tröster et al. 2021), the Dark Energy Survey (DES) cosmic shear, galaxy clustering, and galaxy-galaxy lensing (Y1 and Y3) results (Secco et al. 2021;Porredon et al. 2021), the HSC lensing cosmic shear (Hamana et al. 2020b), the HSC-Y1 power spectrum (Hamana et al. 2020a) and SPT-3G (Dutcher et al. 2021), within the estimated errors.However, we prefer a slightly lower value compared to the CMB results.
In contrast to the usual banana-shaped contour, which characterizes the constraints between the S 8 and Ω m degeneracy from the cosmic shear and shear-tSZ cross correlation (Abbott et al. 2018), the shape of our constraints is consistent with the Planck TT, TE, and EE+ lowE measurements (Planck Collaboration et al. 2020) within 1.1σ, as shown in Figure 10.Indeed, our analysis combined the tSZ with the ISW tracer which probes the late-time Universe.Therefore, such a small difference could arise from the contribution introduced by the late-Universe probe in our work.

Constraints on the gas parameter
The tSZ gas parameter is estimated as the joint constraint on the product of the mean electron density ne , electron temperature T e , and gas bias b gas at redshift z = 0. We parameterize the product as W SZ (via Equation ( 12)) and obtained a value of12 W SZ ≡ b gas T e 0.1 keV ne 1 m −3 = 3.09 +0.320  −0.380 .(27)   This result is consistent with the value estimated by Van Waerbeke et al. ( 2014) within 1.7σ.Let us remark that, as explained above, in our analysis we fixed the n s , τ , and A s parameters.This choice was taken because our parameter space could be too complicated for our likelihood.However, it is worth noting that, whereas the n s and A s uncertainties are very small compared to the ones derived for our parameters, τ is constrained at around 12% by Planck Collaboration et al. ( 2020), similar to the one obtained for W SZ .Taking into account that W SZ is proportional to n e and that this is also proportional to τ , we can derive a rough estimation of how much our uncertainty on W SZ could grow by if τ would be marginalized instead of fixed.In that case, our error will grow from the actual 12% to 18%.Notice that the temperature values derived from different measurements correspond to cosmic structures on different scales; hence, this plot is only intended to provide a comparison between their orders of magnitude, which are found to be in broad agreement.
We can further derive an estimate on the electron temperature for today by adopting the numerical values of b gas and ne from Refregier et al. (2000) and Seljak et al. (2001), as 4 ≤ b gas ≤ 9. Taking b gas = 6 and ne = 0.25 m −3 , we found which is consistent with the values obtained from de Graaff et al. ( 2019), in 0.18σ, with the expected WHIM temperature.Note that that analysis in the work of de Graaff et al. ( 2019) was entirely interpreted using halo modeling and hence was able to resolve individual halos, thereby providing information on filaments in Galaxy clusters up to small scales.Now we provide a comparison of our estimate on the electron temperature with other studies in Figure 11.We found our result in broad agreement with previous studies from numerical simulations and observational analyses.Hydrodynamic simulations (Cen & Ostriker 1999;Davé et al. 2001) showed that the diffuse baryons locked up in the intergalactic medium (IGM) are heated up to T ∼ 10 5 -10 7 K, which suggests that a substantial fraction of the "missing baryons" in the Universe are in the phase of the warm-hot IGM that traces dark matter.We now compare our results with other observational analyses.One should notice that different studies probe the LSS tracers with different densities and environments; hence, it is not very meaningful to directly compare the estimates of gas temperature from different measurements.Nonetheless, the cross-correlation analysis usually receives contributions from the twohalo term, which corresponds to the correlated gas on large scales.In contrast, stacking analyses are usually sensitive to the particular small-scale configurations of individual systems, which usually correspond to highdensity regions.Therefore, in Figure 11, we separate the two-point correlation measurements of the gas temperature from stacking analyses.
In the same category of two-point correlation constraints, Muñoz & Loeb (2018) cross-correlated FRBs with the tSZ map and estimated an IGM temperature of T e = (1.1±0.1)×10 6K.The cross-correlation between the tSZ maps from the Planck data and the lensing convergence map from the Canada-France-Hawaii Lensing Survey (Erben et al. 2013) showed that the diffuse baryons, which are responsible for the lensing-SZ cross-correlation, have temperatures 10 5 -10 7 K (Van Waerbeke et al. 2014;Ma et al. 2015).Génova-Santos et al. (2015) performed a follow-up analysis using the foreground-cleaned CMB maps in cross-correlation with the Two-Micron All Sky Redshift Survey (2MASS) of galaxies and obtained a fraction of the baryons at temperatures in the range 10 4.5 -10 7.5 K, consistent with Van Waerbeke et al. ( 2014) at z ∼ 0.5.The temperature values found from these studies are in broad agreement with our current study.
As for other analyses, Eckert et al. (2015) investigated the A2744 cluster over a physical length of 8 Mpc using X-ray observations and found out that the filamentary structures of gas are in a plasma temperature range of T = 1 × 10 7 K-2 × 10 7 K for the various filaments, after discarding and eliminating coincidental X-ray point sources.Similarly, Bulbul et al. (2016) did a thorough analysis of the A1750 cluster using Suzaku and Chandra X-ray observations and Multiple Mirror Telescope optical observations, to report that in both the filaments and off-filament directions, the gas mass is in agreement with the cosmic baryon fraction at R 200 .In addition, for stacking analyses, Tanimura et al. (2019) stacked the luminous red galaxy samples from the DR12 catalog of SDSS on Planck tSZ map to search for warm gas filamentary structure.They obtained a constraint on the mean electron-density-weighted temperature as ∼ (8.2 ± 0.6) × 10 6 K.With X-ray observations, Tanimura et al. (2022a) obtained an average temperature of (1.2 +0.35  −0.23 ) × 10 7 K. Alvarez et al. (2018) also obtained a measurement of the global projected filament temperature to an estimate of T e = (5.2+1.03 −0.64 )×10 7 K, using measurements from the Chandra and XMM-Newton X-ray observations of the A3391/A3395 intercluster filament, an estimate that is higher than other measurements, including our results, but compatible with Tanimura et al. (2022a) in the 3.8σ CL.These single-system stacking analyses used specific tracers to infer the WHIM properties.We stress that our result is somehow in between the stacking and crosscorrelation results, as shown in Figure 11.We think this fact deserves further studies with numerical simulation and the kinematic Sunyaev-Zeldovich (kSZ) effect (Sunyaev & Zeldovich 1980;Hernández-Monteagudo et al. 2015;Ma 2017;Schaan et al. 2021).In terms of breaking degeneracy in Equation ( 27), the kSZ effect is solely dependent on the electron density and peculiar velocity of the cluster, so if the latter can be inferred by another method (e.g.reconstructing the velocity field from the redshift-space density field; Hernández-Monteagudo et al. (2015); Planck Collaboration et al. (2016d)), the electron density can be measured.In this way, the de-generacy between the electron density and temperature can be broken and one can provide an independent measurement of the gas temperature (see, e.g.Schaan et al. 2021).Other relevant measurements that have probed WHIM properties include but are not limited to Wang (1993), Werner et al. (2008), Nevalainen et al. (2019) and Erciyes et al. (2023).
We now consider the constraints on the foreground parameters.We quantify the level to which the clustered CIB contributes to the observed auto-power spectrum for yy and find a value of We remark that the anticorrelation between A CIB and W SZ seen in Figure 8 is expected, because the tSZ effect and the CIB arise from different physical mechanisms.
The CIB is thought to be produced by star formation in galaxies, which is more likely to occur in regions of high gas density, while the tSZ effect occurs in regions of low gas density.Also, the redshift and mass dependences of the tSZ effect and the CIB are different.While the tSZ comes from more massive clusters, CIB is stronger for less massive objects.Therefore, an anticorrelation is possible between the two parameters.
Similarly, we stress that the clustered CIB could also contribute to the cross-correlation function, as reported in Makiya et al. (2018).Therefore, we parameterized such a contribution in this analysis with B CIB .This parameter gauges the contamination from the clustered CIB to the Ty cross-correlation.We obtained the bestfitting value in the order of ∼ 10 −4 sr MJy −1 and show in Figure 3 its contribution to the correlation of the Planck ISW with the Compton parameter map.This contamination is dominant on low -ℓ, high z (ℓ eff ≤ 30).

Hubble constant measurement
The Hubble constant is a fundamental cosmological parameter that quantifies the expansion rate of the Universe.In recent years, the tension on the value of H 0 has risen between early (global) and late (local) Universe measurements (see Figure 12 for the comparison).The early Universe measurement from the CMB radiation prefers a somewhat lower value, as H 0 = 67.4± 0.5 km s −1 Mpc −1 (Planck Collaboration et al. 2020).The DES clustering, baryon acoustic oscillation (BAO), and Big Bang Nucleosynthesis (BBN) measurements also give a lower value H 0 = 67.4± 1.2 km s −1 Mpc −1 (Abbott et al. 2018).In contrast, late-time observables give a higher value than these.The standard distance ladder from Cepheid variables gives H 0 = 74.0 ± 1.4 km s −1 Mpc −1 (Riess et al. 2019), and the strong-lensing time-delay distances (H0LiCOW team) give H 0 = 73.3± 1.8 km s −1 Mpc −1 (Wong et al. 2020).In addition, a combination of supernova type Ia (SN), cosmic chronometers (CC), and BAO data with the time-delay cosmography from H0LiCOW were used to obtain a value of H 0 = 73.8± 0.8 km s −1 Mpc −1 (Bonilla et al. 2021).However, the H 0 value obtained by the Carnegie-Chicago Hubble Program (CCHP), based on the Tip of the Red Giant Branch (TRGB) measurements in the Large Magellanic Clouds, falls between the early-and late-time measurement, yielding H 0 = 68.9±1.9 km s −1 Mpc −1 (Freedman et al. 2019).
In this work, to study the robustness of our Hubble constant measurement, we carefully examine the H 0 value using individual likelihoods and joint probes.The individual likelihood just uses the specific covariance blocks in Figure 4, calculated via Equation ( 5).By using tSZ data only (the middle block in Figure 4 for its covariance) and ISW data only (the lower left block), we obtain where for tSZ we marginalized over the (Ω m , σ 8 , W SZ , B CIB , A CIB ) parameters, and for ISW we marginalized only the (Ω m , σ 8 ) parameters.One can see that the tSZ-only data make the H 0 value consistent with the Planck result within the error range, whereas the ISW data make the H 0 constraint closer to the "SH0ES" local measurements (Riess et al. 2019), "SN+CC+BAO+H0LiCOW" from Bonilla et al. (2021), and H0LiCOW measurement (Wong et al. 2020;see Figure 12).By using the tSZ-ISW cross-correlation likelihood (using the Cov Ty,Ty covariance; i.e., the upper right block of Figure 4), we obtained which error bar overlaps with both single-likelihood cases and is consistent with the tSZ-only case at the ≲ 0.5σ CL.We now combine the datasets and constrain H 0 .We first used the full cross-correlation data by adding crosscovariance terms in Equation (5), i.e. using all the blocks in Figure 4, and obtained This result is in excellent agreement with the values obtained by Planck (Planck Collaboration et al. 2020) and that obtained by the DES clustering, BAO, and BBN measurement (Abbott et al. 2018) in the ≲ 0.5σ CL.However, this joint constraint sits outside the range of individual constraints, and deviates from the ISW-only constraint by more than 1σ.As a sanity check, we added the offdiagonal cross-covariance matrix one after the other to the diagonal covariance and found that the driving force of this low-H 0 value is the correlated cosmic variance from the TT to Ty statistics.In other words, all the different combinations except the crosscovariance Cov TT,Ty provide consistent results with the two individual constraints.Therefore, instead of combining all covariance blocks, if we only use the diagonal blocks of the covariance matrices (the TT, yy, and Ty autocorrelations), we obtain which sits in between the ISW-only and tSZ-only constraints and is close to the tSZ and ISW crosscorrelation-only constraint (the green data point in Figure 12).This result is also in excellent agreement with the CCHP measurement (< 0.1σ;  33), (30), and (31) (also in Figure 12), we conclude that the low value of the joint probes may be due to some unaccounted systematics in the TT-Ty covariance block.Therefore, instead of quoting the "all-combined" constraint on H 0 , we quote the joint constraint of the H 0 value obtained by using the diagonal blocks of the covariance matrix, i.e.Equation (33).For completeness of different scenarios, we also show the H 0 value obtained using the full covariance case and the individual constraints in Figure 12.We also notice that recently Capozziello et al. (2023) reported that the Hubble tension can be removed in the framework of ΛCDM in terms of the cosmological lookback time (the redshift at which the measurements are performed).No definitive conclusion has been reached about the value of H 0 , which provides a strong reason to suspect the existence of physics beyond the ΛCDM model (Dai et al. 2020;Perivolaropoulos & Skara 2022).Future CMB experiment such as the Simon Observatory (Ade et al. 2019) and CMB-S4 (Abazajian et al. 2016), combined with the new generation of galaxy surveys, such as Euclid and LSST, are expected to reach a precision of ∼ 0.15% in the H 0 estimate (Di Valentino et al. 2021), thus providing a stringent test for the adopted cosmological models.

CONCLUSIONS
In this paper, we have provided a novel constraint on the gas parameter on large scales using the crosscorrelation of the ISW with the tSZ effect.There have been several theoretical studies on the benefits of cross-correlating the ISW with the tSZ effect (Cooray 2002;Taburet et al. 2011;Creque-Sarbinowski et al. 2016), but the actual measurement with the Planck data has not been done before.In this work, we measured the tSZ-ISW cross-correlation power spectrum by using the Planck ISW data and the tSZ Compton-y map from the Planck NILC algorithm.We then use the cross-correlation data to constrain fundamental cosmological parameters (Ω m , h, σ 8 ) and the gas properties on very large scales, via the parameter W SZ = b gas (T e /0.1 keV) ne /1 m −3 .Our estimates on Ω m and σ 8 enabled us to estimate the amplitude of matter density fluctuations S 8 ≡ σ 8 (Ω m /0.5) 0.5 .We also calculated the theoretical power spectrum based on large-scale gas bias models to account for the measured spectra (Van Waerbeke et al. 2014).We used MCMC methods to fit the model parameters and derived large-scale constraints on the cosmological and gas parameters.
For the measurement part, we employed the MAS-TER (pseudo-C ℓ ) approach encoded in the NaMaster code to measure the auto-and cross-correlation angular power spectra.Our measured ISW and tSZ autocorrelation power spectra (C TT ℓ and C yy ℓ ) are consistent with previous measurements (ISW -see Figure 7 of Rahman & Jawed Iqbal 2016; tSZ-see the results in Makiya et al. 2018;Bolliet et al. 2018;Osato et al. 2019;Tanimura et al. 2022b;Ibitoye et al. 2022).The important result is the first detection of the ISW-tSZ cross-correlation power spectrum.Here we simulated 900 ISW maps to quantify the covariance of correlation with tSZ data.With the 100 simulated ISW maps being cross-correlated with true tSZ data, we quantify that the true ISW-tSZ cross-correlation is detected at the 3.6σ CL., with its amplitude and shape compatible with the theoretical prediction in Creque-Sarbinowski et al. ( 2016) (see, e.g., their Figure 1).To maximize the statistical information in our measurement we combined the power spectra into a vector and then fit for the cosmological and astrophysical parameters.We explored the parameter space with an MCMC method, using the Python emcee package (Foreman-Mackey et al. 2013) to determine the maximum likelihood of our data, and obtained posterior probability distributions on all parameters.For a sanity check, we also used another independent package, ultranest (Buchner 2016(Buchner , 2021)), to explore the parameter space and found its results compatible with the ones obtained using emcee package (see Figure 14 for such a comparison).We then adopted the GetDist Python software detailed in Lewis (2019) to plot the final distributions in Figure 8.
We show the constraints obtained on Ω m , and σ 8 in Figure 9. Additionally, we compared different H 0 measurements from the literature and summarized them in Sec.5.4 and Figure 12.Our cross-correlation-data estimated H 0 value is in between the CMB measured value and the local distance ladder measurement value (SH0ES) and is as such more consistent with the result from the TRGB measurement (Freedman et al. 2019).On the other hand, our diagonal constraints yield a value of H 0 = 69.7 +2.0 −1.5 km s −1 Mpc −1 , a value that sits between the global and local measurements of the Hubble constant, while the joint analysis prefers a value of H 0 = 66.5 +2.2 −1.9 km s −1 Mpc −1 , a value that is consistent with Planck.
We obtained a joint fit of the gas bias, electron temperature, and electron number density as W SZ = b gas (T e /0.1 keV) ne /1 m −3 = 3.09 +0.320  −0.380 .We used the model amplitude, and an assumed value of b gas and ne , to estimate the electron temperature as T e = (2.40 +0.250  −0.300 ) × 10 6 K. Our cross-correlation signal would then account for all the missing baryons that lie within the temperature range 10 5 -10 7 K if we assume that they were all located in the WHIM component.
Last, we also estimated the CIB contamination to the cross-correlation function as B CIB , while A CIB dictates the foreground amplitude of the CIB contamination in the Compton parameter map.The results are B CIB × 10 4 = 4.0 +1.200  −1.100 and A CIB = 0.81 +0.086 −0.079 .We also derived a constraint on the amplitude of the matter fluctuations on an 8h −1 Mpc scale, S 8 = σ 8 (Ω m /0.3) 1/2 = 0.755 ± 0.060.Our result is in broad agreement with the values obtained in the KiDS measurement and our error estimate with that obtained in the DES cosmic shear Y1 and Y3 results.In general, we show in Figure 9 a comparison of our result with different data and methods (Troxel et al. 2018;Osato et al. 2019;Hamana et al. 2020a,b;Heymans et al. 2021;Porredon et al. 2021;Secco et al. 2021;Tröster et al. 2021).
Upcoming surveys like LSST (?), Euclid (Amendola et al. 2018), and CMB Stage-4 experiments (Abazajian et al. 2016), with increased redshift depth and sky coverage, will provide a better measurement of the ISW than current observational data.The better angular resolution of CMB-S4 and the higher-precision measurement of the ISW will enable the extension of this correlation up to high-ℓ modes, which have the promise to reach a much higher precision of the cross-correlation measurement (Cooray 2002 forecasted a ∼ 60σ detection).In parallel, high-resolution measurements would require a full halo model to account for the diffuse baryon compo-nents in halos with different masses.Combining sophisticated modeling with accurate data will constrain astrophysical and cosmological parameters to higher precision.and Cosmology with Theoretical Models Confronting Observational Data" of the National Institute for Theoretical and Computational Sciences of South Africa.A.I. acknowledges the support of the National SKA Program of China with grant No. 2022SKA0110100, and the Alliance of International Science Organizations, grant No. ANSO-VF-2022-01.Y.Z.M. acknowledges the support of the National Research Foundation with grant No. 150580.P.V. thanks the Spanish Agencia Estatal de Investigación (AEI, MICIU) for the financial support provided under the projects with references PID2019-110610RB-C21, ESP2017-83921-C2-1-R, and AYA2017-90675-REDC, co-funded with EU FEDER funds, and acknowledges support from Universidad de Cantabria and Consejería de Universidades, Igualdad, Cultura y Deporte del Gobierno de Cantabria, via the "Instrumentación y ciencia de datos para sondear la naturaleza del universo" project, as well as from Unidad de Excelencia María de Maeztu (MDM-2017-0765).D.T. acknowledges financial support from the XJTLU Research Development Fund (RDF) grant with the number RDF-22-02-068.W.M.D. acknowledges the support from the "Big Data for Science and Society" UKZN Research Flagship.X.L. acknowledges the support of the Ministry of Science, and Technology (MoST) inter-government cooperation program China-South Africa Cooperation Flagship Project 2018YFE0120800.We thank R. Belén Barreiro for providing us with the Planck ISW simulations.

SOFTWARE AND DATA
For the analysis presented in this manuscript, we made use of the following software: HEALPix (Górski et al. 2005), MASTER (Hivon et al. 2002), NaMaster (Alonso et al. 2019), emcee (Foreman-Mackey et al. 2013), ultranest (Buchner 2016(Buchner , 2021)), and GetDist (Lewis 2019).The data underlying this article are publicly available.The Planck Compton parameter map, the Integrated Sachs-Wolfe map, the Galactic, and the pointsource masks are available at the Planck Legacy Archive at http://pla.esac.esa.int/pla.The cross-correlation  12)) using the yy-auto power spectrum and the joint analysis.
products and their covariance matrices can be shared on request to the corresponding author.APPENDIX A. GAS PARAMETER The gas parameter, W SZ , which is defined as a product of the mean electron density ne , electron temperature T e , and gas bias b gas , at redshift z = 0, is a quantity that directly sets the amplitude of the yy-auto power spectrum.Figure 13 shows the marginalized likelihoods of W SZ by using only the C yy ℓ -only and the joint spectra (C Tot ℓ , i.e.Equation ( 5)).One can see that W SZ = 3.29 +0.28 −0.24 for C yy ℓ -only and W SZ = 3.09 +0.32 −0.38 for the joint spectra.These two estimates are consistent within 0.5σ and the difference in the constrained value may not be due to the fact that a higher value of the gas parameter is needed to reproduce the yy-auto power spectrum, but that the tSZ is more sensitive to the abundance of cosmic gas (baryons) on large scales.

B. CONSISTENCY CHECK ON PARAMETERS
To ensure the accuracy of our parameter estimation, we employ both the Python emcee package (Foreman-Mackey et al. 2013) package and the ultranest package (Buchner 2016(Buchner , 2021) ) to explore the parameter space.We used the two packages independently to determine the minimum chi-square or the maximum likelihood for our data and obtained posterior probability distributions on all parameters.We then used the GetDist Python software detailed in Lewis (2019) to plot the final distributions in Figure 14.One can see from Figure 14 that both packages yield close results in the parameter distributions, therefore we adopt the Python emcee package in our concluding section.

C. COMPARING YY POWER SPECTRUM ANALYSES
Here we compare and contrast between the measurements and the modeling presented in our analysis and those presented in Planck Collaboration et al. (2016b).The yy angular power spectrum and the associated error bars measured by Planck were computed using XSPECT 13 (Tristram et al. 2005).On the other hand, we extracted the 13 XSPECT is a code that was initially written to measure the CMB temperature angular power spectrum from Archeops data.Given a set of maps, it could also measure the cross-correlation angular power spectrum and the associated error bars computed analytically from the data with no MC simulations involved.It utilizes the standard MASTER-like approach (Hivon et al. 2002) to correct for the beam convolution and the pixelization, as well as the mode-coupling induced by masking sky regions that are foreground-contaminated.

Figure 2 .
Figure 2. All-sky Compton parameter map (y-map) measured by Planck based on the NILC algorithm.The pointsource mask and the 40% Galactic-plane mask have been applied.

Figure 3 .
Figure 3.The measurements of the power spectra, C yy ℓ (y-auto, upper-left), C TT ℓ (ISW-auto, upper-right) and C yT ℓ (ISW-y cross, lower).Black dots and bars in each case are the means of the spectrum measured from the associated map and the error bars are estimated as the square root of the diagonal values of the corresponding covariance matrix.Other curves have been defined in the legend of each plot.

Figure 4 .
Figure 4.The six independent blocks of the full correlation matrix defined in Equation (7).The diagonal blocks show the standard correlations for the TT, yy, and Ty spectra; the offdiagonal blocks show the additional cross-correlations between different spectra.The correlation values are shown for pairs of the effective band power multipoles defined in Table1.
Figure3are the binned-ℓ ones, to reduce the empirical scatter between adjacent multipole amplitudes ℓ.We adopt the same scheme for binning as in PlanckCollaboration et al. (2016c).In addition, this makes the power spectrum coupling matrix invertible to some extent.One may notice that the Planck Collaboration et al. (2016c) used 19 multipole bins in total, but we used 10 bins for the range of ℓ eff = 18.0 to ℓ eff = 181.0 in total.The binning is done by evaluating each bin's average power spectrum values.The resulting values are quoted in Table1.We remind the reader that the y-auto power spectrum presented in the upper left of Figure3was measured in full sky while the one presented on the left of Figure 11 in Planck Collaboration et al. (2016b) was measured in partial sky.As a comparison, we show in Figure 15 in Appendix C the C yy ℓ both from our analysis and Planck's partial-sky result.A closer look at Figure

Figure 5 .
Figure 5. Upper Left-the measurements of the ISW auto-power spectrum using the true Planck ISW map (black solid line) and the 100 simulated ISWs map (colored lines).Upper Right-the respective S/Ns of the 100 simulated ISW maps (red dots; with mean 2.90 and rms 0.24) and the true ISW map (black dashed line; 2.32).Lower Left-the cross-correlation power spectra of the 100 simulated ISWs map with the Planck tSZ map (colored lines) and the true Planck ISW map cross-correlation.Lower Right-the significance of the true ISW-tSZ cross-correlation (black dashed line; 5.63) and the 100 simulated ones (blue dots; with mean 3.12 and rms 0.70.).

Figure 6 .
Figure 6.The effective redshift range driving the autoand cross-power spectra, shown for three reference effective multipoles at ℓ= 10, 50, and 100.

Figure 7 .
Figure 7. Cross-power spectrum of the CIB in three Planck frequencies with the ISW.The spectrum was measured using the public code NaMaster.

Figure 8 .
Figure 8. Posterior distributions from MCMC with the full covariance matrix for the six free parameters in our model.The figure shows the joint constraints for all parameter pairs and the marginalized distributions for each parameter along the table diagonal.W SZ is a dimensionless quantity defined in Equation (12).

Figure 10 .
Figure 10.Base ΛCDM model 68% and 95% posterior distribution constraint contours on the matter density parameter Ωm and on the structure amplitude parameter S8.

Figure 11 .
Figure 11.Comparison of the estimate on baryon temperature from two-point correlation studies (FRBs × Planck tSZ; Muñoz & Loeb 2018), stacking analyses (Planck tSZ-LRGs DR12-Tanimura et al. 2019; Planck tSZ-CMASS-de Graaff et al. 2019; Planck tSZ-Planck lensing-Tanimura et al. 2020), and X-ray observations (X-Ray ROSAT; Wang 1993).Notice that the temperature values derived from different measurements correspond to cosmic structures on different scales; hence, this plot is only intended to provide a comparison between their orders of magnitude, which are found to be in broad agreement.

Figure 13 .
Figure13.Comparison of the estimates on the gas parameter W SZ (Equation (12)) using the yy-auto power spectrum and the joint analysis.

Figure 14 .Figure 15 .
Figure14.Posterior distributions from MCMC using both emcee and ultranest with the full covariance matrix for the six free parameters in our model.The figures show the joint constraints for all parameter pairs and the marginalized distributions for each parameter along the table diagonal.Both results are compatible within uncertainties, and they have no obvious discrepancies.Therefore, for this analysis, we present the results obtained using the emcee package, a more versatile and widely employed MCMC engine.

Table 1 .
Summary of Power Spectra Measurements Averaged over [ℓmin, ℓmax]Multipole Intervals.Note.For each effective multipole ℓ eff , we report the values of the power spectra labeled as C yy ℓ , C TT ℓ , C Ty ℓ .We also report the Gaussian error contribution estimated using the Namaster-based covariance, the non-Gaussian error contribution estimated using the pyccl-based covariance and the simulation-based covariance described in Sections 2.4 and 3.1, labeled as σG(C yy ℓ ), σNG(C yy ℓ ), σTT, and σTy, respectively.