KiDS and Euclid : Cosmological implications of a pseudo angular power spectrum analysis of KiDS-1000 cosmic shear tomography (cid:63)

We present a tomographic weak lensing analysis of the Kilo Degree Survey Data Release 4 (KiDS-1000), using a new pseudo angular power spectrum estimator (pseudo-C (cid:96) ) under development for the ESA Euclid mission. Over 21 million galaxies with shape information are divided into ﬁve tomographic redshift bins, ranging from 0.1 to 1.2 in photometric redshift. We measured pseudo-C (cid:96) using eight bands in the multipole range 76 < (cid:96) < 1500 for auto-and cross-power spectra between the tomographic bins. A series of tests were carried out to check for systematic contamination from a variety of observational sources including stellar number density, variations in survey depth, and point spread function properties. While some marginal correlations with these systematic tracers were observed, there is no evidence of bias in the cosmological inference. B -mode power spectra are consistent with zero signal, with no signiﬁcant residual contamination from E / B -mode leakage. We performed a Bayesian analysis of the pseudo-C (cid:96) estimates by forward modelling the e ﬀ ects of the mask. Assuming a spatially ﬂat Λ CDM cosmology, we constrained the structure growth parameter S 8 = σ 8 ( Ω m / 0 . 3) 1 / 2 = 0 . 754 + 0 . 027 − 0 . 029 . When combining cosmic shear from KiDS-1000 with baryon acoustic oscillation and redshift space distortion data from recent Sloan Digital Sky Survey (SDSS) measurements of luminous red galaxies, as well as the Lyman-α forest and its cross-correlation with quasars, we tightened these constraints to S 8 = 0 . 771 + 0 . 006 − 0 . 032 . These results are in very good agreement with previous KiDS-1000 and SDSS analyses and conﬁrm a ∼ 3 σ tension with early-Universe constraints from cosmic microwave background experiments.


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, mea-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(Riess et al. , 2016(Riess et al. , 2020;;Bernal et al. 2016;Di Valentino et al. 2016;Lin & Ishak 2017;Efstathiou 2020), and the structure growth parameter, S 8 = σ 8 √ Ω m /0.3, when measured using cosmic shear surveys (Heymans et al. 2013(Heymans et al. , 2021;;MacCrann et al. 2015;Joudaki et al. 2017Joudaki et al. , 2020;;Hildebrandt et al. 2017;Abbott et al. 2018;Park & Rozo 2020;Hikage et al. 2019;Lemos et al. 2021;Asgari et al. 2020Asgari et al. , 2021;;Tröster et al. 2021;Amon et al. 2022;Secco et al. 2022;DES Collaboration 2022).
Cosmic shear is the study of the weak gravitational lensing of background galaxies due to the intervening large-scale structure of the Universe.This lensing effect produces correlations in galaxy shapes, and these correlations can be used to determine how the dark and baryonic foreground matter is distributed (see Bartelmann & Schneider 2001;Kilbinger 2015 for comprehensive reviews of cosmic shear).Moreover, by measuring the cosmic shear signal in tomographic redshift bins, we can study how the distribution of matter evolves over time; this enables us to place tight constraints on the structure formation of the Universe and its evolution with redshift (Hall 2021), as well as the dark energy equation-of-state and other standard cosmological parameters.
As more data are independently collected by current galaxy surveys, the slight tension in S 8 between CMB and cosmic shear surveys does not seem to disappear.More specifically, all results coming from the Dark Energy Survey (DES, Abbott et al. 2018;Troxel et al. 2018) and the Hyper Supreme Camera (HSC, Hikage et al. 2019) indicate a lower value for the amplitude of matter density fluctuations, consistent with the Kilo Degree Survey (KiDS, Hildebrandt et al. 2020;Heymans et al. 2021), but also consistent with the Planck projections at a 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 8 with cosmic shear probes when combined with other CMB experiments (Aiola et al. 2020;Han et al. 2021).This also suggests that it is unlikely that the tension we currently observe comes from a systematic contamination in the CMB experiments.If this discrepancy between measurements is indeed induced by systematic contamination, it would have to be shared across several independent datasets with widely different properties and analysis choices.Hopes are that experiments planned for the near future will shine a light on this interesting cosmic puzzle as there is still room to improve the precision on cosmic shear experiments.
Future galaxy surveys such as the Legacy Survey of Space and Time (LSST), carried out at the Vera C. Rubin Observatory (LSST Dark Energy Science Collaboration 2012), the European Space Agency's Euclid Mission (Laureijs et al. 2011), and the National Aeronautics and Space Administration's (NASA) Nancy Grace Roman Space Telescope (Spergel et al. 2015) will strongly rely on cosmic shear as one of their primary probes to understand the large-scale structure of the Universe.These future surveys will obtain cosmological information from billions of galaxy shapes using a plethora of statistical tools including correlation functions, angular power spectra, mass maps, and bi-spectra.Amongst these, the pseudo-C 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 1 .In the presence of a mask, the mask's characteristic scales, from its shapes and holes, will mix nearby modes in harmonic space causing correlations between multipoles to appear.For spin-2 fields, such as the shear field, partial sky observations can also cause ambiguity between E/B-modes.The multipole mixing, however, can be fully accounted for by analytically calculating the mixing matrix from the survey's geometry.
In this work, our approach consists of forward-modelling the mode mixing effect due to masking in the theory.Taking this approach, we are least prone to numerical instability from inverting and deconvolving the mixing matrix from pseudo-C 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 forwardmodelling 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 andTrö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.51Million 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-ofthe-art data from a Stage-III cosmic shear experiment: the public ESO Kilo Degree Survey (KiDS).KiDS is currently at its fourth Data Release (KiDS DR4, Kuijken et al. 2019), having observed around 1000 deg2 of the sky with nine-band matched photometry.Giblin et al. (2021) outlined the KiDS-1000 weak lensing catalogue production and validation, with redshift calibration and measurements presented in Hildebrandt et al. (2021).Using three different 2pt statistics, Asgari et al. (2021) performed a cosmic shear analysis of the KiDS-1000 catalogue, including a series of consistency checks.The main KiDS-1000 analysis is finally presented in Heymans et al. (2021), combining it with galaxy clustering from BOSS DR12 (Reid et al. 2016;Sánchez et al. 2017a), to perform a joint (3 × 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. 2017Alam et al. , 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.

Cosmic shear data from KiDS-1000
Designed as a weak lensing experiment, the Kilo Degree Survey (Kuijken et al. 2015;de Jong et al. 2015de Jong et al. , 2017) is a public sur-vey by the European Southern Observatory, currently at its 4th Data Release (Kuijken et al. 2019) 2 , with an observed area of around 1000 deg 2 .Good seeing nights were prioritised to obtain high-resolution images in the r-band, resulting in an excellent average seeing of 0.7 arc-seconds.Like in previous analyses (Wright et al. 2019;Hildebrandt et al. 2020), KiDS DR4 photometry is processed with Theli (Erben et al. 2005;Schirmer 2013) and AstroWise (Begeman et al. 2013) using five-band NIR photometry (ZY JHK s ) from VIKING 3 (Edge et al. 2013) combined with four-band optical photometry (ugri) from KiDS to estimate photometric redshifts in an accurate and precise way.In this section, we briefly describe the most important and relevant details of the catalogue used in our study.
We perform our main analysis using the latest cosmic shear catalogue from the 4th Data Release (KiDS-1000) 4 , covering an effective unmasked area of 777.4 deg 2 and containing a total of 21 262 011 lensed galaxies.Detailed descriptions of the catalogue's construction, including the LensFit (Miller et al. 2007(Miller et al. , 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 = σ 8 √ Ω 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 ∆z (1−4) B = 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 ∆z (5) B = 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 andHildebrandt et al. 2021 for details).Table A.1 contains more details about the samples' redshift tomographic bins and their multiplicative shear calibration corrections.
Cosmic shear, γ, 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, where the additive bias, c = obs , is assumed to be the weighted average observed ellipticity over all galaxies in a given 0.1 < z B 0.3 0.3 < z B 0.5 0.5 < z B 0.7 0.7 < z B 0.9 0.9 < z B 1.2 Fig. 1.KiDS-1000 redshift distribution for the five tomographic bins used.Solid lines show the redshift distributions obtained via the SOM technique (Hildebrandt et al. 2021).The shaded bands show the corresponding ranges of photometric redshift point estimates, z B .tomographic bin.We note that all these quantities (shear, ellipticity and additive bias) are complex quantities, where = 1 + i 2 .

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

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(Kitching et al. , 2020)).The average shear in a pixel localised by the unit vector n is given as where the summation is taken over all galaxies in the pixel.For each galaxy, corr i is constructed using Eq. ( 1), with weights, w i , provided by shape measurement techniques such as LensFit (Miller et al. 2007(Miller et al. , 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, ±2 Y m ,is defined as (3) where γ( n) is a complex number and * represents its complex conjugate.The spherical harmonic coefficients for each mode, over a region Ω n of the sky, are 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 W( n); where W( n) = 0 for regions that are not observed, masked out due to bright stars or bad seeing, or contain no galaxies; and 1 otherwise -where there are galaxies in the given pixel.In other words, W( n) is an effective binary mask constructed from the catalogue (which already excludes bad A56, page 4 of 23 regions as detailed in Giblin et al. 2021).Although not the case for the KiDS-1000 dataset, future applications to Euclid data W( 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) where we denote the partial sky quantities with a tilde.Here, the pseudo spherical harmonic coefficients, Ẽ m and B m , are related to the full-sky E-and B-modes as where the mixing kernels W ± mm 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 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: 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 ẼR,i m .Finally, the Ñ term from Eq. ( 10) is simply defined as 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 Ñ .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 σ 2 obs = σ 2 γ + σ 2 , where the shear field variance, σ 2 γ , is mixed with the variance from the intrinsic shape and measurement errors, σ 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.

The mixing matrix
The mixing of modes, introduced in Eqs. ( 8) and ( 9) via the W ± mm terms, is a purely geometrical effect caused by the empty pixels where data are not observed.Differently from galaxy clustering analyses, in cosmic shear the absence of data are taken as the absence of information -if we do not observe lensed galaxies in a region of the sky, we are unable to estimate shear information for that region.Therefore, to deal with the mode mixing caused by the effective geometry of the survey, we construct effective tomographic masks from the KiDS-1000 cosmic shear catalogue.As previously mentioned, for each redshift tomographic bin this mask, W( 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): or, assuming no parity violation modes are present The individual elements of the mixing matrix in the equation above are composed of where W ± mm are the same mixing kernels present in Eqs. ( 8) and ( 9).These can be defined for each tomographic bin as (Hikage et al. 2011) 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 -.
where the 2 × 3 objects above are the Wigner 3 j symbols, calculated using UCLWig3j library 5 , and is the spin-0 spherical harmonic decomposition of the effective binary mask with the corresponding spin-0 spherical harmonic basis, Y m ( n).The mode mixing effect introduced by the mask is then forward-modelled into the theory data-vector, with more details presented in Sect.5.1.An example of the full mixing matrix for the third redshift tomographic bin is presented in Fig. 2.
Here we would like to point out two important details regarding using a binary mask as opposed to a mask containing the shear weights or a variable depth mask.First, a mask containing shear weights, used as W( n), would cause an imprint from the clustering of galaxies into the mixing matrix, causing an artificial excess power for W 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 informa-5 https://github.com/LorneWhiteway/UCLWig3j-this library optimises the calculation of Wigner 3 j symbols using the recurrence relation by Schulten & Gordon (1976).tion into the mask.In this case, it is not clear how one can include this effect into the estimator or the data-vector using a non-binary mask.Ideally, to include variable depth effects, one should take an approach similar to Heydenreich et al. ( 2020) and forward model the effect into the likelihood.Yet, Heydenreich et al. (2020) shows that this effect has very little impact on the inferred cosmological parameters for a KiDS-like survey.Although we do not include this effect in the estimator or in the mixing matrix calculation, variable depth effects are modelled into the covariance via the Egretta mocks (see Sect. 3.4), where their impact is more significant (Joachimi et al. 2021b).Therefore, the additional scatter introduced by variable depth is propagated through the inference by including them directly in the covariance matrix.

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 pix = 12N 2 side ), 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) 6 .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  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.
intuition for the k-scales that enter our analysis as a function of redshift and multipole given the bandwidth -bins used.Following Brown et al. (2005) 7 , the measured pseudo-C s are then binned using eight log-spaced multipole bins between 76 ≤ ≤ 1500, where L < ≤ L+1 are the boundaries of the bandwidth bins centred in L, shown in Fig. 3 as the vertical dotted lines.The mea-sured pseudo-C estimates are shown in Fig. 4 for the E-and Bmodes, 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 χ 2 red = 1.07 for a null-detection, indicating that even with the mode-mixing introduced by the mixing matrix, the contribution of E-modes mixed into B-modes is not enough to bias the cosmological measurements.The auto-power spectrum for E 2 and combinations with it, mostly E 2 − E 3 , have a slightly higher amplitude than the best-fit theory lines.A similar effect was found in previous KiDS-1000 studies but Asgari et al. (2021) have shown that this has almost no impact on the cosmological constraints nor the goodness-of-fit, as shown in Fig. 7  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.

Simulations and covariance matrix estimation
While analytical covariance models are widely used in cosmic shear (Efstathiou 2004;Joachimi et al. 2008Joachimi et al. , 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 Flask8 (Xavier et al. 2016) and Salmo9 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 lognormal matter distribution simulations from Flask with postprocessing using Salmo to populate the matter density field with galaxies based on their local redshift distribution and properties.The mock creation begins by generating the dark matter density field using angular power spectra calculated from 18 redshift tomographic bins of similar intervals, in comoving distances of 150-200 Mpc h −1 , using the fiducial cosmology presented in Table A.1 from Joachimi et al. (2021b).Flask then integrates along the line-of-sight to obtain the shear and convergence fields.Using Salmo, galaxies are positioned in the sky according to the underlying matter density field and following a Poisson distribution, with their shapes given by the weak lensing fields generated by Flask.In addition, Salmo allows for spatially variable redshift distributions to be implemented, reproducing variable depth effects in the final simulated catalogues (more details in 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 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) A56, page 8 of 23  ( 1 , 1 ) ( 2 , 3 ) ( 2 , 4 ) ( 2 , 5 ) where σ i m is defined as the standard deviation of the mcorrections shown in Table A.1 and CE 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 σ i m 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.

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  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.
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 11 + Q 22 , where Q i j is defined in terms of the second-order moments of the two-dimensional angular light distribution, I( n); see Sect.3.1 of Giblin et al. (2021).
-Stellar Density: star sample from Gaia DR210 containing all Gaia DR2 objects in a HEALPix grid of N side = 1024.For each of the quantities mentioned above, we create HEALPix maps using the catalogue object positions with N side = 1024 and the mean value for the observables in each pixel.A spherical harmonic transform is applied to the maps -spin-0 for all scalar quantities and spin-2 transforms for the PSF ellipticities -in order to obtain the systematic's spherical harmonic coefficients, a syst m , and the data's spherical harmonic coefficients, X i m -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 To estimate the errors for the cross-correlations between data and systematics, we obtain covariances in two different ways.For the catalogue quantities, we use a jackknife resampling method.We start by subdividing the survey footprint into 60 distinct regions using a K-Means algorithm on the sphere (Kwan et al. 2017) 11 .While removing the kth region, we cal-culate 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 pvalue 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 8 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-tonoise 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 −2 .This means that this excess correlation is subdominant given the data's estimated covariance in this tomographic bin, leading us to draw the same conclusions as Giblin et al. (2021) in that this correlation has a sufficiently low significance that it has no impact on the cosmological parameter inference.
In summary, our conclusion from the examined potential systematics is that no considerable contamination is biasing the measured E-mode angular power spectra shown in Fig. 4. Therefore, the estimated cosmological parameters shown in Sect.6.1 are purely a result of the cosmological and astrophysical signals captured by our measured pseudo-C s.

Cosmological inference
This section describes the final elements necessary for cosmological inference using tomographic cosmic shear angular power  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.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) 12 modified to directly sample from Gaussian priors 13 for the Multinest sampler (Feroz & Hobson 2008;Feroz et al. 2009Feroz et al. , 2019;;Buchner et al. 2014).It has also been modified to sample directly from the S 8 parameter according to the findings reported in Joachimi et al. (2021b).The Monte Python 3 pipeline was compared with KCAP, the pipeline used in Asgari et al. (2021) and Heymans et al. (2021), and results were found to be consistent.
Following Joachimi et al. (2021b), we report our constraints in two different ways: the traditional marginalised constraints, extracted from the 1D marginalised posteriors for each individual parameter, denoted marginal; and the maximum a posteriori (MAP) with the projected joint highest posterior density (PJ-HPD); this is also referred to as MAP+PJ-HPD.This method ensures, by definition, that the multi-dimensional MAP is always within the PJ-HPD region and it has been demonstrated to accurately report the posterior distributions for each parameter, impacting significantly the reported results from skewed distributions.Monte Python originally does not deal with Gaussian priors and its minimisation algorithm acts directly on the likelihoods, making the MAP estimates from minimising the logposterior 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).

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)  14 .Following the other KiDS-1000 analyses, we fix the sum of neutrino masses to m ν = 0.06eVc −2 , using the one massive species approximation of the normal hierarchy.Not only was this the approach taken by Planck Collaboration VI (2020), but it was shown in Loureiro & Cuceu (2019) that it has no impact on the standard flat ΛCDM cosmological parameters compared to more complex neutrino mass models.
The non-linear matter power spectrum, P nl m (k), is calculated using the halo model from HMCode (Mead et al. 2015) included in Class 15 .This approach incorporates AGN baryonic feedback into the matter power spectrum using a model with two parameters: the halo bloating parameter, η 0 , and the halo massconcentration relation amplitude, A bary .These parameters can be constrained by the following relation from Joudaki et al. (2018): η 0 = 0.98−0.12Abary .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): 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 comoving distance χ can be related as k f K (χ) = + 1/2, we can express the right-hand side C s from Eq. ( 22) as 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 comoving 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: The intrinsic alignments kernel is calculated using the nonlinear alignment (NLA) model from Bridle & King (2007) without redshift evolution, with C 1 ρ 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 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 bandpass filter, showing which scales get mixed for each of the logspaced -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.

