Cosmology from weak lensing, galaxy clustering, CMB lensing and tSZ: I. 10 × 2pt Modelling Methodology

The overlap of galaxy surveys and CMB experiments presents an ideal opportunity for joint cosmological dataset analyses. In this paper we develop a halo-model-based method for the first joint analysis combining these two experiments using 10 correlated two-point functions (10 × 2pt) derived from galaxy position, galaxy shear, CMB lensing convergence, and Compton-y fields. We explore this method using the Vera Rubin Observatory Legacy Survey of Space and Time (LSST) and the Simons Observatory (SO) as examples. We find such LSS × CMB joint analyses lead to significant improvement in Figure-of-Merit of Ω m and S 8 over the constraints from using LSS-only probes within Λ CDM. We identify that the shear-y and y - y correlations are the most valuable additions when tSZ is included. We further identify the dominant sources of halo model uncertainties in the small-scale modelling, and investigate the impact of halo self-calibration due to the inclusion of small-scale tSZ information.


INTRODUCTION
Observations of the Universe's large-scale structure (LSS) provide valuable insights into cosmic structure formation and expansion history, enabling tests of theories of gravity and constraints on the mass and number of neutrino species, the nature of dark matter, and dark energy.To address these science questions, several galaxy survey experiments have been developed, including the Kilo-Degree Survey (KiDS, Hildebrandt et al. 2017;Heymans et al. 2021), the Dark Energy Survey (DES, Abbott et al. 2018Abbott et al. , 2022)), the Hyper Suprime-Cam (HSC, Hikage et al. 2019), the Dark Energy Spectroscopic Instrument (DESI, DESI Collaboration et al. 2016), the Vera C. Rubin Observatory Legacy Survey of Space and Time (LSST, Ivezić et al. 2019), the Nancy Grace Roman Space Telescope (Akeson et al. 2019), the Euclid mission (Laureijs et al. 2011), and the Spectro-Photometer for the History of the Universe, Epoch of Reionization, and Ices Explorer (SPHEREx, Doré et al. 2014).Jointly analysing multiple cosmological probes, such as "3×2pt" analyses that combine galaxy clustering and weak lensing statistics, has become the standard in KiDS, DES, and HSC since it increases the overall constraining power and the robustness to systematic effects.This "multiprobe analysis" approach is expected to remain a key strategy in extracting cosmological information from upcoming galaxy survey experiments.
A series of experiments have been dedicated to measuring the ⋆ E-mail: xfang@berkeley.eduenergy composition and structure growth of the Universe at much higher redshifts via the cosmic microwave background (CMB).Whether there exists an S 8 parameter difference within the standard ΛCDM model between high-z CMB measurements from Planck (Planck Collaboration et al. 2020) and several low-z measurements, such as KiDS (Joudaki et al. 2018;Heymans et al. 2021), DES (Abbott et al. 2018(Abbott et al. , 2019(Abbott et al. , 2022)), unWISE and Planck CMB lensing tomography (Krolewski et al. 2020(Krolewski et al. , 2021)), and Atacama Cosmology Telescope (ACT) CMB lensing (Qu et al. 2023;Madhavacheril et al. 2023), has been extensively investigated, yet the issue remains unresolved.The significance of the tension/agreement between these measurements, and potential deviations from the ΛCDM model, will become more clear with future CMB analyses from ACT, (Aiola et al. 2020;Choi et al. 2020), the South Pole Telescope (SPT, Benson et al. 2014;Dutcher et al. 2021;Balkenhol et al. 2022), the Simons Observatory (SO, Galitzki et al. 2018;Ade et al. 2019), and the CMB Stage-4 survey (S4, Abazajian et al. 2016Abazajian et al. , 2019)), along with wider and deeper galaxy surveys.
This paper explores the synergies of a CMB×LSS joint analysis using the Simons Observatory (SO) and Vera C. Rubin Observatory's Legacy Survey of Space and Time (LSST) as examples.We perform the first "10×2pt" simulated analysis, combining two-point functions between fields of galaxy density, shear, CMB lensing convergence, and tSZ Compton-y, which extends our previous work on 6×2pt analyses combining galaxy density, shear, CMB lensing convergence (Fang et al. 2022). 1 We present our theoretical modelling of the combined "10×2pt" probes in Sect.2, including the analysis choices, data vector, and analytic covariance modelling.We carry out the simulated likelihood forecast analysis for LSST + SO in Sect.3. We compare the cosmological constraining power of 10×2pt with 3×2pt and 6×2pt and identify the subset of probes containing most cosmological information.We also identify a subset of halo parameters responsible for improving the self-calibration between probes related to matter distribution (i.e., lensing and galaxy clustering) and probes related to gas distribution (i.e., tSZ).We conclude in Sect. 4.
In a companion paper (Fang et al. in prep, Paper II), we further investigate the 10×2pt synergies between various scenarios of the Roman Space Telescope galaxy survey and SO/S4, including the benefits of choosing an alternative galaxy sample and the impact of halo model uncertainties.

