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

# KiDS-SBI: Simulation-Based Inference Analysis of KiDS-1000 Cosmic Shear

Maximilian von Wietersheim-Kramsta 112233 Kiyam Lin 11 Nicolas Tessore 11
Benjamin Joachimi 11 Arthur Loureiro 44 5 5 Robert Reischke 6677 Angus H.
Wright 77

(Received xx; accepted xx)

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 Levin’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 $\Lambda$CDM
cosmological parameters from 18,000 realisations. We infer a mean marginal of
the growth of structure parameter
$S_{8}\equiv\sigma_{8}(\Omega_{\mathrm{m}}/0.3)^{0.5}=0.731\pm 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 ${\sim}5\%$. Our results are in agreement with previous analysis of
KiDS-1000 and reinforce a $2.9\sigma$ 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.

###### Key Words.:

gravitational lensing: weak – methods: data analysis – methods: observational

##  1 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,
$\Lambda$CDM. Recent cosmic shear analyses of the data taken in Stage-III
galaxy surveys, such as the Kilo-Degree
Survey111https://kids.strw.leidenuniv.nl/ (KiDS; Kuijken et al. 2019; Asgari
et al. 2021; Heymans et al. 2021; Loureiro et al. 2022; van den Busch et al.
2022; Li et al. 2023a), the Hyper Suprime-Cam
survey222https://hsc.mtk.nao.ac.jp/ssp/ (HSC; Sugiyama et al. 2022; Aihara et
al. 2022; Li et al. 2023b), and the Dark Energy
Survey333https://www.darkenergysurvey.org/ (DES; Gatti et al. 2021; Amon et
al. 2022; Secco et al. 2022), have constrained these cosmological parameters
with unprecedented precision. Upcoming Stage-IV galaxy surveys, such as
_Euclid_ 444https://www.euclid-ec.org/ (Laureijs et al., 2011), Rubin
555https://www.lsst.org/ (LSST Science Collaboration et al., 2009) or _Roman_
666https://roman.gsfc.nasa.gov/ (Spergel et al., 2015), 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}$, $\sigma_{8}$, or its analogue
$S_{8}\equiv\sigma_{8}(\Omega_{\mathrm{m}}/0.3)^{0.5}$, agree well with each
other (Hikage et al., 2019; Asgari et al., 2021; van den Busch et al., 2022;
Amon et al., 2022; Secco et al., 2022; Li et al., 2023b; Dalal et al., 2023)
despite using largely independent methodologies and survey volumes. However,
the constraints from observations of the late Universe (Asgari et al., 2021;
Heymans et al., 2021; van den Busch et al., 2022; Amon et al., 2022; Abbott et
al., 2022; Li et al., 2023a; Dark Energy Survey and Kilo-Degree Survey
Collaboration et al., 2023) disagree by up to ${\sim}$3.4$\sigma$ with the
$\sigma_{8}$ value consistent with the early-Universe observations from the
Cosmic Microwave Background (CMB; Planck Collaboration et al. 2020). Recently,
cluster abundances as measured by SRG/eROSITA (Ghirardini et al., 2024) 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 (Mantz et al.,
2015; Garrel et al., 2022; Chiu et al., 2023; Bocquet et al., 2024). 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. Handley & Lemos 2021;
Dark Energy Survey and Kilo-Degree Survey Collaboration et al. 2023). Such
systematic effects may be physical, e.g. baryonic feedback or the intrinsic
alignments of galaxies (Mandelbaum, 2018; Amon & Efstathiou, 2022; Li et al.,
2023b; Miyatake et al., 2023; Dalal et al., 2023) 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 nth-order statistics depends on up to 2nth-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, $\mathbf{d}$, one would require a number of
simulations $>|\mathbf{d}|$ to get an accurate sample covariance matrix. This
is already the approximate amount of forward simulations needed for an SBI
analysis (Alsing et al., 2018; Lin et al., 2023). 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
(Schneider & Hartlap, 2009; Sellentin & Heavens, 2018; Sellentin et al., 2018;
Taylor et al., 2019; Upham et al., 2021; Hall & Taylor, 2022). Additionally,
when observational and physical systematics are included, the likelihood may
become even less Gaussian (Jeffrey et al., 2021). 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 (Fluri et al., 2022; Euclid Collaboration et al., 2023; Lu et al.,
2023). 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. Asgari et al.
2021; Loureiro et al. 2022), while also harnessing the efficiency of log-
normal forward models of matter density fields which are mostly accurate for
two-point statistics (Hall & Taylor, 2022; Piras et al., 2023).

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 (Guzik & Bernstein, 2005; Shirasaki et
al., 2019). 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 (Heydenreich et al., 2020), while
also affecting the standard deviation of the KiDS-1000 cosmic shear signal (by
$\lesssim 20\%$; Joachimi et al. 2021). 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, Alsing et al. 2019), as it offers good performance and
scales better when increasing the dimensionality of the parameter space of
interest than other SBI methods (Leclercq, 2018; Alsing et al., 2018), such as
Approximate Bayesian Computation (ABC; Rubin 1984; Pritchard et al. 1999; Lin
& Kilbinger 2015; Lin et al. 2016; Beaumont 2019). 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
(Papamakarios & Murray, 2016; Alsing et al., 2018; Gerardi et al., 2021;
Jeffrey et al., 2021; Lemos et al., 2021; Legin et al., 2021; Mishra-Sharma &
Cranmer, 2022; Spurio Mancini et al., 2023; Lin et al., 2023; Chen et al.,
2023; Hu et al., 2022; Lemos et al., 2023b, c; Gatti et al., 2024; Moser et
al., 2024). We highlight the methodology of Gatti et al. (2024) which
conducted a similar analysis to the one outlined in this paper with the DES
Year 3 weak lensing data; see Jeffrey et al. (2024) for results. Whereas the
DES analysis focuses on SBI of two-point and higher-order shear statistics, as
well as cosmologies beyond $\Lambda$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
$\Lambda$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 (Lin et al., 2023, L23 hereafter), we introduced a new
inference pipeline for estimating the cosmological parameters from the cosmic
shear data of KiDS-1000 (Kuijken et al., 2019) using DELFI. In L23, 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 (Joachimi et al., 2021) 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. 2 presents the KiDS data used
for this analysis. Section 3 gives a detailed description of the forward
simulations developed for the SBI analysis. Section 4 details how the galaxy
shapes and summary statistics are measured in both the simulations and the
data. Section 5 describes the setup of the simulation-based inference pipeline
using DELFI based on L23. Section 6 presents the novel measure of goodness-of-
fit for SBI that we use in this paper. In Sect. 7, we conduct consistency
tests to ensure the SBI pipeline is unbiased and robust, while Sect. 8
discusses the validity of a Gaussian likelihood assumption. Section 9 presents
the cosmological inference results from the KiDS-1000 weak lensing data.
Lastly, we give concluding remarks in Sect. 10, and provide additional details
about the analysis in the Appendix.

Figure 1: Spatial map of the KiDS-1000 footprint. The top panel shows a
Mollweide projection of the full KiDS-1000 footprint, while the two panels at
the bottom show zoomed-in Cartesian projections of KiDS-North and KiDS-South
fields, respectively. Figure 2: The redshift distributions of the five
KiDS-1000 tomographic bins. The shaded areas show to limits of each
tomographic bin, while the solid lines show the $n(z)$ of the source galaxies
in each tomographic bin as a function of both redshift, $z$, and comoving
distance, $\chi$ (the latter is derived assuming a Planck 2018 cosmology;
Planck Collaboration et al. 2020). The black dashed lines show the limits of
the spherical matter shells in our forward simulations.

##  2 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 (Kuijken,
2011) which is attached to the 2.6 m VLT Survey Telescope (VST). The survey
covers approximately 1,350 deg2 between two distinct fields known as KiDS
North and South which are fields located in the Northern and Southern Galactic
Caps (with DEC $\approx$ 0º and DEC $\approx$ -30º, respectively). This area
is the same as the one covered by the VISTA Kilo-degree INfrared Galaxy survey
(VIKING; Edge et al. 2013), which means that both surveys together observe
every object with a total of nine photometric bands: $ugriZYJHK_{\rm s}$. The
analysis presented in this paper makes use of DR4 (Kuijken et al., 2019), also
known as KiDS-1000, which covers approximately 1,000 deg2. The effective area
covered by KiDS-1000 shape measurements shown in Fig. 1 extends over 773.3
deg2 (Joachimi et al., 2021). 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 Astro-WISE (McFarland et al., 2013) and the
methodology outlined in Wright et al. (2019) to join the VST and VIKING
photometric measurements. To obtain the shape measurements, Theli (Erben et
al., 2005) is used to process the $r$-band images, while the gravitational
shear estimates are then obtained from lensfit (Miller et al., 2007, 2013;
Fenech Conti et al., 2017). Depending on the shape noise variance and the
ellipticity measurement noise variance, lensfit 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 lensfit 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 lensfit 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 (Giblin et al., 2021). 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 Kuijken
et al. (2019); Wright et al. (2019); Giblin et al. (2021).

The photometric redshifts are estimated from the KiDS 9-band photometry using
Bayesian template-fitting as incorporated within the BPZ code (Benítez, 2000;
Wright et al., 2019, 2020). The photometric redshifts are calibrated with
spectroscopic redshift measurements made of galaxies in the KiDS sample
(Hildebrandt et al., 2021) based on a subsample of $99\%$ completeness. The
mapping from photometric to spectroscopic redshift, $P(z|z_{\mathrm{ph}})$, is
learnt through a self-organising map (Wright et al., 2020; Hildebrandt et al.,
2021). 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. 2. 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 (Asgari et al., 2021; Heymans et al., 2021; van den
Busch et al., 2022). This not only avoids biases in the cosmological parameter
estimates but also provides an additional test on the redshift calibration.

##  3 Forward simulations

