The growth of density perturbations in the last $\sim$10 billion years from tomographic large-scale structure data

In order to investigate the origin of the ongoing tension between the amplitude of matter fluctuations measured by weak lensing experiments at low redshifts and the value inferred from the cosmic microwave background anisotropies, we reconstruct the evolution of this amplitude from $z\sim2$ using existing large-scale structure data. To do so, we decouple the linear growth of density inhomogeneities from the background expansion, and constrain its redshift dependence making use of a combination of 6 different data sets, including cosmic shear, galaxy clustering and CMB lensing. We analyze these data under a consistent harmonic-space angular power spectrum-based pipeline. We show that current data constrain the amplitude of fluctuations mostly in the range $0.2<z<0.7$, where it is lower than predicted by Planck. This difference is mostly driven by current cosmic shear data, although the growth histories reconstructed from different data combinations are consistent with each other, and we find no evidence of systematic deviations in any particular experiment. In spite of the tension with Planck, the data are well-described by the $\Lambda$CDM model, albeit with a lower value of $S_8\equiv\sigma_8(\Omega_m/0.3)^{0.5}$. As part of our analysis, we find constraints on this parameter of $S_8=0.7781\pm0.0094$ (68\% confidence level), reaching almost percent-level errors comparable with CMB measurements, and 3.4$\sigma$ away from the value found by Planck.


Introduction
The current era of precision cosmology has seen a tremendous growth in the volume, quality, and variety of astronomical data that have become available. This has allowed us to progressively improve the constraints on the cosmological parameters of the Λ Cold Dark Matter (ΛCDM) model. The ΛCDM model is a remarkable fit to most observations, ranging from the cosmic microwave background (CMB) [1] to large scale structure (LSS) inferred from galaxy clustering and weak lensing [2][3][4], yet a number of intriguing discrepancies or "tensions" have begun to arise. A particularly important tension is related to the amplitude of density fluctuations (or linear matter perturbations) at low redshifts predicted by CMB data in comparison with direct measurements by cosmic shear and galaxy clustering data. This is commonly summarized in the so-called S 8 parameter, defined as where σ 8 is the variance of the linear matter overdensity field in spheres with a 8 h −1 Mpc radius, and Ω m is the fractional energy density of non-relativistic matter, both defined at z = 0. The exponent α is chosen to minimize the correlation between S 8 and Ω m , although the choice α = 0.5, which we use here, is a good approximation in most cases. This quantity has been measured by various weak lensing surveys [4][5][6] , encountering different levels of tension with the value inferred from measurements by the Planck satellite. Overall, cosmic shear data tend to recover somewhat smaller values of S 8 than the CMB and, to date, the strongest disagreement has been reported by the Kilo-Degree Survey collaboration (KiDS) [4], with a significance around 3σ. State-of-the-art measurements now exploit the complementarity between weak gravitational lensing and galaxy clustering [2,4,7]. Weak lensing is a direct tracer of the integrated density perturbations between source and observer, and therefore can be used to measure their amplitude in a relatively clean way. Its cumulative nature, however, washes out most features in the distribution of these perturbations, and makes it difficult to track the evolution of the amplitude in detail [8,9]. Galaxy clustering, on the other hand, is a high signal-tonoise, but biased, tracer of the local (as opposed to integrated) density perturbations. The galaxy bias makes it difficult to extract information about the growth of structure from the projected galaxy distribution alone, although important features in this distribution (e.g. baryon acoustic oscillations) can still be recovered [10,11]. In combination with weak lensing data, however, galaxy clustering improves our ability to reconstruct the history of density perturbations, significantly enhancing the associated cosmological constraints [2,4,[7][8][9]12].
This paper focuses on reconstructing the growth history, i.e. the evolution of the amplitude of density perturbations (characterized by σ 8 or S 8 ) as a function of redshift. The motivation for this endeavour is twofold: on the one hand, reconstructing this evolution from existing data sets will allow us to understand the redshift ranges over which current data are able to constrain the growth of structure, potentially shedding light on the origin of the aforementioned "S 8 tension". On the other hand, the amplitude of density perturbations is, arguably, the natural "observable" to which projected weak lensing surveys are sensitive, in the same way that the density-weighted growth rate, f σ 8 , is for redshift-space distortions and the transverse and longitudinal distance indicators, α ⊥ and α , are for baryonic acoustic oscillations (BAOs) in spectroscopic clustering data sets 1 [3,15,16]. Thus, this analysis will allow us to present current constraints from photometric redshift surveys in terms of observables, as opposed to final parameter constraints, which can then be used to e.g. study deviations with respect to ΛCDM.
To do this, we will make use of various galaxy clustering and weak lensing data sets. In particular, we will use galaxy clustering and cosmic shear data from the first data release of the Dark Energy Survey (DES, [2]), cosmic shear data from the fourth data release of the KiDS collaboration [17], galaxy clustering from the DESI Legacy Survey [18,19], the clustering of high-redshift quasars in the extended Baryon Oscillation Spectroscopic Survey [20], and maps of the CMB weak lensing convergence made by the Planck collaboration [21]. These data will allow us to recover the growth history in the range 0.2 z 2, as well as to compare the histories reconstructed by the two different cosmic shear experiments (DES and KiDS) independently, and by the combination of clustering and CMB lensing in the absence of shear data.
A key aspect of this work is that our analysis will be based on an independent estimation of the relevant two-point correlations between data sets using a consistent harmonic-spacebased framework to estimate both angular power spectra and their covariance matrix. This has the added value of being able to compare the constraints obtained by the different experiments, as well as their combination, under the same analysis pipeline (including the modelling of systematic effects). As we shall show, in some cases this leads to a significant improvement (in terms of goodness of fit, for example) with respect to previous analyses.
For the impatient reader, our main findings are: • Decoupling the linear growth of fluctuations from the background expansion, and modelling the former through quadratic splines, we show that existing large-scale structure data prefer a lower amplitude of fluctuations in the range 0.2 z 0.7 than that predicted by Planck by ∼ 5%. This is also the range of redshifts where current data are most sensitive to structure growth.
• This result is recovered consistently by independent data set combinations, and is driven by existing cosmic shear data.
• The recovered growth history is in good agreement with a ΛCDM model, although with a lower value of S 8 than that predicted from CMB. From a combined analysis of all the data, we obtain a constraint S 8 = 0.7781 ± 0.0094. This is in ∼ 3.4σ tension with current CMB data and with smaller uncertainties.
This work is organized as follows. Section 2 describes the theoretical background used to model the different auto-and cross-correlations, as well as the method used to parametrize the growth history. Section 3 describes the different data sets used in our analysis, together with the methods used to process them into maps tracing the projected mater overdensities. In Section 4 we present the methods used to estimate all power spectra and their covariance matrices, and the likelihood used to connect data and theory. Section 5 presents the main results of this analysis in terms of parameter constraints on ΛCDM and on the growth history. The main results are then discussed and summarized in Section 6.

Projected anisotropies
Our analysis will be based on a set of cross-correlations between fields defined on the celestial sphere u(n) that are related to a three-dimensional quantity U (x, z) in the lightcone through line-of-sight integrals with a radial kernel q u (χ) [8]: u(n) = dχ q u (χ) U (χn, z(χ)), (2.1) where χ is the comoving radial distance. In the Limber approximation [22], appropriate for the broad kernels used here, the angular power spectrum between two projected fields u and v is related to the three-dimensional power spectrum of their associated quantities U and V via: Our analysis will consider three fields: the projected galaxy overdensity δ g , the galaxy shear γ G , and the CMB lensing convergence κ: • δ g is related to the three-dimensional galaxy overdensity ∆ g via a radial kernel proportional to the redshift distribution of sources in the tomographic bin. Assuming a simple linear bias model relating ∆ g and the 3D matter overdensity ∆ M , the effective radial kernel for galaxy clustering is where b g is the linear galaxy bias, and p(z) is the redshift distribution normalized to unit integral.
In the case of eBOSS quasars, we will also take into account the impact of lensing magnification on the observed clustering. Magnification is caused by a combination of the displacement in galaxy positions and the modification in the observed source flux due to gravitational lensing. It is thus a direct tracer of the matter overdensity ∆ M with a radial kernel where s is the slope of the source magnitude distribution, and q L is the lensing kernel is the scale factor, and K is a scale-dependent factor that accounts for the difference between the three-dimensional transverse Laplacian, connecting the matter overdensity and the gravitational potential, and the angular Laplacian on the sphere: which is only significantly (above ∼ 1%) different from 1 on large scales ( 10). We will use s = 0.2 for the quasar sample [23]. The final kernel for the eBOSS quasars is given by the sum of q δg and q µg .
• γ G and κ are directly related to the three-dimensional matter overdensity ∆ M with radial kernels is the equivalent of K for the shear field, and χ * is the comoving distance of the source plane (i.e. the distance to the last scattering surface in the case of CMB lensing). We will also account for the impact of intrinsic alignments on the galaxy shear field, using the non-linear alignment model of [24,25]. In this model, the intrinsic alignment contribution to the observed shear, γ I , is proportional to the local tidal field, and therefore its radial kernel is proportional to the redshift distribution of the source sample: We will parametrize the intrinsic alignment amplitude A(z) as was done for the analysis of the DES first-year data [2,5]: with A IA,0 and η IA two free parameters, z 0 the redshift pivot (which we fix to z 0 = 0.62 as in [2,5]), and D(z) the linear growth factor. As in the case of the eBOSS quasars, the final kernel is given by the sum of q γ G (or q κ ) and q γ I . In our analysis we will also account for residual multiplicative biases in the shear power spectra. These enter the power spectra as an overall multiplicative factor (1 + m) for each shear field being correlated, with m a free parameter of the model.
As described above, the three tracers (galaxy clustering, cosmic shear and CMB lensing) used here can be used to measure the matter density fluctuations. The last remaining ingredient of the model therefore is the matter power spectrum P mm (k, z) entering Eq. 2.2. The model used here is described in the next section. We used the Core Cosmology Library 2 (CCL) [26] to compute all cosmological quantities, making use of the CLASS [27] Boltzmann code to calculate the linear matter power spectrum.

