# KiDS and Euclid: Cosmological implications of a pseudo angular power spectrum analysis of KiDS-1000 cosmic shear tomography

A. Loureiro, L. Whittaker, A. Spurio Mancini, B. Joachimi, A. Cuceu, M.

Asgari, B. Stölzner, T. Tröster, A.H. Wright, M. Bilicki, et al.

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

A. Loureiro, L. Whittaker, A. Spurio Mancini, B. Joachimi, A. Cuceu, et al., KiDS and Euclid: Cosmological implications of a pseudo angular power spectrum analysis of KiDS-1000 cosmic shear tomography. Astronomy and Astrophysics - A&A, 2022, 665, pp.A56. 10.1051/0004-6361/202142481. hal-03839670

HAL Id: hal-03839670

https://hal.science/hal-03839670

Submitted on 10 Mar 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# KiDS and Euclid: Cosmological implications of a pseudo angular power spectrum analysis of KiDS-1000 cosmic shear tomography*

A. Loureiro<sup>1,2,3</sup>, L. Whittaker<sup>1,4</sup>, A. Spurio Mancini<sup>1,5</sup>, B. Joachimi<sup>1</sup>, A. Cuceu<sup>1</sup>, M. Asgari<sup>6,3</sup>, B. Stölzner<sup>1</sup>, T. Tröster<sup>3</sup>, A. H. Wright<sup>7</sup>, M. Bilicki<sup>8</sup>, A. Dvornik<sup>9</sup>, B. Giblin<sup>3</sup>, C. Heymans<sup>3,7</sup>, H. Hildebrandt<sup>7</sup>, H. Shan<sup>10,11</sup>, A. Amara<sup>12</sup>, N. Auricchio<sup>13</sup>, C. Bodendorf<sup>14</sup>, D. Bonino<sup>15</sup>, E. Branchini<sup>16,17</sup>, M. Brescia<sup>18</sup>, V. Capobianco<sup>15</sup>, C. Carbone<sup>19</sup>, J. Carretero<sup>20,21</sup>, M. Castellano<sup>22</sup>, S. Cavuoti<sup>23,18,24</sup>, A. Cimatti<sup>25,26</sup>, R. Cledassou<sup>27,28</sup>, G. Congedo<sup>3</sup>, L. Conversi<sup>29,30</sup>, Y. Copin<sup>31</sup>, L. Corcione<sup>15</sup>, M. Cropper<sup>5</sup>, A. Da Silva<sup>32,33</sup>, M. Douspis<sup>34</sup>, F. Dubath<sup>35</sup>, C. A. J. Duncan<sup>36</sup>, X. Dupac<sup>30</sup>, S. Dusini<sup>37</sup>, S. Farrens<sup>38</sup>, S. Ferriol<sup>31</sup>, P. Fosalba<sup>39,40</sup>, M. Frailis<sup>41</sup>, E. Franceschi<sup>13</sup>, M. Fumana<sup>19</sup>, B. Garilli<sup>19</sup>, B. Gillis<sup>3</sup>, C. Giocoli<sup>42,43</sup>, A. Grazian<sup>44</sup>, F. Grupp<sup>14,45</sup>, S. V. H. Haugan<sup>46</sup>, W. Holmes<sup>47</sup>, F. Hormuth<sup>48,49</sup>, K. Jahnke<sup>49</sup>, M. Kümmel<sup>45</sup>, S. Kermiche<sup>50</sup>, A. Kiessling<sup>47</sup>, M. Kilbinger<sup>38</sup>, T. Kitching<sup>5</sup>, K. Kuijken<sup>51</sup>, M. Kunz<sup>52</sup>, H. Kurki-Suonio<sup>53</sup>, S. Ligori<sup>15</sup>, P.B. Lilje<sup>46</sup>, I. Lloro<sup>54</sup>, O. Mansutti<sup>41</sup>, O. Marggraf<sup>55</sup>, K. Markovic<sup>47</sup>, F. Marulli<sup>13,56,57</sup>, R. Massey<sup>58,59</sup>, M. Meneghetti<sup>13,56</sup>, G. Meylan<sup>60</sup>, M. Moresco<sup>57,13</sup>, B. Morin<sup>38</sup>, L. Moscardini<sup>57,13,56</sup>, E. Munari<sup>41</sup>, S. M. Niemi<sup>61</sup>, C. Padilla<sup>20</sup>, S. Paltani<sup>35</sup>, F. Pasian<sup>41</sup>, K. Pedersen<sup>62</sup>, V. Pettorino<sup>38</sup>, S. Pires<sup>38</sup>, M. Poncet<sup>27</sup>, L. Popa<sup>63</sup>, F. Raison<sup>14</sup>, J. Rhodes<sup>47</sup>, H. Rix<sup>49</sup>, M. Roncarelli<sup>57,13</sup>, R. Saglia<sup>14,45</sup>, P. Schneider<sup>55</sup>, A. Secroun<sup>50</sup>, S. Serrano<sup>39,40</sup>, C. Sirignano<sup>65,37</sup>, G. Sirri<sup>56</sup>, L. Stanco<sup>37</sup>, J. L. Starck<sup>38</sup>, P. Tallada-Crespi<sup>66,21</sup>, A. N. Taylor<sup>3</sup>, I. Tereno<sup>32,67</sup>, R. Toledo-Moreo<sup>68</sup>, F. Torradeflot<sup>21,66</sup>, E. A. Valentijn<sup>69</sup>, Y. Wang<sup>70</sup>, N. Welikala<sup>3</sup>, J. Weller<sup>14,45</sup>, G. Zamorani<sup>13</sup>, J. Zoubian<sup>50</sup>, S. Andreon<sup>71</sup>, M. Baldi<sup>57,13,56</sup>, S. Camera<sup>72,15,73</sup>, R. Farinelli<sup>74</sup>, G. Polenta<sup>75</sup>, and N. Tessore<sup>1</sup>

*(Affiliations can be found after the references)*

Received 19 October 2021 / Accepted 27 June 2022

## ABSTRACT

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_l$) 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_l$ using eight bands in the multipole range $76 < l < 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_l$ estimates by forward modelling the effects of the mask. Assuming a spatially flat ΛCDM cosmology, we constrained the structure growth parameter $S_8 = \sigma_8(\Omega_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.032}^{+0.006}$. These results are in very good agreement with previous KiDS-1000 and SDSS analyses and confirm a ~3$\sigma$ tension with early-Universe constraints from cosmic microwave background experiments.

**Key words.** gravitational lensing: weak – cosmology: observations – large-scale structure of Universe – cosmological parameters

## 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 Λ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 VI 2020), which demonstrated an unprecedented constraining power, mea-

suring some of these parameters with percent-level precision. The *Planck* Mission, however, did not contradict the standard model of cosmology. Instead, it established the ΛCDM model as a difficult paradigm to overturn for other cosmological probes and future surveys (Efstathiou & Gratton 2020, 2021).

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 Λ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

<footnote>* This paper is published on behalf of the KiDS Collaboration and the *Euclid* Consortium.</footnote>
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, 2020; Bernal et al. 2016; Di Valentino et al. 2016; Lin & Ishak 2017; Efstathiou 2020), and the structure growth parameter, $S_8 = \sigma_8\sqrt{\Omega_m}/0.3$, when measured using cosmic shear surveys (Heymans et al. 2013, 2021; MacCrann et al. 2015; Joudaki et al. 2017, 2020; Hildebrandt et al. 2017; Abbott et al. 2018; Park & Rozo 2020; Hikage et al. 2019; Lemos et al. 2021; Asgari et al. 2020, 2021; Tröster et al. 2021; Amon et al. 2022; Secco et al. 2022; DES Collaboration 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 $\ge 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<fnref n="1" />. 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 Collaboration, 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 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–13 redshift tomographic bins (Euclid Collaboration 2020; Tutusaus et al. 2020). This would require dealing with a mixing matrix with 2.1–3.51 Million dimensions, posing potential numerical inaccuracies which might well exceed the very tight accuracy requirements

<footnote n="1">¹ With 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).</footnote>
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 deg<fnref n="2" /> 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 × 2pt) analysis for a flat ΛCDM model, following the methodology described by Joachimi et al. (2021b). Finally, constraints on extensions to ΛCDM using the same 3 × 2pt 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, 2021; Gil-Marín et al. 2020; Bautista et al. 2021; du Mas des Bourboux et al. 2020). 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 VI 2020) 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 sur-

vey by the European Southern Observatory, currently at its 4th Data Release (Kuijken et al. 2019)<fnref n="2" />, with an observed area of around 1000 deg<fnref n="2" />. 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 (ZYJHKs) from VIKING<fnref n="3" /> (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)<fnref n="4" />, covering an effective unmasked area of 777.4 deg<fnref n="2" /> and containing a total of 21 262 011 lensed galaxies. Detailed descriptions of the catalogue's construction, including the LensFit (Miller et al. 2007, 2013; Kitching et al. 2008) 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σ for $S_8 = \sigma_8 \sqrt{\Omega_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_B$. The first four bins have $\Delta z_B^{(1-4)} = 0.2$, ranging from $0.1 < z_B \le 0.9$, while the last tomographic bin ranges from $0.9 < z_B < 1.2$, that is $\Delta z_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 A.1 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^{corr}$, obtained from the observed galaxy ellipticities, $\epsilon^{obs}$. In the weak lensing regime ($|\gamma| \ll 1$), $\epsilon^{corr}$ is estimated by taking into account the additive bias $c$,

$$ \epsilon^{corr} = \epsilon^{obs} - c, \quad (1) $$

where the additive bias, $c = \langle \epsilon^{obs} \rangle$, is assumed to be the weighted average observed ellipticity over all galaxies in a given

<footnote n="2">² http://kids.strw.leidenuniv.nl/DR4/</footnote>

<footnote n="3">³ The VISTA Kilo-degree Infrared Galaxy (VIKING) Survey.</footnote>

<footnote n="4">⁴ http://kids.strw.leidenuniv.nl/DR4/lensing.php</footnote>
<figure><img src="image_3.png" /><figcaption>Fig. 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<sub>B</sub>.</figcaption></figure>

