KiDS-SBI: Simulation-Based Inference Analysis of KiDS-1000 Cosmic Shear
Maximilian von Wietersheim-Kramsta Kiyam Lin Nicolas Tessore Benjamin Joachimi Arthur Loureiro Robert Reischke Angus H. Wright
 


   

   

   Department of Physics and Astronomy, University College London, Gower Street, London, WC1E 6BT, UK
        Department of Physics, Institute for Computational  Cosmology, Durham University, South Road, Durham DH1 3LE, UK
        Department of Physics, Centre for Extragalactic Astronomy, Durham University, South Road, Durham DH1 3LE, UK
         The Oskar Klein Centre, Department of Physics, Stockholm University, AlbaNova University Centre, SE-106 91, Stockholm, Sweden
        Astrophysics Group and Imperial Centre for Inference and Cosmology (ICIC), Blackett Laboratory, Imperial College London, London SW7 2AZ, UK
        Argelander-Institut für Astronomie, Auf dem Hügel 71, D-53121 Bonn, Germany
        Ruhr University Bochum, Faculty of Physics and Astronomy, Astronomical Institute (AIRUB), German Centre for Cosmological Lensing, D-44780 Bochum, Germany
             

   
 
  We present a simulation-based inference (SBI) cosmological analysis of cosmic shear two-point statistics from the fourth weak gravitational lensing data release of the ESO Kilo-Degree Survey (KiDS-1000). KiDS-SBI efficiently performs non-Limber projection of the matter power spectrum via 's method, and constructs log-normal random matter fields on the curved sky for arbitrary cosmologies, including effective prescriptions for intrinsic alignments and baryonic feedback. The forward model samples realistic galaxy positions and shapes based on the observational characteristics of KiDS-1000, incorporating shear measurement and redshift calibration uncertainties, as well as angular anisotropies due to variable survey depth and point-spread function variations. To enable direct comparison with standard inference, we limit our analysis to pseudo-angular power spectra as summary statistics. The SBI is based on neural density estimation of the likelihood with active learning to infer the posterior distribution of spatially-flat ΛCDM cosmological parameters from 18,000 realisations. We infer a mean marginal of the growth of structure parameter S_8≡σ_8 (Ω_m / 0.3)^0.5 = 0.731± 0.033 (68 %). We present a measure of goodness-of-fit for SBI and determine that the forward model fits the data well with a probability-to-exceed of 0.42. For fixed cosmology, the learnt likelihood is approximately Gaussian, while constraints widen compared to a Gaussian likelihood analysis due to cosmology dependence in the covariance. Neglecting variable depth and anisotropies in the point spread function in the model can cause S_8 to be overestimated by ∼5%. Our results are in agreement with previous analysis of KiDS-1000 and reinforce a 2.9 σ tension with early-Universe constraints from cosmic microwave background measurements. This work highlights the importance of forward-modelling systematic effects in upcoming galaxy surveys, such as Euclid, Rubin and Roman.

   

   KiDS-SBI: Simulation-Based Inference Analysis of KiDS-1000 Cosmic Shear
    Maximilian von Wietersheim-Kramsta
          1,2,3
          Kiyam Lin1
          Nicolas Tessore1
          Benjamin Joachimi1
          Arthur Loureiro4, 5
          Robert Reischke6,7
          Angus H. Wright7
          
    Received xx; accepted xx
====================================================================================================================================================================================================================================






§ INTRODUCTION