MULTI-SURVEY MULTI-PROBE ANALYSIS FRAMEWORK
This section describes the ingredients for simulated 10×2pt analyses with focus on the tSZ auto-and cross-correlation modelling, building on the 3×2pt and 6×2pt analyses modelling and inference described in Krause & Eifler (2017) and Fang et al. (2022).

Survey Specifications
Simons Observatory The Simons Observatory (SO, The nominal design of the observatory will include one 6 m largeaperture telescope (LAT, Xu et al. 2021;Parshley et al. 2018) and three small-aperture 0.5 m telescopes (SATs, Ali et al. 2020).The LAT will produce temperature and polarisation maps of the CMB with ∼arcmin resolution over 40% of the sky, with a sensitivity of ∼6 µK•arcmin when combining 90 and 150 GHz bands (Xu et al. 2021).These wide-field maps will be the key input to measure CMB lensing with SO.We also note that a large investment was recently announced that will significantly increase the detector count of the LAT and double the number of SATs, an upgrade known as "Advanced Simons Observatory" (ASO).Since little technical information is publicly available about ASO at present, our forecasts adopt the nominal configuration, and are therefore conservative.Further gains can be expected with ASO.We present the details of the noise levels and component separation techniques used in our analysis in App.B.
Vera C. Rubin Observatory's Legacy Survey of Space and Time will rapidly and repeatedly cover a ∼ 18,000 deg 2 footprint in six optical bands (320 nm-1050 nm).With a single exposure depth of 24.7 r-band magnitude (5σ point source), the 10 years of operations will achieve an overall depth of 27.5 r-band magnitude.
The performance of the LSST dark energy analysis given specific analysis choices is explored in the LSST-Dark Energy Science Collaboration's (LSST-DESC) Science Requirements Document (DESC-SRD, The LSST Dark Energy Science Collaboration et al. 2018).
For galaxy clustering and weak lensing analyses, we adopt the galaxy samples and the survey parameters from Fang et al. (2022), which is largely based on the DESC-SRD.We assume that LSST Y6 will cover the same area as the final SO Y5, i.e., 40% of the sky or Ω s = 16, 500 deg 2 .We also assume an i-band depth i depth = 26.1 mag for the weak lensing analysis, and i-band limiting magnitude i lim = i depth − 1 = 25.1 mag for the clustering analysis.
We parameterise the photometric redshift distributions of both the lens and source samples as normalised by the effective number density nx .N x is the number counts of lens/source galaxies, z ph is the photometric redshift, Ω is the solid angle.The parameter values of nx , z 0 , α are given in Tab. 1, x after redshift cuts.A fiducial Gaussian photo-z error has been convolved with each tomographic bin as described by Eq. (A14).
which are updated values from the LSST-DESC's Observing Strategy2 (Lochner et al. 2018).These values, slightly different from the DESC-SRD, are computed for the more recent optimisation studies of LSST observing strategies (Lochner et al. 2018).Following the redshift cuts in DESC-SRD (z min = 0.2, z max = 1.2 for the lens sample, z min = 0.2 for the source sample), we further divide each galaxy sample into N tomo = 10 equally populated tomographic bins as shown in Fig. 1.

Multi-probe modelling
We extend the existing 6×2pt model of the CosmoLike3 modelling and inference framework (Eifler et al. 2014;Krause & Eifler 2017) to include tSZ (cross-) correlations.To obtain an (internally) consistent model of all 10×2pt statistics, we adopt a halo model (Seljak 2000;Ma & Fry 2000;Peacock & Smith 2000;Cooray & Sheth 2002;Hill & Pajer 2013;Hill & Spergel 2014, c.f. App.A1.1) approach and specify all systematics at the level of the individual observables, i.e., the density contrast of lens galaxies δ g , the lensing convergence of source galaxies κ g , the CMB lensing convergence κ CMB , and the Compton-y field y.We summarise the well-established computation of angular power spectra and our model ingredients for δ g , κ g and κ CMB , including astrophysical systematics linear galaxy bias model, intrinsic alignments (IA) using the "nonlinear linear alignment" model (Hirata & Seljak 2004;Bridle & King 2007), and observational systematics, in App. A. Combined with tSZ and baryonic

Modelling baryons and tSZ cross-correlations
Our implementation of tSZ cross-correlations and baryonic feedback effects on the matter distribution follows closely the empirical HMx halo model parameterisation for matter and pressure (Mead et al. 2020).The parameters of the HMx model are calibrated so as to ensure agreement at the power spectrum level between the prediction of the HMx prescription and the BAHAMAS simulations (Mc-Carthy et al. 2017), which are a set of smooth-particle hydrodynamical simulations with separate dark matter and gas particles, including sub-grid recipes describing effects due to stars, supernovae, AGNs, and gas heating and cooling.As the dominant contribution to pressure (cross-) power spectra comes from massive clusters, an accurate calibration of the pressure power spectrum would require very large simulation volumes.The accuracy of HMx prescription for the matter-pressure cross-spectrum is of order ∼ 15% and worse for the pressure auto-power spectrum (Mead et al. 2020).While Stage-IV data analyses will certainly require refinements of the model parameterisation and parameter calibration, we use the HMx parameterisa-tion only as a toy model to qualitatively explore information content in the non-linear regime as well as the impact of parameter selfcalibration on cosmology constraints.The halo model implementation in CosmoLike uses the Navarro et al. (1997, NFW) profile and the Duffy et al. (2008) halo massconcentration relation to model the halo density profile ũδ (k, M, z), the Tinker et al. (2010) fitting functions to compute the linear halo bias b h,1 (M) and the halo mass function dn(M)/dM.
We use the ∆ ρ = 200 halo definition, so that the halo radius r 200 satisfies M = 4πr 3 200 ρ∆ ρ/3 and the halo scale radius is given by r s = r 200 /c(M, z).We note that several of these halo model ingredients differ from those assumed by Mead et al. (2020).However, we do not expect these differences to qualitatively change the parameter degeneracy structure relevant for self-calibration studies presented here.Our main goal is indeed to investigate, within a coherent modelling choice, the interplay between the different model parameters and estimate how much the inter-calibration between datasets improves the constraints of the parameters of interest.We expect that, even with its limitations, our model is representative enough to allow this study.
Electron pressure The Fourier-transformed y halo profile, ũy (k, M, z), which is the building block for halo model calculations of tSZ (cross-) power spectra (c.f.App.A1.1) is then computed from the electron pressure profile P e (M, r, z), (2) The HMx halo model parameterisation describes the gas distribution with two components, bound and ejected.The bound gas is modelled with the Komatsu-Seljak (KS) profile (Komatsu & Seljak 2002), while the ejected gas is assumed to follow the linear perturbations of the matter field.Therefore, the ejected gas only contributes to the 2-halo term.
The fractions of bound and ejected gas in halos are parameterised as where M 0 and β are an empirical parameterisation of the dependence of the fraction of ejected gas on halo mass, and the stellar fraction f * (M) is assumed to take a log-normal form, parameterised by A * , M * , σ * ).Ω b and Ω m are the cosmic baryon and matter density parameters.
The bound gas density profile ρ bnd can then be written as where the radial dependence follows the KS profile with polytropic index Γ (Komatsu & Seljak 2002).The normalisation is given In the KS profile the bound gas temperature T g is determined by hydrostatic equilibrium, where m p is the proton mass, µ p = 4/(3 + 5 f H ), with f H the hydrogen mass fraction, and α is a free parameter that encapsulates the deviations from the hydrostatic equilibrium due to non-thermal components of the gas.The electron pressure profile of the bound gas can then be written as where µ e = 2/(1+ f H ), and k B is the Boltzmann constant.The density profile for the ejected gas ρ ejc in 2-halo term is approximated as a 3D Dirac delta function distribution with a constant temperature T w and total mass M ejc , following the footnote 7 of Mead et al. (2020), Impact of baryonic feedback on the matter distribution To account for the impact of baryonic feedback on the dark matter distribution, we follow the HMx approach which modifies the dark matter-only halo mass-concentration (c D (M, z)) where ϵ 1 , ϵ 2 are free parameters and ϵ 1 = ϵ 2 = 0 reduces to the unmodified case.
Stellar component Throughout the analyses presented here, the HMx halo model parameters A * , M * , σ * , which describe the stellar mass fraction, are held fixed as the 10×2pt data vector has almost no sensitivity to these parameters.