tomographic bin. We note that all these quantities (shear, ellipticity and additive bias) are complex quantities, where $\epsilon = \epsilon_1 + i\epsilon_2$.

## 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{n}$ is given as

$$ \gamma(\hat{n}) = \frac{1}{(1 + m_z)} \frac{\sum_{i \in \hat{n}} w_i \epsilon_i^{\text{corr}}}{\sum_{i \in \hat{n}} w_i}, \quad (2) $$

where the summation is taken over all galaxies in the pixel. For each galaxy, $\epsilon_i^{\text{corr}}$ is constructed using Eq. (1), with weights, $w_i$, provided by shape measurement techniques such as LensFit (Miller et al. 2007, 2013; Kitching et al. 2008). 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 2Y_{\ell m}$, is defined as

$$ \gamma(\hat{n}) = \gamma_1(\hat{n}) + i\gamma_2(\hat{n}) = \sum_{\ell m} (E_{\ell m} + iB_{\ell m}) + 2Y_{\ell m}(\hat{n}); \quad (3) $$

$$ \gamma^*(\hat{n}) = \gamma_1(\hat{n}) - i\gamma_2(\hat{n}) = \sum_{\ell m} (E_{\ell m} - iB_{\ell m}) - 2Y_{\ell m}(\hat{n}), \quad (4) $$

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

$$ E_{\ell m} = \frac{1}{2} \int d\Omega_{\hat{n}} \{ \gamma(\hat{n}) + 2Y_{\ell m}^* + \gamma^*(\hat{n}) - 2Y_{\ell m}^* \}; \quad (5) $$

