KiDS-1000 cosmology: Combined second- and third-order shear statistics

This paper performs the first cosmological parameter analysis of the KiDS-1000 data with second- and third-order shear statistics. This work builds on a series of papers that describe the roadmap to third-order shear statistics. We derive and test a combined model of the second-order shear statistic, namely the COSEBIs and the third-order aperture mass statistics $\langle M_\mathrm{ap}^3\rangle$ in a tomographic set-up. We validate our pipeline with $N$-body simulations that mock the fourth Kilo Degree survey data release. To model the second- and third-order statistics, we use the latest version of \textsc{HMcode2020} for the power spectrum and \textsc{BiHalofit} for the bispectrum. Furthermore, we use an analytic description to model intrinsic alignments and hydro-dynamical simulations to model the effect of baryonic feedback processes. Lastly, we decreased the dimension of the data vector significantly by considering for the $\langle M_\mathrm{ap}^3\rangle$ part of the data vector only equal smoothing radii, making a data analysis of the fourth Kilo Degree survey data release using a combined analysis of COSEBIs third-order shear statistic possible. We first validate the accuracy of our modelling by analysing a noise-free mock data vector assuming the KiDS-1000 error budget, finding a shift in the maximum-a-posterior of the matter density parameter $\Delta \Omega_m<0.02\, \sigma_{\Omega_m}$ and of the structure growth parameter $\Delta S_8<0.05\, \sigma_{S_8}$. Lastly, we performed the first KiDS-1000 cosmological analysis using a combined analysis of second- and third-order shear statistics, where we constrained $\Omega_m=0.248^{+0.062}_{-0.055}$ and $S_8=\sigma_8\sqrt{\Omega_m/0.3}=0.772\pm0.022$. The geometric average on the errors of $\Omega_\mathrm{m}$ and $S_8$ of the combined statistics increased compared to the second-order statistic by 2.2.