The forward simulations for cosmic shear analyses described in this paper are
fully implemented and available in the KiDS-SBI777Kilo-Degree Survey –
Simulation-Based Inference; https://github.com/mwiet/kids_sbi module. This
code is built within the Cosmological Survey Inference System
(CosmoSIS888https://github.com/joezuntz/cosmosis; Zuntz et al. 2015) and it is
based upon the KiDS Cosmology Analysis Pipeline
(KCAP999https://github.com/KiDS-WL/kcap; Joachimi et al. 2021; Asgari et al.
2021; Heymans et al. 2021; Tröster et al. 2021). In addition, the KiDS-SBI
pipeline is generalisable to other cosmic shear analyses.

Cosmological parameters, $\mathbf{\Theta}$ 3D matter power spectrum,
$P_{\delta,\mathrm{nl}}(\mathbf{k},z)$: Sect. 3.1 Baryonic feedback Split
space into concentric shells: $i,j...$: Sect. 3.2 Resolution along the line-
of-sight Non-Limber 2D projection: Sect. 3.3 2D power spectra,
$C_{\delta\delta}^{(ij)}(\ell)$ Generate log-normal random fields: Sect. 3.4
Matter, $\delta^{(i)}(\bm{\theta})$, & convergence fields,
$\kappa^{(i)}(\bm{\theta})$ Intrinsic alignments,
$\kappa_{\mathrm{IA}}^{(i)}(\bm{\theta})$: Sect. 3.5 Sample galaxies: Sect.
3.6 $\&$ 3.7 Angular galaxy positions, $\bm{\theta}$, redshifts,
$z_{\mathrm{true}}$, observed tomographic bin, $p$, and observed ellipticity,
$\epsilon$ Survey characteristics • Mask • Redshift distributions • Shear
biases • Point spread function • Galaxy density • Intrinsic galaxy shapes •
Depth variability: Sect. 3.8 Measure and correct shear: Sect. 4.1 Observed
shear, $\epsilon^{\mathrm{corr}}$ Measure pseudo-Cls, subtract noise and
binning: Sect. 4.2 Pseudo-Cl, $\tilde{C}_{\epsilon\epsilon}^{(pq)}(\ell)$
Figure 3: Flowchart describing the steps in a single forward simulation of
cosmic shear observables from cosmological parameters. The dark blue rounded
boxes represent the inputs and outputs which are given to the simulation-based
inference pipeline. The green slanted boxes represent relevant quantities
which are calculated during the simulation. The grey rectangular boxes show
steps in the calculations, while the blue slanted boxes show any (systematic)
effects which are included. All variables are defined within the respective
sections quoted in the diagram.

The overall outline of the forward simulations is shown in Fig. 3. Firstly,
the simulation is fed a set of cosmological and astrophysical parameters,
$\mathbf{\Theta}$, from a sampler of our choice. All of these parameters are
passed into CAMB (Lewis et al., 2000; Lewis & Challinor, 2002; Howlett et al.,
2012) which calculates the three-dimensional matter power spectrum,
$P_{\delta,\mathrm{nl}}$, at a given cosmology (see also Sect. 3.1 for more
details). We incorporate baryonic feedback using HMCode (Mead et al., 2016).
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, GLASS (Tessore et al., 2023). 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. 3.2 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_{\delta\delta}$. 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
(Levin, 1996) which allows for the efficient integration of the spherical
Bessel functions found in non-Limber integrals (more details are given in
Sect. 3.3).

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 (GLASS101010Generator for Large Scale Structure;
https://github.com/glass-dev/glass, Tessore et al. 2023), which efficiently
generates correlated log-normal random matter fields, $\delta$, for each
shell. Then, it accurately computes the respective convergence field,
$\kappa$, (see Sect. 3.4) while also incorporating an effective convergence
field due to intrinsic galaxy alignments, $\kappa_{\mathrm{IA}}$ (see Sect.
3.5).

Upon the construction of the matter, $\delta$, and convergence fields,
$\kappa$, 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 SALMO module (Joachimi et al., 2021), 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. 3.6. 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 (Asgari et al.,
2021; van den Busch et al., 2022), the shear fields (Porqueres et al., 2022;
Loureiro et al., 2023) 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 (Loureiro et al., 2022; Tröster et al., 2022), we choose observed
angular power spectra, a.k.a. pseudo-Cls (Peebles, 1973), 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 Loureiro et al.
(2022). 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. 4.2.

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.

###  3.1 Cosmology dependence: three-dimensional matter power spectrum

The basis of the forward simulations is the 3-dimensional matter power
spectrum, $P_{\delta,\mathrm{nl}}(k,z;\mathbf{\Theta})$. In the linear regime,
the equal-time three-dimensional matter power spectrum’s cosmology dependence
comes from the following relation

$P_{\delta,\mathrm{l}}(k,z;\mathbf{\Theta})=T^{2}(k;\mathbf{\Theta})\,D^{2}(z;\mathbf{\Theta})\,P_{\mathcal{R}}(k;\mathbf{\Theta})\,,$
(1)

where $k$ is a wavenumber, $z$ is redshift, $\mathbf{\Theta}$ is the set of
cosmological parameters, $P_{\mathcal{R}}(k;\mathbf{\Theta})$ is the
primordial density fluctuation power spectrum, $D$ is the growth factor of
structure, $T$ is the transfer function between the matter overdensity field,
$\mathbf{\delta}(\mathbf{k},z)$, and the primordial curvature fluctuation
field, $\mathcal{R}(\mathbf{k};\mathbf{\Theta})$, such that
$\mathbf{\delta}(\mathbf{k},z)=T(\mathbf{k};\mathbf{\Theta})D(z;\mathbf{\Theta})\mathcal{R}(\mathbf{k};\mathbf{\Theta})$.
To compute this, we make use of CAMB (Lewis et al., 2000; Lewis & Challinor,
2002; Howlett et al., 2012). To stay in line with the main KiDS-1000 analysis
(Joachimi et al., 2021; Asgari et al., 2021; van den Busch et al., 2022), we
assume a normal neutrino hierarchy, while also assuming a fixed sum of
neutrino masses of $\sum m_{\nu}=0.06$ eV/$c^{2}$.

We compute the non-linear matter power spectrum via a non-perturbative model,
HMCode-2016 (Mead et al., 2015, 2016). Although updated iterations of this
model exist (Mead et al., 2021), we limit ourselves to HMCode-2016, so our
model remains comparable to the KiDS-1000 cosmic shear analysis. HMCode 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_{\mathrm{bary}}$, and the halo bloating parameter,
$\eta_{0}$. We also fix the relation between these two parameters to
$\eta_{0}=0.98-0.12A_{\mathrm{bary}}$ in accordance with Joudaki et al. (2018)
and Joachimi et al. (2021), and treat $A_{\mathrm{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. 3.5). Nevertheless, we do not expect this to cause a large
discrepancy for the scales which we are probing in the cosmic shear signal
($\ell<1500$). Field-level implementations of baryonic feedback already exist
(Schneider & Teyssier, 2015; Schneider et al., 2019, 2022; Porqueres et al.,
2023), and they would constitute interesting avenues for future extensions to
KiDS-SBI.

###  3.2 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 GLASS (Tessore et al., 2023) to model the
underlying large-scale structure. Within GLASS, 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\in\\{1,2,3,...,N_{\mathrm{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 (Tessore et al., 2023),

$W^{(i)}(z;\mathbf{\Theta})=\begin{cases}f_{k}^{2}(z;\mathbf{\Theta})/E(z;\mathbf{\Theta}),&\text{if
$z_{i}\leq z<z_{i+1}$\,,}\\\ 0,&\text{otherwise\,,}\end{cases}$ (2)

where $E(z;\mathbf{\Theta})$ is the dimensionless Hubble function and
$f_{k}(z;\mathbf{\Theta})$ is the transverse comoving distance for a given
redshift and cosmology, which is defined as

$f_{k}(z;\mathbf{\Theta})=\begin{cases}\frac{c}{H_{0}\sqrt{\Omega_{k}}}\mathrm{sinh}\bigg{(}\frac{\sqrt{\Omega_{k}}H_{0}}{c}\chi(z;\mathbf{\Theta})\bigg{)},&\text{if
$\Omega_{k}>0$\,,}\\\ \chi(z;\mathbf{\Theta}),&\text{if $\Omega_{k}=0$}\,,\\\
\frac{c}{H_{0}\sqrt{|\Omega_{k}|}}\mathrm{sin}\bigg{(}\frac{\sqrt{|\Omega_{k}|}H_{0}}{c}\chi(z;\mathbf{\Theta})\bigg{)},&\text{if
$\Omega_{k}<0$\,,}\end{cases}$ (3)

where $H_{0}$ is the Hubble constant, $\Omega_{k}$ is the curvature density
parameter and $c$ is the speed of light in a vacuum. As noted in Tessore et
al. (2023), different choice of window function may induce percent-level
biases in the resulting shear power spectrum. We test for this in Appendix B,
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 D).

A given shell should not be thicker than a few hundred $\mathrm{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
(Xavier et al., 2016; Hall & Taylor, 2022; Tessore et al., 2023). Otherwise,
discretisation effects would smooth away large-scale structure (Tessore et
al., 2023). 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_{\mathrm{shells}}$. Simultaneously, the
simulations’ run-time for a given set of cosmological parameters increases
with $N_{\mathrm{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\lesssim z\lesssim 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. 2).

###  3.3 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. 3.2. 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, $\mathbf{\Theta}$, and a given shell combination $(ij)$ as follows

$C_{\delta\delta,\ell}^{(ij)}(\mathbf{\Theta})=\langle\tilde{\delta}^{(i)}_{\ell
m}(\mathbf{\Theta})\,\tilde{\delta}^{(j)*}_{\ell
m}(\mathbf{\Theta})\rangle\,,$ (4)

where $\tilde{\delta}^{(i)}_{\ell m}(\mathbf{\Theta})$ represents the harmonic
coefficients defined through a spherical harmonic transform of the projection
of the spin-0 matter field, $\delta^{(i)}$, as given by (Zaldarriaga & Seljak,
1997; Reinecke, 2011),

$\tilde{\delta}^{(i)}_{\ell
m}(\mathbf{\Theta})=\int\mathrm{d}^{2}\bm{\theta}\,\delta^{(i)}(\bm{\theta};\mathbf{\Theta})\,\,_{0}\\!\\!\operatorname{Y}^{*}_{\ell
m}(\bm{\theta})\,,$ (5)

where $\,{}_{0}\\!\\!\operatorname{Y}_{\ell m}(\bm{\theta})$ are the spherical
harmonics for a spin value of 0 as a function of the two-dimensional spatial
sky position $\bm{\theta}$. The projected matter field,
$\delta^{(i)}(\bm{\theta};\mathbf{\Theta})$, can in turn be defined with
respect to the underlying three-dimensional matter overdensity field,
$\delta(\bm{\theta},z;\mathbf{\Theta})$, as

$\delta^{(i)}(\bm{\theta};\mathbf{\Theta})=\sum_{\ell
m}\tilde{\delta}^{(i)}_{\ell
m}(\mathbf{\Theta})\,_{0}\\!\\!\operatorname{Y}_{\ell
m}(\bm{\theta})=\int\mathrm{d}z\,W^{(i)}(z;\mathbf{\Theta})\,\delta(\bm{\theta},z;\mathbf{\Theta})\,,$
(6)

where $W^{(i)}(z;\mathbf{\Theta})$ is the weight function of the given shell
as defined in Eq. (2).

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

$\langle\delta(\mathbf{k},z;\mathbf{\Theta})\,\delta^{*}(\mathbf{k}^{\prime},z^{\prime};\mathbf{\Theta})\rangle=(2\pi)^{3}\delta_{\mathrm{D}}(\mathbf{k}-\mathbf{k}^{\prime})\,P_{\delta,\mathrm{nl}}(k,z,z^{\prime};\mathbf{\Theta})\,,$
(7)

where $\delta(\mathbf{k},z;\mathbf{\Theta})$ is the Fourier transform of the
two-dimensional matter overdensity field at $z$,
$\delta(\bm{\theta},z,;\mathbf{\Theta})$, while
$\delta_{\mathrm{D}}(\mathbf{k}-\mathbf{k}^{\prime})$ is the Dirac delta
function.

If we combine the relations in Eqs. (4) to (7) and, for the sake of
simplicity, assume that $\Omega_{k}=0$, i.e.
$f_{k}(z;\mathbf{\Theta})=\chi(z;\mathbf{\Theta})$, we can characterise the
cosmology dependence of $C_{\delta\delta}^{(ij)}(\ell;\mathbf{\Theta})$
directly with respect to the three-dimensional matter power spectrum,
$P_{\delta,\mathrm{nl}}(k,z,z^{\prime};\mathbf{\Theta})$, using the following
relation

$C_{\delta\delta}^{(ij)}(\ell;\mathbf{\Theta})=\frac{2}{\pi}\int\mathrm{d}\chi\,W^{(i)}(z[\chi];\mathbf{\Theta})\int\mathrm{d}\chi^{\prime}\,W^{(j)}(z[\chi^{\prime}];\mathbf{\Theta})\\\
\int\mathrm{d}k\,k^{2}\,P_{\delta,\mathrm{nl}}(k,z[\chi],z[\chi^{\prime}];\mathbf{\Theta})\,j_{\ell}(k\chi)\,j_{\ell}(k\chi^{\prime})\,,$
(8)

where $\chi^{\prime}\equiv z[\chi^{\prime}]$, $j_{\ell}(k\chi)$ are spherical
Bessel functions of order $\ell$, and $\ell\in\mathbb{Z}^{0+}$. We can
simplify Eq. (8) by making the geometric approximation, i.e. taking the
geometric mean of two equal-time power spectra (Castro et al., 2005; Kitching
& Heavens, 2017; Kilbinger et al., 2017; de la Bella et al., 2021), such that

$P_{\delta,\mathrm{nl}}(k,z[\chi],z[\chi^{\prime}];\mathbf{\Theta})=\sqrt{P_{\delta,\mathrm{nl}}(k,z[\chi];\mathbf{\Theta})\,P_{\delta,\mathrm{nl}}(k,z[\chi^{\prime}];\mathbf{\Theta})}\,,$
(9)

which has been shown to be an accurate approximation both in the linear and
non-linear regimes (Kitching & Heavens, 2017).

Figure 4: The two-dimensional angular matter power spectra,
$C_{\delta\delta}^{(ij)}(\ell)$, projected with Levin describing the
correlations between the large-scale structure within a set of 19 concentric
shells spaced along the line-of-sight as shown in Fig. 2. The left panel shows
the autocorrelations of all shells, i.e. $|i-j|=0$. The right panel shows the
correlations of each shell with its nearest neighbour, i.e. $|i-j|=1$
($|i-j|>1$ are not shown as they do not have a large effect on matter fields
within each shell, see Tessore et al. 2023). The colour of each line is given
by the mean redshift of the $i^{\mathrm{th}}$ bin, $\overline{z}^{(i)}$. The
underlying linear three-dimensional matter power spectrum is based on flat
$\Lambda$CDM assuming $\Omega_{\mathrm{c}}=0.05$, $\Omega_{\mathrm{b}}=0.28$,
$\sigma_{8}=0.79$, $S_{8}=0.84$ and $H_{0}=67$ km s-1 Mpc-1 calculated using
CAMB (Lewis et al., 2000; Lewis & Challinor, 2002; Howlett et al., 2012),
while the non-linear contribution is calculated using HMcode-2016 (Mead et
al., 2015, 2016) assuming $A_{\mathrm{bary}}=3.1$. The non-Limber projection
done by Levin assumes $\ell_{\mathrm{max,nL}}=1200$ which causes the
$C_{\delta\delta}^{(ij)}$ with $|i-j|>0$ go to a value of zero for
$\ell>\ell_{\mathrm{max,nL}}$.

Additionally, it is common to make the so-called Limber approximation (Limber,
1953; Kaiser, 1992) which consists in approximating a spherical Bessel
function as a single peak or delta function at $k\chi=\ell+1/2$. This
approximation is reasonably accurate in the auto-correlations for broad window
functions and high angular modes,
$C_{\delta\delta}^{(ii)}(\ell;\mathbf{\Theta})$, 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_{\delta\delta}^{(ij)}$, where $i\neq j$, will be zero as long as their
associated weight functions, $W^{(i)}(z[\chi];\mathbf{\Theta})$, are non-
overlapping, as is the case for the weights constructed in Eq. (2). This can
be seen by looking at the full Taylor expansion of Eq. (8) as given by
(LoVerde & Afshordi, 2008)

$C_{\delta\delta}^{(ij)}(\ell)=\int\frac{\mathrm{d}\chi}{\chi^{2}}\,W^{(i)}(z[\chi])\,W^{(j)}(z[\chi])\,P_{\delta,\mathrm{nl}}\Bigg{(}\frac{\ell+1/2}{\chi},z[\chi]\Bigg{)}\\\
\Bigg{\\{}1-\frac{1}{(\ell+1/2)^{2}}\Bigg{[}\frac{\chi^{2}}{2}\Bigg{(}\frac{W^{(i)^{\prime\prime}}(z[\chi])}{W^{(i)}(z[\chi])}\,+\frac{W^{(j)^{\prime\prime}}(z[\chi])}{W^{(j)}(z[\chi])}\Bigg{)}\,\\\
+\frac{\chi^{3}}{6}\Bigg{(}\frac{W^{(i)^{\prime\prime\prime}}(z[\chi])}{W^{(i)}(z[\chi])}+\frac{W^{(j)^{\prime\prime\prime}}(z[\chi])}{W^{(j)}(z[\chi])}\Bigg{)}\Bigg{]}+\mathcal{O}\Big{(}[\ell+1/2]^{-4}\Big{)}\Bigg{\\}}\,,$
(10)

where each apostrophe, ′, denotes a partial derivative with respect to $\chi$,
i.e. $\partial/\partial\chi$, 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 $\ell$ 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_{\delta\delta}^{(ii)}$ or less. In this case, the simulated
$\delta^{(i)}(\bm{\theta};\mathrm{\Theta})$ 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. (8). 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 Levin module111111Levin;
https://github.com/rreischke/nonLimber_matter_shells. This is a novel Python
module which implements the Levin integration method for oscillatory functions
(Levin, 1996) in the context of non-Limber integrals for weak lensing and
galaxy clustering which is discussed in more detail in Appendix A.

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_{\delta\delta}^{(ij)}(\ell)$ into two domains: non-Limber
($1<\ell\leq\ell_{\mathrm{max,nL}}$) and second-order extended Limber given in
Eq. (10) for $\ell_{\mathrm{max,nL}}<\ell\leq 30,000$. Over the former domain,
Levin numerically integrates the expression shown in Eq. (8) for the given
matter power spectrum and the weights. Over the latter domain, we perform the
numerical integration of Eq. (10) up to the second-order term in
$(\ell+1/2)^{-1}$. This reduces the computational resources needed for the
integration without any substantial loss of accuracy as long as
$\ell_{\mathrm{max,nL}}$ is sufficiently large, since the residuals scale with
$(\ell+1/2)^{-4}$ (LoVerde & Afshordi, 2008). With the shells shown in Fig. 2
weighted with Eq. (2), we obtain angular power spectra,
$C_{\delta\delta}^{(ij)}(\ell;\mathbf{\Theta})$, shown in Fig. 4.

With this configuration, we find that we can obtain computationally efficient
and accurate (consistent with CAMB 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 Levin
module in previous applications (Zieser & Merkel, 2016; Spurio Mancini et al.,
2018a, b; Baleato Lizancos & White, 2023).

###  3.4 Log-normal matter field simulations

With the input angular power spectra for a given cosmology computed as
described in Appendix A, the next step in the forward simulations involves
sampling random matter fields consistent with these correlations. These fields
are generated with GLASS (Tessore et al. 2023). 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 GLASS 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,
$\delta^{(i)}(\bm{\theta};\mathbf{\Theta})$, as well as the associated
convergence fields for each shell,
$\kappa^{(i)}(\bm{\theta};\mathbf{\Theta})$. These are then integrated along
the line-of-sight to give the convergence field for a source within a given
shell, $\kappa^{(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
(Hall & Taylor, 2022; Piras et al., 2023). For this reason, log-normal random
fields are a common approximation used for matter overdensity fields (Coles &
Jones, 1991; Böhm et al., 2017; Abramo et al., 2016, 2022) as well as for
convergence fields (Hilbert et al., 2011; Clerkin et al., 2017; Giocoli et
al., 2017; Gatti et al., 2020).

Log-normal fields within GLASS 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,
$\delta_{\mathrm{log}}^{(i)}(\bm{\theta};\mathbf{\Theta})$, as follows

$\langle\delta_{\mathrm{log}}^{(i)}(\bm{\theta};\mathbf{\Theta})\,\delta_{\mathrm{log}}^{(j)*}(\bm{\theta};\mathbf{\Theta})\rangle=C_{\delta\delta}^{(ij)}(\theta;\mathbf{\Theta})\,.$
(11)

By defining the mean and the variance for
$G_{\delta\delta}^{(ij)}(\theta;\mathbf{\Theta})$ as

$\langle\delta_{\mathrm{G}}^{(i)}(\bm{\theta};\mathbf{\Theta})\,\delta_{\mathrm{G}}^{(j)*}(\bm{\theta};\mathbf{\Theta})\rangle=G_{\delta\delta}^{(ij)}(\theta;\mathbf{\Theta})\,,$
(12)

one can then sample $\delta_{\mathrm{G}}^{(i)}(\bm{\theta};\mathbf{\Theta})$
for each shell. To achieve this with numerical methods, the two-dimensional
spherical shell of each layer is discretised using HEALPix pixels (Górski et
al., 2005). This means that $\bm{\theta}$ is discretised such that it can be
mapped to a linearised variable given by $\theta_{m}$ where
$m\in\\{1,2,3,...,12N_{\mathrm{side}}^{2}\\}$ and $N_{\mathrm{side}}$ is the
HEALPix resolution parameter. To obtain the final matter overdensity fields,
we use the following expression (Coles & Jones, 1991; Kayo et al., 2001;
Hilbert et al., 2011; Xavier et al., 2016),

$\delta_{\mathrm{log}}^{(i)}(\theta_{m};\mathbf{\Theta})=e^{\delta_{\mathrm{G}}^{(i)}(\theta_{m};\mathbf{\Theta})}-\lambda\,,$
(13)

where $\lambda$ is the shift of the log-normal distribution which is assumed
to be $\lambda=1$ for matter fields. For more details about this calculation,
see Tessore et al. 2023.

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. 3.2, if one chooses sufficiently large
$\ell_{\mathrm{max}}$ up to which the input angular power spectra are
calculated, one can accurately sample
$\delta_{\mathrm{log}}^{(i)}(\theta_{m};\mathbf{\Theta})$ up to an
$N_{\mathrm{side}}\sim 0.5\ell_{\mathrm{max}}$ (Leistedt et al., 2013; Alonso
et al., 2019).

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
$\bm{\theta}$ and $z$ along a continuous line-of-sight under the Born
approximation (Bartelmann & Schneider, 2001),

$\begin{split}\kappa(\bm{\theta},\mathit{z};\mathbf{\Theta})=\frac{3\Omega_{\mathrm{m}}}{2}\int_{0}^{\mathit{z}}\mathrm{d}\mathit{z}^{\prime}&\,\frac{\mathit{f}_{\mathit{k}}(\mathit{z}^{\prime};\mathbf{\Theta})\,[\mathit{f}_{\mathit{k}}(\mathit{z}^{\prime};\mathbf{\Theta})-\mathit{f}_{\mathit{k}}(\mathit{z};\mathbf{\Theta})]}{\mathit{f}_{\mathit{k}}(\mathit{z};\mathbf{\Theta})}\\\
&\times\frac{1+\mathit{z}^{\prime}}{E(\mathit{z}^{\prime};\mathbf{\Theta})}\,\delta(\bm{\theta},\mathit{z}^{\prime};\mathbf{\Theta})\,,\end{split}$
(14)

where $\Omega_{\mathrm{m}}$ is the matter density fraction at $z=0$, and
$\Omega_{\mathrm{m}}\in\mathbf{\Theta}$. We then discretise the continuous
two-dimensional matter overdensity field along over $\bm{\theta}$ using
HEALPix pixels, while also discretising along the line-of-sight using Eq. (6)
and defining a mean redshift, $\overline{z}^{(i)}$, for a given weighted shell
as follows

$\overline{\mathit{z}}^{(i)}(\mathbf{\Theta})=\frac{\int\mathrm{d}\mathit{z}\,\mathit{z}\,\mathit{W}^{(i)}(\mathit{z};\mathbf{\Theta})}{\int\mathrm{d}\mathit{z}\,\mathit{W}^{(i)}(\mathit{z});\mathbf{\Theta}}\,.$
(15)

To accelerate the calculation of the discretised convergence fields for each
shell, $\kappa^{(i)}$, within GLASS, 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 (Schneider, 2016). This gives a recurrence
relation for the convergence field, $\kappa^{(i)}$, within shells where $i\geq
2$, which only depends on $\kappa^{(i-1)}$, $\kappa^{(i-2)}$ and
$\delta^{(i-1)}$ (see Tessore et al. 2023 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, $\mathbf{\Theta}$, needed for this
analysis can be obtained with fewer resources when compared to forward models
in previous analysis (see Sect. 5.1).

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.

###  3.5 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 (Heavens et al., 2000a; King & Schneider, 2002;
Heymans & Heavens, 2003; Bridle & King, 2007).

To model this effect, we follow the methodology set out in Tessore et al.
(2023). In this prescription, for a given intrinsic alignment model, we define
an associated effective convergence field, $\kappa^{(i)}_{\mathrm{IA}}$. This
$\kappa^{(i)}_{\mathrm{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
$\kappa^{(i)}_{\mathrm{IA}}\propto\delta^{(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

$\kappa^{(i)}(\theta_{m};\mathbf{\Theta})\rightarrow\kappa^{(i)}(\theta_{m};\mathbf{\Theta})+\kappa_{\mathrm{IA}}^{(i)}(\theta_{m};\mathbf{\Theta})\,.$
(16)

For simplicity and consistency with typical modelling assumptions in weak
lensing surveys such as KiDS-1000 (Joachimi et al., 2021), we choose to model
$\kappa^{(i)}_{\mathrm{IA}}$ with the Non-Linear Alignment (NLA) model
(Catelan et al., 2001; Hirata et al., 2004; Bridle & King, 2007), 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 (Bridle & King 2007,
similar to our approach described in Sect. 3.1). This means that it is also
proportional to the local matter overdensity field, which is given by (Hirata
et al., 2004),

$\kappa_{\mathrm{IA}}^{(i)}(\theta_{m};\mathbf{\Theta})=-A_{\mathrm{IA}}\frac{C_{1}\,\Omega_{\mathrm{m}}\,\overline{\rho}_{\mathrm{cr}}(\overline{z}^{(i)};\mathbf{\Theta})}{D(\overline{z}^{(i)};\mathbf{\Theta})}\,\delta^{(i)}(\theta_{m};\mathbf{\Theta})\,,$
(17)

where $A_{\mathrm{IA}}$ is the intrinsic alignments amplitude (which we are
treating as a nuisance parameter that is also sampled, see Sect. 5), $C_{1}$
is a normalisation constant which we set to $C_{1}=5\times
10^{-14}h^{-2}M_{\odot}^{-1}\mathrm{Mpc}^{3}$ in agreement with the IA
measurements at low redshifts by SuperCOSMOS (Brown et al., 2002),
$\overline{\rho}_{\mathrm{cr}}(\overline{z}^{(i)};\mathbf{\Theta})$ is the
mean critical matter density as a function of redshift and
$D(\overline{z}^{(i)};\mathbf{\Theta})$ is the linear growth factor normalised
to be unity at $z=0$. Note that in Eq. (17), we do not consider an explicit
redshift dependence which is often expressed as a power-law term (Joachimi et
al., 2011). We omit this term, as it is not considered in the fiducial
KiDS-1000 analysis (Asgari et al., 2021; van den Busch et al., 2022), and it
has been found that at least for $z<1$, the IA signals do not vary largely
with redshift (Fortuna et al., 2021).

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 (Blazek et al., 2015). This effect can in fact produce a
signal in the shape measurements comparable to the correction due to non-
linear structure growth (Krause et al., 2016). 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 (Joachimi et al., 2011; Blazek et al.,
2011; Heymans et al., 2013; Krause et al., 2016; Hilbert et al., 2017),
including KiDS-1000 (Fortuna et al., 2021).

###  3.6 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, $\Omega_{\mathrm{survey}}$, the galaxy number density which the
survey can observe, $n_{\mathrm{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,
$\sigma_{\mathrm{\epsilon}}$, 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_{\mathrm{shells}}$
concentric shells, they ought to be modelled across $N_{\mathrm{tomo}}$
tomographic bins. The tomographic bins of a given survey are defined by finite
usually non-overlapping domains in photometric redshift, $z_{\mathrm{ph}}$, in
which all observed galaxies are grouped, such that $z_{p}<z_{\mathrm{ph}}\leq
z_{p+1}$, $p\in\\{1,2,3,...,N_{\mathrm{tomo}}\\}$121212Note 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_{\mathrm{ph}})$ by mapping the measured $z_{\mathrm{ph}}$ to known
spectroscopic redshift measurements which have a comparatively low bias in the
redshift estimate (Hildebrandt et al., 2021; van den Busch et al., 2022).
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_{\mathrm{ph}})$.

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

$n_{\mathrm{eff}}^{(p)}\equiv\frac{1}{|\Omega_{\mathrm{survey}}|}\frac{\Big{[}\sum_{i\in(p)}w_{i}\Big{]}^{2}}{\sum_{i\in(p)}w_{i}^{2}}\,,$
(18)

where $|\Omega_{\mathrm{survey}}|$ is the effective area covered by the survey
footprint, $\overline{M}^{(p)}$ is the mean multiplicative shear bias in
tomographic bin $p$ as calibrated in (Giblin et al., 2021), and $w_{i}$ is the
lensfit 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_{\mathrm{gal}}$, when not correcting for the lensfit weights for each
galaxy. However, as KiDS-SBI does not model galaxy weights, in the forward
simulations $n_{\mathrm{eff}}=n_{\mathrm{gal}}$ since $w_{i}=1\forall i$.

Similarly, we define the galaxy intrinsic ellipticity dispersion,
$\sigma_{\mathrm{\epsilon}}$, within a given tomographic bin as follows

$\sigma_{\mathrm{\epsilon}}^{(p)2}\equiv\frac{\sum_{i\in(p)}w_{i}^{2}\,\Big{(}\epsilon_{\mathrm{obs},i,1}^{2}+\epsilon_{\mathrm{obs},i,2}^{2}\Big{)}}{\Big{(}1+\overline{M}^{(p)}\Big{)}^{2}\sum_{i\in(p)}w_{i}^{2}},$
(19)

where $\epsilon_{\mathrm{obs},1}$ and $\mathrm{i}\epsilon_{\mathrm{obs},2}$
are the two components of the observed galaxy ellipticity. Applying Eqs. (18)
and (19) to the KiDS-1000 gold sample, we obtain the values shown in Table 1.

Bin | $\overline{n}_{\mathrm{eff}}$ | $\overline{\sigma}_{\epsilon}/\sqrt{2}$ | $a_{n_{\mathrm{eff}}}$ | $b_{n_{\mathrm{eff}}}$ | $a_{\sigma_{\epsilon}}$ | $b_{\sigma_{\epsilon}}$
---|---|---|---|---|---|---
| [$\mathrm{arcmin}^{-2}$] |  |  |  | $\times 10^{3}$ |
S1 | 0.62 | 0.27 | -0.035 | 0.72 | 1.81 | 0.267
S2 | 1.18 | 0.26 | -0.042 | 1.30 | 1.30 | 0.257
S3 | 1.85 | 0.28 | -0.243 | 2.58 | -0.86 | 0.280
S4 | 1.26 | 0.25 | -0.250 | 2.01 | -0.62 | 0.255
S5 | 1.31 | 0.27 | -0.416 | 2.56 | 2.39 | 0.261
Table 1: Parameters used to sample galaxies and their shapes in line with the
expectations for the KiDS-1000 gold sample for each tomographic bin (from S1
to S5). $\overline{n}_{\mathrm{eff}}$ is the mean galaxy number density for a
given tomographic bin, $\overline{\sigma}_{\epsilon}/\sqrt{2}$ is the mean
per-component shape dispersion, $a_{n_{\mathrm{eff}}}$ and
$b_{n_{\mathrm{eff}}}$ are the slope and y-intercept, respectively, for the
linear interpolation of the galaxy density as a function of the root-mean-
square of the background noise in the KiDS catalogue, $\sigma_{\mathrm{rms}}$,
according to Eq. (28), while $a_{\sigma_{\epsilon}}$ and
$b_{\sigma_{\epsilon}}$ are the parameters to linearly interpolate
$\sigma_{\epsilon}$ from $\sigma_{\mathrm{rms}}$ according to Eq. (29).

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 SALMO131313Speedy Acquisition for Lensing and Matter
Observables; https://github.com/Linc-tw/salmo framework (Joachimi et al.,
2021) for simulating galaxy catalogues, we Poisson-sample galaxies within each
HEALPix pixel, $m$, on the sky using the following expectation for the galaxy
counts

$\langle
N_{m}^{(i)(p)}\rangle(\mathbf{\Theta})=w_{m}(\Omega_{\mathrm{survey}})\,\Big{[}1+b^{(i)}\delta^{(i)}(\theta_{m};\mathbf{\Theta})\Big{]}\,P_{m}(p|i)\,n_{\mathrm{gal},m}^{(p)}\,A_{\mathrm{pix},m}\,,$
(20)

where $w_{m}(\Omega_{\mathrm{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,
$\Omega_{\mathrm{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\forall 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_{\mathrm{ph}})$),
$n_{\mathrm{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.
3.8), and $A_{\mathrm{pix},m}$ is the area of pixel $m$. As
$n_{\mathrm{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. 3.8).

Taking $w_{m}(\Omega_{\mathrm{survey}})$, $P_{m}(p|i)$ and
$n_{\mathrm{gal},m}^{(p)}$ from the photometric survey’s measurements, we can
Poisson sample $\smash{N_{m}^{(i)(p)}(\mathbf{\Theta})}$ given the expectation
value shown in Eq. (20) from which we can get the expected number of galaxy
counts within each pixel of a tomographic bin $p$ by taking

$N_{m}^{(p)}(\mathbf{\Theta})=\sum_{i=1}^{N_{\mathrm{shells}}}N_{m}^{(i)(p)}(\mathbf{\Theta})\,.$
(21)

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 $\smash{N_{m}^{(i)(p)}(\mathbf{\Theta})}$. 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_{\mathrm{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.

Figure 5: The distributions for each tomographic bin (S1 to S5) from which the
shear bias parameters from KiDS-1000 (Giblin et al., 2021) shown in Eq. (27)
are sampled. The first panel from the left shows the multiplicative shear
bias, $M^{(p)}$. The second and third panel show the real and imaginary part
of the additive shear bias, $c_{1}^{(p)}$ and $c_{2}^{(p)}$, respectively. The
fourth and fifth panel show the real and imaginary part of the amplitude of
the shear bias due to variations in the point-spread function,
$\alpha_{1}^{(p)}$ and $\alpha_{2}^{(p)}$, respectively.

###  3.7 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, $\epsilon$, and any galaxy ellipticity in this work as follows

$\epsilon=\epsilon_{1}+\mathrm{i}\,\epsilon_{2}\,,$ (22)

where $\mathrm{i}\equiv\sqrt{-1}$. The intrinsic $\epsilon_{1}$ and
$\epsilon_{2}$ of each galaxy are then sampled within SALMO as two independent
normal random variates with zero mean and a variance equal to
$\sigma_{\epsilon,m}^{(p)2}$. As with the galaxy density,
$n_{\mathrm{gal},m}^{(p)}$, the shape dispersion, $\sigma_{\epsilon,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. 3.8 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 (Seitz & Schneider,
1996; Hu, 2000),

$\epsilon_{\mathrm{lensed}}(\mathbf{\Theta})=\frac{\epsilon_{\mathrm{int}}+g(\mathbf{\Theta})}{1+g^{*}(\mathbf{\Theta})\epsilon_{\mathrm{int}}}\,,$
(23)

where $\epsilon_{\mathrm{lensed}}(\mathbf{\Theta})$ is the lensed galaxy
ellipticity, $\epsilon_{\mathrm{int}}$ is the intrinsic galaxy ellipticity and
$g(\mathbf{\Theta})$ is the reduced shear ($g\in\mathbb{C}$) which is given by

$g(\mathbf{\Theta})=\frac{\gamma(\mathbf{\Theta})}{1-\kappa(\mathbf{\Theta})}\,,$
(24)

where $\gamma(\mathbf{\Theta})$ is the shear factor ($\gamma\in\mathbb{C}$).
The use of Eq. (23) is an extension to the KiDS-1000 analysis which assumed
the reduced shear approximation, i.e.
$\epsilon_{\mathrm{lensed}}\approx\epsilon_{\mathrm{int}}+g$. In any case,
this should cause a negligible change in the measured shear considering the
sensitivity of KiDS-1000 (Joachimi et al., 2021). To calculate the shear
field, it is useful to decompose the spin-0 convergence field into spherical
harmonics as follows

$\kappa_{\ell
m}^{(i)}(\mathbf{\Theta})=\int\mathbf{d}^{2}\bm{\theta}\,\kappa^{(i)}(\bm{\theta};\mathbf{\Theta})\,_{0}\\!\\!\operatorname{Y}^{*}_{\ell
m}(\bm{\theta})\,,$ (25)

where $\kappa_{\ell m}^{(i)}(\mathbf{\Theta})$ 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, $\gamma_{\ell
m}^{(i)}$, given by (Peebles, 1973),

$\gamma_{\ell
m}^{(i)}(\mathbf{\Theta})=-\Bigg{(}\frac{(\ell+2)(\ell-1)}{\ell(\ell+1)}\Bigg{)}^{1/2}\,\kappa_{\ell
m}^{(i)}(\mathbf{\Theta})\,.$ (26)

The $\gamma_{\ell m}^{(i)}(\mathbf{\Theta})$ coefficients can then be used to
calculate the spin-2 discrete shear field,
$\gamma^{(i)}(\theta_{m},\mathbf{\Theta})$. Having defined the effective
$\kappa^{(i)}(\theta_{m},\mathbf{\Theta})$ in Sect. 3.4 and 3.5, while also
defining $\gamma^{(i)}(\theta_{m},\mathbf{\Theta})$, we can assign a value of
$\kappa$ and $\gamma=\gamma_{1}+\mathrm{i}\,\gamma_{2}$ to each galaxy sampled
according to the procedure described in Sect. 3.6. This is done by taking the
value of $\kappa^{(i)}(\theta_{m},\mathbf{\Theta})$ and
$\gamma^{(i)}(\theta_{m},\mathbf{\Theta})$ within the shell $i$ and pixel $m$
in which a given galaxy is located. Between neighbouring pixels, the values of
$\kappa^{(i)}(\theta_{m},\mathbf{\Theta})$ and
$\gamma^{(i)}(\theta_{m},\mathbf{\Theta})$ are linearly interpolated. This
choice may not be optimal (Tessore et al., 2023), but we find that at the
precision level of KiDS-1000, it does not lead to any significant biases (see
Appendix B). We then use Eq. (24) to calculate the associated reduced shear,
$g(\mathbf{\Theta})$, as well as the lensed ellipticity,
$\epsilon_{\mathrm{lensed}}(\mathbf{\Theta})$, by combining $g$ with the
intrinsic ellipticity using Eq. (23).

However, this lensed ellipticity,
$\epsilon_{\mathrm{lensed}}(\mathbf{\Theta})$, 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. Mandelbaum 2018 for a review). To account for this, it is
common to parametrically estimate any residual systematics which may be
affecting the shape measurements (Hildebrandt et al., 2017, 2021; Zuntz et
al., 2018; Giblin et al., 2021). Within the forward simulations, we use the
first-order parametric expansion of the observed galaxy shapes given by
Heymans et al. (2006) as follows

$\epsilon_{\mathrm{obs},i\in
m}^{(p)}=\Big{(}1+M^{(p)}\Big{)}\epsilon_{\mathrm{lensed},i}+\alpha^{(p)}\epsilon_{\mathrm{PSF},m}+\beta^{(p)}\delta\epsilon_{\mathrm{PSF},m}+c^{(p)}\,,$
(27)

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$, $\alpha^{(p)}$ gives the fraction of the PSF ellipticity
which remains in the shear estimator, $\epsilon_{\mathrm{PSF}}(m)$ is the
local PSF measured within pixel $m$ (see Figure 6 for a map of
$\epsilon_{\mathrm{PSF}}$ in KiDS-1000), $\beta^{(p)}$ gives the amplitude of
the shear bias due to residuals which are not taken into account by the PSF
model, $\delta\epsilon_{\mathrm{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 (Kitching et al., 2019). In
accordance with previous findings (Giblin et al., 2021), we take
$\beta^{(p)}=0\,\forall p$. The amplitudes of the remaining systematics as
calibrated by Giblin et al. (2021) are shown in Fig. 5. 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 (Giblin et al., 2021; Joachimi et al.,
2021; Asgari et al., 2021; Heymans et al., 2021). Additionally, the measured
values of $M^{(p)}$, $\alpha^{(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)}$, $\alpha^{(p)}$ and $c^{(p)}$
for each realisation of the forward simulations from a normal distribution as
shown in Fig. 5. 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, $\mathbf{\Theta}$, 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.

Figure 6: Cartesian spatial map ($N_{\mathrm{side}}=1024$) of the observed
magnitude of the point-spread function ellipticities,
$|\epsilon_{\mathrm{PSF}}|=\sqrt{\epsilon_{\mathrm{PSF},1}^{2}+\epsilon_{\mathrm{PSF},2}^{2}}$,
throughout the KiDS-1000 North field in the upper panel and the KiDS-1000
South field in the lower panel. $\epsilon_{\mathrm{PSF}}$ is added to the
lensed galaxy shapes in the forward simulations within KiDS-SBI in accordance
with Eq. (27).

Figure 7: Cartesian spatial map ($N_{\mathrm{side}}=1024$) of root-mean-square
of the observed background noise, $\sigma_{\mathrm{rms}}$, throughout the
KiDS-1000 North field in the upper panel and the KiDS-1000 South field in the
lower panel. Figure 8: Spatial map ($N_{\mathrm{side}}=4096$) in a Cartesian
projection of root-mean-square of the observed background noise,
$\sigma_{\mathrm{rms}}$, for a $5^{\circ}\times 5^{\circ}$ patch of the
KiDS-1000 North field at a right ascension (RA) of 180∘ and a declination
(DEC) of 0∘.

###  3.8 Variable depth

As alluded to in Sect. 3.6 and 3.7, 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_{\mathrm{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 ${\sim}10\%$ in the standard deviation of cosmic
shear observables (Heydenreich et al., 2020; Joachimi et al., 2021).

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, $\sigma_{\epsilon,m}^{(p)}$, and further exacerbate
the effect of observational depth variability on the cosmic shear signal.

Figure 9: The dependence of the per-component Gaussian shape dispersion,
$\sigma_{\epsilon}/\sqrt{2}$, (left panel) and the galaxy density,
$n_{\mathrm{eff}}$, (right panel) on the root-mean square of the background
noise, $\sigma_{\mathrm{rms}}$ in the KiDS-1000 DR4 data. For both panels, the
data points represent the mean $\sigma_{\epsilon}$ or $n_{\mathrm{eff}}$ of
ten equi-populated bins in $\sigma_{\mathrm{rms}}$ with their boundaries in
$\\{1.70,\,2.33,\,2.57,\,2.80,\,3.04,\,3.28,\,3.52,\,3.76,\,4.00,\,4.23,\,12.96\\}$.
The solid line shows the linear fit to the aforementioned data points of their
respective tomographic bin according to Eqs. (28) and (29). The parameters
obtained from this fit are given in Table 1. The dotted horizontal lines show
the mean values of $\sigma_{\epsilon}$ and $n_{\mathrm{eff}}$ calculated from
the galaxy samples with variable depth per tomographic bin, while the dashed
horizontal lines show the values of $\sigma_{\epsilon}$ and $n_{\mathrm{eff}}$
for the respective galaxy samples without any spatial variations in the
observational depth. Both of these lines agree exceptionally well by
construction, so that for some source bins there is not any observable
difference between them.

To take all these effects into consideration within the forward simulations,
we calibrate the spatial variability of $n_{\mathrm{eff},m}^{(p)}$ and
$\sigma_{\epsilon,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_{\mathrm{eff},m}^{(p)}$ and
$\sigma_{\epsilon,m}^{(p)}$: the root-mean square of the background noise,
$\sigma_{\mathrm{rms}}$. We found this quantity to be a good indicator of
variable depth as it correlates well with the measured
$n_{\mathrm{eff},m}^{(p)}$ and $\sigma_{\epsilon,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 Joachimi et
al. (2021). 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, $\sigma_{\mathrm{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. 7 and 8, a map of $\sigma_{\mathrm{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. 9 shows that
$n_{\mathrm{eff}}^{(p)}$ and $\sigma_{\epsilon}^{(p)}$ vary linearly with
$\sigma_{\mathrm{rms}}$ for KiDS-1000 data. Thanks to this fact, we do not
need to model the variation of $n_{\mathrm{eff},m}^{(p)}$ and
$\sigma_{\epsilon,m}^{(p)}$ with the HEALPix pixel $m$ with individual maps
for each tomographic bin $p$. Instead, using the linear relations shown in
Fig. 9, we may assign a value for $n_{\mathrm{eff},m}^{(p)}$ and
$\sigma_{\epsilon,m}^{(p)}$ for each $\sigma_{\mathrm{rms},m}$ (i.e. the mean
value of level within each pixel $m$ contained within the footprint observed
by KiDS-1000) given that

$\displaystyle n_{\mathrm{eff},m}^{(p)}=$
$\displaystyle\;a_{n_{\mathrm{eff}}}^{(p)}\sigma_{\mathrm{rms},m}+b_{n_{\mathrm{eff}}}^{(p)}\,,$
(28) $\displaystyle\sigma_{\epsilon,m}^{(p)}=$
$\displaystyle\;a_{\sigma_{\epsilon}}^{(p)}\sigma_{\mathrm{rms},m}+b_{\sigma_{\epsilon}}^{(p)}\,,$
(29)

where $a_{n_{\mathrm{eff}}}^{(p)}$ and $a_{\sigma_{\epsilon}}^{(p)}$ are the
slopes of the linear fits shown in Fig. 9 of galaxy density and galaxy shape
dispersion for a given tomographic bin, respectively, and
$b_{n_{\mathrm{eff}}}^{(p)}$ and $b_{\sigma_{\epsilon}}^{(p)}$ are the
associated y-intercepts. The values of these parameters are shown in Table 1.
This drastically helps the performance of the forward simulations, as this
allows to just use a single map of $\sigma_{\mathrm{rms}}$ to model the
variable depth rather than requiring $2N_{\mathrm{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_{\mathrm{eff}}^{(p)}$ and $\sigma_{\epsilon}^{(p)}$ on
$\sigma_{\mathrm{rms}}$ to define an interpolation.

These inhomogeneities in the galaxy selection across the survey footprint also
change the redshift distribution of galaxies. As $\sigma_{\mathrm{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. 10. The redshift distributions
shown in Fig. 10 have been calculated using the same approach as in KiDS-1000
where the photometric redshifts are calibrated from spectroscopic samples
using self-organising maps (Wright et al., 2019). 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\sigma_{\mathrm{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
$\sigma_{\mathrm{rms}}$ in a given pixel, $m$, on the survey footprint, we
sample galaxy redshifts in that pixel from the associated redshift
distribution from Fig. 10. 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 (Heydenreich et al., 2020; Baleato Lizancos &
White, 2023).

Figure 10: In the upper panels, the normalised redshift distributions, $P(z)$,
for each tomographic bin (S1 to S5). The redshift distribution from the entire
KiDS-1000 DR4 galaxy sample, $P_{\mathrm{total}}(z)$, is shown in black, while
the other ten redshift distributions are derived from 10 equi-populated
subsamples of DR4 based on their observational depth (i.e. the mean value of
the root-mean-square of the background noise,
$\overline{\sigma}_{\mathrm{rms}}$) which is shown with its respective colour.
The lower panels show the associated residual change in the redshift
distributions with respect to $P_{\mathrm{total}}(z)$ per unit redshift. It is
apparent that variable depth mostly affects the source distributions at high
redshifts, while the effect tends to decrease the mean of the redshift
distribution with increasing $\sigma_{\mathrm{rms}}$.

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

##  4 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\forall i$, while the weights used in the measurements from
the KiDS-1000 gold sample are given by the lensfit weights (see details in
Sect. 2). 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_{\mathrm{gal}}$, and
intrinsic shape dispersions, $\sigma_{\epsilon}$, where the corrections for
the lensfit weight and the multiplicative shear bias have been applied, as
defined in Eqs. (18) and (19).

###  4.1 Shape measurements

From the KiDS-1000 gold sample data as well as the forward simulated
catalogues, we measure the observed galaxy ellipticity,
$\epsilon_{\mathrm{obs}}$, for each galaxy per tomographic bin. To this end,
we apply the same corrections to the shear measurements as described in Giblin
et al. (2021). The observed shape measurements are corrected as follows

$\epsilon_{\mathrm{obs},i}^{\mathrm{corr}\,(p)}(\mathbf{\Theta})=\frac{1}{1+\overline{M}^{(p)}}\frac{\sum_{i\in
p}{w_{i}(\epsilon_{\mathrm{obs},i}(\mathbf{\Theta})-\langle\epsilon_{\mathrm{obs}}(\mathbf{\Theta})\rangle)}}{\sum_{i\in
p}{w_{i}}}\,,$ (30)

where $w_{i}$ is the galaxy lensfit weight for galaxy $i$ in the KiDS-1000
gold sample (for the case of the simulations, $w_{i}=1\forall i$),
$\langle\epsilon_{\mathrm{obs}}(\mathbf{\Theta})\rangle$ compensates for any
additive biases (Asgari et al., 2019), i.e. the $c^{(p)}$ term and the mean of
$\alpha^{(p)}$ in Eq. (27), and it is defined as the mean observed shape of
all galaxies within tomographic bin $p$, and $\overline{M}^{(p)}$ is the mean
multiplicative shear bias measured for each tomographic bin $p$ as shown in
Fig. 5. Note that $\overline{M}^{(p)}$ is not necessarily the same value as
the value of $M^{(p)}$ applied in Eq. (27). For each instance of the
simulations, a different $M^{(p)}$ is drawn from a Gaussian probability
distribution with mean $\overline{M}^{(p)}$ and the standard deviations shown
in Fig. 5. Regardless of the value of $M^{(p)}$ drawn, the shear measurements
of each simulation are corrected by the same mean $\overline{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.

###  4.2 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,
$\mathbf{\Theta}$, 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 (Hall & Taylor, 2022; Piras et al., 2023). 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 (Loureiro et al., 2022; Asgari et al., 2021) 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,
$\epsilon_{\mathrm{obs}}^{\mathrm{corr}\,(p)}(\theta_{m};\mathbf{\Theta})$,
into its curl-free E-modes and its divergence-free B-modes as follows

$\epsilon_{\mathrm{obs}}^{\mathrm{corr}\,(p)}(\theta;\mathbf{\Theta},\Omega_{\mathrm{survey}})=\sum_{\ell=0}^{\ell_{\mathrm{max}}}\sum_{m=-\ell}^{\ell}(\tilde{E}_{\ell
m}^{(p)}(\mathbf{\Theta})+\mathrm{i}\tilde{B}_{\ell
m}^{(p)}(\mathbf{\Theta}))\,_{2}\\!\\!\operatorname{Y}_{\ell m}(\theta)\,,$
(31)

where $\,{}_{\pm 2}\\!\\!\operatorname{Y}^{\mathstrut}_{\ell m}(\theta)$ are
the spin-2 spherical harmonic functions which define an orthonormal basis for
the observed shear field, such that

$\int\mathrm{d}^{2}\theta\,_{\pm 2}\\!\\!\operatorname{Y}^{\mathstrut}_{\ell
m}(\bm{\theta})\,_{\pm
2}\\!\\!\operatorname{Y}^{*}_{\ell^{\prime}m^{\prime}}(\bm{\theta})=\delta^{\mathrm{K}}_{\ell\ell^{\prime}}\delta^{\mathrm{K}}_{mm^{\prime}}\,,$
(32)

while $\tilde{E}_{\ell m}^{(p)}(\mathbf{\Theta})$ and $\tilde{B}_{\ell
m}^{(p)}(\mathbf{\Theta})$ are spherical harmonic coefficients of the curl-
free and the divergence-free observed shear fields within each tomographic
bin, respectively, and $\delta^{\mathrm{K}}$ is the Kronecker delta. These
coefficients are defined as follows:

$\displaystyle\tilde{E}_{\ell
m}^{(p)}(\mathbf{\Theta})\;=\;\frac{1}{2}\;\int\mathrm{d}^{2}$
$\displaystyle\theta\,[\epsilon_{\mathrm{obs}}^{\mathrm{corr}\,(\mathit{p})}(\bm{\theta};\mathbf{\Theta})\,_{2}\\!\\!\operatorname{Y}^{*}_{\ell
m}(\bm{\theta})$
$\displaystyle+\epsilon_{\mathrm{obs}}^{\mathrm{corr}\,(\mathit{p})*}(\bm{\theta};\mathbf{\Theta})\,_{-2}\\!\\!\operatorname{Y}^{*}_{\ell
m}(\bm{\theta})]\,,$ (33) $\displaystyle\tilde{B}_{\ell
m}^{(p)}(\mathbf{\Theta})=\frac{-\mathrm{i}}{2}\int\mathrm{d}^{2}{\theta}$
$\displaystyle\,[\epsilon_{\mathrm{obs}}^{\mathrm{corr}\,(\mathit{p})}(\bm{\theta};\mathbf{\Theta})\,_{2}\\!\\!\operatorname{Y}^{*}_{\ell
m}(\bm{\theta})$
$\displaystyle-\epsilon_{\mathrm{obs}}^{\mathrm{corr}\,(\mathit{p})*}(\theta;\mathbf{\Theta})\,_{-2}\\!\\!\operatorname{Y}^{*}_{\ell
m}(\bm{\theta})]\,.$ (34)

From these coefficients, we then calculate the pseudo-Cls,
$\tilde{C}^{(pq)}_{\epsilon\epsilon}(\ell;\mathbf{\Theta})$, as follows

$\tilde{C}^{(pq)}_{\epsilon\epsilon,\mu}(\ell;\mathbf{\Theta})=\begin{pmatrix}\tilde{C}_{\epsilon\epsilon}^{EE(pq)}(\ell;\mathbf{\Theta})\\\
\tilde{C}_{\epsilon\epsilon}^{EB(pq)}(\ell;\mathbf{\Theta})\\\
\tilde{C}_{\epsilon\epsilon}^{BB(pq)}(\ell;\mathbf{\Theta})\end{pmatrix}=\frac{1}{2\ell+1}\sum_{m=-\ell}^{\ell}\begin{pmatrix}\tilde{E}^{(p)}_{\ell
m}\tilde{E}^{(q)*}_{\ell m}(\mathbf{\Theta})\\\ \tilde{E}^{(p)}_{\ell
m}\tilde{B}^{(q)*}_{\ell m}(\mathbf{\Theta})\\\ \tilde{B}^{(p)}_{\ell
m}\tilde{B}^{(q)*}_{\ell m}(\mathbf{\Theta})\end{pmatrix}\,,$ (35)

where $\mu\in\\{1,2,3\\}$ such that $1$ stands for the EE component, $2$ for
the EB component and $3$ for the BB component. Note that
$\tilde{C}_{\epsilon\epsilon}^{EB}(\ell;\mathbf{\Theta})=\tilde{C}_{\epsilon\epsilon}^{BE}(\ell;\mathbf{\Theta})$,
since $\tilde{E}_{\ell m}\tilde{B}_{\ell
m}^{*}(\mathbf{\Theta})=\tilde{B}_{\ell m}\tilde{E}_{\ell
m}^{*}(\mathbf{\Theta})$.

Figure 11: Measured pseudo-Cls,
$\tilde{C}^{(pq)}_{\epsilon\epsilon,L}(\bm{\Theta})$, for all combinations for
five tomographic bins (S1 to S5). The green triangles show the pseudo-Cls from
the KiDS-1000 gold sample data. The blue points show the same for a noise-free
realisation single realisation of the anisotropic systematics model in KiDS-
SBI (see Sect. 5) at the maximum-a-posteriori, MAP (see Table 5), while the
magenta crosses show the measurements for noise-free realisation single
realisation of the standard isotropic systematics model in KiDS-SBI (see Sect.
7.2) at the MAP for this model (see Table 5). The uncertainties on the
measurements are derived from the covariance matrix described in Sect. 5.3.
The solid black line shows the pseudo-Cls, $\tilde{C}_{\epsilon\epsilon}$, as
derived from theory (see Appendix D for details on this) at the MAP from KiDS-
SBI based on the anisotropic systematics model. The horizontal dashed line
marks the line of zero signal.

Going forward, similar to other previous analyses (Hikage et al., 2019;
Loureiro et al., 2022; Tröster et al., 2022), we only take into consideration
$\tilde{C}^{EE(pq)}_{\epsilon\epsilon}(\ell;\mathbf{\Theta})$, so any
cosmological signal which may have ended up in
$\tilde{C}^{EB(pq)}_{\epsilon\epsilon}(\ell;\mathbf{\Theta})$ or
$\tilde{C}^{BB(pq)}_{\epsilon\epsilon}(\ell;\mathbf{\Theta})$ because of EE to
EB, or EE to BB mode mixing will be lost.
$\tilde{C}^{EE(pq)}_{\epsilon\epsilon}(\ell;\mathbf{\Theta})$ is expected to
still contain all of the cosmological cosmic shear signal because the shear
field resulting from a scalar gravitational field in $\Lambda$CDM is predicted
to be a curl-free field (Bartelmann & Schneider, 2001; Hikage et al., 2011;
Kilbinger, 2015).

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

$\tilde{C}^{(pq)}_{\epsilon\epsilon}(\ell;\mathbf{\Theta})=\tilde{C}^{EE(pq)}_{\epsilon\epsilon}(\ell;\mathbf{\Theta})-\delta^{\mathrm{K}}_{pq}\langle\tilde{C}^{EE(pq)}_{\mathrm{noise}}(\ell)\rangle\,,$
(36)

where $\tilde{C}^{(pq)}_{\epsilon\epsilon}(\ell;\mathbf{\Theta})$ is the
E-mode pseudo angular power spectrum for cosmic shear with the shape noise
bias subtracted, while
$\langle\tilde{C}^{EE(pq)}_{\mathrm{noise}}(\ell)\rangle$ is the mean of the
curl-free angular power spectrum of the shape noise. The latter is estimated
as follows (Becker et al., 2016; Hikage et al., 2019; Nicola et al., 2021;
Loureiro et al., 2022),

$\tilde{C}^{EE(pp)}_{\mathrm{noise}}(\ell)=\frac{1}{2\ell+1}\sum_{m=-\ell}^{\ell}\tilde{E}^{\mathrm{rand}\,(p)}_{\ell
m}(\mathbf{\Theta})\tilde{E}^{\mathrm{rand}\,(p)*}_{\ell
m}(\mathbf{\Theta})\,,$ (37)

where $\tilde{E}^{\mathrm{rand}\,(p)}_{\ell m}$ is the curl-free spherical
harmonic coefficient of the randomly rotated shear values,
$\epsilon^{\mathrm{rand}\,(p)}_{\mathrm{obs}}(\theta_{m};\mathbf{\Theta})$. In
turn, we define this field from the following galaxy shear values as follows:

$\displaystyle\epsilon^{\mathrm{rand}}_{\mathrm{obs},i}(\mathbf{\Theta})=\;\epsilon^{\mathrm{rand}\,}_{\mathrm{obs},i,1}(\mathbf{\Theta})+\mathrm{i}\epsilon^{\mathrm{rand}\,}_{\mathrm{obs},i,2}(\mathbf{\Theta})\,,$
(38)
$\displaystyle\epsilon^{\mathrm{rand}\,}_{\mathrm{obs},i,1}=\;\epsilon^{\mathrm{corr}\,}_{\mathrm{obs},i,1}\,\mathrm{cos}(\theta_{\mathrm{rand},i})-\epsilon^{\mathrm{corr}\,}_{\mathrm{obs},i,2}\,\mathrm{sin}(\theta_{\mathrm{rand},i})\,,$
(39)
$\displaystyle\epsilon^{\mathrm{rand}\,}_{\mathrm{obs},i,2}=\;\epsilon^{\mathrm{corr}\,}_{\mathrm{obs},i,2}\,\mathrm{cos}(\theta_{\mathrm{rand},i})+\epsilon^{\mathrm{corr}\,}_{\mathrm{obs},i,1}\,\mathrm{sin}(\theta_{\mathrm{rand},i})\,,$
(40)

where $\theta_{\mathrm{rand},i}$ is a randomly drawn angle for each galaxy $i$
from a uniform distribution where $\theta_{\mathrm{rand},i}\in[0,2\pi)$. To
calculate the mean of the curl-free angular power spectrum of the shape noise,
$\langle\tilde{C}^{EE(pq)}_{\mathrm{noise}}(\ell)\rangle$, we take mean of all
modes in $\tilde{C}^{EE(pq)}_{\mathrm{noise}}(\ell)$ as follows

$\langle\tilde{C}^{EE(pq)}_{\mathrm{noise}}(\ell)\rangle=\frac{\sum_{\ell=\ell_{\mathrm{min}}}^{\ell_{\mathrm{max}}}(2\ell+1)\,\tilde{C}^{EE(pq)}_{\mathrm{noise}}(\ell)}{\sum_{\ell=\ell_{\mathrm{min}}}^{\ell_{\mathrm{max}}}(2\ell+1)}\,,$
(41)

where the mean is weighted by a factor of $2\ell+1$ to account for the fact
that each angular scale, $\ell$, contains $2\ell+1$ modes, $m$. We find
setting $\ell_{\mathrm{min}}{\sim}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, $\langle\tilde{C}^{EE(pq)}_{\mathrm{noise}}(\ell)\rangle$, agrees
within $<0.5\%$ with the underlying shape noise signal computed directly from
the known intrinsic galaxy shapes. Additionally, we find that
$\langle\tilde{C}^{EE(pq)}_{\mathrm{noise}}(\ell)\rangle$ is self-consistent
with $\langle\tilde{C}^{BB(pq)}_{\mathrm{noise}}(\ell)\rangle$ within
$<0.01\%$. The introduction of a weighted mean given by Eq. (41) 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
$\langle\tilde{C}^{EE(pq)}_{\mathrm{noise}}(\ell)\rangle$ by taking the
average over all the instances while weighting each angular scale, $\ell$,
equally (see e.g. Loureiro et al. 2022). When weighting each $\ell$ equally,
the low-$\ell$ scales have an equally strong pull on the estimator as the
high-$\ell$ scales which increases the statistical uncertainty on
$\langle\tilde{C}^{EE(pq)}_{\mathrm{noise}}(\ell)\rangle$ 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 ${\sim}$$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,
$\tilde{C}^{(pq)}_{\epsilon\epsilon}(\ell;\mathbf{\Theta})$, into 8 log-spaced
bins between $\ell=76$ and $\ell=1500$ in line with Loureiro et al. (2022).
For this, we choose the pseudo-Cl binning scheme described in Brown et al.
(2005) given by

$\tilde{C}^{(pq)}_{\epsilon\epsilon,L}(\mathbf{\Theta})=\frac{1}{2\pi}\sum_{\ell=\ell_{L}}^{\ell_{L+1}}\frac{\ell(\ell+1)}{(\ell_{L+1}-\ell_{L})}\tilde{C}^{(pq)}_{\epsilon\epsilon}(\ell;\mathbf{\Theta}),$
(42)

where $\ell_{L}$ and $\ell_{L+1}$ are the lower and upper limits of the
$L^{\mathrm{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 D) 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. 11 for a single run of the forward simulations. The
calculation of the data vector calculated from the theory prediction shown in
Fig. 11 is described in Appendix D. As can be seen from Figure 11, 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
B). 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. 11, for the five tomographic bins in KiDS-1000, we obtain
$8\times N_{\mathrm{tomo}}(N_{\mathrm{tomo}}+1)/2$ data points in
$\tilde{C}^{(pq)}_{\epsilon\epsilon,L}(\mathbf{\Theta})$, i.e. a
120-dimensional data vector. The analysis choices described within Sect 3 and
4 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 $\mathbf{\Theta}$:
5 cosmological parameters ($\sigma_{8}$, the root-mean-square matter
fluctuation over 8 Mpc$/h_{0}$; $\omega_{b}$, baryonic matter density;
$\omega_{c}$, cold dark matter density; $n_{s}$, the scalar spectral index of
the primordial density fluctuation power spectrum, $P_{\mathcal{R}}$; and
$h_{0}$141414Note 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_{\mathrm{bary}}$, the baryonic feedback amplitude within the non-linear
three-dimensional matter power spectrum; $A_{\mathrm{IA}}$, the galaxy
intrinsic alignment amplitude of the NLA model; and five correlated
$\delta_{z}$ parameters which define the shift in the mean of the source
redshift distribution, $P(z|z_{\mathrm{ph}})$, of each tomographic bin. All
other parameters on which the simulation depends are fixed from run to run
(e.g. $\omega_{k}=0$, the equation of state of dark energy, $w=-1$, $\sum
m_{\nu}=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. 3.7. 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.

##  5 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 (pyDELFI151515https://github.com/justinalsing/pydelfi,
Alsing et al. 2019, or DELFI henceforth to denote the method and software
package) using the same analysis pipeline outlined in L23.

To implement this, we create a bespoke interface that allows DELFI to be the
sampler for CosmoSIS (Zuntz et al., 2015). Initially, the sampler chooses a
fixed number of points within the considered parameter space based on the
hypercube-generating algorithm from
PyDOE161616https://github.com/tisimst/pyDOE with some minor modifications (see
L23 for more details). Subsequently, as is outlined in Fig. 12, the measured
pseudo-Cls for each evaluation of the parameters on the hypercube are
compressed further using score compression (see Sect. 5.3). 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 DELFI (see Sect. 5.4).

As the ensemble of neural density estimators evaluates the sampling
distribution throughout the hypercube, through the use of active learning,
DELFI finds areas in parameter space where the likelihood is undersampled to
select new parameter vectors (see Alsing et al. 2019 for details). As found in
L23, this can halve the number of simulations needed to accurately train DELFI
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 DELFI.

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. 3 and 4.
It is designed to be consistent with previous analyses of KiDS-1000,
specifically, with Loureiro et al. (2022), while adding additional realism.
Both Loureiro et al. (2022) and this work assume a flat $\Lambda$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 Loureiro et al. (2022) in that it
also includes variable depth in the signal modelling (as described in Sect.
3.8), 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. 3.7). 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 Loureiro et al. (2022), these effects are only
considered separately as they pertain to the signal or the uncertainty, not
both at once.

Cosmological parameters, $\bm{\Theta}$ Forward simulations Mock data: Pseudo-
Cl, $\tilde{C}_{\epsilon\epsilon}^{(pq)}(\ell)$ Score compression Train Neural
Density Estimators Effective likelihood Observed data Priors Posterior Active
learning Figure 12: Flowchart describing the structure of the simulation-based
inference pipeline. The dark blue rounded boxes represent the inputs and
outputs which are given to the simulation-based inference pipeline. The grey
rectangular boxes show steps in the inference pipeline.

###  5.1 Pipeline performance

A single such evaluation of the forward simulations from cosmological
parameters to the shear pseudo-Cls as described in Sect. 3 runs within
${\sim}$20 minutes on a single core (with $N_{\mathrm{shells}}=19$,
$N_{\mathrm{tomo}}=5$ and $N_{\mathrm{side}}=1024$). We choose this spatial
resolution $N_{\mathrm{side}}$, because it implies that the measured two-point
statistics should be representative of the input angular power spectra up to
an $\ell_{\mathrm{max}}{\sim}2\times N_{\mathrm{side}}$ (Leistedt et al.,
2013; Alonso et al., 2019). In addition, for KiDS-1000’s galaxy number
density, the pixel size at an $N_{\mathrm{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 D).

In comparison, as can be seen in the bar chart in Fig. 13, log-normal random
field simulations similar to the ones presented in Joachimi et al. (2021),
take ${\sim}$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 Joachimi et al. (2021) are used exclusively for the
characterisation of a numerical covariance matrix. Hence, only a single
evaluation of the underlying power spectrum from CAMB is necessary for every
simulation. At the same time, the Joachimi et al. (2021) simulations have a
lower resolution of shells along the line-of-sight, so the runtime of FLASK
(Xavier et al., 2016) in the original setting is lower than the runtime
indicated in Fig. 13. We also note that the runtime of TreeCorr (Jarvis et
al., 2004) could be drastically reduced from the single-core runtime shown in
Fig. 13 by increasing the number of cores and/or by pixelising the galaxy
shears prior to calculating the spatial two-point correlation functions,
$\xi(\theta)$.

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 GLASS, the choice of using a quick-to-compute
summary statistic like angular power spectra, $\tilde{C}(\ell)$, 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. 3.6, 3.7 and 3.8, while still being able to compute the
${\sim}$$10^{4}$ forward simulations needed to adequately characterise the
effective likelihood L23.

Figure 13: Bar chart comparing the run-time of a single evaluation of KiDS-SBI
(above) versus a single evaluation of a simulation based on Joachimi et al.
(2021) (below), both on a single core ($N_{\mathrm{shells}}=19$,
$N_{\mathrm{tomo}}=5$ and $N_{\mathrm{side}}=1024$). Both suites of
simulations use CAMB (Lewis et al., 2000; Lewis & Challinor, 2002; Howlett et
al., 2012) to compute the three-dimensional matter power spectrum. For the
reference simulations, we use the non-Limber projection built into CAMB with
limber_phi_lmin $=1200$ rather than Levin with $\ell_{\mathrm{max,nL}}=1200$.
We run FLASK (Xavier et al., 2016) rather than GLASS (Tessore et al., 2023) to
compute the underlying matter and convergence fields of each of the 19 shells.
Subsequently, we sample galaxies using SALMO in both cases, and then calculate
the spatial two-point correlation functions, $\xi_{\pm}(\theta)$, rather than
calculating the angular power spectra, $\tilde{C}(\ell)$. To calculate
$\xi_{\pm}(\theta)$ in the reference simulations, we use TreeCorr (Jarvis et
al., 2004).

###  5.2 Parameters and priors

The parameters which are varied in the simulation-based inference in this work
are shown in Table 2. Their priors match those of the previous KiDS-1000
cosmic shear analyses (Asgari et al., 2021; Heymans et al., 2021; Loureiro et
al., 2022; van den Busch et al., 2022; Tröster et al., 2022).

To avoid overinformative priors, most are flat top-hat functions. For the same
reason, the top-hat priors in $S_{8}$, $h_{0}$, $\omega_{\mathrm{b}}$,
$A_{\mathrm{IA}}$ and $A_{\mathrm{bary}}$ are chosen to be wide. The prior on
$n_{\mathrm{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 $\omega_{\mathrm{c}}$ is defined to be consistent with a
range in $\Omega_{\mathrm{m}}\in[0.188,0.408]$, where the limits are given by
the $\pm 5\sigma$ intervals of the marginal constraints from independent
measurements of luminosity distance to Type Ia Supernovae in Scolnic et al.
(2018). See Joachimi et al. (2021) 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, $\bm{\delta}_{z}$, are given by
a multivariate Gaussian, since the $\delta_{z}$ of a given tomographic bin is
not independent of the shifts in the other bins. This is quantified by the
covariance, $\bm{C}_{z}$, as estimated in Hildebrandt et al. (2021).

Parameter | Symbol | Prior type | Prior range | Initial fiducial | Final fiducial
---|---|---|---|---|---
Density fluctuation amp. | $S_{8}$ | Flat | [0.1, 1.3] | 0.760 | 0.739
Hubble constant | $h_{0}$ | Flat | [0.64, 0.82] | 0.767 | 0.653
Cold dark matter density | $\omega_{\mathrm{c}}$ | Flat | [0.051, 0.255] | 0.118 | 0.135
Baryonic matter density | $\omega_{\mathrm{b}}$ | Flat | [0.019, 0.026] | 0.026 | 0.023
Scalar spectral index | $n_{\mathrm{s}}$ | Flat | [0.84, 1.1] | 0.901 | 0.972
Intrinsic alignment amp. | $A_{\mathrm{IA}}$ | Flat | [-6, 6] | 0.264 | 0.603
Baryon feedback amp. | $A_{\mathrm{bary}}$ | Flat | [2, 3.13] | 3.10 | 2.66
Redshift displacement | $\bm{\delta}_{z}$ | Gaussian | $\mathcal{N}(\bm{0},\bm{C}_{z})$ | $\mathbf{0}$ | $\mathbf{0}$
Multiplicative shear bias | $M^{(p)}$ | Gaussian | $\mathcal{N}(\overline{M}^{(p)},\sigma^{(p)}_{M})$ | $\overline{M}^{(p)}$ | $\overline{M}^{(p)}$
Additive shear bias | $c_{1,2}^{(p)}$ | Gaussian | $\mathcal{N}(\overline{c}_{1,2}^{(p)},\sigma^{(p)}_{c_{1,2}})$ | $\overline{c}_{1,2}^{(p)}$ | $\overline{c}_{1,2}^{(p)}$
PSF variation shear bias | $\alpha_{1,2}^{(p)}$ | Gaussian | $\mathcal{N}(\overline{\alpha}_{1,2}^{(p)},\sigma^{(p)}_{\alpha_{1,2}})$ | $\overline{\alpha}_{1,2}^{(p)}$ | $\overline{\alpha}_{1,2}^{(p)}$
Table 2: Table showing the parameters which are varied within the simulation-
based inference pipeline. The prior ranges are selected to be exactly in line
with previous KiDS-1000 analyses (Asgari et al., 2021; Heymans et al., 2021;
Loureiro et al., 2022; van den Busch et al., 2022; Tröster et al., 2022). The
upper five rows show the cosmological parameters of interest, while the three
rows below show the nuisance parameters which quantify systematic biases. The
last three rows show shear bias parameters for each tomographic bin, $p$, as
given in Fig. 5 which are sampled implicitly within each simulation by drawing
from the prior distribution, so any posterior is pre-marginalised over these
parameters. For flat priors, the lower and upper limits of the normalised
rectangular function defines the prior. For the Gaussian prior on
$\bm{\delta}_{z}$, we use a five-dimensional multivariate Gaussian with its
mean at the zero vector and the covariance, $\bm{C}_{z}$, defined by the one
estimated in (Hildebrandt et al., 2021). Note that for simplicity the
$\bm{\delta}_{z}$ are implicitly marginalised throughout this analysis. The
fiducial values given here constitute the fiducial parameter choice for the
score compression described in Sect. 5.3.

###  5.3 Score compression

To improve the computational efficiency of the analysis and to facilitate the
use of DELFI, we compress each of the 120-dimensional data vectors measured
from each forward simulation using score compression as described in Alsing et
al. (2019).

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 (Zablocki & Dodelson, 2016;
Alsing & Wandelt, 2018; Alsing et al., 2018, 2019). This allows us to compress
the data vector, $\bm{d}$, down to a vector, $\bm{t}$, with the dimensionality
of $|\bm{\Theta}|$ as follows

$\bm{t}=\bm{\nabla}\mathcal{L}(\bm{d}|\bm{\Theta}_{*})\,,$ (43)

where $\mathcal{L}$ is the log-likelihood distribution evaluated at the
fiducial cosmology, $\bm{\Theta}_{*}$.

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 2 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 MOPED
(Heavens et al., 2000b) or a linear compression based on Karhunen-Loéve
eigenvalue decomposition (Tegmark et al., 1997).

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 L23 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 2 are not
necessary to characterise the degrees of freedom of the 120-dimensional cosmic
shear pseudo-Cls. As all $\delta_{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
$\bm{\delta}_{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 $\delta_{z}$ parameters.

###  5.4 Density estimation likelihood-free inference (DELFI)

To estimate the effective likelihood for the anisotropic systematics model, we
train neural density estimators (NDEs) using DELFI (Alsing et al., 2019) 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; Papamakarios et
al. 2017). Each MAF is made up of between three and eight Masked Autoencoders
for Density Estimation (MADEs; Germain et al. 2015) each with two hidden
layers of 50 neurons. For a detailed definition and description of these
network architectures, see Alsing et al. (2019).

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 DELFI 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(\bm{t}|\bm{\Theta},\bm{w})$ is well learnt, recursive sampling
continues until the validation loss stops decreasing for 20 training epochs,
where the loss function is defined as,

$-\mathrm{ln}[U(\bm{w})]=-\sum_{i=1}^{N_{\mathrm{sims}}}\mathrm{ln}[P(\bm{t}_{i}|\bm{\Theta}_{i},\bm{w})],$
(44)

where $N_{\mathrm{sims}}$ is the number of forward simulations and the loss
function, $-\mathrm{ln}[U(\bm{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 Alsing et al. 2019 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 L23 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 nautilus (Lange, 2023) and check for
consistency with emcee (Foreman-Mackey et al., 2013).

As advised in L23, to minimise the loss of information through from the
compression step described in Sect. 5.3, 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. 3). 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. (43) 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 (Lin et al.,
2023), we find that this step only marginally alters the inferred posterior.

##  6 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 $\chi^{2}$ goodness-of-fit test. We follow the work of
Gelman et al. (1996) and implement a SBI goodness-of-fit test from simulations
in the ECP-GF171717https://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
(Nelder & Mead, 1965) optimiser. As the parameters of $n_{s}$ and
$a_{\mathrm{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 Lin et al. (2023) 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 $\chi^{2}$ discrepancy
measure as follows

$\chi^{2}(\bm{t}|\bm{\Theta})\coloneqq(\bm{t}_{i}-\mathrm{E}[\bm{t}_{*}|\bm{\Theta}_{*}])^{\mathrm{T}}{(\mathrm{Cov}(\bm{t}_{*}|\bm{\Theta}_{*})})^{-1}(\bm{t}_{i}-\mathrm{E}[\bm{t}_{*}|\bm{\Theta}_{*}])\,,$
(45)

where $\bm{t}_{*}$ is the compressed data vector at the fiducial cosmology,
$\bm{\Theta}_{*}$ the fiducial cosmology and $\bm{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