11institutetext: Affiliations are listed at the end of the paper.

#  KiDS & Euclid: Cosmological implications of a pseudo angular power spectrum
analysis of KiDS-1000 cosmic shear tomography††thanks: This paper is published
on behalf of the KiDS Collaboration and the Euclid Consortium.

A. Loureiro1,2,3<EMAIL_ADDRESS>L. Whittaker1,4 A. Spurio
Mancini1,5 B. Joachimi1 A. Cuceu1 M. Asgari6,3 B. Stölzner1 T. Tröster3 A.H.
Wright7 M. Bilicki8 A. Dvornik9 B. Giblin3 C. Heymans3,7 H. Hildebrandt7 H.
Shan10,11 A. Amara12 N. Auricchio13 C. Bodendorf14 D. Bonino15 E.
Branchini16,17 M. Brescia18 V. Capobianco15 C. Carbone19 J. Carretero20,21 M.
Castellano22 S. Cavuoti23,18,24 A. Cimatti25,26 R. Cledassou27,28 G. Congedo3
L. Conversi29,30 Y. Copin31 L. Corcione15 M. Cropper5 A. Da Silva32,33 M.
Douspis34 F. Dubath35 C.A.J. Duncan36 X. Dupac30 S. Dusini37 S. Farrens38 S.
Ferriol31 P. Fosalba39,40 M. Frailis41 E. Franceschi13 M. Fumana19 B.
Garilli19 B. Gillis3 C. Giocoli42,43 A. Grazian44 F. Grupp14,45 S.V.H.
Haugan46 W. Holmes47 F. Hormuth48,49 K. Jahnke49 M. Kümmel45 S. Kermiche50 A.
Kiessling47 M. Kilbinger38 T. Kitching5 K. Kuijken51 M. Kunz52 H. Kurki-
Suonio53 S. Ligori15 P.B. Lilje46 I. Lloro54 O. Mansutti41 O. Marggraf55 K.
Markovic47 F. Marulli13,56,57 R. Massey58,59 M. Meneghetti13,56 G. Meylan60 M.
Moresco57,13 B. Morin38 L. Moscardini57,13,56 E. Munari41 S.M. Niemi61 C.
Padilla20 S. Paltani35 F. Pasian41 K. Pedersen62 V. Pettorino38 S. Pires38 M.
Poncet27 L. Popa63 F. Raison14 J. Rhodes47 H. Rix49 M. Roncarelli57,13 R.
Saglia64,45 P. Schneider55 A. Secroun50 S. Serrano39,40 C. Sirignano65,37 G.
Sirri56 L. Stanco37 J.L. Starck38 P. Tallada-Crespí66,21 A.N. Taylor3 I.
Tereno32,67 R. Toledo-Moreo68 F. Torradeflot21,66 E.A. Valentijn69 Y. Wang70
N. Welikala3 J. Weller14,45 G. Zamorani13 J. Zoubian50 S. Andreon71 M.
Baldi57,13,56 S. Camera72,15,73 R. Farinelli74 G. Polenta75 N. Tessore1

(Received October 19, 2021; accepted June 27, 2022)

We present a tomographic weak lensing analysis of the Kilo Degree Survey Data
Release 4 (KiDS-1000), using a new pseudo angular power spectrum estimator
(pseudo-$C_{\ell}$) under development for the ESA Euclid mission. Over 21
million galaxies with shape information are divided into five tomographic
redshift bins, ranging from 0.1 to 1.2 in photometric redshift. We measured
pseudo-$C_{\ell}$ using eight bands in the multipole range $76<\ell<1500$ for
auto- and cross-power spectra between the tomographic bins. A series of tests
were carried out to check for systematic contamination from a variety of
observational sources including stellar number density, variations in survey
depth, and point spread function properties. While some marginal correlations
with these systematic tracers were observed, there is no evidence of bias in
the cosmological inference. B-mode power spectra are consistent with zero
signal, with no significant residual contamination from E/B-mode leakage. We
performed a Bayesian analysis of the pseudo-$C_{\ell}$ estimates by forward
modelling the effects of the mask. Assuming a spatially flat $\Lambda$CDM
cosmology, we constrained the structure growth parameter
$S_{8}=\sigma_{8}(\Omega_{\rm m}/0.3)^{1/2}=0.754_{-0.029}^{+0.027}$. When
combining cosmic shear from KiDS-1000 with baryon acoustic oscillation and
redshift space distortion data from recent Sloan Digital Sky Survey (SDSS)
measurements of luminous red galaxies, as well as the Lyman-$\alpha$ forest
and its cross-correlation with quasars, we tightened these constraints to
$S_{8}=0.771^{+0.006}_{-0.032}$. These results are in very good agreement with
previous KiDS-1000 and SDSS analyses and confirm a $\sim 3\sigma$ tension with
early-Universe constraints from cosmic microwave background experiments.

###### Key Words.:

gravitational lensing: weak, methods: data analysis, methods: statistical,
surveys, cosmology: observations

##  1 Introduction

Since the turn of the millennium, there has been a dramatic increase in the
number and precision of new cosmological observations. Measurements of the
cosmic microwave background (CMB) radiation have achieved extraordinary
precision for the parameters of the standard cosmological model — the
$\Lambda$CDM model, which contains mostly the cosmological constant
($\Lambda$) and cold dark matter (CDM). The most stringent constraints were
obtained by the Planck Mission (Planck Collaboration et al., 2020a), which
demonstrated an unprecedented constraining power, measuring some of these
parameters with percent-level precision. The Planck Mission, however, did not
contradict the standard model of cosmology. Instead, it established the
$\Lambda$CDM model as a difficult paradigm to overturn for other cosmological
probes and future surveys (Efstathiou & Gratton, 2019, 2020).

Measurements of the CMB provide insight into the characteristics of the early
Universe, with information coming from the surface of last scattering, at
redshift $z\approx 1100$. If the $\Lambda$CDM paradigm is correct, one expects
agreement between the cosmological parameters estimated using CMB experiments
and those recovered from late-time probes. In other words, late-time
cosmological parameters such as the amplitude of the matter density
fluctuations, $\sigma_{8}$, predicted by CMB observations, should match what
we observe today with galaxy surveys. Several studies seem to suggest tensions
between late-time probes and the CMB (Mörtsell & Dhawan, 2018; Verde et al.,
2019). These tensions arise in a few parameters, such as the Hubble constant,
when constrained using Cepheids and type Ia Supernovae using distance ladder
measurements (Riess et al., 2011, 2016; Bernal et al., 2016; Di Valentino et
al., 2016; Lin & Ishak, 2017; Efstathiou, 2020; Riess et al., 2020), and the
structure growth parameter, $S_{8}=\sigma_{8}\sqrt{\Omega_{\rm m}/0.3}$, when
measured using cosmic shear surveys (Heymans et al., 2013; MacCrann et al.,
2015; Joudaki et al., 2017; Hildebrandt et al., 2017; Abbott et al., 2018;
Park & Rozo, 2020; Hikage et al., 2019; Lemos et al., 2021; Asgari et al.,
2020; Joudaki et al., 2020; Tröster et al., 2021; Asgari et al., 2021; Heymans
et al., 2021; Amon et al., 2021; Secco et al., 2022; DES Collaboration et al.,
2022).

Cosmic shear is the study of the weak gravitational lensing of background
galaxies due to the intervening large-scale structure of the Universe. This
lensing effect produces correlations in galaxy shapes, and these correlations
can be used to determine how the dark and baryonic foreground matter is
distributed (see Bartelmann & Schneider 2001; Kilbinger 2015 for comprehensive
reviews of cosmic shear). Moreover, by measuring the cosmic shear signal in
tomographic redshift bins, we can study how the distribution of matter evolves
over time; this enables us to place tight constraints on the structure
formation of the Universe and its evolution with redshift (Hall, 2021), as
well as the dark energy equation-of-state and other standard cosmological
parameters.

As more data are independently collected by current galaxy surveys, the slight
tension in $S_{8}$ between CMB and cosmic shear surveys does not seem to
disappear. More specifically, all results coming from the Dark Energy Survey
(DES, Abbott et al. 2018; Troxel et al. 2018) and the Hyper Supreme Camera
(HSC, Hikage et al. 2019) indicate a lower value for the amplitude of matter
density fluctuations, consistent with the Kilo Degree Survey (KiDS,
Hildebrandt et al. 2020; Heymans et al. 2021), but also consistent with the
Planck projections at a $\gtrsim 1.6\sigma$ level. Combining a subset of these
independent surveys, which all use different instruments and observing
strategies, the tension with the early Universe probes is exacerbated. In a
re-analysis of DES first year data combined with KiDS, Asgari et al. (2020)
and Joudaki et al. (2020) demonstrate that the tension can be around
3.2$\sigma$ between cosmic shear and CMB measurements from Planck. These
consistent results from different experiments encourage us to believe that it
is unlikely that this tension arises from an observational systematic error;
however, there is still room for different explanations (Joachimi et al.,
2021a). From a similar perspective, recent results from the Atacama Cosmology
Telescope (ACT), a ground-based CMB experiment, observe a similar tension in
$S_{8}$ with cosmic shear probes when combined with other CMB experiments
(Aiola et al., 2020; Han et al., 2021). This also suggests that it is unlikely
that the tension we currently observe comes from a systematic contamination in
the CMB experiments. If this discrepancy between measurements is indeed
induced by systematic contamination, it would have to be shared across several
independent datasets with widely different properties and analysis choices.
Hopes are that experiments planned for the near future will shine a light on
this interesting cosmic puzzle as there is still room to improve the precision
on cosmic shear experiments.

Future galaxy surveys such as the Legacy Survey of Space and Time (LSST),
carried out at the Vera C. Rubin Observatory (LSST Dark Energy Science
Collaboration, 2012), the European Space Agency’s Euclid Mission (Laureijs et
al., 2011), and the National Aeronautics and Space Administration’s (NASA)
Nancy Grace Roman Space Telescope (Spergel et al., 2015) will strongly rely on
cosmic shear as one of their primary probes to understand the large-scale
structure of the Universe. These future surveys will obtain cosmological
information from billions of galaxy shapes using a plethora of statistical
tools including correlation functions, angular power spectra, mass maps, and
bi-spectra. Amongst these, the pseudo-$C_{\ell}$ technique is a direct
harmonic space approach to estimate the angular power spectra of cosmic shear
and galaxy positions (Peebles, 1973; Brown et al., 2005; Hikage et al., 2011;
Asgari et al., 2018; Alonso et al., 2019; Nicola et al., 2021; García-García
et al., 2021).

The pseudo-$C_{\ell}$ technique has its strengths and weaknesses in comparison
with other two-point (2pt) function estimators (Efstathiou, 2004; Asgari et
al., 2018; Nicola et al., 2021). Among its shared advantages with other 2pt
functions, pseudo-$C_{\ell}$s have a fast and simple application to data
distributed on the sphere, with no need for a flat-sky approximation. The
speed to obtain pseudo-$C_{\ell}$ estimates, compared to 2pt correlation
functions, is particularly convenient when estimating covariances from
simulated catalogues. Furthermore, measuring the 2pt function in harmonic
space increases the speed of theory calculations for cosmological inference as
most Einstein-Boltzmann equation solvers use the advantages of calculating
correlations in harmonic space. Performing the analysis in harmonic space
avoids transforming the theory-vector to real space, as it is done for
correlation functions (Schneider et al., 2002), or transforming the data-
vector into harmonic space, as is the case for band-power analysis (Becker &
Rozo, 2016; van Uitert et al., 2018).

Yet, pseudo-$C_{\ell}$ techniques have one liability regarding partial sky
observations. Naturally, any realistic galaxy survey will contain masked parts
of the sky due to bright foreground objects, survey strategy and other
observational artefacts.111With the increase of extremely bright low-orbit
satellite tracks from mega satellite constellations, this problem will affect
the masks geometry of future ground-based surveys even more, with significant
consequences for cosmology and other areas of astronomy (Gallozzi et al.,
2020a, b). In the presence of a mask, the mask’s characteristic scales, from
its shapes and holes, will mix nearby modes in harmonic space causing
correlations between multipoles to appear. For spin-2 fields, such as the
shear field, partial sky observations can also cause ambiguity between
E/B-modes. The multipole mixing, however, can be fully accounted for by
analytically calculating the mixing matrix from the survey’s geometry.

In this work, our approach consists of forward-modelling the mode mixing
effect due to masking in the theory. Taking this approach, we are least prone
to numerical instability from inverting and deconvolving the mixing matrix
from pseudo-$C_{\ell}$ measurements. Our pseudo-$C_{\ell}$ method is a
prototype for one of the default 2pt functions used in Euclid (Euclid Collab.,
in prep.). The motivation for the forward-modelling of the mask effects is so
that the estimator can achieve the requirements for a pseudo-$C_{\ell}$
analysis in Euclid (Laureijs et al., 2011; Euclid Collaboration et al., 2020;
Tutusaus et al., 2020). Although binning the data-vector and mixing matrix to
perform a more stable deconvolution (as done in Hikage et al. (2019) and
Tröster et al. (2022)) poses no challenges for the data set used in this
analysis, for the future Euclid survey this is challenging. For Euclid, we are
required to measure angular power spectra in 100 log-spaced band-powers for 10
to 13 redshift tomographic bins (Euclid Collaboration et al., 2020; Tutusaus
et al., 2020). This would require dealing with a mixing matrix with 2.1 to
3.51 Million dimensions, posing potential numerical inaccuracies which might
well exceed the very tight accuracy requirements for the Figure-of-Merit for
the equation-of-state of Dark Energy measured by the Euclid Mission (Laureijs
et al., 2011).

Hence, we apply a pseudo-$C_{\ell}$ estimator currently being developed for
the Euclid Science Ground Segment to state-of-the-art data from a Stage-III
cosmic shear experiment: the public ESO Kilo Degree Survey (KiDS). KiDS is
currently at its fourth Data Release (KiDS DR4, Kuijken et al. 2019), having
observed around 1000 deg2 of the sky with nine-band matched photometry. Giblin
et al. (2021) outlined the KiDS-1000 weak lensing catalogue production and
validation, with redshift calibration and measurements presented in
Hildebrandt et al. (2021). Using three different 2pt statistics, Asgari et al.
(2021) performed a cosmic shear analysis of the KiDS-1000 catalogue, including
a series of consistency checks. The main KiDS-1000 analysis is finally
presented in Heymans et al. (2021), combining it with galaxy clustering from
BOSS DR12 (Reid et al., 2016; Sánchez et al., 2017a), to perform a joint
($3\times 2$pt) analysis for a flat $\Lambda$CDM model, following the
methodology described by Joachimi et al. (2021b). Finally, constraints on
extensions to $\Lambda$CDM using the same $3\times 2$pt methodology and data
are presented in Tröster et al. (2021).

Here, we demonstrate that a pseudo-$C_{\ell}$ estimator performs comparatively
accurately to other 2pt estimators, and we show the advantages of forward-
modelling effects of the mask through the use of the mixing matrix as part of
the theory modelling in Bayesian inference. We also demonstrate the
versatility of using this estimator for systematic contamination analysis,
calculating the angular power spectra between our data and some observational
artefacts and characteristics. Finally, we show the constraining power of
combining cosmic shear data with the latest clustering data from the Baryon
Oscillations Spectroscopic Survey and its extension (BOSS and eBOSS,
respectively; Alam et al., 2017; Gil-Marín et al., 2020; Bautista et al.,
2021; du Mas des Bourboux et al., 2020; Alam et al., 2021). We demonstrate
that our independent analysis of the KiDS-1000 cosmic shear catalogue exhibits
a similar level of tension with early Universe measurements from Planck Legacy
(Planck Collaboration et al., 2020a) and the Atacama Cosmology Telescope Data
Release 4 (Aiola et al., 2020).

This study is organised as follows: Section 2 summarises the KiDS-1000 shear
catalogue and shape measurements. Section 3 explains the methodology for
pseudo-$C_{\ell}$ estimators applied to cosmic shear data, the mixing matrix
calculation, the measurements performed on the data for E/B-modes and the
covariance matrix estimation. Section 4 explains our systematic contamination
null tests using the pseudo-$C_{\ell}$ estimator. Section 5 details our
Bayesian inference: theory modelling for cosmic shear, forward-modelling of
the mixing matrix, external galaxy clustering data from baryon acoustic
oscillations (BAO) and redshift space distortions (RSD) used to complement the
cosmic shear data, as well as the priors and the likelihood implemented in the
analysis. Section 6 outlines our main results from the inferred posteriors
with a discussion on tension with cosmic microwave background measurements.
Section 7 summarises our findings and conclusions. Extra figures with
cosmological constraints and nuisance parameters, as well as a complete table
with all parameters probed in our analysis can be found in Appendix A.
Finally, Appendix B discusses the details of the systematics null tests for
B-modes specifically.

##  2 Cosmic shear data from KiDS-1000