External data from SDSS
Since our main objective in this study is to extract cosmological information from late-time probes using galaxy surveys, we have combined our cosmic shear measurements with clustering information.For the current generation of surveys, weak lensing alone has little constraining power beyond the S 8 parameter and it is limited by the degeneracy in the σ 8 − Ω 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 likelihoods16 , which combine pre-reconstruction full-shape17 results and postreconstruction 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 8 , since this parameter is degenerate with the bias.
Using Flask + Salmo mocks, Joachimi et al. (2021b) showed that the BOSS and KiDS-1000 data sets can be considered independent, with their cross-covariance consistent with zero.At the same time, even though the BOSS and KiDS-1000 footprints have a small overlap, Heymans et al. (2021) found that the galaxy galaxy-lensing correlation between the two has very little impact on the cosmological results, mostly due to the necessarily aggressive scale cuts in the data-vector due to bias and IA model limitations.Yet, it does have a small impact on the nuisance parameters, A 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 )σ 8 (z eff ).Here, 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) andTrö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 )σ 8 (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).

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, A56, page 12 of 23 As cosmic shear data are directly sensitive to the amplitude of density fluctuations, σ 8 , or the growth of structure, S 8 , this prior choice is of fundamental importance in our analysis (Hall 2021).In Joachimi et al. (2021b), it was shown that probing this parameter as a derived parameter from the primordial power spectrum amplitude, A s or ln(10 10 A s ), can lead to small biases as the prior becomes informative for S 8 .Therefore, in our analysis, we sample S 8 directly instead of treating it as a derived parameter.This guarantees that our prior on this parameter is flat.
Other parameter priors that are not properly constrained by cosmic shear data are chosen to be ±5σ from independent measurements such as the prior on h, taken from the Riess et al. (2016) measurements, and the prior on ω c = Ω c h 2 , 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 = Ω b h 2 .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.

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+lowE18 (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, σ 8 , Ω m , and S 8 , 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, σ 8 , 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 σ 8 than the band-powers constraints (in orange) but similar to other 2pt estimators presented in Asgari et al. (2021).No degeneracy appears for the SDSS BAO and RSD analysis on this parameter plane, meaning that its combination with cosmic shear considerably breaks the known cosmic shear degeneracy on this plane.This is in line with the 3 × 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.
The clustering posteriors we obtain demonstrate a slightly higher value of σ 8 than the BOSS constraints presented in Heymans et al. (2021), which is more apparent when considering the S 8 = σ 8 √ Ω m /0.3 parameter instead, as shown in  2021) are the addition of highredshift 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 8 measurement robustness from KiDS-1000 cosmic shear analyses.The S 8 constraints with the pseudo-C estimator are slightly larger than the one found with band-powers from Asgari et al. (2021)  increase in comparison.This is in line with the expected increase due to our noisy covariance matrix estimation method, which applies the estimator to simulations instead of using a theoretical covariance.
A summary of S 8 measurements with comparisons to other cosmic shear and CMB probes is displayed in Fig. 10.Overall, weak lensing experiments are consistent between each other.The purple vertical shaded region in Fig. 10 reflects our main constraints from KiDS-1000 pseudo-C 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) 19 are all in agreement with our analysis.Consistency between KiDS-1000 and previous KiDS analyses, such as KV450 (Hildebrandt et al. 2020) and KiDS-450 (Hildebrandt et al. 2017), was explored in detail by Asgari et al. (2021).A lower value for S 8 , when compared to the CMB constraints, is consistently found across a range of different and independent cosmic shear experiments carried out by various research groups.
When comparing the amplitude of intrinsic alignments between the pseudo-C 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: [2PCF] A IA = 0.387 +0.321 −0.374 ; [COSEBIs] A IA = 0.264 +0.424 −0.337 ; 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 χ 2 for the cosmic shear-only analysis is 1.18, reflecting an acceptable p-value of 0.06. 19Here, we quote the mid nominal values as this analysis is closer to the analysis performed for KV-450 and KiDS-1000. .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σ.