Introduction
Gravitational lensing describes the deflection of light by massive objects.It is sensitive to baryonic and dark matter and, therefore, ideal for probing the total matter distribution in the Universe.Since the distribution of matter is highly sensitive to cosmological parameters, it is excellent to test and probe the standard model of cosmology, called the Λ Cold Dark Matter model (ΛCDM).Although the ΛCDM model can describe observations of the early Universe, such as the cosmic microwave background ★ E-mail: pburger@astro.uni-bonn.de(CMB;e.g. Planck Collaboration et al. 2020), or the Local Universe, such as the observed large-scale structure (LSS) of matter and galaxies (Sánchez et al. 2017), with remarkable accuracy, it is being put under stress due to tension observed between early and local probes.A ∼ 2  tension is in the structure growth parameter  8 =  8 √︁ Ω m /0.3, where  8 is the normalisation of the power spectrum and Ω m is the total matter density parameter (Hildebrandt et al. 2017;Planck Collaboration et al. 2020; early-time measurements when extrapolated under the ΛCDM model.If these tensions are not due to systematics, extensions to the ΛCDM model are necessary, which will be tested with the next generation of cosmic shear surveys such as Euclid (Laureĳs et al. 2011) or the Vera Rubin Observatory Legacy Survey of Space and Time (LSST, Ivezic et al. 2008).
Commonly, two-point statistics for weak lensing and galaxy positions are used to infer cosmological parameters since they can be modelled accurately and systematic inaccuracies are well understood (Schneider et al. 1998;Troxel et al. 2018;Hildebrandt et al. 2017;Hikage et al. 2019;Asgari et al. 2019Asgari et al. , 2020;;Hildebrandt et al. 2020).Two-point statistics are excellent for capturing the entire information content of a Gaussian random field.However, non-linear gravitational instabilities created a significant amount of non-Gaussian features during the evolution of the Universe, such that the local matter distribution departed strongly from a Gaussian field.Therefore, higher order statistics are needed to extract all the available information in the local LSS of matter and galaxies.Furthermore, such higher order statistics usually depend differently on cosmological parameters and systematic effects such as intrinsic alignment (IA), meaning that a joint investigation of second-and higher order statistics tightens cosmological parameter constraints (see, e.g.Kilbinger & Schneider 2005;Bergé et al. 2010;Pires et al. 2012;Fu et al. 2014;Pyne & Joachimi 2021).Recently used examples of higher order statistics are the peak count statistics (Martinet et al. 2018;Harnois-Déraps et al. 2021), persistent homology (Heydenreich et al. 2021), density split statistics (Gruen et al. 2018;Burger et al. 2023), and the integrated three-point correlation function used in (Halder et al. 2021;Halder & Barreira 2022), along with a second-and third-order convergence moment analysis (Gatti et al. 2022).
This work considers second-and third-order shear statistics, where the former probes the variance and the latter the skewness of the LSS at various scales.Our chosen second-order statistic is the   -modes of the COSEBIs (Schneider et al. 2010), and the third-order statistic M 3 ap is described in Schneider et al. (2005) and recently measured in the Dark Energy Survey Year (DES) 3 Results (Secco et al. 2022).Furthermore, this work belongs to a series of papers that aim for cosmological parameter analyses using third-order shear statistics.In the first paper, Heydenreich et al. (2023, hereafter H23), we validated the analytical fitting formulae for a non-tomographic analysis and the conversion from three-point correlation function (3PCF) of cosmic shear to M 3 ap .We found that M 3 ap , even though they are combined from the shear 3PCF at different scales, contain a similar amount of information on Ω m and  8 as the 3PCF itself.The fact that M 3 ap can be measured from 3PCF is very convenient since this allows unbiased estimates for any survey geometry.A fast computation method of the aperture mass statistics by measuring the shear 3PCF is tackled in Porth et al. (2023).Lastly, adding third-order statistics increases the dimension of the data vector significantly, and an analytical expression for the covariance is preferred, which is derived and validated for M 3 ap in a non-tomographic setup in Linke et al. (2023).However, as this analysis considers combining second-order statistics with M 3 ap in a tomographic setup, and a joint covariance matrix has not been derived yet (Wielders et al. in prep), we must still rely on numerical simulations to determine the covariance matrix.
This article presents the first cosmological parameter analysis using the fourth data release of KiDS (KiDS-1000), combining second-and third-order shear statistics.We show the cosmological results, preceded by validation of several extensions of the analytical model to allow for a tomographic analysis and to in-clude astrophysical effects such as IA (e.g.Joachimi et al. 2015) or baryonic feedback processes (Chisari et al. 2015).We aim to find the smallest set of scales that retain most of the cosmological information.
The paper is structured as follows: In Sect.2, we review the basics of the second-and third-order shear statistics, extend the modelling to a tomographic analysis and describe our method to model the IA analytically.In Sect.4, we describe the KiDS-1000 data and in Sect.5, we introduce the simulation data used to validate our model and to correct it for baryonic feedback processes.In Sect.6, we briefly review our method to perform cosmological parameter interference.In Sect.7, and 8, we validate our analysis pipeline against several systematics.The final cosmological results are presented in Sect.9, and we present our conclusions in Sect.10.

Theoretical background
This section offers a review of the basics of weak gravitational lensing formalism and aperture statistics.For more detailed reviews, we refer to Bartelmann & Schneider (2001), Hoekstra & Jain (2008), Munshi et al. (2008), Bartelmann (2010), andKilbinger (2015).In this work, we assume a spatially flat universe, such that the comoving angular-diameter distance is expressed as   [ ()] = (), where () is the comoving distance at redshift .Given the matter density ,(  , ), at comoving position,   , and redshift, , the density contrast is (  , ) = (  ,) ρ() − 1, where ρ() is the average matter density at redshift, .The dimensionless surface mass density or convergence, , for sources at redshift, , is determined by the line-of-sight integration as where    is the angular position on the sky,  0 is the Hubble constant, and  is the scale factor.The second argument of  simultaneously describes the radial direction and the cosmological epoch, related through the light-cone condition | d| = () d.

Limber projections of power-and bispectrum
Given the Fourier transform δ of the matter density contrast, the matter power spectrum,    (, ), and bispectrum,     ( 1 ,  2 ,  3 , ), are: where  D is the Dirac-delta distribution.The statistical isotropy of the Universe implies that the power-and bispectrum only depend on the moduli of the -vectors.
The projected power-and bispectrum can then be computed using the Limber approximation (Limber 1954;Kaiser & Jaffe 1997;Bernardeau et al. 1997;Schneider et al. 1998;LoVerde & Afshordi 2008), where ℓ 3 = |ℓ ℓ ℓ 1 + ℓ ℓ ℓ 2 |, and  () ( ) denotes the lensing efficiency and is defined as with  () ( ) being the redshift probability distribution of the th tomographic  ph -bin.Since the Limber approximation breaks down for small values of ℓ (Kilbinger et al. 2017), we consider only scales below 4 • to model M 3 ap .In principle, we could have also used the modified Limber approximation and shift ℓ 1,2,3 → ℓ 1,2,3 +1/2, but since we found a difference at maximum for the largest filter radii of 1.5%, we neglected it here.To model the non-linear matter power spectrum    , we use the revised HMcode2020 model of Mead et al. (2021) and for the matter bispectrum     the BiHalofit of Takahashi et al. (2020).

Non-linear alignment model
The impact of galaxy IA is a known contaminating signal to the cosmic shear measurements that must be accounted for in all weak lensing studies.To model the effects of intrinsic alignment, we use the non-linear alignment (NLA) model (Bridle & King 2007), which is a one-parameter model described as where  + is the linear growth factor at redshift , C1 = 5 × 10 −14  −1 ⊙ ℎ −2 Mpc 3 , as calibrated in Brown et al. (2002), and  IA captures the coupling strength between the matter density and the tidal field.
By considering the tidal alignment field,  I , as a biased tracer of the matter density contrast field , neglecting all higher order bias terms, we get: With this, we find that: where    (, ) is the non-linear matter power spectrum.The projected power spectra then follow to and the total projected power spectrum becomes The term (ℓ) describes the actual lensing signal, which is given in Eq. ( 4) with the weighting kernel in Eq. ( 6).The II contribution describes how two galaxies spatially close together tend to be aligned.The term GI(  >   ) describes the fact that high matter density regions align the lower redshift galaxies but also affect the shear of the background galaxies.While the II term is dominant if galaxies of the same tomographic bin are considered, GI and IG start to dominate if galaxies of separated tomographic bins are considered.With   <   , the term GI is expected to vanish for two tomographic bins with no significant redshift overlap.
Following the ansatz for modelling IA in the power spectrum, we get where     is the non-linear matter bispectrum, which we calculated with BiHalofit.The projected bispectra are The total projected bispectrum is: where the tuple is l = (ℓ 1 , ℓ 2 , ℓ 3 ).The interpretation of all these terms is analogous to the ones from the power spectrum, where ( l) is the actual pure lensing signal given by Eq. ( 5) with the weighting kernel in Eq. (6).

Aperture mass statistics
One of the major problems of weak lensing mass reconstruction techniques is the mass-sheet degeneracy (Falco et al. 1985;Schneider & Seitz 1995, hereafter MSD), which corresponds to adding a uniform surface mass density without affecting lensing observables such as the shear.However, it is possible to define quantities invariant under the MSD, one example being the aperture mass statistics (Schneider 1996;Bartelmann & Schneider 2001).Another advantage of aperture mass statistics is that they separate the signal into so-called E-and B-modes (Schneider et al. 2002), where, to leading order, the weak gravitational lensing effect cannot create B-modes.Lastly, as shown in H23, aperture statistics are an excellent strategy to compress thousands of bins of 3PCF into a few hundred M 3 ap bins.
Article number, page 3 of 19 A&A proofs: manuscript no.main The aperture mass M ap at position    with filter radius  ap is defined through the convergence , as follows: where   ap ( ′ ) is a compensated filter such that ∫ d ′  ′   ap ( ′ ) = 0.The tangential shear component  t of the complex shear in Cartesian coordinates  =  1 + i 2 , is defined as  t = −ℜ(e −2i ), where  is the polar angle of  ′  ′  ′ .Given the tangential shear  t , the aperture mass M ap can also be calculated as where   ap is related to   ap via We define   ap () =  −2 ap (/ ap ), denote by û() the Fourier transform of  and use the filter function introduced in Crittenden et al. (2002), (24)

Modelling aperture mass moments
The expectation value of the aperture mass M ap ( ap ), which approximates the ensemble average over all positions   , vanishes by construction.However, the second-order (variance) of the aperture mass is nonzero and can be calculated as Equivalently, the third-order moment of the aperture statistics M 3 ap , can be computed from the convergence bispectrum via (Jarvis et al. 2004;Schneider et al. 2005) where Later, we differentiate between equal filter radii, for which  ap,1 =  ap,2 =  ap,3 and non-equal filter radii, where the values of  ap, are all allowed to vary.

Modelling COSEBIs
where  ± () are filter functions with support in [ min ,  max ] (Schneider et al. 2010).If where  0,4 are the zeroth and fourth order Bessel function, then the COSEBIs offer the advantage of cleanly separating all welldefined E-and B-modes within the range [ min ,  max ].This is not given, for instance, for second-order aperture mass statistics, as this would require information of  ± over the full space.Analytically, the COSEBIs can be calculated from the E-mode power spectrum (ℓ) defined in Eq. ( 4) and a B-mode angular power spectra Since B-modes cannot be created by gravitational lensing directly, we neglected the modelling of   for this work.Therefore, we focus, henceforth, only on the   modes from the COSEBIs.

Measuring shear statistics
As discussed in H23 M 3 ap can be estimated from data in three ways, which we review below.The first uses the convergence field, , the second uses the shear field, and the third uses correlation functions.The latter is used to measure the real data and compute the covariance matrix used for the real data analysis.

Measuring shear statistics from convergence maps
If the convergence field, , is available (e.g. from simulations), the easiest way to measure M ap is to use Eq. ( 21).If the convergence field does not contain masks, this estimator is also unbiased as long as a border of the size of the filter function is removed from the aperture mass field before calculating the spatial average.However, no boundaries need to be removed for unmasked full-sky convergence fields to get an unbiased estimator.To estimate the aperture mass from (full-sky) convergence maps, the maps are smoothed with the healpy1 (Zonca et al. 2019) function smoothing.This function needs a beam window function created by the function beam2bl, which is determined from the corresponding   ap filter.
In the absence of B-modes the   modes can also be calculated by using convolutions as (Schneider et al. 2010): where with  + introduced in Eq. ( 29).

Measuring shear statistics from galaxy catalogues
In a real survey, the convergence is not observable and can be inferred only from the measured shear.However, estimating M ap from these reconstructed convergence maps is not accurate as the reconstructed convergence is necessarily smoothed and potentially also contains other systematic effects caused by masks or boundaries of the survey (Seitz & Schneider 1997).However, as motivated by Eq. ( 22), the aperture mass can also be estimated from an unmasked shear field (e.g. for simulations).Given a galaxy catalogue, where galaxies are only found at specific positions such that the number of galaxies within an aperture varies on the sky, Eq. ( 22) needs to get modified to where the sum over the filter functions serves as a normalisation,  t are the observed galaxy ellipticities converted into their tangential component, and   are their respective positions.We note that we sampled the galaxies on a grid using a cloud-in-cell method2, so that convolutions can determine the aperture masses.Since we used this approach solely for finite fields, we applied the same cut-off of 4 max ap to all aperture mass maps, where  max ap is the largest filter radius.For both approaches, where the aperture mass is determined, the second-and third-order aperture statistics follow by multiplying the respective aperture mass maps pixel-wise and then taking the average of all pixel values.

Measuring shear statistics from correlation functions
The first two methods to estimate unbiased aperture statistics are not applicable to observed data with masks.Another disadvantage is the removal of the edges, which can be significant for a complex survey footprint, leading to decreased statistical certainty.The best method to estimate aperture shear statistics or the   is by measuring them from two-and three-point shear correlation functions (Schneider et al. 2002;Jarvis et al. 2004;Schneider et al. 2005).The advantage of correlation functions is that they can be estimated for any survey geometry.The measurement of the two-point shear correlation functions  ± ( ap ) is unbiased and easily performed by treecorr (Jarvis et al. 2004) and the conversion to COSEBIs are given in Eq. 27 and Eq.28.
For the measurement of the aperture statistics, M3 ap , we refer to our companion paper, namely, Porth et al. (2023).It describes an efficient estimation procedure of the natural components of the shear 3PCF (Schneider & Lombardi 2003), which is then transformed into aperture statistics.The corresponding equations for a tomographic set-up can be found in their Section 5.2.

Observed data
This analysis explores the fourth data release of KiDS (Kuĳken et al. 2015(Kuĳken et al. , 2019;;de Jong et al. 2015de Jong et al. , 2017)).The weak lensing data observed with the high-quality VST-OmegaCAM camera is public3.It is collectively called 'KiDS-1000' as it covers ∼ 1000 deg2 of images, which is then reduced to an effective area of 777.4 deg 2 after masking.The significant advantage, compared to previous weak lensing surveys and data releases, is its overlap with its partner survey VIKING (VISTA Kilo-degree Infrared Galaxy survey, Edge et al. 2013), which observes galaxy images at infrared wavelength.Therefore, galaxies were observed in nine optical and near-infrared bands, , , , , ,  , , ,  s , allowing for better control over redshift uncertainties (Hildebrandt et al. 2021, hereafter H21).The KiDS-1000 cosmic shear catalogue is divided into five tomographic  ph -bins, whose redshifts are calibrated using the self-organising map (SOM) method4 described in Wright et al. (2020).The redshift distributions of all five tomographic bins are shown in Fig. 1, and were initially presented in H21.The residual systematic uncertainties on the redshift distributions are listed in Table 1 and are included in this work as nuisance parameters which we marginalise over.We note that they are correlated, which we account for by using their correlation matrix for the marginalisation.The galaxy shear ellipticities and their corresponding weights, , are estimated by the lensfit tool (Miller et al. 2013;Fenech Conti et al. 2017;Kannawadi et al. 2019) and are described in more detail in Giblin et al. (2021).These shear-related systematic effects shift the  8 parameter by (at most) 0.1 when measured by cosmic shear two-point functions.The resulting systematics are stated in Table 1, where we marginalise over the shear multiplicative -bias correction in the resulting posteriors.

Simulated data
We use several simulated data sets created to resemble the observed KiDS-1000 data, to validate our inference pipeline, study the impact of key systematic uncertainties, and forecast the expected KiDS-1000 analysis.
In particular, we used the full-sky gravitational lensing simulations described in Takahashi et al. (2017, hereafter T17) to generate data vectors and numerical covariance matrices.Notes.The second column shows the effective number density, which also accounts for the lensfit weights.The third column shows the mean and uncertainty on the redshift bias taken from H21.The redshift bias uncertainties are correlated, which we accounted for by the correlation matrix given in H21.The fourth column displays the measured ellipticity dispersion per component,   , as given in Giblin et al. (2021).The fifth column shows the shear multiplicative -bias correction updated in van den Busch et al. (2022).
whose strength we regulate with a free parameter in the posterior estimation.

Takahashi simulations
Since the T17 simulations are used in this series of previous works, we only reiterate the essential details here.The T17 simulations follow the non-linear evolution of 2048 3 particles evolved in a large series of nested cosmological volumes with side length starting at  = 450 Mpc/ℎ at low redshift, and increasing at higher redshift, resulting in 108 different full-sky realisations.These were produced by the Gadget-3 -body code (Springel 2005) and are publicly available5.The cosmological parameters of the matter and vacuum energy density are fixed to Ω m = 1 − Ω Λ = 0.279, the baryon density parameter to Ω b = 0.046, the dimensionless Hubble constant to ℎ = 0.7, the normalisation of the power spectrum to  8 = 0.82, and the spectral index to  s = 0.97.The shear information of the T17 simulations is given in terms of 108  and  full-sky realisations, where each realisation is divided into 38 ascending redshift slices.
To reproduce the KiDS-1000 data for a given tomographic  phbin shown in Fig. 1, we built the weighted average of the first 30  and  redshift slices.The weight for each redshift slice is measured by integrating the () over the corresponding width of the redshift slice.We decided on two approaches (convergence maps vs. galaxy ellipticity) to measure the data vectors and covariance matrices from the T17 simulations.While the first approach (Sect.5.1.1)has the advantage that a large number of mock data is available to measure a reliable covariance matrix even for large data vectors, the second approach (Sect.5.1.2) has the advantage that the exact galaxy positions are used, which implies that the holes (masks) match the data.

Convergence mock data
The first approach is to convolve the full sky convergence maps with the  ±,n ( ap ) or  ( ap ) filters, then multiply with corresponding other smoothed maps and then take the spatial average.For the covariance matrix, we divided the smoothed convergence maps into 48 sub-patches, where each patch has an area of 860 deg 2 .Since the  maps are first convolved with the filter functions and the sub-patches have common edges, the individ-ual patches are not fully independent from each other, which decreases the covariance matrix.However, we found in our testing that selecting only 18 sub-patches, that have no common edges, gives almost identical results.Shape noise is added to the convergence maps by drawing random numbers from a Gaussian distribution with a vanishing mean and a standard deviation, as follows: with the pixel area as  pix = 0.74 arcmin 2 (nside = 4096),  eff as the effective galaxy number density that included the lensfit weights, and   as the shape noise contribution given in Table 1.With 48 sub-patches and 108 realisations, the covariance matrix is measured from 5184 mock data.As the actual KiDS-1000 area is roughly 777.4 deg 2 , we rescaled the covariance by 860 deg 2 /777.4 deg 2 ≈ 1.1.The reference data vector is measured from one T17 realisation with a resolution of  pix = 0.18 arcmin 2 (nside = 8192) without shape noise.

Galaxy shear catalogue mock data
For real surveys with a complex topology, the convolved maps give biased results (Seitz & Schneider 1997).Therefore, we decided on our second approach to measure our statistics from second-and third-order shear correlation functions.For this, we created galaxy shear catalogues from projected  and  fields by extracting the shear information only at the true positions of the observed galaxies in KiDS-1000.Since the correlation functions need to be measured to very small scales, we used realisations with a pixel resolution of  pix = 0.18 arcmin 2 (nside = 8192).We checked that this resolution is sufficient by comparing the data vectors obtained from catalogues constructed from the same initial conditions on the reference resolution (nside = 8192) and a higher resolution (nside = 16384), finding a maximum deviation of 1% for an aperture filter scale of 4 ′ .To increase the number of mock data, we shifted the galaxy positions 18 times by 20 • along the lines of constant declination and extracted the shear information at the new positions.Afterwards, the shifted galaxy positions and the corresponding lensing information are shifted back to the original footprint.The back shifting is done only for simplicity and has no physical reason.With this procedure, 1944 almost independent mock data were created, from which the covariance and the reference data vector were measured.To add shape noise, we combined the two-component reduced shear  = /(1 − ) of each object with a shape noise contribution,  s , to create observed ellipticities  obs (Seitz & Schneider 1997), as follows: The quantities here are all complex numbers and the asterisk ' * ' indicates complex conjugation.To mock the  s , we use the observed ellipticities  obs KiDS of the KiDS-1000 data and randomly rotate them to erase the underlying correlated shear signal.This procedure offers the advantage that the resulting distribution of  s matches the distribution of  obs KiDS .The estimated mean and dispersion   are given in Table 1.Furthermore, when computing our statistics via correlation functions, we need to consider the corresponding weight  for each shape measurement, which ensures that we use the correct effective number density.Since the lensing weights and the intrinsic ellipticities of source galaxies are correlated, adding the rotated ellipticities to the shear signal preserves this correlation.

Data vector measurement and modelling
We follow the analysis of Asgari et al. (2021, hereafter A21) as it is the fiducial cosmic shear analysis of the KiDS-10006.We only use the first five   -moments determined from two-point correlation functions, measured from 0. ′ 5 to 300 ′ in 400 radial bins.The resulting   are shown in Fig. 2, where we checked that the   are consistent with zero.The   data vector for 7 and  is the mean of all 1944 T17 mock data; and for the  data vector, we used one realisation with resolution  pix = 0.18 arcmin 2 (nside = 8192).The analytical model was computed using Cos-moSIS (Zuntz et al. 2015a).

cosmo-SLICS+IA
The cosmo-SLICS+IA simulations are cosmic shear mock galaxy catalogues infused with the non-linear alignment model of Bridle & King (2007), which is ideally suited for testing and validating our analytical IA model.They are based on -body simulations of identical box size and particle density as the SLICS (Harnois-Déraps et al. 2018), which were already used in previous works of this series; we therefore only summarise the essential details.
The cosmology corresponds to the fiducial cosmo-SLICS model presented in Harnois-Déraps et al. ( 2019), with Ω m = 0.2905, Ω Λ = 0.7095, Ω b = 0.0473, ℎ = 0.6898,  8 = 0.836,  0 = −1.0 and  s = 0.969.We use the full set of 50 simulated galaxy catalogues, each covering 100 deg 2 and reproducing the KiDS-1000 () (see Fig. 1) and galaxy number densities  eff specified in H21.As we use these mock data only to infuse IA effects into the T17 data vector, shape noise is not included.Following the methods described in Harnois-Déraps et al. ( 2022), the intrinsic ellipticity components  IA 1/2 of these galaxies are computed as where  IA is defined in Eq. ( 7), and  is the gravitational potential.
The partial derivatives of the gravitational potential describe the Cartesian components of the projected tidal field tensors.The  IA 1,2 terms were then combined with the noise-free cosmic shear signal using Eq. ( 36), resulting in an IA-contaminated weak lensing sample that is consistent with the NLA model.As these mock data cover a square patch without any masking, we used the methods described in Sect.3.2 to estimate the aperture statistics.Furthermore, they were used only to validate the IA modelling, which we discuss in Appendix A.

Magneticum simulations
The feedback processes due to baryonic matter significantly affect the distribution of the LSS, such that the clustering of the matter is reduced on intra-cluster scales by up to tens of per cent (van Daalen et al. 2011).However, in a quantitative sense, this suppression is not well understood (Chisari et al. 2015).We use the Magneticum simulations (Hirschmann et al. 2014) to investigate the impact of baryonic feedback processes.Magneticum was run using Gadget3 code which is a more efficient version of Gadget 2 (Springel 2005) that includes modern smoothed particle hydrodynamics (Beck et al. 2016).The dark matter particle mass is 6.9 × 108 ℎ −1  ⊙ and gas particle mass 1.4 × 10 8 ℎ −1  ⊙ .The underlying matter fields were constructed from the Magneticum Pathfinder simulations8 in the Run-2 with a comoving volume of side 352 ℎ −1 Mpc and Run-2b with a comoving volume of side 640 ℎ −1 Mpc, and are described in Hirschmann et al. ( 2014) and Ragagnin et al. (2017), respectively.The Magneticum simulations account for radiative cooling, star formation, supernovae and active galactic nuclei.
The Magneticum shear maps used in this work were first presented in Castro et al. (2018).Although the underlying creation of these shear maps is similar to the production of the cosmo-SLICS, meaning that the mock data follow the same redshift distribution and number density as the KiDS-1000 data, there are three main differences.First, only ten instead of 50 pseudoindependent light cones are available.Second, the galaxies are placed at the exact galaxy positions as in the observed data rather than randomly.Third, as these mock data are flat sky simulations and cover only an area of 100 deg 2 , the KiDS-1000 footprint is divided into 18 regions.This results in 18 catalogues for each of the ten pseudo-independent light cones.Due to the non-trivial geometry of the masks, we used correlation functions to estimate the shear statistics.
To incorporate baryonic feedback processes into our modelling pipeline using the Magneticum simulations, we calculated the data vector, d DM , using dark matter only and the data vector, d DM+BA , where dark matter and baryons are included.With them we modified the model vector, m, to where we introduce the parameter  ba .A value of zero for  ba means no baryonic feedback, and a value of one is exactly the strength of baryonic feedback as in the Magneticum.Furthermore, we interpolate between zero and one and extrapolate to KiDS-1000 mock data that are also used to compute the covariance matrix with a resolution  pix = 0.18 arcmin 2 and shape noise.The red dots represent the data vector measured from one full-sky T17 realisation with a resolution  pix = 0.18 arcmin 2 and no shape noise.The grey band indicates the expected uncertainty from the KiDS-1000 survey.
Fig. 3: Same as Fig. 2, but here the measured data and model are the M 3 ap vector for equal-scale aperture filter radii  ap ∈ {4 ′ , 6 ′ , 8 ′ , 10 ′ , 14 ′ , 18 ′ , 22 ′ , 26 ′ , 32 ′ , 36 ′ } in the T17 mock data for some selected  ph -bin combinations.The full data set is shown in five, which we chose arbitrarily.We note that for the joint analysis of second-and third-order statistics, we need at least extrapolation to  ba = 2 (at 95% confidence), which approximately agrees with the baryonic strength of other hydrodynamical simulations (Martinet et al. 2021).The third-order analysis alone is bound by the prior  ba = 5.Although extrapolating to these large values is probably incorrect, we introduced  ba only to give the model some flexibility to account for baryonic feedback effects without having a direct physical meaning such as the ejected gas.

Cosmological parameter inference methodology
In the following sections, we determine multiple posterior distributions using Markov chain Monte Carlo (MCMC) samplings for different model ingredients.
We varied the two cosmological parameters Ω m and  8 while fixing the remaining parameters to the T17 cosmology.Additionally, we varied the intrinsic alignment amplitude,  IA , and the parameter  ba that accounts for the strengths of the response to baryonic feedback9 and was described in Sect.5.3.Lastly, we always marginalise over the shifts of the redshift distributions ⟨⟩ given in H21 and the -bias (Giblin et al. 2021) of all five tomographic bins shown in Table 1.For these nuisance parameters, we assume Gaussian priors and account for the fact that the uncertainties of shifts of the redshift distributions are correlated using its correlation matrix.We give an overview of the priors in Table 2.If the covariance matrix, C, is measured from simulations, it is a random variable.We followed the method from Percival et al. (2022), which leads to credible intervals that can also be interpreted as confidence intervals with approximately the same coverage probability.The posterior distribution of a model vector, m, that depends on   parameters, , if the covariance matrix C is measured from  r mock data, is given by: where d is the measured data vector and The power-law index  is with  d being the number of data points and To give a quantified value for the comparison of different modelling choices, we use the Figure of Merit (FoM), which we calculate as where  is the parameter covariance matrix of the  8 -Ω m plane resulting from the MCMC process.
Depending on the used filter combinations, the modelling of M 3 ap takes around 10 to 20 min, which stands as an obstacle to directly running an MCMC using the model.To circumvent this issue, we use the emulation tool contained in CosmoPower (Spurio Mancini et al. 2022).We trained the emulator on 3000 model points for the analysis based on mock data in the parameter space {Ω m ,  8 ,  IA } and 5000 model points for the real data analysis in the parameter space {Ω m ,  8 ,  IA , ⟨⟩}, which we distributed in a Latin hypercube for the given prior range, and fixed all other parameters to the ones used in the T17 simulations.The -bias and the  ba do not need to be emulated as they follow a simple correction formalism that is applied during the MCMC sampling.The ⟨⟩ of each  ph -bin for the training and testing are distributed uncorrelated between [−0.06, 0.06].The accuracy of the emulator is tested by comparing the emulator prediction with the model at 500 independent points in the same parameter space.The fractional error of each vector element is smaller than 2% for   and smaller than 5% for M 3 ap , which is well within the accuracy of the HMcode2020 or the BiHalofit model itself.We used a Metropolis-Hastings sampler for MCMC10, where we used 1000 walkers running each 20 000 steps and cut the first 2000 steps away to ensure that the posteriors are not biased by the burn-in phase.

Reduction of data vector combinations and filter radii
Generally, as long as the covariance is converged, the more information is used, the higher the constraining power.However, for third-order statistics with four different filter radii and five tomographic bins, the data vector contains 700 elements.To obtain a reliable covariance matrix, roughly ∼ 10 times the dimension of the data vector is needed.This implies that modelling and measuring the model/data vector and covariance matrix are time-consuming and require thousands of simulations.Therefore, it is inevitable for this and for future analyses to compress the data vector without substantial information loss.We note that this section is only concerned with reducing the dimension of the data vector and is not meant as a proper forecast, which we do in the following sections.The spatial resolution of the mock data used to compute the covariance induces an error smaller than < 1% in the   covariance matrix.However, as we only reduce the elements of M 3 ap , while using all combinations of the   , this is unproblematic for the purpose of this section.We note again that the data vector was measured from one realisation with a resolution of  pix = 0.18 arcmin 2 , and that we use the  based covariance for this section to ensure the covariance matrix is converged even for the data vector with the largest dimension.As a first check, we show in Fig. 4 multiple element choices that discard a significant part of the M 3 ap data vector.Using only the five auto-tomographic bins is not a good choice for reduction, as seen in the fact that the FoM is reduced by 32%.However, using only equal-scale aperture radii is a better reduction, reducing the data vector to 28% while reducing the FoM by only ≈ 8%.Using only equal-scale aperture radii also has the advantage of being modelled faster.Therefore, we continue using only the equalscale aperture radii for the rest of this work, which reduces the dimension of the data vector from 775 to 215.
Next, we investigated the choice of aperture filter radii.The aperture filter radii under consideration are  ap ∈ {4 ′ , 6 ′ , 8 ′ , 10 ′ , 14 ′ , 18 ′ , 22 ′ , 26 ′ , 32 ′ , 36 ′ }.The resulting posteriors are displayed in Fig. 5, which reveals that using only the large filter radii above > 22 ′ have the worst constraining power.However, using only the four smallest filter radii is also not the best choice, as they are largely correlated and miss the large-scale information.The best choice of filter function is to use a filter from each range such as  ap ∈ {4 ′ , 8 ′ , 14 ′ , 32 ′ }.
Next, we considered the full data vector, meaning the M 3 ap , and the   for the next compression strategy.The idea is to decrease the number of elements by considering only those with the highest constraining power on  8 .We start with the element with the highest /, then consecutively add those vector elements that maximise the  8 Fisher information content.This Fisher information is calculated as (Tegmark et al. 1997) where the partial derivatives are computed with a five-point stencil beam (Fornberg 1988), where  ± 8 =  T17 8 ± 0.02  T17 8 and  ±± 8 =  T17 8 ± 0.04  T17 8 .For this analysis, we used only the equal-scale filter radii  ap ∈ {4 ′ , 6 ′ , 8 ′ , 10 ′ , 14 ′ , 18 ′ , 22 ′ , 26 ′ , 32 ′ , 36 ′ }.The first ∼ 200 elements are sufficient to get converged posteriors.It is also interesting to see which elements help increase the constraining power.Unsurprisingly, the first elements are all COSEBIs, but among the first 100, approximately half are M 3 ap elements.Furthermore, it is interesting to see that cross-tomographic bin elements  are more likely to be selected by our method, which is expected because they have a higher / than auto-tomographic bin elements.Nevertheless, we also observe that using equal-scale filter radii  ap ∈ {4 ′ , 8 ′ , 14 ′ , 32 ′ } results in a better FoM then using the best 215 elements.This is likely because we optimised only the  8 parameter here while fixing all others.Therefore, we also checked for the equal scale filter case if a principal component analysis (PCA) applied to the covariance matrix performs better in compressing the data.However, it needs more elements to get the same constraining power as our FoM maximiser.This is probably because a PCA considers only the covariance matrix and ignores the derivatives, meaning that a PCA is not necessarily sensitive to cosmology.
Finally, we give in Table 3 an overview of the FoM for the Ω m - 8 plane and the size of the data vector for some more element choices.As the covariance matrix is measured from 5184 mock data, we can assume that for all data vector dimensions under consideration, the covariance matrix is converged.Nevertheless, for the next sections, we use the covariance matrix measured from 1944 galaxy shear catalogues, limiting us to a maximum dimension of our data vector ∼ 200 elements.Given the investigations in this section, we restrict the further sections to the case where we use all   elements and all M 3 ap elements with the equal-scale filter radii  ap ∈ {4 ′ , 8 ′ , 14 ′ , 32 ′ }, as this resulted in the best FoM for the restricted dimension of the data vector.

Validation of data vector estimator
In a real survey analysis, two further difficulties arise.The first issue is that the lensing information is not given in terms of convergence maps but by point estimates (galaxies).The finite area where galaxies are measured implies that no lensing information is available outside that area.It is, therefore, necessary to measure the statistic from correlation functions described in Sect.3.3.Second, these point estimates are the reduced shear  = /(1 − ), which increases the signal and therefore needs to be accounted for.To correct these effects, we measured data vectors without shape noise and with the largest available resolution  pix = 0.05 arcmin 2 .Since these data vectors result only from one realisation without shape noise, we can expect that the ratio is a good approximation for the reduced shear effect.Although the deviations are small (as seen in Figs.A.3 and A.4), we decided to scale the model vectors by the ratio of  and  data vector.
The resulting posteriors are shown in Fig. 6.Our first observation is that similar to the finding of H23, combining second-order with third-order shear statistics significantly improves the constraining power on  8 and Ω m .Compared to the   -only, the  8 -Ω m FoM increases by 93% and for the M 3 ap -only by 233%.For these improvements, we have not considered that the posteriors of the individual statistics are bound by the priors on Ω m , meaning that the improvements are lower bounds.Compared to the analysis based on  mock data (see Fig. 5), the posteriors are broader because the covariance based on  maps also uses information outside the patches since the boundaries were not removed.This decreases the variance between the patches.A further difference is that the  analysis is not subject to masks and uses a slightly larger effective area.Although we rescaled the covariance to correct for the different effective areas, we found in Linke et al. (2023) that this rescaling is not necessarily accurate.Lastly, we notice that our modelling within the KiDS-1000 uncertainty is accurate, which we quantified by measuring the shift of the maximum of the posterior (MP) from the true values in the matter density parameter ΔΩ m < 0.02  Ω m and in the clustering amplitude Δ 8 < 0.05   8 with respect to the averaged noisy mock data vector and the KiDS-1000 uncertainty.We define the MP as the maximum of the one-dimensional marginal distributions.

Cosmological results
Finally, we are ready to present the first cosmological constraints from the KiDS-1000 data using second-and third-order shear statistics, displayed in Figs.7 and 8.The data vectors are described in more detail in A21 for the   , and the M 3 ap in Porth et al. (2023).Given the fact that the covariance matrix is measured only from 1944 realisations using the reduced shear , we decide to build our data and model vector from all   modes and M 3 ap with all equal filter radii  ap ∈ {4 ′ , 8 ′ , 14 ′ , 32 ′ }.To control whether the model accurately fits the data, we minimised the  2 real from the real data and the  2 T17 from each of the 1944 T17 mock data used to compute the covariance matrix.To estimate the probability of measuring a  2 >  2 real (-value), we counted the number of  2 T17 that are greater than  2 real divided by 1944.The resulting -value for the   -only, the M 3 ap -only, and the combination are given in Table 4.The resulting -values for all combinations are better than 0.1, indicating that our covariance is well matched to the observed data and our model is accurate enough to describe the data.Our maximum posterior  2 value increases to 81 if we swap our   covariance matrix to the analytical expression in Joachimi et al. (2021) calculated at the MP parameter values in A21.This is unsurprising as the numerical covariance is computed at the T17 cosmology, giving a signal larger than that computed at the MP of A21.Furthermore, our covariance matrix is measured from reduced shear mock data that slightly increases the covariance matrix.We also notice that the model and the data seem inconsistent for some  ph -bin combinations.However, adjacent COSEBI modes are highly correlated, so visually inspecting the model's goodness of fit to the data is misleading.We refer to A21 for further details on this discrepancy.
We show the resulting posteriors in Fig. 9, where we marginalised over the shift in the redshift distribution and multiplicative shear correction, both stated in Table 1.We improve the constraints on  8 by 23% and on Ω m by 47% compared to the   -only case.This shows how powerful a combined analysis of the second-and third-order shear statistics is.The constraints on   IA and  ba are basically untouched, showing that M 3 ap is not helpful for constraining these nuisance parameters.
Compared to the maximum of the one-dimensional marginal distributions constraints of   measurements given in A21 ( 8 = 0.758 +0.017 −0.026 and Ω m = 0.253 +0.088 −0.074 ), we have slightly larger constraints when using only the   .This is because we use a numerical covariance, which is larger than the analytical one due to the underlying cosmology and the fact that we model the reduced shear effect.Furthermore, we note that A21 varied Ω cdm ℎ 2 , Ω b ℎ 2 and ℎ, while we fixed ℎ and Ω b and varied only Ω m .We also find that if we use the same pipeline as A21 but allow larger Ω cdm ℎ 2 , the posteriors increase towards larger Ω m and therefore get more consistent with our results.Furthermore, we use a different sampler compared to A21 and different baryon feedback Here the   are compared to M 3 ap using all available combinations for the aperture filter radii  ap ∈ {4 ′ , 8 ′ , 14 ′ , 32 ′ } and redshift bin combinations.We note that the prior (see Table 2) range on Ω m is enlarged compared to the previous validation plots.The FoM values of the   -only or M 3 ap -only case should not be compared to Fig. 6, as priors of Ω m bind their contours.
process modelling.Nevertheless, our results from the   analysis are consistent with A21 within 0.03  in Ω m and within 0.14  in  8 .Similar to A21, we also perform an internal consistency check, removing one  ph -bin at a time and finding consistent results.We discuss this in more detail in Appendix A, and we find that (as expected) the fifth  ph -bin is most important to constrain Ω m and  8 .We further discuss in Appendix A modelling checks regarding the infusion of IA, baryon feedback and the reduced shear correction.The baryon feedback and reduced shear correc-tions are always sub-dominant compared to a shift in  8 .The IA, however, is important, especially if lower  ph -bins are included and must be accurately modelled.
Finally, we notice that all statistics are consistent with  IA = 0, although both statistics alone seem to favour positive  IA .Interestingly, the joint analysis is shifted to lower  IA and is slightly more constraining, which we do not observe in the validation in Sect.8 where all posteriors accurately peak at the input value and all constraining power comes from the   .This Notes.The projected MP and their 68% confidence intervals result from the MCMC chains shown in Fig. 9.The values corresponding to the best  2 might differ slightly.We fixed ℎ = 0.6898,  0 = −1 and  s = 0.969 but marginalised over the ⟨⟩ and -bias uncertainties.The -values were estimated by counting the T17 mock data that give larger  2 than the real  2 relative to the total number of mock data.
might indicate that for real data with non-zero  IA , third-order shear statistics can contribute (at least a bit) to constraining IA, in line with predictions by, for instance, Pyne & Joachimi (2021) and Troxel & Ishak (2012).Furthermore, the baryonic parameter can be confined to  ba < 1.4 at 68% confidence with our statistic.
Here, we should note that M 3 ap does not help in constraining  ba , which is probably due to the fact that changes in  8 absorb all  ba effects.

Conclusions
This work validates the combined modelling of second-and thirdorder shear statistics, showing that its accuracy is well-suited for a KiDS-1000 analysis.Our second-order shear statistic of choice is the   -modes of the COSEBIs (Schneider et al. 2010) and the third-order statistic is M 3 ap (Schneider et al. 2005).In particular, we incorporated intrinsic alignment modelling based on the non-linear alignment model of Bridle & King (2007) and validated its accuracy against simulations infused with IA effects.This test is also interesting for other simulation-based analyses for which IA cannot be modelled analytically.We incorporated the impact of the baryonic feedback process by measuring a response function using the Magneticum simulations.Since the amplitude of this response function has no physical meaning, it is considered a nuisance parameter, which does not bias our cosmological parameter predictions.
We investigate which parts of the data vector can be neglected without losing too much cosmological information.This is important because the data vector for third-order shear statistics, due to the possibility it offers of combining three different filters with three different redshift bins, inflating the data vector to several hundred elements easily.Therefore, both a numerical and an analytical covariance are difficult to compute.We find that cross-tomographic redshift bins contain a large amount of cosmological information.Using only equal-scale filter radii but all available tomographic bin combinations was the best data compression strategy, which comes with the advantage that equalscale filters are faster to compute analytically as they require a lower integral accuracy.Next, we investigated the chosen filter radii.For this we measured and modelled the M 3 ap for  ap ∈ {4 ′ , 6 ′ , 8 ′ , 10 ′ , 14 ′ , 18 ′ , 22 ′ , 26 ′ , 32 ′ , 36 ′ }.We find that filter radii above 30 ′ do not contribute much constraining power and that having many small filter radii is unnecessary.The best option is to use filter radii  ap ∈ {4 ′ , 8 ′ , 14 ′ , 32 ′ }.We also tested whether selecting elements that maximise the Fisher information matrix on  8 is more optimal but find no difference compared to selecting all equal-scale filter radii.However, this method can be used to speed up the modelling and measurement of future analyses by discarding irrelevant elements.
Next, we validated the assumption that the correlation function estimator gives accurate results.This is important since only correlation function estimators have the potential to give unbiased results for real data.We used realistic mock data created from the T17 simulations.In particular, we created several galaxy catalogues where the positions of the galaxies are exactly at the KiDS-1000 galaxy positions.We had to rely on correlation functions to measure the second-and third-order statistics, which give unbiased results also if the data has a complex topology.Our first finding is that our modelling and measurement result in unbiased cosmological parameters given the KiDS-1000 uncertainty.Second, we find that using the reduced shear, or the shear itself, does not change the results and can, therefore, be ignored for the KiDS-1000 data analysis.
We conclude this paper with an analysis of the real KiDS-1000 data.Overall, our   constraints are less informative than the original KiDS-1000 analysis.This is mostly because we used a numerical covariance matrix from T17 simulations.However, the chosen sampler, the priors on the cosmological parameters, and the modelling strategy of the baryonic feedback processes impact our constraints, too.We find an  8 = 0.772 ± 0.022 and an Ω m = 0.248 +0.062 −0.055 , which are improved compared to the  only case by 23% and by 47%, respectively.With a -value of 0.25, we also find a good agreement of model and data given the KiDS-1000 uncertainty.This demonstrates that combining second-and third-order statistics is powerful in constraining cosmological parameters.The gain in constraining power in Ω m is also interesting for combined weak lensing and galaxy clustering analysis because the constraining power in Ω m for clustering analysis comes with the issue of further nuisance parameters such as galaxy bias.However, since second-and third-order shear statistics constrain Ω m quite well, combining it with clustering statistics might enable us to learn more about these nuisance parameters.
We leave the optimisation of the IA and baryonic feedback modelling for future analysis.An especially interesting improvement would be a more physically motivated description of the baryon feedback processes, which are identical for power and bispectrum.As all baryon feedback models rely on hydro simulations, we have to use the same simulations for power and bispectrum.Furthermore, although we found that the reduced shear and limber approximation is sufficient for a KiDS-1000 analysis, we probably have to model these effects for future Stage IV surveys.Lastly, we ignore the effect of source clustering for this work.Although Gatti et al. (2023) found it relevant for thirdorder weak lensing statistics based on convergence mass maps, we expect it to be less critical for our analysis, which uses shear catalogues and no mass map reconstruction.For future Stage IV surveys, this needs to be investigated.ap vector for all combinations of aperture filter radii  ap ∈ {4 ′ , 8 ′ , 16 ′ , 32 ′ }.The data vector results from one full-sky convergence realisation without shape noise.The grey band shows the expected uncertainty of KiDS-1000 estimated with the convergence maps.The M 3 ap is scaled by the third root of the product of the corresponding filter radii.

Fig. 1 :
Fig. 1: Redshift distribution of the five tomographic  ph -bins of the KiDS-1000 data.
The cosmo-SLICS+IA simulations, described in Harnois-Déraps et al. (2022), were used to test the modelling of IA; while the Magneticum lensing simulations, first introduced in Hirschmann et al. (2014), were used to infuse Baryon feedback on the model

Fig. 2 :
Fig.2: Measured and modelled   vector for the first five moments in the T17 mock data.The green and blue dots are the mean of all 1944 KiDS-1000 mock data that are also used to compute the covariance matrix with a resolution  pix = 0.18 arcmin 2 and shape noise.The red dots represent the data vector measured from one full-sky T17 realisation with a resolution  pix = 0.18 arcmin 2 and no shape noise.The grey band indicates the expected uncertainty from the KiDS-1000 survey.

Fig. 4 :
Fig. 4: Posterior distribution for data vectors and covariance measured from the T17 convergence maps catalogues.Here only specific parts of the M 3 ap data vector are used.Here 'only equal  ph -bins' means that only the auto tomographic bins are used.Lastly, all non-equal filter radii are discarded for 'only equal filter radii'.

Fig. 5 :
Fig. 5: Posterior distribution for different filter radii combinations if the data vector and covariances are measured from the T17 convergence maps.

Fig. 6 :
Fig. 6: Posterior distribution for data vectors and covariance measured from the T17 galaxy shear catalogues.The covariance matrix and reference data vector were measured from 1944 noisy  mock data.

Fig. 7 :
Fig. 7: Measured and modelled   for the first five components.The blue points show the measurements from the KiDS-1000 data A21.The blue error bars indicate the KiDS-1000 uncertainty.The different dashed lines show analytical descriptions at the MP if the combination of   + M 3 ap or only   is used.

Fig. 8 :
Fig. 8: Measured and modelled M 3 ap with filter radii  ap ∈ {4 ′ , 8 ′ , 14 ′ , 32 ′ }.The blue error bars indicate the KiDS-1000 uncertainty.The different dashed lines show analytical descriptions at the MP if the combination of   + M 3 ap or only M 3 ap is used.

Fig. 9 :
Fig. 9: Posterior distribution for the real KiDS-1000 data vector while the covariance is measured from the 1944 T17 galaxy shear catalogues.

Fig. A. 3 :
Fig. A.3: Illustration of several effects that change the   model vector.The fiducial model  fid  is without reduced shear,  IA = 0,  ba = 0,  8 = 0.78 and Ω m = 0.25.For the different points, one of the effects is changed.We scaled the model differences by the KiDS-1000 uncertainty, indicating the relevance of the change in the context of this work.

Fig. A. 4 :
Fig. A.4: Same as Fig. A.3 but for the M 3 ap model vector.

Fig. B. 2 :
Fig. B.2: Measured data and model M 3ap vector for all combinations of aperture filter radii  ap ∈ {4 ′ , 8 ′ , 16 ′ , 32 ′ }.The data vector results from one full-sky convergence realisation without shape noise.The grey band shows the expected uncertainty of KiDS-1000 estimated with the convergence maps.The M 3 ap is scaled by the third root of the product of the corresponding filter radii.

Table 1 :
Overview of the observational KiDS-1000 data.

Table 2 :
All varied parameters and their flat prior knowledge.Uniformly distributed priors on the parameters used in our cosmological inferences.The priors given in Table1on the multiplicative shear -bias and photometric redshift errors ⟨⟩ are used both for the real data analysis and for the mock data analysis, where for the latter the expectation values are set to zero.The ⟨⟩ for the sources follow a joint normal distribution with covariance matrix   ⟨⟩ shown in figure6of H21.The -bias follows an uncorrelated normal distribution.We increased the upper prior range on Ω m to 0.75 for the real data analysis to have   posteriors that are not bound by the prior.The  IA of the real data analysis are narrower to improve the emulator accuracy.The upper bound of  ba is chosen arbitrarily and does play a role only for M 3 ap -only analysis. Notes.

Table 4 :
MP values with their marginal 68% credible intervals.