Simulated likelihood analysis methodology
To simulate joint analyses of multi-survey multi-probe data, we perform simulated parameter inference assuming a Gaussian likelihood for the data d given a point p in cosmological and nuisance parameter space, where m is the model vector and C is the data covariance matrix.The synthetic (noiseless) data is calculated from the input model, with cosmic variance and shape/shot noise entering only through the covariance; this setup bypasses scatter in the inferred parameters due to realisation noise and allows us to focus on the constraining power of different analysis choices.The analytic Fourier space 10×2pt covariance calculation is summarised in App.B and we make available an implementation CosmoCov_Fourier4 , which is an extension of CosmoCov5 (Fang et al. 2020a).

Scale cuts
We compute all two-point functions in 25 logarithmically spaced Fourier mode bins ranging from ℓ min = 20 to ℓ max = 8000, and impose the following set of ℓ-cuts: • Galaxy clustering δ g scale cuts are driven by the modelling inaccuracy of non-linear galaxy biasing.Following the DESC-SRD, we adopt ℓ i max = k max χ(z i ), where k max = 0.3 h/Mpc and zi is the mean redshift of the lens bin i.
• Weak lensing κ g scale cuts are driven by model misspecifications for the impact of baryonic processes on the non-linear matter power spectrum as well as gravitational non-linearity.
We assume ℓ max = 3000 for all tomography bins, following the DESC-SRD.We also consider a more optimistic case in which we can reliably model shear up to ℓ max = 8000.We note that baryonic effects on the shear power spectrum will be significant for both of these scale cuts, and Stage-IV analyses will likely require modelling beyond the modified halo mass-concentration relation (Eq.11) implemented here to reach these scale cuts, these different ℓ max choices will illustrate the interplay of scale cuts and parameter degeneracies.
• CMB lensing κ CMB measurements are challenging at high ℓ due to lensing reconstruction challenges from foregrounds; similarly, various systematic errors can limit the reconstruction of lowℓ lensing modes (Darwish et al. 2021).We adopt the scale cuts ℓ min = 100, ℓ max = 3000.6More details on the lensing reconstruction are presented in App.B.
• tSZ y measurements are limited by atmospheric noise at lowℓ (in the case of SO) and component separation at high-ℓ (Ade et al. 2019).We adopt the somewhat optimistic scale cuts ℓ min = 80, ℓ max = 8000, noting in particular that contamination from the Cosmic Infrared Background (CIB) will need to be modelled and (likely) marginalised over.
For cross-power spectra, we adopt the more restrictive scale cut combination of the two fields.We further exclude δ g -κ g combinations without cosmological signal, i.e., with the lens tomography bin at higher redshift than the source tomography bin.