Designed as a weak lensing experiment, the Kilo Degree Survey (Kuijken et al.,
2015; de Jong et al., 2015, 2017) is a public survey by the European Southern
Observatory, currently at its 4th Data Release (Kuijken et al.,
2019)222http://kids.strw.leidenuniv.nl/DR4/, with an observed area of around
$1000$ deg2. Good seeing nights were prioritised to obtain high-resolution
images in the r-band, resulting in an excellent average seeing of 0.7 arc-
seconds. Like in previous analyses (Wright et al., 2019; Hildebrandt et al.,
2020), KiDS DR4 photometry is processed with Theli (Erben et al., 2005;
Schirmer, 2013) and AstroWise (Begeman et al., 2013) using five-band NIR
photometry ($ZYJHK_{s}$) from VIKING333The VISTA Kilo-degree Infrared Galaxy
(VIKING) Survey (Edge et al., 2013) combined with four-band optical photometry
($ugri$) from KiDS to estimate photometric redshifts in an accurate and
precise way. In this section, we briefly describe the most important and
relevant details of the catalogue used in our study.

We perform our main analysis using the latest cosmic shear catalogue from the
4th Data Release (KiDS-1000)444http://kids.strw.leidenuniv.nl/DR4/lensing.php,
covering an effective unmasked area of 777.4 deg2 and containing a total of 21
262 011 lensed galaxies. Detailed descriptions of the catalogue’s
construction, including the LensFit (Miller et al., 2007; Kitching et al.,
2008; Miller et al., 2013) shape measurement procedures, inverse variance
weight estimation, with contributions from measurement error and intrinsic
shape noise, as well as systematic contamination tests, can be found in Giblin
et al. (2021). In Giblin et al. (2021), several potential contaminants to the
data were analysed, including, but not limited to, point spread function (PSF)
contamination and the effects of multiplicative and additive bias calibration,
finding an impact smaller than $0.1\sigma$ for
$S_{8}=\sigma_{8}\sqrt{\Omega_{\rm m}/0.3}$ after the calibration is applied.
We perform additional checks in the pseudo-$C_{\ell}$ approach in Sect. 4.

Here, for ease of comparison, we apply the same redshift tomographic binning
as in Asgari et al. (2021), adopting five bins selected using the most
probable redshift assigned by BPZ (Benítez, 2000; Coe et al., 2006; Raichoor
et al., 2014), $z_{\rm B}$. The first four bins have $\Delta z_{\rm
B}^{(1\text{-}4)}=0.2$, ranging from $0.1<z_{\rm B}\leq 0.9$, while the last
tomographic bin ranges from $0.9<z_{\rm B}<1.2$, that is $\Delta z_{\rm
B}^{(5)}=0.3$. The photometric redshift distributions for the lensing sample,
detailed in Hildebrandt et al. (2021), were estimated using the Self-
organising Maps technique (SOM; Maehoenen & Hakala 1995; Naim et al. 1997;
Masters et al. 2015; Kitching et al. 2019b; Wright et al. 2020) to match
galaxies using the nine-band photometry to groups with similar properties
within spectroscopic samples. We select galaxies detected in all nine bands
which have SOM matches and follow extra criteria according to the Gold Sample
presented in Hildebrandt et al. (2021). Figure 1 shows the tomographic
selection and the SOM redshift distributions which have been validated with a
cross-clustering analysis using spectroscopic surveys (see van den Busch et
al. 2020 and Hildebrandt et al. 2021 for details). Table 4 contains more
details about the samples’ redshift tomographic bins and their multiplicative
shear calibration corrections.

Cosmic shear, ${\gamma}$, can be estimated using the galaxy ellipticity
measurements, corrected by the additive bias, $\epsilon^{\rm corr}$, obtained
from the observed galaxy ellipticities, $\epsilon^{\text{obs}}$. In the weak
lensing regime ($|{\gamma}|\ll 1$), $\epsilon^{\rm corr}$ is estimated by
taking into account the additive bias $c$,

$\displaystyle\epsilon^{\rm corr}={\epsilon^{\text{obs}}-{c}}\,,$ (1)

where the additive bias, ${c}=\langle\epsilon^{\text{obs}}\rangle$, is assumed
to be the weighted average observed ellipticity over all galaxies in a given
tomographic bin. We note that all these quantities (shear, ellipticity and
additive bias) are complex quantities, where
$\epsilon=\epsilon_{1}+\rm{i}\epsilon_{2}$.

Figure 1: KiDS-1000 redshift distribution for the five tomographic bins used.
Solid lines show the redshift distributions obtained via the SOM technique
(Hildebrandt et al., 2021). The shaded bands show the corresponding ranges of
photometric redshift point estimates, $z_{\rm B}$.

##  3 Methodology

This section briefly discusses the pseudo-$C_{\ell}$ methodology used to
estimate the cosmic shear angular power spectra. The prototype estimator
applied here is currently under development for Euclid (Laureijs et al.,
2011). A detailed study of the methodology performance for Stage IV surveys
will be presented in future work. Apart from the noise power spectra
estimation, our method is similar to previous works for CMB polarisation power
spectra (Kogut et al., 2003; Brown et al., 2005) and the full-sky formalism
presented by Hikage et al. (2011).

In contrast to the approach taken by Hikage et al. (2019) and Tröster et al.
(2022), we focus on forward-modelling the effects of the mask via the mixing
matrix instead of deconvolving it from the data as done in widely used
pseudo-$C_{\ell}$ implementations such as Hivon et al. (2002) and Alonso et
al. (2019). It is important to note that although not obvious at first, the
forward modelling of the mask does not imply any extra numerical overhead when
modelling the theory vector (discussed in detail in Sect. 5.1). The binning
and mixing matrix deconvolution operations do not commute. Therefore, if one
chooses to take the standard approach of deconvolving the binned data-vector,
similar operations need to be performed for the theoretical modelling during
the likelihood analysis (convolve the theory with the mixing matrix, bin the
theory-vector, bin the mixing matrix, and finally deconvolve it from the
theory vector). All these operations can be written as a single matrix
multiplication (Alonso et al., 2019), which is precisely what the forward
modelling of the mixing matrix requires. Hence, both approaches have similar
computational time with the sole difference that for future galaxy surveys the
binned mixing matrix will be much larger and harder to invert numerically. The
forward model approach we take in this work avoids this potential issue.

Details on the mixing matrix estimation are given in Sect. 3.2, while
modelling of the data-vector is detailed in Sect. 5.1. The methodology we use
for covariance matrix estimation using simulations is detailed in Sect. 3.4.

###  3.1 Pseudo angular power spectrum estimator for cosmic shear

We represent the cosmic shear fields by pixelated maps created from shear
estimates derived from the catalogue using a weighted average and correcting
for the multiplicative bias, $m_{z}$, where z is the redshift tomographic bin
index (Kitching et al., 2019a, 2020). The average shear in a pixel localised
by the unit vector $\hat{\mathbf{n}}$ is given as

${\gamma}\left(\hat{\mathbf{n}}\right)=\frac{1}{(1+m_{z})}\frac{\sum_{i\in\hat{\mathbf{n}}}w_{i}\epsilon^{\rm
corr}_{i}}{\sum_{i\in\hat{\mathbf{n}}}w_{i}}\,,$ (2)

where the summation is taken over all galaxies in the pixel. For each galaxy,
$\epsilon^{\rm corr}_{i}$ is constructed using Eq. (1), with weights, $w_{i}$,
provided by shape measurement techniques such as LensFit (Miller et al., 2007;
Kitching et al., 2008; Miller et al., 2013). While the additive bias from Eq.
(1) can be estimated from the catalogue itself, the multiplicative bias
($m_{z}$, m-correction) needs to be estimated from careful calibration using
pixel-level simulations and is chosen in KiDS to be the average value for all
galaxies in a given redshift tomographic bin ${z}$. These corrections are
averaged for both shear components and, therefore, are applied to both. The
m-correction estimation is detailed in Kannawadi et al. (2019) and Giblin et
al. (2021).

The shear fields can be expressed in terms of a spin-2 field via the E- and
B-modes of spherical harmonic decomposition; the observed shear components are
analogous to the $Q$ and $U$ Stokes parameters, sharing a similar mathematical
structure (Seljak & Zaldarriaga, 1997; Bartelmann & Schneider, 2001). For an
ideal case with a complete coverage of the sky, the shear field decomposition
into spin-2 spherical harmonics, ${}_{\pm 2}Y_{\ell m}$,is defined as

$\displaystyle\gamma(\hat{\mathbf{n}})$
$\displaystyle=\gamma_{1}(\hat{\mathbf{n}})+{\rm
i}\gamma_{2}(\hat{\mathbf{n}})=\sum_{\ell m}\left({E}_{\ell{m}}+{\rm
i}{B}_{\ell{m}}\right)\,_{+2}{Y}_{\ell{m}}(\hat{\mathbf{n}})\,;$ (3)
$\displaystyle\gamma^{*}(\hat{\mathbf{n}})$
$\displaystyle=\gamma_{1}(\hat{\mathbf{n}})-{\rm
i}\gamma_{2}(\hat{\mathbf{n}})=\sum_{\ell{m}}\left({E}_{\ell{m}}-{\rm
i}{B}_{\ell{m}}\right)\,_{-2}{Y}_{\ell{m}}(\hat{\mathbf{n}})\,,$ (4)

where $\gamma(\hat{\mathbf{n}})$ is a complex number and ∗ represents its
complex conjugate. The spherical harmonic coefficients for each mode, over a
region $\Omega_{\hat{\mathbf{n}}}$ of the sky, are

$\displaystyle E_{\ell m}=$ $\displaystyle\frac{1}{2}\int{\rm
d}\Omega_{\hat{\mathbf{n}}}\left\\{\gamma(\bm{\hat{}}{\mathbf{n}})\,_{+2}Y^{*}_{\ell
m}+\gamma^{*}(\hat{\mathbf{n}})\,_{-2}Y^{*}_{\ell m}\right\\}\,;$ (5)
$\displaystyle B_{\ell m}=$ $\displaystyle\frac{-\text{i}}{2}\int{\rm
d}\Omega_{\hat{\mathbf{n}}}\left\\{\gamma(\bm{\hat{}}{\mathbf{n}})\,_{+2}Y^{*}_{\ell
m}-\gamma^{*}(\bm{\hat{}}{\mathbf{n}})\,_{-2}Y^{*}_{\ell m}\right\\}\,.$ (6)

Considering the weak lensing limit, the shear field that emerges from a scalar
gravitational field should be a gradient or a curl-free field, containing
$B_{\ell m}=0$ (Bartelmann & Schneider, 2001; Hikage et al., 2011; Kilbinger,
2015). Although for multiple lenses a B-mode power spectrum can be generated,
its power is expected to be orders of magnitude smaller than the E-mode power
spectrum, meaning that effectively all the cosmological information should
come from the E-mode power spectra and that any significant signals in B-modes
would arise from systematic contamination.

The situation is different for realistic cases when sky coverage is partial.
The survey area is now constrained to a region
$\mathcal{W}(\hat{\mathbf{n}})$; where $\mathcal{W}(\hat{\mathbf{n}})=0$ for
regions that are not observed, masked out due to bright stars or bad seeing,
or contain no galaxies; and 1 otherwise – where there are galaxies in the
given pixel. In other words, $\mathcal{W}(\hat{\mathbf{n}})$ is an effective
binary mask constructed from the catalogue (which already excludes bad regions
as detailed in Giblin et al. 2021). Although not the case for the KiDS-1000
dataset, future applications to Euclid data $\mathcal{W}(\hat{\mathbf{n}})$
need careful consideration on how to deal with a binary mask and shear weights
for high-resolution maps (Euclid Collab. et al., in prep).

For partial sky observations on the celestial sphere, the decomposition of the
observed shear field into spin-2 spherical harmonics is given by (Brown et
al., 2005)

$\displaystyle\tilde{\gamma}_{1}(\hat{\mathbf{n}})\pm\rm{i}\tilde{\gamma}_{2}(\hat{\mathbf{n}})$
$\displaystyle\equiv\mathcal{W}(\hat{\mathbf{n}})\left[\gamma_{1}(\hat{\mathbf{n}})\pm{\rm
i}\gamma_{2}(\hat{\mathbf{n}}))\right]\,$ $\displaystyle=\sum_{\ell
m}\left(\tilde{E}_{\ell m}\pm{\rm i}\tilde{\textit{B}}_{\ell{m}}\right)\,_{\pm
2}Y_{\ell{m}}(\hat{\mathbf{n}})\,,$ (7)

where we denote the partial sky quantities with a tilde. Here, the pseudo
spherical harmonic coefficients, $\tilde{E}_{\ell m}$ and $\tilde{B}_{\ell
m}$, are related to the full-sky E- and B-modes as

$\displaystyle\tilde{E}_{\ell m}$
$\displaystyle=\sum_{\ell^{\prime}m^{\prime}}\left(E_{\ell
m}\@tens{W}^{+}_{\ell\ell^{\prime}mm^{\prime}}+B_{\ell
m}\@tens{W}^{-}_{\ell\ell^{\prime}mm^{\prime}}\right)\,;$ (8)
$\displaystyle\tilde{B}_{\ell m}$
$\displaystyle=\sum_{\ell^{\prime}m^{\prime}}\left(B_{\ell
m}\@tens{W}^{+}_{\ell\ell^{\prime}mm^{\prime}}-E_{\ell
m}\@tens{W}^{-}_{\ell\ell^{\prime}mm^{\prime}}\right)\,,$ (9)

where the mixing kernels $\@tens{W}^{\pm}_{\ell\ell^{\prime}mm^{\prime}}$ are
explained in more detail in the following section.

Now, considering i and j as redshift tomographic bins, we can define an
estimator for the pseudo angular power spectra for cosmic shear as

$\displaystyle\tilde{C}^{E_{i}E_{j}}_{\ell}$
$\displaystyle\approx\frac{1}{2\ell+1}\sum_{m}\tilde{E}^{i}_{\ell
m}\tilde{E}^{j}_{\ell
m}-\left\langle\widetilde{\mathcal{N}}^{i,j}_{\ell}\right\rangle\delta_{ij}\,,$
(10)

with a similar expression for the pseudo B-modes. The last element on the
right-hand side of Eq. (10) is the average shape noise pseudo power spectrum
defined in the following paragraph.

We estimate the noise power spectrum by randomising the orientation of the
galaxy shapes in the catalogue by a random angle $\theta^{R}$, while keeping
the galaxy positions and weights fixed:

$\displaystyle\epsilon_{1}^{R}$
$\displaystyle=\epsilon^{\text{corr}}_{1}\cos(\theta^{R})-\epsilon^{\text{corr}}_{2}\sin(\theta^{R})\,,$
(11) $\displaystyle\epsilon_{2}^{R}$
$\displaystyle=\epsilon^{\text{corr}}_{2}\cos(\theta^{R})+\epsilon^{\text{corr}}_{1}\sin(\theta^{R})\,.$
(12)

This keeps the mask and effective number density of galaxies fixed. Then, we
measured the pseudo E/B-mode angular power spectrum of the randomised
ellipticities catalogue. This step is repeated 100 times and the mean noise
power spectrum is finally subtracted from the data’s pseudo-$C_{\ell}$. Using
Eq. (11) and (12) together with Eq. (2) to define the randomised shear
estimates, $\gamma^{R}$, one can propagate this quantity through Equations (7)
and (8) to obtain $\tilde{E}^{R,i}_{\ell m}$. Finally, the
$\tilde{\mathcal{N}}_{\ell}$ term from Eq. (10) is simply defined as

$\displaystyle\tilde{\mathcal{N}}^{i,j}_{\ell}=\frac{1}{2\ell+1}\sum_{m}\tilde{E}^{R,i}_{\ell
m}\tilde{E}^{R,j}_{\ell m}\,.$ (13)

It is important to note that even if the number of catalogue randomisations
has a very small impact on the measured angular power spectrum, it can have a
significant impact on the covariance estimation (cf. Sect. 3.4).

Although this has been the standard methodology for noise estimates in
harmonic space for cosmic shear surveys (Becker et al., 2016; Hikage et al.,
2019), it could potentially lead to a biased estimate of
$\tilde{\mathcal{N}}_{\ell}$. Starting from Eq. (1) and assuming both the
intrinsic ellipticity and shear are Gaussian distributed, by randomly
orienting galaxies in a shear catalogue one obtains an estimate of
$\sigma^{2}_{\epsilon_{\rm obs}}=\sigma^{2}_{\gamma}+\sigma^{2}_{\epsilon}$,
where the shear field variance, $\sigma^{2}_{\gamma}$, is mixed with the
variance from the intrinsic shape and measurement errors,
$\sigma^{2}_{\epsilon}$. However, using the mocks described in Sect. 3.4, we
have verified that this effect is negligible for the KiDS-1000 cosmic shear
catalogue. This test consisted comparing the estimated noise from
ellipticities and pure shear dispersion from the mocks using the method
described above. This randomised shape catalogue approach for noise power
spectra estimation also proved to be a potential bottleneck for future Euclid
applications; alternatives such as taking the expectation value of the noise
will be implemented for Stage IV surveys.

###  3.2 The mixing matrix

