12 × 2 pt combined probes: pipeline, neutrino mass, and data compression

With the rapid advance of wide-field surveys it is increasingly important to perform combined cosmological probe analyses. We present a new pipeline for simulation-based multi-probe analyses, which combines tomographic large-scale structure (LSS) probes (weak lensing and galaxy clustering) with cosmic microwave background (CMB) primary and lensing data. These are combined at the C ℓ-level, yielding 12 distinct auto- and cross-correlations. The pipeline is based on UFalconv2, a framework to generate fast, self-consistent map-level realizations of cosmological probes from input lightcones, which is applied to the CosmoGridV1 N-body simulation suite. It includes a non-Gaussian simulation-based covariance for the LSS tracers, several data compression schemes, and a neural network emulator for accelerated theoretical predictions. We validate the pipeline by comparing the simulations to these predictions, and our derived constraints to earlier analyses. We apply our framework to a simulated 12×2 pt tomographic analysis of KiDS, BOSS, and Planck, and forecast constraints for a ΛCDM model with a variable neutrino mass. We find that, while the neutrino mass constraints are driven by the CMB data, the addition of LSS data helps to break degeneracies and improves the constraint by up to 35%. For a fiducial Mν = 0.15 eV, a full combination of the above CMB+LSS data would enable a 3σ constraint on the neutrino mass. We explore data compression schemes and find that MOPED outperforms PCA and is made robust using the derivatives afforded by our automatically differentiable emulator. We also study the impact of an internal lensing tension in the CMB data, parametrized by AL , on the neutrino mass constraint, finding that the addition of LSS to CMB data including all cross-correlations is able to mitigate the impact of this systematic. UFalconv2 and a MOPED compressed Planck CMB primary + CMB lensing likelihood are made publicly available.[UFalconv2: https://cosmology.ethz.ch/research/software-lab/UFalcon.html, compressed Planck CMB primary + CMB lensing likelihood: https://github.com/alexreevesy/planck_compressed.]

1 Introduction With the onset of next-generation surveys such as Euclid [1], the Large Synoptic Survey Telescope (LSST) [2], the Simons Observatory [3], Roman Space Telescope [4] and CMB-S4 [5], we will be uniquely positioned to rigorously test our cosmological model.In this context of high-precision data, increasingly dominated by systematics, combined probes analyses offer several advantages.This method uses a common framework to analyze different survey datasets, allowing for a consistent exploration of discrepancies between them [6,7].In addition, combining cosmological probes can break parameter degeneracies inherent in single probe analyses, thus leading to tighter constraints [8,9].Finally, by leveraging cross-correlations, we can calibrate systematic effects that do not correlate between the datasets [10,11].It is no surprise, then, that combined probe analyses have been a central consideration in the planning of recent and future surveys.In the context of large-scale structure (LSS) data, 3 × 2pt analyses, which combine weak lensing (WL) and galaxy clustering data at the 2-point level, are now standard for current and recent optical surveys such as the Dark Energy Survey (DES) [12], Hyper Suprime Cam (HSC) [13] and the Kilo-Degree Survey (KiDS) [14,15].Cosmic microwave background (CMB) measurements have also been jointly analyzed with LSS data to derive strong constraints: There have been a number of recent studies considering the CMB lensing-galaxy clustering cross-correlation [16][17][18][19][20][21], the Integrated Sachs-Wolfe (ISW) cross-correlation [22,23] and the CMB lensing-galaxy weak lensing cross-correlation [24].
The aspiration to create a comprehensive framework for analyzing all datasets, from both the CMB and LSS realms, remains a fundamental goal within combined probes analysis.In recent years a number of studies have emerged that move towards this goal [25][26][27][28][29][30].In particular, Ref. [26] builds upon the map-based approach introduced in Refs.[28,29] to combine CMB and LSS data.Recently, Ref. [25] (see also [31,32]) provides a halo modelbased scheme to combine galaxy clustering, weak lensing, CMB lensing, and thermal Sunyaev Zel'dovich (tSZ) data, resulting in 10 different 2-point functions.
In this work, we introduce a novel combined probes pipeline, designed to incorporate several cosmological datasets at the C ℓ level, inclusive of their respective cross-correlations.Specifically, this pipeline is designed to incorporate CMB primary, galaxy clustering, weak lensing (WL), and CMB lensing data, resulting in a suite of 12 distinct 2-point (2pt) functions once cross-correlations are included.The pipeline includes several advantages compared to previous approaches including a non-Gaussian simulation-based covariance for the LSS probes, methods for data compression, and a neural network emulator for accelerated parameter exploration.The pipeline uses UFalconv2, a new publicly available software that produces consistent map-level realizations of LSS probes (WL, galaxy clustering CMB lensing, and ISW) from an input lightcone.This is applied to the CosmoGridV1 N-body simulation suite to produce a set of correlated map-level realizations from which the LSS mock data vector and covariance matrix are derived.
We apply this framework to perform a forecast tomographic analysis of KiDS [14,15], BOSS [33], and Planck [34].We validate our pipeline by comparing our constraints with those of previous analyses and ensuring that we recover unbiased constraints.We then examine two applications of the framework.In the first, we use the full 12×2pt data vector and forecast its constraining power on the sum of the neutrino mass eigenstates, M ν = i m ν i .We also analyze how the 'A-lens' parameter (A L ), first introduced in [35], impacts this measurement 1 .
In a second application, we explore methods of data compression.This is especially important in the case of multiprobe analyses where uncompressed data vectors are typically large.We specifically compare the performance of two data compression techniques: MOPED [37] and Principal Component Analysis (PCA) [38].We also apply MOPED to the full Planck CMB primary and lensing likelihood, extending the work presented in Ref. [39].

Overview of pipeline
A schematic of our combined probes pipeline is shown in Fig. 1.On the left-hand side, the framework is built upon UFalconv2 which processes N-body simulations into correlated map-level realizations of each probe.These are used to produce both the fully non-Gaussian covariance matrix (bottom left panel in Fig. 1) and in producing realistic mock data vectors (middle panel in Fig. 1).We then optionally compress these data vectors, combine them with theoretical predictions obtained from a differentiable neural-network emulator (bottom right panel in Fig. 1), and eventually derive the likelihood and cosmological parameter constraints (respectively bottom two panels in Fig. 1).In this way, the pipeline is built in a modular fashion and the interested reader can skip to one of the descriptions in the sections shown in Fig. 1 for details.
The paper is structured as follows: We begin with a description of the suite of simulations in § 2, including links to the publicly available UFalconv2 software.Then, in § 3 we describe the theoretical modeling of probes in the pipeline and the training of our multiprobe neural network emulator.In § 4, we outline the data products from real surveys that are used to build our realistic mock observations, and in § 5 how the multiprobe likelihood is computed and interfaced with the inference part of the pipeline.Finally, in § 6 we describe and discuss the results of our analysis before concluding in § 7.
Figure 1: An overview of the combined-probes pipeline.The left-hand channel represents the processing of simulations in this pipeline to produce a Non-Gaussian covariance, whilst the right-hand channel represents the process of training a neural network emulator for rapid theoretical predictions.Finally, the middle section shows how survey data is used to derive a data vector (in this case a mock data vector) before all three elements are combined and optionally compressed (pink panel) to compute the likelihood and eventually derive parameter constraints.

CosmoGridV1
In this section, we describe the simulations that we use for both estimating the covariance matrix and producing realistic mock data vectors.We use the simulations provided in the CosmoGridV12 [40].This is a large suite of simulations produced with the GPU-accelerated N-Body code PkDGraV3 [41].The simulations have a box size of L = 900 Mpc/h, contain 8323 particles to sample the dark matter distribution, and are stored as lightcone output.The lightcone consists of 69 shells stored as HEALPix 3 maps with nside=2048, in the redshift range z = 0 − 3.5, which allows for accurate map projections for various cosmological probes.Whilst the full CosmoGridV1 consists of a suite of 2500 cosmologies with varied [Ω m , σ 8 , h, w 0 , n s , Ω b ], in this pipeline we make use only of the 200 fiducial simulations, run at a fixed set of cosmological parameters shown in Table 1.The effect of massive neutrinos is included in the simulations following the implementation from Ref. [42] available in PkDGraV3.This method accurately models the effect of massive, but light neutrinos, on the matter clustering, using linear perturbation theory (see Ref. [40]

Map making with UFalconv2
To process the simulations described above into cosmological signal maps we created UFalconv2 (Ultra Fast lighcone), an updated version of the UFalcon code, first presented in [43,44].
The source code and documentation (including usage and installation instructions) for the package are available at the following webpage: https://cosmology.ethz.ch/research/software-lab/UFalcon.html.The updated code is configured to output maps of galaxy weak lensing, galaxy clustering, CMB lensing, and the ISW temperature perturbation.Whilst UFalcon was originally configured to process N-body code snapshot output, the updated version directly handles lightcone output, enabling processing of the CosmoGridV1.Processing a CosmoGridV1 lightcone into any of the signal maps takes O(∼ minutes) on a Mac M1 machine.
In the following subsections, we briefly describe the procedures implemented in UFalconv2 to produce full-sky realizations of each of the above probes.We produce all maps in the pipeline at a fiducial HEALPix [45] resolution of nside=1024.Some example maps generated using UFalconv2 are shown in Fig. 2.

Galaxy weak lensing
Galaxy weak lensing (WL) refers to the distortion of the observed shapes of galaxies due to the gravitational deflection of light as it passes through intervening mass distributions.We can indirectly measure the dark matter distribution via analysis of the coherent alignment of galaxy shapes due to gravitational lensing (see [46,47] for reviews).In this context, galaxy shear refers to the differential stretching of the observed shape of a galaxy, while convergence describes its dilation or contraction.In the Born approximation, which ignores the fact that light rays do not follow straight lines, the convergence κ can be related to the line-of-sight projection of the dark matter field δ with a kernel as [47]: where Ω m is the matter density parameter, z s represents the limiting redshift of the galaxy sample, H 0 is the Hubble constant, E(z) = H(z)/H 0 is the dimensionless Hubble parameter, c is the speed of light, and g(z, z s ) is the galaxy lensing radial function.The lensing radial function is defined as [47]: where D(z) = H 0 c χ(z) is the dimensionless comoving distance at redshift z, D(z, z ′ ) = D(z ′ )− D(z) and n(z) is the redshift distribution of source galaxies normalized such that n(z ′ )dz ′ = 1.
To compute a map from the CosmoGridV1 lightcones, we approximate the integral in equation 2.1 as a sum over the redshift shells that comprise a lightcone.Using the relative number of simulation particles in a given pixel as a proxy for the underlying dark matter contrast, the convergence field can be written (up to an overall constant) as [43,48]: In this equation, N pix is the total number of pixels in the map, n p is the number of particles at the pixel position θ pix , N p is the total number of particles, V sim is the volume of the simulation box and ∆z b represents the redshift width of the shell.Note that the user must specify the input cosmological parameters and the redshift distribution function n(z) to compute these maps.
Using equation 2.3, we are able to compute weak lensing convergence maps from an N-body simulation lightcone.However, in survey data, we measure the spin-2 weak lensing shear field by analyzing galaxy shapes.We can relate this to the convergence field via [49]: where 2 γ ℓm are the spin-2 spherical harmonic coefficients of the weak lensing shear observable.This relation is integrated into the UFalconv2 code and the user can specify if they want convergence or shear output.

Galaxy clustering
Galaxy clustering refers to the non-random spatial distribution of galaxies in the universe which act as biased tracers of the underlying dark matter field (see [50] for a recent review).In UFalconv2 galaxy clustering is modeled using the linear bias approach such that the galaxy over-density δ g is related to the underlying dark matter over-density δ via: We can compute the projected galaxy field as an integral over the dark matter distribution: where n g (z) represents the normalized redshift distribution of the galaxy survey.We can then make an analogy to equation 2.3, and approximate this integral as a weighted sum over lightcone slices: . (2.7)

CMB weak lensing
Weak gravitational lensing of the CMB refers to the deflection of CMB photons as they travel through the large-scale structure of the Universe.This induces a signature correlation in CMB temperature and polarization measurements which can be used to reconstruct the lensing potential.This is a probe of the intervening dark matter between the surface of last scattering where the CMB photons are emitted and the observer (see e.g.[51] for a review).
We implement the CMB weak lensing probe in UFalconv2 in an analogous manner to galaxy weak lensing.In this case, the redshift distribution of the source can be effectively modeled as a Dirac delta-function n(z) = δ D (z − z * ) where z * ∼ 1100 is the redshift of recombination, and δ D (z) is the Dirac delta function.Inserting n(z) = δ D (z − z * ) in the galaxy weak lensing window function of equation 2.3 and integrating this out leads to a simplified expression for the CMB lensing window function: . ( N-body simulation lightcones typically contain shells up to a limiting redshift z max .This means the CMB lensing maps produced using UFalcon2 include only the contribution from 0 < z < z max .For sufficiently high z max we can model the remaining CMB weak lensing signal as a Gaussian random field.We can then supplement the UFalconv2-derived CMB lensing map with a Gaussian realization for the z max < z < z * part (note we assume the radial correlation between this additional contribution and the 0 < z < z max piece is negligible)4 .

Integrated Sachs-Wolfe Effect
The ISW effect [53,54] arises from the time evolution of the gravitational potential, which leads to an imprint on the CMB at large physical scales (ℓ < 200).The ISW signal can be detected via cross-correlations with large-scale structure tracers and provides a powerful probe of the late-time dark energy evolution at low redshifts.
The previous computation of the ISW signal in UFalcon relied on snapshots of N-body simulations to compute the time dependence of the gravitational potentials [44].In the new release, we provide an alternative approach that works on the lightcone, enabling crossprobe compatibility and reduced computational requirements.Our pipeline closely follows the method presented in Ref. [55].
The ISW effect leads to a change in the CMB photon temperature according to: with r * the comoving distance to the surface of last scattering.We assume no anisotropic stress and ignore non-linear effects such as the Rees-Sciama effect [54], which have a negligible impact at the low redshifts and large scales analyzed in this work [56].In the linear regime, the time derivative of the potential Φ can be expressed in terms of the Hubble parameter and the linear growth factor: with f (a) = d(ln(D(a))) d(ln(a)) the linear growth rate, D(a) the linear growth factor, and H(a) the Hubble parameter, all computed at scale factor a.
As a final step to enable computations on the lightcone, we rewrite Equation 2.10 in terms of matter constrast δ instead of the potential Φ. Motivated by Poisson's equation, we apply a Spherical Bessel Transform (SBT) decomposition of Φ and δ, since the SBT basis functions are the eigenfunctions of the Laplacian operator on a spherical manifold.The SBT decomposition of a real-valued field in 3D is given by: with k ℓn = q ℓn /r(z max ), z max the highest redshift used in the analysis, q ℓn the n th root of j ℓ (x), and N ℓn = [r(z max )] 3 • [j ℓ+1 (k ℓn • r(z max ))] 2 /2 a normalization factor.We apply the SBT decomposition to Poisson's equation in the linear regime, to obtain a relationship between the coefficients Φ ℓmn and δ ℓmn .Furthermore, we decompose the ISW temperature anisotropies from equation 2.10 into spherical harmonics.Combining both of these equations, we obtain the final result: We note that this SBT algorithm works within a range of physical scales (denoted by k min and k max in the pipeline implementation).These can be overwritten by the user, but are automatically set to reasonable default values by UFalconv2 (related to the physical size of the simulation box and the non-linear physical scale at the relevant redshifts under analysis).

Theory predictions
In this section, we describe the theoretical predictions in our pipeline, which form a key component of the eventual parameter estimation.We begin by describing how the different angular power spectra C ℓ s are modeled theoretically in § 3.1, then how these predictions are computed by a theory code in § 3.2, before finally describing how we build a neural network emulator to rapidly output the predictions for each probe in § 3.3.For the clustering map we use the redshift distribution from the BOSS CMASS selection [33] and a corresponding linear bias value of b = 2.086, and for the shear map we use the first redshift bin from KiDS-1000 [57].All maps plotted above are dimensionless.

Modeling
The theoretical modeling of two-point functions in this pipeline generally separates into two regimes: the LSS probes (WL, clustering, and their respective cross-correlations) and the CMB probes (primary temperature/polarization and CMB lensing).
To compute the LSS correlations we make use of the Limber approximation [58-61] which simplifies the computation of the C ℓ s.This is generally accurate for ℓ ≳ 20 [62].In this approximation, the angular power spectrum can be written as: where P (k) is the matter power spectrum, χ is the comoving distance, z(χ) is the redshift and W X/Y are window functions specific to the probes in the pipeline (X, Y ∈ {γ, δ g , κ CM B , ∆T ISW }).
In the following subsection, we describe the window functions for the different LSS C ℓ s computed in the pipeline as well as any systematics nuisance parameters used in their modeling.

Galaxy weak lensing
The window function for the WL shear signal can be derived by considering the deflection of light rays due to the presence of a gravitational potential [47].In the Born approximation, this takes the following form (see e.g.[63]): where n(χ) represents the weak lensing photometric redshift distribution and χ is the comoving distance5 .We add an intrinsic alignment (IA) component to the pure WL cosmological signal.This is a correlation in the shape of galaxies, in addition to any correlation induced by gravitational lensing, due to local tidal forces which tend to align galaxies (see [64][65][66] for reviews).The total observed signal is computed as the sum of a cosmological shear part C gg ℓ , an IA piece C II ℓ and their cross-correlation C gI ℓ .We use the non-linear alignment model (NLA) to model the IA contribution [67], for which the auto-ad cross-correlation contributions can be computed in the Limber framework with the following window function: where A IA is the NLA amplitude parameter which we assume to be redshift-independent in this work (see [57]).We also include a ∆ z parameter for each WL redshift bin which can shift the photometric redshift distribution to account for redshift miscalibration.

Galaxy clustering
We model galaxy clustering using the scale-independent linear bias approach [68,69].The window function for this model reads: where b(χ) is the redshift dependent linear bias.We allow for one free bias parameter per redshift bin in our analysis.We do not include any parameters to account for redshift distribution systematics.These are assumed to be negligible given we model spectroscopic surveys for galaxy clustering in this pipeline.Furthermore, we do not include the effect of redshift-space distortions (RSDs).

CMB lensing cross-correlations
The CMB lensing probe has a window function that can be derived analogously to that of the WL case but with a single source redshift at z = z * .The window function thus takes the following form (see e.g.[21]): where χ * represents the comoving distance to the last scattering surface.We note that the Limber approximation is used only in the prediction of the CMB-lensing cross-correlations in the pipeline whilst the auto-correlation is computed without making this approximation.

ISW cross-correlations
For the ISW cross-correlation C ℓ s, we make an analogy to the Limber approximation following [29]: ) where z * is the redshift at recombination and W X is as before the window function of an LSS probe in the pipeline.

CMB TTTEEE and lensing auto-correlation
For the CMB probes (C T T ℓ , C T E ℓ , C EE ℓ and C κκ ℓ ), we instead compute the angular power spectra using the CLASS code directly [70].This computation does not rely on the Limber approximation and hence we can safely include larger scales (ℓ < 20) for these spectra in our analysis.

Theory code
There exist many available codes to compute the non-linear matter power spectrum required for the Limber integrals, and the computation of the CMB C ℓ s given a set of cosmological parameters, ⃗ θ (see e.g.CLASS [70] and PyCosmo [63,71,72]).In our pipeline, we adopt a hybrid approach.We use CLASS as our primary engine to compute the matter power spectrum and the CMB C ℓ s.We use the native implementation of HMcode-2016 [73] to compute the non-linear contribution to the matter power spectrum.This spectrum is then passed into PyCosmo, which handles the Limber integrals for the LSS probes in the pipeline [63].

Emulator
Emulators are neural networks that approximate the output of computationally intensive functions.This is done by representing the mapping from input parameters to function output as a series of non-linear transfer function operations.These are increasingly used in cosmology due to the significant speed gains over traditional methods (see [74][75][76][77][78][79][80][81]).Another advantage is that emulators are by default fully differentiable enabling robust numerical predictions of dC ℓ /d ⃗ θ using automatic differentiation (AD).
In this pipeline, we train emulators to learn the mapping between input cosmological parameters and the output spectra of the theory code described in § 3.2.There are 43 distinct spectra to emulate (including all possible auto-and cross-correlation C ℓ s across all redshift bins of the LSS surveys) which can be interpreted as distinct outputs per input parameter set.Instead of training a separate network for each spectrum, we managed to train just six models to cover the entire range of C ℓ s whereby some of the networks are made to output concatenated sets of multiple spectra.This number was found as a compromise between the implementation complexity of training many separate models and the accuracy requirement for each emulated spectrum (see § B for further details).Although each network in the pipeline is created separately, we used a common training data generation framework and identical network architecture for each.
Training data To produce the training data we sample roughly half a million points in parameter space using a Latin Hypercube (LHC) sampling.We train over the parameters of the νΛCDM model and nuisance parameters associated with the target C ℓ (see § B.1 for details on the exact parameters and priors used for each network).At each sampled point of this parameter space, we use our theory code to compute the relevant spectra (the target output of the network).Note we fix the redshift distributions to the respective distributions of the surveys we model for the galaxy weak lensing and clustering spectra (see description of surveys in § 4).Instead of training the models directly on the raw C ℓ , we further compute log(C ℓ ) where possible.This helps to reduce the dynamic range of the training data and vastly improves the emulator fidelity.In the case of C T E ℓ , taking the logarithm is not possible as this spectrum contains negative values across the entire parameter space.In this case, we follow the approach adopted in Ref. [79]: we train on the 512 highest variance PCA modes of the TE training data which similarly helps to reduce the dynamic range and ameliorates the final model accuracy.For C κκ ℓ , we follow Ref. [79] in also training on 512 PCA components as we similarly find this approach gives slightly better performance compared to training on log(C κκ ℓ ).For a subset of the LSS C ℓ s, there are also negative values in the training data (for example WL cross-correlation C ℓ s with large IA amplitudes).In these corner cases, we add an "offset" vector, ⃗ b, to these spectra to enforce positivity before taking the logarithm.We choose the offset vector to be twice the absolute value of the minimum of the concatenated vector of C ℓ s across the emulator training set i.e. ⃗ b = 2 • | min(C ℓ )|.We note that we choose ω cdm as opposed to Ω m in our parameterization of the cosmological model as our initial attempt with the latter consistently gave insufficient accuracy for the C T E ℓ spectrum. 6rchitecture For each emulator in this pipeline, we use the architecture outlined in Ref. [79] 7 .This consists of a neural network, implemented in the Tensorflow framework [82], with 4 hidden layers of 512 neurons each.The activation function for each layer with input ⃗ x is: where the free parameters ⃗ β and ⃗ γ are trained in addition to the weights in each layer.To train the models, we use the gradient descent optimizer ADAM [83] with default parameters.We use a batch size of 1024 and train for ∼ 200 epochs which takes ∼ 3 hours on a Mac M1 machine.We check across 60000 test points that the median relative absolute deviation between the emulated and true spectra is less than 0.1% for all of our models.For further details and emulator validation, the reader is referred to Appendix B.

Survey modeling
To obtain forecasted constraints for mock data we now choose three specific surveys to model.Specifically, we include mock observations of KiDS-1000 WL data [15,57], BOSS galaxy clustering data [33] and Planck CMB data (both primary TTTEEE and reconstructed lensing) [85,86].
We analyze the LSS probes and their cross-correlations using map-level mock observations.In the following, we explain in turn how we create noisy, masked mock observations for each of the LSS surveys.We then describe how we measure the auto-and cross-correlation power spectra and account for the survey mask.Finally, we describe the scale-cuts chosen for this analysis.
and C κκ ℓ are instead included in the framework at the likelihood level using official data products (e.g.covariance matrices).We, therefore, delay the description of these until we describe our combined likelihood in § 5.

Coordinate system
Given that the LSS data from BOSS DR12 and KiDS-1000 is supplied in celestial coordinates, whilst Planck data is in galactic coordinates we must choose a common coordinate choice for each data set.We choose to place all data in galactic coordinates and achieve this by a pre-processing step using the rotator method from healpy to change the coordinates of the observed galaxy positions in the LSS surveys at the catalog level.

KiDS-1000
We create mock observations to mimic the fourth data release [87] of the Kilo-Degree Survey (KiDS).This is a public survey run by the European Southern Observatory (ESO) 8 , which aims to measure the correlations in distortion of galaxy shapes due to weak gravitational lensing.The photometric redshift of each galaxy is measured using 9 bands including the optical ugri as well as the near-infrared ZYJHK s .The addition of the near-infrared bands is possible thanks to the KiDS partner survey VIKING (the VISTA Kilo-degree INfrared Galaxy survey [88]).The data in the fourth release cover 1006 deg 2 on the sky and contain ∼ 2×10 7 individual galaxies across all of the photometric redshift bins for which shape measurements are made [89].The fiducial cosmological parameter constraints from 2-point statistics for this data release are presented in [57], and there have since been many studies including different (higher-order) statistics as well as cross-correlations with other surveys (see [14,90,91]).For the present analysis, we focus on the two-point C ℓ statistic and follow a similar methodology to [91] in analyzing the data using the pseudo-C ℓ approach.In particular, we subdivide the survey galaxies into 5 distinct redshift bins and use the associated photometric redshift distribution functions [84] (γ 1 to γ 5 in Fig. 3).
To produce map-level mock observations of the KiDS-1000 survey, we first use UFalconv2 to produce full-sky cosmological signal maps for each photometric redshift bin.We then add a noise realization generated by randomly rotating the shapes of galaxies in the catalog with a shift ϕ ∈ [0, 2π) but keeping their positions on the sky fixed.This cancels the correlation between galaxy shapes due to lensing and leaves a random shape noise realization.Finally, we create a survey mask using the 'lensing weights' derived from lensfit following Ref.[91] 9 .For this, we bin the lensfit weights into pixels of a HEALPix map to produce a mask that weights the shear pixel value according to the shape measurement uncertainty.

BOSS DR12
We include a mock data analysis of the BOSS 12 th and final public data release10 (hereafter BOSS DR12).BOSS is part of the Sloan Digital Sky Survey (SDSS) III project which took measurements using a 2.5-meter optical imaging telescope.This data release consists of two spectroscopic galaxy catalogs (LOWZ and CMASS).The galaxies within these samples are targeted by applying color-color and color-magnitude cuts to the SDSS photometric catalog [92,93].
We use two broad redshift bins of 0 < z < 0.4 and 0.4 < z < 0.8 corresponding to the LOWZ and CMASS cuts in the BOSS fiducial analysis (see Fig. 3).We produce galaxy over-density maps using UFalconv2 corresponding to the redshift distributions for each cut.To create noisy mock realizations, we then add a random shot-noise realization produced by randomly shuffling the positions of the BOSS DR12 catalog galaxies within the survey mask following the procedure outlined in Ref. [29].
We compute the survey mask from the acceptance and veto masks provided by the BOSS collaboration.These are made available in the MANGLE format [94,95].The acceptance mask for a given region is defined as: i.e. the ratio of the number of observed galaxies with measured redshift to the number of targeted galaxies.The veto masks are instead binary maps that mask areas affected by observational factors [92].We follow the approach taken in Ref. [96] in creating a binary mask from these: first, we convert the acceptance mask to a high resolution nside=16384 HEALPix map and subsequently apply the binary veto masks at this resolution.We then downgrade this combined map to the resolution used in our analysis (nside=1024).Finally, we produce a binary mask by assigning 1 to pixels where the resultant low-resolution map is greater than 0.75 and 0 otherwise following Ref.[19].

Planck CMB lensing cross-correlations
The Planck project, a public survey from the European Space Agency, has provided unprecedented full-sky observations of the CMB and its anisotropies [34].In this pipeline, we create mock observations to simulate the third and final full data release (PR3) [97].The gravitational lensing of the CMB, caused by the deflection of the CMB photons by large-scale structure at late-times, was detected by the Planck team at over 40σ [86].
In this pipeline, the CMB primary C T T ℓ , C T E ℓ and, C EE ℓ as well as the lensing autocorrelation C κκ ℓ are modeled at the likelihood level (see § 5.1.1).However, for the crosscorrelations of Planck CMB lensing, we follow a procedure similar to the LSS probes whereby we create map-level mock observations.To achieve this, we use UFalconv2 to produce a CMB lensing signal realization.To this, we add a CMB lensing noise realization.This is produced by subtracting the input cosmological signal as well as the mean field bias estimate from the Planck Full Focal Plane (FFP10) CMB lensing simulations available from the Planck Legacy Archive 11 .
Finally, we apply the Planck lensing survey mask.We downgrade this mask to our fiducial resolution of nside=1024 following the procedure outlined in [85] and further apodize the mask using a Gaussian smoothing kernel of FWHM 1.0 deg.This apodization scale is chosen as we find there is significant power leakage from noise-dominated small-scale modes until this scale [27,98].

Planck ISW cross-correlations
We make mock observations of the ISW signal on the map-level utilizing the new routine implemented in UFalconv2 (see § 2.2.4).To create a noisy ISW realization we add one Planck temperature FFP10 end-to-end simulation to each of the maps produced by UFalconv2.The simulations contain both a CMB signal and noise component.Adding these to our ISW signal maps is hence an approximation since we are double-counting the ISW signal component.We validate this approximation for our use-case by checking that a covariance matrix built using these noisy ISW realizations has deviations in the diagonal elements of no more than 10% compared to an analytically derived Gaussian ISW covariance.
We use the mask provided by the Planck team again downgrading to our fiducial resolution of nside=1024 and apodizing using a Gaussian smoothing kernel of FWHM 1.0 deg following Ref.[23].

Measuring spherical harmonic power spectra
Once we have produced the masked, noisy map-level mock observations for each of the LSS probes in the pipeline, we can measure the associated auto-and cross-correlation power spectra.The spherical harmonic power spectrum or C ℓ is defined as the ensemble average of the products of the spherical harmonic coefficients of two maps ψ lm and ϕ lm : We can estimate the spherical harmonic power spectra from our maps by computing the spherical harmonic decomposition and using the following estimator: However, we must account for the fact that the data is recorded over the partial sky when comparing the measured C ℓ s to theoretical predictions.We can generally relate the ensemble average of the C ℓ s measured on the cut-sky to the full-sky predictions via a mode coupling matrix: where < CXY ℓ > is the ensemble average of the cut-sky measurement, M ℓℓ ′ is the mode coupling matrix, and C ℓ is the full-sky prediction.Producing an unbiased estimate of the full-sky C ℓ (to compare with theory) then amounts to solving the inverse of equation 4.4.This is in general not possible given the implicit loss of information associated with measuring C ℓ s on the masked sky.To make progress, we use the MASTER algorithm [99] implemented in the code NaMaster [100].This enables an efficient solution to the inverse problem by assuming that the true angular power spectrum is piecewise constant over some bins.It is then possible to invert the associated binned mode coupling matrix and produce an unbiased "bandpower" estimate of the full-sky C ℓ for each bin.
Since the theoretical Cls are not, in reality, constant over the ℓ-range of a given bin, we must make a correction to these before comparing them with the bandpowers of the data.To achieve this we first couple the theoretical C ℓ s with the original mode-coupling matrix (M ℓℓ ′ ).Then we bin the resulting coupled C ℓ s into the same bandpowers as the corresponding data we wish to compare with and decouple the result with the binned mode-coupling matrix [100].This process accounts for the (small) bias induced by the assumption of step-wise full-sky Cls.Note for the cross-correlations we use the intersection of the two associated survey masks when measuring the cross C ℓ .
An important cross-check of the pipeline so far is made at this point in Fig. 4.Here we compare the mean C ℓ of 200 noiseless masked map-level mock observations against the theoretical predictions from the theory code and the emulator including the effects of mask deconvolution.Whilst Fig. 4 shows only a selection of the 43 total spectra in the pipeline, we checked that for each spectrum the deviation between the average of 200 simulated spectra and the emulated spectrum is less than 0.1σ across all bandpowers.The 1 − σ errorbar is determined from the square root of the diagonal of the multiprobe covariance matrix (described in § 5.1.3).

Scale cuts
Scale-cuts generally refer to the restriction of multipole range over which we perform cosmological analysis to the region where we trust the underlying theoretical model.For the CMB probes modeled at the likelihood level (see § 5.1.1),we take the scale-cuts and binning provided by the official analyses.We instead choose the scale-cuts wherever we independently implement the likelihoods (for example for the cross-correlations of the LSS probes).In this section, we describe the motivation behind the scale-cuts in this pipeline for each auto-and cross-correlation C ℓ (see Table 2 for an overview).

KiDS-1000
For the KiDS-1000 auto-correlations, we choose the same scale cuts (100 < ℓ < 1500) as used in [91].These scales represent the range over which the KiDS-1000 methodology is validated.Higher ℓ modes are discarded due to complications in modeling the non-linear matter power spectrum and baryonic processes on small scales.

BOSS DR12
We restrict the analysis of galaxy clustering (both auto-and cross-correlations) to large scales where the linear bias model implemented in the pipeline is accurate.Specifically, we choose ℓ max such that the Fourier modes entering the Limber integral satisfy k < 0.1Mpc −1 following [19].For the low-ℓ cut, we choose a conservative l min = 30 for the autocorrelations since, for lower ℓ values, RSDs become important [101] and the Limber approximation is no longer accurate [62].

KiDS-1000 cross-correlations with BOSS DR12
We generally follow the principle that the cross-correlation scale-cuts are given the two most conservative cuts applied to the respective auto-correlations (i.e., in this case, the high-ℓ cut imposed by the linear bias modeling of BOSS).Furthermore, we follow Ref. [14] in only including cross-correlation bins where there is not a significant overlap between the source galaxies of KiDS-1000 and the galaxies in the BOSS DR12 catalog.This is because in this regime the IA terms become important and it is not clear that the NLA model works sufficiently well to produce unbiased constraints.In particular, we choose to only include cross-correlations for C Planck CMB lensing cross-correlations For the Planck CMB lensing cross-correlations, we consider multipoles in the range 30 < ℓ < 400.The low-ℓ cut is motivated by the use of the Limber approximation used for the theoretical predictions which can be inaccurate for lower ℓ's [62].For high ℓ's, we use the conservative ℓ cut from the Planck 2018 lensing analysis of ℓ max = 400.This was chosen by the Planck team due to potential uncertainties in the lensing bias terms becoming increasingly important at smaller scales.
Planck ISW cross-correlations Whilst for the other probes we choose a conservative ℓ min = 30 due to the Limber approximation, we choose a more aggressive scale-cut of ℓ min = 20 for the ISW cross-correlations.This helps to boost the ISW signal-to-noise ratio (S/N) which peaks at large scales.We expect any bias due to the inaccuracy of the Limber approximation in this regime to be small, especially given the ISW cross-correlations are typically low S/N measurements compared to the other probes in the analysis.We do not include the T γ 1 , T γ 2 cross-correlations in the analysis as we found making the simulationbased estimates of these cross-correlations unstable due to the very low S/N.We leave the possibility of using beyond-Limber techniques [102,103], to push the analysis of the ISW to larger scales and thereby further increase the S/N, to future work.

Table 2:
The ℓ-cuts and binning for the auto-and cross-C ℓ s in this analysis.The CMB and C κκ ℓ are implemented at the likelihood level using Planck official data products.We, therefore, inherit the scale-cuts and binning from these analyses [86,104].Note the number of bins for these probes is not shown as these use binning schemes that are in general more complex than the simple logarithmic binning used for the LSS probes.5 Inference

Likelihood
In this section, we describe how we evaluate the likelihood in the combined probes framework.We begin by describing the CMB probes that are included in the analysis at the 2-point level (i.e.not through the map-making procedure described in § 4.1).We then detail how the mock data vector and covariance matrix are derived for each of the C ℓ s in the pipeline before combining to arrive at the likelihood.

CMB probes
In this pipeline, C T T ℓ , C T E ℓ , C EE ℓ and C κκ ℓ are modeled at the 2-point level.We follow this approach primarily because the CosmoGridV1 simulations go up to z = 3.5 and hence cannot produce maps for the higher redshift probes 12 .
CMB primary For the high-ℓ (ℓ > 30) part of the CMB primary TTTEEE likelihood, we follow the Python implementation of the Plik_lite likelihood, planck-lite-py 13 , described in Ref. [39].This likelihood is already marginalized over the Planck foreground parameters following the procedure first outlined in Ref. [105].The marginalization simplifies the likelihood significantly: this uses just one nuisance parameter, A planck , which accounts for the uncertainty in the overall Planck calibration.For the low-ℓ (ℓ < 30) TT and EE likelihood we use the log-normal binned likelihood, planck-low-py presented in [106]14 .Ref. [106] showed that D T T /EE ℓ = ℓ(ℓ+1)C T T /EE ℓ /2π, once binned into 2 and 3 bins for TT and EE respectively, closely follows an offset-lognormal distribution.This means that one can write a Gaussian likelihood in x = ln(D bin ) as: where x 0 is the offset parameter and µ and σ are the mean and width of the Gaussian distribution respectively.We take the values for these parameters found in Ref. [106] for the Planck PR3 low-ℓ likelihood.We validate our implementation of the Plik_lite likelihood and the low-ℓ lognormal bins by comparing to the official Planck chains from the full likelihood and find excellent agreement (see § A.1).

CMB lensing auto-correlations
Our CMB lensing analysis follows the usual quadraticestimator approach [51,107].Generally, when one measures the lensing potential using a quadratic estimator, the estimated C ℓ s are contaminated with biases which must be subtracted to extract the unbiased C ϕϕ L .Specifically, one must remove: • The N 0 bias.This is a 0th-order bias that accounts for the disconnected signal expected even in the absence of lensing.This can be generally assumed to be fixed with respect to varying cosmology and is typically estimated from simulations.
• The N 1 bias.This enters at first-order and is the non-Gaussian secondary contraction of the trispectrum.This is generally cosmology-dependent.
• Point source correction.This is due to the non-Gaussian signal produced by unresolved point sources in the data and can be assumed fixed and computed from the data.
In this analysis, since we use mock data, we generally do not have to consider the removal of terms that are fixed with respect to cosmological parameters (N 0 and point source correction).However, since the computation of the cosmology-dependent N1 bias relies on sampling the CMB primary C ℓ s at each point in parameter space, this can impact the shape and width of the contours.We hence include the effect of its removal in our likelihood to give realistic forecast constraints [108].To do this we implement the "CMB marginalized" likelihood supplied by the Planck team [86].This uses a minimum variance (MV) reconstruction of the lensing potential based on measurements of both the temperature and polarization data.In this likelihood, the dependence on the CMB primary spectra is marginalized over to obtain a likelihood that depends on cosmological parameters only through the lensing potential power spectrum C ϕϕ L .This procedure boils down to adding a term to the original covariance matrix and applying a linear correction to the theoretical prediction of this quantity.This approach was introduced in the regime of "lensing-only" analyses (i.e.not including the CMB primary , and validated only in this case [86].We, therefore, checked that using this likelihood in the case where we include CMB primary spectra does not introduce any biases compared to the official Planck TTTEEE+ lensing analysis, which is non-marginalized (see § A.2 for further details).

Mock data vector
In this analysis, we use two different mock data vector cosmologies.The baseline case is the fiducial cosmology of the CosmoGridV1, shown in Table 1.For τ reio and the nuisance parameters in the pipeline (which are not included in the simulation cosmology), we choose the baseline values shown in Table 3.The second mock cosmology is termed the "high-M ν " case which has M ν = 0.15 eV and therefore a raised neutrino energy density Ω ν .The increase in Ω ν is compensated by a decrease in Ω cdm such that the overall matter density Ω m = Ω cdm + Ω ν + Ω b is kept the same and all other parameters are identical to the baseline case.We use these two cosmologies to explore the neutrino mass constraints in two settings: (1) where the constraint looks like an upper limit (the baseline case) and ( 2) where there is a possibility of a Gaussian constraint (the high M ν scenario).We additionally vary the A L parameter for some of the mock observations away from the baseline value of A L = 1.

bias cmass 2.086
A planck 1.0 Table 3: Optical depth and nuisance parameter choices for the mock data vectors produced in this pipeline.Note that for the galaxy biases we choose the best-fit values corresponding to [19].
We follow the procedure outlined in § 4 to create mock observations for each of the LSS and associated cross-correlation C ℓ s in the pipeline.For the CMB primary and CMB lensing mock signals, we instead use a C ℓ computed by our theory code described in § 3.2.The C ℓ is then binned using the same procedure from the official analysis pipelines we have re-implemented.We compute our results using mock data vectors that do not additionally include a survey-specific noise draw, though survey noise is included in the covariance matrix (see § 5.1.3) 15.

Simulation-based covariance matrix
We estimate the covariance of the LSS probes in the pipeline from a set of 2000 correlated mock realizations.These are generated by applying both a survey mask and a survey-specific noise realization to the UFalconv2-generated cosmological signal maps (as described in 4).We add 10 random noise realizations for each of the 200 fiducial signal realizations to produce 2000 pseudo-random realizations.We use the NaMaster code to estimate the binned maskdeconvolved C ℓ s for each set of correlated mock observations from which we can estimate the LSS block of the covariance matrix.We add the covariance for C κκ ℓ to the LSS covariance matrix as follows.We We checked that the obtained covariance matrix is converged by analyzing the fractional change in the diagonal elements as a function of the number of realizations used to compute the covariance.We found that above ∼ 600 realizations, the mean change is at the subpercent level (see § C.1).We checked that the diagonal power of our covariance matrix built from simulations matches within 10% of the Gaussian analytic covariance matrix built using the NaMaster framework and further that using this Gaussian covariance as opposed to our simulation-based covariance produces only a small shift in the derived contours(see § C.2).

Combining the CMB and LSS likelihoods
In this pipeline, we treat the CMB primary contribution to likelihood as independent of the LSS likelihood.This is equivalent to summing the individual contributions at the level of the log-likelihood (log(L)): We, therefore, lose the cross-covariance between the CMB correlations (TTTEEE) and the other spectra in our pipeline (meaning we slightly underestimate the error bars in constraints).We expect this to have a negligible impact given that the amplitude of the crosscorrelation with LSS tracers is negligible compared to the CMB primary C T T ℓ , C T E ℓ and C EE ℓ .For the CMB lensing auto-correlation (C κκ ℓ ), where these cross-terms are important, we do include the cross-covariance with the other probes in the simulation-based covariance matrix (see § 5.1.3).
For the LSS probes (including the ISW and CMB lensing cross-correlations) and C κκ ℓ , we compute the likelihood using our full covariance matrix and the mock data vector assuming the Gaussian likelihood approximation: where ⃗ D is the mock data vector, ⃗ f (θ) is the theory prediction (from the emulator) and C −1 is the inverse covariance matrix.Note, we de-bias our estimate of the inverse covariance matrix following Ref.[109].To do this we multiply the naïve inverse by a factor: where n = 2000 is the number of quasi-independent realizations used to estimate the covariance and p is the number of data points (which is set-up dependent but takes a maximum value of p = 243 when all possible LSS correlations are included).

Data Compression
The need for effective data compression in cosmological analysis has been consistently emphasized, especially in the context of multi-probe studies [110,111].For example, compression is important in the case that the covariance is estimated from simulations where for a data vector of size n, typically ∼ n + 3 simulations are required to get an accurate estimate of covariance [111].Even in the case where the covariance is analytically estimated, reducing the dimensionality via data compression can help to improve the numerical stability of the inverse of the covariance matrix.Compression also allows for less memory usage and computational cost in the inference step of analysis due to the reduced size of the covariance matrix & data vector involved.We explore two distinct compression methods: (1) MOPED and ( 2) PCA.The compression algorithms are implemented at the stage of the pipeline just before inference takes place.The covariance matrix and mock data vector are compressed using the methods described below and passed into the inference pipeline.The associated compression vectors are also supplied so that the theoretical predictions can be compressed at each stage of the inference.
The MOPED [37] algorithm is a compression technique that reduces the dimensionality of a dataset (N data ) to the number of parameters being constrained.This works by finding a decomposition of the data vector onto a set of 'MOPED basis vectors' (b i ) which are chosen to maximize the sensitivity of the likelihood to each of the parameters being constrained.The set of compression vectors is found by weighting the derivative of the given summary statistic (here, the C ℓ s) with respect to the input cosmological parameters by the inverse of the covariance matrix and subsequently using Gram-Schmidt orthonormalization: where Y i is the compressed data-vector and C and Ψ are respectively the covariance matrix and theoretical C ℓ at the fiducial cosmology.Ref. [37] showed that MOPED compression is optimal and lossless for Gaussian-distributed data if the input covariance matrix is parameterindependent (see also [112]).The main computational challenge involved in applying this technique is to find accurate derivatives of the data vector at the fiducial cosmology.We compute the derivatives using the AD afforded to us by using an emulator-based approach for the theoretical predictions.As in [78], we found that this approach was consistently more robust than the more rudimentary 'finite differences' method (see also [113] for a recent demonstration of using AD for MOPED compression).
The second technique, PCA, allows us to find a new set of coordinates that can efficiently represent the data by choosing the modes giving rise to the highest variance across a dataset.In the present case, we perform a PCA compression on the emulator training set to find the vectors that capture the most variance as we span the parameter space within the priors (following an approach similar to [77,114]).PCA is a useful comparison to MOPED as one is free to keep as many PCA components (N P CA ) as desired.This means that in the limit of N P CA → N data we will always recover the full information content of a dataset and hence this scheme is guaranteed to be lossless for some choice of N P CA .

Fiducial model
The fiducial model in this study is standard ΛCDM.When exploring neutrino mass constraints we additionally allow the M ν parameter to vary.In this case, we model neutrino masses using three degenerate mass states given this is sufficiently accurate for the precision of current data [115][116][117].In some cases, we also include a free lensing amplitude rescaling A L parameter.This phenomenological parameter, first introduced in [35], rescales the amplitude of the lensing potential in the calculation of the lensed CMB primary anisotropies (i.e. the TTEEEE spectra): (5.6) Note the lensing potential is not rescaled when computing the CMB lensing reconstruction.This is used as an internal consistency check of CMB data where the expected value is A L = 1.Deviations from A L = 1 could indicate internal systematics in CMB data that can have consequences for parameter constraints as will be discussed in more detail in § 6.3.2.
In addition to the extension parameters M ν and A L , we also vary several nuisance parameters specific to each of the probes: Weak lensing For the WL modeling, we allow for a free ∆ z shift in the mean redshift for each of the bins as well as a free IA amplitude (A IA ).Since in this analysis, we rely on DMonly simulations, we make a simplification compared to the fiducial KiDS-1000 analysis by not modeling the effect of baryonic feedback, which can affect the matter power spectrum at scales included in our forecast analysis.Finally, the KiDS-1000 block of the covariance matrix is increased by an additive term C m to account for the uncertainty in the shear calibration as described in Ref. [118].
Galaxy clustering For galaxy clustering, we allow for one linear bias parameter for each bin.Using the linear bias approach limits our analysis (for both auto-and cross-C ℓ s) to large, linear scales where the model is accurate [19].We do not include any nuisance parameters that account for foreground contamination 16 .
CMB For the CMB analysis, we include the A planck nuisance parameter wherever we encounter a CMB temperature/polarization auto-or cross-correlation C ℓ (including ISW crosscorrelations).This parameter is defined via an overall rescaling at the map level, and hence at the C ℓ -level this becomes: where X represents an LSS probe in the ISW cross-correlation.We include this in the likelihood.We do not vary any nuisance parameters for the CMB lensing part of the likelihood.

Inference pipeline
To perform cosmological inference in this pipeline we use the emcee package [119].The emcee algorithm uses an ensemble of walkers in parallel and capitalizes on the affine-invariant ensemble sampling technique to efficiently explore the parameter space, providing accurate Table 4: The parameters varied in the analysis and the associated priors.
and computationally cost-effective parameter space exploration.Using emceee combined with our emulators enables running a full MCMC analysis in O(∼ minutes).
The priors for the cosmological/nuisance parameters in the analysis runs are shown in Table 4.We use the same priors for all results shown in the main text of this paper.These are typically uninformative, flat uniform priors, with only two exceptions.First, we follow [57] in placing a correlated Gaussian prior on the ∆ z shift nuisance parameters for which we use the covariance matrix provided in the KiDS-1000 data release.This encodes the correlated errors from the photometric redshift distribution determination procedure [57].Second, following [104] we place a tight Gaussian prior on the re-scaling parameter A planck .

Baseline ΛCDM constraints
The baseline ΛCDM posterior contours are shown in Fig. 6.These contours, shown in blue, represent the constraints from analyzing the full set of mock data at the baseline cosmology.This consists of a combination of the autocorrelation signals: CMB primary (C T T ℓ , C T E ℓ , C EE ℓ ), CMB lensing (C κκ ℓ ), WL (C γγ ℓ ) and galaxy clustering (C δδ ℓ ) as well as the cross-correlations: . We find that the pipeline produces unbiased contours in the ΛCDM setting: the fiducial model parameters (the purple stars) are within the marginalized contours.We also find this when combining smaller subsets of auto-and cross-correlations in the pipeline (e.g. a mock WL-only analysis).
Whilst this full combination of data is somewhat optimistic given the current levels of tensions between datasets (see e.g.[120]), the constraining power of the 12 × 2pt combination is nevertheless interesting to compute.We find that this combination of data enables a determination of the Hubble constant at the 0.5% level and Ω m at the 1.5% level improving the constraints on these parameters from the official Planck CMB TTTEEE + CMB lensing analysis by ∼ 30% in each case [34].The increase in constraining power comes from the consistent addition of the LSS auto-and cross-correlation C ℓ s.

Data compression for combined probes
In this subsection, we compare the PCA and MOPED compression algorithms.In Fig. 7 we show the comparison of the two techniques for a mock KiDS-1000 analysis (left) and a mock Planck analysis (right).The uncompressed contours are shown in blue and represent the contours derived from the full data vector including all of the extractable information from the 2-point statistics.In red and green are the contours for the same setting but using PCA  compression with different numbers of components.We find that ∼ 40 and ∼ 100 components are required respectively to faithfully reproduce the uncompressed contours (green contours in Fig. 7) whilst fewer components lead to deviations (red contours in Fig. 7).In contrast, the MOPED compression algorithm condenses the information to a data vector of size just 11 in both cases whilst faithfully reproducing the original contours.The compression to 11 data points with MOPED can be explained as follows.For the KiDS-1000 WL case, we have 5 cosmological (ΛCDM parameters but with fixed τ reio as this is not constrained by LSS data) and 6 nuisance parameters (A IA and 5 ∆z shifts).For the Planck CMB primary case we compress the high-ℓ likelihood to just 6 data points (one for each of the 6 cosmological parameters of ΛCDM).We then add the 5 log-normal bins for the low-ℓ T/E likelihood giving overall 11 data points.
In Fig. 6 we apply the MOPED compression scheme to the full combined probes set-up.The uncompressed baseline constraints use 859 data points for all of the CMB and LSS data auto and cross-correlations.In principle, both of the CMB primary (C T T ℓ , C T E ℓ , C EE ℓ ) part and the LSS part of the likelihood could be compressed together to produce a data vector equal the total count of variable parameters (plus 5 lognormal bins for the compressed low-ℓ T/E likelihood).We, however, choose to apply a separate compression of these components.This is necessary as we found that the combined covariance matrix of LSS and CMB is too large to be inverted without numerical instability which induces inaccuracies in the computation of the MOPED compression vectors and ultimately results in biased contours.
Using the separate MOPED compression scheme described above, we demonstrate that the data vector can be drastically reduced from 859 to just 24 data points whilst preserving the information contained in the full combination of CMB + LSS.This consists of 13 data points for the LSS block where we vary 5 cosmological parameters and 8 nuisance, and 11 for the CMB block where we vary 6 cosmological parameters and add 5 log-normal bins for the compressed low-ℓ likelihood.We find that whilst there is a slight loss in constraining power in some of the 2D parameter projections (most notably in the Ω m -σ 8 plane), the 1D marginalized parameter constraints in the MOPED compressed scenario show excellent agreement with the uncompressed case with the mean and standard deviations of each parameter constraint agreeing to within 1%.
In addition to the results on mock data, we checked that the MOPED algorithm is also able to reproduce the full Planck high-ℓ TTTEEE + low-ℓ T/E + lensing results with just 16 data points.Here we again perform a separate compression and have 11 data points for the CMB primary likelihood and 5 for the CMB lensing contribution.using the real survey data (see A.2 for further details).We release, for the first time, this MOPED compressed Planck PR3 CMB primary + CMB lensing likelihood publicly17 .

Constraining power
In this section, we present the constraining power on M ν for various combinations of datasets using a mock data vector in the baseline and high M ν scenarios in Fig. 8 & Fig. 9.A change in the neutrino energy density generally affects both the background evolution of the scale factor and the evolution of the metric potentials at linear-order in perturbation theory.We can therefore measure the neutrino mass by finding the signatures of these effects on the different cosmological probes.One problem commonly encountered with single-probe analyses of neutrino masses is degeneracies in the way that M ν and the other parameters in the ΛCDM model affect the measured spectra.For example, in the CMB alone case, there are clear degeneracies between M ν and σ 8 , h and Ω m (dark blue contours in Fig. 8 & Fig. 9).Whilst the full combination of LSS data (WL, galaxy clustering, CMB lensing, and associated cross-correlations including ISW) is not able to directly constrain M ν (see green contours in Fig. 8 & Fig. 9), the addition of this data enables breaking degeneracies in the constraints from the CMB alone (see also [19,115,121,122]).In particular, the addition of LSS data provides a CMB-independent measure of both σ 8 and Ω m and hence enhances the M ν constraint by breaking the degeneracy with these parameters.In the baseline cosmology with M ν = 0.06eV, we are able to place an upper limit M ν < 0.16eV at the 95% CL when combining CMB+LSS compared to M ν < 0.25eV using CMB TTTEEE data alone representing an improvement of ∼ 35%.The combined limit is slightly less constraining than the Planck +BAO result of M ν < 0.12eV reported in [34].We hypothesize this is because in our case the LSS data helps to exclude the very low M ν < 0.06eV region (see the left-hand side of the dark orange 1D-marginalized M ν constraint in Fig. 8).This biases the 95% CL constraint to larger values meaning it is not a fair measure of the constraining ability of this data combination with respect to M ν .Another reason may be due to the internal A L tension in Planck PR3 data biasing this constraint towards small M ν [36].In the high M ν scenario we are able to detect a neutrino mass at the ∼ 3σ-level (M ν = 0.15 ± 0.06eV).This is similar to the constraining power found in Ref. [123] using a combination of CMB data from Planck and SPT-3G, LSS and supernova data.

Impact of A L
Recently, several studies using data independent from Planck have been able to place weak limits on the neutrino mass at regions of the parameter space (M ν ∼ 0.2eV) that are at face-value strongly ruled out by the Planck -derived upper limits (see for example [36,123]).One postulated reason for this discrepancy is due to the internal lensing tension in the Planck 2018 dataset which may bias the neutrino mass constraints found using this data towards M ν = 0eV.As described in § 5.3, the internal lensing consistency of a CMB dataset can be parameterized via A L .The expected value of this parameter is A L = 1 in the standard ΛCDM model.However, the Planck collaboration PR3 analysis finds A L = 1.18 ± 0.065 from primary anisotropy CMB data (detection of A L > 1 at 99% CL).The effect of lensing on the CMB is anti-correlated with M ν and therefore such an anomaly can bias neutrino mass constraints towards lower values (see discussion in e.g.[36,124]).Other CMB experiments have found values of A L more consistent with the expected value (ACT finds A L = 1.01 ± 0.11 in their DR4 analysis [125]) so it remains unclear whether this is due to unaccounted physics (see for example [126]) or due to a systematic error in the Planck data.It should also be noted that a recent re-analysis with the NPIPE pipeline presented in Ref. [127] finds A L = 1.095 ± 0.056 representing a milder 1.7σ tension.In any case, it is interesting to consider the impact of A L on neutrino mass constraints in the context of this novel combined probes pipeline.
In Figures 10 & 11, we explore the impact of artificially including an A L systematic on the neutrino mass constraint.To achieve this, we create mock CMB primary TTTEEE data vectors with varying A L .However, we do not vary this parameter in the subsequent analysis, thereby treating this modification as an unaccounted-for systematic in the data.Given the anti-correlation between A L and M ν , this leads to a bias in our M ν constraints    which we quantify in this section.The left-hand plots in Figures 10 & 11 show that when considering CMB data alone there is a strong induced bias towards low M ν .For example, when using the baseline cosmology, the CMB alone constraint with A L = 1 is M ν < 0.25 (95 % CL).As the underlying A L in the signal is increased, the upper limit constraint on M ν monotonically decreases.In particular, for the central A L value found with Planck PR3 data (A L = 1.180 [86]), we would find an upper limit of M ν < 0.13eV.Considering the right-hand plot of Fig. 10, we find that the inclusion of LSS data in the full-combined case produces neutrino mass constraints that are more robust to systematic issues of this type.In particular, for A L = 1.18 the upper limit is M ν < 0.15eV compared to the no-systematic (A L = 1.0) baseline of M ν < 0.16eV.The addition of external data also lessens the impact of the A L internal systematic on the σ 8 and h constraints.The corresponding plot in the high M ν scenario (Fig. 11) also demonstrates that the neutrino mass constraint is also stabilized by the addition of LSS data.In particular, the constraint found for A L = 1.18, is M ν = 0.14 ± 0.07 compared to M ν = 0.15 ± 0.06 for the no-systematic baseline.These results suggest that robust neutrino mass constraints can be derived by combining CMB and LSS data even in the presence of an A L -like systematic in the CMB data.

Conclusions
In this work, we introduce a new pipeline to rapidly infer cosmological parameters from a multi-probe analysis including a combination of LSS and CMB data.The pipeline is built upon correlated, map-level realizations of these probes, generated by applying the UFalconv2 package to the CosmoGridV1 suite of simulations.We use these realizations to estimate a non-Gaussian covariance matrix and a mock data vector for the LSS probes.To overcome the challenge of large data vectors for combined probe analyses, we integrate both the MOPED and PCA data compression schemes into the pipeline.We additionally make use of a neural network emulator for accelerated theoretical predictions.We validate our approach by comparing simulations and theoretical predictions, and we also compare our derived constraints to earlier analyses (see Appendix A).We apply our framework to a simulated 12 × 2pt tomographic analysis of KiDS, BOSS, and Planck, and forecast constraints for both a standard ΛCDM model and νΛCDM.In the standard scenario, we are able to determine h, and Ω m to within 0.5% and 1.5% respectively from the full 12 × 2pt combination, representing a ∼ 30% improvement compared to the official Planck CMB TTTEEE + CMB lensing analysis [34].
As part of this analysis, we demonstrate the power of differentiable emulators both in terms of rapid parameter space exploration and in terms of accurately predicting derivatives, dC ℓ /d ⃗ θ, which we find to be crucial for robust MOPED-based data compression.Compared to PCA compression, MOPED provides a more extreme compression in all cases, and using this we demonstrate a lossless compression of the full combined-probes data vector from 856 to just 24 data points.We also release a Planck CMB TTTEEE+lensing likelihood compressed to just 16 data points and show that this produces constraints in excellent agreement with the official results.
We additionally explore neutrino mass constraints in the multiprobe setting.We find that, whilst the majority of the constraining power comes from primary CMB data, the combination of CMB lensing, WL, and galaxy clustering data helps to exclude regions of parameter space by breaking degeneracies.In this optimistic setting with all simulated data sharing the same underlying cosmology, we forecast a constraint of M ν < 0.16eV for the fiducial case where M ν = 0.06eV and determine a ∼ 3σ detection in the high M ν scenario (M ν = 0.15 ± 0.06).We reconfirm that an A L -like internal tension in CMB data can bias determinations of M ν [36,128,129], for example when analyzing CMB data alone we find constraints of M ν < 0.25eV for A L = 1.0 and M ν < 0.13eV for A L = 1.18 (the central value in Planck PR3 data [86]).However, the impact of A L -like systematics on determinations of M ν is significantly mitigated by the addition of LSS data, which helps to anchor the constraints towards the correct cosmological parameters in both the fiducial and high M ν scenarios.This analysis relies on several simplifying assumptions, such as neglecting baryonic feedback effects on WL and restricting to linear bias models for galaxy clustering.In the future, it will be interesting to relax these assumptions and extend our analysis to smaller scales, which we expect to boost the constraining power of the LSS data.Another natural extension of the pipeline is to add additional datasets: On the CMB side, we envisage including ACT CMB primary or CMB lensing data [124], as well as the recently released likelihoods from Planck PR4 maps [127].On the LSS side, it would be interesting to include DES data [12,130,131] and the data from the Dark Energy Spectroscopic Instrument (DESI) [132,133].
This analysis presents a further demonstration of the power of combined probes analyses for deriving strong constraints, robust to systematics.The results and extensive validation tests presented as part of this work strongly motivate applying the combined probes pipeline to current and future survey data, thus deriving robust constraints on ΛCDM and its extensions.

A.1 Planck 2018 primary
In order to validate our approach for the individual surveys, we compare our results with official analyses where possible.For the Planck primary C T T ℓ , C T E ℓ and C EE ℓ case, we compare our emulator implementation with the official chains downloaded from the Planck legacy archive.note these chains are for the full likelihood that is not marginalized over the foreground parameters whereas in the pipeline we implement the Plik_lite likelihood.In addition to this, we also validate and release, a MOPED compressed high-ℓ TTTEEE likelihood consisting of just 6 data points.Combining with the 5 bins of the planck-low-py log-normal Gaussian binned likelihood from [106] the full Planck likelihood can be compressed to just 11 data points that faithfully reproduce the full-likelihood within ΛCDM.Our constraints (both in the uncompressed and compressed case) match with the official analysis to within 1% in both the σ and mean for each 1D marginalized parameter constraint.

A.2 Planck 2018 lensing
For Planck CMB lensing, we compare our emulator implementation with the official chains in two scenarios: firstly on the left-hand plot of Fig. 13 we derive "lensing only" contours.For these, we place the same priors on the varied cosmological parameters as in Ref. [86] (gaussian priors of ω b = 0.0222 ± 0.0005, n s = 0.96 ± 0.02, a uniform prior of 0.4 < h < 1 and fixing τ reio = 0.055) 18 and compare to the official chains.In addition to the cosmological parameters we also derive the posterior distribution of the "lensing parameter" defined by [86]: this parameter describes the degeneracy direction probed by CMB-lensing data [108].We find excellent agreement in each cosmological parameter and, especially in the lensing parameter validating that the likelihood implementation is accurate in the "lensing only" regime.On the right-hand plot of Fig. 13, we compare our implementation of the Planck CMB lensing likelihood in the case that the primary CMB spectra (C T T ℓ , C T E ℓ and C EE ℓ ) are also varied.As stated in § 5.1.1,we are using the "CMB marginalized" likelihood for C κκ ℓ which has only been validated in the case of analyzing CMB lensing alone.We, therefore, check that we are in good agreement (sub percent level match in the mean and σ of the marginalized 1D parameter distributions) with the official CMB + lensing results from Planck which are unmarginalized and so account for the effect of the CMB primary spectra on the lensing normalization and N 1 bias at each step of the inference [34].In addition to this, we also validate and release, a MOPED compressed "CMB marginalized" lensing likelihood consisting of just 5 data points (one for each parameter of the ΛCDM model excluding τ reio which is not varied).Combining with the planck-low-py log-normal Gaussian binned likelihood from [106] and the MOPED compressed Planck high-ℓ TTTEEE likelihood, the full Planck CMB + lensing likelihood can be compressed to just 16 data points that reproduce the full-  19 .Finally, our emulators are trained over a set of fixed parameters that are not identical to those used in the official KiDS-1000 bandpower analysis.We match the priors where possible to make this comparison (see Table 5).The comparison between our implementation and the band power results is shown in Fig. 14.The official results give a slightly lower S8 compared to our implementation which can be 19 The KiDS team tested the impact of this on real data analysis and found sub 1σ shifts in the value of S8.
due to the different covariance matrix, the lack of a baryonic feedback implementation and the sampled parameters/priors not being exactly matched.We note that our results match very closely to the "power spectra" results of [90] which used the same set of simulations in computing the covariance matrix.To evaluate the suitability of a given emulator for use in cosmological inference it is important to consider the impact on the eventual parameter constraints.We follow a procedure similar to Ref. [79] in assessing the emulator accuracy in terms of the relative difference between the true C ℓ derived from a theory code and the emulated spectrum compared with the error associated with the data used in the inference step.This pragmatic approach to emulator accuracy informs our decision to train just six emulators to cover the C ℓ s in the pipeline.The data associated with the CMB probes (C T T ℓ , C EE ℓ ,C T E ℓ ,C κκ ℓ ) has a small error budget and hence we require a high emulator fidelity.We, therefore, train a dedicated emulator for each of these C ℓ s.We also train the models to learn the mapping from cosmological parameters to the full unbinned spectra (from 2 < ℓ < 2508) so that the model does not additionally have to learn a survey-specific binning scheme (we hence apply the required binning at the inference step).The data associated with the LSS probes, on the other hand, has a larger error margin and hence we find that we do not need dedicated models: just two networks are able to predict all of the LSS auto-and cross-correlation C ℓ s to sufficient accuracy.In this case, each emulator outputs a concatenation of multiple binned C ℓ s which already include the effect of convolution and deconvolution with the respective survey mixing matrices (computed using NaMaster).The networks are trained over different parameter sets given we only include the nuisance parameters associated with the particular spectra that the network is trained on.The parameters and associated prior ranges for the emulators are shown in Tables 6 & 7.

B.2 Emulator accuracy
We follow the approach taken in [79] in testing the emulator accuracy by considering the fractional error in the emulation of each of a set of 60000 test spectra divided by the error given approximately by: σ(C ℓ ) = COV (ℓ, ℓ). (B.1) We follow this procedure for each of the emulators in the pipeline and plot the distribution of the errors.We checked that 99% of the testing samples have a relative error < 0.15σ and hence that the emulator accuracy is unlikely to bias our derived constraints.In Fig. 15 we show the distributions for all of the spectra in our pipeline.

C.1 Covariance matrix convergence
To validate that our covariance matrix is converged, we measure the fractional change in the diagonal elements of our covariance matrix as a function of the number of realizations used to estimate it (see Fig. 16).We find that the 2000 realizations we use (200 underlying simulations each with 10 random noise realizations) are sufficient to produce a sub-percent mean fractional change with a negligible spread across the different diagonal components.

C.2 Comparison with Gaussian analytic covariance
A common approach to estimating covariances in cosmology is to assume a negligible non-Gaussian component to total covariance.The covariance matrix can then be computed analytically as a Gaussian piece plus some corrections (e.g. the super-sample covariance).It is possible to write an analytic expression for the Gaussian contribution to the covariance matrix given a) a mock data vector at some fiducial cosmology, b) an estimated noise vector c) the survey mask [100].To compare our simulation-based approach to the analytic method, we first generate an analytic covariance using the NaMaster framework which handles solving for the Gaussian covariance given the inputs described above.We checked that there are deviations of no more than 10% over the diagonal elements of these matrices.Then, in Fig. 17, we compare the results of our full LSS analysis (combining WL, galaxy clustering, CMB lensing, and all associated cross-correlations including ISW).We observe excellent agreement between the constraints, which is a non-trivial result given the completely independent methods used to estimate the two covariance matrices.This agreement provides further evidence that our simulation-based estimation of the covariance matrix works as expected.