Analytical Marginalisation of photo-z and Shear Calibration Parameters
To speed up our forecast analyses, we analytically marginalise the 22 photo-z parameters and 10 shear calibration parameters, assuming that constraints on these parameters are dominated by their (Gaussian) priors listed in Tab. 2. This analytic marginalisation was derived in the context of cosmological surveys in Bridle et al. (2002), which we summarise below (see also Petersen & Pedersen 2012).Suppose that the likelihood L of obtaining the N-dimensional data vector d given the model vector m(p, λ λ λ) is a multivariate Gaussian with covariance C, i.e., d|m(p, λ λ λ) ∼ N(m(p, λ λ λ), C), where p, λ λ λ are parameters of the model, and λ λ λ contains the parameters we want to marginalise over, with a joint prior distribution p(λ λ λ).Marginalising the posterior distribution over λ λ λ has a simple analytic solution if we assume (1) the model is linear, and (2) the prior p(λ λ λ) is a multivariate Gaussian distribution.
The marginalised posterior can be derived as d|m(p, λ λ λ 0 ) ∼ N(m(p, λ λ λ 0 ), C mod ), where Therefore, we can eliminate the parameters λ λ λ and replace the covariance C with the modified version C mod in the likelihood.
For a non-linear model, the linear model is still a good approximation when the prior of λ λ λ is narrow around λ λ λ 0 , such that higher-order expansion can be neglected.We numerically compute the Jacobian matrix J for the 32 nuisance parameters around their fiducial values.
This method reduces the number of sampled parameters from 57 to 25, enabling quick convergence at the cost of our ability to investigate the self-calibration of those 32 nuisance parameters; we apply this method throughout this paper.

SIMULATED LIKELIHOOD ANALYSIS RESULTS
For each shear scale cut choice, ℓ shear max = 500 or 3000 or 8000, we run simulated likelihood analyses of 3×2pt, 6×2pt, 10×2pt, and 8×2pt, where the latter analysis drops the κ CMB -y and δ g -y combinations compared to 10×2pt analysis that includes all possible crosscorrelations.We use emcee (Foreman-Mackey et al. 2013) to sample the parameter space.
In this section, we present a series of simulated analyses that quantify and interpret the gain in constraining unlocked in 10×2pt analyses: We first quantify the measurement signal-to-noise ratio and cosmological constraining power in the context of our baseline halo model for different probe combination.We then design simulated analyses that explore the degradation of cosmology parameter constraints due to halo model parameter uncertainties, which limit the cosmological interpretation of small-scale measurements, and isolate the importance of baryonic feedback parameter selfcalibration.