Growth reconstruction
Given the ongoing debate around the value of the amplitude of matter fluctuations at late time measured by large-scale structure data [2,[4][5][6]14] compared with that extrapolated by CMB data in the context of ΛCDM [1,28], our main objective in this paper is to explore: a) whether this tension can be directly mapped into the time evolution of inhomogeneities within the range of redshifts covered by existing data sets, and b) whether the growth histories recovered by different data sets are compatible with each other. To do so, we will decouple the linear growth factor D(z) from the background cosmological parameters and instead treat it as a free function that we will constrain directly from the data.
In a ΛCDM Universe with no massive neutrinos, in which the only relevant density inhomogeneities are those of pressureless matter, the linear matter overdensity field ∆ L M grows in a "self-similar" fashion, in which time dependence is factorizable: is the linear growth factor normalized as D(0) = 1. This result is easily understandable since, in this case, the equation for the evolution of matter overdensities (conservation of the energy momentum tensor) is scale independent. This factorizability then maps directly into the linear matter power spectrum: The dependence of the growth factor on redshift as a time variable can be a rich observable to constrain the main energy components of the Universe. Unfortunately, the non-linear nature of the gravitational interaction causes the matter overdensity to quickly depart from this linear behaviour, causing smaller, non-linear scales to grow faster than larger, linear ones.
Fortunately, for a wide range of cosmological models, the power spectrum of the observable non-linear overdensity can be expressed as a functional of the linear power spectrum: (2.13) with little dependence on the specific ingredients of the cosmological model beyond those that determine the form of P L and its evolution in time [29][30][31][32][33][34]. Here we will make use of one such parametrization to connect the linear and non-linear power spectra: the popular HALOFIT model as implemented in [33].
There are different options to parametrize D(z), such as expanding it as a linear combination of basis functions (e.g. as a polynomial in z or a) or, more ideally, modelling it as a Gaussian process with its hyperparameters and conditional distribution determined from the data. Existing off-the-shelf parameter inference frameworks for the cosmological analysis of large-scale structure data are not efficient enough yet to deal with the high-dimensional parameter spaces associated with Gaussian processes (although the community is moving fast in that direction [35,36]), and therefore we choose a middle ground. In our case, D(z) is determined by its value at a set of fixed redshift nodesD z i ≡ D(z i ). Each node is treated as an additional free parameter in the likelihood. The growth factor D(z) is then calculated at any z by interpolating over the values ofD z i , extrapolating beyond the range covered by the fixed redshift nodes z i . Specifically, we use a quadratic spline interpolation in the space (log(1 + z), log(D(z))). In order to avoid allowing for unphysically large or negative values of D(z) at high redshifts where our data have no constraining power, we fix D(z) to the ΛCDM prediction with the best-fit Planck cosmological parameters beyond z = 5 3 .
We choose four redshift nodes centered at the mean redshifts of some of the galaxy clustering tracers used in our analysis. The logic behind this choice is that, on large scales, a combination of galaxy clustering and weak lensing data can be used to effectively measure the galaxy bias and the amplitude of matter fluctuations at the mean redshift of the clustering sample, and thus the galaxy clustering bins act as natural anchors at which growth is effectively measured. In particular, we choose the position of the first, third and fifth redshift bins for the DES clustering sample (see Section 3.1.1), as well as the mean redshift of the eBOSS-QSO sample. The nodes are thus located at z i ∈ {0.24, 0.53, 0.83, 1.5}. As stated above, this is in addition to a fixed node at z = 5 that matches Planck's best fit ΛCDM cosmology.
There are a few caveats associated with the method chosen to quantify structure growth. First, the extrapolation above and below the redshift range covered by the nodes may lead to biases in the final constraints. This is more relevant in the range z < 0.24 which the cosmic shear data are sensitive to, rather than at high redshifts. Additionally, the choice of a quadratic spline instead of other interpolation methods will likely have an impact on the final constraints on the spline parameters (e.g. in terms of their final uncertainty or correlation between modes). As we discuss in Section 5.3, neither of these effects has a strong impact on our final result. Additionally, the self-similar linear growth assumed here (Eq. 2.11) is not valid in general (e.g. in the presence of massive neutrinos, or through scale-dependent growth in modified gravity models e.g. [37]). Current constraints on neutrino mass from particle physics experiments and CMB observations place any signature of neutrinos below the level of detectability of the data and scales used here, and therefore any significant departure from eBOSS DELS KiDS DES Figure 1: Sky footprint of the galaxy surveys used in this analysis. The Planck CMBκ map overlaps with all the data sets and is not shown. We carry out separate analyses of the DES (green) and KiDS +DELS (cyan and yellow) 3 × 2pt data, which we label SD and ND respectively. When analysed separately, we include their combination with Planck CMBκ and eBOSS-QSO (dark blue). We will also consider the combination of all the data sets shown.
ΛCDM in terms of growth history we observe is unlikely to be due to the scale-dependent signature of massive neutrinos. Finally, the use of HALOFIT to relate the linear and nonlinear power spectra, while sufficiently accurate for a wide family of cosmological models, is not guaranteed to be appropriate for models with an arbitrary growth history. Since the deviations with respect to ΛCDM found in our analysis are relatively mild, we do not expect this to have a significant impact on our conclusions.

Data
This work is based on the analysis of 6 different data sets. These are: the cosmic shear and galaxy clustering samples used in the cosmological analysis of the first-year data release of the Dark Energy Survey (DESg and DESγ respectively, [2], Section 3.1), the cosmic shear sample used in the fourth data release of the Kilo-Degree Survey (KiDS1000, [38], Section 3.4), a galaxy clustering sample extracted from the DESI Legacy Survey for the analysis presented in [19] (DELS, Section 3.5), the clustering of quasars in the extended Baryon Oscillation Spectroscopic Survey (eBOSS-QSO, [20], Section 3.2), and the lensing convergence of the CMB measured by Planck (CMBκ, [21], Section 3.3).
Our work will conceptually divide these data sets into two groups of data, determined by the data combinations for which a "3 × 2-point" analysis (i.e. the combination of twopoint correlations involving cosmic shear and galaxy clustering) can be carried out over  non-overlapping sky regions. First, the "South data set" (SD) will comprise the DESY1 galaxy clustering and cosmic shear samples. The "North data set" (ND) instead focuses on the combination of cosmic shear from KiDS1000 and galaxy clustering from DELS. The ND and SD sets also include the cross-correlation of their clustering and shear samples with the Planck convergence map in their respective footprints. Finally, when analyzed separately, both data sets also include auto-correlations of the eBOSS-QSO sample, and its cross-correlations with the CMB lensing map, in order to provide a high-redshift lever arm for the growth reconstruction. We will also consider the combination of all six data sets (which we will label FD). The specific auto-and cross-correlations between the different probes and their associated scale cuts are described in Section 4.1. Figure 1 shows the sky footprints covered by each of the data sets used here. The redshift distributions and associated radial kernels for each tracer used in our analysis (as defined in Eq. 2.1) are shown in Figure 2. In combination, the data used in this analysis allow us to cover the range of redshifts z 2 as well as a significant fraction of the celestial sphere.
The next subsections describe the procedure used to process these data sets and extract  Table 1: Summary table describing the methods used to generate the main map-level data products needed to estimate power spectra and covariance matrices (signal maps, masks and noise power spectra) for the 6 data sets used in this analysis. The cells reading "PDR" correspond to quantities that are directly provided in the public data release associated with these data.
four key data products: signal maps, sky masks, noise power spectra and redshift distributions. A summary is provided in Table 1.