Figure 2: KiDS-1000 full mixing matrix, $\@tens{M}_{\ell\ell^{\prime}}$ from
Eqs. (14) and (15), for the third redshift tomographic bin. For this case,
$\@tens{W}^{++}$ in the diagonal is responsible for most of the mixing for the
E-modes. In a scenario with non-zero but small B-modes, the off-diagonal terms
are still subdominant to the E-mode signal. For each individual block,
$\@tens{W}^{++}$ and $\@tens{W}^{--}$, we show the matrices calculated in the
range $2\leq\ell\leq 3070$. From Eqs. (16) and (17) one can see that these
objects are not symmetric in $\ell$-$\ell^{\prime}$.

The mixing of modes, introduced in Eqs. (8) and (9) via the
$\@tens{W}^{\pm}_{\ell\ell^{\prime}mm^{\prime}}$ terms, is a purely
geometrical effect caused by the empty pixels where data are not observed.
Differently from galaxy clustering analyses, in cosmic shear the absence of
data are taken as the absence of information – if we do not observe lensed
galaxies in a region of the sky, we are unable to estimate shear information
for that region. Therefore, to deal with the mode mixing caused by the
effective geometry of the survey, we construct effective tomographic masks
from the KiDS-1000 cosmic shear catalogue. As previously mentioned, for each
redshift tomographic bin this mask, $\mathcal{W}(\hat{\mathbf{n}})$, will be 1
if there are galaxies assigned to that pixel or zero in case the pixel is
empty. In other words, in our case the mask is dependent only on the empty
pixels for a given tomographic bin.

For a given mask, the expectation value of the pseudo-$C_{\ell}$ estimates
taken over all realisations of noise and cosmic variance can be related to the
underlying full-sky angular power spectra through the mixing matrix
$\@tens{M}_{\ell\ell^{\prime}}$, which contains information about the survey
geometry in harmonic space (Peebles, 1973; Brown et al., 2005; Hikage et al.,
2011):

$\displaystyle\left\langle\tilde{\bm{C}}_{\ell}^{ij}\right\rangle=\sum_{\ell^{\prime}}\@tens{M}_{\ell\ell^{\prime}}^{ij}\bm{C}_{\ell^{\prime}}^{ij}\,,$
(14)

or, assuming no parity violation modes are present
$\left(\tilde{C}_{\ell}^{E_{i}B_{j}}=\tilde{C}_{\ell}^{B_{i}E_{j}}\right)$,

$\displaystyle\begin{pmatrix}\left\langle\tilde{C}_{\ell}^{E_{i}E_{j}}\right\rangle\\\
\left\langle\tilde{C}_{\ell}^{B_{i}B_{j}}\right\rangle\end{pmatrix}=\sum_{\ell^{\prime}}\begin{pmatrix}\@tens{W}^{++}_{\ell\ell^{\prime}}(i,j)&\@tens{W}^{--}_{\ell\ell^{\prime}}(i,j)\\\
\@tens{W}^{--}_{\ell\ell^{\prime}}(i,j)&\@tens{W}^{++}_{\ell\ell^{\prime}}(i,j)\\\
\end{pmatrix}\begin{pmatrix}{C}_{\ell^{\prime}}^{E_{i}E_{j}}\\\
{C}_{\ell^{\prime}}^{B_{i}B_{j}}\end{pmatrix}\,.$ (15)

The individual elements of the mixing matrix in the equation above are
composed of

$\displaystyle\@tens{W}^{\pm\pm}_{\ell\ell^{\prime}}(i,j)=\frac{1}{2\ell+1}\sum_{mm^{\prime}}\@tens{W}^{\pm,i}_{\ell\ell^{\prime}mm^{\prime}}\left(\@tens{W}^{\pm,j}_{\ell\ell^{\prime}mm^{\prime}}\right)^{*}\,,$
(16)

where $\@tens{W}^{\pm}_{\ell\ell^{\prime}mm^{\prime}}$ are the same mixing
kernels present in Eq. (8) and (9). These can be defined for each tomographic
bin as (Hikage et al., 2011)

$\displaystyle\@tens{W}^{\pm,i}_{\ell\ell^{\prime}mm^{\prime}}\equiv$
$\displaystyle\int{\rm d}\Omega_{\hat{\mathbf{n}}}\,{}_{\pm
2}Y_{\ell^{\prime}m^{\prime}}(\hat{\mathbf{n}})\mathcal{W}^{i}(\hat{\mathbf{n}})\,_{\pm
2}Y^{*}_{\ell m}(\hat{\mathbf{n}})\,;$ $\displaystyle=$
$\displaystyle\sum_{\ell^{\prime\prime}m^{\prime\prime}}(-1)^{m}\sqrt{\frac{(2\ell+1)(2\ell^{\prime}+1)(2\ell^{\prime\prime}+1)}{4\pi}}\mathcal{W}^{i}_{\ell^{\prime\prime}m^{\prime\prime}}$
$\displaystyle\times\begin{pmatrix}\ell&\ell^{\prime}&\ell^{\prime\prime}\\\
\pm 2&\mp
2&0\end{pmatrix}\begin{pmatrix}\ell&\ell^{\prime}&\ell^{\prime\prime}\\\
m&m^{\prime}&m^{\prime\prime}\end{pmatrix}\,,$ (17)

where the $2\times 3$ objects above are the Wigner $3j$ symbols, calculated
using UCLWig3j library555https://github.com/LorneWhiteway/UCLWig3j – this
library optimises the calculation of Wigner $3j$ symbols using the recurrence
relation by Schulten & Gordon (1976)., and

$\displaystyle\mathcal{W}_{\ell m}^{i}=\oint{\rm
d}\Omega_{\hat{\mathbf{n}}}\mathcal{W}^{i}(\hat{\mathbf{n}})Y^{*}_{\ell
m}(\hat{\mathbf{n}})\,$ (18)

is the spin-0 spherical harmonic decomposition of the effective binary mask
with the corresponding spin-0 spherical harmonic basis, $Y_{\ell
m}(\hat{\mathbf{n}})$. The mode mixing effect introduced by the mask is then
forward-modelled into the theory data-vector, with more details presented in
Sect. 5.1. An example of the full mixing matrix for the third redshift
tomographic bin is presented in Fig. 2.

Here we would like to point out two important details regarding using a binary
mask as opposed to a mask containing the shear weights or a variable depth
mask. First, a mask containing shear weights, used as
$\mathcal{W}(\hat{\mathbf{n}})$, would cause an imprint from the clustering of
galaxies into the mixing matrix, causing an artificial excess power for
$\mathcal{W}_{\ell m}$ and biasing the forward-modelling of the mixing matrix
into the likelihood (see Sect. 5.1). Although this is not a problem for
KiDS-1000 given the multipole range, HEALPix resolution, and effective galaxy
density ($n_{\rm eff}$), it could be a problem for Stage-IV surveys such as
Euclid. The pixelisation scheme smooths out the weights information if the
resolution is too low for a denser survey. Meanwhile, increasing the
resolution too much can cause many pixels containing a single galaxy only,
which, given Eq. (2), would cause the weights information to be lost and bias
the analysis. Details regarding how to implement the weights information
without causing the aforementioned bias will be investigated in forthcoming
work.

When considering variable depth effects via a non-binary mask, one should
carefully avoid imprinting clustering information into the mask. In this case,
it is not clear how one can include this effect into the estimator or the
data-vector using a non-binary mask. Ideally, to include variable depth
effects, one should take an approach similar to Heydenreich et al. (2020) and
forward model the effect into the likelihood. Yet, Heydenreich et al. (2020)
shows that this effect has very little impact on the inferred cosmological
parameters for a KiDS-like survey. Although we do not include this effect in
the estimator or in the mixing matrix calculation, variable depth effects are
modelled into the covariance via the Egretta mocks (see Sect. 3.4), where
their impact is more significant (Joachimi et al., 2021b). Therefore, the
additional scatter introduced by variable depth is propagated through the
inference by including them directly in the covariance matrix.

###  3.3 Pseudo-$C_{\ell}$ measurements

Figure 3: Binned mixing matrices as band-pass filters and the physical scales
mixed for each redshift tomographic bin edge. Top: redshift as a function of
multipole for different wave-numbers. The coloured horizontal dashed-lines are
the centre of the redshift tomographic bins used in the analysis and shown in
Fig. 1; the black lines are curves of constant wavenumber, $k=\ell/\chi(z)$,
in units of $h$ Mpc-1 for a given $\ell$ and $z$, where $\chi(z)$ is the co-
moving distance. Bottom: the KiDS-1000 pseudo-$C_{\ell}$s binned mixing
matrices (see Fig. 2, for example). These are convolved with the theory
$C_{\ell}$s (see Eq. 23) to model the effect of multipole mixing introduced by
the survey mask. The filters are calculated using all the multipoles inside
the eight log-spaced bandwidths from $76\leq\ell\leq 1500$. Figure 4:
KiDS-1000 E-mode (upper triangle) and B-mode (lower triangle)
pseudo-$C_{\ell}$ measurements with the best-fit model from our cosmological
analysis in Sect. 6.1 convolved with the data’s mixing matrix (Eq. 14). The
auto- and cross-angular power spectra have been measured in eight log-spaced
bandwidths from $76\leq\ell\leq 1500$. The error bars are estimated from the
covariance (described in Sect. 3.4), calculated with Flask (Xavier et al.,
2016) and Salmo (Joachimi et al., 2021b) simulations. We include the best-fit
theory for the pseudo B-modes using Eq. (15) (red line); the amplitude is too
small to be distinguishable from the zero line with a reduced $\chi^{2}\sim 1$
for a null detection.

Our pseudo-$C_{\ell}$ implementation uses pixelised shear maps for each
redshift bin using a HEALPix grid (Górski et al., 2005). For the map
generation, we use $N_{\text{side}}=1024$
($N_{\text{pix}}=12N_{\text{side}}^{2}$), measuring multipoles between
$76\leq\ell\leq 1500$ with the method described above. The lower bound we used
in this work contains larger scales than the band-power measurements from
Asgari et al. (2021)666For comparison, the scales used in Asgari et al. (2021)
for band-power measurements were $100\leq\ell\leq 1500$.. Observing the binned
mixing matrices shown in Fig. 3, we note that by cutting at an bellow an
$\ell_{\rm min}=100$, we find that just a small number of lower multipoles get
mixed into the scales of interest. Therefore, we decreased the $\ell$ lower
bound towards larger scales. Although this choice is driven by survey geometry
as below these scales practically no information is in the data, the specific
value of $\ell_{\rm min}=76$ was chosen arbitrarily. We were able to push to
slightly lower $\ell$ values than the band-powers analysis from (Asgari et
al., 2021) as the latter suffers strong mixing bellow $\ell=100$. Figure 3
also gives us an intuition for the $k$-scales that enter our analysis as a
function of redshift and multipole given the bandwidth $\ell$-bins used.

Following Brown et al. (2005)777We note that this choice of weight could be
optimised for cosmic shear analysis, which will be investigated in future
work., the measured pseudo-$C_{\ell}$s are then binned using eight log-spaced
multipole bins between $76\leq\ell\leq 1500$,

$\displaystyle\tilde{C}^{ij}_{L}=\frac{1}{2\pi}\sum_{\ell}\frac{\ell(\ell+1)}{(\ell^{L+1}-\ell^{L})}\tilde{C}_{\ell}^{ij}\,,$
(19)

where $\ell^{L}<\ell\leq\ell^{L+1}$ are the boundaries of the bandwidth bins
centred in $L$, shown in Fig. 3 as the vertical dotted lines. The measured
pseudo-$C_{\ell}$ estimates are shown in Fig. 4 for the E- and B-modes,
together with the best-fit theory line from the cosmological constraints shown
in Sect. 6.1. The pseudo-$C_{\ell}$ B-modes shown in the lower triangle part
of Fig. 4 are consistent with zero with a $\chi_{\rm red}^{2}=1.07$ for a
null-detection, indicating that even with the mode-mixing introduced by the
mixing matrix, the contribution of E-modes mixed into B-modes is not enough to
bias the cosmological measurements. The auto-power spectrum for $E_{2}$ and
combinations with it, mostly $E_{2}-E_{3}$, have a slightly higher amplitude
than the best-fit theory lines. A similar effect was found in previous
KiDS-1000 studies but Asgari et al. (2021) have shown that this has almost no
impact on the cosmological constraints nor the goodness-of-fit, as shown in
Fig. 7 of Asgari et al. (2021).

The number of multipole bins was arbitrarily selected so that the overlaps
between the binned mixing matrices (filters), shown in Fig. 3, did not exceed
10% of their integrated areas. This was done to prevent excessive correlations
between adjacent multipole bins. However, it should be pointed out that since
we estimate the covariance matrix used for our cosmological analysis using
simulations (see Section 3.4 for details), we expect any residual correlations
caused by these overlaps to be well described in our analysis.

###  3.4 Simulations and covariance matrix estimation

While analytical covariance models are widely used in cosmic shear
(Efstathiou, 2004; Joachimi et al., 2008; Krause & Eifler, 2017; Alonso et
al., 2019; Hikage et al., 2019; García-García et al., 2019; Li et al., 2019;
Krause et al., 2021; Asgari et al., 2021; Joachimi et al., 2021b; Heymans et
al., 2021; Nicola et al., 2021) and it will be used for future Euclid
applications (Upham et al., 2022), the forward-modelling approach for the
mixing matrix makes it convenient to estimate the covariance matrix from
simulations. We decided to use the Egretta validated suite of 1000 log-normal
simulations generated with Flask888http://www.astro.iag.usp.br/~flask/ (Xavier
et al., 2016) and Salmo999https://github.com/Linc-tw/salmo presented in
Joachimi et al. (2021b). This suite of mocks includes a number of known
observational complexities such as variable depth source redshift
distributions, and spatially varying galaxy sample properties. Flask + Salmo
mocks like Egretta were used to validate the covariance matrix modelling in
Joachimi et al. (2021b), as well as cosmic shear signal recovery, demonstrated
in Heydenreich et al. (2020) using semi-analytic models. Contrary to N-Body
simulations, log-normal simulations are much faster to obtain, allowing for
quick survey simulations that accurately reproduce the desired 2pt functions
in just a few CPU hours. Figure D.9 from Joachimi et al. (2021b) shows the
agreement between the covariances obtained by Egretta mocks and the analytical
covariances. Both methods agree to a few percent for cosmic shear 2pt
functions, with the largest deviations coming from the first redshift bin. In
other words, we are confident that our log-normal simulation-based covariance
is accurate and that it reproduces, up to a few percent difference, what is
expected by theoretical covariance calculations for cosmic shear analysis.

Table 1: Comparison between the original mocks, the sub-sampled mocks and the
Gold Sample SOM

Bin |  $\bar{z}$ (J21) |  $\bar{z}$ (Sub-sampled) |  $\bar{z}$ (Gold Sample) |  $n_{\text{eff}}$ (J21) |  $n_{\text{eff}}$ (Sub-sampled) |  $n_{\text{eff}}$ (Gold Sample)
---|---|---|---|---|---|---
1 | 0.39 | 0.40 | 0.26 | 0.85 | 0.62 | 0.62
2 | 0.49 | 0.47 | 0.40 | 1.56 | 1.19 | 1.18
3 | 0.67 | 0.62 | 0.56 | 2.23 | 1.86 | 1.85
4 | 0.83 | 0.80 | 0.79 | 1.52 | 1.27 | 1.26
5 | 1.00 | 0.96 | 0.98 | 1.39 | 1.32 | 1.31

101010Here we present a comparison between characteristics from the sub-
sampled and original redshift distributions from Egretta mocks with the Gold
Sample SOM $N(z)$s from (Hildebrandt et al., 2021). The original mock values
from Joachimi et al. (2021b) are labelled as J21 and values labelled as sub-
sampled are the ones used in this work.

Egretta mocks were generated using a combination of log-normal matter
distribution simulations from Flask with post-processing using Salmo to
populate the matter density field with galaxies based on their local redshift
distribution and properties. The mock creation begins by generating the dark
matter density field using angular power spectra calculated from 18 redshift
tomographic bins of similar intervals, in comoving distances of 150 – 200
Mpc$\,h^{-1}$, using the fiducial cosmology presented in Table A.1 from
Joachimi et al. (2021b). Flask then integrates along the line-of-sight to
obtain the shear and convergence fields. Using Salmo, galaxies are positioned
in the sky according to the underlying matter density field and following a
Poisson distribution, with their shapes given by the weak lensing fields
generated by Flask. In addition, Salmo allows for spatially variable redshift
distributions to be implemented, reproducing variable depth effects in the
final simulated catalogues (more details in Section 4 of Joachimi et al.
2021b).

Even though the Egretta mocks were created with KiDS-1000 analysis in mind,
they were originally generated using the KiDS+VIKING-450 (KV450) DIR redshift
distributions (Hildebrandt et al., 2017, 2020). In order to use these mocks
for our purposes, we have sub-sampled them to match the KiDS-1000 SOM redshift
distributions (Hildebrandt et al., 2021). The KV450 galaxy sample uses all the
galaxies available, whereas in KiDS-1000 we remove galaxies using the SOM-
based Gold sample. For this reason, we can sub-sample the Egretta mocks; we
have a lower $n_{\rm eff}$ for KiDS-1000 than for KV450. Table 1 shows the
sub-sampling from the original mocks to match the KiDS-1000 gold catalogue
galaxy density. Although the mean redshift, $\bar{z}$, is slightly higher than
the SOM’s redshift distributions for the first two tomographic bins, the
effective number of galaxies, $n_{\text{eff}}$, matches quite well. As higher
mean redshifts lead to a higher variance in the sample’s pseudo-$C_{\ell}$s,
the differences in the mean redshifts yield more conservative errors and,
therefore, are not a concern.