Cosmic shear tension with early Universe probes
We now turn our discussion towards the known inconsistencies between cosmic shear surveys and the cosmic microwave background in the amplitude of matter density fluctuations and, consequently, the growth of structure as measured by S 8 .By assuming Gaussianity and adopting the marginal constraints to quantify tension, we initially employ the conventional marginal tension estimate where S A 8 and σ 2 (S A 8 ) are the mean and variance for a given probe A. Under the Gaussianity assumption, τ measures how consistent the difference between S A 8 and S B 8 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 earlytime 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 8 include information from the final WMAP public data release (Bennett et al. 2013).The tension between our constraints from cosmic shear and clustering increases to 2.2σ 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 8 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 With the metric above, we measured d H = 0.95 between our KiDS-1000 pseudo-C and Planck Legacy's S 8 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.

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, √ Ω m /0.3), σ 8 as a derived parameter, the physical baryonic matter density (Ω b h 2 ), the dimensionless Hubble constant (h), and the spectral index (n s ).Included are also two nuisance parameters: the amplitude of intrinsic alignments (A IA ) and the baryonic feedback amplitude (A bary ) as these correlate with some cosmological parameters.When combining the KiDS-1000 pseudo-C 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).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 available20 , 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 2 .
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 A56, page 16 of 23 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 9 galaxies); hence we will implement a noise power spectra subtraction based on the theoretical expectation of the noise.Second, the treatment used for the shear weights in the analysis is sufficient for KiDS, but it may not be the case for Euclid.For higher resolution shear maps such as those required by the Euclid Science Ground Segment (N 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 8 , quite well but it has a known degeneracy in the amplitude of matter density fluctuations and density of matter plane of the parameter space.To break this degeneracy, we combined our pseudo-C 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 exam-ple.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.027 −0.029 , or 9.8 ± 3.3% lower than the S 8 nominal constraints from Planck Collaboration VI (2020).When adding clustering data from SDSS, we increased the precision on the growth of structure parameter, obtaining a measurement of S 8 = 0.771 +0.006 −0.032 , reflecting a 7.8 ± 2.3% lower value when comparing the marginal constraints with Planck's measurement of S 8 .Both constraints are in very good agreement with other KiDS-1000 measurements from Asgari et al. (2021), Heymans et al. (2021) and other cosmic shear surveys such as the Dark Energy Survey (Troxel et al. 2018), the Hyper Supreme Camera survey (Hikage et al. 2019) and CFHTLenS (Joudaki et al. 2017).
Using two different metrics to evaluate tension in a marginalised space, we assessed the tension in the growth of structure parameter between cosmic shear and probes of the cosmic microwave background.Tension with the Planck Legacy TTTEEElowE analysis is around 2.8σ 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.the redshift displacements δ z i to D z i using a Cholesky decomposition of the covariance of the redshift distributions.Therefore, the D z i are uncorrelated.Naturally, this is transformed back into the original correlated space before using the redshift distributions in the forward-modelling described in Sect.5.1.The results reported by Table A.2 are for δ z i , that is the correlated redshifts.Figure A.2 shows good agreement between our analysis and band-powers from Asgari et al. (2021) for all nuisance parameters, except for the amplitude of intrinsic alignments, as discussed in detail in Sect.6.1. Finally,in Fig. A.3 we show a comparison between the SDSS BAO and RSD likelihood and data we present in Sect.5.2 with the clustering likelihood used in Heymans et al. (2021) from BOSS DR12 (Sánchez et al. 2017b).The constraints are in agreement but differ due to containing different data sets and approaches.The SDSS likelihood we use extracts BAO and RSD information from BOSS and eBOSS LRGs, Ly-α autocorrelations 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)    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 E 5 and 2.51 (p-value = 0.010) for B 1 .The S/N is below 5 × 10 −2 for all the bandwidths, meaning that the estimated data covariance is much larger than this correlation.
shows that the pseudo-C B-modes are consistent with zero, indicating that any leakage from any significant correlations between B-modes and systematics are unlikely to interfere with the measured E-modes used for our cosmological analysis.
The highest correlations with B-modes manifest in the crosscorrelations with the PSF E-mode and the stellar density from Gaia DR2, with a reduced χ 2 = 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 where σ X i ,X i is the standard deviation for the auto-power spectra X i , where X = E, B. Figure B.2 shows the cross-correlations between PSF E-mode and B 1 (E 5 ), as well as the S/N for these cross-correlations.From the S/N analysis, we can see that even with a relatively high χ 2 red , the cross-correlations signal is subdominant for both B 1 × PSF-E and E 5 × 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 −2 .
For the case shown in Fig. B.3, we cannot perform a S/N analysis as the units between the CB 2 ,syst and the data covariance do not match.It is possible that we are observing a stellar contamination for the second tomographic bin B-mode.Nonetheless, this possible contamination has no impact in the E-modes used in our cosmological analysis.We know that even by considering E/B-mode leakage, the auto and cross angular power spectra with B 2 shown in Fig. 4 are consistent with zero.Another piece of evidence that this possible contamination does not affect our cosmological signal is the cross-correlation between E 2 and stellar overdensity, with a reduced χ 2 ≈ 1.5 (p-value = 0.13), also shown in Fig. B.3.Had this potential contamination leaked from the B-modes to the E-modes in the second tomographic bin, we would have observed a higher correlation between E 2 and stellar density, which does not seem to be the case.We conclude that any significant correlations between B-modes and systematics do not affect any of our cosmic shear data-vectors in a way that could potentially bias S 8 and other cosmological parameters.

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

