# Cosmology from weak lensing, galaxy clustering, CMB lensing, and tSZ – I. 10 × 2pt modelling methodology

Xiao Fang<sup>1,2*</sup> Elisabeth Krause,<sup>2,3</sup> Tim Eifler<sup>1,2</sup> Simone Ferraro,<sup>1,4</sup> Karim Benabed,<sup>5</sup> Pranjal R. S.<sup>1,2</sup> Emma Ayçobery<sup>1,5</sup> Yohan Dubois<sup>5</sup> and Vivian Miranda<sup>1,6</sup>

Accepted 2023 December 7. Received 2023 November 14; in original form 2023 August 9

## ABSTRACT

The overlap of galaxy surveys and cosmic microwave background (CMB) experiments presents an ideal opportunity for joint cosmological data set 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 Ω<sub>m</sub> and S<sub>8</sub> over the constraints from using LSS-only probes within Λ cold dark matter (ΛCDM) model. We identify that the shear–y and y–y correlations are the most valuable additions when thermal Sunyaev–Zel'dolvich (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.

**Key words:** cosmic background radiation – cosmological parameters – large-scale structure of Universe – cosmology: theory.

## 1 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. 2018, 2022), the Hyper Suprime-Cam (HSC; Hikage et al. 2019), the Dark Energy Spectroscopic Instrument (DESI; DESI Collaboration 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 (SPHERE; 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 energy composition and structure growth of the Universe at much higher redshifts via the cosmic microwave background (CMB). Whether there exists an S<sub>8</sub> parameter difference within the standard ΛCDM model between high-z CMB measurements from *Planck* (Planck Collaboration 2020) and several low-z measurements, such as KiDS (Joudaki et al. 2018; Heymans et al. 2021), DES (Abbott et al. 2018, 2019, 2022), unWISE and *Planck* CMB lensing tomography (Krolewski et al. 2020; Krolewski, Ferraro & White 2021), and *Atacama Cosmology Telescope (ACT)* CMB lensing (Madhavacheril et al. 2023; Qu 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. 2016, 2019), along with wider and deeper galaxy surveys.

CMB secondary effects, such as CMB lensing (Lewis & Challinor 2006) and thermal Sunyaev–Zel’dolvich (tSZ) effect (Carlstrom, Holder & Reese 2002), trace the foreground gravitational potential and baryon distributions, hence are complementary probes of LSS.

<footnote n="1"><sup>1</sup>Berkeley Center for Cosmological Physics, UC Berkeley, CA 94720, USA</footnote>

<footnote n="2"><sup>2</sup>Department of Astronomy and Steward Observatory, University of Arizona, 933 North Cherry Avenue, Tucson, AZ 85721, USA</footnote>

<footnote n="3"><sup>3</sup>Department of Physics, University of Arizona, 1118 E. Fourth Street, Tucson, AZ 85721, USA</footnote>

<footnote n="4"><sup>4</sup>Lawrence Berkeley National Laboratory, One Cyclotron Road, Berkeley, CA 94720, USA</footnote>

<footnote n="5"><sup>5</sup>Sorbonne Université, CNRS, UMR7095, Institut d'Astrophysique de Paris, 98 bis Boulevard Arago, F-75014 Paris, France</footnote>

<footnote n="6"><sup>6</sup>C. N. Yang Institute for Theoretical Physics, Stony Brook University, Stony Brook, NY 11794, USA</footnote>

<footnote>*E-mail: xfang@berkeley.edu</footnote>
As the sensitivities of CMB experiments improve, these probes become increasingly powerful and relevant to precision cosmology measurements. For example, the constraints on Ω_m (the matter density parameter) and σ_8 (the amplitude of the linear power spectrum on the scale of 8 h⁻<fnref n="1" /> Mpc) from the ACT CMB lensing alone have achieved precision comparable to current weak lensing results (DES Y3/KiDS-1000/HSC-Y3; see section 3 of Madhavacheril et al. 2023). Besides the cosmological information carried by the CMB secondary effects, cross-correlating these signals with LSS probes, such as galaxy positions and shapes, provides measurements with an independent set of systematic errors and improves the calibration of systematics in galaxy surveys. The role of CMB lensing cross-correlations has been emphasized in the context of SO (Ade et al. 2019) and CMB-S4 (Abazajian et al. 2016). There is a wealth of literature studying the synergy of CMB lensing and galaxy fields and the synergy of tSZ with galaxy fields. The theory of the former has recently been studied in the context of upcoming galaxy surveys and CMB experiments (e.g. Schaan et al. 2017; Schmittfull & Seljak 2018; Schaan, Ferraro & Seljak 2020; Fang et al. 2022; Wenzl et al. 2022), and the synergies have been carried out in experiments, including SDSS + *Planck* (e.g. Singh, Mandelbaum & Brownstein 2017), DES + *Planck* (e.g. Giannantonio et al. 2016; Abbott et al. 2019; Omori et al. 2019a, b), DES + *Planck* + SPT (e.g. Abbott et al. 2023; Chang et al. 2023; Omori et al. 2023), unWISE + *Planck* (e.g. Krolewski, Ferraro & White 2021), HSC + *Planck* (e.g. Miyatake et al. 2022), KiDS + ACT + *Planck* (e.g. Robertson et al. 2021), DESI/DESI-like + *Planck* (e.g. Hang et al. 2021; Kitanidis & White 2021; White et al. 2022). The synergies of tSZ with galaxy fields have also been explored both in theory (e.g. Pandey, Baxter & Hill 2020a; Nicola et al. 2022) and in observations, such as Canada France Hawaii Lensing Survey + *Planck* (Ma et al. 2015), 2MASS + *Planck* (Makiya, Ando & Komatsu 2018), 2MASS + WISE × SuperCOSMOS + *Planck* (Koukoufilippas et al. 2020), SDSS + *Planck* (e.g. Vikram, Lidz & Jain 2017; Hill et al. 2018; Chiang et al. 2020; Tanimura et al. 2020), SDSS + ACT (e.g. Schaan et al. 2021), KiDS + *Planck* + ACT (e.g. Yan et al. 2021; Tröster et al. 2022), and DES + *Planck* + ACT (e.g. Gatti et al. 2022; Pandey et al. 2022a). However, there has not been any fully joint analysis combining all these probes, and simultaneously constraining cosmological and baryonic properties with a consistent halo model.

This paper explores the synergies of a CMB × LSS joint analysis using the SO and Vera C Rubin Observatory's 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).<fnref n="1" /> We present our theoretical modelling of the combined '10 × 2pt' probes in Section 2, including the analysis choices, data vector, and analytic covariance modelling. We carry out the simulated likelihood forecast analysis for LSST + SO in Section 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 Section 4.