The mocks were originally created at $N_{\text{side}}=4096$; however, we apply
the pseudo-$C_{\ell}$ estimator at $N_{\text{side}}=1024$ to match the process
applied to the data (see Sect. 3.1 for details), measuring angular power
spectra using the same eight log-spaced $\ell$-bins. The difference in
resolution is taken into account by correcting for the pixel window function
at the target $N_{\rm side}=1024$ (Górski et al., 2005). Lastly, we calculate
the covariance of the noise subtracted pseudo-$C_{\ell}$s for all 1000 Egretta
mocks. It is fundamental that the noise estimates from each mock realisation
are calculated in the same way as in the data: for each mock, we re-orient
galaxies in the catalogue, calculate the noise angular power spectrum to
obtain a noise realisation, and subtract the average of 100 noise realisations
from the measured pseudo-$C_{\ell}$ for that mock catalogue. We also note that
the number of noise realisations used to estimate the noise for each of the
mocks can have a significant impact on the covariance matrix, which directly
impacts the cosmological constraints since we have a noisy estimate of the
covariance. We performed tests using only 20 realisations and the impact on
cosmological parameters was significant with an increase as high as 28% on the
$S_{8}$ error-bars. Therefore, we stress that a high number of noise power
spectra estimates is crucial for accurate estimation of the covariance matrix.

We also verified that, for the mocks, the impact of the mixing matrix
variation between each mock realisation is subdominant in the
pseudo-$C_{\ell}$s. We are forward-modelling the mixing matrix into the
theory, meaning that we include the partial sky effects into the mock’s
measurements. Although the mask varies slightly from mock to mock, the
pseudo-$C_{\ell}$s effects due to the masking vary within the mock’s cosmic
variance and noise. We calculated the variation of the mixing matrix from a
sub-sample of the mocks and verified that, overall, they are similar to the
data’s mixing matrix, with a maximum of a few percent difference. Moreover, we
also verified that the variation is subdominant when convolving theory-vectors
with the mocks’ mixing matrix and comparing with theory-vectors convolved with
the data’s mixing matrix, showing a few percent difference.

The Egretta mocks have no multiplicative bias, being equivalent to the
KiDS-1000 data after the m-correction is applied. The m-correction has an
uncertainty associated with it, which needs to be included in the covariance
as an additional term. To propagate correctly the effects of the
multiplicative shear calibration into our covariance matrix we take an
analytical approach (Joachimi et al., 2021b), modelling the m-correction
contribution to the covariance matrix as (Bridle et al., 2002; Taylor &
Kitching, 2010; Doux et al., 2021; Joachimi et al., 2021b)

$\displaystyle{\rm
Cov}_{m}\left(\tilde{C}^{i,j}_{\ell};\tilde{C}^{i^{\prime},j^{\prime}}_{\ell^{\prime}}\right)=$
$\displaystyle\,\tilde{C}^{i,j}_{\ell}\tilde{C}^{i^{\prime},j^{\prime}}_{\ell^{\prime}}\left[\sigma_{m}^{i}\sigma_{m}^{i^{\prime}}+\sigma_{m}^{i}\sigma_{m}^{j^{\prime}}\right.$
$\displaystyle\left.\quad+\,\sigma_{m}^{j}\sigma_{m}^{i^{\prime}}+\sigma_{m}^{j}\sigma_{m}^{j^{\prime}}\right]\,,$
(20)

where $\sigma_{m}^{i}$ is defined as the standard deviation of the
m-corrections shown in Table 4 and $\tilde{C}^{E_{i}E_{j}}_{\ell}$ are theory
angular power spectra convolved with the mixing matrix using Eq. (14). This
expression assumes that the expectation value of m is zero after the
correction, with $\sigma^{i}_{m}\ll 1$. Since the $\sigma_{m}$ part of Eq.
(20) is independent of $\ell$, the mixing of modes does not affect this part
of the covariance. Asgari et al. (2021) demonstrated that this approach gives
equivalent results to having the $m_{z}$ correction as a free parameter, for
each redshift tomographic bin z, during the cosmological inference analysis.

Finally, our total covariance matrix is the combination of the simulated
covariance from the Egretta mocks and the analytic m-correction covariance:
$\rm{Cov}_{\rm tot}=\rm{Cov}_{\rm mocks}+{\rm Cov}_{\textit{m}}$. The
resulting covariance, $\rm{Cov}_{\rm tot}$, is presented as a correlation
matrix in Fig. 5 and shows a similar structure as previously found by Joachimi
et al. (2021b) for the band-powers analysis.

Figure 5: Correlation matrix for the E-mode pseudo-$C_{\ell}$ covariance. The
covariance matrix was calculated using 1000 mocks from the Egretta suite of
Flask + Salmo simulations (Joachimi et al., 2021b). The indices $(i,j)$ in the
labels are the redshift bin numbers for each pair of $C^{i,j}_{\ell}$ E-modes.

##  4 Systematic contamination analysis

Figure 6: P-values for the cross-correlations between the data’s tomographic
pseudo E-modes and different systematics for a null detection ($p\leq 0.05$
would indicate a detection). Here we show the systematics motivated in Sect.
4: object detection threshold ($\mu$-threshold), extinction for the r-band,
magnitude limit for the r-band, PSF ellipticities components in harmonic space
(E/B-mode), the PSF size, and stellar overdensity from Gaia DR2 (Gaia
Collaboration et al., 2018). The highest correlation occurs for the PSF E-mode
in the fifth tomographic bin, this is shown in Fig. 7 to be subdominant to the
cosmological signal given the estimated statistical error in this particular
tomographic bin.

Photometric galaxy surveys are affected by observational effects related to
instrumental, astrophysical, and atmospheric phenomena. When measuring the
shapes and positions of galaxies, calibration and control of such systematic
effects are at the core of a successful analysis. Like any galaxy survey, the
KiDS-1000 data set experiences systematic contaminations given the complex
nature of combining optical and infrared data spread through nine different
bands. For example, instrumental effects such as the telescope’s point spread
function can vary due to telescope inclination, impacting galaxy shape
measurements resulting in the introduction of spurious and artificial
correlations in the final data-vector. Other observational effects like a
spatially varying magnitude limit in a given band also have the potential to
contaminate or bias the sample, inducing an artificial selection, for example.
Finally, astrophysical phenomena like Galactic extinction and stellar
contamination in the galaxy sample can impact the quality of data, resulting
in correlations that alter the data’s measured angular power spectra leading
to a potential bias in the inferred cosmology.

An effort to select, clean, and validate this sample was undertaken by Giblin
et al. (2021), including a range of contamination checks. In this section, our
goal is to analyse any possible residual systematic contamination in the
KiDS-1000 catalogue that could affect the measured angular power spectra of
galaxy ellipticities. We perform null tests by cross-correlating data with
relevant observational quantities in the catalogue and stellar density from
Gaia Data Release 2 (Gaia Collaboration et al., 2018). We then observe the
goodness-of-fit to a zero-correlation case. The approach taken here is similar
to null tests in Leistedt & Peiris (2014), but without the use of mode-
projection.

A summary of observational, astrophysical, and instrumental systematics
analysed is as follows:

  * •

Extinction ($ugriZYJHK_{s}$): Galactic extinction for each of the observed
nine bands derived from the $E(B-V)$ maps from Schlegel et al. (1998) combined
with the extinction coefficients, $R_{V}=3.1$, from Schlafly & Finkbeiner
(2011). The extinction coefficients for each filter, $R_{f}$, follow a linear
relation between the observed bands: $R_{f}=A_{f}/E(B-V)$. Therefore, we have
only shown results for the r-band since this is the detection band for
KiDS-1000 and other values are just linearly scaled based in this band. The
values for $R_{f}$ used for KiDS-1000 can be found in Table 4 of Kuijken et
al. (2019).

  * •

Magnitude limit ($ugriZYJHK_{s}$): magnitude limit in each of the observed
nine bands. This quantity follows a forced photometry based on the r-band
detection. We only show results for the r-band as it is the detection band
used for forced photometry.

  * •

$\mu$-Threshold: object detection threshold above the background.

  * •

PSF ellipticities: mean PSF ellipticity components, treated as a spin-2 field,
converted to PSF E/B-modes in harmonic space (see Eqs. 8 and 9). The value
quoted in Fig. 6 has the PSF $\alpha$-correction term from Giblin et al.
(2021) applied to the PSF ellipticities before converting them to harmonic
space. The subsequent text explores and explains this in detail.

  * •

PSF Size: The trace of the PSF magnitude moments, $T=Q_{11}+Q_{22}$, where
$Q_{ij}$ is defined in terms of the second-order moments of the two-
dimensional angular light distribution, $I(\hat{\mathbf{n}})$; see Sect. 3.1
of Giblin et al. (2021).

  * •

Stellar Density: star sample from Gaia DR2111111Obtained from
https://cds.u-strasbg.fr/gaia containing all Gaia DR2 objects in a HEALPix
grid of $N_{\rm side}=1024$.

For each of the quantities mentioned above, we create HEALPix maps using the
catalogue object positions with $N_{\rm side}=1024$ and the mean value for the
observables in each pixel. A spherical harmonic transform is applied to the
maps — spin-0 for all scalar quantities and spin-2 transforms for the PSF
ellipticities — in order to obtain the systematic’s spherical harmonic
coefficients, $a_{\ell m}^{\text{syst}}$, and the data’s spherical harmonic
coefficients, $X_{\ell m}^{i}$ – where $\text{X}=\\{E,\,B\\}$ and i is the
tomographic bin index. Using the spherical harmonic coefficients, we estimate
the cross-power spectra between data and the systematic quantities of interest
as

$\tilde{C}^{\text{X}_{i},\text{syst}}_{\ell}=\frac{1}{(2\ell+1)}\sum_{m}\text{X}_{\ell
m}^{i}\left(a_{\ell m}^{\text{syst}}\right)^{*}\,.$ (21)

To estimate the errors for the cross-correlations between data and
systematics, we obtain covariances in two different ways. For the catalogue
quantities, we use a jackknife re-sampling method. We start by subdividing the
survey footprint into 60 distinct regions using a K-Means algorithm on the
sphere (Kwan et al., 2017)121212Our K-Means code uses the code kmeans_radec,
developed by Erin Sheldon which can be found here:
https://github.com/esheldon/kmeans_radec.. While removing the kth region, we
calculate the pseudo-$C_{\ell}$ (Eq. 21) of each sub-sample for all 60 regions
in the footprint. For the stellar density Gaia map, we take a different
approach. The KiDS-1000 footprint is naturally over a low stellar density
region in the sky, meaning that the jackknife covariance method leads to very
small variations in the cross-power spectra. This is due to most of the
jackknife regions containing a very small number of stellar objects in the
Gaia DR2 map, if any. Therefore, the covariances for the cross-power spectra
between stellar density and the data are calculated using the cross-
correlations between the Gaia DR2 map and 200 noise realisations of the galaxy
ellipticities in the catalogue (created by randomly orienting the galaxies in
the KiDS-1000 catalogue for each realisation). Given the stochastic nature of
this process, it yields a somewhat underestimated covariance, which leads to a
conservative null test.

Figure 7: PSF E-mode and $E_{5}$ cross-power spectrum. (Top panel) Cross-
angular power spectra between PSF E-mode component and $E_{5}$ with the
$\alpha$-correction from Giblin et al. (2021) applied to the PSF
ellipticities. (Bottom panel) Signal-to-noise ratio as defined by Eq. (32),
showing that the largest S/N multipole bandwidth has $\text{S/N}\sim-2.5\times
10^{-2}$, meaning that it is well within the estimated covariance for our
analysis.

Subsequently, we calculate the goodness-of-fit using the p-value for a fit to
zero correlations between systematics and measured ellipticities. This is
shown for the E-modes in Fig. 6 and in Fig. 15 for the B-modes as heat-maps.
From Fig. 4 we do not detect any pseudo-$C_{\ell}$ B-modes; therefore, the
B-modes systematics analysis are shown in Appendix B.

Initially, we find a significant correlation with the PSF ellipticity and the
5th tomographic E-mode, with a p-value of $0.004$ for a cross-correlation
detection. This result is to be expected from the analysis of Giblin et al.
(2021), who quantified the level of contamination to the cosmic shear signal
from PSF leakage. Giblin et al. (2021) measured the $\alpha$ parameter per
tomographic bin, where $\alpha$ quantifies the average fraction of PSF
ellipticity in the shear estimator. As the KiDS PSF has a very low ellipticity
on average, and as $\alpha$ was found to be small (at the level of 4% for
$z_{\rm B}>0.7$), or consistent with zero (for $z_{\rm B}<0.7$), Giblin et al.
(2021) concluded that the statistically significant systematic PSF leakage
introduced less than a $0.1\sigma$ in the inferred $S_{8}$ constraints for
KiDS-1000, inducing no significant biases.

In this analysis, we recover the angular dependency on this systematic using
the $\alpha$ correction for the PSF ellipticity measurements, finding a cross-
correlation between the fifth tomographic bin E-mode with a p-value of $0.027$
for a non-zero signal. Although the value is not worryingly low, the fifth
tomographic bin contains most of our constraining power (Asgari et al., 2021;
Heymans et al., 2021). Therefore, it is crucial to understand if such
correlations could indeed mimic a cosmological signal and bias our final
cosmological analysis. We show this cross-angular power spectrum in detail in
Fig. 7 together with its signal-to-noise ratio (S/N, Eq. 32) for the cross-
correlations between data and PSF E-mode and the data’s covariance estimated
in Sect. 3.4. The S/N demonstrates that this cross-correlation is subdominant
in the E-mode signal from the fifth tomographic bin, with the bandwidth with
the highest S/N around $5\times 10^{-2}$. This means that this excess
correlation is subdominant given the data’s estimated covariance in this
tomographic bin, leading us to draw the same conclusions as Giblin et al.
(2021) in that this correlation has a sufficiently low significance that it
has no impact on the cosmological parameter inference.

In summary, our conclusion from the examined potential systematics is that no
considerable contamination is biasing the measured E-mode angular power
spectra shown in Fig. 4. Therefore, the estimated cosmological parameters
shown in Section 6.1 are purely a result of the cosmological and astrophysical
signals captured by our measured pseudo-$C_{\ell}$s.

##  5 Cosmological inference

This section describes the final elements necessary for cosmological inference
using tomographic cosmic shear angular power spectra data: theory forward
modelling, the likelihood and the priors we use to obtain the posteriors
presented in Section 6.1. Here, we also describe the external clustering data
included to obtain combined constraints using galaxy clustering information
from BOSS & eBOSS with Lyman-$\alpha$ forest data from eBOSS.

Cosmological inference is performed using a version of the Monte Python 3
software (Audren et al., 2013; Brinckmann & Lesgourgues,
2019)131313https://github.com/brinckmann/montepython_public modified to
directly sample from Gaussian priors141414The modified version sampling from
Gaussian priors can be found at:
https://github.com/BStoelzner/montepython_public/tree/gaussian_prior for the
Multinest sampler (Feroz & Hobson, 2008; Feroz et al., 2009; Buchner et al.,
2014; Feroz et al., 2019). It has also been modified to sample directly from
the $S_{8}$ parameter according to the findings reported in Joachimi et al.
(2021b). The Monte Python 3 pipeline was compared with KCAP, the pipeline used
in Asgari et al. (2021) and Heymans et al. (2021), and results were found to
be consistent.

Following Joachimi et al. (2021b), we report our constraints in two different
ways: the traditional marginalised constraints, extracted from the 1D
marginalised posteriors for each individual parameter, denoted marginal; and
the maximum a posteriori (MAP) with the projected joint highest posterior
density (PJ-HPD); this is also referred to as MAP+PJ-HPD. This method ensures,
by definition, that the multi-dimensional MAP is always within the PJ-HPD
region and it has been demonstrated to accurately report the posterior
distributions for each parameter, impacting significantly the reported results
from skewed distributions. Monte Python originally does not deal with Gaussian
priors and its minimisation algorithm acts directly on the likelihoods, making
the MAP estimates from minimising the log-posterior unreliable. We take a
different approach than the one proposed by Joachimi et al. (2021b) by finely
sampling the posterior using a high number of live-points and obtaining the
MAP estimates from the MultiNest chains instead of minimising the log-
posterior. Tests show that the MAP values we recover have a better goodness-
of-fit from this approach than the minimisation approach. This means our
estimate of the MAP is not as accurate as it could be; however, there is no
impact for the PJ-HPD boundaries. For more details, see Section 6.4 from
Joachimi et al. (2021b).

###  5.1 Theoretical modelling