Cosmic shear, the weak gravitational lensing effect on distant galaxies due to matter in the foreground, is a powerful tool to study the distribution of matter in the Universe and to probe its large-scale structure. By measuring the distortions in the shapes of galaxy images caused by the gravitational influence of intervening matter, we can infer a combination of the matter density and the amplitude of the matter power spectrum within the framework of a spatially flat cold dark matter cosmology with a cosmological constant, ΛCDM. Recent cosmic shear analyses of the data taken in Stage-III galaxy surveys, such as the Kilo-Degree Survey[<https://kids.strw.leidenuniv.nl/>] (KiDS; ), the Hyper Suprime-Cam survey[<https://hsc.mtk.nao.ac.jp/ssp/>] (HSC; ), and the Dark Energy Survey[<https://www.darkenergysurvey.org/>] (DES; ), have constrained these cosmological parameters with unprecedented precision. Upcoming Stage-IV galaxy surveys, such as Euclid[<https://www.euclid-ec.org/>] <cit.>, Rubin [<https://www.lsst.org/>] <cit.> or Roman[<https://roman.gsfc.nasa.gov/>] <cit.>, will further improve upon these constraints, as these surveys will be wider, deeper and more precise. However, many analysis challenges still remain.

Current Stage-III galaxy survey measurements of the root-mean square of the matter overdensity field at 8 Mpc h^-1, σ_8, or its analogue S_8 ≡σ_8 (Ω_m / 0.3)^0.5, agree well with each other <cit.> despite using largely independent methodologies and survey volumes. However, the constraints from observations of the late Universe <cit.> disagree by up to ∼3.4σ with the σ_8 value consistent with the early-Universe observations from the Cosmic Microwave Background (CMB; ). Recently, cluster abundances as measured by SRG/eROSITA <cit.> have provided another measurement of the fluctuation amplitude even higher than the ones of the CMB, although similar recent cluster studies have obtained S_8 constraints more consistent with weak lensing and CMB analyses <cit.>. This discrepancy might be pointing at new physical phenomena, but it could also be caused by unconsidered systematic effects when modelling the low-redshift large-scale structure or CMB signal and noise (see e.g. ). Such systematic effects may be physical, e.g. baryonic feedback or the intrinsic alignments of galaxies <cit.> or observational.

We aim to shed some light onto this by conducting a full cosmic shear analysis using neural likelihood estimation simulation-based inference (SBI) with the same complexity as current Stage-III analyses. SBI, also known as Likelihood-Free Inference or Implicit Likelihood Inference, is a Bayesian inference method which does not require an explicit formulation for the likelihood function of the data given the parameters of interest. Instead, the likelihood is implicitly calculated by evaluating the joint probability of data and parameters of interest from forward simulations which map the parameters to the corresponding mock data vectors. This comes with multiple advantages with respect to other standard approaches which require an explicit form for the likelihood. Firstly, the likelihood is allowed to take an arbitrary form, so one can avoid the typical assumption of a Gaussian likelihood or avoid having to define a complex analytical expression for the likelihood. Secondly, for some models and measurements, it may not even be possible to define an analytical likelihood or it would be too resource-intensive, e.g. the covariance of n^th-order statistics depends on up to 2n^th-order correlation functions which become geometrically more expensive to compute. In such cases, as long as the observables can be simulated, an effective/implicit likelihood may be found using SBI. Lastly, in standard analyses, it often is required to validate the signal and noise modelling with forward simulations, so one can make use of existing tools to perform SBI.

When running a Markov-Chain Monte Carlo (MCMC) chain to sample the posterior distributions from a Gaussian likelihood, one either has to compute a numerical covariance matrix from forward simulations or create an analytical model for the covariance which ought to be validated with a numerical one. In both cases, for a data vector, 𝐝, one would require a number of simulations >|𝐝| to get an accurate sample covariance matrix. This is already the approximate amount of forward simulations needed for an SBI analysis <cit.>. Therefore, SBI allows one to do full Bayesian uncertainty propagation from data to parameters for any model which can be simulated, without substantial additional computational costs when compared to standard approaches at comparable accuracy.

In cosmic shear, we know that the commonly used two-point statistics to quantify the shear-shear correlations of galaxies on the sky have an approximately Gaussian likelihood, but this assumption has its limitations <cit.>. Additionally, when observational and physical systematics are included, the likelihood may become even less Gaussian <cit.>. In such cases, like the one presented in this paper, SBI allows one to capture potential non-Gaussanities and propagate them to the posterior distributions.

SBI can also compute effective likelihoods for summary statistics for which analytical methods may be intractable, such as higher-order cosmic shear statistics <cit.>. In this work, we limit the scope to only angular power spectra of cosmic shear and leave higher-order statistics or field-level statistics within SBI as an avenue for future work. This allows us to make a direct comparison to previous standard cosmic shear analyses (e.g. ), while also harnessing the efficiency of log-normal forward models of matter density fields which are mostly accurate for two-point statistics <cit.>.

A potential cause for deviations from Gaussianity in the likelihood of two-point statistics/angular power spectra is the anisotropy in the observational depth of a galaxy survey across the sky <cit.>. As a survey such as KiDS images the sky, each pointing and each pixel are subject to variable environmental/background conditions which may affect the selection of galaxies (e.g. atmospheric conditions, thermal expansion of the telescope, moonlight, zodiacal light, Galactic extinction, etc.).

Due to the given survey strategy and/or dithering to compensate for gaps in the telescope's charge-coupled device (CCD), it is possible that there is spatial variation in the overlap between pointings. When the images of each pointing are subsequently combined, some galaxies are observed more often than others which can cause variations in the signal-to-noise ratio which are spatially correlated.

In effect, all these variations lead to fluctuations in the effective galaxy density which are observed across the footprint, while also selecting galaxies according to their shape. This has been found to be a percent-level effect on the cosmic shear signal measured by KiDS <cit.>, while also affecting the standard deviation of the KiDS-1000 cosmic shear signal (by ≲ 20%; ). Although for KiDS-1000 this is within the statistical uncertainty, it will become important for upcoming galaxy shear surveys such as Euclid, Rubin or Roman. Additionally, the effect of spatial variability on the non-Gaussianity of the likelihood has not been explored in previous work. To address these needs, in this work, our forward model includes the effects of the observational depth variations in KiDS, so that we may propagate the effects of spatial variations in the depth to the uncertainty of inferred cosmological parameters.

 To conduct the SBI, we choose to make use of Density Estimation Likelihood-Free Inference (DELFI, ), as it offers good performance and scales better when increasing the dimensionality of the parameter space of interest than other SBI methods <cit.>, such as Approximate Bayesian Computation (ABC; ). DELFI discards less information from the forward simulations than ABC. ABC only keeps the forward simulations which have a certain level of agreement with the desired data vector. In contrast, DELFI learns a probability density distribution of the data vectors as a function of the model parameters based on all the forward simulations. To achieve this, DELFI employs ensembles of neural density estimators in order to learn the sampling distribution of the data and the parameters from the forward simulations.

 Of course, as with any inference method, one needs to be mindful of any intrinsic biases which might come with SBI, how to mitigate them and avoid them. As the likelihood is not known analytically, it is important that the chosen SBI method, sampling scheme and number of simulations are sufficient to ensure that the learnt likelihood has converged to the model's true likelihood function. To avoid biasing the inference, this requires careful robustness and convergence tests of the implicit likelihood which may not always be necessary for a standard Gaussian likelihood analysis. Hence, we dedicate a considerable effort in this work to ensure the inference pipeline is validated, the neural density estimators of DELFI have converged, the inferred posterior is accurate, and the inferred parameters give a good fit to the data.

SBI is increasing in popularity, among other reasons, because it requires about an order of magnitude fewer evaluations of the likelihood than a standard Gaussian likelihood MCMC analysis to constrain the posterior <cit.>. We highlight the methodology of <cit.> which conducted a similar analysis to the one outlined in this paper with the DES Year 3 weak lensing data; see <cit.> for results. Whereas the DES analysis focuses on SBI of two-point and higher-order shear statistics, as well as cosmologies beyond ΛCDM based on a fixed number of N-body simulations while incorporating standard systematics, the analysis in this paper exclusively considers two-point statistics within the context of ΛCDM. We highlight that in our analysis the data's likelihood is learned by actively sampling as many efficient forward simulations as needed which capture many physical and observational systematics, including anisotropic systematics which are typically not considered in weak lensing analyses. With this, we aim to give more insight into the robustness against systematics of standard weak lensing analyses.

In a previous paper <cit.>, we introduced a new inference pipeline for estimating the cosmological parameters from the cosmic shear data of KiDS-1000 <cit.> using DELFI. In , the SBI pipeline for KiDS is shown to be robust, accurate and efficient, even when constraining a 12-dimensional posterior distribution. However, in that work, the simulated vectors were based on random samples from the covariance matrix used in the fiducial Gaussian KiDS-1000 analysis, so the likelihood was Gaussian by construction to allow for direct comparison with the standard formalism. In this paper, we present a novel suite of physically motivated and realistic forward simulations of the cosmic shear observables as seen by KiDS, so the neural density estimators learn any non-Gaussianities in the likelihood induced by our model, plus any complex relations between the data and the model parameters. We then use these simulations to create two inference pipelines: a fiducial one which makes the same modelling choices as the KiDS-1000 analysis <cit.> and a more realistic one which explicitly models systematic effects that would be prohibitively difficult to incorporate into a standard likelihood analysis. From each pipeline, we obtain a separate posterior estimate for S_8 which we then scrutinise with a novel measure of goodness-of-fit for SBI posteriors.

This paper has the following structure: Sect. <ref> presents the KiDS data used for this analysis. Section <ref> gives a detailed description of the forward simulations developed for the SBI analysis. Section <ref> details how the galaxy shapes and summary statistics are measured in both the simulations and the data. Section <ref> describes the setup of the simulation-based inference pipeline using DELFI based on . Section <ref> presents the novel measure of goodness-of-fit for SBI that we use in this paper. In Sect. <ref>, we conduct consistency tests to ensure the SBI pipeline is unbiased and robust, while Sect. <ref> discusses the validity of a Gaussian likelihood assumption. Section <ref> presents the cosmological inference results from the KiDS-1000 weak lensing data. Lastly, we give concluding remarks in Sect. <ref>, and provide additional details about the analysis in the Appendix.








§ KIDS-1000 DATA
 
The Kilo-Degree Survey (KiDS) is a large public galaxy survey conducted by the European Southern Observatory using the OmegaCAM CCD mosaic camera <cit.> which is attached to the 2.6 m VLT Survey Telescope (VST). The survey covers approximately 1,350 deg^2 between two distinct fields known as KiDS North and South which are fields located in the Northern and Southern Galactic Caps (with DEC ≈ 0º and DEC ≈ -30º, respectively). This area is the same as the one covered by the VISTA Kilo-degree INfrared Galaxy survey (VIKING; ), which means that both surveys together observe every object with a total of nine photometric bands: ugriZYJHK_ s. The analysis presented in this paper makes use of DR4 <cit.>, also known as KiDS-1000, which covers approximately 1,000 deg^2. The effective area covered by KiDS-1000 shape measurements shown in Fig. <ref> extends over 773.3 deg^2 <cit.>. This area is calculated directly from the mosaic mask defined a the native OmegaCAM pixel scale of 0.213 arcsec. 

KiDS images are processed using  <cit.> and the methodology outlined in <cit.> to join the VST and VIKING photometric measurements. To obtain the shape measurements,  <cit.> is used to process the r-band images, while the gravitational shear estimates are then obtained from  <cit.>. Depending on the shape noise variance and the ellipticity measurement noise variance,  assigns a weight to the shear measurement of each galaxy, w_i, which scales the shear signal such that the signal-to-noise ratio is optimal. Although KiDS-SBI currently does not model  weights (all galaxies are weighted equally by default), the weights are carried through the measurement pipeline when measuring the cosmic shear summary statistics from the real KiDS-1000 data. If the  weight is mischaracterised, this may lead to shear bias. For this purpose, the shape measurements in the KiDS-1000 gold sample are calibrated to determine any shear biases due to selection biases, noise, weight bias, point-spread function residuals or otherwise are calibrated <cit.>. The resulting estimate of the linear multiplicative shear bias per tomographic bin is then used to unbias the galaxy shear measurements in KiDS-1000. All details on the reduction of the images and the shape calibration are given in <cit.>.

The photometric redshifts are estimated from the KiDS 9-band photometry using Bayesian template-fitting as incorporated within the  code <cit.>. The photometric redshifts are calibrated with spectroscopic redshift measurements made of galaxies in the KiDS sample <cit.> based on a subsample of 99% completeness. The mapping from photometric to spectroscopic redshift, P(z|z_ph), is learnt through a self-organising map <cit.>. As is the case for previous KiDS-1000 analyses, we use the KiDS-1000 "gold-sample" which only considers galaxies with reliable redshift and shape measurements. The maximum posterior photometric redshift is subsequently used to split the source galaxy catalogue into five tomographic bins at the following boundaries: (0.1, 0.3, 0.5, 0.7, 0.9, 1.2]. These then give the redshift distributions shown in Fig. <ref>. As the cosmological constraints are sensitive to shifts in the redshift distributions, the KiDS-1000 analysis also allows the shift in the mean of each tomographic bin to vary freely in the cosmological inference as a nuisance parameter while incorporating for the correlations between them <cit.>. This not only avoids biases in the cosmological parameter estimates but also provides an additional test on the redshift calibration. 




§ FORWARD SIMULATIONS

The forward simulations for cosmic shear analyses described in this paper are fully implemented and available in the KiDS-SBI[Kilo-Degree Survey – Simulation-Based Inference; <https://github.com/mwiet/kids_sbi>] module. This code is built within the Cosmological Survey Inference System ([<https://github.com/joezuntz/cosmosis>]; ) and it is based upon the KiDS Cosmology Analysis Pipeline ([<https://github.com/KiDS-WL/kcap>]; ). In addition, the KiDS-SBI pipeline is generalisable to other cosmic shear analyses. 





The overall outline of the forward simulations is shown in Fig. <ref>. Firstly, the simulation is fed a set of cosmological and astrophysical parameters, Θ, from a sampler of our choice. All of these parameters are passed into  <cit.> which calculates the three-dimensional matter power spectrum, P_δ, nl, at a given cosmology (see also Sect. <ref> for more details). We incorporate baryonic feedback using  <cit.>. Next, the simulation splits the light-cone along the line-of-sight into concentric spherical shells centred at the observer, as required by our large-scale structure simulator of choice,  <cit.>. Such an approach is computationally efficient as all observables are simulated with a coarse resolution along the line-of-sight. Within each shell, we sample log-normal random matter fields to efficiently model the two-point statistics with high accuracy, while also incorporating some higher-order fluctuations. As an added benefit, this simplifies the implementation of effects such as redshift evolution as well as survey characteristics, as they usually also follow a spherical geometry. On the other hand, the concentric shells also allow to split the Universe into shells, which are of a scale large enough such that a log-normal field becomes an accurate description of the matter over-densities within it (see Sects. <ref> for details).

Once we know the geometry of our matter shells, we can project the three-dimensional matter power spectrum using the redshift distributions of these shells in order to obtain their two-dimensional angular power spectra, C_δδ. The latter describes the correlations between the matter over-densities in each shell. Since we also require accurate cross-correlations between different shells to compute the log-normal fields, the Limber approximation is not applicable here, as the cross-correlations are zero in that approximation with the shells being non-overlapping along the line-of-sight. To avoid this, we project the three-dimensional matter power spectrum with a non-Limber integral. This is usually an expensive computation, so we show here a novel implementation of the Levin method of integration <cit.> which allows for the efficient integration of the spherical Bessel functions found in non-Limber integrals (more details are given in Sect. <ref>).

We can then construct log-normal matter fields which are consistent with these correlations and their associated cosmology. To do this, we use the Generator for Large Scale Structure ([Generator for Large Scale Structure; <https://github.com/glass-dev/glass>], ), which efficiently generates correlated log-normal random matter fields, δ, for each shell. Then, it accurately computes the respective convergence field, κ, (see Sect. <ref>) while also incorporating an effective convergence field due to intrinsic galaxy alignments, κ_IA (see Sect. <ref>).

Upon the construction of the matter, δ, and convergence fields, κ, we use them to Poisson-sample galaxies that reside within the matter fields as well as their properties, i.e. intrinsic ellipticities and shear. During this step, with the aid of the  module <cit.>, we also consider survey characteristics, such as the survey footprint, its spatial variability in depth, in its redshift distributions, as well as in the shape and shot noise as described in Sect. <ref>. Thus, we obtain a full galaxy catalogue with shape information consistent with the input cosmology as well as with the previously mentioned survey properties.

The last steps of the forward simulation pipeline involve the post-processing and compression of the simulated data catalogue. In principle, any data vector could be obtained here, such as two-point shear statistics <cit.>, the shear fields <cit.> or even the full catalogues. However, in order to reduce the dimensionality of the data vector efficiently, while still retaining most information about the cosmology, similar to previous KiDS-1000 analyses <cit.>, we choose observed angular power spectra, a.k.a. pseudo-Cls <cit.>, as our data vector of choice. This allows us to compare our SBI analysis to the results from the standard analysis of KiDS-1000 data using pseudo-Cls in <cit.>. The KiDS-SBI pipeline is, in principle, set up to make cosmological inferences from cosmic shear based on any statistic which could be derived from a galaxy catalogue, but the underlying statistical random fields are currently only capable of achieving percent-level accuracy for two-point statistics. Any higher-order statistics measured from the current forward simulations are not necessarily accurate enough to conduct inference. The statistics and their post-processing are defined in Sect. <ref>.

The forward simulations are designed to model weak gravitational lensing observations on the level of galaxy catalogues. The galaxy populations are sampled such that they trace the underlying log-normal random matter fields, while the galaxy shapes are lensed in accordance with the lensing potential of the matter fields along the line-of-sight. As the simulations do not model images, any image-level systematic effects, such as shear biases or variable depth, are included as spatially varying probability density functions.




 §.§ Cosmology dependence: three-dimensional matter power spectrum

The basis of the forward simulations is the 3-dimensional matter power spectrum, P_δ, nl(k, z; Θ). In the linear regime, the equal-time three-dimensional matter power spectrum's cosmology dependence comes from the following relation

    P_δ, l(k, z; Θ) =  T^2(k; Θ)   D^2(z; Θ)   P_ℛ(k; Θ) ,

where k is a wavenumber, z is redshift, Θ is the set of cosmological parameters, P_ℛ(k; Θ) is the primordial density fluctuation power spectrum, D is the growth factor of structure, T is the transfer function between the matter overdensity field, δ(𝐤, z), and the primordial curvature fluctuation field, ℛ(𝐤; Θ), such that δ(𝐤, z) = T(𝐤; Θ)  D(z; Θ) ℛ(𝐤; Θ). To compute this, we make use of  <cit.>. To stay in line with the main KiDS-1000 analysis <cit.>, we assume a normal neutrino hierarchy, while also assuming a fixed sum of neutrino masses of ∑ m_ν = 0.06 eV/c^2.

We compute the non-linear matter power spectrum via a non-perturbative model,  <cit.>. Although updated iterations of this model exist <cit.>, we limit ourselves to , so our model remains comparable to the KiDS-1000 cosmic shear analysis.  uses a halo model approach to incorporate the effects of baryonic feedback on the matter distribution. The main driving factor in this is the baryonic matter expelled by Active Galactic Nuclei (AGN). This suppresses the power at small scales, i.e. large k, as a function of the amplitude of the halo mass-concentration relation, A_bary, and the halo bloating parameter, η_0. We also fix the relation between these two parameters to η_0 = 0.98 - 0.12 A_bary in accordance with <cit.> and <cit.>, and treat A_bary as the only free parameter related to baryonic feedback.

We note that this treatment may lead to systematic biases at small scales, as the galaxies will be sampled from a matter power spectrum which is already modified by baryonic feedback. It would be more physically accurate to sample the galaxies from a matter power spectrum without baryonic non-linearities and subsequently add a perturbation to the matter density contrast fields due to baryonic feedback (similar to the approach used for the intrinsic alignments described in Sect. <ref>). Nevertheless, we do not expect this to cause a large discrepancy for the scales which we are probing in the cosmic shear signal (ℓ < 1500). Field-level implementations of baryonic feedback already exist <cit.>, and they would constitute interesting avenues for future extensions to KiDS-SBI. 



 §.§ Working on the sphere


Astronomical observations can be viewed as images of a three-dimensional space projected onto the inside of a two-dimensional surface of a sphere at a given distance with the observer at its centre. Therefore, galaxy position measurements, most observational biases and some systematics are naturally expressed in spherical coordinates on the sky. For this reason, the forward model presented in this work is expressed in the same geometry.

The KiDS-SBI forward simulations use  <cit.> to model the underlying large-scale structure. Within , we choose to model the large-scale structure through log-normal random matter fields within non-overlapping concentric shells centred at the observer. Therefore, a given shell is a comoving volume which spans the full sky and has a finite width along the line-of-sight from one redshift, z_i, to another, z_i+1, such that i ∈{1, 2, 3, ..., N_shells}. We define the matter weight function W^(i)(z) describing the distribution along the line-of-sight of the comoving volume spanned by a given shell as <cit.>,

    W^(i)(z; Θ) =  f_k^2(z; Θ)/E(z; Θ),    if z_i≤ z < z_i+1 ,
     0,    otherwise ,

where E(z; Θ) is the dimensionless Hubble function and f_k(z; Θ) is the transverse comoving distance for a given redshift and cosmology, which is defined as

    f_k(z; Θ) = c/H_0√(Ω_k)sinh(√(Ω_k) H_0/cχ(z; Θ) ),    if Ω_k > 0 ,
    χ(z; Θ),    if Ω_k = 0 , 
    c/H_0√(|Ω_k|)sin(√(|Ω_k|) H_0/cχ(z; Θ) ),    if Ω_k < 0 ,

where H_0 is the Hubble constant, Ω_k is the curvature density parameter and c is the speed of light in a vacuum. As noted in <cit.>, different choice of window function may induce percent-level biases in the resulting shear power spectrum. We test for this in Appendix <ref>, and find sub-percent level agreement between the measured shear from the simulations with this choice of window function and the shear signal expected from theory (see Appendix <ref>). 

A given shell should not be thicker than a few hundred Mpc of comoving distance along the line-of-sight for log-normal random fields to be an accurate enough representation of the distribution of large-scale structure <cit.>. Otherwise, discretisation effects would smooth away large-scale structure <cit.>. Additionally, one needs to be mindful of computational resources, because the number required of angular power spectra between the matter fields within each shell scales as 2N_shells. Simultaneously, the simulations' run-time for a given set of cosmological parameters increases with N_shells. To balance these factors, we aim to have as many shells along the line-of-sight as computationally feasible. 

For an accurate weak lensing simulation at the relevant redshifts, 0.6 ≲ z ≲ 1.2, we require shells that finely sample the lowest redshifts. The lensing signal in a given matter shell depends on the weighted sum of the matter fields of the shells with lower redshifts. At redshifts where lensing becomes more important as the lensing efficiency increases, prioritising thinner shells helps to reduce discretisation effects. We find that the following set of redshift limits for 19 shells are sufficient for the resolution needed for KiDS-1000, while still being efficient: {0,   0.04,   0.08,   0.12,   0.16,   0.2,   0.27,   0.34,   0.42,   0.5,   0.58,   0.66, 0.75,   0.88,   1.03,   1.19,   1.36,   1.55,   1.76,   2 } (see Fig. <ref>).




 §.§ Non-Limber projection

Although the three-dimensional matter power spectrum characterises all cosmological dependence of the matter fields in our simulation, for it to define the cosmology dependence in the simulated matter fields, it has to be projected according to the geometry of the shells considered in Sect. <ref>. To determine the correlations between the matter fields of comoving spherical shells, we define the two-dimensional angular matter power spectra for a given cosmology, Θ, and a given shell combination (i j) as follows

    C_δδ, ℓ^(i j)(Θ) = ⟨δ^(i)_ℓ m(Θ)  δ^(j)*_ℓ m(Θ) ⟩ ,

where δ^(i)_ℓ m(Θ) represents the harmonic coefficients defined through a spherical harmonic transform of the projection of the spin-0 matter field, δ^(i), as given by <cit.>,

    δ^(i)_ℓ m(Θ) = ∫d^2θ δ^(i)(θ; Θ) 0*ℓ m(θ) ,

where 0ℓ m (θ) are the spherical harmonics for a spin value of 0 as a function of the two-dimensional spatial sky position θ. The projected matter field, δ^(i)(θ; Θ), can in turn be defined with respect to the underlying three-dimensional matter overdensity field, δ(θ, z; Θ), as

    δ^(i)(θ; Θ) = ∑_ℓ mδ^(i)_ℓ m(Θ)0ℓ m(θ) = ∫dz   W^(i)(z; Θ)  δ(θ, z; Θ) ,

where W^(i)(z; Θ) is the weight function of the given shell as defined in Eq. (<ref>).

The three-dimensional matter overdensity field's cosmology dependence can be quantified as follows

    ⟨δ(𝐤, z; Θ)  δ^*(𝐤', z'; Θ) ⟩ = (2π)^3δ_D(𝐤-𝐤')   P_δ, nl(k, z, z'; Θ) ,

where δ(𝐤, z; Θ) is the Fourier transform of the two-dimensional matter overdensity field at z, δ(θ, z,; Θ), while δ_D(𝐤-𝐤') is the Dirac delta function.

If we combine the relations in Eqs. (<ref>) to (<ref>) and, for the sake of simplicity, assume that Ω_k=0, i.e. f_k(z; Θ) = χ(z; Θ), we can characterise the cosmology dependence of C_δδ^(i j)(ℓ; Θ) directly with respect to the three-dimensional matter power spectrum, P_δ, nl(k, z, z'; Θ), using the following relation

    C_δδ^(i j)(ℓ; Θ) = 2/π∫dχ  W^(i)(z[χ]; Θ) ∫dχ'   W^(j)(z[χ']; Θ) 
    ∫dk   k^2   P_δ, nl(k, z[χ], z[χ']; Θ)   j_ℓ(k χ)   j_ℓ(k χ') ,

where χ'≡ z[χ'], j_ℓ(kχ) are spherical Bessel functions of order ℓ, and ℓ∈ℤ^0+. We can simplify Eq. (<ref>) by making the geometric approximation, i.e. taking the geometric mean of two equal-time power spectra <cit.>, such that

    P_δ, nl(k, z[χ], z[χ']; Θ) = √(P_δ, nl(k, z[χ]; Θ)   P_δ, nl(k, z[χ']; Θ)) ,

which has been shown to be an accurate approximation both in the linear and non-linear regimes <cit.>.



Additionally, it is common to make the so-called Limber approximation <cit.> which consists in approximating a spherical Bessel function as a single peak or delta function at kχ = ℓ + 1/2. This approximation is reasonably accurate in the auto-correlations for broad window functions and high angular modes, C_δδ^(i i)(ℓ; Θ),  but depending on the thickness of the shells, the Limber approximation can lead to substantial biases, especially, at large scales. However, when making the Limber approximation, all off-diagonal angular power spectra, C_δδ^(i j), where i ≠ j, will be zero as long as their associated weight functions, W^(i)(z[χ]; Θ), are non-overlapping, as is the case for the weights constructed in Eq. (<ref>). This can be seen by looking at the full Taylor expansion of Eq. (<ref>) as given by <cit.>

    C_δδ^(i j)(ℓ) = ∫dχ/χ^2  W^(i)(z[χ])   W^(j)(z[χ])   P_δ, nl(ℓ +1/2/χ, z[χ] ) 
    { 1 - 1/(ℓ + 1/2)^2[ χ^2/2( W^(i)”(z[χ])/W^(i)(z[χ])  + W^(j)”(z[χ])/W^(j)(z[χ]))  
     + 
        χ^3/6( W^(i)”'(z[χ])/W^(i)(z[χ])  + W^(j)”'(z[χ])/W^(j)(z[χ])) ] + 𝒪([ℓ + 1/2]^-4) } ,

where each apostrophe, ', denotes a partial derivative with respect to χ, i.e. ∂/ ∂χ, while we also omit the explicit cosmological dependence of the power spectrum and the weight functions here for the sake of clarity. The non-Limber integral therefore also contributes to the autocorrelations, where |i-j|=0, particularly for small ℓ and the narrow weight functions considered here. If not considered this can induce a large bias in the angular power spectra of the shells which could significantly bias cosmological constraints from KiDS-1000. At the same time, we find that for the cross-correlations of adjacent shells, such that |i-j|=1, the angular power can be of the order of 10^-1 C_δδ^(i i) or less. In this case, the simulated δ^(i)(θ; Θ) would have a percent-level bias while being void of any large-scale structure along the line-of-sight beyond the scales of individual shells. To avoid this, we choose to calculate the full non-Limber projection as given by Eq. (<ref>). The correlations between shell pairs which are not immediate neighbours, i.e. |i-j|>1, contribute less than 0.1% to the overall angular power of a shell, so their effect is undetectable for KiDS.

Calculating the full non-Limber integral is a computationally expensive endeavour as it involves three numerical integrals over highly oscillatory Bessel functions. To address this, we present the  module[; <https://github.com/rreischke/nonLimber_matter_shells>]. This is a novel Python module which implements the Levin integration method for oscillatory functions <cit.> in the context of non-Limber integrals for weak lensing and galaxy clustering which is discussed in more detail in Appendix <ref>.

The module allows one to divide the domain of the angular power spectra into three regions: a non-Limber domain, a second-order extended Limber domain and a Limber domain. In this work, we choose to split the calculation of C_δδ^(i j)(ℓ) into two domains: non-Limber (1 < ℓ≤ℓ_max, nL) and second-order extended Limber given in Eq. (<ref>) for ℓ_max, nL < ℓ≤ 30,000. Over the former domain,  numerically integrates the expression shown in Eq. (<ref>) for the given matter power spectrum and the weights. Over the latter domain, we perform the numerical integration of Eq. (<ref>) up to the second-order term in (ℓ+1/2)^-1. This reduces the computational resources needed for the integration without any substantial loss of accuracy as long as ℓ_max, nL is sufficiently large, since the residuals scale with (ℓ + 1/2)^-4 <cit.>. With the shells shown in Fig. <ref> weighted with Eq. (<ref>), we obtain angular power spectra, C_δδ^(i j)(ℓ; Θ), shown in Fig. <ref>.

With this configuration, we find that we can obtain computationally efficient and accurate (consistent with  within 0.1%) angular power spectra which describe the correlations between the matter fields within each spherical shell of the simulation. This is in line with testing conducted on the  module in previous applications <cit.>. 



 §.§ Log-normal matter field simulations

With the input angular power spectra for a given cosmology computed as described in Appendix <ref>, the next step in the forward simulations involves sampling random matter fields consistent with these correlations. These fields are generated with  (). This framework allows one to create efficient and accurate random field simulations of large-scale structure matter fields as well as the associated weak gravitational lensing signals. Within the  suite, many different choices can be made, while it also allows one to sample the angular positions, redshifts, shapes and shears of galaxies. Nevertheless, in this analysis, we only choose to use it to simulate the underlying matter overdensity fields, δ^(i)(θ; Θ), as well as the associated convergence fields for each shell, κ^(i)(θ; Θ). These are then integrated along the line-of-sight to give the convergence field for a source within a given shell, κ^(i).

It has been found that a good approximation to such fields is a log-normally distributed random field (which can be derived from a Gaussian random field). In fact, log-normal fields are even able to produce reasonable approximations to the three-point and four-point statistics measured from N-body simulations <cit.>. For this reason, log-normal random fields are a common approximation used for matter overdensity fields <cit.> as well as for convergence fields <cit.>. 

Log-normal fields within  are based on Gaussian random fields. One can obtain these by taking the angular power spectra of the shells and assuming that they describe the correlations between log-normal fields, δ_log^(i)(θ; Θ), as follows

    ⟨δ_log^(i)(θ; Θ)   δ_log^(j)*(θ; Θ)  ⟩ =  C_δδ^(i j)(θ; Θ) .

By defining the mean and the variance for G_δδ^(i j)(θ; Θ) as 

    ⟨δ_G^(i)(θ; Θ)   δ_G^(j)*(θ; Θ)  ⟩ =  G_δδ^(i j)(θ; Θ) ,

one can then sample δ_G^(i)(θ; Θ) for each shell. To achieve this with numerical methods, the two-dimensional spherical shell of each layer is discretised using  pixels <cit.>. This means that θ is discretised such that it can be mapped to a linearised variable given by θ_m where m ∈{1, 2, 3, ..., 12 N_side^2} and N_side is the  resolution parameter. To obtain the final matter overdensity fields, we use the following expression <cit.>,

    δ_log^(i)(θ_m; Θ) = e^δ_G^(i)(θ_m; Θ) - λ ,

where λ is the shift of the log-normal distribution which is assumed to be λ = 1 for matter fields. For more details about this calculation, see .

The layers of concentric log-normal matter overdensity fields within each shell give a full description of the large-scale structure. The resolution of this structure will of course be limited by the size of the shells along the line-of-sight, while being limited by the resolution of the discretised two-dimensional pixels on the sphere. However, as discussed with regards to the resolution of shells in Sect. <ref>, if one chooses sufficiently large ℓ_max up to which the input angular power spectra are calculated, one can accurately sample δ_log^(i)(θ_m; Θ) up to an N_side∼ 0.5 ℓ_max <cit.>.

As we intend to make use of these simulations for weak gravitational lensing, the next step in the forward simulations is to model how the image of galactic sources within a given shell will be distorted by weak gravitational lensing due to the matter overdensities the light encounters along the line-of-sight to the observer. Here, the geometry of the simulations aids us again, since the concentric volumes in comoving distance that each shell makes up intrinsically allow one to trace all possible light-cones emanating from the observer.

We start from the definition of the convergence for a source located at θ and z along a continuous line-of-sight under the Born approximation <cit.>,

    κ(θ, 𝑧; Θ) = 3 Ω_m/2∫_0^𝑧d𝑧'    𝑓_𝑘(𝑧'; Θ)   [𝑓_𝑘(𝑧'; Θ)-𝑓_𝑘(𝑧; Θ)]/𝑓_𝑘(𝑧; Θ)
       ×1+𝑧'/E(𝑧'; Θ) δ(θ, 𝑧'; Θ) ,

where Ω_m is the matter density fraction at z=0, and Ω_m∈Θ. We then discretise the continuous two-dimensional matter overdensity field along over θ using  pixels, while also discretising along the line-of-sight using Eq. (<ref>) and defining a mean redshift, z^(i), for a given weighted shell as follows

    𝑧^(i)(Θ) =  ∫d𝑧 𝑧 𝑊^(i)(𝑧; Θ)/∫d𝑧 𝑊^(i)(𝑧); Θ .


To accelerate the calculation of the discretised convergence fields for each shell, κ^(i), within , it makes use of the fact that for all Robertson-Walker space-times, transverse comoving distances will scale in such a way that one can write down a recurrence relation from one interval to the next along the line-of-sight <cit.>. This gives a recurrence relation for the convergence field, κ^(i), within shells where i ≥ 2, which only depends on κ^(i-1), κ^(i-2) and δ^(i-1) (see  for details). Such a relation allows one to calculate the convergence fields for all shells, while only holding three fields in memory at a given time. This greatly reduces the amount of computational resources needed per realisation, so the large amount of simulations at different cosmologies, Θ, needed for this analysis can be obtained with fewer resources when compared to forward models in previous analysis (see Sect. <ref>).

We now have all the ingredients to sample galaxies and to shear their shapes, while simultaneously having the flexibility to add any field-level systematics which may be relevant.



 §.§ Intrinsic alignments: non-linear alignment model

An important effect to consider when modelling weak gravitational lensing is intrinsic alignments (IAs). This refers to the fact that correlations between source galaxy shapes in different parts of the sky may not only be caused by weak gravitational lensing due to the foreground matter distribution. This happens due to two local processes which occur irrespective of weak gravitational lensing:  tidal alignments and tidal torquing. The latter is thought to arise when galaxy discs form perpendicular to the angular momentum axis of the tidal field, while tidal alignments occur when the gravitational tidal forces from neighbouring matter distributions cause the galaxies to align with the surrounding matter density field. This alignment leads to a coherent stretching or compression of the galaxy shapes as a function of the local matter overdensity field. Hence, the intrinsic ellipticities of source galaxies are systematically aligned with the underlying large-scale structure within which they form. Thus, the measured cosmic shear signal will be biased by these intrinsic alignments <cit.>.

To model this effect, we follow the methodology set out in <cit.>. In this prescription, for a given intrinsic alignment model, we define an associated effective convergence field, κ^(i)_IA. This κ^(i)_IA is a useful construct which describes the contribution to the observed weak lensing signal from IA under the Born approximation. Alternatively, an IA model could also directly be implemented into the dark matter distribution of the matter fields which could be an interesting avenue for future work.

Based on the fact that for the IA model, we take into consideration κ^(i)_IA∝δ^(i), we can assume that, for a given shell, this field can be added linearly to the underlying convergence field due to weak lensing, as given by

    κ^(i)(θ_m; Θ) →κ^(i)(θ_m; Θ) + κ_IA^(i)(θ_m; Θ) .

For simplicity and consistency with typical modelling assumptions in weak lensing surveys such as KiDS-1000 <cit.>, we choose to model κ^(i)_IA with the Non-Linear Alignment (NLA) model <cit.>, so we consider tidal alignments and do not consider tidal torquing. The NLA model assumes that the bias in the shear signal from IAs is linearly dependent on the projected local tidal field.  The "non-linear" part of the NLA model then simply refers to the fact that it has been found that the modelling of IAs is more accurate when modelling the underlying large-scale structure using a non-linear matter power spectrum, rather than a linear one (, similar to our approach described in Sect. <ref>). This means that it is also proportional to the local matter overdensity field, which is given by <cit.>,

    κ_IA^(i)(θ_m; Θ) = - A_IAC_1 Ω_m ρ_cr(z^(i); Θ)/D(z^(i); Θ) δ^(i)(θ_m; Θ) ,

where A_IA is the intrinsic alignments amplitude (which we are treating as a nuisance parameter that is also sampled, see Sect. <ref>), C_1 is a normalisation constant which we set to C_1 = 5 × 10^-14 h^-2 M_⊙^-1Mpc^3 in agreement with the IA measurements at low redshifts by SuperCOSMOS <cit.>, ρ_cr(z^(i); Θ) is the mean critical matter density as a function of redshift and D(z^(i); Θ) is the linear growth factor normalised to be unity at z = 0. Note that in Eq. (<ref>), we do not consider an explicit redshift dependence which is often expressed as a power-law term <cit.>. We omit this term, as it is not considered in the fiducial KiDS-1000 analysis <cit.>, and it has been found that at least for z < 1, the IA signals do not vary largely with redshift <cit.>. 

Nevertheless, the NLA model still comes with some caveats. Mainly, it neglects non-linearities which may not be captured in the underlying matter power spectrum, such as source density weighing, i.e. the next-to-leading order perturbations proportional to the source density considered in a full tidal alignment model <cit.>. This effect can in fact produce a signal in the shape measurements comparable to the correction due to non-linear structure growth <cit.>. Despite these limitations, the NLA model has been found to be accurate enough to model the IA signal seen by past and current weak lensing surveys <cit.>, including KiDS-1000 <cit.>.




 §.§ Galaxy positions and redshifts

With the cosmological dependence of the underlying large-scale structure characterised and astrophysical effects such as intrinsic alignments and baryonic feedback included, the forward simulations still need to include observational biases from our instruments. We choose to apply these on the level of galaxies, i.e. by biasing the sampled galaxy positions, redshifts and shapes/shears in accordance with our models for the observational systematic effects.

When sampling the galaxies' positions on the sky, we consider the survey footprint, Ω_survey, the galaxy number density which the survey can observe, n_gal, and its variability as a function of redshift and observational depth across the sky, i.e. as a function of the local observing conditions. When sampling the galaxies' redshifts, we take into consideration the redshift distributions calibrated from photometry for a given survey, n(z). Lastly, when determining the observed galaxy shapes, we consider the intrinsic ellipticity dispersion of the galaxies, σ_ϵ, the multiplicative and additive shear bias, while also including the effects on the measured galaxy shapes caused by the variation in the point spread function of the instrument as a function of the position in the sky.

Since we are interested in simulating a photometric survey, we must model the resolution along the line-of-sight realistically. This means that rather than modelling large-scale structure observables across N_shells concentric shells, they ought to be modelled across N_tomo tomographic bins. The tomographic bins of a given survey are defined by finite usually non-overlapping domains in photometric redshift, z_ph, in which all observed galaxies are grouped, such that z_p < z_ph≤ z_p+1, p ∈{ 1, 2, 3, ..., N_tomo}[Note that we use the indices p and q for tomographic bins and the indices i and j for underlying shells.]. The measured photometric redshift of a galaxy does not necessarily agree with its true/spectroscopic redshift, z. This means that galaxies which are grouped within a given photometric redshift bin, p, may not be located in this redshift range. To account for this, one can characterise a conditional redshift distribution for each tomographic bin P(z|z_ph) by mapping the measured z_ph to known spectroscopic redshift measurements which have a comparatively low bias in the redshift estimate <cit.>. Therefore, a given galaxy is in an associated pixel, m, shell, i, and tomographic bin, p. The probability of the galaxy being in tomographic bin, p, while it is sampled in shell, i, then depends on P(z|z_ph).

In line with the KiDS-1000 analysis, we define the effective galaxy number density as follows

    n_eff^(p)≡1/|Ω_survey|[∑_i ∈ (p) w_i]^2/∑_i ∈ (p) w_i^2 ,

where |Ω_survey| is the effective area covered by the survey footprint, M^(p) is the mean multiplicative shear bias in tomographic bin p as calibrated in <cit.>, and w_i is the  weight for a given galaxy i. Note that the effective galaxy density within the KiDS-1000 gold sample is not the same as the galaxy density, n_gal, when not correcting for the  weights for each galaxy. However, as KiDS-SBI does not model galaxy weights, in the forward simulations n_eff = n_gal since w_i = 1 ∀ i.

Similarly, we define the galaxy intrinsic ellipticity dispersion, σ_ϵ, within a given tomographic bin as follows

    σ_ϵ^(p)2≡∑_i ∈ (p) w_i^2  ( ϵ_obs, i, 1^2 + ϵ_obs, i, 2^2)/(1+M^(p))^2∑_i ∈ (p) w_i^2 ,

where ϵ_obs, 1 and iϵ_obs, 2 are the two components of the observed galaxy ellipticity. Applying Eqs. (<ref>) and (<ref>) to the KiDS-1000 gold sample, we obtain the values shown in Table <ref>.



Although it does not have detectable effects on the cosmic shear signal, the KiDS-SBI forward simulations sample galaxy positions such that they trace the underlying matter fields, so spatially varying systematics may be accurately applied. Using the [Speedy Acquisition for Lensing and Matter Observables; <https://github.com/Linc-tw/salmo>] framework <cit.> for simulating galaxy catalogues, we Poisson-sample galaxies within each  pixel, m, on the sky using the following expectation for the galaxy counts

    ⟨ N_m^(i)(p)⟩ (Θ) = w_m(Ω_survey)  [1+b^(i)δ^(i)(θ_m; Θ)]   P_m(p|i)   n_gal, m^(p)  A_pix, m ,

where w_m(Ω_survey) is the weight for a given pixel m which we take to be unity if a pixel is within the survey's mask/footprint, Ω_survey, and zero if it is not, b^(i) is galaxy bias parameter within a given shell which we assume to be b^(i) = 1 ∀ i for the purposes of this cosmic shear analysis, P_m(p|i) is the probability of a galaxy in pixel m being detected within tomographic bin p while being within shell i (i.e. a discretised version of P_m(z|z_ph)), n_gal, m^(p) is the observed galaxy number per unit area for a given tomographic bin p and a given pixel m (here the pixel dependence accounts for the spatial variability of the detection limit of galaxies as a function of observing conditions, the number of exposures, etc.; see Sect. <ref>), and A_pix, m is the area of pixel m. As n_gal, m^(p) changes due to variations in observational depth, it is more or less likely to observe high-redshift galaxies. Hence, a change in depth from one pixel to the next can change the observed redshift distribution (see Sect. <ref>).

Taking w_m(Ω_survey), P_m(p|i) and n_gal, m^(p) from the photometric survey's measurements, we can Poisson sample N_m^(i)(p)(Θ) given the expectation value shown in Eq. (<ref>) from which we can get the expected number of galaxy counts within each pixel of a tomographic bin p by taking

    N_m^(p)(Θ) = ∑_i=1^N_shells N_m^(i)(p)(Θ) .

With that, we then randomly sample each galaxy's right ascension and declination within each pixel m assuming a uniform distribution. This means that below the scales of pixels, the simulations do not have any information on galaxy clustering.

To sample each galaxy's redshift, we can use the fact that we know exactly how many galaxies are in each shell i for a given tomographic bin p and pixel m from N_m^(i)(p)(Θ). We can therefore randomly assign an index i and p to each sampled galaxy in accordance with that number. We can then randomly sample a specific z_spec for each galaxy by assuming that the galaxies are uniformly distributed within a given shell. Again, this means that weak lensing and galaxy clustering measurements from the simulated galaxies do not contain any information about the large-scale structure along the line-of-sight below the scales of shells.





 §.§ Galaxy shapes


In order to determine the galaxy shapes and shears, we begin by assigning intrinsic ellipticities. We define a given galaxy's shape through a complex ellipticity, ϵ, and any galaxy ellipticity in this work as follows

    ϵ = ϵ_1 + i ϵ_2 ,

where i≡√(-1). The intrinsic ϵ_1 and ϵ_2 of each galaxy are then sampled within  as two independent normal random variates with zero mean and a variance equal to σ_ϵ, m^(p)2. As with the galaxy density, n_gal, m^(p), the shape dispersion, σ_ϵ, m^(p), can also vary with tomographic bin and as a function of the position on the sky due to the effects of the spatial variability of observational depth. See Sect. <ref> for a more detailed discussion. The ellipticity of each source galaxy is then altered through weak gravitational lensing by the matter along the line-of-sight. The exact shape distortion is given by <cit.>,

    ϵ_lensed(Θ) = ϵ_int + g(Θ)/1+ g^*(Θ) ϵ_int ,

where ϵ_lensed(Θ) is the lensed galaxy ellipticity, ϵ_int is the intrinsic galaxy ellipticity and g(Θ) is the reduced shear (g ∈ℂ) which is given by

    g(Θ) = γ(Θ)/1-κ(Θ) ,

where γ(Θ) is the shear factor (γ∈ℂ). The use of Eq. (<ref>) is an extension to the KiDS-1000 analysis which assumed the reduced shear approximation, i.e. ϵ_lensed≈ϵ_int + g. In any case, this should cause a negligible change in the measured shear considering the sensitivity of KiDS-1000 <cit.>. To calculate the shear field, it is useful to decompose the spin-0 convergence field into spherical harmonics as follows

    κ_ℓ m^(i)(Θ) = ∫𝐝^2θ κ^(i)(θ; Θ)0*ℓ m(θ) ,

where κ_ℓ m^(i)(Θ) are the harmonic coefficients of the convergence field of a given shell i. One can then use the following relation to define the shear field's harmonic coefficients, γ_ℓ m^(i), given by <cit.>,

    γ_ℓ m^(i)(Θ) = - ((ℓ+2)(ℓ-1)/ℓ(ℓ+1))^1/2 κ_ℓ m^(i)(Θ) .

The γ_ℓ m^(i)(Θ) coefficients can then be used to calculate the spin-2 discrete shear field, γ^(i)(θ_m, Θ). Having defined the effective κ^(i)(θ_m, Θ) in Sect. <ref> and <ref>, while also defining γ^(i)(θ_m, Θ), we can assign a value of κ and γ = γ_1 +  i γ_2 to each galaxy sampled according to the procedure described in Sect. <ref>. This is done by taking the value of κ^(i)(θ_m, Θ) and γ^(i)(θ_m, Θ) within the shell i and pixel m in which a given galaxy is located. Between neighbouring pixels, the values of κ^(i)(θ_m, Θ) and γ^(i)(θ_m, Θ) are linearly interpolated. This choice may not be optimal <cit.>, but we find that at the precision level of KiDS-1000, it does not lead to any significant biases (see Appendix <ref>). We then use Eq. (<ref>) to calculate the associated reduced shear, g(Θ), as well as the lensed ellipticity, ϵ_lensed(Θ), by combining g with the intrinsic ellipticity using Eq. (<ref>).

However, this lensed ellipticity, ϵ_lensed(Θ), is not the observed shape measurement in a weak gravitational lensing survey such as KiDS-1000. Additionally, the shape measurement may be distorted by instrumental effects or shape modelling inaccuracies. Some relevant effects are selection biases, noise biases from a low signal-to-noise ratio, biases in the galaxy weights, artefacts, non-linear CCD responses, asymmetries in the point spread function, etc. (see e.g.  for a review). To account for this, it is common to parametrically estimate any residual systematics which may be affecting the shape measurements <cit.>. Within the forward simulations, we use the first-order parametric expansion of the observed galaxy shapes given by <cit.> as follows

    ϵ_obs, i ∈ m^(p)= (1+M^(p))ϵ_lensed, i + α^(p)ϵ_PSF, m + β^(p)δϵ_PSF, m + c^(p) ,

where i is the index for a single galaxy within tomographic bin p and pixel m, M^(p) is the multiplicative shear bias as measured for tomographic bin p, α^(p) gives the fraction of the PSF ellipticity which remains in the shear estimator, ϵ_PSF(m) is the local PSF measured within pixel m (see Figure <ref> for a map of ϵ_PSF in KiDS-1000), β^(p) gives the amplitude of the shear bias due to residuals which are not taken into account by the PSF model, δϵ_PSF represents the residuals in question, and c^(p) is the additive shear bias within a given tomographic bin p. Note that spatially varying multiplicative shear biases, M^(p), can also arise, but previous work has found that the mean bias across all positions is sufficient to accurately model the effect <cit.>. In accordance with previous findings <cit.>, we take β^(p) = 0  ∀ p. The amplitudes of the remaining systematics as calibrated by <cit.> are shown in Fig. <ref>. It becomes apparent that these effects are relatively small which explains why only the multiplicative and additive bias is taken into consideration in the modelling on the fiducial KiDS-1000 analysis <cit.>. Additionally, the measured values of M^(p), α^(p) and c^(p) have associated uncertainties. To take this uncertainty into account and make sure it is represented in the effective likelihood, rather than using a fixed value, we randomly sample different values for  M^(p), α^(p) and c^(p) for each realisation of the forward simulations from a normal distribution as shown in Fig. <ref>. Note that we do not explicitly marginalise over these parameters, but as they are sampled randomly from predetermined probability density distributions, so that their associated uncertainty is propagated into the effective likelihood. This is similar to how the seed for each log-normal random field, and the galaxy positions which populate the associated matter fields are varied from one forward simulation to the next, such that cosmic variance enters the effective likelihood as well.

As a result, for given cosmological parameters, Θ, the KiDS-SBI forward simulations produce a galaxy catalogue containing the galaxies' position on the sky, their spectroscopic redshift, the tomographic bin in which they are detected and their observed ellipticities, while taking into account many relevant systematics which are considered in the modelling of weak gravitational lensing measurements.









 §.§ Variable depth

As alluded to in Sect. <ref> and <ref>, when randomly sampling galaxies and their intrinsic ellipticities, the parameters which charaterise the probability distributions depend on the location on the sky, i.e. the pixel m in question. Specifically, the number of galaxies which is sampled (i.e. "observed") within a given pixel depends on the average galaxy density, n_gal, m^(p), measured within that pixel. The galaxy density varies from pixel to pixel as it depends on the observational depth of the survey in different parts of the sky, and on anisotropies and time variations in the atmospheric seeing (particularly, for surveys such as KiDS which only visit each pointing in the sky once at a given time). Since many parts of the sky are observed at different times, as is the case for ground-based telescopes such as the VLT Survey Telescope and the Visible and the Infrared Survey Telescope for Astronomy (VISTA) used for KiDS-1000, the atmospheric conditions and background light will be different for each pointing for a given exposure time. Thus, the depth of observations may change as seeing, the point spread functions and/or the background flux varies (for example, due to moonlight, zodiacal light or galactic absorption). 

In addition, surveys such as KiDS often have overlapping pointings in order to ensure that the footprint is observed without gaps. However, this implies that galactic sources which happen to be located near the edge of a pointing will be observed more often in different overlapping fields than a source located in the centre of a pointing. Consequently, near the edge of a pointing, there will tend to be a higher signal-to-noise ratio which allows for deeper observations. Since these systematic effects in the observational depth can happen at fixed scales and/or have certain periodicities, they can induce significant systematic effects into galaxy clustering and weak gravitational lensing signals. For cosmic shear in KiDS-1000, the bias can be near 1% in the signal and an average of ∼ 10 % in the standard deviation of cosmic shear observables <cit.>.

As the systematics modify the local observational depth, the galaxies which can be observed within a given pixel are selected according to their magnitude. This will systematically bias the dispersion of intrinsic ellipticities observed, σ_ϵ, m^(p), and further exacerbate the effect of observational depth variability on the cosmic shear signal.



To take all these effects into consideration within the forward simulations, we calibrate the spatial variability of n_eff, m^(p) and σ_ϵ, m^(p) from the KiDS-1000 measurements directly. Rather than defining maps over the entire survey footprint of both quantities for each tomographic bin, we use a map of a direct estimator of observational depth which also correlates with n_eff, m^(p) and σ_ϵ, m^(p): the root-mean square of the background noise, σ_rms. We found this quantity to be a good indicator of variable depth as it correlates well with the measured n_eff, m^(p) and σ_ϵ, m^(p), while being uncorrelated with the measured photometric redshifts and magnitudes in KiDS-1000.

A similar approach to modelling variable depth has been taken in <cit.>. However, in that analysis, the variable depth along the KiDS-1000 footprint was modelled with a direct estimate of the magnitude limit in the r-band (the band in which galaxy shape measurements are made in KiDS). We have found that the magnitude limit associated with a given galaxy in the KiDS-1000 gold sample is correlated with the galaxy's observed r-band magnitude as well as the photometric redshift estimate for the galaxy. This means that a selection according to the magnitude limit in the r-band can bias cosmological estimates, as the selection is not independent of the cosmic shear signal. In contrast, σ_rms is uncorrelated with the r-band magnitude measurements and the photometric redshifts in the KiDS-1000 gold sample, while still being highly correlated with the local magnitude limit in the r-band. Thus, we conclude it is a more direct parametrisation of the variations across the survey footprint of the galaxy selection.

As can be seen in Fig. <ref> and  <ref>, a map of σ_rms across the KiDS North and South fields clearly shows variations over different pointings as well as with areas of increased overlap. In addition, Fig. <ref> shows that n_eff^(p) and σ_ϵ^(p) vary linearly with σ_rms for KiDS-1000 data. Thanks to this fact, we do not need to model the variation of n_eff, m^(p) and σ_ϵ, m^(p) with the  pixel m with individual maps for each tomographic bin p. Instead, using the linear relations shown in Fig. <ref>, we may assign a value for n_eff, m^(p) and σ_ϵ, m^(p) for each σ_rms, m (i.e. the mean value of level within each pixel m contained within the footprint observed by KiDS-1000) given that

    n_eff, m^(p) =      a_n_eff^(p)σ_rms, m + b_n_eff^(p) ,
    σ_ϵ, m^(p) =     a_σ_ϵ^(p)σ_rms, m + b_σ_ϵ^(p) ,

where a_n_eff^(p) and a_σ_ϵ^(p) are the slopes of the linear fits shown in Fig. <ref> of galaxy density and galaxy shape dispersion for a given tomographic bin, respectively, and b_n_eff^(p) and b_σ_ϵ^(p) are the associated y-intercepts. The values of these parameters are shown in Table <ref>. This drastically helps the performance of the forward simulations, as this allows to just use a single map of σ_rms to model the variable depth rather than requiring 2N_tomo maps for it. Such a linear relation may not exist in other weak lensing surveys, but it should still be possible to numerically calibrate the relations dependence of n_eff^(p) and σ_ϵ^(p) on σ_rms to define an interpolation.

These inhomogeneities in the galaxy selection across the survey footprint also change the redshift distribution of galaxies. As σ_rms decreases, the local magnitude limit increases and fainter galaxies can be detected by the survey. As more distant galaxies with high redshifts also tend to be fainter, variable depth can shift the redshift distributions of the local galaxy population as can be seen in Fig. <ref>. The redshift distributions shown in Fig. <ref> have been calculated using the same approach as in KiDS-1000 where the photometric redshifts are calibrated from spectroscopic samples using self-organising maps <cit.>. The self-organising map has been applied to ten equi-populated subsamples of each tomographic bin of the KiDS-1000 gold sample which are binned in 100σ_rms along the following boundaries: { 1.7,   2.33,   2.57,   2.8,    3.04,   3.28,   3.52,   3.76,   4,   4.23,   12.96 }. To account for variable depth in our forward model, depending on the value of σ_rms in a given pixel, m, on the survey footprint, we sample galaxy redshifts in that pixel from the associated redshift distribution from Fig. <ref>. This is effect is important to model as an unaccounted shift in the redshift distribution can lead to significant biases in the cosmic shear signal <cit.>.



To summarise, in KiDS-SBI, we account for the effects of spatial variations in the observational depth on galaxy number density, galaxy shape dispersion and redshift distributions. Moreover, we also consider anisotropies in the shear bias by including distortions to the lensed galaxy ellipticities calibrated from the spatial variation in the point-spread function ellipticity in the KiDS-1000 gold sample (see Sect. <ref>).



§ MEASUREMENTS

Upon sampling galaxy shapes as well as shear, and applying all relevant instrumental systematics within KiDS-SBI, we measure the galaxy shapes, and the cosmic shear summary statistics thereof in the same manner as the measurements made from the KiDS-1000 gold sample's weak lensing catalogue. The only difference between the measurements of the forward-modelled catalogues and the KiDS-1000 gold sample is the values of the galaxy weights, w_i. In the forward simulations, the galaxy weights are set to unity for all sampled galaxies, w_i = 1 ∀ i, while the weights used in the measurements from the KiDS-1000 gold sample are given by the  weights (see details in Sect. <ref>). To avoid biasing the statistics of the galaxy populations in the forward simulations with respect to the real data, we sample the galaxy positions and shapes from the galaxy densities, n_gal, and intrinsic shape dispersions, σ_ϵ, where the corrections for the  weight and the multiplicative shear bias have been applied, as defined in Eqs. (<ref>) and (<ref>).



 §.§ Shape measurements

From the KiDS-1000 gold sample data as well as the forward simulated catalogues, we measure the observed galaxy ellipticity, ϵ_obs, for each galaxy per tomographic bin. To this end, we apply the same corrections to the shear measurements as described in <cit.>. The observed shape measurements are corrected as follows

    ϵ_obs, i^corr  (p)(Θ) = 1/1 + M^(p)∑_i ∈ pw_i (ϵ_obs, i(Θ) - ⟨ϵ_obs(Θ) ⟩ )/∑_i ∈ pw_i ,

where w_i is the galaxy  weight for galaxy i in the KiDS-1000 gold sample (for the case of the simulations, w_i = 1 ∀ i), ⟨ϵ_obs(Θ) ⟩ compensates for any additive biases <cit.>, i.e. the c^(p) term and the mean of α^(p) in Eq. (<ref>), and it is defined as the mean observed shape of all galaxies within tomographic bin p, and M^(p) is the mean multiplicative shear bias measured for each tomographic bin p as shown in Fig. <ref>. Note that M^(p) is not necessarily the same value as the value of M^(p) applied in Eq. (<ref>). For each instance of the simulations, a different M^(p) is drawn from a Gaussian probability distribution with mean M^(p) and the standard deviations shown in Fig. <ref>. Regardless of the value of M^(p) drawn, the shear measurements of each simulation are corrected by the same mean M^(p). Any discrepancy between these values is going to introduce noise in the simulations which accounts for the uncertainty on the shear bias measurements.




 §.§ Observed angular power spectra


With a full simulated cosmic shear catalogue akin to the KiDS-1000 gold sample, we can map any set of cosmological and astrophysical parameters, Θ, to a corresponding set of cosmic shear measurements from the KiDS survey as required for simulation-based inference. Although theoretically possible, it is computationally impractical to conduct cosmological inference at the level of catalogues. For this reason, it is useful to compress the catalogues down to useful statistics which still retain most of the relevant cosmological information about the underlying large-scale structure.

We choose to compress the catalogues down to two-point statistics. The main reason for this is that the forward simulation pipeline is designed to be only percent-level accurate in two-point statistics. Higher-order correlations of the shear field are captured to a certain degree within the simulations, but their accuracy is limited <cit.>. In addition, the use of two-point statistics assures that the only distinguishing feature of our analysis when compared to the fiducial cosmic shear KiDS-1000 analyses <cit.> is the dropping of the assumption of a Gaussian likelihood, so we can carry out a direct comparison between the analyses. This way we can ensure that any non-Gaussianities in our likelihood are either attributed to inherent non-Gaussianities in the likelihood of two-point statistics or non-Gaussianities induced by systematic effects.

Due to its computational efficiency when compared to spatial two-point correlation functions, we choose pseudo-Cls as our two-point statistic. We decompose the observed shear field, ϵ_obs^corr  (p)(θ_m; Θ), into its curl-free E-modes and its divergence-free B-modes as follows

    ϵ_obs^corr  (p)(θ; Θ, Ω_survey) = ∑_ℓ=0^ℓ_max∑_m=-ℓ^ℓ (Ẽ_ℓ m^(p) (Θ) +  iB̃_ℓ m^(p)(Θ)) 2ℓ m(θ)  ,

where ± 2ℓ m(θ) are the spin-2 spherical harmonic functions which define an orthonormal basis for the observed shear field, such that 

    ∫d^2θ± 2ℓ m(θ) ± 2*ℓ' m'(θ) = δ^K_ℓℓ'δ^K_m m' ,

while Ẽ_ℓ m^(p) (Θ) and B̃_ℓ m^(p) (Θ) are spherical harmonic coefficients of the curl-free and the divergence-free observed shear fields within each tomographic bin, respectively, and δ^K is the Kronecker delta. These coefficients are defined as follows:

    Ẽ_ℓ m^(p) (Θ) =  1/2 ∫d^2    θ  [ϵ_obs^corr  (𝑝)(θ; Θ)2*ℓ m(θ)
        + ϵ_obs^corr  (𝑝)*(θ; Θ)-2*ℓ m(θ)] ,
    B̃_ℓ m^(p) (Θ) = -i/2∫d^2θ    [ϵ_obs^corr  (𝑝)(θ; Θ)2*ℓ m(θ) 
        - ϵ_obs^corr  (𝑝)*(θ; Θ)-2*ℓ m(θ)] .


From these coefficients, we then calculate the pseudo-Cls, C̃^(pq)_ϵϵ(ℓ; Θ), as follows


    C̃^(pq)_ϵϵ, μ(ℓ; Θ) = [ C̃_ϵϵ^EE (pq)(ℓ; Θ); C̃_ϵϵ^EB (pq)(ℓ; Θ); C̃_ϵϵ^BB (pq)(ℓ; Θ) ] = 1/2ℓ + 1∑_m = -ℓ^ℓ[   Ẽ^(p)_ℓ mẼ^(q)*_ℓ m  (Θ);  Ẽ^(p)_ℓ mB̃^(q)*_ℓ m  (Θ); B̃^(p)_ℓ mB̃^(q)*_ℓ m  (Θ) ] ,

where μ∈{1, 2, 3} such that 1 stands for the EE component, 2 for the EB component and 3 for the BB component. Note that C̃_ϵϵ^EB(ℓ; Θ) = C̃_ϵϵ^BE(ℓ; Θ), since Ẽ_ℓ mB̃_ℓ m^*  (Θ)=  B̃_ℓ mẼ_ℓ m^*  (Θ). 



Going forward, similar to other previous analyses <cit.>, we only take into consideration C̃^EE(pq)_ϵϵ(ℓ; Θ), so any cosmological signal which may have ended up in  C̃^EB(pq)_ϵϵ(ℓ; Θ) or  C̃^BB(pq)_ϵϵ(ℓ; Θ) because of EE to EB, or EE to BB mode mixing will be lost. C̃^EE(pq)_ϵϵ(ℓ; Θ) is expected to still contain all of the cosmological cosmic shear signal because the shear field resulting from a scalar gravitational field in ΛCDM is predicted to be a curl-free field <cit.>. 

To account for the EE signal from the correlations between the intrinsic galaxy shapes, i.e. the shape noise bias, we subtract the mean shape noise power spectrum as follows

    C̃^(pq)_ϵϵ(ℓ; Θ) = C̃^EE(pq)_ϵϵ(ℓ; Θ) - δ^K_pq⟨C̃^EE(pq)_noise(ℓ) ⟩ ,

where C̃^(pq)_ϵϵ(ℓ; Θ) is the E-mode pseudo angular power spectrum for cosmic shear with the shape noise bias subtracted, while ⟨C̃^EE(pq)_noise(ℓ) ⟩ is the mean of the curl-free angular power spectrum of the shape noise. The latter is estimated as follows <cit.>,

    C̃^EE(pp)_noise(ℓ) = 1/2ℓ + 1∑_m = -ℓ^ℓẼ^rand  (p)_ℓ m (Θ) Ẽ^rand  (p)*_ℓ m  (Θ) ,

where Ẽ^rand  (p)_ℓ m is the curl-free spherical harmonic coefficient of the randomly rotated shear values, ϵ^rand  (p)_obs(θ_m; Θ). In turn, we define this field from the following galaxy shear values as follows:

    ϵ^rand_obs, i(Θ) =  ϵ^rand _obs, i, 1(Θ) +  iϵ^rand _obs, i, 2(Θ) ,
       ϵ^rand _obs, i, 1 =  ϵ^corr _obs, i, 1 cos(θ_rand, i) - ϵ^corr _obs, i, 2 sin(θ_rand, i) ,
       ϵ^rand _obs, i, 2 =  ϵ^corr _obs, i, 2 cos(θ_rand, i) + ϵ^corr _obs, i, 1 sin(θ_rand, i) ,

where θ_rand, i is a randomly drawn angle for each galaxy i from a uniform distribution where θ_rand, i∈ [0, 2π). To calculate the mean of the curl-free angular power spectrum of the shape noise, ⟨C̃^EE(pq)_noise(ℓ) ⟩, we take mean of all modes in C̃^EE(pq)_noise(ℓ) as follows

    ⟨C̃^EE(pq)_noise(ℓ) ⟩ =  ∑_ℓ = ℓ_min^ℓ_max(2 ℓ + 1)  C̃^EE(pq)_noise(ℓ)/∑_ℓ = ℓ_min^ℓ_max (2 ℓ + 1) ,

where the mean is weighted by a factor of 2 ℓ + 1 to account for the fact that each angular scale, ℓ, contains 2 ℓ + 1 modes, m. We find setting ℓ_min∼ 200 yields optimal results, as high angular scales are typically shape noise-dominated, and thus put better constraints on the average shape noise bias. In fact, we find that the measured average shape noise bias, ⟨C̃^EE(pq)_noise(ℓ) ⟩, agrees within <0.5% with the underlying shape noise signal computed directly from the known intrinsic galaxy shapes. Additionally, we find that ⟨C̃^EE(pq)_noise(ℓ) ⟩ is self-consistent with ⟨C̃^BB(pq)_noise(ℓ) ⟩ within <0.01%. The introduction of a weighted mean given by Eq. (<ref>) allows us to achieve this level of precision with only a single random rotation of the galaxies' shapes. Instead, previous approaches relied on creating many different instances of random rotations and then computed ⟨C̃^EE(pq)_noise(ℓ) ⟩ by taking the average over all the instances while weighting each angular scale, ℓ, equally (see e.g. ). When weighting each ℓ equally, the low-ℓ scales have an equally strong pull on the estimator as the high-ℓ scales which increases the statistical uncertainty on ⟨C̃^EE(pq)_noise(ℓ) ⟩ as the large scales are mostly dominated by the large-scale structure signal. For this reason, it is common to take the mean over many realisations of the shape noise to reduce the variance of the mean. Although such approaches achieve similar levels of precision when estimating the shape noise bias as the approach shown in this work, computing the angular power spectra for hundreds of realisations of randomly rotated galaxy shape catalogues can become computationally expensive and time-consuming. Hence, it is not feasible for an SBI analysis which requires ∼10^4 realisations and we opt to use the weighted mean of a single realisation of a random rotation to estimate the shape noise bias.

Furthermore, we bin the noise-free observed angular power spectrum, C̃^(pq)_ϵϵ(ℓ; Θ), into 8 log-spaced bins between ℓ = 76 and ℓ = 1500 in line with <cit.>. For this, we choose the pseudo-Cl binning scheme described in <cit.> given by

    C̃^(pq)_ϵϵ, L(Θ) = 1/2 π∑_ℓ = ℓ_L^ℓ_L+1ℓ(ℓ+1)/(ℓ_L+1 - ℓ_L)C̃^(pq)_ϵϵ(ℓ; Θ),

where ℓ_L and ℓ_L+1 are the lower and upper limits of the L^th bin, respectively.

Note that we do not choose to deconvolve the pseudo-Cls to estimate the full-sky E-mode angular power spectra for cosmic shear. We find that with the sky-coverage of KiDS-1000 and the complexity of the geometry of the KiDS-1000 footprint, the mixing matrix (see Appendix <ref>) is not necessarily invertible, and the deconvolution is not single-valued.

We therefore obtain for a given forward simulation a cosmic shear measurement like the one shown in Fig. <ref> for a single run of the forward simulations. The calculation of the data vector calculated from the theory prediction shown in Fig. <ref> is described in Appendix <ref>. As can be seen from Figure <ref>, the measured cosmic shear signal is consistent with theory predictions which we find is the case at all cosmologies throughout the prior volume (see Appendix <ref>). This shows that the log-normal random fields accurately recover the two-point statistics of the galaxy populations throughout parameter space, as expected.

As seen in Fig. <ref>, for the five tomographic bins in KiDS-1000, we obtain 8 × N_tomo(N_tomo + 1)/2 data points in C̃^(pq)_ϵϵ, L(Θ), i.e. a 120-dimensional data vector. The analysis choices described within Sect <ref> and <ref> in conjunction constitute the forward model assumed within this SBI analysis, and it is henceforth labelled as the anisotropic systematics model, where the plus sign is in reference to the fact that the model includes novel features such as variable depth, PSF variation, and dropping of the reduced shear approximation.

To summarise, each data vector depends on 12 parameters in Θ: 5 cosmological parameters (σ_8, the root-mean-square matter fluctuation over 8 Mpc/h_0; ω_b, baryonic matter density; ω_c, cold dark matter density; n_s, the scalar spectral index of the primordial density fluctuation power spectrum, P_ℛ; and h_0[Note that throughout we deviate from standard notation and denote the dimensionless Hubble constant as  h_0.], the normalised Hubble constant), and 7 astrophysical/nuisance parameters related to systematics: A_bary, the baryonic feedback amplitude within the non-linear three-dimensional matter power spectrum; A_IA, the galaxy intrinsic alignment amplitude of the NLA model; and five correlated δ_z parameters which define the shift in the mean of the source redshift distribution, P(z|z_ph), of each tomographic bin. All other parameters on which the simulation depends are fixed from run to run (e.g. ω_k = 0, the equation of state of dark energy, w = - 1, ∑ m_ν = 0.06 eV, etc.). The exceptions to this are the seed used to sample the matter overdensities and the galaxies, and the amplitudes of the shear biases as described in Sect. <ref>. These values are varied from run to run in order to simulate cosmic variance, shape noise as well as the uncertainty in the shear bias, respectively.




§ SIMULATION-BASED INFERENCE (SBI)


The final data vector as defined by the forward simulations is used to determine the effective likelihood of the data using sequential neural likelihood estimation based on the Density Estimation Likelihood-Free Inference package ([<https://github.com/justinalsing/pydelfi>], , or  henceforth to denote the method and software package) using the same analysis pipeline outlined in .

To implement this, we create a bespoke interface that allows  to be the sampler for  <cit.>. Initially, the sampler chooses a fixed number of points within the considered parameter space based on the hypercube-generating algorithm from PyDOE[<https://github.com/tisimst/pyDOE>] with some minor modifications (see  for more details). Subsequently, as is outlined in Fig. <ref>, the measured pseudo-Cls for each evaluation of the parameters on the hypercube are compressed further using score compression (see Sect. <ref>). This reduces the dimensionality of the data vector to the size of the parameter vector (in this case, 7 dimensions). The compressed data is then used to train neural density estimators through  (see Sect. <ref>). 

As the ensemble of neural density estimators evaluates the sampling distribution throughout the hypercube, through the use of active learning,  finds areas in parameter space where the likelihood is undersampled to select new parameter vectors (see  for details). As found in , this can halve the number of simulations needed to accurately train  by drawing more simulations in the parameter region of interest when compared to a Latin hypercube alone. Once the ensemble of neural density estimators has reached the stopping criterion, controlled for over-fitting by comparison to a hold-out validation data set, we can calculate an effective likelihood from the learnt sampling distribution condition on the observed data vector (which can be a mock or measured from real data). We can then obtain a posterior distribution by sampling the product of the likelihood and the prior. This approach of learning the sampling distribution conditioned on the observed compressed data to yield an effective likelihood has the advantage that if the prior or data vector were to be changed, there is nothing new that must be re-trained within .

We refer to the model choices that we make within KiDS-SBI as the anisotropic systematics model. The model follows the choices outlined in Sects. <ref> and <ref>. It is designed to be consistent with previous analyses of KiDS-1000, specifically, with <cit.>, while adding additional realism. Both <cit.> and this work assume a flat ΛCDM cosmology to model the cosmic shear signal, both use pseudo-Cls as their data vector of choice, and both consider systematics such as multiplicative shear bias, variable depth in the uncertainty and intrinsic alignments in the signal. The anisotropic systematics model differs from <cit.> in that it also includes variable depth in the signal modelling (as described in Sect. <ref>), it considers the effect of intrinsic alignments on the likelihood, and it takes into account the variance in the additive and PSF shear biases (as described in Sect. <ref>). This is a consequence of the fact that in KiDS-SBI the effects which are modelled in the signal are intrinsically considered in the uncertainty model. In <cit.>, these effects are only considered separately as they pertain to the signal or the uncertainty, not both at once.





 §.§ Pipeline performance


A single such evaluation of the forward simulations from cosmological parameters to the shear pseudo-Cls as described in Sect. <ref> runs within ∼20 minutes on a single core (with N_shells = 19, N_tomo = 5 and N_side = 1024). We choose this spatial resolution N_side, because it implies that the measured two-point statistics should be representative of the input angular power spectra up to an ℓ_max∼ 2 × N_side <cit.>. In addition, for KiDS-1000's galaxy number density, the pixel size at an N_side = 1024 is sufficient for almost all pixels to contain at least one galaxy. This suppresses mode mixing in the pseudo-Cls due to random masking of empty pixels which were observed, as at this resolution all pixels are on average populated by at least one galaxy (see Appendix <ref>). 

In comparison, as can be seen in the bar chart in Fig. <ref>, log-normal random field simulations similar to the ones presented in <cit.>, take ∼280 minutes to compute a single forward simulation when run with the same accuracy and precision settings on a single core. In practice, the simulations described in <cit.> are used exclusively for the characterisation of a numerical covariance matrix. Hence, only a single evaluation of the underlying power spectrum from  is necessary for every simulation. At the same time, the <cit.> simulations have a lower resolution of shells along the line-of-sight, so the runtime of  <cit.> in the original setting is lower than the runtime indicated in Fig. <ref>. We also note that the runtime of  <cit.> could be drastically reduced from the single-core runtime shown in Fig. <ref> by increasing the number of cores and/or by pixelising the galaxy shears prior to calculating the spatial two-point correlation functions, ξ(θ).

These gains in speed are driven by four main factors: the fast non-Limber integration using the Levin method, the use of recurrence relations to calculate the convergence within , the choice of using a quick-to-compute summary statistic like angular power spectra, C̃(ℓ), and the efficient subtraction of the shape noise bias from a single random rotation of the galaxy shapes. This is a substantial improvement with respect to other typical codes used for such simulations, while also improving accuracy with respect to previous models by including systematics, improving the resolution along the line-of-sight, and dropping common approximations like the Limber and reduced shear approximation. Thanks to this, it becomes feasible to increase the precision of our simulations as necessary and add the realism discussed in Sect. <ref>, <ref> and <ref>, while still being able to compute the ∼10^4 forward simulations needed to adequately characterise the effective likelihood .





 §.§ Parameters and priors


The parameters which are varied in the simulation-based inference in this work are shown in Table <ref>. Their priors match those of the previous KiDS-1000 cosmic shear analyses <cit.>.

To avoid overinformative priors, most are flat top-hat functions. For the same reason, the top-hat priors in S_8, h_0, ω_b, A_IA and A_bary are chosen to be wide. The prior on n_s spans a smaller range around the theoretical value of unity for scale-invariant primordial fluctuations. This avoids artefacts within the prior volume as this parameter is not well constrained by weak gravitational lensing. The prior on ω_c is defined to be consistent with a range in Ω_m∈ [0.188, 0.408], where the limits are given by the ± 5 σ intervals of the marginal constraints from independent measurements of luminosity distance to Type Ia Supernovae in <cit.>. See <cit.> for more details on the motivations for the chosen priors. The priors on the nuisance parameters capturing any shifts in the mean of each of the five tomographic bins, δ_z, are given by a multivariate Gaussian, since the δ_z of a given tomographic bin is not independent of the shifts in the other bins. This is quantified by the covariance, C_z, as estimated in <cit.>.





 §.§ Score compression

To improve the computational efficiency of the analysis and to facilitate the use of , we compress each of the 120-dimensional data vectors measured from each forward simulation using score compression as described in <cit.>.

If the likelihood is known a priori, one can compress a given data vector to a summary of the same dimensionality as the degrees of freedom in the assumed model, such that Fisher information is conserved <cit.>. This allows us to compress the data vector, d, down to a vector, t, with the dimensionality of |Θ| as follows

    t = ∇ℒ(d | Θ_*) ,

where ℒ is the log-likelihood distribution evaluated at the fiducial cosmology, Θ_*. 

In this work, the exact form of the likelihood is not known a priori, as the main motivation is to characterise the form of the effective likelihood. However, it is still possible to perform score compression on the data by assuming an analytical form for the likelihood. To this end, we assume a Gaussian likelihood with a fixed covariance matrix from 1,600 forward simulations at the fiducial cosmology given in Table <ref> and evaluate its gradient through a five-point stencil near the fiducial cosmology. In this form with a Gaussian likelihood, this compression is equivalent to  <cit.> or a linear compression based on Karhunen-Loéve eigenvalue decomposition <cit.>.

This compression is optimal if the true likelihood is Gaussian, which it is asymptotically near the peak of the likelihood and the chosen fiducial set of parameters is equal to the true parameter values. The downside of this is that the compression can lose information if this is not the case. Nevertheless, it was found in  that, in KiDS-SBI, the score compression is robust to suboptimal choices of fiducial parameter values and data covariance.

Additionally, we find that the 12 parameters shown in Table <ref> are not necessary to characterise the degrees of freedom of the 120-dimensional cosmic shear pseudo-Cls. As all δ_z are broadly consistent with zero, we may compress all data vectors to a seven-dimensional summary which is still capable of capturing all the complexity in the data. In any case, the δ_z parameters are still explicitly varied within the simulations to propagate any uncertainties on the mean of the tomographic bins. For simplicity, henceforth all posteriors are implicitly marginalised over the five δ_z parameters.



 §.§ Density estimation likelihood-free inference ()
 
To estimate the effective likelihood for the anisotropic systematics model, we train neural density estimators (NDEs) using  <cit.> such that the NDEs learn the sampling distribution between the compressed data and the input cosmological parameters. To this end, we initialise an ensemble of six independent conditional Masked Autoregressive Flows (MAFs;  ). Each MAF is made up of between three and eight Masked Autoencoders for Density Estimation (MADEs; ) each with two hidden layers of 50 neurons. For a detailed definition and description of these network architectures, see <cit.>.

As the first step in the exploration of the prior space, the ensemble of these neural density estimators (NDEs) initially learns the sampling distribution from a set 2000 points on a Latin hypercube. Based on this, we then make use of the active learning feature within  to sample additional points in parameter space which efficiently contribute to learning the effective likelihood within prior space. To ensure that the sampling distribution density, P(t| Θ, w) is well learnt, recursive sampling continues until the validation loss stops decreasing for 20 training epochs, where the loss function is defined as,


    - ln [U(w)] = - ∑_i=1^N_simsln[P(t_i|Θ_i, w)],

where N_sims is the number of forward simulations and the loss function, - ln [U(w)], is defined such that it is a Monte Carlo estimate of the Kullback-Leibler divergence between the learnt effective likelihood and the true effective likelihood (see  for a more detailed explanation). This process is repeated over an ensemble of many independent neural density estimator networks, where each network learns the conditional likelihood independently. The final estimate of the effective likelihood is then given by the weighted sum of the learnt distributions from each neural density estimator.

We find that for the anisotropic systematics model, this occurs already after 5,000 realisations. We note that the NDEs require fewer realistic KiDS-SBI forward simulations to minimise the loss function as is the case in  where the forward simulations were idealised random samples from a Gaussian distribution. In any case, to ensure that the prior space is densely sampled, we choose to train the NDEs on 18,000 realisations. This is facilitated by the computational efficiency of the forward simulations making additional realisations relatively inexpensive.

Upon training each of the six MAF NDEs using forward simulations independently, the NDEs are stacked while weighting each according to their relative validation losses. From the stacked NDE, we obtain the final posterior distributions by sampling with  <cit.> and check for consistency with  <cit.>.

As advised in , to minimise the loss of information through from the compression step described in Sect. <ref>, once the simulation-based inference has been conducted the first time on the KiDS-1000 data for a given forward-model, we repeat the inference again a second time while taking the fiducial model parameters assumed in the score compression to be given by the maximum a posteriori (MAP) from the initial inference (see Tab. <ref>). Note that in the second inference, we do not change the forward model in any way, and keep the exact same realisations of the model as in the first analysis. Only the fiducial data vector used to characterise the likelihood in Eq. (<ref>) are updated to re-compress all simulated data vectors. Considering that the initial fiducial parameters based on previous KiDS-1000 analyses are mostly consistent with the MAP which we infer in the initial analysis, while score compression is also robust against the biased fiducial parameters <cit.>, we find that this step only marginally alters the inferred posterior.



§ GOODNESS-OF-FIT IN SBI
 
As the inference is performed using SBI, we cannot perform a goodness of fit test the traditional way when we have access to an analytical likelihood. However, we still wish to keep this metric and so create a novel implementation of the χ^2 goodness-of-fit test. We follow the work of <cit.> and implement a SBI goodness-of-fit test from simulations in the [<https://github.com/Kiyam/ecp_gf_tests>] module.

The methodology is as follows, initially, we find the maximum a posteriori (MAP) values from the posterior sampled in the SBI using the Nelder-Mead <cit.> optimiser. As the parameters of n_s and a_bary are unconstrained by our analysis, and the MAP can therefore vary freely across the prior, instead of making use of the MAP cosmology, we run simulations at the cosmology chosen to perform score compression for those parameters, which is already a pseudo maximum likelihood estimation (see <cit.> for more details).

At the MAP cosmological parameter values, we simulate the sampling distribution by running a further 500 forward simulations while randomly varying noise realisations. We can thus define a χ^2 discrepancy measure as follows

    χ^2(t|Θ) (t_i - E[t_*|Θ_*])^T(Cov(t_*|Θ_*))^-1(t_i - E[t_*|Θ_*]) ,

where t_* is the compressed data vector at the fiducial cosmology, Θ_* the fiducial cosmology and t_i one of n randomly varied noise realisations of the compressed data vector at the MAP parameter values. A fully Bayesian implementation would simply require forward simulations to be realised from posterior parameter draws as opposed to only at the inferred MAP parameter value. We have stuck to making use of simulations at the MAP to be analogous to the traditional KiDS-1000 analysis. Furthermore, as we are making use of linear score compression, we have linearised our model, and thus our goodness of fit is analogous to the classic χ^2 test. The main difference is that the degrees of freedom of the maximum χ^2 does not trivially reflect the degrees of freedom of our model due to the effects of optimal score compression in performing data dimensionality reduction.



For our setup, we evaluate a mean model data vector and a covariance from a set of 1,600 simulations run at the fiducial cosmology. It should be noted that these simulations are the same as the ones required for the score compression step described in Sect. <ref>, so this comes without any additional computational costs. Using this newly obtained model data vector and covariance, and the χ^2 discrepancy measure defined above, we calculate the discrepancy measure for the 100 separate forward simulations from the anisotropic systematics model, and the KiDS-1000 gold sample data, which are depicted in the left panel of Fig. <ref>. As we obtain a probability to exceed (PTE) of 0.42, we find the fit for the data to be good.

Note that, in this test, the degrees of freedom are not recovered directly due to the effects of normalisation and optimal score compression in doing dimensionality reduction.




§ VALIDATION
 
The SBI pipeline presented in this work has previously been extensively tested in an idealised case. From the findings in , we conclude that the SBI pipeline based on  implemented within KiDS-SBI robustly and accurately recovers the posterior distribution of ΛCDM cosmological parameters from KiDS-1000 cosmic shear data, even if the data compression loses information. However,  only considered simulated data whose likelihood is by construction Gaussian. For this reason, it becomes necessary to test the robustness of the SBI pipeline again for the anisotropic systematics forward model. We find in Appendix <ref> that accurately and robustly recovers the underlying posterior distribution of ΛCDM parameters.



 §.§ Coverage






To determine whether the learnt posterior distribution from the anisotropic systematics model is representative of the true underlying posterior and is unbiased, we conduct the Tests of Accuracy with Random Points (TARP) shown in Fig. <ref>. TARP <cit.> measures the expected coverage probability of random posterior samples within a given credibility level of the learnt posterior empirically. Whilst the original authors have their implementation of the coverage test, we made use of our implementation (note:gf) to interface with the outputs from our simulations and DELFI in a more straightforward manner.  

We did not perform the coverage test over our entire parameter space, as there are many cosmological parameters that are entirely unconstrained by our inference. Instead, the coverage test was only performed on S_8 and Ω_m as those are the parameters of greatest interest for cosmic shear.

To perform this coverage test, we first ran 500 new simulations at the fiducial cosmology for A_IA, n_s, h_0, A_bary, freely varying the correlated δ_z [5] values, but varied S_8 ∈ [0.5, 0.9], ω_b∈ [0.019, 0.026], ω_c∈ [0.051, 0.255] and h_0 ∈ [0.64, 0.82], with a latin hypercube. A slightly narrower range was chosen for S_8 than the prior for the full inference so that the testing points were positioned closer to the final parameter posterior range of interest. For each of these simulations, the noise realisation was allowed to freely vary. We also performed a coverage test with A_IA also varied, and obtained similar results to what is shown in Fig. <ref>.

We then used our trained  model to perform inference on the new suite of 500 simulations, each acting as a mock data vector and thus obtaining 1000 posterior points for all 500 simulations. From the posterior samples, the coverage probability was calculated empirically for a given credibility level (see  for a detailed discussion of the algorithm). Due to the inherently stochastic nature of the coverage probability test, the test was bootstrapped 100 times for each credibility level. Fig. <ref> depicts both the mean and two standard deviations of the spread of the bootstrapped TARP coverage test.

For a well-estimated posterior that is not biased and representative of the true underlying posterior, we expect the coverage probability to be directly proportional to the credibility level. Any biases in the learnt posterior would lead to deviations from linearity in the TARP.

As can be seen from Fig. <ref>, the learnt posterior is very close to being ideal. As the mean coverage is seen to be slightly lower than the credibility level but always within one standard deviation of the ideal case. This means that the learnt posterior is slightly underconfident. This is acceptable as it implies that the KiDS-SBI posterior estimates are conservative but not biased. We repeat this test while also varying A_IA, and find similar results.



 §.§ Impact of variable depth and shear bias

To assess the impact of the additional observational systematics considered in the anisotropic systematics model, we re-train a set of NDE ensembles on another 18,000 simulations of a different model which does not consider variable depth or additional shear biases while following the same procedure as outlined in Sect. <ref>. We refer to this model as the standard isotropic systematics model as it is more in line with previous "standard" KiDS-1000 analyses such as <cit.>, <cit.>, <cit.>, or <cit.>. This model considers the cosmological signal, intrinsic alignments, and baryonic feedback in the same way as the anisotropic systematics model, but it has two main differences: it does not consider variable depth, and it does not consider the additive or PSF shear biases. Firstly, when sampling galaxies as described in Eq. (<ref>) the galaxy density, n_gal, and the photometric redshift distribution, P(p|i), are assumed to be isotropic for a given tomographic bin. Secondly, when biasing the observed galaxy shapes according to the shear biases measured in KiDS-1000, the standard isotropic systematics model only considers the multiplicative shear bias, i.e. only M^(p)≠ 0 while α^(p) = 0,   c^(p) = 0  ∀ p.

This simplification of the forward model impacts the measured cosmic shear signal in a few ways. Firstly, from just comparing the measured two-point statistics, as can be seen in more detail in Appendix <ref>, we find that variable depth can contribute up to 1% to the measured two-point shear signal in KiDS-1000, which is consistent with the predictions from <cit.>. At the same time, variable depth can significantly alter the uncertainty of the two-point statistic (by ≲ 20%; ). As the intrinsic galaxy shape dispersion and the galaxy density vary anisotropically, the shape noise of the sample is no longer determined by the mean value per tomographic bin. Instead, the shape noise becomes an anisotropic distribution too which in the case of KiDS-1000 is skewed below the mean (see Appendix <ref>). In addition, as can be seen by eye from Fig. <ref> and <ref>, variable depth can add angular correlations to the data at the scales of pixels, pointing overlaps, whole pointings or even over the whole footprint, while also changing the signal along the line-of-sight anisotropically as shown in Fig. <ref>. These additional correlations can significantly alter the covariance of the two-point statistic or may even lead to non-Gaussian noise.

In fact, these non-cosmological correlations are significant enough for us to find that the overall uncertainty on the two-point statistic of a model including variable depth is larger than the uncertainty in a model which does not consider the effect. Hence, variable depth contributes approximately 1% (∼ 0.2 σ) additional non-cosmological two-point signal.



The other distinguishing systematic between the anisotropic systematics model and the standard isotropic systematics model is the shear bias. When including additive and PSF shear biases in the anisotropic systematics model, that can also add a per-cent level (up to ∼ 0.3 σ, but less than 0.1 σ at most angular scales) non-cosmological signal to the pseudo-Cls, mostly driven by the PSF variation shear bias (see Appendix <ref> for details). This is consistent with the negligible impact on the measurement found in previous studies <cit.>. At the same time, the inclusion of the α_1^(p), α_2^(p),  c_1^(p) and c_2^(p) parameters for each tomographic bin and their associated uncertainties adds additional variance to the anisotropic systematics model when we pre-marginalise over these 20 parameters. Therefore, the additive and PSF variation shear biases are yet another source of additional non-cosmological signal and variance.

These effects clearly propagate to the posteriors as can be seen when comparing the full simulation-based inference based on a likelihood learnt from 18,000 realisations of the standard isotropic systematics model to a likelihood based on 18,000 realisations of the anisotropic systematics model in Fig. <ref>.

To assess whether the inclusion or exclusion of variable depth and shear biases can bias the SBI of cosmological parameters, we repeat the inference of a mock data vector, but this time we generate a random realisation from the standard isotropic systematics model and analyse it using the likelihood learnt from the forward simulation which assumes the anisotropic systematics model, and vice versa. This allows us to determine whether the inclusion of the effects in the forward model would bias our inference if they were not actually present in the data that is analysed. We also show these posteriors in Fig. <ref> where we see that when the data model and the likelihood model are mismatched, the constraints are consistent within a given model. This confirms that, given the learned likelihood, the bias in the two-point shear signal from both variable depth and PSF shear bias is negligible in KiDS-1000. These findings are consistent with previous work which found this to be the case for KiDS-1000 for variable depth <cit.> and for the PSF shear bias <cit.>.

Nonetheless, systematics such as variable depth and the PSF shear bias do not only affect the measured signal, but also its uncertainty. As shown in the right panel of Fig. <ref>, when analysing the same noisy mock data vector simulated with the anisotropic systematics model with the SBI trained on both models (anisotropic systematics and standard isotropic systematics), the marginal in S_8 shifts considerably by about 0.7σ. At the same time, the constraints based on the anisotropic systematics model give ∼5% smaller 1σ confidence intervals in the S_8 marginals than standard isotropic systematics model. As each model characterises a different effective likelihood, when analysing the same data vector with both models, each data point is weighted differently depending on the model assumed during the inference. Depending on the noise realisation of the mock data and the effective likelihood of a given model, the posterior may vary considerably. This shows that when assessing the impact of a systematic a measurement, it is not sufficient only checking its contribution to the signal with respect to the noise. Simultaneously, it is important to assess the impact of the systematic on the likelihood itself.

Of course, the noise realisation depicted in the right panel of Fig. <ref> is only one possible realisation of many. We find that for other noise realisations from mocks, the difference in the cosmological constraints between anisotropic systematics and standard isotropic systematics models varies as each posterior scatters as shown in Fig. <ref>. Consequently, neglecting variable depth and the PSF shear bias in the forward modelling, while it is present in the data that is analysed, can bias the inferred value of S_8 by up to ∼1σ. As discussed in Sect. <ref>, this is consistent with what we find for the KiDS-1000 data.





§ IMPACT OF THE GAUSSIAN LIKELIHOOD ASSUMPTION

One of the main aims of this work is to determine whether the standard assumption of the likelihood of cosmic shear two-point statistics being Gaussian holds under realistic conditions while including systematics. To test this assumption, we perform the same analysis of a mock data vector using SBI described in Sect. <ref> and Appendix <ref> again with different inference methods.

As shown in Fig. <ref>, we compare the posterior from the SBI analysis learnt from forward simulations based on the anisotropic systematics model to a standard Gaussian likelihood analysis and an SBI analysis which forces the learnt likelihood to be Gaussian. The standard analysis assumes a Gaussian likelihood characterised by the same numerical covariance from 1,600 realisations of the anisotropic systematics model fixed at the fiducial cosmology in Table <ref>. The model evaluations are calculated from analytical theory for the cosmic shear signal as described in Appendix <ref>, and the sampling is done with the  sampler <cit.> over the same 12-dimensional prior space defined in Table <ref>. In contrast, the Gaussian SBI analysis makes the exact same modelling choices as the main SBI analysis described Sect. <ref> and Appendix <ref> with one change: instead of using ensembles of MAFs made up of MADEs as the architecture of the neural density estimators when learning the effective likelihood, we use an ensemble of three single-component Mixture Density Networks (MDNs; ) to learn the likelihood from the 18,000 realisation of the anisotropic systematics model. As a consequence, the ensemble of MDNs forces the learnt likelihood to be described by a single multi-variate Gaussian distribution. Notably, this is different from a standard Gaussian likelihood analysis with a Gaussian likelihood in that the Gaussian SBI analysis allows the covariance of the likelihood to vary with cosmology.

From Fig. <ref>, we can conclude that the full non-Gaussian SBI is in exceptionally good agreement with the SBI analysis which forces the likelihood to be Gaussian but allows for cosmology dependence in the uncertainty (particularly, when the mock data vector includes realistic noise). Thus, the uncertainty as modelled by the anisotropic systematics model can be considered to be approximately Gaussian at a given cosmology, since a Gaussian NDE network is fully capable of learning an unbiased likelihood. We conclude that, for a fixed cosmology, the assumption of a Gaussian likelihood for the pseudo-Cl cosmic shear signal for ℓ>76 holds well. This is in agreement with previous studies that found this to be the case in angular power spectra for ℓ>50 <cit.>.

At the same time, we find that the standard Gaussian likelhood analysis has a noticeably different posterior to both SBI analyses despite the underlying data vector being identical, and the covariance being sampled from the same anisotropic systematics model. In particular, the S_8 marginal is appreciably narrower in the Gaussian likelihood analysis than in the SBI analyses. In the noiseless case shown in Fig. <ref>, we find that the SBI gives S_8 = 0.743^+0.034_-0.031, while the standard Gaussian likelihood analysis gives S_8 =0.733^+0.034_-0.024, i.e. the 1σ confidence interval in S_8 is ∼11% narrower than in the SBI analysis. When adding realistic noise to the data vector from the same seed, we find that the SBI recovers S_8 = 0.762^+0.032_-0.027, while the Gaussian likelihood analysis recovers S_8 = 0.744^+0.031_-0.022, i.e. the 1σ confidence interval in S_8 is also ∼10% narrower than in the SBI analysis. Additionally, when noise is introduced, the Gaussian likelihood analysis appears more biased with respect to the truth (S_8 = 0.756) than the SBI.

The good agreement between the Gaussian likelihood SBI based on MDNs and the non-Gaussian likelihood SBI based on MAFs indicates that both analyses appear to be sensitive to some form of cosmology dependence in the uncertainty model, as that is the only effect which is not accounted for in the standard Gaussian likelihood analysis.

To investigate the origin of this effect further, we sample the learnt likelihood from the full SBI analysis directly at different values of S_8 in Fig. <ref>. In this plot, we show the seven-dimensional score-compressed data space and label each element of the data vector according to the model parameter with which they are most correlated. Then, we vary the value of S_8 while keeping all other parameters fixed at the fiducial cosmology, and plot the likelihood of the data given the chosen model parameters. We find that as S_8 increases the marginal likelihoods of the elements of the score-compressed data vector increase linearly. In particular, the data elements most correlated with ω_c, ω_b, n_s, and h_0 show the strongest linear correlation with S_8. The implication of this is that the anisotropic systematics model entails that data  at higher assumed S_8 values has larger uncertainties than data modelled for low S_8 values, and therefore the learnt likelihood is different from a Gaussian likelihood with the covariance evaluated at a fixed cosmology. This effect appears to drive the increase in the uncertainty on S_8 of the SBI-based posteriors with respect to the Gaussian likelihood analysis in Fig. <ref>.

In Appendix <ref>, we conduct a detailed investigation into the physical origin of this observed cosmology dependence in the learnt likelihood of cosmic shear two-point statistics. We find that, as has been noted in previous work <cit.>, the cosmic variance of the cosmic shear two-point signal (at least, for 76 ≤ℓ≤ 1500) is measurably cosmology-dependent, and mostly S_8-dependent. We also determine that this effect is expected to be important enough to become detectable in most tomographic bin combinations of the KiDS-1000 gold sample in the regime of ℓ∼ 10^2. Since the cosmic variance is imprinted into the forward simulations at the step of sampling log-normal random matter fields (see Sect. <ref>), the cosmology dependence of the cosmic variance is present in any forward model presented in this work irrespective of the observational systematic effects considered (i.e. both the anisotropic systematics and the standard isotropic systematics model are affected by this).

With this being the case, the noise of the measured cosmic shear pseudo-Cls may contain some information which may impact the inferred cosmology from the SBI. Hence, it is important to validate whether the variance realised by log-normal random simulations does not bias the cosmological inference. Firstly, it has been found in previous studies that log-normal simulations can incorporate accurate higher-order statistics into the sampled matter, galaxy and shear fields when compared to N-body dark matter simulations <cit.>, so log-normal fields can incorporate three- or four-point statistics to model the variance of two-point statistics. Additionally, we find in our own testing that a sample covariance of cosmic shear two-point statistics in KiDS-1000 as generated from log-normal simulations agrees within 10% with analytical covariance matrices computed from analytical considerations (, Reischke et al. in prep.), so the variance of the cosmic shear signal as modelled within KiDS-SBI is consistent with analytical theory.



§ COSMOLOGICAL INFERENCE FROM KIDS-1000 DATA



 §.§ Blinding

Throughout this analysis, we followed a procedure to avoid any form of unconscious bias when making choices within the KiDS-SBI models. Although the KiDS-1000 two-point shear statistics have been measured before, we fixed the forward model used in the analysis before ever evaluating learnt likelihood function for the KiDS data.

All authors did not conduct any measurements on the KiDS-1000 data nor was it analysed with the SBI pipeline until the forward model as well as all validation tests were finalised and fixed. Finally, upon evaluating the learnt likelihood for the KiDS data, we initially inferred an anomalously low S_8 with the SBI and the standard MCMC caused by an inconsistency between the measurement code for the forward simulations and the real data. Upon re-measuring the KiDS data in a consistent fashion, we obtain the pseudo-Cls from the KiDS-1000 gold sample as shown in Fig. <ref> and ran it through the entire KiDS-SBI pipeline without making any adjustments to the methodology described in this paper.



 §.§ Cosmological parameter constraints

With the learnt likelihood found to be unbiased, and with an understanding of the sensitivity of the likelihood to systematics and cosmology, we proceed to analyse the cosmic shear pseudo-Cls measured from the KiDS-1000 gold sample described in Sect. <ref> (see Fig. <ref> for the full measurement).





Our main results are based on the likelihood learnt from the anisotropic systematics model in KiDS-SBI which incorporates variable depth and additional shear biases when compared to previous KiDS analyses (e.g. ). The marginalised posterior and the constraints on the main parameters are shown in Fig. <ref> and Table <ref>, respectively (see Appendix <ref> for the full posteriors and all parameter best estimates). We report the maximum a posteriori (MAP) ± the projected joint highest posterior density (PJ-HPD; ) of 0.743^+0.015_-0.051, and a mean marginal with 68% confidence intervals of 0.731± 0.033. We note that Fig. <ref> also shows that the SBI significantly broadens the posterior in Ω_m with respect to a Gaussian likelihood analysis, while the mean value of Ω_m from SBI is highly consistent with the value measured from the early Universe with Planck 2018 TT,TE,EE+lowE <cit.>.

As shown in Fig. <ref>, we find that the MAP best-fit cosmology given the anisotropic systematics model provides a good fit to the data with a probability-to-exceed (PTE) of 0.42. Additionally, we plot the line of best fit at the MAP cosmology in Fig. <ref> and find that it is in good agreement with the measured data vector and the analytical theory at the MAP cosmology.

In addition to our main results based on the anisotropic systematics model, we re-analyse the KiDS-1000 gold sample data with the SBI pipeline based on the standard isotropic systematics model which makes similar assumptions to previous KiDS-1000 analyses and does not consider the effects of variable depth and some shear biases (see Sect. <ref> for details). When considering this model instead, we find a MAP±PJ-HPD of S_8 = 0.780^+0.020_-0.048, and a mean marginal with 68% confidence intervals of 0.772^+0.038_-0.032, while the fit to the data at the MAP is similarly good with a χ^2 that gives a PTE of 0.52. With the standard isotropic systematics model, the marginal in S_8 is shifted by 0.9 σ upwards. This is broadly consistent with the expectation from the validation discussed in Sect. <ref>, i.e. that the extra systematics in the anisotropic systematics model can bias the inference when the data is noisy despite the bias in the signal itself being negligible at the precision of KiDS-1000. This is driven by the uncertainty modelling being considerably different in both models (see Sect. <ref>): the shape noise model is different, variable depth changes the correlation functions due to the anisotropic selection, and the inclusion of additional shear biases adds statistical noise from 20 more parameters. Consequently, both models respond differently to the noise measured in the KiDS-1000 data.

We also conduct a Bayesian model comparison between the two models by evaluating the Bayesian evidence, Z, for each posterior using  <cit.>. We find that for the anisotropic systematics model, ln(Z) = -14.18, while for the standard isotropic systematics model, ln(Z) = -15.07. Hence, assuming equal prior probability of both models, the log of the Bayes factor, |ln(B)| = 1.11 which implies no significant preference for either model when compared to the other.

As also shown in Fig. <ref>, we find that when re-analysing the KiDS-1000 data assuming a Gaussian likelihood instead, we find results consistent with the testing conducted on mocks (see Fig. <ref>): the Gaussian likelihood assumption causes S_8 to be overconstrained by ∼10% with respect to the SBI irrespective of how the observational systematics are modelled as it neglects to consider the cosmology dependence of cosmic variance. Additionally, the posterior mean marginal in S_8 assuming the Gaussian likelihood shifts from 0.764^+0.031_-0.025 when deriving the sample covariance from the standard isotropic systematics model to a value of 0.725^+0.034_-0.023 when assuming the anisotropic systematics model in the covariance instead. Note that both assume the same modelling in the signal described in Appendix <ref>.

We highlight that the Gaussian likelihood analysis based on the standard isotropic systematics model agrees  well with the constraints from <cit.> where S_8 = 0.758^+0.017_-0.026. We note that the constraints in our work are broader despite making consistent modelling assumptions. This is caused by the use of a different MCMC sampler. Here we have employed the  sampler <cit.>, while in <cit.> the  sampler <cit.> was used. It has been found that  can be overconfident by approximately ∼10% <cit.> which we also find in our analysis and show in Appendix <ref>. When accounting for the difference in sampler used, we find that the Gaussian likelihood constraints from the isotropic systematics model agree exceptionally well with the results from <cit.>.

Likewise the inferred S_8 value when assuming a Gaussian likelihood derived from the anisotropic systematics model is consistent with the result from <cit.> where the mean marginal is measured to be 0.742^+0.034_-0.023. Despite both analyses being based on pseudo-Cls, the one presented in this work finds a slightly lower value of S_8 (∼0.27σ) which can be driven by the additional signal and uncertainty induced by the PSF variations into the anisotropic systematics model (see Appendix <ref>) which was not considered in the modelling of <cit.>.

When comparing the constraints from the anisotropic systematics model SBI to previous cosmic shear analyses in Fig. <ref>, we find that our results are in agreement within 1σ of S_8 best-fits from all previous KiDS-1000 analyses <cit.>, the HSC-Y3 pseudo-Cl analysis <cit.>, and the DES-Y3 two-point analysis <cit.>. The one exception is the anisotropic systematics DES-Y3 analysis <cit.> whose marginal from the main hybrid analysis is different by 1.65σ from our marginal estimate of S_8 based on the anisotropic systematics model. However, in the same work, all re-analyses of the KiDS-1000 data alone (with KiDS-like, DES-like, and hybrid analysis) are consistent within less than 1σ with the best-fit for S_8 presented in this work. These discrepancies may be partially attributed to cosmic variance, different modelling choices for certain systematics, and the impact of the cosmology dependence in the learnt likelihood which is included in our SBI analysis, but it is not considered in any other previous analyses.

In fact, as discussed in detail in Sect. <ref>, since the uncertainty is cosmology-dependent and scales with S_8 in KiDS-SBI, the SBI infers a larger uncertainty in S_8 on top of the added uncertainty from variable depth and shear biases. Additionally, variable depth and the PSF shear bias can impact the likelihood to significantly re-weight the data vector despite their impact on the signal being negligible. This may explain why the best estimate for S_8 from the anisotropic systematics model is about 0.5σ lower than previous standard analyses of KiDS-1000.

When comparing the results from the SBI with the anisotropic systematics model to early-Universe probes of large-scale structure, such as the results from the Planck collaboration shown in Fig. <ref>, we find that the main KiDS-SBI S_8 marginal is in 2.9 σ tension with the best estimate of S_8 = 0.834 ± 0.016 from Planck 2018 TT,TE,EE+lowE <cit.>.





§ CONCLUSIONS

We performed a novel simulation-based inference (SBI) analysis, KiDS-SBI, of the cosmic shear two-point statistics as measured by the Kilo-Degree Survey's fourth data release (KiDS-1000; ).

We presented a new suite of stochastic forward simulations of cosmic shear observables within KiDS-SBI which take into consideration all the systematic effects which are typically modelled in a weak lensing analysis plus additional systematics which have the potential to significantly bias the signal. At the same time, KiDS-SBI introduces considerable improvements in accuracy and performance with respect to similar simulation environments in previous work <cit.>: non-Limber projections of angular power spectra (with the new  module), higher resolution along the line-of-sight, inclusion of field-level intrinsic alignments, consideration of variable depth and anisotropic shear biases, and significant computational performance improvements.

On the basis of these forward simulations, we implement the first simulation-based inference analysis of large-scale structure observables as measured by KiDS-1000 with the same level of complexity as other standard analyses in the field. Thus, this enables us to fully propagate all uncertainties from the data vector all the way to the inferred cosmological parameters in a Bayesian way. The analysis is based on neural density estimators which learn the effective likelihood through the mapping between the model parameters and the forward-simulated data vectors. We conduct inference of the cosmological parameters in flat ΛCDM from cosmic shear two-point statistics, specifically, pseudo-Cls. We assume a forward model consistent with previous KiDS-1000 cosmic shear analyses <cit.> which additionally includes variable depth and anisotropic shear biases: the anisotropic systematics model. Variable depth accounts for the anisotropy in the galaxy selection function in KiDS-1000 and how it impacts the measured galaxy densities, shapes and redshifts. The inclusion of anisotropic shear biases models the angular variation of the residual bias in the galaxy shapes due to variations in the point-spread function of the telescope from the mean model used in the shape calibration.

To assess the quality of any best fit obtained from KiDS-SBI, we also introduce a Bayesian goodness-of-fit measure based on <cit.>. From a small number of realisations of the forward model, we characterise a distribution of χ^2 values for the given model. By determining where the χ^2 of the best fit on the KiDS-1000 data lies, we evaluate the quality of the fit, and whether the data is a feasible realisation of the forward model. From this, we conclude that the KiDS-1000 data is fitted well by the anisotropic systematics model at the maximum a posteriori cosmology (with a probability-to-exceed of 0.42).

Through the testing conducted in this work and in a previous paper <cit.>, we determine that KiDS-SBI recovers accurate and robust posteriors of the 12-dimensional space of cosmological and astrophysical parameters. At the same time, we find that the likelihood learnt from KiDS-SBI forward simulations is well approximated by a Gaussian distribution at a fixed cosmology. However, we determine that the uncertainty on the cosmic shear signal measurably varies with the assumed cosmological parameters which are driven by the cosmic variance scaling with the underlying value of S_8 as has been noted in previous work <cit.>. Consequently, we find that the uncertainty of the two-point statistics in KiDS-1000 is measurably cosmology-dependent which increases the nominal uncertainty of S_8 by ∼10%. Hence, the standard assumption of the covariance of a Gaussian likelihood being fixed for all parameter points in the prior space does not necessarily hold. Therefore, such an assumption may not be valid for current Stage-III galaxy surveys, and it may be also be relevant for upcoming Stage-IV galaxy surveys, such as Euclid <cit.> or Rubin <cit.>, which will further improve the precision of cosmic shear measurements.

Given the anisotropic systematics model in KiDS-SBI, we obtain a mean marginal with 68% confidence intervals of 0.731± 0.033. These constraints are in agreement within 1σ with previous cosmic shear analyses of KiDS-1000 <cit.>, HSC-Y3 <cit.> and DES-Y3 <cit.>. The tension in S_8 from the KiDS-1000 analysis in this work and measurements from the cosmic microwave background from Planck 2018 TT,TE,EE+lowE <cit.> is at 2.9σ, which is similar to the 3.0σ tension found in an equivalent Gaussian likelihood analysis of KiDS-1000 <cit.> or 2.8σ tension found the analysis of KiDS-1000 which incorporated variable depth in the Gaussian likelihood <cit.> or the 3.1σ tension found in the KiDS-1000 beyond ΛCDM <cit.>. This establishes that the σ_8-tension between late-Universe and early-Universe probes of large-scale structure still persists even when one drops the assumption of a Gaussian likelihood for the cosmic shear signal and considers additional contaminant systematics.

To evaluate and isolate the impact of the inclusion of variable depth and shear biases in the anisotropic systematics model, we reproduce the KiDS-SBI analysis with a second model which makes the same assumptions, but does not consider these specific systematic effects: the standard isotropic systematics model (which is more in line with previous KiDS-1000 analyses, e.g. ). We find that the inclusion of both systematics in the forward model increases the S_8 marginal by 0.9σ or 5% and decreases the marginal 1σ confidence intervals by 6%. This reflects the fact that variable depth and PSF shear variation bias can significantly alter the uncertainty modelling despite having a negligible impact on the shear signal at the precision levels of KiDS-1000. The KiDS-1000 data appears to be consistent with the presence of these effects since the anisotropic systematics model provides a good fit to the data, but when examining the Bayesian evidence, we do not find a clear preference between the two models by the data. Nevertheless, both systematic effects have been independently measured in the KiDS-1000 weak lensing data <cit.>.

This work stresses the importance of testing the common assumption of Gaussian likelihoods in realistic conditions. At the same time, we also highlight the potential of simulation-based inference to rigorously perform Bayesian uncertainty propagation even when the models of the uncertainty of the signal or the systematics are analytically intractable. SBI with likelihood-based neural density estimation comes with the added advantage that it does not require the expensive re-evaluation of the posterior distribution if the input measurements or priors are updated, as would be the case for a standard MCMC analysis. Since the likelihood for a given forward model is learnt independently of the input data and the priors, SBI can significantly improve the efficiency of inference analysis where the forward model is not changed between measurements. All these advantages showcase the potential of simulation-based inference becoming the method of choice in future parameter inference and model testing analyses in observational cosmology.

We conclude that SBI is a powerful tool which may be able to tackle many of the physical and statistical challenges which future galaxy surveys will face, as their observations become more precise, less limited by statistical noise, and more limited by the accurate modelling of systematics. To address this need, future work should expand upon the forward model in KiDS-SBI to also include other probes of large-scale structure such as galaxy clustering and galaxy-galaxy lensing, while incorporating the modelling for the relevant systematic effects, such as magnification bias, non-linear galaxy bias and field-level baryonic feedback.


We thank Stephen Feeney, Justin Alsing, Hendrik Hildebrandt, Marika Asgari and Alan Heavens for helpful discussions. MvWK thanks the Science and Technology Facilities Council for support. KL thanks the UCL Cosmoparticle Initiative. NT and BJ are supported by the ERC-selected UKRI Frontier Research Grant EP/Y03015X/1 and the UK Space Agency grants ST/W002574/1 and ST/X00208X/1. BJ is also supported by STFC Consolidated Grant ST/V000780/1. AL acknowledges support from the research project grant `Understanding the Dynamic Universe' funded by the Knut and Alice Wallenberg Foundation under Dnr KAW 2018.0067. RR is supported by the European Research Council (Grant No. 770935). This work used facilities provided by the UCL Cosmoparticle Initiative. We thank Dr Edd Edmonson for technical support at the UCL HPC facilities. The results in this paper are based on observations made with ESO Telescopes at the La Silla Paranal Observatory under programme IDs: 088.D-4013, 092.A-0176, 092.D-0370, 094.D-0417, 177.A-3016, 177.A-3017, 177.A-3018 and 179.A-2004, and on data products produced by the KiDS consortium. The KiDS production team acknowledges support from: Deutsche Forschungsgemeinschaft, ERC, NOVA and NWO-M grants; Target; the University of Padova, and the University Federico II (Naples). Contributions to the data processing for VIKING were made by the VISTA Data Flow System at CASU, Cambridge and WFAU, Edinburgh. The work presented in this paper made use of the following software packages:  <cit.>,  <cit.>,  <cit.>,  <cit.>,  <cit.>,  <cit.>,  <cit.>, and  <cit.>.














aa





§ NUMERICAL NON-LIMBER INTEGRATION
 
Levin's method <cit.> casts a quadrature problem of an oscillatory integral into the solution of a system of ordinary and linear differential equations. Assuming a set basis for the solution the solution to the set of differential equations can be cast into a simple linear algebra problem. It has been used in a variety of cosmic shear applications <cit.>. The method relies on integrals of the type:

    I [h(x,k)] = ∫_x_1^x_2d x  ⟨f( x, k ), w(k,x)⟩ ,

where ⟨·,·⟩ denotes a scalar product, which, in this case, is taken between the vector of non-oscillatory functions f and the oscillatory part w, which has to satisfy:

    dw(x)/dx = Aw(x) ,

with a matrix A. Therefore any oscillatory function with recursion relations is very well suited for this formalism.

The integral in Eq. <ref> is approximated by finding a vector p such that:

    ⟨ p,  w⟩' = ⟨ p' +  A^T p,  w⟩≈⟨ F,  w⟩ .

Thus
⟨ p,  w⟩' = ⟨ F,  w⟩, at n collocation points x_j , j = 1, 2, . . . , n, is subject to the following set of equations:

    ⟨ p' +  A^T  p -  F,  w⟩ (x_j) = 0,    j = 1, ..., n .

Here the trivial solution is the null vector:

    p'(x_j) +  A^T (x_j)  p(x_j) =  F(x_j) .

In a last step a solution is constructed for p using n differentiable basis functions u_m(x) that p can be expanded in:

    p_i(x) = ∑_m c_i^(m) u_m(x),    i= 1,...,d;    m = 1,...,n .

Implying:

    ∑_m c_i^(m)u_m'(x_j) + ∑_m,qA_q i c_q^(m) u_m(x_j) = F_i(x_j) ,

where i,q = 1, ...,d; j,m = 1, ...,n . 
Any basis is suitable, here we chose equidistant collocation points x_j
and the n lowest-order polynomials as basis functions. Finally, the vector w depends on the specific type of integral considered. Eq. <ref> allows for two possibilities: (i) first  integrate over the co-moving distance χ: (ii) first integrate over the wave-vector k. Both approaches have advantages and disadvantages: due to the separation of the χ and χ^' dependence, integrating over the co-moving distance first will be quicker. However, the residual integral over k can still be very oscillatory and must be sampled accordingly, especially with narrow window functions W(χ). The second option is to integrate over k first, this smooths out the oscillations from the Bessel function. However, there are two remaining integrals, thus slowing down the computation, especially for many radial bins. If, however, the radial bins do not overlap, the resulting spectrum for bins further apart than a single bin index will be very small and can be neglected. A more detailed discussion on this is provided in <cit.>.
For integrals of the type

    I_ℓ,k =     ∫dx f_k(x) j_ℓ(kx) , 
    
        I_ℓ,k_1,k_2 =     ∫dx f_k_1,k_2(x) j_ℓ(k_1x)j_ℓ(k_2x) ,

we choose

    w (x) =     [     j_ℓ(xk); j_ℓ + 1(xk) ] ,   A= 
        [      ℓ/x       -k;        k -ℓ + 2/x ] ; 
    
        w (x) =    [     j_ℓ(xk_1)j_ℓ(xk_2);   j_ℓ+1(xk_1)j_ℓ(xk_2);   j_ℓ(xk_1)j_ℓ+1(xk_2); j_ℓ+1(xk_1)j_ℓ+1(xk_2) ] ,  A= 
        [     2ℓ/x     -k_1     -k_2        0;      k_1     -2/x        0    - k_2;      k_2        0     -2/x     -k_1;        0      k_2      k_1 -2ℓ +2/x ] ,

respectively.



§ ACCURACY OF THE COSMIC SHEAR SIGNAL

To test the robustness of the cosmic shear signal as predicted by the log-normal simulations generated in KiDS-SBI in addition to the self-consistency tests already conducted within Sect. <ref>, we perform a comparison of the noise-free pseudo-Cls as measured in KiDS-SBI assuming the standard isotropic systematics model (see Sect. <ref> and <ref>) to the pseudo-Cls predicted from analytical theory (see Appendix <ref>). This comparison is shown for the KiDS-SBI MAP cosmology in Fig. <ref>. When assuming 19 shells and a comoving-volume weighted kernel as discussed in Sect. <ref>, we find per-cent or sub-percent level agreement for all tomographic bin combinations and scales (within 0.5 σ or less). The agreement between the models is also shown to be good in Fig. <ref> where the posterior from the SBI analysis based on KiDS-SBI is unbiased with respect to the posterior sampled with an MCMC while assuming analytical theory. Overall, we find that for small ℓ the forward simulations predict additional power with respect to the pseudo-Cls predicted from theory. This can be explained by discretisation effects along the line-of-sight due to the division of the observable Universe into concentric shells weighted by a window function (see Sect. <ref>), and by the linear interpolation of convergence and shear across  pixels which are built into . Both effects have been shown to produce a similar minor excess in angular power in testing conducted in <cit.>.

When increasing the resolution along the line-of-sight by increasing the number of shells, and when changing the kernel to a uniform window function as suggested in <cit.>, we find that for 26 shells the agreement with the signal predicted by analytical theory is marginally better, but not enough to impact the SBI. When increasing the number of shells further to 31 shells, such that each shell has a width of ∼100 Mpc h^-1, the agreement with analytical theory is poorer for large ℓ and for low ℓ in the autocorrelations, as discretisation effects become important. Consequently, we conclude that the choice of 19 shells allows to accurately forward model the cosmic shear signal in , while still allowing for computational efficiency (see Sect. <ref>). 

These effects are not significant enough to bias the cosmological inference at the precision of KiDS-1000, particularly, since they are more prevalent in tomographic bin combinations with a lower signal-to-noise ratio. However, future simulation-based inference analyses of galaxy surveys with smaller uncertainties (e.g. Stage-IV surveys) based on similar simulations should aim to account for such effects.

Moreover, Fig. <ref> shows that in the bin combinations at higher redshifts (with the highest signal-to-noise ratio) the angular power in the high ℓ modes tends to be slightly underestimated by the standard isotropic systematics model when compared to the analytical theory. We find that this is driven by more cosmic shear signal leaking from the EE modes to the BB modes than predicted by the mixing matrix discussed in Appendix <ref>. This should not bias the SBI analysis, since the KiDS-1000 data and the forward simulations share the same measurement procedure, so mode mixing effects should be intrinsically ingrained in both the real and the forward-modelled data vectors.





§ EFFECTS OF VARIABLE DEPTH AND SHEAR BIASES

To quantify the physical effect of the inclusion of variable depth and shear biases into the anisotropic systematics model as discussed in Sect. <ref> and <ref>, we isolate the bias each systematic causes in the observed two-point cosmic shear signal. First, we forward model variable depth in KiDS-1000 as described in Sect. <ref>, and compare it to a set of forward simulations where the galaxy density, intrinsic galaxy shape dispersion and the photometric redshift distributions are isotropic for each tomographic bin, as is the case in the standard isotropic systematics model. When comparing the cosmic shear signal through the spatial two-point correlations functions (2PCF) measured from  <cit.> in Fig. <ref>, we find that variable depth adds additional two-point correlation signal at all scales and particularly at scales where θ < 100 arcmin. In this case, we measure the cosmic shear signal with 2PCF instead of pseudo-Cls as it simplifies the comparison of the results with previous work. This additional signal is equivalent to about 1% (∼0.2σ) additional non-cosmological shear signal which is in agreement with predictions made for KiDS-1000 from <cit.>. The effect is similar, but slightly smaller in the lower redshift tomographic bin combinations of KiDS-1000 where the shear signal is weaker.



 Additionally, when modelling variable depth, the shape noise becomes anisotropic across the sky as a function of the variations in the effective galaxy density and the intrinsic galaxy shape dispersion. When considering a Gaussian likelihood, this implies that the shape noise scales with σ_ϵ^2/n_eff <cit.>. As discussed in Sect. <ref> both quantities scale linearly with the variable depth measure, σ_rms, which implies that the shape noise in a given tomographic bin either scales with σ_rms^3 in S1, S2 and S5, or with σ_rms^-1 in S3 and S5. Since in KiDS-1000 the distribution of σ_rms is skewed below the mean (see Fig. <ref>), variable depth in fact reduces the overall contribution of the shape noise to the uncertainty of the cosmic shear signal. Nevertheless, as shown in  Figs. <ref> and <ref>, variable depth also causes additional non-cosmological angular correlations in the data at the scales of pixels, pointing overlaps, whole pointings or even over the whole footprint, while also changing the signal along the line-of-sight anisotropically as shown in Fig. <ref>. These additional correlations can significantly alter the covariance of the two-point statistic or may even lead to non-Gaussian noise. Such correlations can exacerbate the uncertainty contributions from cosmic variance and super-sample covariance, which causes the noise predicted in KiDS-1000 by a model including variable depth to ultimately exceed the noise predicted by a model which does not consider the systematic (Reischke et al. in prep.)).



When isolating the effect of additive and PSF variation shear bias on the observed cosmic shear signal as described in Sect. <ref> in Fig. <ref>, we find that, while the additive shear bias has a negligible effect on the cosmic shear signal, the shear bias originating from PSF variations can bias the pseudo-Cls by up to 0.3σ particularly at large scales. Hence, the inclusion of the PSF shear bias can contribute a percentage level non-cosmological signal on certain scales.





§ THEORETICAL SIGNAL MODELLING
 
To provide consistency checks, we also theoretically model the expected signal for a given mock. We model the shape-shape angular power spectrum, C_ϵϵ^(pq)(ℓ), as follows


    C_ϵϵ^(pq)(ℓ) = C_g g^(pq)(ℓ) + C_gI^(pq)(ℓ) + C_Ig^(pq)(ℓ) + C_II^(pq)(ℓ) ,

where ϵ indicates correlations with a galaxy shape field, g indicates correlations with a cosmic shear field, and I indicates correlations between galaxy shapes due to intrinsic alignments. To model each of these terms, we calculate the matter power spectrum using  and proceed to make a Limber projection to obtain the following expression

    C_a b^(p q)(ℓ) = ∫_0^∞dχ/f_k^2(χ)  W_a^(p)(χ)   W_b^(q)(χ)   P_δ(ℓ +1/2/f_k(χ), χ) ,

where W_a is the kernel for a given field, such that a,b ∈{g, I}. For cosmic shear, the kernel is given by 

    W_g^(p)(χ) = 
            3 H_0^2Ω_m/2f_k(χ)/a(χ)∫^χ_hor_χdχ'
            n_S^(p)(χ')  𝑓_k(χ' - χ)/𝑓_k(χ') ,

where n_S^(p) is the redshift distribution of the objects in p^th tomographic bin as shown in Figure <ref>. For consistency with the forward modelling described in Sect. <ref>, we choose an intrinsic alignments kernel in accordance with the non-linear alignment model (NLA; ) given by
    
    W_I^(p)(χ) =
            - A_IAC_1ρ_cr(χ)  Ω_m/D(χ)  n_S^(p)(χ) .


Knowing the full-sky angular power spectra, C^(pq)_ϵϵ, ν(ℓ'; Θ), we can define the pseudo-Cl, C̃^(pq)_ϵϵ, μ(ℓ; Θ), as follows

    C̃^(pq)_ϵϵ, μ(ℓ; Θ) = ∑_ℓ' = 0 ^ℓ'_max∑_ℓ” = 0 ^ℓ”_max∑_ν = 1 ^3∑_ν' = 1 ^3 M_μν', ℓℓ”
        M'^(pq)_ν' ν, ℓ”ℓ' C^(pq)_ϵϵ, ν(ℓ'; Θ) ,


where μ, ν, ν' ∈{1, 2, 3} such that 1 stands for the EE component, 2 for the BB component and 3 for the EB component, M_ℓℓ” is the mixing matrix of the survey mask, W(θ), which does not vary between tomographic bins, and M'^(pq)_ℓ”ℓ' is the mixing matrix of the effective mask imposed by the average random variations in the observed galaxies in a given tomographic bin over many realisations.



The first mixing matrix from the left in Eq. (<ref>), M_ℓℓ”, is the standard block mixing matrix to account for the partial sky coverage in a survey <cit.>. The second mixing matrix, M'^(pq)_ℓ”ℓ', accounts for an additional effect: the mode mixing induced by a mask of uncorrelated noise originating from the random sampling described in Sect. <ref>. In other words, when randomly Poisson sampling galaxies within a discrete pixel k given a probability of P^(m)(N), there will be a non-zero probability of observing zero counts in that given pixel. In essence, all the pixels which are sampled to have zero galaxy counts impose an additional masking on our sample. This can become significant at the resolutions and galaxy densities which we probe in this work. For example, for an  N_side= 2048 and a galaxy density of n=1  arcmin^-2 for a given tomographic bin k (similar to KiDS-1000), P^(m)(N=0) ≈ 0.07. This means that nearly a tenth of all pixels could be "masked" in a given realisation of simulations at a higher resolution. To allow for these cases, we decide to also model the mode mixing induced by the sampling variance. However, for the settings discussed in this paper for the main pipeline of forwards simulations where N_side= 1024, this effect is not so relevant. Even for the sparsest source bin, S1, with an n_gal^(1) = 0.62  arcmin^-2, P^(m)(N=0) ≈ 0.0005, so only a negligible fraction of pixels are masked due to random sample variance.

To provide a general overview, the distinguishing factor between the two mixing matrices is how we define the mask W(θ) which characterises them. For M_μν, ℓℓ”, W(θ) is defined as the mask of the survey footprint, Ω_survey, such that

    W(θ) =  1,    θ∈Ω_survey
     0,    θ∉Ω_survey.


For M^'(pq)_μν, ℓℓ”, rather than defining it using an explicit spatial mask from the pixels masked by each random realisation, W^'(𝑝)(θ), we can define the expectation value of the mixing matrix due to the random sampling of galaxies directly from the angular power spectrum of W^'(𝑝)(θ) as

    𝒲^(pq)'(ℓ)/2ℓ+1=     1/2ℓ+1∑_m=-ℓ^ℓ w^(p)_ℓ m w^(q)*_ℓ m
     =     4 πσ^(pq)/N_pix-1 + 4 π[μ^(pq) - σ^(pq)/N_pix-1] δ^K_ℓ 0,

where μ^(pq) = P^(p)(N>0) P^(q)(N>0), i.e. the product of the mean probabilities of success within each tomographic bin (the probability that a given pixel will be populated by at least one galaxy), while 2σ^(pq)2 = P^(p)2(N>0)(1-P^(p)(N>0))^2+ P^(q)2(N>0)(1-P^(q)(N>0))^2, i.e. the associated standard deviation of the mean, and N_pix is the number of pixels.

As described in Sect. <ref>, the galaxies are Poisson sampled according to the  probability of success is given in Eq. (<ref>). As we are interested in the mean over many iterations, we can assume that the matter overdensities are small, 1+b^(p)δ^(p)≈ 1, so we can rewrite the probability as

    P^(p)(N>0) =     1 - P^(p)(N=0) 
    
        =      1 - e^-⟨ N_m^(p)⟩ (Θ)
    ≈    1 - e^-n^(p)_gal A_pix ,

where n^(p)_gal is the mean galaxy density per tomographic bin p and A_pix is the mean pixel size for a given  resolution. From this, we can calculate both mixing matrices as shown in Fig. <ref>. In practice, we can see that from Fig. <ref> that the mode mixing M^'(pq)_μν, ℓ”ℓ' is negligible when compared to M_μν, ℓℓ” at an N_side=1024. However, we still include the term in the signal modelling to allow the flexibility to choose higher resolutions within KiDS-SBI. 

With Eqs. (<ref>) and (<ref>) defining the selections which define  M_μν, ℓℓ” and M^'(pq)_μν, ℓ”ℓ', respectively, we can compute the matrices as follows <cit.>:

    M_μν, ℓℓ” = [               W_ℓℓ”^++   (W_ℓℓ”^-+ +W_ℓℓ”^+-)                W_ℓℓ”^–;              -W_ℓℓ”^+-   (W_ℓℓ”^++ - W_ℓℓ”^–)               W_ℓℓ”^-+;                W_ℓℓ”^– -(W_ℓℓ”^-+ + W_ℓℓ”^+-)               W_ℓℓ”^++ ] ,


    M^'(pq)_μν, ℓ”ℓ' = [                 W_ℓ”ℓ'^'++   (W_ℓ”ℓ'^'-+ +W_ℓ”ℓ'^'+-)                  W_ℓ”ℓ'^'–;                -W_ℓ”ℓ'^'+-   (W_ℓ”ℓ'^'++ - W_ℓ”ℓ'^'–)                 W_ℓ”ℓ'^'-+;                  W_ℓ”ℓ'^'– -(W_ℓ”ℓ'^'-+ + W_ℓ”ℓ'^'+-)                 W_ℓ”ℓ'^'++ ]^(pq) ,

where each of the Wigner symbols is given by

    W_ℓℓ'^±± = 1/2ℓ+ 1∑_m=-ℓ^ℓ∑_m'=-ℓ'^ℓ' W_ℓℓ' m m'^± (W_ℓ' ℓ m' m^± )^* .


Since the shear field is a spin 2 field

    W_ℓℓ' m m'^+ = 1/2 (_2W_ℓℓ'^m m' + _-2W_ℓℓ' ^m m') ,


    W_ℓℓ' m m'^- = i/2 (_2W_ℓℓ'^m m' - _-2W_ℓℓ' ^m m') ,

where

    _sW_ℓℓ' m m'^+ = ∫dθsℓ' m'(θ) W(θ) s*ℓ m(θ) .


The latter is difficult to compute, so it is more convenient to compute in phase space such that

    ∑_m=-ℓ^ℓ∑_m'=-ℓ'^ℓ'  _sW_ℓℓ'^m m' (_s'W_ℓ' ℓ^m' m)^*
     = (2ℓ+1)(2ℓ' + 1)/4 π∑_ℓ”=0 ^ℓ”_max𝒲_ℓ”[  ℓ ℓ' ℓ”; -s  s  0 ][   ℓ  ℓ'  ℓ”; -s'  s'   0 ] ,

where s ∈{-2, 2} and

    𝒲_ℓ = ∑_m=-ℓ^ℓ w_ℓ m w_ℓ m^* ,

where w_ℓ m are the coefficients of the spherical harmonic transform of the mask, W(θ), such that

    w_ℓ m = ∫dθ W(θ)*ℓ m(θ)
         .







§ PARAMETER RECOVERY


Fig. <ref> shows the posterior distribution of the seven main parameters varied in the analysis as learnt from 18,000 realisations of the anisotropic systematics model. The data vector used for this posterior is a noisy mock data vector generated with the anisotropic systematics model while assuming the input parameters given in Table <ref>. We find that the posteriors from KiDS-SBI accurately recover the true input parameters (see Table <ref> for the inferred parameter values). As shown in Fig. <ref>, the anisotropic systematics posterior mostly only constrains S_8 and A_IA. At the same time, it becomes apparent that the posterior distributions over ω_b, n_s, A_bary and h_0 are mainly prior driven, i.e. the posterior distribution is mostly flat throughout prior space. As the cosmic shear signal is largely degenerate to small changes in these parameters, this is expected.



In summary, we find that the SBI inference pipeline within KiDS-SBI accurately and robustly recovers the underlying posterior distribution of ΛCDM parameters from two-point statistic measurements of the cosmic shear signal in KiDS-1000 when considering a forward model-based log-normal random fields and biased by relevant systematic effects.




§ COSMOLOGY DEPENDENCE OF THE COSMIC SHEAR TWO-POINT STATISTIC COVARIANCE






To determine the origin of the observed cosmology dependence in the likelihood in Fig. <ref>, we investigate whether there is a theoretical justification for the uncertainty in the cosmic shear signal in KiDS-1000 to vary significantly with S_8. First, we assume that the likelihood of a full-sky shear angular power spectrum, C(ℓ), at a given cosmology, is Gaussian, i.e. ℒ∝exp(-1/2χ^2), where χ is the goodness of fit given by

    χ^(p q) = [Ĉ^(pq)(ℓ) - C_fid^(pq)(ℓ)]/σ_fid^(p q)(ℓ) ,

where Ĉ^(pq) is the observed data vector, C_fid^(pq) is the modelled data vector at the fiducial cosmology, and σ_fid^(p q) is the standard deviation of the data at the fiducial cosmology. Assuming that the underlying shear field is described by a Gaussian random field,  σ_fid^(p q) is given by

    ( σ_fid^(p q))^2 = [(C_fid^(pq) (ℓ) ]^2  + C_fid^(pp)(ℓ)   C_fid^(q q) (ℓ)/2 ℓ + 1 + [ N^(p q)(ℓ) ]^2 ,

where the first term on the right-hand side defined the contribution due to cosmic variance, while N^(p q) is the shape noise contribution for a given tomographic bin combination <cit.>. To assess the impact of varying S_8 on such a likelihood distribution, we introduce a toy parameter, a, which linearly scales the model's fiducial data vector similar to S_8 as follows

    C_fid,  a^(pq) (ℓ) = a   C_fid^(pq) (ℓ) .


Hence, we can define χ^(p q)_a as the goodness-of-fit after rescaling the model by combining Eqs. (<ref>), (<ref>) and (<ref>) as follows

    χ_a^(pq) = 1 - a ±χ^(pq) √(1 + [R^(p q)(ℓ) ]^-2/2 ℓ + 1 + [ N^(pq) (ℓ)/C_fid^(pq) (ℓ)]^2)/√(a^2  1 + [R^(p q)(ℓ)]^-2/2l + 1 + [ N^(pq) (ℓ)/C_fid^(pq) (ℓ)]^2) ,

where R^(p q)(ℓ) ≡ C_fid^(pq)(ℓ)/√(C_fid^(pp)(ℓ)   C_fid^(qq) (ℓ)). By calculating the distribution of exp(-1/2χ_a^2) over χ for different values of ℓ and a, we can assess the impact of rescaling the shear angular power spectrum directly on the Gaussian likelihood across different angular scales, ℓ. To facilitate this, one can analytically compute the standard deviation of this distribution, σ_χ_a, as follows

    σ_χ_a^2 = 1 +  (a^2 - 1) (1 + (2ℓ + 1)  [ N^(pq) (ℓ)/C_fid^(pq) (ℓ)]^2/1 + [R^(p q)(ℓ)]^-2)^-1 .


Applying Eq. (<ref>) to the five KiDS-1000 tomographic bins, while assuming a shape noise contribution consistent with the values given in Table <ref>, we obtain Fig. <ref>. We find that, as expected, at large angular scales (ℓ = 500), the likelihood is dominated by the shape noise contribution, so its width does not scale with a under any circumstance. However, at ℓ = 70 and ℓ =100, in some tomographic bin combinations where the signal-to-noise ratio is large enough (i.e. the higher redshift bins), the width of the likelihood distribution scales linearly with the scaling applied to the shear angular power spectrum, a. This implies that the uncertainty of the KiDS-1000 cosmic shear measurements is dependent on the amplitude of the measurement itself, and therefore the underlying cosmology. In fact, as shown in the upper right panel of Fig. <ref>, the broadening of the analytical likelihood with a appears to be generally consistent with the widening observed in the likelihood marginals from the anisotropic systematics model in KiDS-SBI shown in Fig. <ref>.

-0.4em




§ POSTERIOR VARIATIONS WITH NOISE REALISATIONS

As discussed in Sect. <ref>, different uncertainty realisations of the anisotropic systematics model in KiDS-SBI at the same cosmology produces a significant random scatter in the measured final pseudo-Cls. Fig. <ref> shows how the choice of noise realisation impacts the inferred posteriors at a given cosmology. We find that overall the variations in the posteriors between uncertainty realisations are small with the mean marginal of S_8 fluctuating by up to ±0.5σ and the A_IA marginal can fluctuate up to ± 0.75σ.





§ IMPACT OF THE SAMPLER ON THE GAUSSIAN LIKELIHOOD POSTERIORS

To show that the Gaussian likelihood constraints shown in Fig. <ref> and Fig. <ref> obtained from the  sampler <cit.> are consistent with previous KiDS-1000 analyses (e.g. ) which used the  <cit.>, we compare the posteriors obtained from both samplers for the same data vector and covariance matrix in Fig. <ref>. As has been noted in previous work <cit.>, we find that  consistently gives ∼15% narrower S_8 posterior marginals than the ones sampled  irrespective of the model assumed to obtain the sample covariance.





§ KIDS-SBI POSTERIORS AND PARAMETER CONSTRAINTS

To expand on the results of the KiDS-SBI analysis of the KiDS-1000 gold sample presented in Fig. <ref> and Table <ref>, we present the full posterior learnt in KiDS-SBI in  Fig. <ref> for both the anisotropic systematics and standard isotropic systematics models. Table <ref> also provides the parameter best estimates for all cosmological and astrophysical parameters varied in the forward models. As shown in Fig. <ref>, both models leave A_bary, n_s and h_0 unconstrained, while constraining A_IA. The posteriors of both models strongly prefer values of A_IA greater than zero and less than 1.5. Despite the wide range of the A_IA priors, the posterior stays within a physical regime for the parameter which indicates a low chance of residual systematics strongly biasing the parameter constraints.