Figure 2 :
Figure 2: Example maps produced by post-processing the CosmoGridV1 with the UFalconv2 software.For the clustering map we use the redshift distribution from the BOSS CMASS selection [33] and a corresponding linear bias value of b = 2.086, and for the shear map we use the first redshift bin from KiDS-1000 [57].All maps plotted above are dimensionless.

Figure 4 :
Figure 4: A comparison of the mean of 200 UFalconv2-derived mock (noiseless) signal observations at the fiducial cosmology C ℓ s with theoretical predictions from the base theory code and the multiprobe emulator for a selection of spectra.The C ℓ s measured from the mock observations are mask-deconvolved and the theory predictions are treated to include the effects of mask deconvolution following the procedure outlined in 4.2.The 1 − σ errorbars are determined from the square root of the diagonal of the multiprobe covariance matrix described in § 5.1.3.
use the CMB lensing mocks to estimate the C κκ ℓ − C XY ℓ -type (where X and Y stand for any of the other LSS probes) terms of the covariance matrix.We then insert the 'CMB marginalized' official covariance matrix into the C κκ ℓ − C κκ ℓ sub-block of the covariance matrix.The full LSS combined probes correlation matrix (including the C κκ ℓ contributions) is shown in Fig 5.

Figure 5 :
Figure 5: The correlation matrix for the LSS probes including the contribution from CMB-κ and the associated cross-correlations.This is derived from 2000 quasi-random mock realizations of fiducial cosmology.Each mock realization consists of a set of noisy mask-deconvolved C ℓ s which are derived from the same underlying dark matter simulation.This means the signals have the correct cross-correlation power.