DES Y1 galaxy clustering and weak lensing
The Dark Energy Survey is a 5 year survey that will cover 5000 deg 2 in 5 filter bands (grizY ) and has mapped hundreds of millions of galaxies and thousands of galaxy clusters [39]. These observations are taken from the Cerro Tololo Inter-American Observatory (CTIO) with the 4 m Blanco Telescope, using the 570-megapixel Dark Energy Camera (DECam [40]). In this paper we use the first year of data products 4 , which cover 1786 deg 2 before masking [41,42]. In this analysis we use the same fiducial galaxy samples used in the DES Y1 3 × 2pt analysis [2], including their associated redshift distributions. Our fiducial analysis will also employ the same models used in [2] to describe various systematic effects, as discussed in Section 2.1.

Galaxy clustering
We use the clustering sample presented in [43]. The sample was constructed using the red-MaGiC algorithm, which selects red luminous galaxies with excellent photometric redshift accuracy (σ z = 0.017(1+z) in this case, [44]). We split this sample into the same five redshift bins used in the analysis of [2,43,43]. The three lower redshift bins are populated with galaxies from the redMaGiC high-density sample. These have a comoving densityρ n 10 −3 and minimum luminosity, L min 0.5L * . The other two redshift bins contain galaxies from the redMaGiC high-luminosity (ρ n 4 × 10 −4 , L min = L * ) and higher-luminosity (ρ n = 10 −4 , L min = 1.5L * ) samples respectively.
To track the survey geometry, we use the mask publicly available with the Y1 release, in the form of a HEALPix 5 [45] map at resolution N side = 4096 containing the effective fractional area w p of each pixel p. In order to avoid inaccuracies from strongly masked pixels, we set the value of all pixels with w p < 0.5 to zero. This leaves a total unmasked area of 1321 deg 2 .
All galaxies in the redMaGiC tomographic bins are assigned weights to correct for the impact of observational systematics. In order to create the associated overdensity map for each redshift bin we must therefore use both these weights and the geometric information contained in the sky mask. To do so we first create a map containing the weighted number of galaxies lying in a given pixel p, n p = i∈p v i , where v i the weight of the i-th galaxy. The overdensity field is then given by δ p = n p nw p − 1, (3.1) wheren = p n p / p w p is the weighted average galaxy number per pixel and w p the value of the mask on that pixel. The auto-correlation galaxy clustering power spectra has a shot noise contribution, N , that has to be subtracted. Assuming perfect Poisson sampling, this can be done analytically [46]. In short, the "mode-coupled" noise power spectrum (i.e. before multiplying by the inverse mode-coupling matrix in Eq. 4.3) is given bỹ where w is the mean value of the mask across the full sky, andn 2 is the effective mean number density, given byn Here, A pix is the pixel area in steradians. Note thatn 2 reduces ton in the case of equal weights. The DES Y1 release provides an estimate of the redshift distribution for the five clustering redshift bins. We use these in our analysis marginalizing over a parameter corresponding to the mean of each distribution.