Fig. 4 .
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.
of Asgari et al. (2021).A56, page 7 of 23 Notes.Here we present a comparison between characteristics from the sub-sampled and original redshift distributions from Egretta mocks with the Gold Sample SOM N(z)s fromHildebrandt et al. (2021).The original mock values fromJoachimi et al. (2021b) are labelled as J21 and values labelled as sub-sampled are the ones used in this work.

Fig. 5 .
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 i, j E-modes.

Fig. 6 .
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.7to be subdominant to the cosmological signal given the estimated statistical error in this particular tomographic bin.

Fig. 7 .
Fig.7.PSF E-mode and E 5 cross-power spectrum.Top panel: crossangular power spectra between PSF E-mode component and E 5 with the α-correction fromGiblin 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.

8 =
Fig. A.3.Here, we find S SDSS 8 = 0.812 +0.050 −0.049 for the SDSS BAO and RSD data set which is in agreement with both KiDS-1000 pseudo-C , S PCL 0.754 +0.027 −0.029 , and the Planck 2018 results.The two main reasons why our clustering constraints differ from those in Heymans et al. (

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

Fig. 11 .
Fig. 11.Extended marginalised 1D and 2D posteriors for relevant ΛCDM and astrophysical parameters: the matter density parameter (Ω m ), the amplitude of matter density fluctuations (S 8 = σ 8√ Ω m /0.3), σ 8 as a derived parameter, the physical baryonic matter density (Ω b h 2 ), the dimensionless Hubble constant (h), and the spectral index (n s ).Included are also two nuisance parameters: the amplitude of intrinsic alignments (A IA ) and the baryonic feedback amplitude (A bary ) as these correlate with some cosmological parameters.When combining the KiDS-1000 pseudo-C 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) fromAsgari et al. (2021).