Figure 6 :
Figure 6: The baseline 12 × 2pt ΛCDM results from the pipeline (blue) and comparison to MOPED compressed likelihood (red).The fiducial cosmology of the input mock data vector is shown in the purple stars and lines.The MOPED compressed data vector contains just 24 data points compared to the original 859 and preserves the information contained in the full data vector.

Figure 7 :
Figure 7: A comparison of the MOPED and PCA compression schemes: left KiDS-1000 WL, right Planck CMB primary TTTEEE example.The fiducial cosmology of the input mock data vector is shown in the purple stars and lines.MOPED outperforms PCA in both cases in terms of reducing the size of the data vector whilst preserving the information in the likelihood and avoiding biases.

Figure 8 :
Figure 8: The constraining power of different combinations of data at the baseline cosmology (M ν = 0.06eV).The fiducial cosmology of the input mock data vector is shown in the purple stars and lines.The lss line in green corresponds to the constraint from WL, galaxy clustering, CMB lensing, and all associated cross-correlations (including ISW).The parameter constraint values shown above the contour plots are obtained from the full 12×2pt combination (cmb+lss, in dark orange).

Figure 10 :
Figure 10: Exploration of the impact of an A L -like systematic for a data vector at the fiducial cosmology (M ν = 0.06eV): left CMB alone right full 12 × 2pt combination.