Galaxy shear
The DES shear analysis was carried out using two different shape-measurement algorithms, IM3SHAPE [47] and Metacalibration [48]. In this work we make use of the Metacalibration catalog. Metacalibration fits a 2D Gaussian model for each galaxy to the pixel data in the r, i and z bands, convolved with their corresponding point-spread function (PSF). This process is repeated with artificially sheared images to calibrate the shear estimator, which allows for the calculation of shear-dependent selection effects that could bias the statistics a few percents [5,[48][49][50].
The raw ellipticities measured by Metacalibration,ê i , must be corrected as where e i is the i-th component of the resulting calibrated ellipticity. The multiplicative correction isR = (R 11 + R 22 )/2, where R is the 2 × 2 response tensor calculated by Metacalibration. The response tensor contains two additive terms: the shear response R γ , obtained from artificially sheared images as described above, and the selection response R S , which accounts for the selection bias that appears when applying a set of selection criteria on the sheared galaxy sample. The total response tensor is R = R γ + R S which is almost diagonal with R 11 ∼ R 22 , and thus its effect can be well approximated byR. Note that Metacalibration computes the shear response tensor R γ for each galaxy but, as done in [5,49], we average it over the whole sample in each redshift bin to calculateR. Finally, the DES analysis found a non-negligible mean residual ellipticity (ē i ∼ O(10 −4 )) in each redshift bin. Following [5], we subtracted this mean ellipticity per redshift bin after calibration.
We generate maps of the two shear components γ 1,2 as the per-pixel weighted average of galaxy ellipticities where i = 1, 2 is the shear component, v n is the weight associated with the n-th galaxy.
Metacalibration assigns unit weights to all galaxies v n = 1 in it. As discussed in [51], we use the sum of weights in each pixel as the mask associated with the resulting shear maps. This should be a close-to-optimal choice assuming the weights are close to inverse-variance: The mode-coupled noise power spectrum is calculated analytically as [51]: where σ 2 e,i = (e 2 i,1 + e 2 i,2 )/2 is an estimate of the shape-measurement noise rms per galaxy. As shown in [51], this estimate is equivalent to averaging over the power spectra of a large number of catalogs with randomly rotated ellipticities, at a much lower computational cost. This approach is exact in the noise-dominated regime, where the contribution from the cosmic shear signal to each galaxy is negligible. Note that, as a spin-2 field, the power spectrum of the shear field is zero at < 2.
The Metacalibration sample is divided into 4 tomographic redshift bins. We use the fiducial redshift distributions for each bin provided with the Y1 data release. These were estimated by stacking the per-galaxy probability distributions derived by the BPZ [52] photoz algorithm for all sources in a given bin [53]. They were further validated by cross-matching against the COSMOS 30-band catalog [54], and via cross-correlations (see [53] for further details). To a large extent, most of the uncertainty on the redshift distribution for cosmic shear samples can be well-described by an uncertainty in the mean of said distribution, since a shift in this mean impacts the width of the associated redshift kernel significantly [55]. Thus, as done in [5], we marginalize over four parameters corresponding to linear shifts in the mean of the redshift distributions of each bin.

eBOSS quasars
In order to extend the range of redshifts over which we reconstruct the growth history, we use quasar (QSO) clustering measurements from the extended Baryon Oscillation Spectroscopic Survey (eBOSS) derived the from Sloan Digital Sky Survey Data Release 16 (DR16). In particular we use the homogeneous quasar sample used for the cosmological power spectrum analysis of [20,56], and presented in [57]. The catalog comprises 343,708 objects with measured redshifts in the range 0.8 ≤ z ≤ 2.2, covering over 4,800 deg 2 . Note that we will not use the existing redshift-space distortion measurements from eBOSS-QSO to constrain growth here. Instead, we use the quasar sample as another two-dimensional projected tracer of the large-scale structure which we correlate with all the other datasets.
The footprint of this sample is described by a set of random sources, covering the same area in the absence of clustering. Furthermore, in order to correct for the modulation in the observed number density of objects caused by various systematics, observing conditions and Galactic systematics, objects in both the random and data catalogs are assigned weights. Since we carry out a projected 2D analysis, we include all systematic weights (accounting for redshift failures, fiber collisions, sky systematics), but omit the so-called "FKP" weights that maximize the signal-to-noise ratio of the three-dimensional power spectrum.
We combine eBOSS measurements from the North and South Galactic Caps into one single catalog that we split into two different bins with redshifts above and below z = 1.5 respectively. The redshift distribution of each bin is estimated directly from the data as a histogram of the measured spectroscopic redshifts. We do not account for any systematic uncertainty in the redshift distribution thus constructed.
In order to calculate the quasar overdensity map we make use of the random catalog to track the survey geometry. Specifically, the overdensity in pixel p is calculated as where v d,i and v r,i are the data and random weights for the i-th object lying in pixel p. α = p v d,p / p v r,p accounts for the fact that the random catalog is significantly larger than the data to minimize the impact of its associated shot noise. w p is the survey mask, which we compute as the scaled sum of random weights in each pixel: In order to minimize the impact of shot noise in the random catalog when constructing this mask, we do so at a relatively low resolution (N side = 512), ensuring a sufficiently high average number density of random points in non-empty pixels. This is then upgraded to our target resolution (N side = 4096) correcting for the different pixel area. The eBOSS-QSO catalog is by far the sparsest of the clustering samples used in this analysis, and the angular quasar power spectrum is dominated by shot noise over a large range of scales. A careful treatment of the noise bias in the auto-correlation is therefore crucial in order to obtain reliable constraints from it. This is further complicated by the fact that we need to account for the impact of shot noise in the random catalog, which affects both the numerator and denominator of Eq. 3.8. Fortunately, since the auto-correlations are noise-dominated on small scales, we can obtain a reasonable estimate of the noise bias from the data. An uncorrelated noise component would appear in the mode-coupled power spectrum as a scale-independent contribution, which we can thus estimate by averaging the value of the power spectrum calculated from the masked overdensity map in the range 2000 ≤ < 2N side = 8192. We find this method to be a better approximation than analytically calculating the noise contribution via e.g. a generalization of Eq. 3.2 accounting for data and random weights. Nevertheless, given the importance of the shot-noise contribution, we additionally marginalize over a constant noise power spectrum with a free amplitude and 10% Gaussian prior centered on the power spectrum amplitude estimated as we describe in Section 4.3.

CMB lensing from Planck
The CMB weak lensing convergence field (CMBκ) is produced by the matter overdensities between the last scattering surface and us, with most of the contribution coming from 0.5 z 3 [58,59]. In this work, we use the convergence map made available as part of the Planck 2018 data release [21]. The map covers a sky fraction f sky = 0.671 and thus overlaps spatially with all surveys considered here. Various photometric redshift surveys have made use of cross-correlations with CMB lensing maps from different collaborations to extract cosmological constraints [7,19,28,[60][61][62][63].
Planck is a third-generation space mission, following COBE and WMAP, dedicated to measure the CMB anistropies, offering full sky maps of temperature and polarization anisotropies with micro-Kelvin sensitivity per resolution element. The specifications of the 2018 Planck CMB lensing data release are described in detail in [21].
In this analysis we make use of the "Minimum Variance" (MV) lensing convergence harmonic coefficients, which we transform into a HEALPix map at N side = 4096 resolution. The harmonic coefficients are provided in the range < 4096. Following [21,64] we remove the smallest multipoles < 8, which are too sensitive to the mean-field subtraction in the lensing reconstruction process, and go up to = 2000, falling inside the aggressive scale ranges of [21]. We use the binary sky mask made available with this map. Since the lensing reconstruction noise power spectrum rises sharply with on small scales, we apodize this mask with a 0.2 • "C1" kernel [65] in order to minimize the leakage from noise-dominated small-scale modes. The final usable sky fraction is f sky 0.66.
Our analysis does not include the CMB lensing auto-correlation, and therefore we do not need to include a rigorous modelling of the various noise bias terms that enter the CMB lensing likelihood. The CMB lensing noise, however, enters the covariance matrix for any power spectra involving the CMB lensing map. For this, we use the estimate of the noise power spectrum N κ provided with the data release.

KiDS-1000 weak lensing
The Kilo Degree Survey (KiDS) is a large optical survey that has mapped 1350 deg 2 of the sky in four tomographic bands (ugri) using the VLT Survey Telescope (VST) located in the ESO Paranal Observatory. Its main objective is the measurement of the cosmic shear signal. In this paper, we use the Gold Sample from the public data release 4 (DR4) [66]. This data release includes both images covering a total area of 1006 deg 2 , which reduces to 777.4 deg 2 after masking, as well as forced photometry data in five additional infrared bands from the VIKING survey. The Gold Sample is described in [17,66] and it targeted a sample of galaxies with reliable shapes and redshift distributions.
We follow the analysis of [38]. We split the gold sample in 5 tomographic bins assigning each galaxy to a bin based on its best-fitting photometric redshift, z B . As described in [17], in each redshift bin we subtract the residual weighted mean ellipticity, and correct for the multiplicative bias factors listed in Table 1 of [38], estimated by KiDS from image simulations. We then produce shear maps and weight maps (mask) for each redshift bin following Eqs. 3.5 and 3.6 respectively, making use of the shear measurement weights assigned by the lensfit algorithm. The coupled noise power spectrum is estimated analytically as in Eq. 3.7.
We use the redshift distributions provided with the DR4. These redshift distributions were constructed using the self-organizing map (SOM) method described in [67,68]. As in the case of DES, we will marginalize over the mean of each redshift distribution in the likelihood.
Our analysis will make use of the galaxy sample selected by [19] to obtain cosmological constraints from its cross-correlation with CMB temperature and lensing convergence 6 . We also follow their analysis choices closely. Each galaxy was assigned a redshift based on a multi-dimensional matching in colour space with a set of spectroscopic samples. We use these redshifts to separate the sample into the four tomographic bins used in [19].
Our main objective in using these data is to combine with the KiDS cosmic shear data to carry out a 3×2pt analysis that can then be compared with the results found with the DES Y1 data. In order to facilitate this comparison, as well as the combination of both data sets, we remove all data with declination δ < −36 • from the DELS sample, ensuring no area overlap between DELS +KiDS and DES.
For each redshift bin, we compute a first estimate of the galaxy overdensity map using Eq. 3.1, where n p is the number of galaxies in pixel p,n is the mean number density, and w p is the completeness of the pixel (understood as its effective fractional area).n is estimated as n = p∈G n p / p∈G w p , where G is the set of "good" pixels with completeness above 95% and a star density lower than N star = 8515 deg −2 [19]. We make use of the completeness map and star map made available by [19]. The latter corresponds to a smoothed version of the ALLWISE total density map [75]. Additionally, all pixels with completeness w p < 0.86 or star density N star > 1.29 × 10 4 deg −2 were masked. The overdensity field thus created contains residual star contamination, affecting the galaxy auto-correlation on large scales. We correct for this at the map level by subtracting a systematic overdensity map whose pixel values are estimated by evaluating a 5th-order polynomial fit to the mean galaxy overdensity with respect to the local logarithmic number of stars in each pixel, as done in [19]. As in the case of DES, we estimate the coupled noise power spectrum through Eq. 3.2.
We model the redshift distribution of each tomographic bin using the same approach described in [19]. The true redshift distribution can be related to the photo-z distribution through a convolution with the conditional distribution p(z t |z p ), where z t and z p are the true and photometric redshifts respectively. This conditional distribution is modeled to be stationary (i.e. only dependent on z t − z p ) and given by a Lorentzian distribution of the form (3.10) The parameters x 0 , σ and a are determined from the existing spectroscopic data and given in Table 1 of [19]. Additionally, the authors of [19] marginalized over x 0 and a in their cosmological analysis, showing that they can be self-calibrated to a large extent through the use of clustering cross-correlations between different redshift bins. In order to keep a consistent model among different probes, we simply marginalize over linear shifts in the mean of each redshift distribution, and fix x 0 and a to the values found by [19] after self-calibration.

Power spectra
We analyze all data using a common harmonic-space power spectrum framework based on the pseudo-C method as implemented in NaMaster [65], including the approximations described in [51,76] to estimate the power spectrum covariance. We describe the method briefly below and direct the reader to these references for further details on the implementation 7 . An observed mapã is modelled as the product of the true map a and a weights map w: Although, for simplicity, we will limit our discussion to scalar (spin-0) fields such as δ g or κ, the methods used are directly generalized to spin-2 (or, in fact arbitrary spin) fields such as γ. This is described in detail in [51,65,76]. The weights map can be understood both as a mask, i.e., a map defining whether a given pixel has been observed or not (w = 1 and 0 respectively), as well as an inverse-variance local weight, down-weighting regions of high noise.
Following the convolution theorem, the spherical harmonic coefficients ofã are a convolution of the spherical harmonic coefficients of a and w, a fact that leads to mode coupling between the power spectra of the observed and true fields. Defining the "coupled pseudo-C " C ab between fields a and b with weight maps v and w respectively as its relation to the true underlying C ab is given by where denotes averaging over realizations of a and b. M vw is the so-called mode-coupling matrix (MCM), and depends solely on the weight maps of the two fields being correlated. As shown in [77], M vw can be computed efficiently and analytically thanks to the orthogonality of the Wigner 3j symbols. Explicit expressions for the coupling matrix of different combinations of spin-0 and 2 fields can be found in [65]. In the limit of full-sky data, M vw is the identity, and departures from this limit give rise to a statistical off-diagonal coupling between neighbouring multipoles. In many practical cases the MCM is non-invertible, and thusC ab cannot be turned into an unbiased estimator of C ab . The pseudo-C estimator then proceeds along the following three steps: 1. Binning. As a way to regularize the MCM, we bin the power spectra into bandpowers containing weighted sums of different s: We can then define the binned MCM which is usually invertible for sufficiently broad bandpowers.
2. Inverting the MCM. The decoupled bandpowersĈ ab q are then defined by inverting the binned MCM: where we have explicitly removed the (binned) mode-coupled noise power spectrum N ab q . The form ofÑ ab q depends on the maps being correlated, and their calculation has been described in Section 3 for the different data sets used here.
3. Bandpower convolution.Ĉ ab q would be an unbiased estimator of the true power spectrum evaluated at, e.g., the central multipole of each bandpower, if the latter was exactly constant within each bandpower. Since this is not the case, the effects of binning must be propagated through, in order to connect the observed bandpowers with a given theoretical prediction 8 . The theoretical prediction for the bandpowers C ab q is related to the theory power spectrum C ab by convolving the latter with the bandpower window functions F vw q through a fast matrix-vector multiplication where F ab encodes the three pseudo-C linear operations: mode-coupling, binning, and binned mode-decoupling.
We use a common HEALPix resolution of N side = 4,096 for all sky maps, which allows us to compute all power spectra up to a maximum multipole max = 3N side − 1 = 12,287. We bin all power spectra into a common set of bandpowers with the following binning scheme: we use linear bins with width ∆ = 30 between 0 ≤ ≤ 240, and logarithmic bins until max with ∆ log 10 ( ) = 0.055. The resulting bandpower edges are listed in Table 2 of [51].

Covariance matrix
The power spectrum covariance matrix can be decomposed in three contributions [78,79]: the disconnected "Gaussian" part, equivalent to assuming that all fields involved are Gaussianly distributed, the connected non-Gaussian part (cNG), corresponding to the intrinsic connected trispectrum of the fields, and the super-sample covariance (SSC), caused by the coupling of small-scale modes sourced by modes larger than the survey footprint. As we show in appendix A, the cNG and SSC contributions are subdominant on the scales used here, and our fiducial analysis will employ only the Gaussian contribution.
The Gaussian covariance can be computed analytically and depends on the survey geometry, which causes statistical correlations between different bandpowers. Accounting for this effect is particularly critical for the case of cosmic shear, given the high complexity of its associated weights map (which follows the distribution of the source galaxies). However, the exact calculation of these correlations sales as O( 6 max ), making it intractable for our data. In this work, instead, we use the Narrow Kernel Approximation (NKA) presented in [76,80], and improved in [51] for the case of cosmic shear. The method scales as O( 3 max ) and has been shown to be accurate for all spin-0 and spin-2 quantities used here. The core assumption behind the NKA is that the map-level mode-coupling matrix is close to diagonal, such that the power spectra can be treated as constant in the calculation. This is usually an excellent approximation for well-behaved masks, but fails catastrophically in the case of cosmic shear, or for steep power spectra. As shown in [51], the NKA can be improved significantly in these cases if one simply substitutes the power spectra used in the calculation by their mode-coupled versions scaled by the overlapping sky fraction: If the noise properties vary across the footprint, the effective masks of the signal and noise components of the maps are not the same, and thus the signal-signal, signal-noise and noise-noise contributions to the covariance matrix should in principle be computed using different mode-coupling coefficients [81]. As shown in [51], this additional complication can be avoided by simply adding to the signal power spectrum (Eq. 4.9), the equivalent coupled noise power spectrum scaled by the overlapping sky fraction. We have described the calculation of the coupled noise power spectrum for the different data sets used here in Section 3. While this is not an exact result, the impact of this approximation on the covariance matrix is negligible for our current analysis (although its validity must be reassessed in for future, more sensitive cosmic shear data sets). As a technical note, we find that the "spin-0" approximation discussed in [76], which treats the E and B components as independent scalar fields, yields a better estimate of the shear covariance matrix for the KiDS data, given its significant depth variations across the footprint, and thus we make use of this approximation here. In addition to this, we make use of the "Toeplitz approximation" proposed in [82] to accelerate the calculation of the mode-coupling coefficients used in the covariance calculation. The impact of this approximation is at the sub-percent level for our data.
This calculation of the covariance matrix requires an estimate of the power spectra of the different fields. We calculate these using CCL assuming the Planck 2018 best-fit cosmological parameters: ( . We verified that these parameters are a reasonable fit to the measured power spectra on the scales used in this analysis. For the CMBκ autocorrelation, we use the signal+noise power spectra provided with the public Planck data release.

Likelihood
We use a Gaussian power spectrum likelihood to derive constraints on the free parameters describing the measured power spectra. To do so we made use of the MontePython sampler [83,84], which we modified to interact with CCL as the main code for theory calculations 9 . We sample with the Metropolis-Hasting algorithm [85,86] and address the chains convergence requiring the Gelman-Rubin parameter R − 1 0.01 [87].

Parameter
Prior 0.2116 Table 2: Prior distributions for the cosmological parameters. When reconstructing the growth history, we fix the normalization of the linear power spectrum template to σ fid 8 = 0.8111 and do not vary A s . This is taken into account in our growth reconstruction described in Section 2.2. Furthermore, we fix the highest redshift node,D 5 , at its value for the best fit of Planck (Eq. 4.10), ensuring that we recover our fiducial growth at z ≥ 5. U (a, b) and N (µ, σ) are a uniform distribution with boundaries (a, b), and a Gaussian distribution with mean µ and variance σ, respectively. The priors largely follow the choices of [60].
We consider two types of cosmological models: i.e. the amplitude of scalar perturbations, the primordial spectral index, the present value of the matter and baryonic density parameters, and the dimensionless Hubble parameter, respectively. The priors used for each parameter mostly follow the choices made by the DES collaboration [60], and are listed in Table 2. Four particular choices must be noted. We fix the optical depth τ = 0.08, since our data are insensitive to its value, and the redshift of the last scattering surface to a fiducial value z * = 1100. We impose a flat prior on A s , as opposed to other common choices such as log A s , σ 8 or S 8 . The impact of this choice on the final constraints on S 8 is discussed in detail in [88]. Finally, we consider only cosmologies with massless neutrinos. Since these data are not able to place strong constraints on the sum of neutrino masses, this mostly allows us to accelerate the calculation of the matter power spectrum. Although these choices may have an effect on our final constraints, at the level of a few fractions of a σ, we emphasize that our ΛCDM constraints are mainly aimed at validating our analysis pipeline by comparing them with those found in the literature for subsets of our data, rather than performing a thorough study of ΛCDM and its extensions, which we leave for future work.
• Growth reconstruction. In this case we retain 4 ΛCDM cosmological parameters: Ω m , Ω b , h and n s . The amplitude of matter fluctuations is defined by the free growth factor parametersD z (see Section 2.2), which define the value of the linear growth factors at redshifts z = 0.24, 0.53, 0.83, and 1.5. Note that the overall normalization of theD z parameters is degenerate with that of the fiducial linear power spectrum at z = 0 used in Eq. 2.12. Thus, when generating P L (k, 0) in that equation, we fix its normalization to the best-fit value of σ 8 measured by Planck, σ fid 8 = 0.8111. The actual value of σ 8 for that particular model is then given by σ 8 = D(0)σ fid 8 . We are able to constrain the value of all theD z nodes for most of the data combinations explored. However, in order to avoid unphysical (e.g. negative) values of D i in cases where some of the parameters are unconstrained (e.g. in the absence of high-redshift data), we impose a flat prior on the growth parametersD z ∈ [0, 2.5]. This prior is uninformative and sufficiently broad not to bias our result when the data are able to constrain these parameters.
Our model also contains a large number of free nuisance parameters characterizing different sources of astrophysical uncertainties and systematics. These are (see Table 3 for their priors): • Galaxy bias. We use a linear bias model, assigning a different bias parameter for each galaxy clustering sample and redshift bin (i.e. 5, 4 and 2 parameters for DESg, DELS and eBOSS-QSO respectively).
• Intrinsic alignments. We use an evolving alignment amplitude (Eq. 2.10), with free amplitude and redshift evolution parameters A IA,0 , η IA .
• Photo-zs. We characterize the uncertainties on the redshift distributions of the different galaxy samples in terms of a shift in the mean redshift ∆z. We marginalize over one such parameter in each redshift bin for both clustering and shear samples (except in the case of eBOSS-QSO), with the same Gaussian prior used by the corresponding collaborations. Note that the analysis of the DELS sample by [19] used a different parametrization for the redshift distribution uncertainties, stated in terms of the conditional photo-z distribution. The relatively broad priors on ∆z we use should encompass the small uncertainties found by [19] after using cross-bin correlations to self-calibrate the redshift distributions, although a more thorough analysis of all possible modes of uncertainty (e.g. in the width of the distribution [46]) would be desirable both in the case of DES and DELS.
• Quasar shot noise. As noted in Section 3.2, our estimate of the noise bias for the eBOSS-QSO sample is not exact. Therefore we marginalize over two parameters characterizing the amplitude of the noise power spectrum for the two quasar redshift bins with a 10% Gaussian prior. Since these parameters are linear in the power spectra, this marginalization can be done analytically by modifying the covariance matrix of the eBOSS-QSO auto-correlations.
• Multiplicative bias. We marginalize over a multiplicative bias parameter m i in each cosmic shear redshift bin, with Gaussian priors derived by DES and KiDS. Note that, as done in the KiDS analysis [38], it is possible to marginalize over these parameters analytically by linearizing the impact of their uncertainty on the power spectra [89]. We verified that our constraints are not sensitive to the choice of marginalizing over these parameters exactly in the likelihood instead, and thus we do so for consistency with the DES analysis.

Nuiscance parameters priors
Parameter Prior Parameter Prior  Table 3: Prior distributions for the nuisance parameters entering our analysis for each tracer. U (a, b) and N (µ, σ) describe a uniform distribution with boundaries (a, b) and a Gaussian distribution with mean µ and variance σ, respectively. The index i in b i g and m i runs over the different redshift bins. The DES and intrinsic alignment priors have been taken from Table 1 of [60] and the KiDS priors follow [38].

DES lens photo-z bias eBOSS QSO bias
Following the DES Y1 analysis we do not model the impact of baryonic effects on the matter power spectrum. Although this should have a subdominant effect on the scales used here [4,6], a more rigorous study of the cosmological constraints on ΛCDM parameters from the combination of all data sets studied here would require a more careful assessment of this source of uncertainty.
Finally, our data vector will contain a combinaion of the following power spectra: • Clustering auto-correlations. We will only consider auto-correlations between galaxies in the same redshift bins. Although the cosmological signal is concentrated in these auto-correlations, cross-bin correlations can be used to self-calibrate photometric uncertainties [19,46]. Nevertheless, we choose to discard them in order to mimic the choices made in the DES analysis [2]. We impose strict scale cuts on clustering, using multipoles smaller than max = k max /χ, whereχ is the comoving radial distance  Table 4: Scale cuts for the different tracers in our data vector. Scales without units are angular modes; whereas those in Mpc −1 are comoving wavenumbers k, translated into angular multipoles as = kχ, whereχ is the comoving distance at the center of each redshift bin. We exclude bandpowers with effective s lower or larger than the scale cuts.
at the mean redshift of the bin, and k max = 0.15 Mpc −1 . This is done in order to ensure the validity of the linear bias model used here. In addition to this, we will impose large-scale cuts to avoid the impact of large-scale observational systematics in the galaxy auto-correlations. Following [19,20], these correspond to k min = 0.02 Mpc −1 for eBOSS-QSO, and min = 30 for the DELS sample. As done in [43] we do not include any large-scale cut on the redMaGiC sample (and a visual inspection of the power spectra did not reveal any obvious sign of additional large-scale power due to systematics).
• Clustering-lensing cross-correlations. We will use all available cross-correlations between clustering redshift bins and lensing probes, including both tomographic cosmic shear samples and CMB lensing. The only exception to this is the cross-correlation between eBOSS-QSO and any of the cosmic shear samples, given their null spatial overlap (see Figure 1). In each cross-spectrum we impose the high-scale cut of the corresponding galaxy clustering sample as described above.
• Shear-shear correlations. We use all available cross-correlations between different tomographic bins corresponding to the same cosmic shear sample (ignoring DES-KiDS cross-correlations given their zero area overlap). For the DES power spectra we use all bandpowers in the range 30 < < 2000, which were shown in [51] to be free from systematics. For KiDS we use a more strict large-scale cut, following [4,38], with 100 < < 2000.
• Shear-CMBκ correlations. We use all cross-correlations between the Planck CMBκ map and the KiDS and DES samples, using the same scale cuts used for the analysis of the corresponding shear-shear correlations.
The choice of scale cuts used in our analysis is summarized in Table 4. Note that our data vector does not include the auto-correlation of the Planck CMBκ map. Although this would provide a valuable constraint on the integrated evolution of cosmic structures over the range of redshifts studied, we choose to exclude it in order to avoid the complications of modelling the different cosmology-dependent lensing biases [21,90] as well as accurately  Table 4 have been applied as they enter into the MCMC. The yellow points correspond to eBOSS-QSO auto-correlations, and are caused by the analytical noise marginalization. We see that the Gaussian covariance is mainly described by its diagonal elements, although there are some non-negligible correlations between nearby elements. describing the covariance of this power spectrum with all other components of our data vector, given the spatial overlap of the lensing map with all our data sets.
The resulting data vector contains 665 and 662 elements for the SD and ND data sets respectively, and 1275 elements in the case of the full combination. Figure 3 shows the correlation matrices for the ND and SD data sets from their Gaussian covariance matrices after applying the scale cuts.

Validation
Since our constraints will be based on a re-analysis of the different data sets presented in Section 3 using a common harmonic-space framework, we carry out a set of basic validation tests to ensure the robustness of the resulting power spectra and covariance matrices.
NaMaster, the pseudo-C power spectrum estimator we use, has been extensively validated [65,76], and the specific case of cosmic shear data, given the significantly more complex survey geometry, was studied in detail in [51]. The main object of our validation is therefore to diagnose the presence of systematics in the maps that could affect the estimated power spectra.
Since the cosmological shear signal is dominated by a pure E mode (although see e.g. [91] ), we study the significance of any power spectra involving shear B-modes as a signature of systematic contamination. We do so by calculating, for each such power spectrum, the probability-to-exceed (p-value) of its χ 2 with respect to a null signal. The results are shown in Table 5 for all the combinations involving shear tracers explored here. In the vast majority of cases the resulting p-values are acceptable, showing no evidence of B-mode contamination at more than 2σ (p > 0.05). Although we find 4 cases with p < 0.05, such a small number  is compatible with the look-elsewhere effect. We quantify this by conducting a Kolmogorov-Smirnov test on the full set of χ 2 values with respect to a χ 2 distribution, finding probabilities p = 0.72 and p = 0.35 for the ND and SD data sets respectively. Therefore we find no evidence of B-modes in either of the shear samples within the range of scales used here. Another source of map-level systematics is the modification of the observed number density of galaxies due to various observational systematics (e.g. dust absorption, star contamination, completeness variations). The impact of these systematics has, to some extent, already been taken into account in the three galaxy clustering samples used in this analysis. The DES redMaGiC galaxies and eBOSS-QSO sample have galaxy weights associated to them aimed at correcting for the effect of known systematics, and the overdensity maps for the DELS sample used here are corrected by completeness variations as well as a fifth-order polynomial of the local star density [19]. Additionally, we implement large-scale cuts on the auto-correlations of the eBOSS-QSO and DELS data sets following the prescription presented in [20] and [19] respectively, removing scales significantly affected by these systematics (see Table 4). In order to test for the impact of residual contamination, we recompute these power spectra making use of the linear deprojection method implemented in NaMaster. The method removes the impact of specific systematics by projecting the data on the subspace orthogonal to a set of known contaminant maps, correcting for the resulting bias to the power spectrum analytically (see [65,92]). We then compare the resulting power spectra with those calculated without deprojection and calculate their relative χ 2 in order to detect significant deviations. For eBOSS-QSO we deproject a set of 15 systematic templates corresponding to observing conditions in the SDSS survey as well as a dust template based on [93]. For DELS, we deproject the completeness and star density maps made publicly available by [19]. For the redMaGiC sample, we deproject the set of 13 most relevant survey property maps identified in [43]. In all cases we do not observe any significant deviation within the range of scales used here, with a relative reduced χ 2 below 0.1 in all cases.
The methods used to estimate the Gaussian covariance matrix have been thoroughly validated in [51,76] for the type of data used here. This, together with the acceptable χ 2 values we find in Section 5.2, as well as a visual inspection of the scatter in the power spectrum residuals with respect to the ΛCDM best-fit prediction (see Section 5.2) reassure us that the covariance matrix used here is not significantly over-or under-estimated.

ΛCDM constraints
In this Section we discuss how these data sets constrain the ΛCDM model. On the one hand, this is a validating exercise through which we show that our results are in agreement with those found by others in the literature using different combinations of our data. On the other hand, it also allows us to present ΛCDM constraints from a novel combination of data sets. Figure 4 shows the constraints on Ω m and S 8 for the 3 × 2pt analysis of the DES and KiDS +DELS samples (blue and pink contours, respectively), together with the results found in the official 3 × 2pt analysis published by the DES and KiDS collaborations. We find an overall good agreement between both sets of constraints, although our analysis recovers visibly larger uncertainties in both cases. These differences are understandable, since the analyses are not equivalent. The official DES analysis [2] made use of real-space correlation functions, making it difficult to match their choice of scale cuts. The KiDS 3×2pt analysis [4] made use of spectroscopic clustering samples (BOSS and 2dFLenS), as well as a variety of two-point function estimators. Furthermore, unlike these official analyses, we considered only massless neutrinos and did not marginalize over their mass. The poorer constraints we find on Ω m in comparison with the KiDS analysis are expected, since our use of photometric clustering prevents us from taking advantage of the BAOs as a standard ruler. As shown in [94] this could be improved significantly by increasing the range of scales over which galaxy clustering can be used through a perturbative bias expansion. The method used here to estimate the power spectrum covariance matrix, as well as the choice of parameter priors, also differ from those used by both collaborations. Nevertheless, the comparison with these published results for the ΛCDM model, together with the validation tests described in the previous section, does not reveal any significant issues in our analysis.
The left panel of Figure 5 shows in blue. The right panel of Figure 5 then breaks these constraints down by tracer combination, showing the constraints found in the absence of shear (gray) and galaxy clustering (green) in addition to the full data set (red) and Planck (blue). These constraints on S 8 and Ω m are also listed in Table 6 for the different experiment combinations explored here, as well as pictorially represented in Figure 6. Although in most cases we find an increasing tension with the value of S 8 , it is interesting to note that this tension is driven by the shear data, and is not evident through the combination of galaxy clustering and CMB lensing. Note, however, that this tension has also been reported from this probe combination (although using different data sets) by other groups [7,19]. In spite of this tension, the different data sets used here are in reasonable agreement with each other. Since there is no obvious sign of tension between them, we combine them to find a constraint on S 8 given by Compared with the constraints found by Planck on this parameter, S Planck 8 = 0.832 ± 0.013, and assuming Gaussian errors added in quadrature, the level of tension is ∼ 3.4σ. This is 0.4σ larger than the tension found by the KiDS collaboration [4,14], and in agreement with previous results. As noted in [4,95] tension between experiments in their measurement of one particular parameter is not necessarily indicative of tension between their data sets when the full multi-dimensional parameter space is taken into account. Nevertheless, since similar constraints have been consistently obtained by various groups, using different data sets, this parameter tension must be analyzed further. Although more insight is expected from the ongoing analysis of new data from DES and other collaborations, this motivates our reconstruction of the growth history as a means to better understand the origin of this tension. Table 6 also lists the values of the χ 2 as a metric for goodness of fit for each data set combination. In each case we calculate the corresponding p-value assuming an effective number of degrees of freedom given by where N d is the total number of data points, and N b is the number of free galaxy bias parameters. See Figure 6 for a visual representation of these results.  Table 6 for a quantitative representation of these results. of data points, and N b is the number of free galaxy bias parameters. This was found to be a reasonable rule of thumb in [94], and we expect it to be still the case since only two cosmological parameters, Ω m and S 8 , are completely constrained by the data. In all cases we find reasonable p-values above 10%. The ΛCDM model therefore is able to describe the data reasonably well, and we find no signatures of strong tension between data sets. Note that achieving this satisfactory goodness of fit is remarkable given the difficulties that have traditionally been involved in the calculation of the covariance matrix for two-point correlations [2,38,96,97]. This is likely the result of a careful accounting for the impact of survey geometry on all signal and noise contributions to the covariance matrix combined with the relatively simpler covariance calculation needed in a harmonic-space pipeline. In order to ensure that this is not the result of an over-estimation of the power spectrum uncertainties, we visually inspect all the power spectrum residuals with respect to the bestfit ΛCDM model using all data, revealing no obvious trend or signs of error overestimation. More quantitatively, a Komogorov-Smirnov test reveals that the distribution of the residuals normalized by their estimated standard deviation agrees well (p = 0.04) with a standard normal distribution with zero mean and unit variance (see also Figure 7). Thus we find no evidence of misestimation in the power spectrum uncertainties.

Growth reconstruction
We now move on to the constraints on the linear growth of structure found from the different data combinations explored here. As a reminder, we model the redshift-dependence of the linear growth factor in terms of its value at 4 redshift nodes, z = {0.24, 0.53, 0.83, 1.5}, from which we interpolate and extrapolate to other redshifts using a quadratic spline (see Section 2.2 for further details).
The top panel of Figure 8 shows the amplitude of matter fluctuations as a function of redshift, parametrized by S 8 (z) defined as where σ 8 (z) = D(z) σ fid 8 , σ fid 8 = 0.8111 and Ω m is the present value of the matter density parameter. This allows us to visualize the evolution of structure growth as directly measured by the data while setting the amplitude to the parameter that these data sets are able to constrain best. Results are shown for the ND (purple) and SD (orange) data sets, as well as the full FD combination (red), together with the growth history predicted by Planck in the ΛCDM model (blue). The bottom panel of the same figure shows the same reconstructed growth for different probe combinations, excluding shear data (gray) and galaxy clustering (green). As the figure shows, most of the constraining power from these projected data sets comes from the redshift range 0.25 z 0.7, through the combination of photometric clustering and cosmic shear, and, to a lesser extent, at z ∼ 2 from the combination of QSO : Reconstructed growth history from the SD data using our fiducial quadratic spline (orange) and a linear spline (gray) showing 68% C.L. errors. Both methods recover compatible growth histories with slight variations that become most prominent at low redshifts beyond the interpolation range.
clustering and CMB lensing. By comparison with the Planck ΛCDM prediction, we see that the tension between current large-scale structure and CMB data originates around z ∼ 0.4, where the data display a significantly lower amplitude of fluctuations. Furthermore, this feature is recovered consistently by both the ND and SD samples independently. Furthermore, this tension seems to be driven by the cosmic shear data, particularly KiDS, whereas the combination of galaxy clustering and CMB lensing seems to be in reasonable agreement with Planck (although see [7]). The current uncertainties, however, are too large to show any tension between these two probes.
The complementarity between weak lensing and galaxy clustering data allows us to reconstruct the growth history from z ∼ 2. Nevertheless, more data at low (z < 0.2) and high (z > 1) redshifts would be necessary in order to fully recover the evolution of matter fluctuations, and better understand the origin of the slower growth potentially preferred by current data at intermediate redshifts.
When interpreting the reconstructed growth history it is important to bear in mind the potential impact of the choice of interpolation scheme (quadratic spline) on the results. This could have an effect both on the recovered growth parameters and on the extrapolation of the growth factor in redshift ranges where we do not have direct measurements. To quantify this, we have repeated our analysis of the SD data set using linear interpolation between nodes. The result, shown in Figure 9, demonstrates that the choice of interpolation scheme does not have a significant effect on the recovered growth constraints in the range where the galaxy clustering and weak lensing data lie. However, the amplitude of perturbations predicted at z ∼ 0, as well as its statistical uncertainties can be significantly affected by this choice, and therefore care should be taken when interpreting these constraints beyond the range of redshifts covered by the data. In addition, we have checked that we recover the same χ 2 for the best fit ΛCDM SD, ND and FD cases if we compute it with the reconstructed growth model using theD(z i ) values found in ΛCDM as the spline nodes. More precisely, the relative deviations of the χ 2 between both cases for the SD, ND and FD datasets are -0.06%, -0.02%, -0.005%, respectively. Finally, we also checked that we are able to recover the fiducial cosmology when using simulated data in our MCMC.
It is interesting to understand in better detail the redshift ranges over which the different data sets studied here constrain structure growth. To this end, the left and right panels of Figure 10 show the constraints on S 8 (z) from different combinations of the ND and SD experiments respectively. As the figure shows, the constraints at low redshifts (z < 0.5) are dominated by the cosmic shear data sets, while the intermediate redshift regime 0.5 z 2 is dominated by the combination of galaxy clustering and CMB lensing. Both data sets are thus highly complementary and their combination is able to significantly improve the accuracy of the reconstruction over the full redshift range. This is understandable given the cumulative nature of weak lensing, which weights it more heavily towards low redshifts. At high redshifts, the overlap between the galaxy clustering and shear kernels is reduced, whereas the overlap with the CMB lensing kernel is equally significant at all redshifts. The figure also shows that the slight kink at z 0.8 visible in Figure 8 originates from the KiDS cosmic shear data, likely due to a statistical fluctuation.
To further explore the complementarity between different probes, Figure 11 shows the two-dimensional confidence intervals for the growth reconstruction parametersD z in the case of the SD experiments for different data set combinations. Besides the complementarity between shear and clustering data sets in terms of the final statistical uncertainties onD z we just described, we can also interpret the level of correlation between different parameters. While the constraints from the combination of clustering and CMB lensing are approximately uncorrelated between different growth parameters, the constraints from cosmic shear show a clear positive correlation between adjacentD z s. This is also easy to understand in terms of the radial kernel associated to the different tracers. While galaxy clustering traces the amplitude of inhomogeneities locally (i.e. in the redshift range where the galaxies are selected), weak lensing traces the cumulative distribution, and thus the measurements of the shear power spectrum in different redshift bins are visibly correlated.
The final quantitative constraints onD z for the different data set combinations are listed in Table 7, together with the expected values of the growth factor estimated from the FD no FD no g FD Figure 11: Posterior contours at 68 and 95% C.L. for the reconstructed growth nodes. Results are shown for different data combinations of the FD data set. The full constraints are shown in orange, while the green and red contours show the constraints in the absence of galaxy clustering or cosmic shear respectively. Due to its cumulative nature, constraints from cosmic shear produce a visible correlation between adjacent redshift nodes. The constraints on the highest-redshift node are completely dominated by the combination of eBOSS-QSO and CMB lensing.
Planck 18 ΛCDM chains [1]. As is evident from the previous figures, the evidence of tension with respect to CMB data is concentrated on the first two redshift nodes. We find only a mild improvement in the goodness of fit between ΛCDM and this more flexible model. For instance, when using the FD data we find ∆χ 2 = −7, which does not represent strong evidence against ΛCDM given the 3 additional parameters of this more flexible model. We therefore conclude that the linear growth history preferred by current large-scale structure data is in good agreement with the expectation from ΛCDM 10 , and that the main source of  Table 7: Constraints on the growth reconstruction nodesD z for different combinations of data sets. Note that a fifth nodeD 5 = 0.212 was used as anchor to match Planck at z ≥ 5 when doing the interpolation. The Planck 18 ΛCDM row was obtained reanalyzing the official chains [1].
tension with CMB data is in the overall amplitude of fluctuations. In the future, it will be interesting to incorporate RSD measurements in this formalism to simultaneously constrain the growth factor and its derivative with respect to redshift.

Discussion
In this paper we set out to answer the question: "What do current large-scale structure data say about the growth of density fluctuations in the late Universe?". To do so, we have carried out a combined analysis of 6 different projected probes of the large-scale structure, separating the linear growth from the background expansion and parametrizing the former in terms of its value at a set of four redshift nodes covering the range 0.2 z 2. The main motivation for this exercise is to shed light on the tension in the value of the current amplitude of linear fluctuations as measured by low-redshift weak lensing data [4][5][6] and as extrapolated from CMB data in the context of ΛCDM. Answering this question has allowed us to a) evaluate the consistency between different low-redshift probes in their predicted growth history and b) determine at what point in time we have evidence of a slower growth than that predicted by Planck. Our work is based on an independent harmonic-space-based analysis, making use of state-of-the-art methods to obtain unbiased and precise measurements of power spectra for different probes and the associated covariance matrix. This has allowed us to carry out a consistent analysis of galaxy clustering from DES, DELS and eBOSS-QSO, cosmic shear from KiDS and DES, and CMB lensing from Planck. The methods used have been extensively validated and, as we have shown, we find no significant evidence of unaccounted systematic uncertainties in our data vector and covariances within the range of scales we use.
Our main result, shown in Figure 8, is the reconstructed growth history at late times. Our analysis has shown that the tension in the measured value of S 8 is due to the data surements ofDz alone.
preferring a lower amplitude of fluctuations in the range 0.25 z 0.7 than the Planck prediction with a significance higher than 2σ. The data used here are not able to place strong constraints at lower or higher redshifts and thus show no deviation with respect to the Planck preferred value in this regime. As we have shown, this tension is driven by the cosmic shear data sets, whereas the combination of galaxy clustering and CMB lensing used here shows no significant deviation with respect to Planck. Our uncertainties are not small enough to reveal any tension between these probes, however, although the results found by other groups with different data [7] suggest that the tension is also revealed by CMB lensing tomography.
In order to validate our analysis pipeline, we have verified that the constraints on ΛCDM from our data agree well with those found by the different collaborations. Doing so has also allowed us to explore the constraints our full complement of data sets is able to place on this model when combined, measuring the S 8 parameter to be S 8 = 0.7781 ± 0.0094. This is in good agreement with current constraints from other groups, and in tension with Planck at the ∼ 3.4σ level. This also shows that, in combination, current large-scale structure data are able to constrain these parameters with an uncertainty that is significantly (∼ 25%) better than that achieved by primary CMB data.
The results presented here are subject to some caveats. First, the specific properties of the method used to parametrize the linear growth as a free function of redshift (a quadratic spline), has an impact on the results in terms of e.g. the correlation length of this function, or its behaviour outside of the range of redshifts constrained by the data. A more ideal choice would have been, for example, to make use of Gaussian process regression, in order to let the data constrain the smoothness of this function, and to robustly quantify its uncertainty in poorly-constrained regimes. This approach would significantly increase the complexity of the model and most likely require the use of differentiable sampling methods, and therefore we leave this for future work.
Since the focus of this paper is not to derive constraints on ΛCDM, our analysis of this model has overlooked a number of issues that should be explored in more detail. Although the ND and SD data sets are spatially disjoint, the systematic uncertainties in their redshift distribution are, to some extent, correlated due to the use of similar spectroscopic calibration samples. Furthermore, as shown in [88], the potential biases in these calibrating samples could give rise to systematic effects in the derived redshift distributions that are not sufficiently well captured by the mean shifts. This is in addition to other unknown observational systematics correlated across overlapping probes. We have ignored the impact of baryons on our final constraints. Although this effect was found to be subdominant within the scales used here by individual collaborations, the increased sensitivity of their combination could make it more relevant. The same can be said about the level of complexity used to parametrize the impact of intrinsic alignments. Finally, our analysis has made use of the simplest linear bias model to describe galaxy clustering, which significantly limits the range of scales over which this tracer can be used. As shown in [94] the use of recent advances in the modelling of galaxy bias (e.g. [98,99]) would allow us to significantly increase this range, which would have a direct impact on the final cosmological constraints. This could be particularly relevant in the case of growth reconstruction, since the galaxy clustering samples act as redshift anchors enabling a more precise recovery of the redshift dependence.
At this point it is worth asking ourselves what possible explanations could resolve this tension in the light of the growth reconstruction results presented here. The three possible scenarios are the presence of unknown systematic uncertainties in the CMB analysis, the presence of systematic uncertainties in the analysis of cosmic shear (which, as we have shown, is the driving source of this tension), or the presence of new physics giving rise to a slower growth of perturbations at late times. The results found by the Atacama Cosmology Telescope (ACTPol) collaboration independently of Planck data [100], obtain S 8 = 0.840 ± 0.030 in combination with WMAP data, in perfect agreement with Planck and in 2σ tension with our results. Therefore there is currently no significant evidence of systematic uncertainties which have been unaccounted for in the Planck measurement. As we have shown, both DES and KiDS independently recover the same lower amplitude of fluctuations at z ∼ 0.4, and therefore the possibility of experiment-specific systematic uncertainties being the cause of the tension is not strongly favoured. The other possibility would be a mis-modelling of unknown astrophysical systematic uncertainties specific to cosmic shear and therefore common to different experiments. However, similar results have been found by others without cosmic shear data [7].
We must therefore turn to the third scenario: new physics. A number of extensions to ΛCDM have been explored in [101,102] in the context of DES and KiDS respectively. While some simple extensions, such as allowing for departures of the dark energy equation of state from a perfect cosmological constant, are able to reconcile the measured values of S 8 , significant tension still remains in the regions of the multi-dimensional parameter space preferred by the CMB and large-scale structure data sets [102]. Furthermore, the data are not yet able to detect departures from a flat Universe, as well as any signatures of modified gravity [103]. A less exotic proposal (see e.g. [7,19,94]) is the possibility of a value of Ω m slightly lower than that found by CMB experiments, which would also agree with other probes of structure growth.
A final scenario, although arguably a less interesting one, is that the current tension is simply a statistical fluke which sits at the ∼ 3σ level for this particular combination of parameters due to a slight mis-estimation of subdominant systematic uncertainties. Fortunately, current constraints are far from sample-variance dominated. Most immediately, the imminent release of the DES third-year analysis is expected to shed light on this issue, either reinforcing or alleviating the existing tension. Further down the line, data from nextgeneration photometric surveys such as the Vera C. Rubin Observatory [104] or Euclid [105], spectroscopic surveys like DESI [106], and CMB experiments such as the Simons Observatory [107] and CMB Stage-4 [108] will allow us to resolve and exploit the source of this tension. These data will also allow us to map the growth history in a model-independent way with a much higher precision than current data allow over a larger range of redshifts. The question of whether these forthcoming data will be able to strengthen or alleviate the S 8 tension remains open.
Significant progress can still be made on the constraints presented here with currentlyavailable data. At low redshifts, where the data used here are not able to strongly constrain growth, our measurements could be improved by including information from all-sky photometric surveys [109,110], and, potentially, from maps of the thermal Sunyaev-Zel'dovich effect as a probe of structure [62,[111][112][113] assuming astrophysical uncertainties can be kept under control [114,115]. At redshifts z 1 further information could be gained by combining these data with infrared data such as unWISE [7,116] or radio continuum data [117,118]. Finally, other high-precision probes of growth, such as redshift-space distortions and peculiar velocity surveys [119,120], should allow us to constrain the growth history at least at the same level as tomographic data, and therefore it is worth developing methods to combine these measurements in a consistent way. The Legacy Survey team makes use of data products from the Near-Earth Object Widefield Infrared Survey Explorer (NEOWISE), which is a project of the Jet Propulsion Laboratory/California Institute of Technology. NEOWISE is funded by the National Aeronautics and Space Administration.
The   Figure 12: Posterior distribution for the growth reconstruction nodes for the cases using only the Gaussian part of the covariance matrix and the additional non-Gaussian terms. We use only DES data as a worst case scenario where the last node remains unconstrained. On the top right corner, we show the constraints for the (Ω m , S 8 ) parameters for ΛCDM using DES + CMBκ. As can be seen, the non-Gaussian part does not affect the posterior distributions. The small differences are likely to be originated by different chain lengths.

A Effect of non-Gaussian covariances
In order to verify that our results are insensitive to the non-Gaussian contributions to the covariance matrix over the range of scales used here, we have repeated the analysis for both ΛCDM and the growth reconstruction in two different ways that cover different situations. In the first case, we use the DES 5 × 2pt data including the connected trispectrum and supersample terms computed as described in [46]. In comparison, for the reconstruction of growth we only use the DES 3 × 2pt data as a worse case scenario, where there is no cosmological information to constrain the last nodeD 1.5 .. The results are shown in Figure 12 for the growth rate parametersD z and the ΛCDM parameters (Ω m , S 8 ). The non-Gaussian contributions have a negligible effect on the final best-fit parameters and their uncertainties within the range of scales used here. The impact on the goodness of fit is also small, with the best-fit χ 2 changing by less than 1% in both cases. Given these results and for simplicity, we use only the dominant Gaussian part of the covariance matrix, calculated as described in 4.2.