Constraining power of multi-probe combinations
The signal-to-noise ratios (S/N), defined by S/N = d C −1 d ⊤ ⊤ ⊤ , of different measurements can be used as a simple proxy for their constraining power.For a single-parameter model, S/N is directly related to the constraining power on this parameter; for complex models, degeneracies between parameters may degrade the constraining power on parameters compared to simple S/N expectations.We evaluate the S/N with the data computed at the fiducial parameter values and show the results for different multi-probe combination with different ℓ shear max in Fig. 2. With the scale cuts described in Sect.2.1, we find the S/N to improve by ∼ 10 − 15% when very small-scale shear information is included (ℓ shear max = 8000 vs. ℓ shear max = 3000), with only minor variations between different probe combinations.The addition of all CMB lensing and tSZ (cross-)correlations increases the S/N by ∼ 20% compared to 3×2pt data.Notably, there is limited information gain from the κ CMB -y and δ g -y cross-correlations, with the S/N increasing only ∼ 5% from 8×2pt to 10×2pt.
To illustrate the cosmological constraining power of different probe combinations, we focus on Ω m and S 8 ≡ σ 8 (Ω m /Ω fid m ) 0.4 , where the power index 0.4 is chosen such that S 8 captures the well-measured combination, as parameters of interest.We use the 3×2pt, 6×2pt, and 10×2pt analyses are shown in Fig. 3, with the corresponding FoM values shown in Fig. 2. The constraints improve by ∼ 50% (in FoM) as CMB lensing and tSZ are included in the analysis while the corresponding gain in S/N is only ∼ 20%, which illustrates the importance of parameter self-calibration and breaking of degeneracies for parametric constraints.Fig. 2 also shows the limited gain in constraining power from more aggressive shear scale cuts, which is consistent with the corresponding S/N gains.
We also find negligible difference in the Ω m − S 8 constraints between 8×2pt and 10×2pt probes, which suggests that κ CMB -y and δ gy combinations may be dropped with limited impact on cosmology constraints.However, these combinations may contribute to selfcalibration of additional systematic parameters in extended models.

Impact of halo model uncertainties
To illustrate the degradation of cosmological constraints due to the systematic uncertainties, we run two additional simulated analyses: (1) "fix halo", fixing all the eight free halo parameters at their fiducial values, and (2) "cosmo only", further fixing the IA and galaxy bias parameters and sampling only the five cosmological parameters.Note that in all the scenarios, the photo-z and shear calibration parameters are analytically marginalised, which does not allow for internal calibration of these parameters and may thus underestimate the gains from decreasing model complexity.The Ω m −S 8 constraints are shown in Fig. 4 and the corresponding FoMs for "fiducial" and "fix halo" cases are 3.8 × 10 5 and 5.6 × 10 5 , respectively.
One way to identify the halo model parameter(s) most responsible for degrading the cosmological constraints is to run chains similar to the "fix halo" analyses, whereby we fix one or more parameters and compute the impact on the FoM.Since this is a computationally expensive step, we instead utilise the full parameter covariance from the fiducial case.We first approximate the posterior distribution of the fiducial case as a multivariate Gaussian with parameter covariance C p .Fixing a combination of one or more parameters p, we analytically compute the conditional covariance for the remaining parameters (Petersen & Pedersen 2012).Suppose that the fixed parameters p correspond to the covariance block C p p p , then the remaining parameters' joint posterior (i.e., the conditional distribution) will be a multivariate Gaussian with its covariance given by the Schur complement of block C p p p ,7 from which we can calculate the conditional FoMs.We find that in the limiting case when all eight halo parameters are fixed, the estimated FoM of 5.4 × 10 5 is consistent with the FoM measured directly from the "fix halo" chain.
We vary the number of fixed halo parameters N p and compute the conditional FoMs for all possible parameter combinations.We identify the combination that maximises the conditional FoM, which are shown in Tab. 3. We find that N p = 3 is sufficient for the maximal conditional FoM to reach 90% of the limiting FoM.These results identify that the parameters describing the halo mass-concentration relation (ϵ 1 , ϵ 2 , M 0 , β, Eq. 11) fractions of bound and ejected gas (M 0 , β, Eqs.3-4) contribute the most to the degradation of constraining power on Ω m and S 8 parameters.Additionally, the concentration of the gas profile, parameterised by the gas polytropic index Γ, and non-thermal contributions to gas temperature, α, can have a relatively high impact on constraining power as well.This test suggests that improved priors on the halo mass-concentration (from independent observations or/and hydrodynamic simulations) will have the highest impact on constraining power on cosmology.
As the posterior distribution of parameters is not a multivariate Gaussian, we show a simulated 10×2pt analysis with the three highest-ranked halo parameters from Tab. 3 (ϵ 1 , ϵ 2 , Γ) fixed at their fiducial values in Fig. 5 (dashed orange line).The Ω m − S 8 constraints approach the limiting case where all halo parameters are perfectly known ("fix halo", solid red contour).As these halo parameters were selected by their impact on FoM Ω m ,S 8 , the rank ordering may be different for other parameters of interest, and our example performs poorly in recovering H 0 constraining power.The "fiducial" model (blue line) samples over five cosmological parameters, 12 astrophysical nuisance parameters, and eight halo parameters (c.f.Tab. 2).The "fix halo" scenario (red shaded) fixes all halo parameters at their fiducial values, while the "cosmo only" scenario further fixes the remaining 12 nuisance parameters at their fiducial values. .Degradation of cosmological parameter constraints due to halo uncertainties 10×2pt analyses with ℓ shear max = 8000.In addition to "fiducial" model (blue line) and "fix halo" scenario (red shaded) repeated from Fig. 4, the orange dotted contour shows an intermediate model fixing only the three halo model parameters ϵ 1 , ϵ 2 , Γ.These three parameters most degrade the constraining power on Ω m and S 8 and are thus high-reward targets for external calibration.
3.3 Self-calibration of small-scale modelling Within the HMx-like halo model parameterization adopted here, the cosmological information gain from small scales is limited by (halo) modelling uncertainties, as demonstrated in Sect.3.2.When including tSZ information, the halo model parameters describing the matter field are self-calibrated from y (cross-) correlations through the halo parameters p × = {ϵ 1 , ϵ 2 , M 0 , β} that are shared between the matter and pressure model.
To isolate the impact of halo parameter self-calibration, we extend the halo parameterisation with a second copy of parameters p ′ × to decouple the halo parameters for matter and pressure.This allows us to design two scenarios that fall between the fiducial 6×2pt and 8×2pt (6×2pt, y-y, shear-y) analyses as limiting cases: • 8×2pt without halo parameter self-calibration: independent halo parameters for matter and pressure, i.e. y field modelled by p ′ × .• 8×2pt with partial halo parameter self-calibration: y field in shear-y correlations shares p × with the matter field, while the y field in y-y uses p ′ × .
The corresponding constraints on p × are shown in Fig. 6.Without self-calibration between matter field and y field (black dashed-dotted line), the 8×2pt constraints are very similar to those of a 6×2pt analysis (red line), indicating that there is no significant information gain from including tSZ.Partial self-calibration noticeably improves constraints on the halo parameters (blue dashed-dotted line), but falls short of the fiducial 8×2pt analysis with full self-calibration (orange shaded contours).
These results demonstrate the importance of joint parameterisation for matter and pressure to maximise the cosmology constraining power from tSZ cross-correlation analyses.