In a companion paper (Fang et al. in preparation, 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.

## 2 MULTISURVEY MULTIPROBE 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).

### 2.1 Survey specifications

#### 2.1.1 SO

The SO (Galitzki et al. 2018; Ade et al. 2019) is a CMB experiment under construction in the Atacama desert in Chile, at an altitude of 5200 m. It is designed to observe the microwave sky in six frequency bands centred around 30, 40, 90, 150, 230, and 290 GHz, in order to separate the CMB from Galactic and extragalactic foregrounds.

The nominal design of the observatory will include one 6-m large-aperture telescope (LAT; Parshley et al. 2018; Xu et al. 2021) and three small-aperture 0.5-m telescopes (SAT s; Ali et al. 2020). The LAT will produce temperature and polarization maps of the CMB with ~arcmin resolution over 40 per cent 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 SAT s, 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 Appendix B.

#### 2.1.2 Vera C Rubin Observatory's LSST

Vera C Rubin Observatory's LSST will rapidly and repeatedly cover a ~18 000 deg² footprint in six optical bands (320–1050 nm). With a single exposure depth of 24.7 r-band magnitude (5σ point source), the 10 yr 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 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 per cent of the sky or Ω_s = 16 500 deg². 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.

<footnote n="1">¹We note that the term N×2pt analysis is ambiguous, as N = 3/6/10 different two-point correlations can be constructed from any two/three/four different tracer fields; e.g. other 6 × 2pt combinations are possible, such as replacing CMB lensing with cluster density (To et al. 2021).</footnote>
<table><caption>**Table 1.** Parameters for generating the lens and source galaxy samples of LSST Y6.</caption><thead><tr><th>Parameter</th><th>x = lens</th><th>x = source</th></tr></thead><tbody><tr><td>n&#773;<sub>x</sub> (arcmin<sup>-2</sup>)</td><td>41.1</td><td>23.2</td></tr><tr><td>z<sub>0</sub></td><td>0.274</td><td>0.178</td></tr><tr><td>&alpha;</td><td>0.907</td><td>0.798</td></tr></tbody></table>

<figure><img src="image_1.png" /><figcaption>**Figure 1.** The equally binned true redshift distributions $n_x(z) = \frac{dN_x}{dz d\Omega}$ of the source and lens samples of LSST Year 6, normalized by the corresponding binned effective number density $\bar{n}_x^i$ after redshift cuts. A fiducial Gaussian photo-z error has been convolved with each tomographic bin as described by equation (A14).</figcaption></figure>

We parametrize the photometric redshift distributions of both the lens and source samples as

$$n_x(z_{ph}) = \frac{dN_x}{dz_{ph} d\Omega} \propto z_{ph}^2 \exp[-(z_{ph}/z_0)^\alpha], \quad x \in \{\text{lens, source}\} \quad (1)$$

normalized by the effective number density $\bar{n}_x$. $N_x$ is the number counts of lens/source galaxies, $z_{ph}$ is the photometric redshift, $\Omega$ is the solid angle. The parameter values of $\bar{n}_x$, $z_0$, $\alpha$ are given in Table 1, which are updated values from the LSST-DESC's Observing Strategy<fnref n="2" /> (Lochner et al. 2018). These values, slightly different from the DESC-SRD, are computed for the more recent optimization 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.

## 2.2 Multiprobe modelling

We extend the existing 6 × 2pt model of the CosmoLike<fnref n="3" /> 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 (Ma & Fry 2000; Peacock & Smith 2000; Seljak 2000; Cooray & Sheth 2002; Hill & Pajer 2013; Hill & Spergel 2014, c.f. Appendix A1.1) approach and specify all systematics at the level of the individual

<table><caption>**Table 2.** A list of the parameters characterizing the surveys, cosmology, and systematics.</caption><thead><tr><th>Parameters</th><th>Fiducial</th><th>Prior</th></tr></thead><tbody><tr><td><strong>Survey</strong></td><td></td><td></td></tr><tr><td>Ω<sub>s</sub> (deg<sup>2</sup>)</td><td>16 500</td><td>Fixed</td></tr><tr><td>σ<sub>e</sub></td><td>0.26/component</td><td>Fixed</td></tr><tr><td><strong>Cosmology</strong></td><td></td><td></td></tr><tr><td>Ω<sub>m</sub></td><td>0.3156</td><td>U[0.05, 0.6]</td></tr><tr><td>σ<sub>8</sub></td><td>0.831</td><td>U[0.5, 1.1]</td></tr><tr><td>n<sub>s</sub></td><td>0.9645</td><td>U[0.85, 1.05]</td></tr><tr><td>Ω<sub>b</sub></td><td>0.0492</td><td>U[0.04, 0.055]</td></tr><tr><td>h<sub>0</sub></td><td>0.6727</td><td>U[0.4, 0.9]</td></tr><tr><td><strong>Galaxy bias</strong></td><td></td><td></td></tr><tr><td>b<sup>i</sup> (i=1,...,N<sub>tomo</sub>)</td><td>0.95/G(z<sup>i</sup>)</td><td>U[0.4, 3]</td></tr><tr><td><strong>Photo-z</strong></td><td></td><td></td></tr><tr><td>Δ<sup>i</sup><sub>z,lens</sub> (i=1,...,N<sub>tomo</sub>)</td><td>0</td><td>N(0, 0.001<sup>2</sup>)</td></tr><tr><td>σ<sub>z,lens</sub></td><td>0.03</td><td>N(0.03, 0.003<sup>2</sup>)</td></tr><tr><td>Δ<sup>i</sup><sub>z,source</sub> (i=1,...,N<sub>tomo</sub>)</td><td>0</td><td>N(0, 0.001<sup>2</sup>)</td></tr><tr><td>σ<sub>z,source</sub></td><td>0.05</td><td>N(0.05, 0.003<sup>2</sup>)</td></tr><tr><td><strong>Shear calibration</strong></td><td></td><td></td></tr><tr><td>m<sup>i</sup> (i=1,...,N<sub>tomo</sub>)</td><td>0</td><td>N(0, 0.003<sup>2</sup>)</td></tr><tr><td><strong>IA</strong></td><td></td><td></td></tr><tr><td>A<sub>IA</sub></td><td>0.5</td><td>U[-5, 5]</td></tr><tr><td>η<sub>IA</sub></td><td>0</td><td>U[-5, 5]</td></tr><tr><td><strong>Halo and gas parameters</strong></td><td></td><td></td></tr><tr><td>Γ</td><td>1.17</td><td>U[1.05, 1.35]</td></tr><tr><td>β</td><td>0.6</td><td>U[0.2, 1.0]</td></tr><tr><td>log<sub>10</sub>M<sub>0</sub></td><td>14.0</td><td>U[12.5, 15.0]</td></tr><tr><td>α</td><td>1.0</td><td>U[0.5, 1.5]</td></tr><tr><td>log<sub>10</sub>T<sub>w</sub></td><td>6.5</td><td>U[6.0, 7.0]</td></tr><tr><td>ε<sub>1</sub></td><td>0</td><td>U[-0.8, 0.8]</td></tr><tr><td>ε<sub>2</sub></td><td>0</td><td>U[-0.8, 0.8]</td></tr><tr><td>f<sub>H</sub></td><td>0.752</td><td>U[0.7, 0.8]</td></tr><tr><td>M<sub>min</sub>, M<sub>max</sub></td><td>10<sup>6</sup>, 10<sup>17</sup> M<sub>☉</sub>/h</td><td>Fixed</td></tr><tr><td>A<sub>*</sub></td><td>0.03</td><td>Fixed</td></tr><tr><td>log<sub>10</sub>M<sub>*</sub></td><td>12.5</td><td>Fixed</td></tr><tr><td>σ<sub>*</sub></td><td>1.2</td><td>Fixed</td></tr></tbody></table>

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[\text{minimum}, \text{maximum}]$, and Gaussian priors are described through the normal distribution $N(\mu, \sigma^2)$.

observables, i.e. the density contrast of lens galaxies $\delta_g$, the lensing convergence of source galaxies $\kappa_g$, the CMB lensing convergence $\kappa_{CMB}$, and the Compton-y field $y$. We summarize the well-established computation of angular power spectra and our model ingredients for $\delta_g$, $\kappa_g$, and $\kappa_{CMB}$, including astrophysical systematics linear galaxy bias model, intrinsic alignments (IA) using the ‘non-linear linear alignment’ model (Hirata & Seljak 2004; Bridle & King 2007), and observational systematics, in Appendix A. Combined with tSZ and baryonic feedback modelling described in this subsection, our 10 × 2pt model has 57 free parameters, which are listed in Table 2.

### 2.2.1 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 parametrization 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 (McCarthy

<footnote n="2">²https://github.com/LSSTDESC/ObsStrat/tree/static/static</footnote>

<footnote n="3">³https://github.com/CosmoLike</footnote>
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, active galactic nuclei, 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 per cent and worse for the pressure auto-power spectrum (Mead et al. 2020). While Stage-IV data analyses will certainly require refinements of the model parametrization and parameter calibration, we use the HMX parametrization only as a toy model to qualitatively explore information content in the non-linear regime as well as the impact of parameter self-calibration on cosmology constraints.

The halo model implementation in CosmoLike uses the Navarro, Frenk & White (1997) profile and the Duffy et al. (2008) halo mass-concentration relation to model the halo density profile $\tilde{u}_{\delta}(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 $\Delta\bar{\rho} = 200$ halo definition, so that the halo radius $r_{200}$ satisfies $M = 4\pi r_{200}^3 \bar{\rho} \Delta\bar{\rho}/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 data sets improves the constraints of the parameters of interest. We expect that, even with its limitations, our model is representative enough to allow this study.

### 2.2.2 Electron pressure

The Fourier-transformed y halo profile, $\tilde{u}_y(k, M, z)$, which is the building block for halo model calculations of tSZ (cross-)power spectra (c.f. Appendix A1.1) is then computed from the electron pressure profile $P_e(M, r, z)$,

$$ \tilde{u}_y(k, M, z) = 4\pi \int_0^\infty dr \, r^2 j_0(kr) P_e(M, r, z). \quad (2) $$

The HMX halo model parametrization 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 (2h) term.

The fractions of bound and ejected gas in halos are parametrized as

$$ f_{\text{bnd}}(M) = \frac{\Omega_b}{\Omega_m} \frac{1}{1 + (M_0/M)^\beta}, \quad (3) $$

$$ \begin{aligned} f_{\text{ejc}}(M) &= \frac{\Omega_b}{\Omega_m} - f_{\text{bnd}}(M) - f_*(M) \\ &= \frac{\Omega_b}{\Omega_m} \frac{1}{1 + (M_0/M)^{-\beta}} - A_* \exp\left[-\frac{\log_{10}(M/M_*)}{2\sigma_*^2}\right], \end{aligned} \quad (4) $$

where $M_0$ and $\beta$ are an empirical parametrization 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, parametrized by $A_*, M_*, \sigma_*$). $\Omega_b$ and $\Omega_m$ are the cosmic baryon and matter density parameters.

The bound gas density profile $\rho_{\text{bnd}}$ can then be written as

$$ \rho_{\text{bnd}}(M, r, z) = A(M, z) M f_{\text{bnd}}(M) \left[ \frac{\ln(1 + r/r_s)}{r/r_s} \right]^{1/(\Gamma - 1)}, \quad (5) $$

where the radial dependence follows the KS profile with polytropic index $\Gamma$ (Komatsu & Seljak 2002). The normalization is given by $\int_0^{r_v} dr \, 4\pi r^2 \rho_{\text{bnd}}(M, r, z) = M f_{\text{bnd}}(M)$, leading to

$$ [A(M, z)]^{-1} = \int_{0}^{r_v} dr \, 4\pi r^2 \left[ \frac{\ln(1 + r/r_s)}{r/r_s} \right]^{1/(\Gamma - 1)}. \quad (6) $$

In the KS profile the bound gas temperature $T_g$ is determined by hydrostatic equilibrium,

$$ T_g(M, r, z) = T_v(M, z) \frac{\ln(1 + r/r_s)}{r/r_s}, \quad (7) $$

with $T_v(M, z)$ the virial temperature

$$ \frac{3}{2} k_B T_v(M, z) = \alpha \frac{G M m_p \mu_p}{r_v} (1 + z), \quad (8) $$

where $m_p$ is the proton mass, $\mu_p = 4/(3 + 5f_H)$, with $f_H$ the hydrogen mass fraction, and $\alpha$ 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

$$ \begin{aligned} P_{e,bnd}(M, r, z) &= \frac{\rho_{bnd}(M, r, z)}{m_p \mu_e} k_B T_g(M, r, z) \\ &= \frac{\rho_{bnd}(M, r, z)}{m_p \mu_e} k_B T_v(M, z) \frac{\ln(1 + r/r_s)}{r/r_s}, \end{aligned} \quad (9) $$

where $\mu_e = 2/(1 + f_H)$, and $k_B$ is the Boltzmann constant. The density profile for the ejected gas $\rho_{ejc}$ in 2h 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),