Fig
Fig. A.1.Marginalised 1D and 2D posterior distributions for 68% (darker) and 95% (lighter) C.I. in the Ω m , σ 8 , Ω b h 2 , and h plane.The degeneracy broken by the introduction of clustering data to the cosmic shear constraints can be understood as better constraining power from the combination with clustering in the h × Ω b h 2 plane, as well as the BBN-inspired prior.These yield Hubble constant values that are consistent with Planck 2018, in line with what was found by Cuceu et al. (2019) for constraints containing SDSS BAO and RSD information.
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 D z i parameters are shown here as they are sampled, de-correlated displacements calculated from a linear transformation on δ z i using the Cholesky decomposition of the redshift covariance, C z .
used inHeymans et al. (2021) andTröster et al. (2020) fits cosmology directly to measurements of the anisotropic galaxy clustering correlation functions.
Fig. B.1.P-values for a null detection from cross-correlations between the measured B-modes and different systematics.The cross-correlations between B 1 and PSF E-mode are shown in detail in Fig. B.2, while those between B 2 and stellar density are shown in Fig. B.3.

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

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

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

Table 2 .
Hildebrandt et al. (2021)d nuisance parameters.Notes.The Ω b h 2 prior, motivated by BBN constraints fromOlive & Particle Data Group (2014), is a conservative (5σ) equivalent to combining BBN data when external data from SDSS is used, similar to what is shown inCuceu et al. (2019).Here, N(µ, C z ) is a normal distribution with mean µ and covariance C z as estimated byHildebrandt et al. (2021).

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