SUMMARY AND DISCUSSION
The forthcoming overlap of large galaxy photometric surveys and CMB experiments offers the opportunity to jointly analyse these datasets, thereby enhancing constraints on cosmological physics.In this paper, we present the first joint "10×2pt" analysis of the galaxy position and lensing shear fields from galaxy surveys like LSST, and the reconstructed CMB lensing convergence and Compton-y fields from CMB experiments like SO.These simulated analysis are based on MCMC chains with a multi-component halo model, Figure 6.The effect of halo parameter self-calibration on halo parameters β, M 0 , ϵ 1 , ϵ 2 on 6×2pt and 8×2pt analyses with ℓ shear max = 8000.Starting from the fiducial 8×2pt analysis (orange shaded) with a unified set of halo parameters for both δ m and y fields, we decouple halo parameters for matter and pressure to reduce the halo parameter self-calibration partially (blue dasheddotted) and turning it off completely (black dashed-dotted) as detailed in Sect 3.3.The 8×2pt analysis without self-calibration approaches the constraining power of a 6×2pt analysis (red solid).Note that the orange, blue and black contours all are 8×2pt analyses, i.e.have the same data S/N.non-Gaussian covariances, and extensive modelling of observational (shear calibration and photo-z uncertainties) and astrophysical (galaxy bias, intrinsic alignment, halo) systematics.
(2) 8×2pt (excluding galaxy-y and CMB lensing-y correlations) analysis results in cosmology constraints similar to those from 10×2pt, suggesting that shear-y and y-y correlations are the most valuable additions when the tSZ information is included (Sect.3.1).
(3) The small-scale modelling is limited by the halo model uncertainties.By eliminating those uncertainties, one may see another ∼50% improvement in FoM of the fiducial 10×2pt.These uncertainties are mostly attributed to parameters ϵ 1 , ϵ 2 , Γ, while M 0 and β contribute slightly less (Sect.3.2).(4) We demonstrate that the Compton-y field provides crucial halo self-calibration, which reduces the uncertainties in the small-scale matter power spectrum, hence indirectly improving the cosmological constraints (Sect.3.3).
More generally, the gain in constraining power from LSS×CMB joint analyses is caused by both (i) the increase in S/N from additional measurements and (ii) the improved conversion of S/N to parameter constraints due to enhanced self-calibration of nuisance parameters.To illustrate the latter point, we consider different multiprobe analyses that match the S/N of our 10×2pt, ℓ shear max = 8000 analysis, which amounts to artificially rescaling the survey area for different analyses.For example, a conservative, ℓ shear max = 500 3×2pt analysis would require an (impossible) survey area of f sky = 1.18 to reach the same S/N as the 10×2pt, ℓ shear max = 8000 analysis with f sky = 0.4.
Figure 7 shows the constraining power of these S/N-matched analyses, with the upper horizontal axis indicating the rescaled survey area.The different 3×2pt analyses and 6×2pt analysis employ the same set of model parameters; adding qualitatively different information in steps 1 and 2 breaks degeneracies between nuisance parameters thus improves the constraining power (at fixed S/N).Adding the tSZ auto-power spectrum to the 6×2pt analysis (step 3) requires additional nuisance parameters, that are partially degenerate with other model parameters, which at fixed S/N leads to a reduction in constraining power (i.e., in this case the tSZ model parameter are constrained at the expense of the cosmological model).Including all tSZ cross-correlations in the 10×2pt analysis then enables improved self-calibration (step 4).We illustrate the potential gain from improved priors on halo model parameters in step 5, where the extent of the blue vertical bar corresponds to the "fix halo" analysis in Fig. 4 as a limiting case.We discussed in the previous sections how fixing only three of the halo parameters recovers a large fraction of this gain.We note that the similar constraining power of the S/Nmatched 6×2pt and 10×2pt analyses is likely a coincidence, and different priors on the pressure halo model parameters would change the relative constraining power of these two analyses.The survey area required for these different S/N-matched analyses further illustrates the potential of 10×2pt analyses and the impact of further research into joint modelling approaches: even a full-sky LSS survey is insufficient to reach comparable S/N with 3×2pt analyses restricted to moderately non-linear scales (ℓ shear max = 500).With ℓ shear max = 8000, 3×2pt and 6×2pt analyses would still require f sky ∼ 0.6 surveys.The latter is in principle achievable from the ground with multiple facilities, but the associated financial and environmental costs motivate the optimisation of analysis strategies and information extraction.While the accuracy of current modelling prescriptions (like HMx) are insufficient for ambitious 10×2pt analyses, these forecasts demonstrate the outsized impact of joint analyses on constraining power and motivate future research on joint multi-probe modelling and multi-wavelength calibration of halo properties.
It is important to stress that beyond the increase of constraining power showcased in this paper, multi-probe analysis programs will be essential to assess model accuracy by enabling consistency checks between different data splits.For example, one can measure halo model parameters separately from two subsets of the 10×2pt data, namely y × {y, δ g , κ CMB }, which measures halo model parameters from y and is most sensitive to cluster-sized halos, and κ g × {κ g , δ g , κ CMB }, which measures halo model parameters from the matter power spectrum and is most sensitive to group-sized halos.Any inconsistency in these two sets of halo model parameter constraints would indicate incompleteness of the model, e.g., in the parameterisation for the halo mass-dependence of baryonic effects.
While these conclusion inherently depend on the survey and analysis choices, we demonstrate the potential of future joint analyses of galaxy surveys and CMB experiments.In the companion paper (Fang et al, in prep), we apply the 10×2pt methodology and analysis framework to Roman Space Telescope and SO as well as CMB-S4, with various survey and analysis choices.
We do not consider higher-order tidal alignment, tidal torquing models (see e.g., Blazek et al. 2015Blazek et al. , 2019)), or effective field theory models (Vlah et al. 2020), or more complicated IA modelling as a function of galaxy colour (Samuroff et al. 2019), or IA halo model (Fortuna et al. 2021).Similar to non-linear galaxy bias models, the degrees of freedom that are opened up by these models may degrade the constraining power of the galaxy survey and enhance the importance of the information carried by secondary CMB effects.2) is the galaxy-galaxy lensing (κ g -δ g ) tomography correlation with (13) being the galaxy-galaxy combinations of the 2nd lens bin with all the non-overlapping source bins at higher redshifts.( 12) zooms into the matrix of a part of the correlation between cosmic shear and galaxy-galaxy lensing.(3) is the clustering (δ g -δ g ) auto-probe matrix with 10 tomographic bins.(4) corresponds to the galaxy-CMB convergence (δ g -κ CMB ) combinations auto-probe matrix, which is comprised of 10 lens sample redshift bins.( 5) is the auto-probe correlation of the κ g -κ CMB part of the data vector, which uses the 10 source sample redshift bins.( 6) is the auto-probe correlation of the CMB lensing power spectrum, κ CMB -κ CMB .( 7) is the auto-probe correlation of the galaxy-tSZ (δ g -y) combinations, which uses the 10 lens redshift bins.( 14) zooms into the correlation matrix of all 10 galaxy-tSZ combinations with the last 3 of galaxy-galaxy lensing power spectra (κ g -δ g ), all the 10 galaxy clustering (δ g -δ g ), and all 10 δ g -y spectra, respectively from left to right.( 8) is the auto-probe correlation of the galaxy lensing convergence-tSZ (κ g -y) combinations, which uses the 10 source redshift bins.( 9) is the auto-probe correlation of the CMB lensing convergence-tSZ (κ CMB -y).( 10) is the auto-probe correlation of the tSZ power spectrum (y-y).