$$ P_{e, ejc}(M, r, z) = \frac{\rho_{ejc}(M, r, z)}{m_p \mu_e} k_B T_w = \frac{M_{ejc}(M, z) \delta_D(r)}{m_p \mu_e} k_B T_w. \quad (10) $$

### 2.2.3 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)$)

$$ \begin{aligned} c(M, z) &= c_D(M, z) \left[ 1 + \epsilon_1 + (\epsilon_2 - \epsilon_1) \frac{f_{\text{bnd}}(M)}{\Omega_b / \Omega_m} \right] \\ &= c_D(M, z) \left[ 1 + \epsilon_1 + \frac{\epsilon_2 - \epsilon_1}{1 + (M_0/M)^\beta} \right], \end{aligned} \quad (11) $$

where $\epsilon_1, \epsilon_2$ are free parameters and $\epsilon_1 = \epsilon_2 = 0$ reduces to the unmodified case.

### 2.2.4 Stellar component

Throughout the analyses presented here, the HMX halo model parameters $A_*, M_*, \sigma_*$, which describe the stellar mass fraction, are held fixed as the 10 × 2pt data vector has almost no sensitivity to these parameters.

## 2.3 Simulated likelihood analysis methodology

To simulate joint analyses of multisurvey multiprobe data, we perform simulated parameter inference assuming a Gaussian likelihood
for the data d given a point p in cosmological and nuisance parameter space,

$$L(d|p) \propto e^{-\chi^2/2}, \quad \text{where} \quad \chi^2 = [d - m(p)]^\top C^{-1} [d - m(p)], \quad (12)$$

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 realization noise and allows us to focus on the constraining power of different analysis choices. The analytic Fourier space $10 \times 2pt$ covariance calculation is summarized in Appendix B and we make available an implementation CosmoCov_Fourier,<fnref n="4" /> which is an extension of CosmoCov<fnref n="5" /> (Fang, Eifler & Krause 2020a).

### 2.3.1 Scale cuts

We compute all two-point functions in 25 logarithmically spaced Fourier mode bins ranging from $\ell_{min} = 20$ to $\ell_{max} = 8000$, and impose the following set of $\ell$-cuts,

(i) **Galaxy clustering $\delta_g$** scale cuts are driven by the modelling inaccuracy of non-linear galaxy biasing. Following the DESC-SRD, we adopt $\ell_{max}^i = k_{max} \chi(z^i)$, where $k_{max} = 0.3 h/Mpc$ and $\bar{z}^i$ is the mean redshift of the lens bin i.