$$ B_{\ell m} = \frac{-i}{2} \int d\Omega_{\hat{n}} \{ \gamma(\hat{n}) + 2Y_{\ell m}^* - \gamma^*(\hat{n}) - 2Y_{\ell m}^* \}. \quad (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 $W(\hat{n})$; where $W(\hat{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, $W(\hat{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 $W(\hat{\mathbf{n}})$ need careful consideration on how to deal with a binary mask and shear weights for high-resolution maps (Euclid Collaboration 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)

$$ \begin{aligned} \tilde{\gamma}_1(\hat{\mathbf{n}}) \pm i \tilde{\gamma}_2(\hat{\mathbf{n}}) &\equiv W(\hat{\mathbf{n}}) [\gamma_1(\hat{\mathbf{n}}) \pm i \gamma_2(\hat{\mathbf{n}})] \\ &= \sum_{\ell m} (\tilde{E}_{\ell m} \pm i \tilde{B}_{\ell m}) \pm 2 Y_{\ell m}(\hat{\mathbf{n}}), \end{aligned} \quad (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

$$ \tilde{E}_{\ell m} = \sum_{\ell' m'} (E_{\ell m} \tilde{W}_{\ell \ell' mm'}^{+} + B_{\ell m} \tilde{W}_{\ell \ell' mm'}^{-}); \quad (8) $$

$$ \tilde{B}_{\ell m} = \sum_{\ell' m'} (B_{\ell m} \tilde{W}_{\ell \ell' mm'}^{+} - E_{\ell m} \tilde{W}_{\ell \ell' mm'}^{-}), \quad (9) $$

where the mixing kernels $\tilde{W}_{\ell \ell' mm'}^{\pm}$ 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

$$ \tilde{C}_{\ell}^{E_i E_j} \approx \frac{1}{2\ell + 1} \sum_m \tilde{E}_{\ell m}^i \tilde{E}_{\ell m}^j - \langle \tilde{\mathcal{N}}_{\ell}^{i,j} \rangle \delta_{i,j}, \quad (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:

$$ \epsilon_1^R = \epsilon_1^{\text{corr}} \cos(\theta^R) - \epsilon_2^{\text{corr}} \sin(\theta^R), \quad (11) $$

$$ \epsilon_2^R = \epsilon_2^{\text{corr}} \cos(\theta^R) + \epsilon_1^{\text{corr}} \sin(\theta^R). \quad (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 Eqs. (11) and (12) together with Eq. (2) to define the randomised shear estimates, $\gamma^R$, one can propagate this quantity through Eqs. (7) and (8) to obtain $\tilde{E}_{\ell m}^{R,i}$. Finally, the $\tilde{\mathcal{N}}_{\ell}$ term from Eq. (10) is simply defined as

$$ \tilde{\mathcal{N}}_{\ell}^{R,i} = \frac{1}{2\ell + 1} \sum_m \tilde{E}_{\ell m}^{R,i} \tilde{E}_{\ell m}^{R,j}. \quad (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_{\epsilon_{\text{obs}}}^2 = \sigma_{\gamma}^2 + \sigma_{\epsilon}^2$, where the shear field variance,

$\sigma_{\gamma}^2$, is mixed with the variance from the intrinsic shape and measurement errors, $\sigma_{\epsilon}^2$. 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

The mixing of modes, introduced in Eqs. (8) and (9) via the $\tilde{W}_{\ell \ell' mm'}^{\pm}$ 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, $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 $M_{\ell\ell'}$, which contains information about the survey geometry in harmonic space (Peebles 1973; Brown et al. 2005; Hikage et al. 2011):

$$ \langle \tilde{\mathcal{C}}_{\ell}^{ij} \rangle = \sum_{\ell'} \mathcal{M}_{\ell \ell'}^{ij} \mathcal{C}_{\ell'}^{ij}, \quad (14) $$

or, assuming no parity violation modes are present ($\tilde{C}_{\ell}^{E_i E_j} = \tilde{C}_{\ell}^{B_i B_j}$),

$$ \langle \tilde{\mathcal{C}}_{\ell}^{E_i E_j} \rangle = \sum_{\ell'} \left( \begin{matrix} \tilde{W}_{\ell \ell'}^{++}(i, j) & \tilde{W}_{\ell \ell'}^{--}(i, j) \\ \tilde{W}_{\ell \ell'}^{--}(i, j) & \tilde{W}_{\ell \ell'}^{++}(i, j) \end{matrix} \right) \left( \begin{matrix} C_{\ell \ell'}^{E_i E_j} \\ C_{\ell \ell'}^{B_i B_j} \end{matrix} \right). \quad (15) $$

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

$$ \tilde{W}_{\ell \ell'}^{\pm \pm}(i, j) = \frac{1}{2\ell + 1} \sum_{mm'} \tilde{W}_{\ell \ell' mm'}^{\pm, i} (\tilde{W}_{\ell \ell' mm'}^{\pm, j})^*, \quad (16) $$

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

$$ \begin{aligned} \tilde{W}_{\ell \ell' mm'}^{\pm, i} &\equiv \int d\Omega_{\hat{\mathbf{n}}} \pm 2 Y_{\ell' m'}(\hat{\mathbf{n}}) W^{i}(\hat{\mathbf{n}}) \pm 2 Y_{\ell m}(\hat{\mathbf{n}}); \\ &= \sum_{\ell'' m''} (-1)^m \sqrt{\frac{(2\ell + 1)(2\ell' + 1)(2\ell'' + 1)}{4\pi}} \tilde{W}_{\ell \ell'' m''}^{i} \\ &\quad \times \begin{pmatrix} \ell & \ell' & \ell'' \\ \pm 2 & \mp 2 & 0 \end{pmatrix} \begin{pmatrix} \ell & \ell' & \ell'' \\ m & m' & m'' \end{pmatrix}, \end{aligned} \quad (17) $$
<figure><img src="image_4.png" /><figcaption>Fig. 2. KiDS-1000 full mixing matrix, $M_{\ell\ell'}$ from Eqs. (14) and (15), for the third redshift tomographic bin. For this case, $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, $W^{++}$ and $W^{--}$, we show the matrices calculated in the range $2 \le \ell \le 3070$. From Eqs. (16) and (17) one can see that these objects are not symmetric in $\ell-\ell'$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 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 \text{ Mpc}^{-1}$ for a given $\ell$ and $z$, where $\chi(z)$ is the co-moving distance. Bottom: the KiDS-1000 pseudo-$C_\ell$ binned mixing matrices (see Fig. 2, for example). These are convolved with the theory $C_\ell$ (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 \le \ell \le 1500$.</figcaption></figure>

where the $2 \times 3$ objects above are the Wigner 3j symbols, calculated using UCLWig3j library<fnref n="5" />, and

$$W_{\ell m}^{i} = \oint d\Omega_{\hat{n}} W^{i}(\hat{n}) Y_{\ell m}^{*}(\hat{n}) \qquad (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{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 $W(\hat{n})$, would cause an imprint from the clustering of galaxies into the mixing matrix, causing an artificial excess power for $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_{\text{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 informa-

tion 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

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 \le \ell \le 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)<fnref n="6" />. Observing the binned mixing matrices shown in Fig. 3, we note that by cutting at an $\ell_{\text{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_{\text{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

<footnote n="5">⁵ https://github.com/LorneWhiteway/UCLWig3j – this library optimises the calculation of Wigner 3j symbols using the recurrence relation by Schulten & Gordon (1976).</footnote>

<footnote n="6">⁶ For comparison, the scales used in Asgari et al. (2021) for band-power measurements were $100 \le \ell \le 1500$.</footnote>
<figure><img src="image_6.png" /><figcaption>Fig. 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 \le \ell \le 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.</figcaption></figure>

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)<fnref n="7" />, the measured pseudo-$C_\ell$ are then binned using eight log-spaced multipole bins between $76 \le \ell \le 1500$,

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

where $\ell^L < \ell \le \ell^{L+1}$ are the boundaries of the bandwidth bins centred in $L$, shown in Fig. 3 as the vertical dotted lines. The mea-

sured 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^2_{red} = 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).

<footnote n="7">⁷ We note that this choice of weight could be optimised for cosmic shear analysis, which will be investigated in future work.</footnote>
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 Sect. 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, 2021b; 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; 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 *Flask*<fnref n="8" /> (Xavier et al. 2016) and *Salmo*<fnref n="9" /> 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.

*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⁻¹, 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 Sect. 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_{\text{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_{\text{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)

$$ \text{Cov}_m(\tilde{C}_\ell^{i,j}; \tilde{C}_{\ell'}^{i',j'}) = \tilde{C}_\ell^{i,j} \tilde{C}_{\ell'}^{i',j'} \left[ \sigma_m^i \sigma_m^{i'} + \sigma_m^j \sigma_m^{j'} + \sigma_m^i \sigma_m^{j'} + \sigma_m^j \sigma_m^{i'} \right], \quad (20) $$

<footnote n="8">⁸ http://www.astro.iag.usp.br/~flask/</footnote>

<footnote n="9">⁹ https://github.com/Linc-tw/salmo</footnote>
<table><caption>**Table 1.** Comparison between the original mocks, the sub-sampled mocks and the Gold Sample SOM.</caption><thead><tr><th>Bin</th><th>z̄ (J21)</th><th>z̄ (Sub-sampled)</th><th>z̄ (Gold Sample)</th><th>n<sub>eff</sub> (J21)</th><th>n<sub>eff</sub> (Sub-sampled)</th><th>n<sub>eff</sub> (Gold Sample)</th></tr></thead><tbody><tr><td>1</td><td>0.39</td><td>0.40</td><td>0.26</td><td>0.85</td><td>0.62</td><td>0.62</td></tr><tr><td>2</td><td>0.49</td><td>0.47</td><td>0.40</td><td>1.56</td><td>1.19</td><td>1.18</td></tr><tr><td>3</td><td>0.67</td><td>0.62</td><td>0.56</td><td>2.23</td><td>1.86</td><td>1.85</td></tr><tr><td>4</td><td>0.83</td><td>0.80</td><td>0.79</td><td>1.52</td><td>1.27</td><td>1.26</td></tr><tr><td>5</td><td>1.00</td><td>0.96</td><td>0.98</td><td>1.39</td><td>1.32</td><td>1.31</td></tr></tbody></table>

**Notes.** Here 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.

<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Correlation matrix for the E-mode pseudo-$C_ℓ$ 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_{\ell}^{i,j}$ E-modes.</figcaption></figure>

where $\sigma_m^i$ is defined as the standard deviation of the m-corrections shown in Table A.1 and $\tilde{C}_\ell^{E_i E_j}$ 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_m^i \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: $C_{\text{tot}} = C_{\text{mocks}} + C_{\text{m}}$. The resulting covariance, $C_{\text{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.

## 4. Systematic contamination analysis

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 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 (ugriZYJHKs):** 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 (ugriZYJHKs):** 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.

- **μ-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
<figure><img src="image_8.png" /><figcaption>**Fig. 6.** P-values for the cross-correlations between the data's tomographic pseudo E-modes and different systematics for a null detection ($p \le 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 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.</figcaption></figure>

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{n})$; see Sect. 3.1 of Giblin et al. (2021).

- **Stellar Density:** star sample from Gaia DR2¹⁰ containing all Gaia DR2 objects in a HEALPix grid of $N_{side} = 1024$.

For each of the quantities mentioned above, we create HEALPix maps using the catalogue object positions with $N_{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 $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}_{\ell}^{X_{i, \text{syst}}} = \frac{1}{(2\ell + 1)} \sum_{m} X_{\ell m}^{i} (a_{\ell m}^{\text{syst}})^{*}. \qquad (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 resampling method. We start by subdividing the survey footprint into 60 distinct regions using a K-Means algorithm on the sphere (Kwan et al. 2017)<fnref n="11" />. While removing the kth region, we cal-

culate 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.

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. B.1 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_B > 0.7$), or consistent with zero (for $z_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. B.1) 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 Sect. 6.1 are purely a result of the cosmological and astrophysical signals captured by our measured pseudo-$C_\ell$.

## 5. Cosmological inference

This section describes the final elements necessary for cosmological inference using tomographic cosmic shear angular power

<footnote n="10">¹⁰ Obtained from https://cds.u-strasbg.fr/gaia</footnote>

<footnote n="11">¹¹ Our K-Means code uses the code kmeans_radec, developed by Erin Sheldon which can be found here: https://github.com/esheldon/kmeans_radec.</footnote>
<figure><img src="image_9.png" /><figcaption>**Fig. 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. (B.1), showing that the largest S/N multipole bandwidth has $S/N \sim -2.5 \times 10^{-2}$, meaning that it is well within the estimated covariance for our analysis.</figcaption></figure>

spectra data: theory forward modelling, the likelihood and the priors we use to obtain the posteriors presented in Sect. 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)<fnref n="12" /> modified to directly sample from Gaussian priors<fnref n="13" /> for the Multinest sampler (Feroz & Hobson 2008; Feroz et al. 2009, 2019; Buchner et al. 2014). 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 Sect. 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)<fnref n="14" />. Following the other KiDS-1000 analyses, we fix the sum of neutrino masses to $\sum m_\nu = 0.06\,\text{eV}c^{-2}$, using the one massive species approximation of the normal hierarchy. Not only was this the approach taken by Planck Collaboration VI (2020), but it was shown in Loureiro & Cuceu (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_m^{\text{nl}}(k)$, is calculated using the halo model from HMCode (Mead et al. 2015) included in Class<fnref n="15" />. 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_{\text{bary}}$. These parameters can be constrained by the following relation from Joudaki et al. (2018): $\eta_0 = 0.98-0.12A_{\text{bary}}$. Therefore, $A_{\text{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):

$$C_{\ell}^{\epsilon_i, \epsilon_j} = C_{\ell}^{\gamma_i, \gamma_j} + C_{\ell}^{\gamma_i, I_j} + C_{\ell}^{I_i, \gamma_j} + C_{\ell}^{I_i, I_j}. \quad (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 $k f_K(\chi) = \ell + 1/2$, we can express the right-hand side $C_{\ell s}$ from Eq. (22) as

$$C_{\ell}^{X_i, Y_j} = \int_{0}^{\chi^{\text{H}}} d\chi \frac{W_{X}^{i}(\chi) W_{Y}^{j}(\chi)}{f_{K}^{2}(\chi)} P_{m}^{\text{nl}}\left(\frac{\ell + 1/2}{f_{K}(\chi)}, z(\chi)\right), \quad (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_H$ is the co-moving distance to the horizon and $f_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_K(\chi) = \chi$. The weak lensing kernel is defined as:

$$W_{\gamma}^{i}(\chi) = \frac{3 H_{0}^{2} \Omega_{m}}{2 c^{2} a(\chi)} f_{K}(\chi) \int_{\chi}^{\chi^{\text{H}}} d\chi' n_{s}^{i}(\chi') \frac{f_{K}(\chi' - \chi)}{f_{K}(\chi')}. \quad (24)$$

<footnote n="14">¹⁴ https://lesgourg.github.io/class_public/class.html</footnote>

<footnote n="15">¹⁵ Mead 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.</footnote>

<footnote n="12">¹² https://github.com/brinckmann/montepython_public</footnote>

<footnote n="13">¹³ The modified version sampling from Gaussian priors can be found at: https://github.com/BStoelzner/montepython_public/tree/gaussianamanuscript</footnote>
The intrinsic alignments kernel is calculated using the non-linear alignment (NLA) model from Bridle & King (2007) without redshift evolution,

$$W_1^i(\chi) = -A_{1A} \frac{C_1 \rho_{cr} \Omega_m}{D[a(\chi)]} n_s^i(\chi), \qquad (25)$$

with $C_1\rho_{cr} \approx 0.0134$ (Joachimi et al. 2011) and where $D[a(\chi)]$ is the growth function and $A_{1A}$, 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 \le \ell \le 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_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 likelihoods<fnref n="16" />, which combine pre-reconstruction full-shape<fnref n="17" /> 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_{\text{bary}}$ and $A_{\text{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_M(z_{\text{eff}})/r_d$, $D_H(z_{\text{eff}})/r_d$ and the rate of structure growth, $f(z_{\text{eff}})\sigma_8(z_{\text{eff}})$. Here, $D_M(z) = (1+z)^{-1}f_K[\chi(z)]$ is the comoving angular diameter distance, $D_H(z) = c/H(z)$, $r_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_M(z_{\text{eff}})/r_d$ and $D_H(z_{\text{eff}})/r_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_M(z_{\text{eff}})/r_d$, $D_H(z_{\text{eff}})/r_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

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,

<footnote n="16">¹⁶ https://svn.sdss.org/public/data/eboss/DR16cosmo/tags/v1_0_0/likelihoods/BAO-plus</footnote>

<footnote n="17">¹⁷ In 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.</footnote>
<table><caption>**Table 2.** Priors on cosmological and nuisance parameters.</caption><thead><tr><th>Parameter name</th><th>Symbol</th><th>Prior</th></tr></thead><tbody><tr><td>Density fluctuation amp.</td><td>S<sub>8</sub></td><td>[0.1, 1.3]</td></tr><tr><td>Cold dark matter density</td><td>Ω<sub>c</sub>h<sup>2</sup></td><td>[0.051, 0.255]</td></tr><tr><td>Baryonic matter density</td><td>Ω<sub>b</sub>h<sup>2</sup></td><td>[0.019, 0.026]</td></tr><tr><td>Hubble constant</td><td>h</td><td>[0.64, 0.82]</td></tr><tr><td>Spectral index</td><td>n<sub>s</sub></td><td>[0.84, 1.1]</td></tr><tr><td>Intrinsic alignment amp.</td><td>A<sub>IA</sub></td><td>[-6, 6]</td></tr><tr><td>Baryonic feedback amp.</td><td>A<sub>bary</sub></td><td>[2, 3.13]</td></tr><tr><td>Redshift displacement</td><td>δ<sub>z</sub></td><td>N(μ, C<sub>z</sub>)</td></tr></tbody></table>

**Notes.** The $\Omega_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, $N(\mu, C_z)$ is a normal distribution with mean $\mu$ and covariance $C_z$ as estimated by Hildebrandt et al. (2021).

<table><caption>**Table 3.** Summary of the main parameters constrained by cosmic shear data and combinations with galaxy clustering.</caption><thead><tr><th>Parameter</th><th>KiDS-1000 PC<sub>ℓ</sub></th><th>KiDS-1000 PC<sub>ℓ</sub> + SDSS</th></tr></thead><tbody><tr><td>S<sub>8</sub></td><td>0.754<sup>+0.027</sup><sub>-0.029</sub></td><td>0.771<sup>+0.006</sup><sub>-0.032</sub></td></tr><tr><td>σ<sub>8</sub></td><td>0.920<sup>+0.159</sup><sub>-0.215</sub></td><td>0.756<sup>+0.029</sup><sub>-0.044</sub></td></tr><tr><td>Ω<sub>m</sub></td><td>0.202<sup>+0.127</sup><sub>-0.053</sub></td><td>0.312<sup>+0.011</sup><sub>-0.018</sub></td></tr></tbody></table>

**Notes.** Here, we are considering the 68% C.I. ($1\sigma$) MAP+PJ-PHD. It is important to note that $\sigma_8$ and $\Omega_m$ are not properly constrained by cosmic shear alone.

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_s$ or $\ln(10^{10}A_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_c = \Omega_c h^2$, derived from Supernova Type Ia constraints on $\Omega_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_b = \Omega_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 $\text{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 $\mu = (0.0001, 0.0021, 0.0129, 0.0110, -0.0060)$ and a covariance $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. A.2.

# 6. Results

## 6.1. Cosmological posteriors

We present our flat $\Lambda$CDM constraints for the pseudo-$C_ℓ$ 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 TTTEE+lowE<fnref n="18" /> (Planck Collaboration VI 2020; Planck Collaboration V 2020), 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_m$, and $S_8$, using MAP+PJ-HPD constraints are presented in Table 3 for KiDS-1000 pseudo-$C_ℓ$ and its combination with SDSS BAO and RSD data. A more complete summary of results, including the traditional marginal constraints, can be found in Table A.2. Unless stated otherwise, the quoted credible intervals in this section are the 68% credible intervals (CI) using the MAP+PJ-HPD 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_m$. Results between the KiDS-1000 estimators are consistent, with the pseudo-$C_ℓ$ 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_m/0.3}$ parameter instead, as shown in Fig. A.3. Here, we find $S_8^{\text{SDSS}} = 0.812^{+0.050}_{-0.049}$ for the SDSS BAO and RSD data set which is in agreement with both KiDS-1000 pseudo-$C_ℓ$, $S_8^{\text{PCL}} = 0.754^{+0.027}_{-0.029}$, 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 $\text{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_m$ parameter. This demonstrates the $S_8$ measurement robustness from KiDS-1000 cosmic shear analyses. The $S_8$ constraints with the pseudo-$C_ℓ$ estimator are slightly larger than the one found with band-powers from Asgari et al. (2021), reflecting a 3.7%

<footnote n="18">¹⁸ Here we used the plikann_TTTEE+lowL+lowE likelihood installed within Monte Python 3.</footnote>
<figure><img src="image_10.png" /><figcaption>**Fig. 8.** Marginalised 2D posterior distributions for 68% (darker) and 95% (lighter) C.I. in the $\Omega_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 TTTEE low-E combination (green). We also compare it with the 3 × 2pt analysis from Heymans et al. (2021) (black-dashed).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 9.** Same as Fig. 8, but for the marginalised $\Omega_m$ and $S_8$ constraints.</figcaption></figure>

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α 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)<fnref n="19" /> 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:

$$[\text{PCL}] A_{\text{IA}} = 0.396^{+0.251}_{-0.931}, \qquad (26)$$

$$[\text{2PCF}] A_{\text{IA}} = 0.387^{+0.321}_{-0.374}, \qquad (27)$$

$$[\text{COSEBIs}] A_{\text{IA}} = 0.264^{+0.424}_{-0.337}, \qquad (28)$$

$$[\text{BP}] A_{\text{IA}} = 0.973^{+0.292}_{-0.383}, \qquad (29)$$

We speculate that this could be the reason why our method finds a slightly lower peak for the matter density parameter, $\Omega_m$, as a similar effect appears in Asgari et al. (2021). Despite this, $\Omega_m$ is not strongly constrained in our analysis, meaning that the uncertainty is large enough that the constraints remain consistent for both $\Omega_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 ~1.2σ 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_{\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_{\text{bary}}$ and $A_{\text{IA}}$, are shown in Fig. 11, and the redshift nuisance parameters are shown in the Appendix in Fig. A.2 with a comparison to the band-power constraints from Asgari et al. (2021). Table A.2 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.

<footnote n="19">¹⁹ Here, we quote the mid nominal values as this analysis is closer to the analysis performed for KV-450 and KiDS-1000.</footnote>
<figure><img src="image_12.png" /><figcaption>**Fig. 10.** Comparison between different cosmological probe constraints on the growth of structure parameter, $S_8 = \sigma_8(\Omega_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<sub>ℓ</sub> 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σ.</figcaption></figure>

## 6.2. Cosmic shear tension with early Universe probes

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

$$ \tau = \frac{|\langle S_8^A \rangle - \langle S_8^B \rangle|}{\sqrt{\sigma^2(S_8^A) + \sigma^2(S_8^B)}}, \qquad (30) $$

where $\langle S_8^A \rangle$ and $\sigma^2(S_8^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^A \rangle$ and $\langle S_8^B \rangle$ is with a case where there is no difference between measurements.

Comparisons between KiDS-1000 pseudo-C<sub>ℓ</sub> measurements and Planck 2018 TTTEE+lowE nominal results (Planck Collaboration VI 2020; Planck Collaboration V 2020) lead to a $2.8\sigma$ tension between late- and early-Universe constraints. When combining the pseudo-C<sub>ℓ</sub> 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 × 2pt 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<sub>ℓ</sub> 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<sub>ℓ</sub> 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<sub>ℓ</sub> measurements and the Planck Legacy results, we also adopt the Hellinger distance measure, $d_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_H$ is defined as

$$ d_H^2[p, q] = \frac{1}{2} \int d\theta [\sqrt{p(\theta)} - \sqrt{q(\theta)}]^2. \qquad (31) $$

With the metric above, we measured $d_H = 0.95$ between our KiDS-1000 pseudo-C<sub>ℓ</sub> 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_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_H$ between KiDS-1000, ACT and combinations with ACT.

## 7. Conclusions

We performed a pseudo angular power spectra (pseudo-C<sub>ℓ</sub>) analysis of the current state-of-the-art galaxy shape measurements from the Kilo Degree Survey’s fourth data release (KiDS-1000,
<figure><img src="image_13.png" /><figcaption>**Fig. 11.** Extended marginalised 1D and 2D posteriors for relevant ΛCDM and astrophysical parameters: the matter density parameter (Ω_m), the amplitude of matter density fluctuations (S_8 = σ_8√Ω_m/0.3), σ_8 as a derived parameter, the physical baryonic matter density (Ω_bh^2), the dimensionless Hubble constant (h), and the spectral index (n_s). Included are also two nuisance parameters: the amplitude of intrinsic alignments (A_{IA}) and the baryonic feedback amplitude (A_{bary}) as these correlate with some cosmological parameters. When combining the KiDS-1000 pseudo-C_l 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 TTTEE low E (Planck Collaboration V 2020) constrains (green) and the KiDS-1000 band-powers constraints (orange) from Asgari et al. (2021).</figcaption></figure>

Kuijken et al. 2019), obtaining cosmological constraints for the current concordance model, the flat ΛCDM cosmology. The KiDS-1000 shear catalogue (Giblin et al. 2021) is now publicly available<fnref n="20" />, containing nine-band photometry with shape and photometric redshift estimates for more than 21 million lensed objects distributed over an area of around 1000 deg².

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_l 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_l estimator performed as well as other 2pt functions, with no flat-sky

<footnote n="20">²⁰ http://kids.strw.leidenuniv.nl/DR4/KiDS-1000_shearcatalogue.php
https://github.com/KiDS-WL/Cat_to_Obs_K1000_P1</footnote>
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 (~2 × 10⁹ 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_{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_{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_{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_{min}$ = 100, while we used $\ell_{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_{IA}$ between the different 2pt functions as measurements agree up to ~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.027}_{-0.029}$, or 9.8 ± 3.3% lower than the $S_8$ nominal constraints from Planck Collaboration VI (2020). 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 ± 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 ~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 systems 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 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. M.A. and T.T. acknowledge support from the European Research Council under grant number 647112. A.H.W. and A.D. acknowledge support from the European Research Council Consolidator Grant (No. 770935). M.B. 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. B.G. acknowledges support from the European Research Council under grant number 647112 and from the Royal Society through an Enhancement Award (RGF/EA/181006). C.H. 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). H.Y.S. 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 Economía 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 (A.L., L.W., A.S.M., B.J., A.C.), 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, T. M. C., Abdalla, F. B., Alarcon, A., et al. 2018, Phys. Rev. D, 98, 043526

Aiola, S., Calabrese, E., Maurin, L., et al. 2020, JCAP, 2020, 047

Alam, S., Ata, M., Bailey, S., et al. 2017, MNRAS, 470, 2617

Alam, S., Aubert, M., Avila, S., et al. 2021, Phys. Rev. D, 103

Alonso, D., Sanchez, J., Slosar, A., & LSST Dark Energy Science Collaboration 2019, MNRAS, 484, 4127

Amon, A., Gruen, D., Troxel, M. A., et al. 2022, Phys. Rev. D, 105, 023514

Asgari, M., Taylor, A., Joachimi, B., & Kitching, T. D. 2018, MNRAS, 479, 454

Asgari, M., Tröster, T., Heymans, C., et al. 2020, A&A, 634, A127

Asgari, M., Lin, C.-A., Joachimi, B., et al. 2021, A&A, 645, A104

Astropy Collaboration (Price-Whelan, A. M., et al.) 2018, AJ, 156, 123

Audren, B., Lesgourgues, J., Benabed, K., & Prunet, S. 2013, JCAP, 1302, 001

Bartelmann, M., & Schneider, P. 2001, Phys. Rep., 340, 291

Bautista, J. E., Paviot, R., Vargas Magaña, M., et al. 2021, MNRAS, 500, 736

Becker, M. R., & Rozo, E. 2016, MNRAS, 457, 304

Becker, M. R., Troxel, M. A., MacCrann, N., et al. 2016, Phys. Rev. D, 94, 022002

Begeman, K., Belikov, A. N., Boxhoorn, D. R., & Valentijn, E. A. 2013, Exp. Astron., 35, 1

Benítez, N. 2000, ApJ, 536, 571

Bennett, C. L., Larson, D., Weiland, J. L., et al. 2013, ApJS, 208, 20

Beran, R. 1977, Ann. Stat., 5, 445

Bernal, J. L., Verde, L., & Riess, A. G. 2016, JCAP, 2016, 019

Bernstein, G. M. 2009, ApJ, 695, 652

Blas, D., Lesgourgues, J., & Tram, T. 2011, JCAP, 2011, 034

Bridle, S., & King, L. 2007, New J. Phys., 9, 444

Bridle, S. L., Crittenden, R., Melchiorri, A., et al. 2002, MNRAS, 335, 1193

Brinckmann, T., & Lesgourgues, J. 2019, Phys. Dark Universe, 24, 100260

Brown, M. L., Castro, P. G., & Taylor, A. N. 2005, MNRAS, 360, 1262

Buchner, J., Georgakakis, A., Nandra, K., et al. 2014, A&A, 564, A125

Coe, D., Benítez, N., Sánchez, S. F., et al. 2006, AJ, 132, 926

Cuceu, A., Farr, J., Lemos, P., & Font-Ribera, A. 2019, JCAP, 2019, 044

de Jong, J. T. A., Verdoes Kleijn, G. A., Boxhoorn, D. R., et al. 2015, A&A, 582, A62

de Jong, J. T. A., Verdoes Kleijn, G. A., Erben, T., et al. 2017, A&A, 604, A134

DES Collaboration (Abbott, T. M. C., et al.) 2022, Phys. Rev. D, 105, 023520

Di Valentino, E., Melchiorri, A., & Silk, J. 2013, Phys. Lett. B, 761, 242

Doux, C., Chang, C., Jain, B., et al. 2021, MNRAS, 503, 3796

du Mas des Bourboux, H., Le Goff, J. M., Blomqvist, M., et al. 2017, A&A, 608, A130

du Mas des Bourboux, H., Rich, J., Font-Ribera, A., et al. 2020, ApJ, 901, 153

Edge, A., Sutherland, W., Kuijken, K., et al. 2013, The Messenger, 154, 32

Efstathiou, G. 2004, MNRAS, 349, 603

Efstathiou, G. 2020, arXiv e-prints [arXiv:2007.10716]

Efstathiou, G., & Gratton, S. 2020, MNRAS, 496, L91

Efstathiou, G., & Gratton, S. 2021, Open J. Astrophys., 4, 8

Erben, T., Schirmer, M., Dietrich, J. P., et al. 2005, Astron. Nachr., 326, 432

Euclid Collaboration (Blanchard, A., et al.) 2020, A&A, 642, A191

Feroz, F., & Hobson, M. P. 2008, MNRAS, 384, 449

Feroz, F., Hobson, M. P., & Bridges, M. 2009, MNRAS, 398, 1601

Feroz, F., Hobson, M. P., Cameron, E., & Pettitt, A. N. 2019, Open J. Astrophys., 2, 10

Fortuna, M. C., Hoekstra, H., Joachimi, B., et al. 2021, MNRAS, 501, 2983

Gaia Collaboration (Brown, A. G. A., et al.) 2018, A&A, 616, A1

Gallozzi, S., Paris, D., Scardia, M., & Dubois, D. 2020a, ArXiv e-prints [arXiv:2003.05472]

Gallozzi, S., Scardia, M., & Maris, M. 2020b, ArXiv e-prints [arXiv:2001.10952]

García-García, C., Alonso, D., & Bellini, E. 2019, JCAP, 2019, 043

García-García, C., Ruiz-Zapatero, J., Alonso, D., et al. 2021, JCAP, 2021, 030

Giblin, B., Heymans, C., Asgari, M., et al. 2021, A&A, 645, A105

Gil-Marín, H., Bautista, J. E., Paviot, R., et al. 2020, MNRAS, 498, 2492

Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759

Hall, A. 2021, MNRAS, 505, 4935

Han, D., Sehgal, N., MacInnis, A., et al. 2021, JCAP, 2021, 031

Handley, W., & Lemos, P. 2021, Phys. Rev. D, 103, 063529

Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357

Heydenreich, S., Schneider, P., Hildebrandt, H., et al. 2020, A&A, 634, A104

Heymans, C., Van Waerbeke, L., Miller, L., et al. 2012, MNRAS, 427, 146

Heymans, C., Grocutt, E., Heavens, A., et al. 2013, MNRAS, 432, 2433

Heymans, C., Tröster, T., Asgari, M., et al. 2021, A&A, 646, A140

Hikage, C., Takada, M., Hamana, T., & Spergel, D. 2011, MNRAS, 412, 65

Hikage, C., Oguri, M., Hamana, T., et al. 2019, PASJ, 71, 43

Hildebrandt, H., Viola, M., Heymans, C., et al. 2017, MNRAS, 465, 1454

Hildebrandt, H., Köhlinger, F., van den Busch, J. L., et al. 2020, A&A, 633, A69

Hildebrandt, H., van den Busch, J. L., Wright, A. H., et al. 2021, A&A, 647, A124

Hivon, E., Górski, K. M., Netterfield, C. B., et al. 2002, ApJ, 567, 2

Hotelling, H. 1931, Ann. Math. Stat., 2, 360

Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90

Joachimi, B., & Bridle, S. L. 2010, A&A, 523, A1

Joachimi, B., Schneider, P., & Eifler, T. 2008, A&A, 477, 43

Joachimi, B., Mandelbaum, R., Abdalla, F. B., et al. 2011, A&A, 527, A26

Joachimi, B., Köhlinger, F., Handley, W., & Lemos, P. 2021a, A&A, 647, L5

Joachimi, B., Lin, C. A., Asgari, M., et al. 2021b, A&A, 646, A129

Joudaki, S., Blake, C., Heymans, C., et al. 2017, MNRAS, 465, 2033

Joudaki, S., Blake, C., Johnson, A., et al. 2018, MNRAS, 474, 4894

Joudaki, S., Hildebrandt, H., Traykova, D., et al. 2020, A&A, 638, L1

Kaiser, N. 1992, ApJ, 388, 272

Kannawadi, A., Hoekstra, H., Miller, L., et al. 2019, A&A, 624, A92

Kilbinger, M. 2015, Rep. Prog. Phys., 78, 086901

Kilbinger, M., Heymans, C., Asgari, M., et al. 2017, MNRAS, 472, 2126
Kitching, T. D., Miller, L., Heymans, C. E., van Waaerbeke, L., & Heavens, A. F. 2008, MNRAS, 390, 149

Kitching, T. D., Alsing, J., Heavens, A. F., et al. 2017, MNRAS, 469, 2737

Kitching, T. D., Taylor, P. L., Capak, P., Masters, D., & Hoekstra, H. 2019a, Phys. Rev. D, 99, 063536

Kitching, T. D., Paykari, P., Hoekstra, H., & Cropper, M. 2019b, Open J. Astrophys., 2, 5

Kitching, T. D., Deshpande, A. C., & Taylor, P. L. 2020, Open J. Astrophys., 3, 14

Kluyver, T., Ragan-Kelley, B., Pérez, F., et al. 2016, in Positioning and Power in Academic Publishing: Players, Agents and Agendas, eds. F. Loizides, & B. Schmidt, (Netherlands: IOS Press), 87

Kogut, A., Spergel, D. N., Barnes, C., et al. 2003, ApJS, 148, 161

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

Krause, E., Fang, X., Pandey, S., et al. 2021, ArXiv e-prints [arXiv:2105.13548]

Kuijken, K., Heymans, C., Hildebrandt, H., et al. 2015, MNRAS, 454, 3500

Kuijken, K., Heymans, C., Dvornik, A., et al. 2019, A&A, 625, A2

Kwan, J., Sánchez, C., Clampitt, J., et al. 2017, MNRAS, 464, 4045

Laureijs, R., Amiaux, J., Arduini, S., et al. 2011, arXiv e-prints [arXiv:1110.3193]

Leistedt, B., & Peiris, H. V. 2014, MNRAS, 444, 2

Lemos, P., Challinor, A., & Efstathiou, G. 2017, JCAP, 2017, 014

Lemos, P., Köhlinger, F., Handley, W., et al. 2020, MNRAS, 496, 4647

Lemos, P., Raveri, M., Campos, A., et al. 2021, MNRAS, 505, 6179

Lesgourgues, J. 2011, arXiv e-prints [arXiv:1104.2932]

Lewis, A. 2019, arXiv e-prints [arXiv:1910.13970]

Li, Y., Singh, S., Yu, B., Feng, Y., & Seljak, U. 2019, JCAP, 2019, 016

Lin, W., & Ishak, M. 2017, Phys. Rev. D, 96, 083532

Lin, C.-H., Harnois-Déraps, J., Eifler, T., et al. 2020, MNRAS, 499, 2977

Loureiro, A., Cuceu, A., et al. 2019, Phys. Rev. Lett., 123, 081301

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

LSST Dark Energy Science Collaboration 2012, arXiv e-prints [arXiv:1211.0310]

MacCrann, N., Zuntz, J., Bridle, S., Jain, B., & Becker, M. R. 2015, MNRAS, 451, 2877

Maehoenen, P. H., & Hakala, P. J. 1995, ApJ, 452, L77

Masters, D., Capak, P., Stern, D., et al. 2015, ApJ, 813, 53

Mead, A. J., Peacock, J. A., Heymans, C., Joudaki, S., & Heavens, A. F. 2015, MNRAS, 454, 1958

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

Miller, L., Kitching, T. D., Heymans, C., Heavens, A. F., & van Waaerbeke, L. 2007, MNRAS, 382, 315

Miller, L., Heymans, C., Kitching, T. D., et al. 2013, MNRAS, 429, 2858

Mörtsell, E., & Dhawan, S. 2018, JCAP, 2018, 025

Naim, A., Ratnatunga, K. U., & Griffiths, R. E. 1997, ArXiv e-prints [arXiv:astro-ph/9704012]

Nicola, A., García-García, C., Alonso, D., et al. 2021, JCAP, 2021, 067

Olive, K. A., & Particle Data Group 2014, Chin. Phys. C, 38, 090001

Park, Y., & Rozo, E. 2020, MNRAS, 499, 4638

Peebles, P. J. E. 1973, ApJ, 185, 413

Planck Collaboration V. 2020, A&A, 641, A5

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

Raichoor, A., Mei, S., Erben, T., et al. 2014, ApJ, 797, 102

Reid, B., Ho, S., Padmanabhan, N., et al. 2016, MNRAS, 455, 1553

Riess, A. G., Macri, L., Casertano, S., et al. 2011, ApJ, 730, 119; Erratum: ApJ 732, 129 (2011)

Riess, A. G., Macri, L. M., Hoffmann, S. L., et al. 2016, ApJ, 826, 56

Riess, A. G., Yuan, W., Casertano, S., Macri, L. M., & Scolnic, D. 2020, ApJ, 896, L43

Ross, A. J., Samushia, L., Burden, A., et al. 2014, MNRAS, 437, 1109

Sánchez, A. G., Grieb, J. N., Salazar-Albornoz, S., et al. 2017a, MNRAS, 464, 1493

Sánchez, A. G., Scoccimarro, R., Crocce, M., et al. 2017b, MNRAS, 464, 1640

Schirmer, M. 2013, ApJS, 209, 21

Schlafly, E. F., & Finkbeiner, D. P. 2011, ApJ, 737, 103

Schlegel, D. J., Finkbeiner, D. P., & Davis, M. 1998, ApJ, 500, 525

Schneider, P., van Waaerbeke, L., Kilbinger, M., & Mellier, Y. 2002, A&A, 396, 1

Schulten, K., & Gordon, R. 1976, Comput. Phys. Commun., 11, 269

Scolnic, D. M., Jones, D. O., Rest, A., et al. 2018, ApJ, 859, 101

Secco, L. F., Samuroff, S., Krause, E., et al. 2022, Phys. Rev. D, 105, 023515

Seljak, U., & Zaldarriaga, M. 1997, Phys. Rev. Lett., 78, 2054

Sellentin, E., & Heavens, A. F. 2006, MNRAS, 456, L132

Spergel, D., Gehrels, N., Baltay, C., et al. 2015, ArXiv e-prints [arXiv:1503.03757]

Taylor, A. N., & Kitching, T. D. 2010, MNRAS, 408, 865

Taylor, P. L., Kitching, T. D., Alsing, J., et al. 2019, Phys. Rev. D, 100, 023519

Tröster, T., Sánchez, A. G., Asgari, M., et al. 2020, A&A, 633, L10

Tröster, T., Asgari, M., Blake, C., et al. 2021, A&A, 649, A88

Tröster, T., Mead, A. J., Heymans, C., et al. 2022, A&A, 660, A27

Troxel, M. A., MacCrann, N., Zuntz, J., et al. 2018, Phys. Rev. D, 98, 043528

Tutusaus, I., Martinelli, M., Cardone, V. F., et al. 2020, A&A, 643, A70

Upham, R. E., Whittaker, L., & Brown, M. L. 2020, MNRAS, 491, 3165

Upham, R. E., Brown, M. L., & Whittaker, L. 2021, MNRAS, 503, 1999

Upham, R. E., Brown, M. L., Whittaker, L., et al. 2022, A&A, 660, A114

van den Busch, J. L., Hildebrandt, H., Wright, A. H., et al. 2020, A&A, 642, A200

van Uitert, E., Joachimi, B., Joudaki, S., et al. 2018, MNRAS, 476, 4662

Verde, L., Treu, T., & Riess, A. G. 2019, Nat. Astron., 3, 891

Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nat. Methods, 17, 261

Wes McKinney, 2010, in Proceedings of the 9th Python in Science Conference, eds. S. van der Walt, & J. Millman, 56

Wright, A. H., Hildebrandt, H., Kuijken, K., et al. 2019, A&A, 632, A34

Wright, A. H., Hildebrandt, H., van den Busch, J. L., & Heymans, C. 2020, A&A, 637, A100

Xavier, H. S., Abdalla, F. B., & Joachimi, B. 2016, MNRAS, 459, 3693

¹ Department of Physics and Astronomy, University College London, Gower Street, London WC1E 6BT, UK
e-mail: arthur.loureiro@ucl.ac.uk

² Astrophysics Group, Blackett Laboratory, Imperial College London, SW7 2AZ, UK

³ Institute for Astronomy, University of Edinburgh, Royal Observatory, Blackford Hill, Edinburgh EH9 3HJ, UK

⁴ Jodrell Bank Centre for Astrophysics, Department of Physics and Astronomy, University of Manchester, Oxford Road, Manchester M13 9PL, UK

⁵ Mullard Space Science Laboratory, University College London, Holmbury St Mary, Dorking, Surrey RH5 6NT, UK

⁶ E. A. Milne Centre, University of Hull, Cottingham Road, Hull HU6 7RX, UK

⁷ Ruhr-Universität Bochum, Astronomisches Institut, German Centre for Cosmological Lensing (GCCL), Universitätsstr. 150, 44801 Bochum, Germany

⁸ Center for Theoretical Physics, Polish Academy of Sciences, al. Lotników 32/46, 02-668 Warsaw, Poland

⁹ Astronomisches Institut, Ruhr-Universität Bochum, Universitätsstr. 150, 44801 Bochum, Germany

¹⁰ Shanghai Astronomical Observatory (SHAO), Nandan Road 80, Shanghai 200030, PR China

¹¹ University of Chinese Academy of Sciences, Beijing 100049, PR China

¹² Institute of Cosmology and Gravitation, University of Portsmouth, Portsmouth PO1 3FX, UK

¹³ INAF-Osservatorio di Astrofisica e Scienza dello Spazio di Bologna, Via Piero Gobetti 93/3, 40129 Bologna, Italy

¹⁴ Max Planck Institute for Extraterrestrial Physics, Giessenbachstr. 1, 85748 Garching, Germany

¹⁵ INAF-Osservatorio Astrofisico di Torino, Via Osservatorio 20, 10025 Pino Torinese, TO, Italy

¹⁶ Department of Mathematics and Physics, Roma Tre University, Via della Vasca Navale 84, 00146 Rome, Italy

¹⁷ INFN-Sezione di Roma Tre, Via della Vasca Navale 84, 00146 Roma, Italy

¹⁸ INAF-Osservatorio Astronomico di Capodimonte, Via Moiariello 16, 80131 Napoli, Italy

¹⁹ INAF-IASF Milano, Via Alfonso Corti 12, 20133 Milano, Italy

²⁰ Institut de Física d'Altes Energies (IFAE), The Barcelona Institute of Science and Technology, Campus UAB, 08193 Bellaterra, Barcelona, Spain

²¹ Port d'Informació Científica, Campus UAB, C. Albareda s/n, 08193 Bellaterra, Barcelona, Spain

²² INAF-Osservatorio Astronomico di Roma, Via Frascati 33, 00078 Monteporzio Catone, Italy

²³ Department of Physics “E. Pancini”, University Federico II, Via Cinthia 6, 80126 Napoli, Italy

²⁴ INFN section of Naples, Via Cinthia 6, 80126 Napoli, Italy
25 Dipartimento di Fisica e Astronomia "Augusto Righi" - Alma Mater Studiorum Università di Bologna, Viale Berti Pichat 6/2, 40127 Bologna, Italy

26 INAF-Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5, 50125 Firenze, Italy

27 Centre National d'Etudes Spatiales, Toulouse, France

28 Institut national de physique nucléaire et de physique des particules, 3 rue Michel-Ange, 75794 Paris Cedex 16, France

29 European Space Agency/ESRIN, Largo Galileo Galilei 1, 00044 Frascati, Roma, Italy

30 ESAC/ESA, Camino Bajo del Castillo, s/n., Urb. Villafranca del Castillo, 28692 Villanueva de la Cañada, Madrid, Spain

31 Univ Lyon, Univ Claude Bernard Lyon 1, CNRS/IN2P3, IP2I Lyon, UMR 5822, 69622 Villeurbanne, France

32 Departamento de Física, Faculdade de Ciências, Universidade de Lisboa, Edifício C8, Campo Grande, 1749-016 Lisboa, Portugal

33 Instituto de Astrofísica e Ciências do Espaço, Faculdade de Ciências, Universidade de Lisboa, Campo Grande, 1749-016 Lisboa, Portugal

34 Université Paris-Saclay, CNRS, Institut d'astrophysique spatiale, 91405 Orsay, France

35 Department of Astronomy, University of Geneva, ch. d'Ecogia 16, 1290 Versoix, Switzerland

36 Department of Physics, Oxford University, Keble Road, Oxford OX1 3RH, UK

37 INFN-Padova, Via Marzolo 8, 35131 Padova, Italy

38 AIM, CEA, CNRS, Université Paris-Saclay, Université de Paris, 91191 Gif-sur-Yvette, France

39 Institut d'Estudis Espacials de Catalunya (IEEC), Carrer Gran Capità 2-4, 08034 Barcelona, Spain

40 Institute of Space Sciences (ICE, CSIC), Campus UAB, Carrer de Can Magrans, s/n, 08193 Barcelona, Spain

41 INAF-Osservatorio Astronomico di Trieste, Via G. B. Tiepolo 11, 34131 Trieste, Italy

42 Istituto Nazionale di Astrofisica (INAF) – Osservatorio di Astrofisica e Scienza dello Spazio (OAS), Via Gobetti 93/3, 40127 Bologna, Italy

43 Istituto Nazionale di Fisica Nucleare, Sezione di Bologna, Via Irne-rio 46, 40126 Bologna, Italy

44 INAF-Osservatorio Astronomico di Padova, Via dell'Osservatorio 5, 35122 Padova, Italy

45 Universitäts-Sternwarte München, Fakultät für Physik, Ludwig-Maximilians-Universität München, Scheinerstrasse 1, 81679 München, Germany

46 Institute of Theoretical Astrophysics, University of Oslo, PO Box 1029, Blindern 0315, Oslo, Norway

47 Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA 91109, USA

48 von Hoerner & Sulger GmbH, SchloßPlatz 8, 68723 Schwetzingen, Germany

49 Max-Planck-Institut für Astronomie, Königstuhl 17, 69117 Heidelberg, Germany

50 Aix-Marseille Univ, CNRS/IN2P3, CPPM, Marseille, France

51 Leiden Observatory, Leiden University, Niels Bohrweg 2, 2333 CA Leiden, The Netherlands

52 Université de Genève, Département de Physique Théorique and Centre for Astroparticle Physics, 24 quai Ernest-Ansermet, 1211 Genève 4, Switzerland

53 Department of Physics and Helsinki Institute of Physics, Gustaf Hällströminkatu 2, 00014 University of Helsinki, Finland

54 NOVA optical infrared instrumentation group at ASTRON, Oude Hoogeveensedijk 4, 7991 PD Dwingeloo, The Netherlands

55 Argelander-Institut für Astronomie, Universität Bonn, Auf dem Hügel 71, 53121 Bonn, Germany

56 INFN-Sezione di Bologna, Viale Berti Pichat 6/2, 40127 Bologna, Italy

57 Dipartimento di Fisica e Astronomia "Augusto Righi" - Alma Mater Studiorum Università di Bologna, via Piero Gobetti 93/2, 40129 Bologna, Italy

58 Centre for Extragalactic Astronomy, Department of Physics, Durham University, South Road, Durham DH1 3LE, UK

59 Institute for Computational Cosmology, Department of Physics, Durham University, South Road, Durham DH1 3LE, UK

60 Observatoire de Sauverny, Ecole Polytechnique Fédérale de Lausanne, 1290 Versoix, Switzerland

61 European Space Agency/ESTEC, Keplerlaan 1, 2201 AZ Noordwijk, The Netherlands

62 Department of Physics and Astronomy, University of Aarhus, Ny Munkegade 120, 8000 Aarhus C, Denmark

63 Institute of Space Science, Bucharest 077125, Romania

64 Max-Planck-Institut für Astrophysik, Karl-Schwarzschild Str. 1, 85741 Garching, Germany

65 Dipartimento di Fisica e Astronomia "G.Galilei", Università di Padova, Via Marzolo 8, 35131 Padova, Italy

66 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas (CIEMAT), Avenida Complutense 40, 28040 Madrid, Spain

67 Instituto de Astrofísica e Ciências do Espaço, Faculdade de Ciências, Universidade de Lisboa, Tapada da Ajuda, 1349-018 Lisboa, Portugal

68 Universidad Politécnica de Cartagena, Departamento de Electrónica y Tecnología de Computadoras, 30202 Cartagena, Spain

69 Kapteyn Astronomical Institute, University of Groningen, PO Box 800, 9700 AV Groningen, The Netherlands

70 Infrared Processing and Analysis Center, California Institute of Technology, Pasadena, CA 91125, USA

71 INAF-Osservatorio Astronomico di Brera, Via Brera 28, 20122 Milano, Italy

72 Dipartimento di Fisica, Università degli Studi di Torino, Via P. Giuria 1, 10125 Torino, Italy

73 INFN-Sezione di Torino, Via P. Giuria 1, 10125 Torino, Italy

74 INAF-IASF Bologna, Via Piero Gobetti 101, 40129 Bologna, Italy

75 Space Science Data Center, Italian Space Agency, via del Politecnico snc, 00133 Roma, Italy
## Appendix A: Extra figures and tables

We provide some extra figures and tables in this appendix to complement the relevant information regarding our analysis. Table A.1 shows details about the five redshift tomographic bins presented in Fig. 1 including the m-correction applied to the shear estimates for each tomographic bin. Following that, we present tables and figures related to the Bayesian inference results from Sect. 6.1. Table A.2 shows all the model parameters from our pseudo-C<sub>ℓ</sub> analysis (KiDS PCL) and its combination with clustering from SDSS baryon acoustic oscillations and redshift space distortions (KiDS PCL + SDSS) with the constrained parameters shown in bold. We also show both the maximum a posteriori results with the projected joint highest probability density (PJ-HPD) and the marginal constraints.

When examining a broader parameter sub-space, considering the dimensionless Hubble parameter (Fig. A.1), we can see

that two main factors play a crucial role in breaking the $\sigma_8 - \Omega_m$ degeneracy when combining SDSS and KiDS-1000. The BAO and RSD data help constraining significantly the $h - \Omega_b h^2$ plane, but there is a second factor playing a part here: the BBN-inspired prior on the baryonic matter density. As was highlighted in Sect. 5.3, the $5\sigma$ prior on $\Omega_b h^2$ based on BBN constraints from Olive & Particle Data Group (2014) was a conservative equivalent to what has been considered in Cuceu et al. (2019). However, it is a prior that has an impact when probing SDSS BAO and RSD data. This has a significant role in constraining the Hubble parameter, which impacts the $h - \Omega_b h^2$ plane and consequently the $h - \sigma_8$ plane, resulting in tighter constraints when combining clustering data with cosmic shear information. Here, however, we advocate that the inclusion of BBN data would be more constraining than our prior choice (Cuceu et al. 2019), therefore raising the point that such a choice was, from this perspective, conservative.

<table><caption>**Table A.1.** Details on each tomographic bin.</caption><thead><tr><th>Bin</th><th>z<sub>B</sub> range</th><th>mean z</th><th>Number of objects</th><th>n<sub>eff</sub>[arcmin<sup>-2</sup>]</th><th>m-correction</th></tr></thead><tbody><tr><td>1</td><td>0.1 &lt; z<sub>B</sub> ≤ 0.3</td><td>0.26</td><td>1 792 136</td><td>0.62</td><td>-0.009 ± 0.019</td></tr><tr><td>2</td><td>0.3 &lt; z<sub>B</sub> ≤ 0.5</td><td>0.40</td><td>3 681 319</td><td>1.18</td><td>-0.011 ± 0.020</td></tr><tr><td>3</td><td>0.5 &lt; z<sub>B</sub> ≤ 0.7</td><td>0.56</td><td>6 148 102</td><td>1.85</td><td>-0.015 ± 0.017</td></tr><tr><td>4</td><td>0.7 &lt; z<sub>B</sub> ≤ 0.9</td><td>0.79</td><td>4 544 395</td><td>1.26</td><td>0.002 ± 0.012</td></tr><tr><td>5</td><td>0.9 &lt; z<sub>B</sub> ≤ 1.2</td><td>0.98</td><td>5 096 059</td><td>1.31</td><td>0.007 ± 0.010</td></tr><tr><td>total</td><td>0.1 &lt; z<sub>B</sub> ≤ 1.2</td><td>N/A</td><td>21 262 011</td><td>6.95</td><td>N/A</td></tr></tbody></table>

**Notes.** The effective number of galaxies, $n_{\text{eff}}$, is defined in the same way as in Heymans et al. 2012. Details on the m-corrections estimation can be found in Giblin et al. (2021).

<table><caption>**Table A.2.** Summary of marginalised parameter constraints with 68% confidence intervals. Priors and parameter definitions are shown in Table 2.</caption><thead><tr><th>Parameter</th><th>KiDS-1000 PCL (PJ-HPD)</th><th>KiDS-1000 PCL (Marginal)</th><th>KiDS-1000 PCL + SDSS (PJ-HPD)</th><th>KiDS-1000 PCL + SDSS (Marginal)</th></tr></thead><tbody><tr><td>S<sub>8</sub></td><td><strong>0.754</strong><sup>+0.027</sup><sub>-0.029</sub></td><td><strong>0.742</strong><sup>+0.034</sup><sub>-0.023</sub></td><td><strong>0.771</strong><sup>+0.006</sup><sub>-0.032</sub></td><td><strong>0.762</strong><sup>+0.017</sup><sub>-0.023</sub></td></tr><tr><td>Ω<sub>c</sub>h<sup>2</sup></td><td>0.064<sup>+0.087</sup><sub>-0.013</sub></td><td>0.091<sup>+0.070</sup><sub>-0.034</sub></td><td><strong>0.115</strong><sup>+0.011</sup><sub>-0.008</sub></td><td><strong>0.114</strong><sup>+0.010</sup><sub>-0.009</sub></td></tr><tr><td>Ω<sub>b</sub>h<sup>2</sup></td><td>0.022<sup>+0.003</sup><sub>-0.001</sub></td><td>0.021<sup>+0.003</sup><sub>-0.001</sub></td><td>0.022<sup>+0.002</sup><sub>-0.003</sub></td><td>0.024<sup>+0.001</sup><sub>-0.003</sub></td></tr><tr><td>h</td><td>0.657<sup>+0.118</sup><sub>-0.003</sub></td><td>0.757<sup>+0.049</sup><sub>-0.059</sub></td><td><strong>0.664</strong><sup>+0.014</sup><sub>-0.022</sub></td><td><strong>0.674</strong><sup>+0.010</sup><sub>-0.025</sub></td></tr><tr><td>n<sub>s</sub></td><td>1.000<sup>+0.080</sup><sub>-0.074</sub></td><td>1.000<sup>+0.074</sup><sub>-0.080</sub></td><td>1.003<sup>+0.075</sup><sub>-0.073</sub></td><td>1.050<sup>+0.032</sup><sub>-0.118</sub></td></tr><tr><td>Ω<sub>m</sub></td><td>0.202<sup>+0.127</sup><sub>-0.053</sub></td><td>0.2149<sup>+0.125</sup><sub>-0.062</sub></td><td><strong>0.312</strong><sup>+0.011</sup><sub>-0.018</sub></td><td><strong>0.310</strong><sup>+0.015</sup><sub>-0.016</sub></td></tr><tr><td>σ<sub>8</sub></td><td>0.920<sup>+0.159</sup><sub>-0.215</sub></td><td>0.853<sup>+0.103</sup><sub>-0.226</sub></td><td><strong>0.756</strong><sup>+0.029</sup><sub>-0.044</sub></td><td><strong>0.750</strong><sup>+0.037</sup><sub>-0.039</sub></td></tr><tr><td>A<sub>1A</sub></td><td><strong>0.396</strong><sup>+0.251</sup><sub>-0.941</sub></td><td><strong>0.221</strong><sup>+0.639</sup><sub>-0.431</sub></td><td><strong>0.720</strong><sup>+0.231</sup><sub>-0.655</sub></td><td><strong>0.471</strong><sup>+0.483</sup><sub>-0.440</sub></td></tr><tr><td>A<sub>bary</sub></td><td>3.113<sup>+0.002</sup><sub>-0.663</sub></td><td>2.788<sup>+0.245</sup><sub>-0.422</sub></td><td>2.553<sup>+0.566</sup><sub>-0.103</sub></td><td>2.713<sup>+0.369</sup><sub>-0.333</sub></td></tr><tr><td>δz<sub>1</sub></td><td>0.002<sup>+0.010</sup><sub>-0.009</sub></td><td>0.002<sup>+0.011</sup><sub>+0.008</sub></td><td>0.003<sup>+0.011</sup><sub>-0.008</sub></td><td>0.005<sup>+0.012</sup><sub>-0.009</sub></td></tr><tr><td>δz<sub>2</sub></td><td>0.001<sup>+0.009</sup><sub>-0.010</sub></td><td>0.002<sup>+0.010</sup><sub>-0.010</sub></td><td>0.005<sup>+0.012</sup><sub>-0.009</sub></td><td>0.002<sup>+0.011</sup><sub>-0.011</sub></td></tr><tr><td>δz<sub>3</sub></td><td>-0.011<sup>+0.013</sup><sub>-0.005</sub></td><td>-0.016<sup>+0.012</sup><sub>-0.008</sub></td><td>-0.014<sup>+0.010</sup><sub>-0.010</sub></td><td>-0.016<sup>+0.012</sup><sub>-0.009</sub></td></tr><tr><td>δz<sub>4</sub></td><td>-0.013<sup>+0.008</sup><sub>-0.009</sub></td><td>-0.013<sup>+0.009</sup><sub>-0.008</sub></td><td>-0.015<sup>+0.008</sup><sub>-0.008</sub></td><td>-0.014<sup>+0.008</sup><sub>-0.008</sub></td></tr><tr><td>δz<sub>5</sub></td><td>0.006<sup>+0.009</sup><sub>-0.008</sub></td><td>0.009<sup>+0.007</sup><sub>-0.011</sub></td><td>0.002<sup>+0.006</sup><sub>-0.012</sub></td><td>0.007<sup>+0.009</sup><sub>-0.009</sub></td></tr></tbody></table>

**Notes.** Here we show the marginalised constraints for using the maximum a posteriori values with the projected joint highest posterior density (PJ-HPD) intervals and the usual one-dimensional marginalised constraints (Marginal). As in Asgari et al. (2021), we show the constrained parameters in bold, while parameters that are not constrained by the current data or data combination are mostly just reflecting their priors. Although *h* is constrained by the SDSS data, the value we find is a reflection of our BBN-inspired prior on $\Omega_b h^2$.
<figure><img src="image_14.png" /><figcaption>**Fig. A.1.** Marginalised 1D and 2D posterior distributions for 68% (darker) and 95% (lighter) C.I. in the $\Omega_m, \sigma_8, \Omega_b h^2$, and $h$ plane. The degeneracy broken by the introduction of clustering data to the cosmic shear constraints can be understood as better constraining power from the combination with clustering in the $h \times \Omega_b h^2$ plane, as well as the BBN-inspired prior. These yield Hubble constant values that are consistent with Planck 2018, in line with what was found by Cuceu et al. (2019) for constraints containing SDSS BAO and RSD information.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Fig. A.2.** Marginalised 1D and 2D posterior distributions for 68% (darker) and 95% (lighter) C.I. in the nuisance parameter space, comparing our pseudo-$C_\ell$ results (red) with the band-powers constraints (orange) from Asgari et al. (2021). The $D_{z_i}$ parameters are shown here as they are sampled, de-correlated displacements calculated from a linear transformation on $\delta_{z_i}$ using the Cholesky decomposition of the redshift covariance, $C_z$.</figcaption></figure>

In Fig. A.2 we show a comparison between pseudo-$C_\ell$ and band-power measurements (Asgari et al. 2021) for the sub-space of nuisance parameters, including the amplitude of intrinsic alignments ($A_{\text{IA}}$), the amplitude of baryonic feedback ($A_{\text{bary}}$), and the uncorrelated redshift displacements. In order to properly implement correlated Gaussian priors, we linearly transform

<figure><img src="image_16.png" /><figcaption>**Fig. A.3.** Marginalised 1D and 2D posterior distributions for 68% (darker) and 95% (lighter) C.I. in the cosmological parameter space, comparing the results from the SDSS BAO and RSD (blue) likelihoods we introduce in Sect. 5.2 with the BOSS DR12 likelihood used in the main KiDS-1000 analysis (grey, Sánchez et al. 2017b).</figcaption></figure>

the redshift displacements $\delta_{z_i}$ to $D_{z_i}$ using a Cholesky decomposition of the covariance of the redshift distributions. Therefore, the $D_{z_i}$ are uncorrelated. Naturally, this is transformed back into the original correlated space before using the redshift distributions in the forward-modelling described in Sect. 5.1. The results reported by Table A.2 are for $\delta_{z_i}$, that is the correlated redshifts. Figure A.2 shows good agreement between our analysis and band-powers from Asgari et al. (2021) for all nuisance parameters, except for the amplitude of intrinsic alignments, as discussed in detail in Sect. 6.1.

Finally, in Fig. A.3 we show a comparison between the SDSS BAO and RSD likelihood and data we present in Sect. 5.2 with the clustering likelihood used in Heymans et al. (2021) from BOSS DR12 (Sánchez et al. 2017b). The constraints are in agreement but differ due to containing different data sets and approaches. The SDSS likelihood we use extracts BAO and RSD information from BOSS and eBOSS LRGs, Ly-$\alpha$ auto-correlations and its cross-correlations with quasars. Cosmology is then fitted to BAO and RSD quantities. In contrast to this, the BOSS DR12 approach from Sánchez et al. (2017b) used in Heymans et al. (2021) and Tröster et al. (2020) fits cosmology directly to measurements of the anisotropic galaxy clustering correlation functions.

## Appendix B: B-modes systematics analysis

Here, we present the complete systematic analysis for the B-modes in the data. In the context of cosmic shear in a standard ΛCDM cosmology, the pure B-modes do not contain cosmological information. This means that any measurements of B-modes that are inconsistent with zero could be a potential red flag for systematic contamination. Partial sky observations cause E/B-mode leakage, where part of the power measured from potential B-modes leaks into the E-modes; such modes should be excluded from the analysis. This could be potentially problematic; however, from the scenario we are seeing in Fig. B.1, our B-modes have small correlations with systematics. Figure 4
<figure><img src="image_17.png" /><figcaption>**Fig. B.1.** P-values for a null detection from cross-correlations between the measured B-modes and different systematics. The cross-correlations between $B_1$ and PSF E-mode are shown in detail in Fig. B.2, while those between $B_2$ and stellar density are shown in Fig. B.3.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Fig. B.2.** Top panel: Cross-angular power spectra between PSF E-mode component and $E_5$ (purple circles) or $B_1$ (red triangles). Bottom panel: Signal-to-noise ratio (S/N) as defined by Eq. (B.1). The reduced $\chi^2$ for a null signal is 2.16 (p-value = 0.027) for $E_5$ and 2.51 (p-value = 0.010) for $B_1$. The S/N is below $5 \times 10^{-2}$ for all the bandwidths, meaning that the estimated data covariance is much larger than this correlation.</figcaption></figure>

shows that the pseudo-$C_ℓ$ B-modes are consistent with zero, indicating that any leakage from any significant correlations between B-modes and systematics are unlikely to interfere with the measured E-modes used for our cosmological analysis.

The highest correlations with B-modes manifest in the cross-correlations with the PSF E-mode and the stellar density from Gaia DR2, with a reduced $\chi^2 = 2.51$ and 2.60, respectively. For the PSF ellipticities, it is reasonable, however not guaranteed,

<figure><img src="image_19.png" /><figcaption>**Fig. B.3.** Cross-angular power spectra between stellar density and $E_2$ (dark circles) or $B_2$ (yellow triangles). For these cases the reduced $\chi^2$ for a null signal is 1.54 and 2.60, respectively.</figcaption></figure>

to assume that the contamination is linear in ellipticity space. Therefore, we can compare these with the errors measured from the mocks (see Sect. 3.4). The signal-to-noise ratio is defined as

$$S/N = \left( \frac{\tilde{C}_{\ell}^{X_i, \text{syst}}}{\sigma_{\ell}^{X_i, X_i}} \right), \qquad (B.1)$$

where $\sigma_{\ell}^{X_i, X_i}$ is the standard deviation for the auto-power spectra $X_i$, where $X = E, B$. Figure B.2 shows the cross-correlations between PSF E-mode and $B_1(E_5)$, as well as the S/N for these cross-correlations. From the S/N analysis, we can see that even with a relatively high $\chi^2_{\text{red}}$, the cross-correlations signal is sub-dominant for both $B_1 \times \text{PSF-E}$ and $E_5 \times \text{PSF-E}$ (more details in Sect. 4). For most bandwidths, the S/N is very close to zero, with the exception of one with S/N $\approx 5 \times 10^{-2}$.

For the case shown in Fig. B.3, we cannot perform a S/N analysis as the units between the $\tilde{C}_{\ell}^{B_2, \text{syst}}$ and the data covariance do not match. It is possible that we are observing a stellar contamination for the second tomographic bin B-mode. Nonetheless, this possible contamination has no impact in the E-modes used in our cosmological analysis. We know that even by considering E/B-mode leakage, the auto and cross angular power spectra with $B_2$ shown in Fig. 4 are consistent with zero. Another piece of evidence that this possible contamination does not affect our cosmological signal is the cross-correlation between $E_2$ and stellar overdensity, with a reduced $\chi^2 \approx 1.5$ (p-value = 0.13), also shown in Fig. B.3. Had this potential contamination leaked from the B-modes to the E-modes in the second tomographic bin, we would have observed a higher correlation between $E_2$ and stellar density, which does not seem to be the case. We conclude that any significant correlations between B-modes and systematics do not affect any of our cosmic shear data-vectors in a way that could potentially bias $S_8$ and other cosmological parameters.