The modelling approach for the cosmic shear angular power spectra is very
similar to the one described in Joachimi et al. (2021b) with a few technical
differences and additional steps related to the forward-modelling of the
mixing matrix. Differently from Joachimi et al. (2021b), Monte Python 3
obtains the underlying matter power spectra from Class (Lesgourgues, 2011;
Blas et al., 2011).151515https://lesgourg.github.io/class_public/class.html
Following the other KiDS-1000 analyses, we fix the sum of neutrino masses to
$\sum m_{\nu}=0.06{\rm eV}c^{-2}$, using the one massive species approximation
of the normal hierarchy. Not only was this the approach taken by Planck
Collaboration et al. (2020a), but it was shown in Loureiro et al. (2019) that
it has no impact on the standard flat $\Lambda$CDM cosmological parameters
compared to more complex neutrino mass models.

The non-linear matter power spectrum, $P_{\rm m}^{\rm nl}(k)$, is calculated
using the halo model from HMCode (Mead et al., 2015) included in
Class.161616Mead et al. (2021) introduce a new version of the HMCode, with
improved accuracy. Tröster et al. (2021) compared this new HMCode against the
version used in this work and found that it introduces a $0.26\sigma$ shift in
$S_{8}$. This exemplifies how the non-linear modelling is the limiting
theoretical systematic in cosmic shear analysis. This approach incorporates
AGN baryonic feedback into the matter power spectrum using a model with two
parameters: the halo bloating parameter, $\eta_{0}$, and the halo mass-
concentration relation amplitude, $A_{\rm bary}$. These parameters can be
constrained by the following relation from Joudaki et al. (2018):
$\eta_{0}=0.98-0.12A_{\rm bary}$. Therefore, $A_{\rm bary}$ can be used as a
single free parameter to model baryonic effects.

Subsequently, the tomographic cosmic shear angular power spectra are modelled
using projections along the line-of-sight of the 3D matter power spectrum in
redshift bins i and j. The observed $C_{\ell}$s are a combination of
contributions from cosmic shear, indexed by $\gamma$, and the intrinsic
alignment of galaxies, indexed by I (Bernstein, 2009; Joachimi & Bridle,
2010):

$\displaystyle
C^{\epsilon_{i},\epsilon_{j}}_{\ell}=C^{\gamma_{i},\gamma_{j}}_{\ell}+C^{\gamma_{i},\text{I}_{j}}_{\ell}+C^{\text{I}_{i},\gamma_{j}}_{\ell}+C^{\text{I}_{i},\text{I}_{j}}_{\ell}\,.$
(22)

Using the Limber approximation (Kaiser, 1992; LoVerde & Afshordi, 2008;
Kitching et al., 2017; Kilbinger et al., 2017; Lemos et al., 2017), where the
wavenumber k and the radial co-moving distance $\chi$ can be related as
$kf_{\rm K}(\chi)=\ell+1/2$, we can express the right-hand side $C_{\ell}$s
from Eq. (22) as

$\displaystyle C_{\ell}^{X_{i},Y_{j}}=\int_{0}^{\chi_{\rm H}}{\rm
d}\chi\frac{W_{X}^{i}(\chi)W_{Y}^{j}(\chi)}{f_{\rm K}^{2}(\chi)}P_{\rm m}^{\rm
nl}\left(\frac{\ell+1/2}{f_{\rm K}(\chi)},z(\chi)\right)\,,$ (23)

where the $X,Y$ indices are $\gamma$ and/or I, the $i,j$ indices are the
tomographic bins indices for the different kernels, $\chi_{\rm H}$ is the co-
moving distance to the horizon and $f_{\rm K}(\chi)$ is the Friedmann-
Lemaître-Robertson-Walker (FLRW) co-moving angular diameter distance. Note,
also, that we only consider flat $\Lambda$CDM models in this study, therefore,
$f_{\rm K}(\chi)=\chi$. The weak lensing kernel is defined as:

$\displaystyle W_{\gamma}^{i}(\chi)=\frac{3}{2}\frac{H_{0}^{2}\Omega_{\rm
m}}{c^{2}a(\chi)}f_{\rm K}(\chi)\int_{\chi}^{\chi_{\rm H}}{\rm
d}\chi^{\prime}n^{i}_{s}(\chi^{\prime})\frac{f_{\rm
K}(\chi^{\prime}-\chi)}{f_{\rm K}(\chi^{\prime})}\,.$ (24)

The intrinsic alignments kernel is calculated using the non-linear alignment
(NLA) model from Bridle & King (2007) without redshift evolution,

$\displaystyle W^{i}_{\rm I}(\chi)=-A_{\rm IA}\frac{C_{1}\rho_{\rm
cr}\Omega_{\rm m}}{D[a(\chi)]}n_{s}^{i}(\chi)\,,$ (25)

with $C_{1}\rho_{\rm cr}\approx 0.0134$ (Joachimi et al., 2011) and where
$D[a(\chi)]$ is the growth function and $A_{\rm IA}$, the amplitude of
intrinsic alignments, kept as a free parameter in the analysis but the same
for all redshift tomographic bins. Although not physically motivated in the
way non-linearities are accounted for, the NLA model was found by Fortuna et
al. (2021) to be precise enough for the case we are currently analysing.

The theory angular power spectra from Eq. (22) are then corrected for the
pixel window function (Górski et al., 2005), convolved with the mixing matrix
using Eq. (14), assuming $C_{\ell}^{B_{i},B_{j}}=0$, and binned in eight log-
space bandwidths between $76\leq\ell\leq 1500$ using Eq. (19). Although not
individually binned for the analysis, the binned mixing matrix can be seen as
a band-pass filter, showing which scales get mixed for each of the log-spaced
$\ell$-bins used in the data as in Fig. 3. This final convolved and binned
theory vector is then compared to the measurements in the likelihood for
cosmological inference.

###  5.2 External data from SDSS

Since our main objective in this study is to extract cosmological information
from late-time probes using galaxy surveys, we have combined our cosmic shear
measurements with clustering information. For the current generation of
surveys, weak lensing alone has little constraining power beyond the $S_{8}$
parameter and it is limited by the degeneracy in the $\sigma_{8}\,-\Omega_{\rm
m}$ plane (Hall, 2021). In order to compete with the extremely precise cosmic
microwave background measurements from Planck, this synergy between late-time
probes is crucial. By combining late-time probes we can have a better
understanding of the growth of structure tension with early-time probes.

We use the publicly available SDSS
likelihoods171717https://svn.sdss.org/public/data/eboss/DR16cosmo/tags/v1_0_0/likelihoods/BAO-
plus, which combine pre-reconstruction full-shape181818In this context, the
so-called full shape means that the relevant BAO and RSD quantities are fitted
to the full shape of the correlation function or power spectrum and then used
to constrain cosmology. results and post-reconstruction BAO results, both from
the BOSS and eBOSS LRG samples (Alam et al., 2017; Gil-Marín et al., 2020;
Bautista et al., 2021). We also include BAO results from the Lyman-$\alpha$
(Ly$\alpha$) forest auto-correlation and its cross-correlation with quasar
positions (du Mas des Bourboux et al., 2020). We chose to leave the quasar
auto-correlations out, as sampling their biases could have an impact on the
growth of structure parameter, $S_{8}$, since this parameter is degenerate
with the bias.

Using Flask + Salmo mocks, Joachimi et al. (2021b) showed that the BOSS and
KiDS-1000 data sets can be considered independent, with their cross-covariance
consistent with zero. At the same time, even though the BOSS and KiDS-1000
footprints have a small overlap, Heymans et al. (2021) found that the galaxy
galaxy-lensing correlation between the two has very little impact on the
cosmological results, mostly due to the necessarily aggressive scale cuts in
the data-vector due to bias and IA model limitations. Yet, it does have a
small impact on the nuisance parameters, $A_{\rm bary}$ and $A_{\rm IA}$, as
shown in Fig. 6 from Heymans et al. (2021). Therefore, we also consider the
two data sets, KiDS-1000 and BOSS, as independent and do not consider any
galaxy-galaxy lensing correlations in this work. We use the already-existing
Monte Python BOSS full-shape likelihood, but update it with eBOSS data.
Following Alam et al. (2021), we use BOSS results from the two low-redshift
($0.2<z<0.6$) LRG bins, and the higher redshift ($0.6<z<1$) eBOSS LRG bin that
combines both BOSS and eBOSS data. For each of these bins we use measurements
of $D_{\rm M}(z_{\text{eff}})/r_{\rm d}$, $D_{\rm H}(z_{\text{eff}})/r_{\rm
d}$ and the rate of structure growth,
$f(z_{\text{eff}})\sigma_{8}(z_{\text{eff}})$. Here, $D_{\rm
M}(z)=(1+z)^{-1}f_{\rm K}\left[\chi(z)\right]$ is the comoving angular
diameter distance, $D_{\rm H}(z)=c/H(z)$, $r_{\rm d}$ is the size of the sound
horizon at the end of the drag epoch, $f(z_{\text{eff}})$ is the linear growth
rate, and $z_{\text{eff}}$ is the effective redshift of the measurement (as
defined by Eq. 7, Ross et al. 2014).

For the Ly$\alpha$ forest measurements we use the Monte Python DR14
likelihoods (Cuceu et al., 2019), but updated with DR16 data (du Mas des
Bourboux et al., 2020). In this case we have measurements of $D_{\rm
M}(z_{\text{eff}})/r_{\rm d}$ and $D_{\rm H}(z_{\text{eff}})/r_{\rm d}$, for
both the auto-correlation and cross-correlation. We treat the two as
independent, following a study on synthetic data sets (du Mas des Bourboux et
al., 2017) that found negligible correlations. From now on, we denote this
combined data set as SDSS BAO and RSD.

Finally, we note that in contrast to the approach for the BOSS data used in
Heymans et al. (2021) and Tröster et al. (2020), the SDSS likelihoods and data
we use here do not assume a flat $\Lambda$CDM model, that is they are not
constraining cosmology using the full galaxy power spectrum. Instead, the SDSS
likelihoods we are using in this work measure the quantities mentioned above –
$D_{\rm M}(z_{\text{eff}})/r_{\rm d}$, $D_{\rm H}(z_{\text{eff}})/r_{\rm d}$,
$f(z_{\text{eff}})\sigma_{8}(z_{\text{eff}})$ – fitted from the measured
correlation function or power spectrum estimates, to constrain cosmology. This
approach is considered to be ‘model-independent’ as these quantities do not
assume a flat $\Lambda$CDM power spectrum when measured. Nonetheless, this
‘model-independent’ approach sacrifices precision as pointed out by Tröster et
al. (2020).

###  5.3 Likelihood and priors

Table 2: Priors on cosmological and nuisance parameters.

Parameter name | Symbol | Prior
---|---|---
Density fluctuation amp. | $S_{8}$ | [0.1, 1.3]
Cold dark matter density | $\Omega_{\text{c}}h^{2}$ | [0.051, 0.255]
Baryonic matter density | $\Omega_{\rm b}h^{2}$ | [0.019, 0.026]
Hubble constant | $h$ | [0.64, 0.82]
Spectral index | $n_{\rm s}$ | [0.84, 1.1]
Intrinsic alignment amp. | $A_{\text{IA}}$ | [$-6,6$]
Baryonic feedback amp. | $A_{\rm bary}$ | [2,3.13]
Redshift displacement | $\delta_{z}$ | $\mathcal{N}(\bm{\mu},\@tens{C}_{z})$
191919The $\Omega_{\rm b}h^{2}$ prior, motivated by BBN constraints from Olive
& Particle Data Group (2014), is a conservative (5$\sigma$) equivalent to
combining BBN data when external data from SDSS is used, similar to what is
shown in Cuceu et al. (2019). Here, $\mathcal{N}(\bm{\mu},\@tens{C}_{z})$ is a
normal distribution with mean $\bm{\mu}$ and covariance $\@tens{C}_{z}$ as
estimated by Hildebrandt et al. (2021).

In this section we discuss the likelihood and prior choices made for our
analysis. Consistency with the main KiDS-1000 analysis (Asgari et al., 2021;
Heymans et al., 2021) is at the core of our choices. One of our objectives in
this work is to show the accuracy of a pseudo-$C_{\ell}$ methodology with
forward-modelling which could be applied to future Stage-IV galaxy surveys. We
want to extend the comparison with other 2pt statistics estimators from Asgari
et al. (2021). This means that making similar prior choices is important as it
allows for an easy and fair comparison.

Our covariance matrix is partially estimated from simulations with an
analytical part for the m-correction error propagation. Therefore, we opt for
a Gaussian likelihood instead of a multivariate t-distribution (Hotelling,
1931) as suggested by Sellentin & Heavens (2016) or the exact
pseudo-$C_{\ell}$ likelihood, presented in Upham et al. (2020), as this
approach is more computationally expensive. This is motivated by the fact that
with a sufficiently high number of simulations and for sufficiently high
multipoles, the Gaussian approximation is a satisfactory representation of the
underlying distribution for the cosmological parameters of interest (Taylor et
al., 2019; Lin et al., 2020; Upham et al., 2021).

The priors on all cosmological and nuisance parameters, with the exception of
the redshift displacement parameters, are chosen to be flat and are detailed
in Table 2. The following paragraphs motivate and explain all prior choices.
For more details, please refer to Appendix B from Heymans et al. (2021) and
Sect. 6.1 from Joachimi et al. (2021b).

As cosmic shear data are directly sensitive to the amplitude of density
fluctuations, $\sigma_{8}$, or the growth of structure, $S_{8}$, this prior
choice is of fundamental importance in our analysis (Hall, 2021). In Joachimi
et al. (2021b), it was shown that probing this parameter as a derived
parameter from the primordial power spectrum amplitude, $A_{\rm s}$ or
$\ln(10^{10}A_{\rm s})$, can lead to small biases as the prior becomes
informative for $S_{8}$. Therefore, in our analysis, we sample $S_{8}$
directly instead of treating it as a derived parameter. This guarantees that
our prior on this parameter is flat.

Other parameter priors that are not properly constrained by cosmic shear data
are chosen to be $\pm 5\sigma$ from independent measurements such as the prior
on $h$, taken from the Riess et al. (2016) measurements, and the prior on
$\omega_{\rm c}=\Omega_{\rm c}h^{2}$, derived from Supernova Type Ia
constraints on $\Omega_{\rm m}$ from Scolnic et al. (2018). Nevertheless,
there is one prior choice that has a significant impact when combining cosmic
shear data with external galaxy clustering data from SDSS (see Sect. 5.2 for
more details): the baryonic matter density, $\omega_{\rm b}=\Omega_{\rm
b}h^{2}$. This prior is chosen to reflect $\pm 5\sigma$ constraints from big
bang nucleosynthesis (BBN) presented by Olive & Particle Data Group (2014).
This is a conservative choice for cosmic shear-only analysis as justified by
Heymans et al. (2021) and Asgari et al. (2021). When combining cosmic shear
data with clustering information from BOSS/eBOSS galaxy clustering and the
eBOSS DR16 Ly$\alpha$ forest, this prior is a conservative equivalent to
adding BBN data to our analysis, as shown in Cuceu et al. (2019). Therefore,
we would like to stress that the combined constraints between KiDS-1000 weak
lensing and SDSS BAO and RSD shown in Sect. 6.1 should be interpreted with the
BBN-‘inspired’ priors in mind.

Finally, the redshift distributions for each of the five tomographic bins are
allowed to shift according to a correlated Gaussian prior with
$\bm{\mu}=(0.0001,0.0021,0.0129,0.0110,-0.0060)$ and a covariance
$\@tens{C}_{z}$ calibrated in Hildebrandt et al. (2021) using mock galaxy
surveys, according to the method presented by Wright et al. (2020).
Constraints for the nuisance parameters are presented in Fig. 13.

##  6 Results

###  6.1 Cosmological posteriors

We present our flat $\Lambda$CDM constraints for the pseudo-$C_{\ell}$
analysis of the KiDS-1000 cosmic shear catalogue, as well as combinations with
SDSS BAO and RSD (see Sect. 5.2), and comparisons with Planck 2018
TTTEEE+lowE202020Here we used the plik_lite_TTTEEE+lowl+lowE likelihood
installed within Monte Python 3. (Planck Collaboration et al., 2020a, b), the
band-powers constraints for KiDS-1000 obtained by Asgari et al. (2021) and the
$3\times 2$pt analysis from Heymans et al. (2021). A summary of the main
parameters constrained by cosmic shear data, $\sigma_{8}$, $\Omega_{\rm m}$,
and $S_{8}$, using MAP+PJ-HPD constraints are presented in Table 3 for
KiDS-1000 pseudo-$C_{\ell}$s and its combination with SDSS BAO and RDS data. A
more complete summary of results, including the traditional marginal
constraints, can be found in Table 5. Unless stated otherwise, the quoted
credible intervals in this section are the $68\%$ credible intervals (CI)
using the MAP+PJ-HPD constraints.

Table 3: Summary of the main parameters constrained by cosmic shear data and
combinations with galaxy clustering.

