© A. Loureiro et al. 2022

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 (Λ) 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, measuring some of these parameters with percent-level precision. The Planck Mission, however, did not contradict the standard model of cosmology. Instead, it established the Λ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 ≈ 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, σ₈, predicted by CMB observations, should match what we observe today with galaxy surveys. Several studies seem to suggest tensions between late-time probes and the CMB (Mörtsell & Dhawan 2018; Verde et al. 2019). These tensions arise in a few parameters, such as the Hubble constant, when constrained using Cepheids and type Ia Supernovae using distance ladder measurements (Riess et al. 2011, 2016, 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{{\mathrm{\Omega }}_{\mathrm{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₈ 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 ≳1.6σ 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σ 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₈ 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_ℓ 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_ℓ 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_ℓ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_ℓ 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_ℓ 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¹. 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_ℓ measurements. Our pseudo-C_ℓ 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_ℓ 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 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_ℓ 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² 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_ℓ 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_ℓ 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_ℓ estimator. Section 5 details our Bayesian inference: theory modelling for cosmic shear, forward-modelling of the mixing matrix, external galaxy clustering data from baryon acoustic oscillations (BAO) and redshift space distortions (RSD) used to complement the cosmic shear data, as well as the priors and the likelihood implemented in the analysis. Section 6 outlines our main results from the inferred posteriors with a discussion on tension with cosmic microwave background measurements. Section 7 summarises our findings and conclusions. Extra figures with cosmological constraints and nuisance parameters, as well as a complete table with all parameters probed in our analysis can be found in Appendix A. Finally, Appendix B discusses the details of the systematics null tests for B-modes specifically.

2. Cosmic shear data from KiDS-1000

Designed as a weak lensing experiment, the Kilo Degree Survey (Kuijken et al. 2015; de Jong et al. 2015, 2017) is a public survey by the European Southern Observatory, currently at its 4th Data Release (Kuijken et al. 2019)², with an observed area of around 1000 deg². Good seeing nights were prioritised to obtain high-resolution images in the r-band, resulting in an excellent average seeing of 0.7 arc-seconds. Like in previous analyses (Wright et al. 2019; Hildebrandt et al. 2020), KiDS DR4 photometry is processed with Theli (Erben et al. 2005; Schirmer 2013) and AstroWise (Begeman et al. 2013) using five-band NIR photometry (ZYJHK_s) from VIKING³ (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)⁴, covering an effective unmasked area of 777.4 deg² 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{{\mathrm{\Omega }}_{\mathrm{m}}/0.3}$ after the calibration is applied. We perform additional checks in the pseudo-C_ℓ 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 $\mathrm{\Delta }{z}_{\mathrm{B}}^{(1-4)}=0.2$, ranging from 0.1 < z_B ≤ 0.9, while the last tomographic bin ranges from 0.9 < z_B < 1.2, that is $\mathrm{\Delta }{z}_{\mathrm{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.

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, zB.

Cosmic shear, γ, can be estimated using the galaxy ellipticity measurements, corrected by the additive bias, ϵ^corr, obtained from the observed galaxy ellipticities, ϵ^obs. In the weak lensing regime (|γ|≪1), ϵ^corr is estimated by taking into account the additive bias c,

$$\begin{array}{c}\hfill {ϵ}^{\mathrm{corr}}={ϵ}^{\mathrm{obs}}-c,\end{array}$$(1)

where the additive bias, c = ⟨ϵ^obs⟩, is assumed to be the weighted average observed ellipticity over all galaxies in a given tomographic bin. We note that all these quantities (shear, ellipticity and additive bias) are complex quantities, where ϵ = ϵ₁ + iϵ₂.

3. Methodology

This section briefly discusses the pseudo-C_ℓ 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_ℓ 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 $\widehat{\mathbf{n}}$ is given as

$$\begin{array}{c}\hfill \gamma \left(\widehat{\mathbf{n}}\right)=\frac{1}{(1+m_z)}\frac{{\sum }_{i\in \widehat{\mathbf{n}}}{w}_i{ϵ}_i^{\mathrm{corr}}}{{\sum }_{i\in \widehat{\mathbf{n}}}{w}_i},\end{array}$$(2)

where the summation is taken over all galaxies in the pixel. For each galaxy, ${ϵ}_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, _±2Y_ℓm, is defined as

$$\begin{array}{cc}\hfill \gamma (\widehat{\mathbf{n}})& ={\gamma }_1(\widehat{\mathbf{n}})+\mathrm{i}{\gamma }_2(\widehat{\mathbf{n}})={\displaystyle \sum _{\ell m}(E_{\ell m}+\mathrm{i}{B}_{\ell m}){}_{+2}{Y}_{\ell m}(\widehat{\mathbf{n}});}\hfill \end{array}$$(3)

$$\begin{array}{cc}\hfill {\gamma }^{\ast }(\widehat{\mathbf{n}})& ={\gamma }_1(\widehat{\mathbf{n}})-\mathrm{i}{\gamma }_2(\widehat{\mathbf{n}})={\displaystyle \sum _{\ell m}(E_{\ell m}-\mathrm{i}{B}_{\ell m}){}_{-2}{Y}_{\ell m}(\widehat{\mathbf{n}}),}\hfill \end{array}$$(4)

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

$$\begin{array}{cc}& E_{\ell m}=\frac{1}{2}{\displaystyle \int }\mathrm{d}{\mathrm{\Omega }}_{\widehat{\mathbf{n}}}\{\gamma (\widehat{\mathit{n}}){}_{+2}{Y}_{\ell m}^{\ast }+{\gamma }^{\ast }(\widehat{\mathbf{n}}){}_{-2}{Y}_{\ell m}^{\ast }\};\hfill \end{array}$$(5)

$$\begin{array}{cc}\hfill & B_{\ell m}=\frac{-\mathrm{i}}{2}{\displaystyle \int }\mathrm{d}{\mathrm{\Omega }}_{\widehat{\mathbf{n}}}\{\gamma (\widehat{\mathit{n}}){}_{+2}{Y}_{\ell m}^{\ast }-{\gamma }^{\ast }(\widehat{\mathit{n}}){}_{-2}{Y}_{\ell m}^{\ast }\}.\hfill \end{array}$$(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_ℓm = 0 (Bartelmann & Schneider 2001; Hikage et al. 2011; Kilbinger 2015). Although for multiple lenses a B-mode power spectrum can be generated, its power is expected to be orders of magnitude smaller than the E-mode power spectrum, meaning that effectively all the cosmological information should come from the E-mode power spectra and that any significant signals in B-modes would arise from systematic contamination.

The situation is different for realistic cases when sky coverage is partial. The survey area is now constrained to a region $\mathcal{W}(\widehat{\mathbf{n}})$; where $\mathcal{W}(\widehat{\mathbf{n}})=0$ for regions that are not observed, masked out due to bright stars or bad seeing, or contain no galaxies; and 1 otherwise – where there are galaxies in the given pixel. In other words, $\mathcal{W}(\widehat{\mathbf{n}})$ is an effective binary mask constructed from the catalogue (which already excludes bad regions as detailed in Giblin et al. 2021). Although not the case for the KiDS-1000 dataset, future applications to Euclid data $\mathcal{W}(\widehat{\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{array}{cc}& {\stackrel{\sim }{\gamma }}_1(\widehat{\mathbf{n}})\pm \mathrm{i}{\stackrel{\sim }{\gamma }}_2(\widehat{\mathbf{n}\left\}\right\}})\equiv \mathcal{W}(\widehat{\mathbf{n}\left\}\right\}})[{\gamma }_1(\widehat{\mathbf{n}\left\}\right\}})\pm \mathrm{i}{\gamma }_2(\widehat{\mathbf{n}\left\}\right\}}))]\hfill \\ \hfill & ={\displaystyle \sum _{\ell m}({\stackrel{\sim }{E}}_{\ell m}\pm \mathrm{i}{\stackrel{\sim }{\mathit{B}}}_{\ell m}){}_{\pm 2}{Y}_{\ell m}(\widehat{\mathbf{n}}),}\hfill \end{array}$$(7)

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