Figure 4 .
Figure 4. Degradation of the Ω m − S 8 constraining power for 10×2pt analyses with ℓ shear max = 8000 as a function of the model complexity.The "fiducial" model (blue line) samples over five cosmological parameters, 12 astrophysical nuisance parameters, and eight halo parameters (c.f.Tab. 2).The "fix halo" scenario (red shaded) fixes all halo parameters at their fiducial values, while the "cosmo only" scenario further fixes the remaining 12 nuisance parameters at their fiducial values.
Figure5.Degradation of cosmological parameter constraints due to halo uncertainties 10×2pt analyses with ℓ shear max = 8000.In addition to "fiducial" model (blue line) and "fix halo" scenario (red shaded) repeated from Fig.4, the orange dotted contour shows an intermediate model fixing only the three halo model parameters ϵ 1 , ϵ 2 , Γ.These three parameters most degrade the constraining power on Ω m and S 8 and are thus high-reward targets for external calibration.

Figure 7 .
Figure 7.The Figure-of-Merit (FoM) of different S/N matched-analyses.The survey area required for each analysis to achieve the same S/N as the 10×2pt, ℓ shear max = 8000 analysis with f sky = 0.4 is listed on the upper horizontal axis.See text for details.

Figure B1 .
Figure B1.The diagonal elements of the tSZ power spectrum covariance.The blue, orange and green lines show the Gaussian covariance, the supersample covariance (SSC), the connected non-Gaussian (cNG) covariance, respectively, with cNG dominates over other two components for the entire range of angular scales considered.The cNG component computed here only includes the 1-halo contributions.