Parameter | KiDS-1000 P$C_{\ell}$ | KiDS-1000 P$C_{\ell}$ \+ SDSS
---|---|---
$S_{8}$ | $0.754_{-0.029}^{+0.027}$ | $0.771_{-0.032}^{+0.006}$
$\sigma_{8}$ | $0.920_{-0.215}^{+0.159}$ | $0.756_{-0.044}^{+0.029}$
$\Omega_{\rm m}$ | $0.202_{-0.053}^{+0.127}$ | $0.312_{-0.018}^{+0.011}$
212121Here, we are considering the 68% C.I. (1$\sigma$) MAP+PJ-PHD. It is
important to note that $\sigma_{8}$ and $\Omega_{\rm m}$ are not properly
constrained by cosmic shear alone.

Figure 8: Marginalised 2D posterior distributions for 68% (darker) and 95%
(lighter) C.I. in the $\Omega_{\rm m}-\sigma_{8}$ plane. Here we exhibit the
pseudo-$C_{\ell}$ constraints (red) in comparison with band-powers results
from Asgari et al. (2021) (orange). SDSS BAO and RSD (blue) and its
combination with the KiDS-1000 pseudo-$C_{\ell}$ measurements (purple) have a
completely different degeneracy than the Planck 2018 constraints from the
TTTEEE-lowE combination (green). We also compare it with the $3\times 2$pt
analysis from Heymans et al. (2021) (black-dashed). Figure 9: Same as Fig. 8,
but for the marginalised $\Omega_{\rm m}$ and $S_{8}$ constraints.

In Fig. 8, we display the known cosmic shear degeneracy between the amplitude
of matter fluctuations, $\sigma_{8}$, and the matter density parameters,
$\Omega_{\rm m}$. Results between the KiDS-1000 estimators are consistent,
with the pseudo-$C_{\ell}$ results (in red) demonstrating a slightly lower
value for $\sigma_{8}$ than the band-powers constraints (in orange) but
similar to other 2pt estimators presented in Asgari et al. (2021). No
degeneracy appears for the SDSS BAO and RSD analysis on this parameter plane,
meaning that its combination with cosmic shear considerably breaks the known
cosmic shear degeneracy on this plane. This is in line with the $3\times 2$pt
analysis performed by Heymans et al. (2021), shown as the black dashed
contours in Fig. 8, with the difference that the SDSS BAO and RSD data do not
assume flat $\Lambda$CDM, that is it does not use the full galaxy power
spectrum, yielding slightly larger constraints.

The clustering posteriors we obtain demonstrate a slightly higher value of
$\sigma_{8}$ than the BOSS constraints presented in Heymans et al. (2021),
which is more apparent when considering the $S_{8}=\sigma_{8}\sqrt{\Omega_{\rm
m}/0.3}$ parameter instead, as shown in Fig. 14. Here, we find
$S^{\text{SDSS}}_{8}=0.812_{-0.049}^{+0.050}$ for the SDSS BAO and RSD data
set which is in agreement with both KiDS-1000 pseudo-$C_{\ell}$,
$S_{8}^{\text{PCL}}=0.754_{-0.029}^{+0.027}$, and the Planck 2018 results. The
two main reasons why our clustering constraints differ from those in Heymans
et al. (2021) are the addition of high-redshift eBOSS LRGs and the Ly-$\alpha$
measurements as well as the impact of not using the full measured galaxy power
spectrum for the SDSS measurements, like in Tröster et al. (2020). Yet, when
combining the clustering information with cosmic shear, we can see from Fig. 9
that the difference is mostly in the $\Omega_{\rm m}$ parameter. This
demonstrates the $S_{8}$ measurement robustness from KiDS-1000 cosmic shear
analyses. The $S_{8}$ constraints with the pseudo-$C_{\ell}$ estimator are
slightly larger than the one found with band-powers from Asgari et al. (2021),
reflecting a $3.7\%$ increase in comparison. This is in line with the expected
increase due to our noisy covariance matrix estimation method, which applies
the estimator to simulations instead of using a theoretical covariance.

A summary of $S_{8}$ measurements with comparisons to other cosmic shear and
CMB probes is displayed in Fig. 10. Overall, weak lensing experiments are
consistent between each other. The purple vertical shaded region in Fig. 10
reflects our main constraints from KiDS-1000 pseudo-$C_{\ell}$ combined with
BAO and RSD from BOSS and eBOSS LRGs, and BAO from Ly$\alpha$ forest and its
cross-correlations with quasars (SDSS BAO and RSD). Results from the Hyper
Supreme Camera (HSC) first year data (Hikage et al., 2019), Dark Energy Survey
Year 1 (DES-Y1) analysis (Troxel et al., 2018), and the CFHTLenS constraints
from Joudaki et al. (2017)222222Here, we quote the mid nominal values as this
analysis is closer to the analysis performed for KV-450 and KiDS-1000. are all
in agreement with our analysis. Consistency between KiDS-1000 and previous
KiDS analyses, such as KV450 (Hildebrandt et al., 2020) and KiDS-450
(Hildebrandt et al., 2017), was explored in detail by Asgari et al. (2021). A
lower value for $S_{8}$, when compared to the CMB constraints, is consistently
found across a range of different and independent cosmic shear experiments
carried out by various research groups.

When comparing the amplitude of intrinsic alignments between the
pseudo-$C_{\ell}$ and band-power constraints, we find lower values in the
pseudo-$C_{\ell}$ analysis, in line with what has been found for the other 2pt
statistics in Asgari et al. (2021), COSEBIs and correlation functions:

[PCL] $\displaystyle A_{\text{IA}}=0.396_{-0.931}^{+0.251}\,;$ (26) [2PCF]
$\displaystyle A_{\text{IA}}=0.387_{-0.374}^{+0.321}\,;$ (27) [COSEBIs]
$\displaystyle A_{\text{IA}}=0.264_{-0.337}^{+0.424}\,;$ (28) [BP]
$\displaystyle A_{\text{IA}}=0.973_{-0.383}^{+0.292}\,.$ (29)

We speculate that this could be the reason why our method finds a slightly
lower peak for the matter density parameter, $\Omega_{\rm m}$, as a similar
effect appears in Asgari et al. (2021). Despite this, $\Omega_{\rm m}$ is not
strongly constrained in our analysis, meaning that the uncertainty is large
enough that the constraints remain consistent for both $\Omega_{\rm m}$ and
$A_{\text{IA}}$, as investigated by Asgari et al. (2021). Since both band-
powers and pseudo-$C_{\ell}$s are in harmonic space, one could expect them to
match better than pseudo-$C_{\ell}$s and correlation functions, but the
difference between both measurements is $\sim 1.2\sigma$ away (see Eq. 30). We
note that in contrast to the other 2pt functions presented in Asgari et al.
(2021), band-powers filter out most of the scales below $\ell=100$. Meanwhile,
correlation functions and COSEBIs have filters in harmonic space that capture
power below $\ell=100$. Comparing this with the pseudo-$C_{\ell}$s
$\ell_{\text{min}}=76$, Fig. 3 shows that larger scales also leak through the
mixing matrix filters into our analysis. This slight difference in $\ell$
range potentially explains the differences between band-powers and
pseudo-$C_{\ell}$, and their similarity with the correlation function and
COSEBIs results since intrinsic alignments have a significant impact in larger
scales as seen in Fig. 4 from Asgari et al. (2021).

The full posteriors for cosmological parameters and the lensing nuisance
parameters, $A_{\rm bary}$ and $A_{\text{IA}}$, are shown in Fig. 11, and the
redshift nuisance parameters are shown in the Appendix in Fig. 13 with a
comparison to the band-power constraints from Asgari et al. (2021). Table 5
presents a summary of our cosmological inference results, including the
traditional marginalised constraints. For our best-fit parameters, the reduced
$\chi^{2}$ for the cosmic shear-only analysis is 1.18, reflecting an
acceptable p-value of 0.06.

###  6.2 Cosmic shear tension with early Universe probes

Figure 10: Comparison between different cosmological probe constraints on the
growth of structure parameter, $S_{8}=\sigma_{8}(\Omega_{\rm m}/0.3)^{1/2}$.
Here, we present our results with both the fiducial MAP+PJ-HPD constraints
(solid) and the traditional marginal posterior mode (dashed). The purple
shaded region reflects the MAP+PJ-HPD constraints for KiDS-1000
pseudo-$C_{\ell}$ plus SDSS BAO and RSD. For external probes quoted from the
literature, we quote the marginal posterior mode with tail credible intervals
(dotted). Our analysis is not only consistent with other KiDS-1000 analysis
but also consistent with other independent cosmic shear probes. Tensions
between cosmic shear experiments and Planck 2018 (green shaded region) or
ACT+Planck are between 2.8 to 3.2$\sigma$.  Figure 11: Extended marginalised
1D and 2D posteriors for relevant $\Lambda$CDM and astrophysical parameters:
the matter density parameter ($\Omega_{\rm m}$), the amplitude of matter
density fluctuations ($S_{8}=\sigma_{8}\sqrt{\Omega_{\rm m}/0.3}$),
$\sigma_{8}$ as a derived parameter, the physical baryonic matter density
($\Omega_{\rm b}h^{2}$), the dimensionless Hubble constant ($h$), and the
spectral index ($n_{\rm s}$). Included are also two nuisance parameters: the
amplitude of intrinsic alignments ($A_{\rm IA}$) and the baryonic feedback
amplitude ($A_{\rm bary}$) as these correlate with some cosmological
parameters. When combining the KiDS-1000 pseudo-$C_{\ell}$ constraints (red)
with the SDSS BAO and RSD constraints (blue), we obtain constraints (purple)
that break a few degeneracies that cosmic shear alone cannot. For comparison,
we added the Planck 2018 TTTEEE-lowE (Planck Collaboration et al., 2020b)
constrains (green) and the KiDS-1000 band-powers constraints (orange) from
Asgari et al. (2021).

We now turn our discussion towards the known inconsistencies between cosmic
shear surveys and the cosmic microwave background in the amplitude of matter
density fluctuations and, consequently, the growth of structure as measured by
$S_{8}$. By assuming Gaussianity and adopting the marginal constraints to
quantify tension, we initially employ the conventional marginal tension
estimate

$\displaystyle\tau=\frac{|\langle S_{8}^{\rm A}\rangle-\langle S_{8}^{\rm
B}\rangle|}{\sqrt{\sigma^{2}(S_{8}^{\rm A})+\sigma^{2}(S_{8}^{\rm B})}}\,,$
(30)

where $\langle S_{8}^{\rm A}\rangle$ and $\sigma^{2}(S_{8}^{\rm A})$ are the
mean and variance for a given probe A. Under the Gaussianity assumption,
$\tau$ measures how consistent the difference between $\langle S_{8}^{\rm
A}\rangle$ and $\langle S_{8}^{\rm B}\rangle$ is with a case where there is no
difference between measurements.

Comparisons between KiDS-1000 pseudo-$C_{\ell}$ measurements and Planck 2018
TTTEEE+lowE nominal results (Planck Collaboration et al., 2020a, b) lead to a
$2.8\sigma$ tension between late- and early-Universe constraints. When
combining the pseudo-$C_{\ell}$ measurements with SDSS BAO and RSD data, the
tension with Planck 2018 increases to 2.9$\sigma$, in line with the tension
found in the KiDS-1000 $3\times 2$pt analysis from Heymans et al. (2021).

To further assess the tension between late-time and early-time probes, we
verify the tension between our KiDS-1000 pseudo-$C_{\ell}$ constraints and the
most recent results from the Atacama Cosmology Telescope (ACT), a ground-based
CMB experiment (Aiola et al., 2020). We compare our results with ACT alone,
ACT combined with WMAP, and ACT combined with Planck 2018. These nominal
constraints can also be seen in Fig. 10.

As expected, tensions with ACT constraints alone are much smaller given that
the data set is not as precise as space-based CMB observations; we find
$\tau=1.4\sigma$ when combining cosmic shear and clustering data. However, ACT
does not measure the TT power spectrum below $\ell=600$, and the TE and EE
power spectra below multipole 350, losing the first two acoustic peaks in the
temperature auto spectrum and the entirety of the first peak in the TE/EE
spectra. As a consequence, ACT’s nominal constraints for $S_{8}$ include
information from the final WMAP public data release (Bennett et al., 2013).
The tension between our constraints from cosmic shear and clustering increases
to $2.2\sigma$ for the nominal ACT plus WMAP results.

The largest tension between cosmic microwave background and large-scale
structure measurements appears when combining ACT and Planck. For
pseudo-$C_{\ell}$ alone, we find $3.2\sigma$ tension and combining cosmic
shear with clustering, we find $3.3\sigma$ tension. However, using the
suspiciousness statistic (Lemos et al., 2020), Handley & Lemos (2021) pointed
out that Planck 2018 and ACT have a $2.6\sigma$ tension between each other.

It is clear from Fig. 11 and Table 3, that the marginal distribution of
$S_{8}$ is not Gaussian. To further assess the tension between our KiDS-1000
pseudo-$C_{\ell}$ measurements and the Planck Legacy results, we also adopt
the Hellinger distance measure, $d_{\text{H}}$, as a tension metric (Beran
1977, and more details in Appendix F from Heymans et al. 2021). The Hellinger
distance is a stable metric for the comparison of arbitrary one-dimensional
distributions. For two distinct distributions, $p(\theta)$ and $q(\theta$),
$d_{\text{H}}$ is defined as

$\displaystyle d_{\text{H}}^{2}[p,q]=\frac{1}{2}\int{\rm
d}\theta\left[\sqrt{p(\theta)}-\sqrt{q(\theta})\right]^{2}\,.$ (31)

With the metric above, we measured $d_{\text{H}}=0.95$ between our KiDS-1000
pseudo-$C_{\ell}$ and Planck Legacy’s $S_{8}$ constraints. Using the
methodology described by Heymans et al. (2021), this can be expressed as a
$\sim 3.1\sigma$ tension. For our combination of cosmic shear from KiDS and
clustering from SDSS, we obtain $d_{\text{H}}=0.94$ – corresponding to a $\sim
3\sigma$ tension with Planck. Since this tension metric uses the marginalised
distribution for a given parameters, one needs the MCMC chains to calculate
it. The ACT DR4 chains were not publicly available at the time. Therefore, we
were not able to verify $d_{\text{H}}$ between KiDS-1000, ACT and combinations
with ACT.

##  7 Conclusions

We performed a pseudo angular power spectra (pseudo-$C_{\ell}$) analysis of
the current state-of-the-art galaxy shape measurements from the Kilo Degree
Survey’s fourth data release (KiDS-1000, Kuijken et al. 2019), obtaining
cosmological constraints for the current concordance model, the flat
$\Lambda$CDM cosmology. The KiDS-1000 shear catalogue (Giblin et al., 2021) is
now publicly
available232323http://kids.strw.leidenuniv.nl/DR4/KiDS-1000_shearcatalogue.php
https://github.com/KiDS-WL/Cat_to_Obs_K1000_P1, containing nine-band
photometry with shape and photometric redshift estimates for more than $21$
million lensed objects distributed over an area of around 1000 deg2.

With this work, we added an additional 2pt function to the family of
measurements used in Asgari et al. (2021). With its highly accurate calibrated
shear and photometric redshift measurements, the KiDS-1000 shear catalogue is
an excellent proxy for future surveys, including their early data-releases.
The performance of our pseudo-$C_{\ell}$ estimator on KiDS-1000 data marks an
important road-test for its future application on data from Euclid cosmic
shear catalogue. Even for a comparatively small survey area and a disjointed
footprint, the pseudo-$C_{\ell}$ estimator performed as well as other 2pt
functions, with no flat-sky approximation — an important and necessary
characteristic that estimators of the next generation of galaxy surveys need.
The forward-modelling approach allowed us to maintain accuracy and avoid
numerical instabilities from inverting and deconvolving a mixing matrix from
our estimates with no significant computational cost when calculating our
theory-vector for cosmological inference.

Through this re-analysis in harmonic space of KiDS-1000 data as a road-test of
our prototype estimator for the Euclid Science Ground Segment we find that,
with KiDS data, we are operating at the limit of this methodology. First,
noise realisations of the galaxy shape catalogue will not be feasible for the
size of the expected Euclid shape catalogue ($\sim 2\times 10^{9}$ galaxies);
hence we will implement a noise power spectra subtraction based on the
theoretical expectation of the noise. Second, the treatment used for the shear
weights in the analysis is sufficient for KiDS, but it may not be the case for
Euclid. For higher resolution shear maps such as those required by the Euclid
Science Ground Segment ($N_{\rm side}=4096$), the shear weights will need to
be incorporated into the mask, requiring our formalism to change. Finally, we
will also investigate the possibility of deconvolving the effects of the
mixing matrix while maintaining the level of accuracy and precision required
by Euclid outlined in Laureijs et al. (2011), which does not allow for mask
apodisation. We have learned valuable lessons for Euclid which we will explore
in a subsequent paper (Euclid Consortium et al., in prep).