$$\begin{array}{cc}& {\stackrel{\sim }{E}}_{\ell m}={\displaystyle \sum _{{\ell }^{\prime }{m}^{\prime }}(E_{\ell m}{\mathtt{W}}_{\ell {\ell }^{\prime }m{m}^{\prime }}^{+}+B_{\ell m}{\mathtt{W}}_{\ell {\ell }^{\prime }m{m}^{\prime }}^{-});}\hfill \end{array}$$(8)

$$\begin{array}{cc}\hfill & {\stackrel{\sim }{B}}_{\ell m}={\displaystyle \sum _{{\ell }^{\prime }{m}^{\prime }}(B_{\ell m}{\mathtt{W}}_{\ell {\ell }^{\prime }m{m}^{\prime }}^{+}-E_{\ell m}{\mathtt{W}}_{\ell {\ell }^{\prime }m{m}^{\prime }}^{-}),}\hfill \end{array}$$(9)

where the mixing kernels ${\mathtt{W}}_{\mathcal{l}{\mathcal{l}}^{\prime }m{m}^{\prime }}^{\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

$$\begin{array}{c}\hfill {\stackrel{\sim }{C}}_{\ell }^{E_i{E}_j}\approx \frac{1}{2\ell +1}{\displaystyle \sum _m{\stackrel{\sim }{E}}_{\ell m}^i{\stackrel{\sim }{E}}_{\ell m}^j-\langle {\stackrel{\sim }{\mathcal{N}}}_{\ell }^{i,j}\rangle {\delta }_{\mathit{ij}},}\end{array}$$(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 θ^R, while keeping the galaxy positions and weights fixed:

$$\begin{array}{cc}& {ϵ}_1^R={ϵ}_1^{\mathrm{corr}}cos({\theta }^R)-{ϵ}_2^{\mathrm{corr}}sin({\theta }^R),\hfill \end{array}$$(11)

$$\begin{array}{cc}\hfill & {ϵ}_2^R={ϵ}_2^{\mathrm{corr}}cos({\theta }^R)+{ϵ}_1^{\mathrm{corr}}sin({\theta }^R).\hfill \end{array}$$(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_ℓ. Using Eqs. (11) and (12) together with Eq. (2) to define the randomised shear estimates, γ^R, one can propagate this quantity through Eqs. (7) and (8) to obtain ${\stackrel{\sim }{E}}_{\ell m}^{R,i}$. Finally, the ${\stackrel{\sim }{\mathcal{N}}}_{\ell }$ term from Eq. (10) is simply defined as

$$\begin{array}{c}\hfill {\stackrel{\sim }{\mathcal{N}}}_{\ell }^{i,j}=\frac{1}{2\ell +1}{\displaystyle \sum _m{\stackrel{\sim }{E}}_{\ell m}^{R,i}{\stackrel{\sim }{E}}_{\ell m}^{R,j}.}\end{array}$$(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 ${\stackrel{\sim }{\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 }_{{ϵ}_{\text{obs}}}^2={\sigma }_{\gamma }^2+{\sigma }_{ϵ}^2$, where the shear field variance, ${\sigma }_{\gamma }^2$, is mixed with the variance from the intrinsic shape and measurement errors, ${\sigma }_{ϵ}^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 ${\mathtt{W}}_{\mathcal{l}{\mathcal{l}}^{\prime }m{m}^{\prime }}^{\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, $\mathcal{W}(\widehat{\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_ℓ 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_ℓℓ′, which contains information about the survey geometry in harmonic space (Peebles 1973; Brown et al. 2005; Hikage et al. 2011):

$$\begin{array}{c}\hfill \langle {\stackrel{\sim }{\mathit{C}}}_{\ell }^{\mathit{ij}}\rangle ={\displaystyle \sum _{{\ell }^{\prime }}{\mathtt{M}}_{\ell {\ell }^{\prime }}^{\mathit{ij}}{\mathit{C}}_{{\ell }^{\prime }}^{\mathit{ij}},}\end{array}$$(14)

or, assuming no parity violation modes are present $({\stackrel{\sim }{C}}_{\ell }^{E_i{B}_j}={\stackrel{\sim }{C}}_{\ell }^{B_i{E}_j})$,

$$\begin{array}{c}\hfill \left(\begin{array}{c}\langle {\stackrel{\sim }{C}}_{\ell }^{E_i{E}_j}\rangle \\ \langle {\stackrel{\sim }{C}}_{\ell }^{B_i{B}_j}\rangle \end{array}\right)={\displaystyle \sum _{{\ell }^{\prime }}}\left(\begin{array}{cc}{\mathtt{W}}_{\ell {\ell }^{\prime }}^{++}(i,j)& {\mathtt{W}}_{\ell {\ell }^{\prime }}^{--}(i,j)\\ {\mathtt{W}}_{\ell {\ell }^{\prime }}^{--}(i,j)& {\mathtt{W}}_{\ell {\ell }^{\prime }}^{++}(i,j)\\ \end{array}\right)\left(\begin{array}{c}{C}_{{\ell }^{\prime }}^{E_i{E}_j}\\ C_{{\ell }^{\prime }}^{B_i{B}_j}\end{array}\right).\end{array}$$(15)

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

$$\begin{array}{c}\hfill {\mathtt{W}}_{\ell {\ell }^{\prime }}^{\pm \pm }(i,j)=\frac{1}{2\ell +1}{\displaystyle \sum _{m{m}^{\prime }}{\mathtt{W}}_{\ell {\ell }^{\prime }m{m}^{\prime }}^{\pm ,i}{\left({\mathtt{W}}_{\ell {\ell }^{\prime }m{m}^{\prime }}^{\pm ,j}\right)}^{\ast },}\end{array}$$(16)

where ${\mathtt{W}}_{\mathcal{l}{\mathcal{l}}^{\prime }m{m}^{\prime }}^{\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{array}{cc}& {\mathtt{W}}_{\ell {\ell }^{\prime }m{m}^{\prime }}^{\pm ,i}\equiv {\displaystyle \int }\mathrm{d}{\mathrm{\Omega }}_{\widehat{\mathbf{n}}}{}_{\pm 2}{Y}_{{\ell }^{\prime }{m}^{\prime }}(\widehat{\mathbf{n}}){\mathcal{W}}^i(\widehat{\mathbf{n}}){}_{\pm 2}{Y}_{\ell m}^{\ast }(\widehat{\mathbf{n}});\hfill \\ \hfill & ={\displaystyle \sum _{{\ell }^{″}{m}^{″}}{(-1)}^m\sqrt{\frac{(2\ell +1)(2{\ell }^{\prime }+1)(2{\ell }^{″}+1)}{4\pi }}{\mathcal{W}}_{{\ell }^{″}{m}^{″}}^i}\hfill \\ \hfill & \times \left(\begin{array}{ccc}\ell & {\ell }^{\prime }& {\ell }^{″}\\ \pm 2& \mp 2& 0\end{array}\right)\left(\begin{array}{ccc}\ell & {\ell }^{\prime }& {\ell }^{″}\\ m& m^{\prime }& m^{″}\end{array}\right),\hfill \end{array}$$(17)

where the 2 × 3 objects above are the Wigner 3j symbols, calculated using UCLWig3j library⁵, and

$$\begin{array}{c}\hfill {\mathcal{W}}_{\ell m}^i={\displaystyle \oint }\mathrm{d}{\mathrm{\Omega }}_{\widehat{\mathbf{n}}}{\mathcal{W}}^i(\widehat{\mathbf{n}})Y_{\ell m}^{\ast }(\widehat{\mathbf{n}})\end{array}$$(18)

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

Fig. 2. KiDS-1000 full mixing matrix, Mℓℓ′ 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 ≤ ℓ ≤ 3070. From Eqs. (16) and (17) one can see that these objects are not symmetric in ℓ-ℓ′.

Here we would like to point out two important details regarding using a binary mask as opposed to a mask containing the shear weights or a variable depth mask. First, a mask containing shear weights, used as $\mathcal{W}(\widehat{\mathbf{n}})$, would cause an imprint from the clustering of galaxies into the mixing matrix, causing an artificial excess power for 𝒲_ℓ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_eff), it could be a problem for Stage-IV surveys such as Euclid. The pixelisation scheme smooths out the weights information if the resolution is too low for a denser survey. Meanwhile, increasing the resolution too much can cause many pixels containing a single galaxy only, which, given Eq. (2), would cause the weights information to be lost and bias the analysis. Details regarding how to implement the weights information without causing the aforementioned bias will be investigated in forthcoming work.

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

3.3. Pseudo-C_ℓ measurements

Our pseudo-C_ℓ 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_side = 1024 ($N_{\text{pix}}=12N_{\text{side}}^2$), measuring multipoles between 76 ≤ ℓ ≤ 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)⁶. Observing the binned mixing matrices shown in Fig. 3, we note that by cutting at an bellow an ℓ_min = 100, we find that just a small number of lower multipoles get mixed into the scales of interest. Therefore, we decreased the ℓ 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 ℓ_min = 76 was chosen arbitrarily. We were able to push to slightly lower ℓ values than the band-powers analysis from Asgari et al. (2021) as the latter suffers strong mixing bellow ℓ = 100. Figure 3 also gives us an intuition for the k-scales that enter our analysis as a function of redshift and multipole given the bandwidth ℓ-bins used.

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 = ℓ/χ(z), in units of h Mpc−1 for a given ℓ and z, where χ(z) is the co-moving distance. Bottom: the KiDS-1000 pseudo-Cℓs binned mixing matrices (see Fig. 2, for example). These are convolved with the theory Cℓs (see Eq. (23)) to model the effect of multipole mixing introduced by the survey mask. The filters are calculated using all the multipoles inside the eight log-spaced bandwidths from 76 ≤ ℓ ≤ 1500.

Following Brown et al. (2005)⁷, the measured pseudo-C_ℓs are then binned using eight log-spaced multipole bins between 76 ≤ ℓ ≤ 1500,

$$\begin{array}{c}\hfill {\stackrel{\sim }{C}}_L^{\mathit{ij}}=\frac{1}{2\pi }{\displaystyle \sum _{\ell }\frac{\ell (\ell +1)}{({\ell }^{L+1}-{\ell }^L)}{\stackrel{\sim }{C}}_{\ell }^{\mathit{ij}},}\end{array}$$(19)

where ℓ^L < ℓ ≤ ℓ^L + 1 are the boundaries of the bandwidth bins centred in L, shown in Fig. 3 as the vertical dotted lines. The measured pseudo-C_ℓ 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_ℓB-modes shown in the lower triangle part of Fig. 4 are consistent with zero with a ${\chi }_{\text{red}}^2=1.07$ for a null-detection, indicating that even with the mode-mixing introduced by the mixing matrix, the contribution of E-modes mixed into B-modes is not enough to bias the cosmological measurements. The auto-power spectrum for E₂ and combinations with it, mostly E₂ − E₃, 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).

Fig. 4. KiDS-1000 E-mode (upper triangle) and B-mode (lower triangle) pseudo-Cℓ 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 ≤ ℓ ≤ 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 χ2 ∼ 1 for a null detection.

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⁸ (Xavier et al. 2016) and Salmo⁹ presented in Joachimi et al. (2021b). This suite of mocks includes a number of known observational complexities such as variable depth source redshift distributions, and spatially varying galaxy sample properties. Flask + Salmo mocks like Egretta were used to validate the covariance matrix modelling in Joachimi et al. (2021b), as well as cosmic shear signal recovery, demonstrated in Heydenreich et al. (2020) using semi-analytic models. Contrary to N-Body simulations, log-normal simulations are much faster to obtain, allowing for quick survey simulations that accurately reproduce the desired 2pt functions in just a few CPU hours. Figure D.9 from Joachimi et al. (2021b) shows the agreement between the covariances obtained by Egretta mocks and the analytical covariances. Both methods agree to a few percent for cosmic shear 2pt functions, with the largest deviations coming from the first redshift bin. In other words, we are confident that our log-normal simulation-based covariance is accurate and that it reproduces, up to a few percent difference, what is expected by theoretical covariance calculations for cosmic shear analysis.

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_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, $\overline{z}$, is slightly higher than the SOM’s redshift distributions for the first two tomographic bins, the effective number of galaxies, n_eff, matches quite well. As higher mean redshifts lead to a higher variance in the sample’s pseudo-C_ℓs, the differences in the mean redshifts yield more conservative errors and, therefore, are not a concern.

The mocks were originally created at N_side = 4096; however, we apply the pseudo-C_ℓ estimator at N_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 ℓ-bins. The difference in resolution is taken into account by correcting for the pixel window function at the target N_side = 1024 (Górski et al. 2005). Lastly, we calculate the covariance of the noise subtracted pseudo-C_ℓ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_ℓ 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₈ 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_ℓ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_ℓ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)

$$\begin{array}{cc}& {\mathrm{Cov}}_m({\stackrel{\sim }{C}}_{\ell }^{i,j};{\stackrel{\sim }{C}}_{{\ell }^{\prime }}^{i^{\prime },j^{\prime }})={\stackrel{\sim }{C}}_{\ell }^{i,j}{\stackrel{\sim }{C}}_{{\ell }^{\prime }}^{i^{\prime },j^{\prime }}[{\sigma }_m^i{\sigma }_m^{i^{\prime }}+{\sigma }_m^i{\sigma }_m^{j^{\prime }}\hfill \\ \hfill & +{\sigma }_m^j{\sigma }_m^{i^{\prime }}+{\sigma }_m^j{\sigma }_m^{j^{\prime }}],\hfill \end{array}$$(20)

where ${\sigma }_m^i$ is defined as the standard deviation of the m-corrections shown in Table A.1 and ${\stackrel{\sim }{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 σ_m part of Eq. (20) is independent of ℓ, 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: Cov_tot = Cov_mocks + Cov_m. The resulting covariance, Cov_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.

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.

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 (ugriZYJHK_s): Galactic extinction for each of the observed nine bands derived from the E(B − V) maps from Schlegel et al. (1998) combined with the extinction coefficients, R_V = 3.1, from Schlafly & Finkbeiner (2011). The extinction coefficients for each filter, R_f, follow a linear relation between the observed bands: R_f = A_f/E(B − V). Therefore, we have only shown results for the r-band since this is the detection band for KiDS-1000 and other values are just linearly scaled based in this band. The values for R_f used for KiDS-1000 can be found in Table 4 of Kuijken et al. (2019).

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

• μ-Threshold: object detection threshold above the background.

• PSF ellipticities: mean PSF ellipticity components, treated as a spin-2 field, converted to PSF E/B-modes in harmonic space (see Eqs. (8) and (9)). The value quoted in Fig. 6 has the PSF α-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₁₁ + Q₂₂, where Q_ij is defined in terms of the second-order moments of the two-dimensional angular light distribution, $I(\widehat{\mathbf{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.

Fig. 6. P-values for the cross-correlations between the data’s tomographic pseudo E-modes and different systematics for a null detection (p ≤ 0.05 would indicate a detection). Here we show the systematics motivated in Sect. 4: object detection threshold (μ-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.

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}^{\mathrm{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

$$\begin{array}{c}\hfill {\stackrel{\sim }{C}}_{\ell }^{X_i,\mathrm{syst}}=\frac{1}{(2\ell +1)}{\displaystyle \sum _m{X}_{\ell m}^i{\left(a_{\ell m}^{\mathrm{syst}}\right)}^{\ast }.}\end{array}$$(21)

To estimate the errors for the cross-correlations between data and systematics, we obtain covariances in two different ways. For the catalogue quantities, we use a jackknife re-sampling method. We start by subdividing the survey footprint into 60 distinct regions using a K-Means algorithm on the sphere (Kwan et al. 2017)¹¹. While removing the kth region, we calculate the pseudo-C_ℓ (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_ℓ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 α parameter per tomographic bin, where α quantifies the average fraction of PSF ellipticity in the shear estimator. As the KiDS PSF has a very low ellipticity on average, and as α 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σ in the inferred S₈ constraints for KiDS-1000, inducing no significant biases.

In this analysis, we recover the angular dependency on this systematic using the α 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 × 10⁻². 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.

Fig. 7. PSF E-mode and E5 cross-power spectrum. Top panel: cross-angular power spectra between PSF E-mode component and E5 with the α-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 ∼ −2.5 × 10−2, meaning that it is well within the estimated covariance for our analysis.

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_ℓs.

5. Cosmological inference

This section describes the final elements necessary for cosmological inference using tomographic cosmic shear angular power spectra data: theory forward modelling, the likelihood and the priors we use to obtain the posteriors presented in 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-α forest data from eBOSS.

Cosmological inference is performed using a version of the Monte Python 3 software (Audren et al. 2013; Brinckmann & Lesgourgues 2019)¹² modified to directly sample from Gaussian priors¹³ 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₈ 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)¹⁴. Following the other KiDS-1000 analyses, we fix the sum of neutrino masses to ∑m_ν = 0.06eVc⁻², 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 ΛCDM cosmological parameters compared to more complex neutrino mass models.

The non-linear matter power spectrum, $P_{\mathrm{m}}^{\mathrm{nl}}(k)$, is calculated using the halo model from HMCode (Mead et al. 2015) included in Class¹⁵. This approach incorporates AGN baryonic feedback into the matter power spectrum using a model with two parameters: the halo bloating parameter, η₀, and the halo mass-concentration relation amplitude, A_bary. These parameters can be constrained by the following relation from Joudaki et al. (2018): η₀ = 0.98 − 0.12A_bary. Therefore, A_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_ℓs are a combination of contributions from cosmic shear, indexed by γ, and the intrinsic alignment of galaxies, indexed by I (Bernstein 2009; Joachimi & Bridle 2010):

$$\begin{array}{c}\hfill C_{\ell }^{{ϵ}_i,{ϵ}_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}.\end{array}$$(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 χ can be related as kf_K(χ)=ℓ+1/2, we can express the right-hand side C_ℓs from Eq. (22) as

$$\begin{array}{c}\hfill C_{\ell }^{X_i,Y_j}={\displaystyle {\int }_0^{{\chi }_{\mathrm{H}}}\mathrm{d}\chi \frac{W_X^i(\chi )W_Y^j(\chi )}{f_{\mathrm{K}}^2(\chi )}{P}_{\mathrm{m}}^{\mathrm{nl}}(\frac{\ell +1/2}{f_{\mathrm{K}}(\chi )},z(\chi )),}\end{array}$$(23)

where the X, Y indices are γ and/or I, the i, j indices are the tomographic bins indices for the different kernels, χ_H is the co-moving distance to the horizon and f_K(χ) is the Friedmann-Lemaître-Robertson-Walker (FLRW) co-moving angular diameter distance. Note, also, that we only consider flat ΛCDM models in this study, therefore, f_K(χ)=χ. The weak lensing kernel is defined as:

$$\begin{array}{c}\hfill W_{\gamma }^i(\chi )=\frac{3}{2}\frac{H_0^2{\mathrm{\Omega }}_{\mathrm{m}}}{c^{2}a(\chi )}{f}_{\mathrm{K}}(\chi ){\displaystyle {\int }_{\chi }^{{\chi }_{\mathrm{H}}}\mathrm{d}{\chi }^{\prime }{n}_s^i({\chi }^{\prime })\frac{f_{\mathrm{K}}({\chi }^{\prime }-\chi )}{f_{\mathrm{K}}({\chi }^{\prime })}.}\end{array}$$(24)

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

$$\begin{array}{c}\hfill W_{\mathrm{I}}^i(\chi )=-A_{\mathrm{IA}}\frac{C_1{\rho }_{\mathrm{cr}}{\mathrm{\Omega }}_{\mathrm{m}}}{D[a(\chi )]}{n}_s^i(\chi ),\end{array}$$(25)

with C₁ρ_cr ≈ 0.0134 (Joachimi et al. 2011) and where D[a(χ)] is the growth function and A_IA, the amplitude of intrinsic alignments, kept as a free parameter in the analysis but the same for all redshift tomographic bins. Although not physically motivated in the way non-linearities are accounted for, the NLA model was found by Fortuna et al. (2021) to be precise enough for the case we are currently analysing.

The theory angular power spectra from Eq. (22) are then corrected for the pixel window function (Górski et al. 2005), convolved with the mixing matrix using Eq. (14), assuming $C_{\mathcal{l}}^{B_i,B_j}=0$, and binned in eight log-space bandwidths between 76 ≤ ℓ ≤ 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 ℓ-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₈ parameter and it is limited by the degeneracy in the σ₈  − Ω_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¹⁶, which combine pre-reconstruction full-shape¹⁷ 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-α (Lyα) 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₈, 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_bary and A_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_eff)/r_d, D_H(z_eff)/r_d and the rate of structure growth, f(z_eff)σ₈(z_eff). Here, D_M(z) = (1 + z)⁻¹f_K[χ(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_eff) is the linear growth rate, and z_eff is the effective redshift of the measurement (as defined by Eq. (7), Ross et al. 2014).

For the Lyα 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_eff)/r_d and D_H(z_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 Λ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_eff)/r_d, D_H(z_eff)/r_d, f(z_eff)σ₈(z_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 Λ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_ℓ 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_ℓ likelihood, presented in Upham et al. (2020), as this approach is more computationally expensive. This is motivated by the fact that with a sufficiently high number of simulations and for sufficiently high multipoles, the Gaussian approximation is a satisfactory representation of the underlying distribution for the cosmological parameters of interest (Taylor et al. 2019; Lin et al. 2020; Upham et al. 2021).

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

As cosmic shear data are directly sensitive to the amplitude of density fluctuations, σ₈, or the growth of structure, S₈, 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¹⁰A_s), can lead to small biases as the prior becomes informative for S₈. Therefore, in our analysis, we sample S₈ 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 ±5σ from independent measurements such as the prior on h, taken from the Riess et al. (2016) measurements, and the prior on ω_c = Ω_ch², derived from Supernova Type Ia constraints on Ω_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, ω_b = Ω_bh². This prior is chosen to reflect ±5σ constraints from big bang nucleosynthesis (BBN) presented by Olive & Particle Data Group (2014). This is a conservative choice for cosmic shear-only analysis as justified by Heymans et al. (2021) and Asgari et al. (2021). When combining cosmic shear data with clustering information from BOSS/eBOSS galaxy clustering and the eBOSS DR16 Lyα 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 μ = (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 Λ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 TTTEEE+lowE¹⁸ (Planck Collaboration VI 2020; Planck Collaboration V 2020), the band-powers constraints for KiDS-1000 obtained by Asgari et al. (2021) and the 3 × 2pt analysis from Heymans et al. (2021). A summary of the main parameters constrained by cosmic shear data, σ₈, Ω_m, and S₈, using MAP+PJ-HPD constraints are presented in Table 3 for KiDS-1000 pseudo-C_ℓs and its combination with SDSS BAO and RDS data. A more complete summary of results, including the traditional marginal constraints, can be found in Table 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, σ₈, and the matter density parameters, Ω_m. Results between the KiDS-1000 estimators are consistent, with the pseudo-C_ℓ results (in red) demonstrating a slightly lower value for σ₈ 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 × 2pt 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 ΛCDM, that is it does not use the full galaxy power spectrum, yielding slightly larger constraints.

Fig. 8. Marginalised 2D posterior distributions for 68% (darker) and 95% (lighter) C.I. in the Ωm − σ8 plane. Here we exhibit the pseudo-Cℓ 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ℓ measurements (purple) have a completely different degeneracy than the Planck 2018 constraints from the TTTEEE-lowE combination (green). We also compare it with the 3 × 2pt analysis from Heymans et al. (2021) (black-dashed).

The clustering posteriors we obtain demonstrate a slightly higher value of σ₈ than the BOSS constraints presented in Heymans et al. (2021), which is more apparent when considering the $S_8={\sigma }_8\sqrt{{\mathrm{\Omega }}_{\mathrm{m}}/0.3}$ parameter instead, as shown in Fig. A.3. Here, we find $S_8^{\mathrm{SDSS}}=0.{812}_{-0.049}^{+0.050}$ for the SDSS BAO and RSD data set which is in agreement with both KiDS-1000 pseudo-C_ℓ, $S_8^{\mathrm{PCL}}=0.{754}_{-0.029}^{+0.027}$, and the Planck 2018 results. The two main reasons why our clustering constraints differ from those in Heymans et al. (2021) are the addition of high-redshift eBOSS LRGs and the Ly-α 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 Ω_m parameter. This demonstrates the S₈ measurement robustness from KiDS-1000 cosmic shear analyses. The S₈ 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% 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.

Fig. 9. Same as Fig. 8, but for the marginalised Ωm and S8 constraints.

A summary of S₈ 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_ℓ 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)¹⁹ 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₈, 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.

Fig. 10. Comparison between different cosmological probe constraints on the growth of structure parameter, S8 = σ8(Ω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ℓ 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σ.

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

$$\begin{array}{cc}\hfill [\mathrm{PCL}]A_{\mathrm{IA}}& =0.{396}_{-0.931}^{+0.251};\hfill \end{array}$$(26)

$$\begin{array}{cc}\hfill [2\mathrm{PCF}]A_{\mathrm{IA}}& =0.{387}_{-0.374}^{+0.321};\hfill \end{array}$$(27)

$$\begin{array}{cc}\hfill [\mathrm{COSEBIs}]A_{\mathrm{IA}}& =0.{264}_{-0.337}^{+0.424};\hfill \end{array}$$(28)

$$\begin{array}{cc}\hfill [\mathrm{BP}]A_{\mathrm{IA}}& =0.{973}_{-0.383}^{+0.292}.\hfill \end{array}$$(29)

We speculate that this could be the reason why our method finds a slightly lower peak for the matter density parameter, Ω_m, as a similar effect appears in Asgari et al. (2021). Despite this, Ω_m is not strongly constrained in our analysis, meaning that the uncertainty is large enough that the constraints remain consistent for both Ω_m and A_IA, as investigated by Asgari et al. (2021). Since both band-powers and pseudo-C_ℓs are in harmonic space, one could expect them to match better than pseudo-C_ℓ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 ℓ = 100. Meanwhile, correlation functions and COSEBIs have filters in harmonic space that capture power below ℓ = 100. Comparing this with the pseudo-C_ℓs ℓ_min = 76, Fig. 3 shows that larger scales also leak through the mixing matrix filters into our analysis. This slight difference in ℓ range potentially explains the differences between band-powers and pseudo-C_ℓ, 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_bary and A_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 χ² for the cosmic shear-only analysis is 1.18, reflecting an acceptable p-value of 0.06.

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={\sigma }_8\sqrt{{\mathrm{\Omega }}_{\mathrm{m}}/0.3}$), σ8 as a derived parameter, the physical baryonic matter density (Ωbh2), the dimensionless Hubble constant (h), and the spectral index (ns). Included are also two nuisance parameters: the amplitude of intrinsic alignments (AIA) and the baryonic feedback amplitude (Abary) as these correlate with some cosmological parameters. When combining the KiDS-1000 pseudo-Cℓ constraints (red) with the SDSS BAO and RSD constraints (blue), we obtain constraints (purple) that break a few degeneracies that cosmic shear alone cannot. For comparison, we added the Planck 2018 TTTEEE-lowE (Planck Collaboration V 2020) constrains (green) and the KiDS-1000 band-powers constraints (orange) from Asgari et al. (2021).

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₈. By assuming Gaussianity and adopting the marginal constraints to quantify tension, we initially employ the conventional marginal tension estimate

$$\begin{array}{c}\hfill \tau =\frac{|⟨S_8^{\mathrm{A}}⟩-⟨S_8^{\mathrm{B}}⟩|}{\sqrt{{\sigma }^2(S_8^{\mathrm{A}})+{\sigma }^2(S_8^{\mathrm{B}})}},\end{array}$$(30)

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

Comparisons between KiDS-1000 pseudo-C_ℓ measurements and Planck 2018 TTTEEE+lowE nominal results (Planck Collaboration VI 2020; Planck Collaboration V 2020) lead to a 2.8σ tension between late- and early-Universe constraints. When combining the pseudo-C_ℓ measurements with SDSS BAO and RSD data, the tension with Planck 2018 increases to 2.9σ, 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_ℓ 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 τ = 1.4σ when combining cosmic shear and clustering data. However, ACT does not measure the TT power spectrum below ℓ = 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₈ 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σ 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_ℓ alone, we find 3.2σ tension and combining cosmic shear with clustering, we find 3.3σ tension. However, using the suspiciousness statistic (Lemos et al. 2020), Handley & Lemos (2021) pointed out that Planck 2018 and ACT have a 2.6σ tension between each other.

It is clear from Fig. 11 and Table 3, that the marginal distribution of S₈ is not Gaussian. To further assess the tension between our KiDS-1000 pseudo-C_ℓ 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(θ) and q(θ), d_H is defined as

$$\begin{array}{c}\hfill d_{\mathrm{H}}^2[p,q]=\frac{1}{2}{\displaystyle \int }\mathrm{d}\theta {[\sqrt{p(\theta )}-\sqrt{q(\theta })]}^2.\end{array}$$(31)

With the metric above, we measured d_H = 0.95 between our KiDS-1000 pseudo-C_ℓ and Planck Legacy’s S₈ constraints. Using the methodology described by Heymans et al. (2021), this can be expressed as a ∼3.1σ tension. For our combination of cosmic shear from KiDS and clustering from SDSS, we obtain d_H = 0.94 – corresponding to a ∼3σ 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_ℓ) analysis of the current state-of-the-art galaxy shape measurements from the Kilo Degree Survey’s fourth data release (KiDS-1000, Kuijken et al. 2019), obtaining cosmological constraints for the current concordance model, the flat ΛCDM cosmology. The KiDS-1000 shear catalogue (Giblin et al. 2021) is now publicly available²⁰, 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_ℓ 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_ℓ estimator performed as well as other 2pt functions, with no flat-sky approximation – an important and necessary characteristic that estimators of the next generation of galaxy surveys need. The forward-modelling approach allowed us to maintain accuracy and avoid numerical instabilities from inverting and deconvolving a mixing matrix from our estimates with no significant computational cost when calculating our theory-vector for cosmological inference.

Through this re-analysis in harmonic space of KiDS-1000 data as a road-test of our prototype estimator for the Euclid Science Ground Segment we find that, with KiDS data, we are operating at the limit of this methodology. First, noise realisations of the galaxy shape catalogue will not be feasible for the size of the expected Euclid shape catalogue (∼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_ℓ 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_ℓs contained only cosmological and astrophysical signals and could be used for Bayesian inference of the relevant flat ΛCDM parameters.

Cosmic shear alone constrains the growth of structure parameter, S₈, 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_ℓ 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-α 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_ℓ, 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 ℓ_min = 100, while we used ℓ_min = 76 for pseudo-C_ℓ, 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σ (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_ℓ analysis alone yielded $S_8=0.{754}_{-0.029}^{+0.027}$, or 9.8  ±  3.3% lower than the S₈ 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.032}^{+0.006}$, reflecting a 7.8 ± 2.3% lower value when comparing the marginal constraints with Planck’s measurement of S₈. 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σ 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σ, 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σ. 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σ. 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.

Acknowledgments

The authors would like to thank Dr Lorne Whiteway for useful comments and discussions. We would also like to thank the reviewer for their helpful comments that helped enrich this work. The mixing matrix calculation used a publicly available Wigner 3-j symbol code developed by Dr Whiteway and it can be found in https://github.com/LorneWhiteway/UCLWig3j. The K-Means code used for the jackknife covariances for systematics was built on top of a code developed by Dr Erin Sheldon, publicly available in https://github.com/esheldon/kmeans_radec. The analysis performed in this work also made use of the following packages and software: NumPy (Harris et al. 2020), SciPy (Virtanen et al. 2020), Matplotlib (Hunter 2007), GetDist (Lewis 2019), Pandas (Wes McKinney 2010), AstroPy (Astropy Collaboration 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 Economia y Competitividad, the National Aeronautics and Space Administration, the National Astronomical Observatory of Japan, the Netherlandse Onderzoekschool Voor Astronomie, the Norwegian Space Agency, the Romanian Space Agency, the State Secretariat for Education, Research and Innovation (SERI) at the Swiss Space Office (SSO), and the United Kingdom Space Agency. A complete and detailed list is available on the Euclid web site (http://www.euclid-ec.org). Author Contributions: All authors contributed to the development and writing of this paper. The authorship list is given in several groups: the lead authors (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

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_ℓ 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 σ₈ − Ω_m degeneracy when combining SDSS and KiDS-1000. The BAO and RSD data help constraining significantly the h − Ω_bh² 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σ prior on Ω_bh² 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 − Ω_bh² plane and consequentially the h − σ₈ 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.

Fig. A.1. Marginalised 1D and 2D posterior distributions for 68% (darker) and 95% (lighter) C.I. in the Ωm, σ8, Ωbh2, 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 × Ωbh2 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.

In Fig. A.2 we show a comparison between pseudo-C_ℓ and band-power measurements (Asgari et al. 2021) for the sub-space of nuisance parameters, including the amplitude of intrinsic alignments (A_IA), the amplitude of baryonic feedback (A_bary), and the uncorrelated redshift displacements. In order to properly implement correlated Gaussian priors, we linearly transform the redshift displacements δ_zi to D_zi using a Cholesky decomposition of the covariance of the redshift distributions. Therefore, the D_zi 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 δ_zi, 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.

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ℓ results (red) with the band-powers constraints (orange) from Asgari et al. (2021). The Dzi parameters are shown here as they are sampled, de-correlated displacements calculated from a linear transformation on δzi using the Cholesky decomposition of the redshift covariance, Cz.

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

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

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

Fig. B.1. P-values for a null detection from cross-correlations between the measured B-modes and different systematics. The cross-correlations between B1 and PSF E-mode are shown in detail in Fig. B.2, while those between B2 and stellar density are shown in Fig. B.3.

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 χ² = 2.51 and 2.60, respectively. For the PSF ellipticities, it is reasonable, however not guaranteed, 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

$$\begin{array}{c}\hfill \text{S/N}=\left(\frac{{\stackrel{\sim }{C}}_{\ell }^{{\text{X}}_i,\text{syst}}}{{\sigma }_{\ell }^{{\text{X}}_i,{\text{X}}_i}}\right),\end{array}$$(B.1)

where ${\sigma }_{\ell }^{{\text{X}}_i,{\text{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₁(E₅), 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 }_{\text{red}}^2$, the cross-correlations signal is subdominant for both B₁ × PSF-E and E₅ × 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 ≈ 5 × 10⁻².

Fig. B.2. Top panel: Cross-angular power spectra between PSF E-mode component and E5 (purple circles) or B1 (red triangles). Bottom panel: Signal-to-noise ratio (S/N) as defined by Eq. (B.1). The reduced χ2 for a null signal is 2.16 (p-value = 0.027) for E5 and 2.51 (p-value = 0.010) for B1. The S/N is below 5 × 10−2 for all the bandwidths, meaning that the estimated data covariance is much larger than this correlation.

For the case shown in Fig. B.3, we cannot perform a S/N analysis as the units between the ${\stackrel{\sim }{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₂ 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₂ and stellar overdensity, with a reduced χ² ≈ 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₂ 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₈ and other cosmological parameters.

Fig. B.3. Cross-angular power spectra between stellar density and E2 (dark circles) or B2 (yellow triangles). For these cases the reduced χ2 for a null signal is 1.54 and 2.60, respectively.

All Tables

All Figures

	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, zB.
In the text

Fig. 2. KiDS-1000 full mixing matrix, Mℓℓ′ 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 ≤ ℓ ≤ 3070. From Eqs. (16) and (17) one can see that these objects are not symmetric in ℓ-ℓ′.

In the text

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 = ℓ/χ(z), in units of h Mpc−1 for a given ℓ and z, where χ(z) is the co-moving distance. Bottom: the KiDS-1000 pseudo-Cℓs binned mixing matrices (see Fig. 2, for example). These are convolved with the theory Cℓs (see Eq. (23)) to model the effect of multipole mixing introduced by the survey mask. The filters are calculated using all the multipoles inside the eight log-spaced bandwidths from 76 ≤ ℓ ≤ 1500.

In the text

Fig. 4. KiDS-1000 E-mode (upper triangle) and B-mode (lower triangle) pseudo-Cℓ 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 ≤ ℓ ≤ 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 χ2 ∼ 1 for a null detection.

In the text

	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.
In the text

Fig. 6. P-values for the cross-correlations between the data’s tomographic pseudo E-modes and different systematics for a null detection (p ≤ 0.05 would indicate a detection). Here we show the systematics motivated in Sect. 4: object detection threshold (μ-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.

In the text

Fig. 7. PSF E-mode and E5 cross-power spectrum. Top panel: cross-angular power spectra between PSF E-mode component and E5 with the α-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 ∼ −2.5 × 10−2, meaning that it is well within the estimated covariance for our analysis.

In the text

Fig. 8. Marginalised 2D posterior distributions for 68% (darker) and 95% (lighter) C.I. in the Ωm − σ8 plane. Here we exhibit the pseudo-Cℓ 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ℓ measurements (purple) have a completely different degeneracy than the Planck 2018 constraints from the TTTEEE-lowE combination (green). We also compare it with the 3 × 2pt analysis from Heymans et al. (2021) (black-dashed).

In the text

	Fig. 9. Same as Fig. 8, but for the marginalised Ωm and S8 constraints.
In the text

Fig. 10. Comparison between different cosmological probe constraints on the growth of structure parameter, S8 = σ8(Ω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ℓ 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σ.

In the text

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={\sigma }_8\sqrt{{\mathrm{\Omega }}_{\mathrm{m}}/0.3}$), σ8 as a derived parameter, the physical baryonic matter density (Ωbh2), the dimensionless Hubble constant (h), and the spectral index (ns). Included are also two nuisance parameters: the amplitude of intrinsic alignments (AIA) and the baryonic feedback amplitude (Abary) as these correlate with some cosmological parameters. When combining the KiDS-1000 pseudo-Cℓ constraints (red) with the SDSS BAO and RSD constraints (blue), we obtain constraints (purple) that break a few degeneracies that cosmic shear alone cannot. For comparison, we added the Planck 2018 TTTEEE-lowE (Planck Collaboration V 2020) constrains (green) and the KiDS-1000 band-powers constraints (orange) from Asgari et al. (2021).

In the text

Fig. A.1. Marginalised 1D and 2D posterior distributions for 68% (darker) and 95% (lighter) C.I. in the Ωm, σ8, Ωbh2, 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 × Ωbh2 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.

In the text

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ℓ results (red) with the band-powers constraints (orange) from Asgari et al. (2021). The Dzi parameters are shown here as they are sampled, de-correlated displacements calculated from a linear transformation on δzi using the Cholesky decomposition of the redshift covariance, Cz.

In the text

	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).
In the text

	Fig. B.1. P-values for a null detection from cross-correlations between the measured B-modes and different systematics. The cross-correlations between B1 and PSF E-mode are shown in detail in Fig. B.2, while those between B2 and stellar density are shown in Fig. B.3.
In the text

Fig. B.2. Top panel: Cross-angular power spectra between PSF E-mode component and E5 (purple circles) or B1 (red triangles). Bottom panel: Signal-to-noise ratio (S/N) as defined by Eq. (B.1). The reduced χ2 for a null signal is 2.16 (p-value = 0.027) for E5 and 2.51 (p-value = 0.010) for B1. The S/N is below 5 × 10−2 for all the bandwidths, meaning that the estimated data covariance is much larger than this correlation.

In the text

	Fig. B.3. Cross-angular power spectra between stellar density and E2 (dark circles) or B2 (yellow triangles). For these cases the reduced χ2 for a null signal is 1.54 and 2.60, respectively.
In the text