Figure 11 :
Figure 11: Same as Fig. 10 but in the high M ν scenario.

Figure 12 :
Figure 12: Comparison of the Planck official high-ℓ TTTEEE + low-ℓ T/E chains [34] with the chains derived from the implementation in the combined probes pipeline and the MOPED compressed version.

Figure 14 :
Figure 14: Comparison of the official KiDS-1000 bandpower results [57] with the pseudo-C ℓ s implementation in the combined probes pipeline.

Figure 15 :
Figure 15: The distribution over 60000 testing points of the relative difference between the emulator predictions and the theory code predictions, scaled by the 1σ error bars of the modeled surveys.The different colors represent the 68/95/99% limits of the distributions of each spectrum.

Figure 16 :
Figure 16: The mean fractional change in the diagonals of the covariance matrix (M ean(∆C N ii /C N −1 ii )) as a function of the number of realizations used to generate the matrix.The error bars in purple represent the standard deviation of this fractional change across the different diagonal covariance matrix elements.

Figure 17 :
Figure 17: Comparison of constraints using an analytic covariance vs simulation-based covariance for a combination of WL, galaxy clustering, CMB lensing, and associated crosscorrelations including ISW.

Table 1 :
for further details).Fiducial Cosmological parameters used in the CosmoGridV1.

Table 6 :
B Emulator set-up and accuracyB.1 Distribution of spectra over networks and prior ranges Priors used for the three dedicated CMB primary networks which output C T T ℓ ,

Table 7 :
Priors used for the LSS emulators in this pipeline.Identical priors are used for the CMB lensing auto correlation C κκ ℓ except that the WL/clustering nuisance parameters are not varied.