Following, our systematic contamination null-test analysis also used the
pseudo-$C_{\ell}$ estimator at its core, demonstrating the versatility of this
method. Using 2pt statistics in harmonic space to assess systematic
contamination yields insights regarding the validity of cosmological
information on the scales used to infer cosmology. This approach allows us to
be confident when using a broad range of multipoles from measured angular
power spectra as any spurious correlations arising solely from systematic
contamination will be picked-up by it. We find no significant correlations
between systematics and data for any relevant observational catalogue
quantities nor for stellar density from Gaia DR2. These null detection
analyses, together with the null detection of B-modes in our data, made us
confident that the measured E-mode pseudo-$C_{\ell}$s contained only
cosmological and astrophysical signals and could be used for Bayesian
inference of the relevant flat $\Lambda$CDM parameters.

Cosmic shear alone constrains the growth of structure parameter, $S_{8}$,
quite well but it has a known degeneracy in the amplitude of matter density
fluctuations and density of matter plane of the parameter space. To break this
degeneracy, we combined our pseudo-$C_{\ell}$ measurements with baryon
acoustic oscillations and redshift space distortion measurements from luminous
red galaxies from the Sloan Digital Sky Survey’s BOSS DR12 and eBOSS DR16 data
sets, as well as baryon acoustic oscillations from the Lyman-$\alpha$ forest
and its cross-correlations with quasar positions from eBOSS DR16 (Alam et al.,
2021). The addition of clustering data allowed us to break the aforementioned
degeneracy.

Our cosmological inference takes a few cosmic shear and redshift distribution-
related nuisance parameters into account. From these, we can observe an
interesting trend for the amplitude of intrinsic alignments, $A_{\text{IA}}$,
when compared to previous results from Asgari et al. (2021). In contrast to
band-power measurements, pseudo-$C_{\ell}$, as well as correlation functions
and COSEBIs, found a lower value of $A_{\text{IA}}$. We believe that these
results are related to the larger scales accessible by the different analyses.
Band-powers have a higher maximum scale cut at $\ell_{\text{min}}=100$, while
we used $\ell_{\text{min}}=76$ for pseudo-$C_{\ell}$, for example. We
emphasise that, given the uncertainty in these measurements, there is no
evidence for inconsistencies from $A_{\text{IA}}$ between the different 2pt
functions as measurements agree up to $\sim 1.2\sigma$ (using the same metrics
as in Eq. 30).

Quoting the maximum a posteriori (MAP) value with the projected joint highest
posterior density (PJ-HPD) credible intervals (Joachimi et al., 2021b), our
pseudo-$C_{\ell}$ analysis alone yielded $S_{8}=0.754_{-0.029}^{+0.027}$, or
$9.8\pm 3.3\%$ lower than the $S_{8}$ nominal constraints from Planck
Collaboration et al. (2020a). When adding clustering data from SDSS, we
increased the precision on the growth of structure parameter, obtaining a
measurement of $S_{8}=0.771^{+0.006}_{-0.032}$, reflecting a $7.8\pm 2.3\%$
lower value when comparing the marginal constraints with Planck’s measurement
of $S_{8}$. Both constraints are in very good agreement with other KiDS-1000
measurements from Asgari et al. (2021); Heymans et al. (2021) and other cosmic
shear surveys such as the Dark Energy Survey (Troxel et al., 2018), the Hyper
Supreme Camera survey (Hikage et al., 2019) and CFHTLenS (Joudaki et al.,
2017).

Using two different metrics to evaluate tension in a marginalised space, we
assessed the tension in the growth of structure parameter between cosmic shear
and probes of the cosmic microwave background. Tension with the Planck Legacy
TTTEEElowE analysis is around 2.8$\sigma$ for cosmic shear alone, in line with
what was found by Asgari et al. (2021). When comparing the Planck Legacy
results to our combination of cosmic shear with clustering data from SDSS, the
tension increased to $\sim 2.9\sigma$, also in line with findings from Heymans
et al. (2021).

We further verify that the tensions we observed did not arise exclusively from
comparisons with the Planck Mission results and there is indeed a tension
between the large-scale structure of the late Universe and the cosmic
microwave background radiation from the early Universe. We analysed this by
comparing our results to the nominal constraints from a ground-based CMB
experiment, the Atacama Cosmology Telescope Data Release 4 (ACT, Aiola et al.
2020). The nominal results from ACT also included data from nine years of
observations from WMAP (Bennett et al., 2013); the tension between our
constraints from cosmic shear combined with clustering and ACT+WMAP were not
as high as with Planck Legacy: $2.2\sigma$. However, ACT was also combined
with Planck instead of WMAP. In this case, tensions between large-scale
structure probes and the cosmic microwave background were as high as
$3.3\sigma$. The impact of Planck on this tension value, however apparent,
cannot be properly assessed as evidence points to internal tensions between
ACT, WMAP, and Planck (Handley & Lemos, 2021). The tensions between the late
and early Universe seem to be independent of Planck Legacy data, but are
exacerbated by it.

So far, the tension has remained at a level that is inconclusive as to whether
it is arising from unknown systematics in the data, noise in the realisation
of the portions of the Universe that we are observing, or hints of new
physics. The Kilo Degree Survey still has a final future data release (DR5)
and our hopes are that, with it, we will obtain a better understanding of what
is causing this tension with the growth of structure in the Universe.
Nonetheless, the next decade will see the light of new and more powerful
cosmic shear surveys with the launch of the Euclid Space Telescope and the
start of the survey carried out by the Vera C. Rubin Observatory. Both
surveys, independently and combined, will produce an almost complete sky map
of the large-scale structure of the Universe, using cosmic shear as one of
their main probes of understanding cosmology, and with the potential of
finally solving this cosmic puzzle.

###### Acknowledgements.