(ii) **Weak lensing $\kappa_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 $\ell_{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 $\ell_{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 (equation 11) implemented here to reach these scale cuts, these different $\ell_{max}$ choices will illustrate the interplay of scale cuts and parameter degeneracies.

(iii) **CMB lensing $\kappa_{CMB}$** measurements are challenging at high $\ell$ due to lensing reconstruction challenges from foregrounds; similarly, various systematic errors can limit the reconstruction of low-$\ell$ lensing modes (Darwish et al. 2021). We adopt the scale cuts $\ell_{min} = 100$, $\ell_{max} = 3000$.<fnref n="6" /> More details on the lensing reconstruction are presented in Appendix B.

(iv) **tSZ y** measurements are limited by atmospheric noise at low-$\ell$ (in the case of SO) and component separation at high-$\ell$ (Ade et al. 2019). We adopt the somewhat optimistic scale cuts $\ell_{min} = 80$, $\ell_{max} = 8000$, noting in particular that contamination from the cosmic infrared background (CIB) will need to be modelled and (likely) marginalized over.

For cross-power spectra, we adopt the more restrictive scale cut combination of the two fields. We further exclude $\delta_g - \kappa_g$ combinations without cosmological signal, i.e. with the lens tomography bin at higher redshift than the source tomography bin.

### 2.3.2 Analytical marginalization of photo-z and shear calibration parameters

To speed up our forecast analyses, we analytically marginalize the 22 photo-z parameters and 10 shear calibration parameters, assuming that constraints on these parameters are dominated by their (Gaussian) priors listed in Table 2. This analytic marginalization was derived in the context of cosmological surveys in Bridle et al. (2002), which we summarize 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, \lambda)$ is a multivariate Gaussian with covariance $C$, i.e. $d|m(p, \lambda) \sim N(m(p, \lambda), C)$, where $p, \lambda$ are parameters of the model, and $\lambda$ contains the parameters we want to marginalize over, with a joint prior distribution $p(\lambda)$. Marginalizing the posterior distribution over $\lambda$ has a simple analytic solution if we assume (1) the model is linear, and (2) the prior $p(\lambda)$ is a multivariate Gaussian distribution.

Assuming that $p(\lambda)$ follows $N(\lambda_0, C_\lambda)$, and a linear model

$$m(p, \lambda) = m(p, \lambda_0) + J \cdot (\lambda - \lambda_0), \quad (13)$$

where the Jacobian matrix $J$ is defined by $J = \nabla_{\lambda} m(p, \lambda)|_{\lambda=\lambda_0}$. The marginalized posterior can be derived as $d|m(p, \lambda_0) \sim N(m(p, \lambda_0), C_{mod})$, where

$$C_{\text{mod}} = C + J C_{\lambda} J^{\top}. \quad (14)$$

Therefore, we can eliminate the parameters $\lambda$ 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 $\lambda$ is narrow around $\lambda_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.

## 3 SIMULATED LIKELIHOOD ANALYSIS RESULTS

For each shear scale cut choice, $\ell_{max}^{\text{shear}} = 500$ or 3000 or 8000, we run simulated likelihood analyses of $3 \times 2pt$, $6 \times 2pt$, $10 \times 2pt$, and $8 \times 2pt$, where the latter analysis drops the $\kappa_{CMB}-y$ and $\delta_g-y$ combinations compared to $10 \times 2pt$ analysis that includes all possible cross-correlations. 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 \times 2pt$ analyses: we first quantify the measurement signal-to-noise ratio (S/N) 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 self-calibration.

## 3.1 Constraining power of multiprobe combinations

The S/N, defined by $S/N = \sqrt{d C^{-1} d^\top}$, 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

<footnote n="4">⁴https://github.com/CosmoLike/CosmoCov_Fourier</footnote>

<footnote n="5">⁵https://github.com/CosmoLike/CosmoCov</footnote>

<footnote n="6">⁶While the interpretation of CMB lensing (cross-)spectra is of course also limited by the same high-$\ell$ processes that set the galaxy weak lensing scale cuts, the fractional impact of these processes on the power spectrum decreases with redshift. Hence $\ell_{max} = 3000$ for $\kappa_g$ would correspond to a less restrictive scale cut for $\kappa_{CMB}$.</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 2.** The S/N and $\Omega_m-S_8$ Figure-of-Merit (FoM) of different combined probes for three shear scale cut choices $\ell_{\text{max}}^{\text{shear}} = 500, 3000$, and 8000. The dotted line shows the FoM $\propto (\text{S/N})^2$ scaling.</figcaption></figure>

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 multiprobe combination with different $\ell_{\text{max}}^{\text{shear}}$ in Fig. 2. With the scale cuts described in Section 2.1, we find the S/N to improve by ~10–15 per cent when very small-scale shear information is included ($\ell_{\text{max}}^{\text{shear}} = 8000$ vs. $\ell_{\text{max}}^{\text{shear}} = 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 per cent compared to 3 × 2pt data. Notably, there is limited information gain from the $\kappa_{\text{CMB}}-y$ and $\delta_g-y$ cross-correlations, with the S/N increasing only ~5 per cent from 8 × 2pt to 10 × 2pt.

To illustrate the cosmological constraining power of different probe combinations, we focus on $\Omega_m$ and $S_8 \equiv \sigma_8(\Omega_m/\Omega_m^{\text{fid}})^{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 Figure-of-Merit (FoM) in the $\Omega_m-S_8$ subspace, $\text{FoM}_{\Omega_m, S_8} = 1/\sqrt{\det \text{Cov}(\Omega_m, S_8)}$, where the covariances $\text{Cov}(\Omega_m, S_8)$ are computed using the last 500k samples of the Markov chain Monte Carlo (MCMC) chains, as a simple metric for constraining power. The cosmological constraints from 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 per cent (in FoM) as CMB lensing and tSZ are included in the analysis while the corresponding gain in S/N is only ~20 per cent, 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 $\Omega_m-S_8$ constraints between 8 × 2pt and 10 × 2pt probes, which suggests that $\kappa_{\text{CMB}}-y$ and $\delta_g-y$ combinations may be dropped with limited impact on cosmology constraints. However, these combinations may contribute to self-calibration of additional systematic parameters in extended models.

## 3.2 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

<figure><img src="image_3.png" /><figcaption>**Figure 3.** Increase in LSST Y6 + SO Y5 constraining power as the data vector increases from 3 × 2pt (black dashed) to 6 × 2pt (red dashed) to 10 × 2pt (blue shaded), for galaxy shear scale cuts $\ell_{\text{max}}^{\text{shear}} = 3000$.</figcaption></figure>

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 marginalized, which does not allow for internal calibration of these parameters and may thus underestimate the gains from decreasing model complexity. The $\Omega_m-S_8$ constraints are shown in Fig. 4 and the corresponding FoMs for ‘fiducial’ and ‘fix halo’ cases are $3.8 \times 10^5$ and $5.6 \times 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 utilize 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 $\tilde{p}$, we analytically compute the conditional covariance for the remaining parameters (Petersen & Pedersen 2012). Suppose that the fixed parameters $\tilde{p}$ correspond to the covariance block $C_{\tilde{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_{\tilde{p}}$,<fnref n="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 \times 10^5$ is consistent with the FoM measured directly from the ‘fix halo’ chain.

We vary the number of fixed halo parameters $N_{\tilde{p}}$ and compute the conditional FoMs for all possible parameter combinations. We identify the combination that maximizes the conditional FoM, which are shown in Table 3. We find that $N_{\tilde{p}} = 3$ is sufficient for the maximal conditional FoM to reach 90 per cent of the limiting FoM. These results identify that the parameters describing the halo mass–

<footnote n="7">⁷For a block matrix $A = (A_{11} \ A_{12})$, the Schur complement of block $A_{11}$ is given by $A_{22} - A_{21}A_{11}^{-1}A_{12}$.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 4.** Degradation of the $\Omega_m-S_8$ constraining power for 10 × 2pt analyses with $\ell_{\text{max}}^{\text{shear}} = 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. Table 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.</figcaption></figure>

<table><caption>**Table 3.** Estimated $\Omega_m-S_8$ FoM of 10 × 2pt with galaxy shear scale cuts $\ell_{\text{max}}^{\text{shear}} = 8000$ when the number $N_{\tilde{p}}$ of fixed halo parameters $\tilde{p}$ increases.</caption><thead><tr><th>$N_{\tilde{p}}$</th><th>Fixed parameters $\tilde{p}$</th><th>FoM<sub>Ω<sub>m</sub>,S<sub>8</sub></sub>/FoM<sub>limit</sub></th></tr></thead><tbody><tr><td>0</td><td>None</td><td>70 per cent</td></tr><tr><td rowspan="3">1</td><td>$\epsilon_1$</td><td>76 per cent</td></tr><tr><td>$\Gamma$</td><td>75 per cent</td></tr><tr><td>$x \in \{\alpha, \beta, \epsilon_2, f_H, M_0\}$</td><td>71 per cent</td></tr><tr><td rowspan="3">2</td><td>$\epsilon_1, M_0$</td><td>84 per cent</td></tr><tr><td>$\epsilon_1, \beta$</td><td>84 per cent</td></tr><tr><td>$\epsilon_1, \epsilon_2$</td><td>82 per cent</td></tr><tr><td rowspan="7">3</td><td>$\epsilon_1, \epsilon_2, \Gamma$</td><td>91 per cent</td></tr><tr><td>$\epsilon_1, \epsilon_2, \beta$</td><td>87 per cent</td></tr><tr><td>$\epsilon_1, M_0, \Gamma$</td><td>87 per cent</td></tr><tr><td>$\epsilon_1, \beta, \Gamma$</td><td>87 per cent</td></tr><tr><td>$\epsilon_1, \beta, \alpha$</td><td>87 per cent</td></tr><tr><td>$\epsilon_1, M_0, \alpha$</td><td>87 per cent</td></tr><tr><td>$\epsilon_1, \epsilon_2, M_0$</td><td>86 per cent</td></tr></tbody></table>

concentration relation ($\epsilon_1, \epsilon_2, M_0, \beta$, equation 11) fractions of bound and ejected gas [$M_0, \beta$, equations (3)–(4)] contribute the most to the degradation of constraining power on $\Omega_m$ and $S_8$ parameters. Additionally, the concentration of the gas profile, parametrized by the gas polytropic index $\Gamma$, and non-thermal contributions to gas temperature, $\alpha$, 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

<figure><img src="image_5.png" /><figcaption>**Figure 5.** Degradation of cosmological parameter constraints due to halo model uncertainties 10 × 2pt analyses with $\ell_{\text{max}}^{\text{shear}} = 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 $\epsilon_1, \epsilon_2, \Gamma$. These three parameters most degrade the constraining power on $\Omega_m$ and $S_8$ and are thus high-reward targets for external calibration.</figcaption></figure>

highest-ranked halo parameters from Table 3 ($\epsilon_1, \epsilon_2, \Gamma$) fixed at their fiducial values in Fig. 5 (dashed orange line). The $\Omega_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$_{\Omega_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.

## 3.3 Self-calibration of small-scale modelling

Within the HMX-like halo model parametrization adopted here, the cosmological information gain from small scales is limited by (halo) modelling uncertainties, as demonstrated in Section 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_x = \{\epsilon_1, \epsilon_2, M_0, \beta\}$ that are shared between the matter and pressure model.

To isolate the impact of halo parameter self-calibration, we extend the halo parametrization with a second copy of parameters $p'_x$ 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:

(i) 8 × 2pt without halo parameter self-calibration: independent halo parameters for matter and pressure, i.e. y field modelled by $p'_x$.

(ii) 8 × 2pt with partial halo parameter self-calibration: y field in shear-y correlations shares $p_x$ with the matter field, while the y field in y-y uses $p'_x$.

The corresponding constraints on $p_x$ are shown in Fig. 6. Without self-calibration between matter field and y field (black dash-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
<figure><img src="image_6.png" /><figcaption>**Figure 6.** The effect of halo parameter self-calibration on halo parameters $\beta$, $M_0$, $\epsilon_1$, $\epsilon_2$ on 6 × 2pt and 8 × 2pt analyses with $\ell_{\text{max}}^{\text{shear}} = 8000$. Starting from the fiducial 8 × 2pt analysis (orange shaded) with a unified set of halo parameters for both $\delta_m$ and $\gamma$ fields, we decouple halo parameters for matter and pressure to reduce the halo parameter self-calibration partially (blue dash-dotted) and turning it off completely (black dash-dotted) as detailed in Section 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.</figcaption></figure>

gain from including tSZ. Partial self-calibration noticeably improves constraints on the halo parameters (blue dash-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 parametrization for matter and pressure to maximize the cosmology constraining power from tSZ cross-correlation analyses.

# 4 SUMMARY AND DISCUSSION

The forthcoming overlap of large galaxy photometric surveys and CMB experiments offers the opportunity to jointly analyse these data sets, 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 multicomponent halo model, non-Gaussian covariances, and extensive modelling of observational (shear calibration and photo-z uncertainties) and astrophysical (galaxy bias, IA, halo) systematics.

As an example, we simulate the LSST Y6 + SO Y5 joint analysis. For this specific analysis setup, our main findings are summarized as follows:

(i) Within ΛCDM, the 10 × 2pt constraints of $\Omega_m$ and $S_8$ improve by around 70 per cent in FoM from 3 × 2pt, or around 30 per cent from 6 × 2pt (Section 3.1).

(ii) 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

<figure><img src="image_7.png" /><figcaption>**Figure 7.** The FoM of different S/N matched analyses. The survey area required for each analysis to achieve the same S/N as the 10 × 2pt, $\ell_{\text{max}}^{\text{shear}} = 8000$ analysis with $f_{\text{sky}} = 0.4$ is listed on the upper horizontal axis. See text for details.</figcaption></figure>

most valuable additions when the tSZ information is included (Section 3.1).

(iii) The small-scale modelling is limited by the halo model uncertainties. By eliminating those uncertainties, one may see another ~50 per cent improvement in FoM of the fiducial 10 × 2pt. These uncertainties are mostly attributed to parameters $\epsilon_1$, $\epsilon_2$, $\Gamma$, while $M_0$ and $\beta$ contribute slightly less (Section 3.2).

(iv) 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 (Section 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, $\ell_{\text{max}}^{\text{shear}} = 8000$ analysis, which amounts to artificially rescaling the survey area for different analyses. For example, a conservative, $\ell_{\text{max}}^{\text{shear}} = 500$ 3 × 2pt analysis would require an (impossible) survey area of $f_{\text{sky}} = 1.18$ to reach the same S/N as the 10 × 2pt, $\ell_{\text{max}}^{\text{shear}} = 8000$ analysis with $f_{\text{sky}} = 0.4$.

Fig. 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/N-matched 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 ($\ell_{\text{shear}}^{\text{max}} = 500$). With $\ell_{\text{shear}}^{\text{max}} = 8000$, 3 × 2pt and 6 × 2pt analyses would still require $f_{\text{sky}} \sim 0.6$ surveys. The latter is in principle achievable from the ground with multiple facilities, but the associated financial and environmental costs motivate the optimization 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 multiprobe modelling and multiwavelength calibration of halo properties.

It is important to stress that beyond the increase of constraining power showcased in this paper, multiprobe 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 \times \{y, \delta_g, \kappa_{\text{CMB}}\}$, which measures halo model parameters from y and is most sensitive to cluster-sized halos, and $\kappa_g \times \{\kappa_g, \delta_g, \kappa_{\text{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 parametrization 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 preparation), 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.

## ACKNOWLEDGEMENTS

We thank Shivam Pandey, Emmanuel Schaan, Blake Sherwin, and Chun-Hao To for helpful discussions. We thank the anonymous referee for helpful comments that improved this manuscript.

XF is supported by the BCCP fellowship at the Berkeley Center for Cosmological Physics. EK is supported in part by Department of Energy grant DE-SC0020247, the David and Lucile Packard Foundation and an Alfred P Sloan Research Fellowship. TE is supported in part by Department of Energy grant DE-SC0020215. SF is supported by the Physics Division of Lawrence Berkeley National Laboratory and by the U.S. Department of Energy (DOE), Office of Science, under contract No. DE-AC02-05CH11231. KB, EA, and YD acknowledge support by the CNES, the French national space agency. EA and PRS are supported in part by a PhD Joint Programme between the CNRS and the University of Arizona. Calculations in this paper use High Performance Computing (HPC) resources supported by the University of Arizona TRIF, UITS, and RDI and maintained by the UA Research Technologies department.

## DATA AVAILABILITY

The data underlying this paper will be shared on reasonable request to the corresponding author.

## REFERENCES

* Abazajian K. N. et al., 2016, preprint (arXiv:1610.02743)

* Abazajian K. et al., 2019, preprint (arXiv:1907.04473)

* Abbott T. M. C. et al., 2018, Phys. Rev. D, 98, 43526

* Abbott T. M. C. et al., 2019, Phys. Rev. D, 100, 23541

* Abbott T. M. C. et al., 2022, Phys. Rev. D, 105, 23520

* Abbott T. M. C. et al., 2023, Phys. Rev. D, 107, 23531

* Ade P. et al., 2019, J. Cosmol. Astropart. Phys., 2019, 56

* Aiola S. et al., 2020, J. Cosmol. Astropart. Phys., 2020, 47

* Akeson R. et al., 2019, preprint (arXiv:1902.05569)

* Ali A. M. et al., 2020, Journal of Low Temperature Physics, 200, 461

* Balkenhol L. et al., 2023, Phys. Rev. D 108 023510

* Benson B. A. et al., 2014, in Holland W. S., Zmuidzinas J. eds, Proc. SPIE Conf. Ser. Vol. 9153, Millimeter, Submillimeter, and Far-Infrared Detectors and Instrumentation for Astronomy VII. SPIE, Bellingham. p. 91531P

* Blazek J., Vlah Z., Seljak U., 2015, J. Cosmol. Astropart. Phys., 2015, 15

* Blazek J. A., MacCrann N., Troxel M. A., Fang X., 2019, Phys. Rev. D, 100, 103506

* Bridle S., King L., 2007, New Journal of Physics, 9, 444

* Bridle S. L., Crittenden R., Melchiorri A., Hobson M. P., Kneissl R., Lasenby A. N., 2002, MNRAS, 335, 1193

* Carlstrom J. E., Holder G. P., Reese E. D., 2002, ARA&A, 40, 643

* Chang C. et al., 2023, Phys. Rev. D, 107, 23530

* Chiang Y.-K., Makiya R., Ménard B., Komatsu E., 2020, ApJ, 902, 56

* Choi S. K. et al., 2020, J. Cosmol. Astropart. Phys., 2020, 45

* Cooray A., Sheth R., 2002, Phys. Rep., 372, 1

* DESI Collaboration, 2016, preprint (arXiv:1611.00036)

* Darwish O. et al., 2021, MNRAS, 500, 2250

* Doré O. et al., 2014, preprint (arXiv:1412.4872)

* Duffy A. R., Schaye J., Kay S. T., Dalla Vecchia C., 2008, MNRAS, 390, L64

* Dutcher D. et al., 2021, Phys. Rev. D, 104, 22003

* Eifler T., Krause E., Schneider P., Honscheid K., 2014, MNRAS, 440, 1379

* Eifler T. et al., 2021, MNRAS, 507, 1514

* Fang X., Eifler T., Krause E., 2020a, MNRAS, 497, 2699

* Fang X., Krause E., Eifler T., MacCrann N., 2020b, J. Cosmol. Astropart. Phys., 2020, 10

* Fang X., Eifler T., Schaan E., Huang H.-J., Krause E., Ferraro S., 2022, MNRAS, 509, 5721

* Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306

* Fortuna M. C., Hoekstra H., Joachimi B., Johnston H., Chisari N. E., Georgiou C., Mahony C., 2021, MNRAS, 501, 2983

* Galitzki N. et al., 2018, in Zmuidzinas J., Gao J.-R. eds, Proc. SPIE Conf. Ser. Vol. 10708, Millimeter, Submillimeter, and Far-Infrared Detectors and Instrumentation for Astronomy IX. SPIE, Bellingham. p. 1070804

* Gatti M. et al., 2022, Phys. Rev. D, 105, 123525

* Giannantonio T. et al., 2016, MNRAS, 456, 3213

* Hang Q., Alam S., Peacock J. A., Cai Y.-C., 2021, MNRAS, 501, 1481

* Heymans C. et al., 2021, A&A, 646, A140

* Hikage C. et al., 2019, PASJ, 71, 43

* Hildebrandt H. et al., 2017, MNRAS, 465, 1454

* Hill J. C., Pajer E., 2013, Phys. Rev. D, 88, 63526

* Hill J. C., Spergel D. N., 2014, J. Cosmol. Astropart. Phys., 2014, 30

* Hill J. C., Baxter E. J., Lidz A., Greco J. P., Jain B., 2018, Phys. Rev. D, 97, 83501

* Hirata C. M., Seljak U., 2004, Phys. Rev. D, 70, 63526

* Hu W., Jain B., 2004, Phys. Rev. D, 70, 43009

* Huang H.-J., Eifler T., Mandelbaum R., Dodelson S., 2019, MNRAS, 488, 1652

* Huang H.-J. et al., 2021, MNRAS, 502, 6010

* Ivezić Ž. et al., 2019, ApJ, 873, 111

* Joudaki S. et al., 2018, MNRAS, 474, 4894

* Kitanidis E., White M., 2021, MNRAS, 501, 6181

* Komatsu E., Seljak U., 2002, MNRAS, 336, 1256
Koukoufilippas N., Alonso D., Bilicki M., Peacock J. A., 2020, MNRAS, 491, 5464

Krause E., Eifler T., 2017, MNRAS, 470, 2100

Krause E., Eifler T., Blazek J., 2016, MNRAS, 456, 207

Krolewski A., Ferraro S., Schlafly E. F., White M., 2020, J. Cosmol. Astropart. Phys., 2020, 47

Krolewski A., Ferraro S., White M., 2021, J. Cosmol. Astropart. Phys., 2021, 28

Laureijs R. et al., 2011, preprint (arXiv:1110.3193)

Lewis A., Challinor A., 2006, Phys. Rep., 429, 1

Limber D. N., 1953, ApJ, 117, 134

LoVerde M., Afshordi N., 2008, Phys. Rev. D, 78, 123506

Lochner M. et al., 2018, preprint (arXiv:1812.00515)

Ma C.-P., Fry J. N., 2000, ApJ, 543, 503

Ma Y.-Z., Waerbeke L. V., Hinshaw G., Hojjati A., Scott D., Zuntz J., 2015, J. Cosmol. Astropart. Phys., 2015, 46

McCarthy I. G., Schaye J., Bird S., Le Brun A. M. C., 2017, MNRAS, 465, 2936

Madhavacheril M. S. et al., 2023, preprint (arXiv:2304.05203)

Makiya R., Ando S., Komatsu E., 2018, MNRAS, 480, 3928

Mead A. J., Tröster T., Heymans C., Van Waerbeke L., McCarthy I. G., 2020, A&A, 641, A130

Mead A. J., Frieden S., Tröster T., Heymans C., 2021, MNRAS, 502, 1401

Miyatake H. et al., 2022, Phys. Rev. Lett., 129, 61301

Modi C., White M., Vlah Z., 2017, J. Cosmol. Astropart. Phys., 2017, 9

Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ, 490, 493

Nicola A. et al., 2022, J. Cosmol. Astropart. Phys., 2022, 46

Omori Y. et al., 2019a, Phys. Rev. D, 100, 43501

Omori Y. et al., 2019b, Phys. Rev. D, 100, 43517

Omori Y. et al., 2023, Phys. Rev. D, 107, 23529

Osato K., Takada M., 2021, Phys. Rev. D, 103, 63501

Pandey S., Baxter E. J., Hill J. C., 2020a, Phys. Rev. D, 101, 43525

Pandey S. et al., 2020b, Phys. Rev. D, 102, 123522

Pandey S. et al., 2022a, Phys. Rev. D, 105, 123526

Pandey S. et al., 2022b, Phys. Rev. D, 106, 43520

Parshley S. C. et al., 2018, in Marshall H. K., Spyromilio J. eds, Proc. SPIE Conf. Ser. Vol. 10700, Ground-based and Airborne Telescopes VII. SPIE, Bellingham. p. 1070041

Peacock J. A., Smith R. E., 2000, MNRAS, 318, 1144

Petersen K. B., Pedersen M. S., 2012, The Matrix Cookbook. http://www2.compute.dtu.dk/pubdb/pubs/3274-full.html, (last accessed on 15 March 2023)

Planck Collaboration, 2020, A&A, 641, A6

Qu F. J. et al., 2023, preprint (arXiv:2304.05202)

Robertson N. C. et al., 2021, A&A, 649, A146

Samuroff S. et al., 2019, MNRAS, 489, 5453

Schaan E., Krause E., Eifler T., Doré O., Miyatake H., Rhodes J., Spergel D. N., 2017, Phys. Rev. D, 95, 123512

Schaan E., Ferraro S., Seljak U., 2020, J. Cosmol. Astropart. Phys., 2020, 1

Schaan E. et al., 2021, Phys. Rev. D, 103, 63513

Schmittfull M., Seljak U., 2018, Phys. Rev. D, 97, 123540

Seljak U., 2000, MNRAS, 318, 203

Singh S., Mandelbaum R., Brownstein J. R., 2017, MNRAS, 464, 2120

Takada M., Hu W., 2013, Phys. Rev. D, 87, 123504

Takada M., Jain B., 2009, MNRAS, 395, 2065

Tanimura H. et al., 2020, MNRAS, 491, 2318

The LSST Dark Energy Science Collaboration, 2018, preprint (arXiv:1809.01669)

Tinker J. L., Robertson B. E., Kravtsov A. V., Klypin A., Warren M. S., Yepes G., Gottlöber S., 2010, ApJ, 724, 878

To C. et al., 2021, Phys. Rev. Lett., 126, 141301

Tröster T. et al., 2022, A&A, 660, A27

Vikram V., Lidz A., Jain B., 2017, MNRAS, 467, 2315

Vlah Z., Chisari N. E., Schmidt F., 2020, J. Cosmol. Astropart. Phys., 2020, 25

Wenzl L., Doux C., Heinrich C., Bean R., Jain B., Doré O., Eifler T., Fang X., 2022, MNRAS, 512, 5311

White M. et al., 2022, J. Cosmol. Astropart. Phys., 2022, 7

Xu Z. et al., 2021, Res. Notes Am. Astron. Soc., 5, 100

Yan Z. et al., 2021, A&A, 651, A76

## APPENDIX A: TWO-POINT FUNCTION MODELLING

This appendix summarizes the computation of angular (cross) power spectra for 10 × 2pt analyses. We use capital Roman subscripts to denote observables, $A, B \in \{\delta_g, \kappa_g, \kappa_{CMB}, y\}$. We assume general relativity (GR) and a flat $\Lambda$CDM cosmology throughout.

### A1 Angular power spectra

Adopting the Limber approximation [Limber 1953; but see LoVerde & Afshordi (2008) and Fang et al. (2020b) for potential impact in current and near future surveys], we write the angular power spectrum between redshift bin $i$ of observable A and redshift bin $j$ of observable B at Fourier mode $\ell$ as

$$C_{AB}^{ij}(\ell) = \int d\chi \frac{q_A^j(\chi)q_B^j(\chi)}{\chi^2} P_{AB} \left( \frac{\ell + 1/2}{\chi}, z(\chi) \right), \quad (A1)$$

where $\chi$ is the comoving distance, $P_{AB}(k, z)$ is the 3D cross-probe power spectra, $q_A^j(\chi)$, $q_B^j(\chi)$ are weight functions of the observables A, B given by

$$q_{\delta_g}^j(\chi) = \frac{n_{\text{lens}}^i(z(\chi))}{n_{\text{lens}}^i} \frac{dz}{d\chi}, \quad (A2)$$

$$q_{kg}^j(\chi) = \frac{3H_0^2 \Omega_m}{2c^2} \frac{\chi}{a(\chi)} \int_{\chi_{\text{min}}}^{\chi_{\text{max}}} d\chi' \frac{n_{\text{source}}^i(z(\chi'))}{n_{\text{source}}^i} \frac{dz}{d\chi'} \frac{\chi' - \chi}{\chi}, \quad (A3)$$

$$q_{\kappa_{CMB}}(\chi) = \frac{3H_0^2 \Omega_m}{2c^2} \frac{\chi}{a(\chi)} \frac{\chi^* - \chi}{\chi^*}, \quad (A4)$$

$$q_y(\chi) = \frac{\sigma_T}{m_e c^2} \frac{1}{[a(\chi)]^2}, \quad (A5)$$

where $\chi_{\text{min/max}}^i$ is the minimum/maximum comoving distance of the redshift bin $i$, $a(\chi)$ is the scale factor, $\Omega_m$ is the matter density fraction at present, $H_0$ is the Hubble constant, $c$ is the speed of light, $\chi^*$ is the comoving distance to the surface of last scattering, $\sigma_T$ is the Thomson scattering cross section, and $m_e$ is the electron mass. Note that the weight functions of $y$ and $\kappa_{CMB}$ do not depend on redshift bins.

We relate the 3D cross-power spectra $P_{AB}(k, z) = P_{BA}(k, z)$ to the non-linear matter power spectrum $P_{\delta\delta}(k, z)$, matter-electron pressure power spectrum $P_{\delta\kappa}(k, z)$, and electron pressure power spectrum $P_{\kappa\kappa}(k, z)$, where $\delta$ is the non-linear matter density contrast and $P_e$ is the electron pressure,

$$P_{\delta g}(k, z) = b_g(z) P_{\delta B}(k, z), \quad (A6)$$

$$P_{\kappa g}(k, z) = P_{\kappa CMB}(k, z) = P_{\delta B}(k, z), \quad (A7)$$

$$P_{yB}(k, z) = P_{\kappa B}(k, z), \quad (A8)$$

where we have assumed that the galaxy density contrast is proportional to the non-linear matter density contrast, weighted by an effective galaxy bias parameter $b_g(z)$. We impose conservative scale cuts on $\delta_g$-related probes to ensure the validity of this assumption in our analysis. Modelling beyond linear galaxy bias would be required to robustly extract cosmological information from smaller scales (e.g. Krause & Eifler 2017; Modi, White & Vlah 2017; Pandey et al. 2020b, 2022b; Hang et al. 2021; Kitanidis & White 2021; Krolewski, Ferraro & White 2021). Due to equation (A8), we will use $P_{yB}$ to
represent $P_{PeB}$, which simplifies the subscripts. Thus, $P_{\delta Pe} \equiv P_{\delta y}$ and $P_{PePe} \equiv P_{yy}$.

We calculate the 3D power spectra $P_{AB}$ ($A, B \in \{\delta, y\}$) using a halo model implementation, which is based on the HMx parametrization from Mead et al. (2020).

### A1.1 Halo model

The halo model (Ma & Fry 2000; Peacock & Smith 2000; Seljak 2000; Cooray & Sheth 2002) describes the distribution of matter, and biased tracers, based on the assumption that all matter is distributed within halos. Our implementation of the gas component largely follows the HMx code, while the other ingredients of the halo model are adapted to the existing CosmoLike implementation (Krause & Eifler 2017), including the halo definition, critical density contrast, and halo bias and mass function. While more accurate models and consistent implementations will be required for data analyses, we do not expect these modelling approximations to significantly impact the exploration of information content in different probe combinations presented here.

The 3D power spectrum of two fields A and B can then be written as the sum of a 2h and 1-halo (1h) term,

$$ P_{AB}(k, z) = P_{AB}^{2h}(k, z) + P_{AB}^{1h}(k, z) \\ = P_{\text{lin}}(k, z) I_A^1(k, z) I_B^1(k, z) + I_{AB}^0(k, z), \quad (A9) $$

where $P_{\text{lin}}(k, z)$ is the linear matter power, and

$$ I_A^\alpha(k, z) = \langle b_{h,\alpha}(M, z) \tilde{u}_A(k, M, z) \rangle_M, \quad (A10) $$

$$ I_{AB}^{\alpha}(k, k', z) = \langle b_{h,\alpha}(M, z) \tilde{u}_A(k, M, z) \tilde{u}_B(k', M, z) \rangle_M, \quad (A11) $$

where $b_{h, \alpha}$ is the $\alpha$-th order halo biasing, with $b_{h, 0} = 1$, and where we neglect higher order biasing, i.e. $b_{h, \ge 2} = 0$. The functions $\tilde{u}_X(k, M, z)$ ($X \in \{\delta, y\}$) are related to the Fourier transforms of the radial profiles of matter density $\rho_m(M, r, z)$ and electron pressure $P_e(M, r, z)$ within a halo of mass $M$, and where we have used a shorthand notation for the halo mass function-weighted average of quantity $X$,

$$ \langle X(M) \rangle_M \equiv \int_{M_{\text{min}}}^{M_{\text{max}}} dM \frac{dn(M)}{dM} X(M), \quad (A12) $$

where $M$ is the halo mass and $dn(M)/dM$ is the halo mass function.

This standard halo model calculation leads to unphysical behaviour of the 1h term as $k \to 0$. We follow the treatment in the public halo model code HMCODE-2020 (equation 17 of Mead et al. 2021) for all 1h terms at low-$k$,

$$ P_{AB}^{1h}(k, z) \to P_{AB}^{1h}(k, z) \frac{(k/k_*)^4}{1 + (k/k_*)^4}, \quad (A13) $$

and take $k_*$ as the fitted functional form in their Table 2.

### A2 Systematics

Systematic uncertainties are parametrized through nuisance parameters, whose fiducial values and priors are summarized in Table 2. The nuisance parameters included in our simulated analysis are very similar to those used in the joint Rubin Observatory-SO forecast in Fang et al. (2022) and the Rubin Observatory-Roman Space Telescope forecast in Eifler et al. (2021), except that here we do not use a principal component analysis approach (Huang et al. 2019, 2021) for the baryonic feedback. Instead, we have treated the baryonic feedback in a consistent halo model as described in the previous subsection. The systematic effects are summarized as follows:

#### A2.1 Photometric redshift uncertainties

We parametrize photometric redshift uncertainties by a Gaussian scatter $\sigma_{z,x}$ for lens and source sample each, and a shift parameter $\Delta_{z,x}^i$ for each redshift bin $i$ of the lens and source samples. The binned true redshift distribution $n_x^i(z)$ is related to the binned photometric redshift distribution $n_x^i(z_{ph})$ (equation 1) by

$$ n_x^i(z) = \int_{z_{\text{min},x}}^{z_{\text{max},x}} \frac{dz_{\text{ph}} n_x^i(z_{\text{ph}})}{\sqrt{2\pi}\sigma_{z,x}(1+z_{\text{ph}})} \exp\left[-\frac{(z-z_{\text{ph}}-\Delta_{z,x}^i)^2}{2[\sigma_{z,x}(1+z_{\text{ph}})]^2}\right]. \quad (A14) $$

We set the fiducial values of $\Delta_{z,x}^i$ as zero, $\sigma_{z, \text{lens}}$ as 0.03, and $\sigma_{z, \text{source}}$ as 0.05, following the DESC-SRD. The resulting distributions are shown in Fig. 1. In total, we have 22 photo-z parameters (10 shift parameters and 1 scatter parameter for lens and source sample each). In this paper, these parameters are analytically marginalized over using Gaussian priors, as described in Section 2.3.2.

For the source samples, the $\Delta_{z, \text{source}}^i$ and $\sigma_{z, \text{source}}$ priors for LSST Y6 are chosen to be the same as the Y10 requirements given in DESC-SRD. For the lens samples, we choose these priors to be same as the corresponding source sample priors.<fnref n="8" />

#### A2.2 Linear galaxy bias

We assume one linear bias parameter $b_g^i$ per lens redshift bin, whose fiducial values follow the simple relation as given in the DESC-SRD: $b_g^i = 0.95/G(z^i)$, where $G(z)$ is the growth function. We independently marginalize over the 10 linear bias parameters with a conservative flat prior [0.4,3.0]. We note that this model may be oversimplified for going into $k_{\text{max}} = 0.3h/Mpc$. More complex models may lead to less constraining power of the galaxy survey and enhance the importance of the secondary CMB information.

#### A2.3 Multiplicative shear calibration

We assume one parameter $m^i$ per source redshift bin, which affects $\kappa_g-X$ ($X \in \{\kappa_g, \delta_g, \kappa_{\text{CMB}}, y\}$) via

$$ q_{kg}^i \to (1 + m^i)q_{kg}^i \quad (A15) $$

The total of 10 $m^i$ parameters are independently marginalized over with Gaussian priors shown in Table 2. The priors are chosen to be the same as the LSST Y10 requirements given in DESC-SRD. In this paper, we analytically marginalize these parameters as described in Section 2.3.2.

#### A2.4 IA

We adopt the ‘non-linear linear alignment’ model (Hirata & Seljak 2004; Bridle & King 2007), which considers only the ‘linear’ response of the elliptical (red) galaxies’ shapes to the tidal field sourced by the underlying ‘non-linear’ matter density field. Our implementation follows Krause, Eifler & Blazek (2016) for cosmic shear, Krause & Eifler (2017) for galaxy–galaxy lensing, Fang et al. (2022) for $\kappa_g-\kappa_{\text{CMB}}$ power spectra, and extend it to $\kappa_g-y$ power spectra. Using the notation in Section A1, we can encapsulate the effect as

$$ q_{kg}^j \to q_{kg}^j + q_{j}^i, \quad (A16) $$

<footnote n="8">⁸This choice follows Fang et al. (2022) and is justified in their footnote 4.</footnote>
where

$$q_I^i = -A(z) \frac{n_{\text{source}}^i(z(\chi))}{\bar{n}_{\text{source}}^i} \frac{dz}{d\chi}. \qquad (A17)$$

$A(z)$ is the IA amplitude at a given redshift $z$, computed by

$$A(z) = \frac{C_1 \rho_{\text{cr}}}{G(z)} A_{\text{IA}} \left( \frac{1+z}{1+z_0} \right)^{\eta_{\text{IA}}}, \qquad (A18)$$

with pivot redshift $z_0 = 0.3$, and with $C_1 \rho_{\text{cr}} = 0.0134$ derived from SuperCOSMOS observations (Hirata & Seljak 2004; Bridle & King 2007). The fiducial values and priors for the nuisance parameters $A_{\text{IA}}$ and $\eta_{\text{IA}}$ are given in Table 2. We neglect luminosity dependence and additional uncertainties in the luminosity function, which can be significant and are discussed in Krause, Eifler & Blazek (2016).

Together with the multiplicative shear calibration, $q_{\kappa_g}^i$ is altered as

$$q_{\kappa_g}^i(\chi) \rightarrow (1+m^i)[q_{\kappa_g}^i(\chi) + q_I^i(\chi)]. \qquad (A19)$$

We do not consider higher-order tidal alignment, tidal torquing models (see e.g. Blazek, Vlah & Seljak 2015; Blazek et al. 2019), or effective field theory models (Vlah, Chisari & Schmidt 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.

## APPENDIX B: ANALYTIC COVARIANCES

The Fourier $10 \times 2pt$ covariance matrix $C(C(\ell_1), C(\ell_2))$ includes the Gaussian part $C^G$ (Hu & Jain 2004), the non-Gaussian part from connected 4-point functions (in short, ‘connected non-Gaussian’ or cNG) in the absence of survey window effect $C^{\text{cNG}}$ (e.g. Cooray &

<figure><img src="image_8.png" /><figcaption>**Figure B1.** The $10 \times 2pt$ correlation matrix for LSST Y6 × SO Y5, plotted as $\log_{10}|R|$, where $R = [\text{diag}(C)]^{-1/2}C[\text{diag}(C)]^{-1/2}$. We have highlighted some parts of the matrix to illustrate the correlation structure: (1) depicts the cosmic shear ($\kappa_g-\kappa_{\text{g}}$) correlation matrix, comprised of 55 tomographic combinations of source bins, each with 25 Fourier $\ell$-bins. (11) shows one of the tomographic combinations, and the individual $\ell_1$, $\ell_2$ elements are clearly visible. (2) is the galaxy–galaxy lensing ($\kappa_g-\delta_g$) tomography correlation with (13) being the galaxy–galaxy combinations of the second 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 ($\delta_g-\delta_g$) auto-probe matrix with 10 tomographic bins. (4) corresponds to the galaxy-CMB convergence ($\delta_g-\kappa_{\text{CMB}}$) combinations auto-probe matrix, which is comprised of 10 lens sample redshift bins. (5) is the auto-probe correlation of the $\kappa_g-\kappa_{\text{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, $\kappa_{\text{CMB}}-\kappa_{\text{CMB}}$. (7) is the auto-probe correlation of the galaxy-tSZ ($\delta_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 three of galaxy–galaxy lensing power spectra ($\kappa_g-\delta_g$), all the 10 galaxy clustering ($\delta_g-\delta_g$), and all 10 $\delta_g-y$ spectra, respectively from left to right. (8) is the auto-probe correlation of the galaxy lensing convergence-tSZ ($\kappa_g-y$) combinations, which uses the 10 source redshift bins. (9) is the auto-probe correlation of the CMB lensing convergence-tSZ ($\kappa_{\text{CMB}}-y$). (10) is the auto-probe correlation of the tSZ power spectrum ($y-y$).</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure B2.** The diagonal elements of the tSZ power spectrum covariance. The dashed, dash-dotted, and solid lines show the Gaussian covariance, the SSC, the 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 1h contributions.</figcaption></figure>

Sheth 2002; Takada & Jain 2009), and the super-sample covariance (SSC) $\mathbf{C}^{\text{SSC}}$ (Takada & Hu 2013), i.e. $\mathbf{C} = \mathbf{C}^{\text{G}} + \mathbf{C}^{\text{cNG}} + \mathbf{C}^{\text{SSC}}$.

## B0.5 Gaussian covariance

We extending the modelling and implementation described in Krause & Eifler (2017), Schaan et al. (2017), and Fang et al. (2022) to include the y-related probes. The Gaussian covariance includes the probe-specific shot noise terms, i.e.

$$ \mathbf{C}^{\text{G}}(C_{\text{AB}}^{ij}(\ell_1), C_{\text{CD}}^{kl}(\ell_2)) = \frac{\delta_{\ell_1 \ell_2}}{f_{\text{sky}}(2\ell_1 + 1)\Delta\ell_1} \left[ \hat{C}_{\text{AC}}^{ik}(\ell_1)\hat{C}_{\text{BD}}^{jl}(\ell_2) + \hat{C}_{\text{AD}}^{il}(\ell_1)\hat{C}_{\text{BC}}^{jk}(\ell_2) \right], \quad (\text{B1}) $$

where $\delta_{ij}$ is the Kronecker delta function, $f_{\text{sky}} = \Omega_{\chi}/(4\pi)$ is the survey's sky coverage fraction, $\Delta\ell$ is the $\ell$ bin width, $\hat{C}_{\text{AC}}^{ik}(\ell_1) = C_{\text{AC}}^{ik}(\ell_1) + \delta_{ik}\delta_{\text{AC}}N_{\text{A}}^{i}(\ell_1)$. The noise term $N_{\text{A}}^{i}(\ell)$ are given by

$$ N_{\kappa_g}^i(\ell) = \sigma_{\epsilon}^2 / \bar{n}_{\text{source}}^i, \quad N_{\delta_g}^i(\ell) = 1 / \bar{n}_{\text{lens}}^i, \quad (\text{B2}) $$

where $\sigma_{\epsilon} = 0.26$ is the shape noise per component given by the DESC-SRD. $N_{\kappa_{\text{CMB}}}^i(\ell)$, $N_{\gamma}^i(\ell)$ are the CMB lensing reconstruction noise and the Compton-y reconstruction noise.

The SO noise models are based on component separated maps for the 'goal' sensitivity (SENS-2), and include the expected level of residual foregrounds as well as atmospheric noise. The noise for the Compton-y (tSZ) map,<fnref n="9" /> is obtained from the standard 'Internal Linear Combination' (ILC, deprojo) without any additional deprojection on 40 per cent of the sky. For the noise properties of the CMB lensing map<fnref n="10" /> we take the same sensitivity and area as the tSZ map, with the minimum variance ILC, including both temperature and polarization data (with $\ell_{\text{min}, T} = \ell_{\text{min}, P} = 30$, $\ell_{\text{max}, T} = 3000$, and

$\ell_{\text{max}, P} = 5000)$ and no additional deprojection (deprojo). The lensing noise includes a small improvement expected from iterative lensing reconstruction. These noise curves were obtained with the methodology explained in Section 2 of Ade et al. (2019) and are publicly available on GitHub at the address linked in the footnotes below.

## B0.6 cNG covariance

The cNG term is computed from the projection of the trispectra. We neglect the mode-coupling due to the survey window and adopt the Limber approximation,

$$ \mathbf{C}^{\text{cNG}}(C_{\text{AB}}^{ij}(\ell_1), C_{\text{CD}}^{kl}(\ell_2)) = \frac{1}{4\pi f_{\text{sky}}} \int_{\ell_1}^{\ell_2} \frac{d^2\ell}{A(\ell_1)} \int_{\ell_2}^{\ell_1} \frac{d^2\ell'}{A(\ell_2)} \\ \times \int d\chi \frac{q_A^i(\chi)q_B^i(\chi)q_C^i(\chi)q_D^i(\chi)}{\chi^6} T_{\text{ABCD}}^{ijkl} \\ \times \left( \frac{\ell}{\chi}, -\frac{\ell}{\chi}, \frac{\ell'}{\chi}, -\frac{\ell'}{\chi}; z(\chi) \right). \quad (\text{B3}) $$

In the case where none of A,B,C,D fields is y field, we approximate the trispectrum as

$$ T_{\text{ABCD}}^{ijkl}(\mathbf{k}, -\mathbf{k}, \mathbf{k}', -\mathbf{k}'; z) \approx T_{\text{ABCD}}^{1h}(\mathbf{k}, \mathbf{k}, \mathbf{k}', \mathbf{k}'; z) \\ + b_{\text{A}}^{i}(z)b_{\text{B}}^{i}(z)b_{\text{C}}^{k}(z)b_{\text{D}}^{l}(z)T_{\text{m}}^{2h+3h+4h} \\ \times (\mathbf{k}, -\mathbf{k}, \mathbf{k}', -\mathbf{k}'; z), \quad (\text{B4}) $$

where the bias $b_{\text{A}}^{i}(z)$ is equal to the galaxy bias $b_{\text{g}}^{i}(z)$ if $A = \delta_{\text{g}}$, and $b_{\text{A}}^{i}(z) = 1$ otherwise. For (cross-)covariances involving at least one y field, we only consider the 1h contribution (Komatsu & Seljak 2002; Makiya, Ando & Komatsu 2018)

$$ T_{\text{ABCY}}(\mathbf{k}, -\mathbf{k}, \mathbf{k}', -\mathbf{k}'; z) = \langle \tilde{u}_{\text{A}}(\mathbf{k}, \mathbf{M}, z) \tilde{u}_{\text{B}}(\mathbf{k}, \mathbf{M}, z) \rangle \\ \times \tilde{u}_{\text{C}}(\mathbf{k}', \mathbf{M}, z) \tilde{u}_{\text{Y}}(\mathbf{k}', \mathbf{M}, z) \rangle_{\text{M}}. \quad (\text{B5}) $$

## B0.7 SSC

We extend the CosmoCov SSC implementation to include y cross-correlations, following the formalism of (Takada & Hu 2013) applied to multiprobe observables (Krause & Eifler 2017) and assuming a polar cap survey footprint to evaluate the variance of super-survey density modes.

We show the full correlation matrix $\mathbf{R}$, given by $\mathbf{R} = [\text{diag}(\mathbf{C})]^{-1/2}\mathbf{C}[\text{diag}(\mathbf{C})]^{-1/2}$, of the $10 \times 2pt$ data vector in Fig. B1, evaluated at the fiducial parameter values. The y-y and shear-y probes are highly correlated across tomographic bins and $\ell$-bins due to the significant contributions from the cNG term. In Fig. B2, we show the Gaussian, cNG, and SSC components of the diagonal elements of the $\mathbf{C}(C_{yy}, C_{yy})$ covariance block. The dominance of the cNG component over the SSC and the Gaussian part is consistent with Osato & Takada (2021). As shown in their paper, applying a cluster mask for massive clusters can significantly suppress the noise level. We leave the studies of its improvement of our $10 \times 2pt$ analysis for future work.

This paper has been typeset from a TeX/LA TeX file prepared by the author.

<footnote n="9">⁹https://github.com/simonsobs/so-noise_models/blob/master/LAT_comp.s ep noise/v3.1.0/SO_LAT_Nell_T_atmv1_goal_fsky0p4_LIC_tSZ.txt.</footnote>

<footnote n="10">¹⁰https://github.com/simonsobs/so-noise_models/blob/master/LAT lensing _noise/lensing_v3.1.1/nlkk_v3.1.0.deprojo_SENS2_fsky0p4_it_1T30-3000_1 P30-5000.dat.</footnote>