Figure B2 .
Figure B2.The 10×2pt correlation matrix for LSST Y6 ×SO Y5, plotted as log 10 |R|, where R = [diag(C)] −1/2 C[diag(C)] −1/2 .We have highlighted some parts of the matrix to illustrate the correlation structure: (1) depicts the cosmic shear (κ g -κ g ) correlation matrix, comprised of 55 tomographic combinations of source bins, each with 25 Fourier ℓ-bins.(11) shows one of the tomographic combinations, and the individual ℓ 1 , ℓ 2 elements are clearly visible.(2) is the galaxy-galaxy lensing (κ g -δ g ) tomography correlation with (13) being the galaxy-galaxy combinations of the 2nd lens bin with all the non-overlapping source bins at higher redshifts.(12) zooms into the matrix of a part of the correlation between cosmic shear and galaxy-galaxy lensing.(3) is the clustering (δ g -δ g ) auto-probe matrix with 10 tomographic bins.(4) corresponds to the galaxy-CMB convergence (δ g -κ CMB ) combinations auto-probe matrix, which is comprised of 10 lens sample redshift bins.(5) is the auto-probe correlation of the κ g -κ CMB part of the data vector, which uses the 10 source sample redshift bins.(6) is the auto-probe correlation of the CMB lensing power spectrum, κ CMB -κ CMB .(7) is the auto-probe correlation of the galaxy-tSZ (δ g -y) combinations, which uses the 10 lens redshift bins.(14) zooms into the correlation matrix of all 10 galaxy-tSZ combinations with the last 3 of galaxy-galaxy lensing power spectra (κ g -δ g ), all the 10 galaxy clustering (δ g -δ g ), and all 10 δ g -y spectra, respectively from left to right.(8) is the auto-probe correlation of the galaxy lensing convergence-tSZ (κ g -y) combinations, which uses the 10 source redshift bins.(9) is the auto-probe correlation of the CMB lensing convergence-tSZ (κ CMB -y).(10) is the auto-probe correlation of the tSZ power spectrum (y-y).
Galitzki et al. 2018;Ade et al. 2019) is a CMB experiment under construction in the Atacama desert in Chile, at an altitude of 5,200 m.It is designed to observe the microwave sky in six frequency bands centred around30, 40, 90, 150, 230, and 290GHz, in order to separate the CMB from Galactic and extragalactic foregrounds.

Table 1 .
Parameters for generating the lens and source galaxy samples of LSST Y6.

Table 2 .
A list of the parameters characterising the surveys, cosmology and systematics.The fiducial values are used for generating the simulated data vectors, and the priors are used in the sampling.Uniform priors are described by U[minimum, maximum], and Gaussian priors are described through the normal distribution N(µ, σ 2 ).

Table 3 .
Estimated Ω m − S 8 Figure-of-Merit (FoM) of 10×2pt with galaxy shear scale cuts ℓ shear max = 8000 when the number N p of fixed halo parameters p increases.