The authors would like to thank Dr Lorne Whiteway for useful comments and
discussions. We would also like to thank the reviewer for their helpful
comments that helped enrich this work. The mixing matrix calculation used a
publicly available Wigner 3-j symbol code developed by Dr Whiteway and it can
be found in https://github.com/LorneWhiteway/UCLWig3j. The K-Means code used
for the jackknife covariances for systematics was built on top of a code
developed by Dr Erin Sheldon, publicly available in
https://github.com/esheldon/kmeans_radec. The analysis performed in this work
also made use of the following packages and software: NumPy (Harris et al.,
2020), SciPy (Virtanen et al., 2020), Matplotlib (Hunter, 2007), GetDist
(Lewis, 2019), Pandas (Wes McKinney, 2010), AstroPy (Astropy Collaboration et
al., 2018) and Jupyter Lab (Kluyver et al., 2016). The Flask + Salmo mocks
were generated and provided by Dr Chieh-An Lin. We thank Dr Edd Edmonson for
technical support at the UCL HPC facilities. This work used computing
equipment funded by the Research Capital Investment Fund (RCIF) provided by
UKRI, and is partially funded by the UCL Cosmoparticle Initiative. MA and TT
acknowledge support from the European Research Council under grant number
647112. AHW and AD acknowledge support from the European Research Council
Consolidator Grant (No. 770935). MB is supported by the Polish National
Science Center through grants no. 2020/38/E/ST9/00395, 2018/30/E/ST9/00698 and
2018/31/G/ST9/03388, and by the Polish Ministry of Science and Higher
Education through grant DIR/WK/2018/12. BG acknowledges support from the
European Research Council under grant number 647112 and from the Royal Society
through an Enhancement Award (RGF/EA/181006). CH acknowledges support from the
European Research Council under grant number 647112, and support from the Max
Planck Society and the Alexander von Humboldt Foundation in the framework of
the Max Planck-Humboldt Research Award endowed by the Federal Ministry of
Education and Research. H. Hildebrandt is supported by a Heisenberg grant of
the Deutsche Forschungsgemeinschaft (Hi 1495/5-1) as well as an ERC
Consolidator Grant (No. 770935). HYS acknowledges the support from NSFC of
China under grant 11973070, the Shanghai Committee of Science and Technology
grant No.19ZR1466600 and Key Research Program of Frontier Sciences, CAS, Grant
No. ZDBS-LY-7013.
The KiDS-1000 data products in this paper are based on observations made with
ESO Telescopes at the La Silla Paranal Observatory under programme IDs
177.A-3016, 177.A-3017 and 177.A-3018, and on data products produced by
Target/OmegaCEN, INAF-OACN, INAF-OAPD and the KiDS production team, on behalf
of the KiDS consortium.
The authors acknowledge the Euclid Consortium, the European Space Agency, and
a number of agencies and institutes that have supported the development of
Euclid, in particular the Academy of Finland, the Agenzia Spaziale Italiana,
the Belgian Science Policy, the Canadian Euclid Consortium, the French Centre
National d’Etudes Spatiales, the Deutsches Zentrum für Luft- und Raumfahrt,
the Danish Space Research Institute, the Fundação para a Ciência e a
Tecnologia, the Ministerio de Economia y Competitividad, the National
Aeronautics and Space Administration, the National Astronomical Observatory of
Japan, the Netherlandse Onderzoekschool Voor Astronomie, the Norwegian Space
Agency, the Romanian Space Agency, the State Secretariat for Education,
Research and Innovation (SERI) at the Swiss Space Office (SSO), and the United
Kingdom Space Agency. A complete and detailed list is available on the Euclid
web site (http://www.euclid-ec.org).
Author Contributions: All authors contributed to the development and writing
of this paper. The authorship list is given in several groups: the lead
authors (AL, LW, ASM, BJ, AC), followed by an alphabetical group that includes
those who are key contributors to both the scientific analysis and the KiDS
data products, and a further alphabetical group covering those who have either
made a significant contribution to the KiDS data products or to the scientific
analysis. Additional alphabetical groups correspond to different authorship
levels of the Euclid Consortium.

## References

  * Abbott et al. (2018) Abbott, T. M. C., Abdalla, F. B., Alarcon, A., et al. 2018, Phys. Rev. D, 98, 043526
  * Aiola et al. (2020) Aiola, S., Calabrese, E., Maurin, L., et al. 2020, J. Cosmology Astropart. Phys., 2020, 047
  * Alam et al. (2017) Alam, S., Ata, M., Bailey, S., et al. 2017, MNRAS, 470, 2617
  * Alam et al. (2021) Alam, S., Aubert, M., Avila, S., et al. 2021, Phys. Rev. D, 103, 083533
  * Alonso et al. (2019) Alonso, D., Sanchez, J., Slosar, A., & LSST Dark Energy Science Collaboration. 2019, MNRAS, 484, 4127
  * Amon et al. (2021) Amon, A., Gruen, D., Troxel, M. A., et al. 2021, arXiv e-prints, arXiv:2105.13543
  * Asgari et al. (2021) Asgari, M., Lin, C.-A., Joachimi, B., et al. 2021, A&A, 645, A104
  * Asgari et al. (2018) Asgari, M., Taylor, A., Joachimi, B., & Kitching, T. D. 2018, MNRAS, 479, 454
  * Asgari et al. (2020) Asgari, M., Tröster, T., Heymans, C., et al. 2020, A&A, 634, A127
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123
  * Audren et al. (2013) Audren, B., Lesgourgues, J., Benabed, K., & Prunet, S. 2013, J. Cosmology Astropart. Phys., 1302, 001
  * Bartelmann & Schneider (2001) Bartelmann, M. & Schneider, P. 2001, Phys. Rep, 340, 291
  * Bautista et al. (2021) Bautista, J. E., Paviot, R., Vargas Magaña, M., et al. 2021, MNRAS, 500, 736
  * Becker & Rozo (2016) Becker, M. R. & Rozo, E. 2016, MNRAS, 457, 304
  * Becker et al. (2016) Becker, M. R., Troxel, M. A., MacCrann, N., et al. 2016, Phys. Rev. D, 94, 022002
  * Begeman et al. (2013) Begeman, K., Belikov, A. N., Boxhoorn, D. R., & Valentijn, E. A. 2013, Experimental Astronomy, 35, 1
  * Benítez (2000) Benítez, N. 2000, ApJ, 536, 571
  * Bennett et al. (2013) Bennett, C. L., Larson, D., Weiland, J. L., et al. 2013, ApJS, 208, 20
  * Beran (1977) Beran, R. 1977, Ann. Statist., 5, 445
  * Bernal et al. (2016) Bernal, J. L., Verde, L., & Riess, A. G. 2016, J. Cosmology Astropart. Phys., 2016, 019
  * Bernstein (2009) Bernstein, G. M. 2009, ApJ, 695, 652
  * Blas et al. (2011) Blas, D., Lesgourgues, J., & Tram, T. 2011, J. Cosmology Astropart. Phys., 2011, 034
  * Bridle & King (2007) Bridle, S. & King, L. 2007, New Journal of Physics, 9, 444
  * Bridle et al. (2002) Bridle, S. L., Crittenden, R., Melchiorri, A., et al. 2002, MNRAS, 335, 1193
  * Brinckmann & Lesgourgues (2019) Brinckmann, T. & Lesgourgues, J. 2019, Physics of the Dark Universe, 24, 100260
  * Brown et al. (2005) Brown, M. L., Castro, P. G., & Taylor, A. N. 2005, MNRAS, 360, 1262
  * Buchner et al. (2014) Buchner, J., Georgakakis, A., Nandra, K., et al. 2014, A&A, 564, A125
  * Coe et al. (2006) Coe, D., Benítez, N., Sánchez, S. F., et al. 2006, AJ, 132, 926
  * Cuceu et al. (2019) Cuceu, A., Farr, J., Lemos, P., & Font-Ribera, A. 2019, J. Cosmology Astropart. Phys., 2019, 044
  * de Jong et al. (2015) de Jong, J. T. A., Verdoes Kleijn, G. A., Boxhoorn, D. R., et al. 2015, A&A, 582, A62
  * de Jong et al. (2017) de Jong, J. T. A., Verdoes Kleijn, G. A., Erben, T., et al. 2017, A&A, 604, A134
  * DES Collaboration et al. (2022) DES Collaboration, Abbott, T. M. C., Aguena, M., et al. 2022, Phys. Rev. D, 105, 023520
  * Di Valentino et al. (2016) Di Valentino, E., Melchiorri, A., & Silk, J. 2016, Phys. Lett. B, 761, 242
  * Doux et al. (2021) Doux, C., Chang, C., Jain, B., et al. 2021, MNRAS, 503, 3796
  * du Mas des Bourboux et al. (2017) du Mas des Bourboux, H., Le Goff, J.-M., Blomqvist, M., et al. 2017, A&A, 608, A130
  * du Mas des Bourboux et al. (2020) du Mas des Bourboux, H., Rich, J., Font-Ribera, A., et al. 2020, ApJ, 901, 153
  * Edge et al. (2013) Edge, A., Sutherland, W., Kuijken, K., et al. 2013, The Messenger, 154, 32
  * Efstathiou (2004) Efstathiou, G. 2004, MNRAS, 349, 603
  * Efstathiou (2020) Efstathiou, G. 2020, arXiv e-prints, arXiv:2007.10716
  * Efstathiou & Gratton (2019) Efstathiou, G. & Gratton, S. 2019, arXiv e-prints, arXiv:1910.00483
  * Efstathiou & Gratton (2020) Efstathiou, G. & Gratton, S. 2020, MNRAS, 496, L91
  * Erben et al. (2005) Erben, T., Schirmer, M., Dietrich, J. P., et al. 2005, Astronomische Nachrichten, 326, 432
  * Euclid Collaboration et al. (2020) Euclid Collaboration, Blanchard, A., Camera, S., et al. 2020, A&A, 642, A191
  * Feroz & Hobson (2008) Feroz, F. & Hobson, M. P. 2008, MNRAS, 384, 449
  * Feroz et al. (2009) Feroz, F., Hobson, M. P., & Bridges, M. 2009, MNRAS, 398, 1601
  * Feroz et al. (2019) Feroz, F., Hobson, M. P., Cameron, E., & Pettitt, A. N. 2019, The Open Journal of Astrophysics, 2, 10
  * Fortuna et al. (2021) Fortuna, M. C., Hoekstra, H., Joachimi, B., et al. 2021, MNRAS, 501, 2983
  * Gaia Collaboration et al. (2018) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, A&A, 616, A1
  * Gallozzi et al. (2020a) Gallozzi, S., Paris, D., Scardia, M., & Dubois, D. 2020a, arXiv e-prints, arXiv:2003.05472
  * Gallozzi et al. (2020b) Gallozzi, S., Scardia, M., & Maris, M. 2020b, arXiv e-prints, arXiv:2001.10952
  * García-García et al. (2019) García-García, C., Alonso, D., & Bellini, E. 2019, J. Cosmology Astropart. Phys., 2019, 043
  * García-García et al. (2021) García-García, C., Ruiz-Zapatero, J., Alonso, D., et al. 2021, J. Cosmology Astropart. Phys., 2021, 030
  * Giblin et al. (2021) Giblin, B., Heymans, C., Asgari, M., et al. 2021, A&A, 645, A105
  * Gil-Marín et al. (2020) Gil-Marín, H., Bautista, J. E., Paviot, R., et al. 2020, MNRAS, 498, 2492
  * Górski et al. (2005) Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759
  * Hall (2021) Hall, A. 2021, MNRAS, 505, 4935
  * Han et al. (2021) Han, D., Sehgal, N., MacInnis, A., et al. 2021, J. Cosmology Astropart. Phys., 2021, 031
  * Handley & Lemos (2021) Handley, W. & Lemos, P. 2021, Phys. Rev. D, 103, 063529
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357
  * Heydenreich et al. (2020) Heydenreich, S., Schneider, P., Hildebrandt, H., et al. 2020, A&A, 634, A104
  * Heymans et al. (2013) Heymans, C., Grocutt, E., Heavens, A., et al. 2013, MNRAS, 432, 2433
  * Heymans et al. (2021) Heymans, C., Tröster, T., Asgari, M., et al. 2021, A&A, 646, A140
  * Heymans et al. (2012) Heymans, C., Van Waerbeke, L., Miller, L., et al. 2012, MNRAS, 427, 146
  * Hikage et al. (2019) Hikage, C., Oguri, M., Hamana, T., et al. 2019, PASJ, 71, 43
  * Hikage et al. (2011) Hikage, C., Takada, M., Hamana, T., & Spergel, D. 2011, MNRAS, 412, 65
  * Hildebrandt et al. (2020) Hildebrandt, H., Köhlinger, F., van den Busch, J. L., et al. 2020, A&A, 633, A69
  * Hildebrandt et al. (2021) Hildebrandt, H., van den Busch, J. L., Wright, A. H., et al. 2021, A&A, 647, A124
  * Hildebrandt et al. (2017) Hildebrandt, H., Viola, M., Heymans, C., et al. 2017, MNRAS, 465, 1454
  * Hivon et al. (2002) Hivon, E., Górski, K. M., Netterfield, C. B., et al. 2002, ApJ, 567, 2
  * Hotelling (1931) Hotelling, H. 1931, The Annals of Mathematical Statistics, 2, 360
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90
  * Joachimi & Bridle (2010) Joachimi, B. & Bridle, S. L. 2010, A&A, 523, A1
  * Joachimi et al. (2021a) Joachimi, B., Köhlinger, F., Handley, W., & Lemos, P. 2021a, A&A, 647, L5
  * Joachimi et al. (2021b) Joachimi, B., Lin, C. A., Asgari, M., et al. 2021b, A&A, 646, A129
  * Joachimi et al. (2011) Joachimi, B., Mandelbaum, R., et al. 2011, A&A, 527, A26
  * Joachimi et al. (2008) Joachimi, B., Schneider, P., & Eifler, T. 2008, A&A, 477, 43
  * Joudaki et al. (2017) Joudaki, S., Blake, C., Heymans, C., et al. 2017, MNRAS, 465, 2033
  * Joudaki et al. (2018) Joudaki, S., Blake, C., Johnson, A., et al. 2018, MNRAS, 474, 4894
  * Joudaki et al. (2020) Joudaki, S., Hildebrandt, H., Traykova, D., et al. 2020, A&A, 638, L1
  * Kaiser (1992) Kaiser, N. 1992, ApJ, 388, 272
  * Kannawadi et al. (2019) Kannawadi, A., Hoekstra, H., Miller, L., et al. 2019, A&A, 624, A92
  * Kilbinger (2015) Kilbinger, M. 2015, Reports on Progress in Physics, 78, 086901
  * Kilbinger et al. (2017) Kilbinger, M., Heymans, C., Asgari, M., et al. 2017, MNRAS, 472, 2126
  * Kitching et al. (2017) Kitching, T. D., Alsing, J., Heavens, A. F., et al. 2017, MNRAS, 469, 2737
  * Kitching et al. (2020) Kitching, T. D., Deshpande, A. C., & Taylor, P. L. 2020, The Open Journal of Astrophysics, 3, 14
  * Kitching et al. (2008) Kitching, T. D., Miller, L., Heymans, C. E., van Waerbeke, L., & Heavens, A. F. 2008, MNRAS, 390, 149
  * Kitching et al. (2019a) Kitching, T. D., Paykari, P., Hoekstra, H., & Cropper, M. 2019a, The Open Journal of Astrophysics, 2, 5
  * Kitching et al. (2019b) Kitching, T. D., Taylor, P. L., Capak, P., Masters, D., & Hoekstra, H. 2019b, arXiv e-prints, arXiv:1901.06495
  * Kluyver et al. (2016) Kluyver, T., Ragan-Kelley, B., Pérez, F., et al. 2016, in Positioning and Power in Academic Publishing: Players, Agents and Agendas, ed. F. Loizides & B. Scmidt (Netherlands: IOS Press), 87–90
  * Kogut et al. (2003) Kogut, A., Spergel, D. N., Barnes, C., et al. 2003, ApJS, 148, 161
  * Krause & Eifler (2017) Krause, E. & Eifler, T. 2017, MNRAS, 470, 2100
  * Krause et al. (2021) Krause, E., Fang, X., Pandey, S., et al. 2021, arXiv e-prints, arXiv:2105.13548
  * Kuijken et al. (2019) Kuijken, K., Heymans, C., Dvornik, A., et al. 2019, A&A, 625, A2
  * Kuijken et al. (2015) Kuijken, K., Heymans, C., Hildebrandt, H., et al. 2015, MNRAS, 454, 3500
  * Kwan et al. (2017) Kwan, J., Sánchez, C., Clampitt, J., et al. 2017, MNRAS, 464, 4045
  * Laureijs et al. (2011) Laureijs, R., Amiaux, J., Arduini, S., et al. 2011, arXiv e-prints, arXiv:1110.3193
  * Leistedt & Peiris (2014) Leistedt, B. & Peiris, H. V. 2014, MNRAS, 444, 2
  * Lemos et al. (2017) Lemos, P., Challinor, A., & Efstathiou, G. 2017, J. Cosmology Astropart. Phys., 2017, 014
  * Lemos et al. (2020) Lemos, P., Köhlinger, F., Handley, W., et al. 2020, MNRAS, 496, 4647
  * Lemos et al. (2021) Lemos, P., Raveri, M., Campos, A., et al. 2021, MNRAS, 505, 6179
  * Lesgourgues (2011) Lesgourgues, J. 2011, arXiv e-prints, arXiv:1104.2932
  * Lewis (2019) Lewis, A. 2019, arXiv e-prints, arXiv:1910.13970
  * Li et al. (2019) Li, Y., Singh, S., Yu, B., Feng, Y., & Seljak, U. 2019, J. Cosmology Astropart. Phys., 2019, 016
  * Lin et al. (2020) Lin, C.-H., Harnois-Déraps, J., Eifler, T., et al. 2020, MNRAS, 499, 2977
  * Lin & Ishak (2017) Lin, W. & Ishak, M. 2017, Phys. Rev. D, 96, 083532
  * Loureiro et al. (2019) Loureiro, A., Cuceu, A., et al. 2019, Phys. Rev. Lett., 123, 081301
  * LoVerde & Afshordi (2008) LoVerde, M. & Afshordi, N. 2008, Phys. Rev. D, 78, 123506
  * LSST Dark Energy Science Collaboration (2012) LSST Dark Energy Science Collaboration. 2012, arXiv e-prints, arXiv:1211.0310
  * MacCrann et al. (2015) MacCrann, N., Zuntz, J., Bridle, S., Jain, B., & Becker, M. R. 2015, Mon. Not. Roy. Astron. Soc., 451, 2877
  * Maehoenen & Hakala (1995) Maehoenen, P. H. & Hakala, P. J. 1995, ApJ, 452, L77
  * Masters et al. (2015) Masters, D., Capak, P., Stern, D., et al. 2015, ApJ, 813, 53
  * Mead et al. (2021) Mead, A. J., Brieden, S., Tröster, T., & Heymans, C. 2021, MNRAS, 502, 1401
  * Mead et al. (2015) Mead, A. J., Peacock, J. A., Heymans, C., Joudaki, S., & Heavens, A. F. 2015, MNRAS, 454, 1958
  * Miller et al. (2013) Miller, L., Heymans, C., Kitching, T. D., et al. 2013, MNRAS, 429, 2858
  * Miller et al. (2007) Miller, L., Kitching, T. D., Heymans, C., Heavens, A. F., & van Waerbeke, L. 2007, MNRAS, 382, 315
  * Mörtsell & Dhawan (2018) Mörtsell, E. & Dhawan, S. 2018, J. Cosmology Astropart. Phys., 2018, 025
  * Naim et al. (1997) Naim, A., Ratnatunga, K. U., & Griffiths, R. E. 1997, arXiv e-prints, arXiv:9704012
  * Nicola et al. (2021) Nicola, A., García-García, C., Alonso, D., et al. 2021, J. Cosmology Astropart. Phys., 2021, 067
  * Olive & Particle Data Group (2014) Olive, K. A. & Particle Data Group. 2014, Chinese Physics C, 38, 090001
  * Park & Rozo (2020) Park, Y. & Rozo, E. 2020, MNRAS, 499, 4638
  * Peebles (1973) Peebles, P. J. E. 1973, ApJ, 185, 413
  * Planck Collaboration et al. (2020a) Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2020a, A&A, 641, A6
  * Planck Collaboration et al. (2020b) Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2020b, A&A, 641, A5
  * Raichoor et al. (2014) Raichoor, A., Mei, S., Erben, T., et al. 2014, ApJ, 797, 102
  * Reid et al. (2016) Reid, B., Ho, S., Padmanabhan, N., et al. 2016, MNRAS, 455, 1553
  * Riess et al. (2011) Riess, A. G., Macri, L., Casertano, S., et al. 2011, ApJ, 730, 119, [Erratum: ApJ 732, 129 (2011)]
  * Riess et al. (2016) Riess, A. G., Macri, L. M., Hoffmann, S. L., et al. 2016, ApJ, 826, 56
  * Riess et al. (2020) Riess, A. G., Yuan, W., Casertano, S., Macri, L. M., & Scolnic, D. 2020, ApJ, 896, L43
  * Ross et al. (2014) Ross, A. J., Samushia, L., Burden, A., et al. 2014, MNRAS, 437, 1109
  * Sánchez et al. (2017a) Sánchez, A. G., Grieb, J. N., Salazar-Albornoz, S., et al. 2017a, MNRAS, 464, 1493
  * Sánchez et al. (2017b) Sánchez, A. G., Scoccimarro, R., Crocce, M., et al. 2017b, MNRAS, 464, 1640
  * Schirmer (2013) Schirmer, M. 2013, ApJS, 209, 21
  * Schlafly & Finkbeiner (2011) Schlafly, E. F. & Finkbeiner, D. P. 2011, ApJ, 737, 103
  * Schlegel et al. (1998) Schlegel, D. J., Finkbeiner, D. P., & Davis, M. 1998, ApJ, 500, 525
  * Schneider et al. (2002) Schneider, P., van Waerbeke, L., Kilbinger, M., & Mellier, Y. 2002, A&A, 396, 1
  * Schulten & Gordon (1976) Schulten, K. & Gordon, R. 1976, Computer Physics Communications, 11, 269
  * Scolnic et al. (2018) Scolnic, D. M., Jones, D. O., Rest, A., et al. 2018, ApJ, 859, 101
  * Secco et al. (2022) Secco, L. F., Samuroff, S., Krause, E., et al. 2022, Phys. Rev. D, 105, 023515
  * Seljak & Zaldarriaga (1997) Seljak, U. & Zaldarriaga, M. 1997, Phys. Rev. Lett., 78, 2054
  * Sellentin & Heavens (2016) Sellentin, E. & Heavens, A. F. 2016, MNRAS, 456, L132
  * Spergel et al. (2015) Spergel, D., Gehrels, N., Baltay, C., et al. 2015, arXiv e-prints, arXiv:1503.03757
  * Taylor & Kitching (2010) Taylor, A. N. & Kitching, T. D. 2010, MNRAS, 408, 865
  * Taylor et al. (2019) Taylor, P. L., Kitching, T. D., Alsing, J., et al. 2019, Phys. Rev. D, 100, 023519
  * Tröster et al. (2021) Tröster, T., Asgari, M., Blake, C., et al. 2021, A&A, 649, A88
  * Tröster et al. (2022) Tröster, T., Mead, A. J., Heymans, C., et al. 2022, A&A, 660, A27
  * Tröster et al. (2020) Tröster, T., Sánchez, A. G., Asgari, M., et al. 2020, A&A, 633, L10
  * Troxel et al. (2018) Troxel, M. A., MacCrann, N., Zuntz, J., et al. 2018, Phys. Rev. D, 98, 043528
  * Tutusaus et al. (2020) Tutusaus, I., Martinelli, M., Cardone, V. F., et al. 2020, A&A, 643, A70
  * Upham et al. (2021) Upham, R. E., Brown, M. L., & Whittaker, L. 2021, MNRAS, 503, 1999
  * Upham et al. (2022) Upham, R. E., Brown, M. L., Whittaker, L., et al. 2022, A&A, 660, A114
  * Upham et al. (2020) Upham, R. E., Whittaker, L., & Brown, M. L. 2020, MNRAS, 491, 3165
  * van den Busch et al. (2020) van den Busch, J. L., Hildebrandt, H., Wright, A. H., et al. 2020, A&A, 642, A200
  * van Uitert et al. (2018) van Uitert, E., Joachimi, B., Joudaki, S., et al. 2018, MNRAS, 476, 4662
  * Verde et al. (2019) Verde, L., Treu, T., & Riess, A. G. 2019, Nature Astronomy, 3, 891
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261
  * Wes McKinney (2010) Wes McKinney. 2010, in Proceedings of the 9th Python in Science Conference, ed. Stéfan van der Walt & Jarrod Millman, 56 – 61
  * Wright et al. (2019) Wright, A. H., Hildebrandt, H., Kuijken, K., et al. 2019, A&A, 632, A34
  * Wright et al. (2020) Wright, A. H., Hildebrandt, H., van den Busch, J. L., & Heymans, C. 2020, A&A, 637, A100
  * Xavier et al. (2016) Xavier, H. S., Abdalla, F. B., & Joachimi, B. 2016, MNRAS, 459, 3693

## Affiliations

1 Department of Physics and Astronomy, University College London, Gower
Street, London WC1E 6BT, UK
2 Astrophysics Group, Blackett Laboratory, Imperial College London, London SW7
2AZ, UK
3 Institute for Astronomy, University of Edinburgh, Royal Observatory,
Blackford Hill, Edinburgh EH9 3HJ, UK
4 Jodrell Bank Centre for Astrophysics, Department of Physics and Astronomy,
University of Manchester, Oxford Road, Manchester M13 9PL, UK
5 Mullard Space Science Laboratory, University College London, Holmbury St
Mary, Dorking, Surrey RH5 6NT, UK
6 E. A. Milne Centre, University of Hull, Cottingham Road, Hull, HU6 7RX, UK
7 Ruhr-Universität Bochum, Astronomisches Institut, German Centre for
Cosmological Lensing (GCCL), Universitätsstr. 150, 44801 Bochum, Germany
8 Center for Theoretical Physics, Polish Academy of Sciences, al. Lotników
32/46, 02-668 Warsaw, Poland
9 Astronomisches Institut, Ruhr-Universität Bochum, Universitätsstr. 150,
44801 Bochum, Germany
${}^{1}0$ Shanghai Astronomical Observatory (SHAO), Nandan Road 80, Shanghai
200030, China
${}^{1}1$ University of Chinese Academy of Sciences, Beijing 100049, China
${}^{1}2$ Institute of Cosmology and Gravitation, University of Portsmouth,
Portsmouth PO1 3FX, UK
${}^{1}3$ INAF-Osservatorio di Astrofisica e Scienza dello Spazio di Bologna,
Via Piero Gobetti 93/3, I-40129 Bologna, Italy
${}^{1}4$ Max Planck Institute for Extraterrestrial Physics, Giessenbachstr.
1, D-85748 Garching, Germany
${}^{1}5$ INAF-Osservatorio Astrofisico di Torino, Via Osservatorio 20,
I-10025 Pino Torinese (TO), Italy
${}^{1}6$ Department of Mathematics and Physics, Roma Tre University, Via
della Vasca Navale 84, I-00146 Rome, Italy
${}^{1}7$ INFN-Sezione di Roma Tre, Via della Vasca Navale 84, I-00146, Roma,
Italy
${}^{1}8$ INAF-Osservatorio Astronomico di Capodimonte, Via Moiariello 16,
I-80131 Napoli, Italy
${}^{1}9$ INAF-IASF Milano, Via Alfonso Corti 12, I-20133 Milano, Italy
${}^{2}0$ Institut de Física d’Altes Energies (IFAE), The Barcelona Institute
of Science and Technology, Campus UAB, 08193 Bellaterra (Barcelona), Spain
${}^{2}1$ Port d’Informació Científica, Campus UAB, C. Albareda s/n, 08193
Bellaterra (Barcelona), Spain
${}^{2}2$ INAF-Osservatorio Astronomico di Roma, Via Frascati 33, I-00078
Monteporzio Catone, Italy
${}^{2}3$ Department of Physics ”E. Pancini”, University Federico II, Via
Cinthia 6, I-80126, Napoli, Italy
${}^{2}4$ INFN section of Naples, Via Cinthia 6, I-80126, Napoli, Italy
${}^{2}5$ Dipartimento di Fisica e Astronomia ”Augusto Righi” - Alma Mater
Studiorum Universitá di Bologna, Viale Berti Pichat 6/2, I-40127 Bologna,
Italy
${}^{2}6$ INAF-Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5, I-50125,
Firenze, Italy