# KiDS-1000 and DES-Y1 combined: Cosmology from peak count statistics

Joachim Harnois-Déraps1, Sven Heydenreich2, Benjamin Giblin3, Nicolas
Martinet4, Tilman Tröster5, Marika Asgari1,6,7, Pierre Burger8,9,10,Tiago
Castro11,12,13,14, Klaus Dolag15, Catherine Heymans3,16, Hendrik
Hildebrandt16, Benjamin Joachimi17 & Angus H. Wright16
1School of Mathematics, Statistics and Physics, Newcastle University, Herschel
Building, NE1 7RU, Newcastle-upon-Tyne, UK
2Department of Astronomy and Astrophysics, University of California, Santa
Cruz, 1156 High Street, Santa Cruz, CA 95064 USA
3Institute for Astronomy, University of Edinburgh, Blackford Hill, Edinburgh,
EH9 3HJ, UK
4Université d’Aix-Marseille, CNRS, CNES, LAM, Marseille, France
5Institute for Particle Physics and Astrophysics, ETH Zürich, 8092 Zürich,
Switzerland
6E.A Milne Centre, University of Hull, Cottingham Road, Hull, HU6 7RX, UK
7Excellence for Data Science, AI, and Modelling (DAIM), University of Hull,
Cottingham Road, Kingston-upon-Hull, HU6 7RX
8Waterloo Centre for Astrophysics, University of Waterloo, Waterloo, ON N2L
3G1, Canada
9Department of Physics and Astronomy, University of Waterloo, Waterloo, ON N2L
3G1, Canada
10Argelander-Institut für Astronomie, Auf dem Hügel 71, 53121 Bonn, Germany
11INAF - Osservatorio Astronomico di Trieste, via Tiepolo 11, I-34131 Trieste,
Italy
12INFN - Sezione di Trieste, Via Valerio 2, 34100 Trieste TS, Italy
13IFPU - Institute for Fundamental Physics of the Universe, via Beirut 2,
34151, Trieste, Italy
14 ICSC - Centro Nazionale di Ricerca in High Performance Computing, Big Data
e Quantum Computing, Via Magnanelli 2, Bologna, Italy
15Max-Planck-Institut fur Astrophysik, Karl-Schwarzschild Strasse 1, 85748
Garching, Germany
16Ruhr University Bochum, Faculty of Physics and Astronomy, Astronomical
Institute (AIRUB), German Centre for Cosmological Lensing, 44780 Bochum,
Germany
17Department of Physics and Astronomy, University College London, Gower
Street, London WC1E 6BT, UK
E-mail<EMAIL_ADDRESS>
(Accepted XYZ, 2023. Received XYZ, 2023; in original form XYZ, 2023)

###### Abstract

We analyse the fourth data release of the Kilo Degree Survey (KiDS-1000) and
extract cosmological parameter constraints based on the cosmic shear peak
count statistics. Peaks are identified in aperture mass maps in which the
filter is maximally sensitive to angular scales in the range 2-4 arcmin,
probing deep into the non-linear regime of structure formation. We interpret
our results with a simulation-based inference pipeline, sampling over a broad
$w$CDM prior volume and marginalising over uncertainties on shape calibration,
photometric redshift distribution, intrinsic alignment and baryonic feedback.
Our measurements constrain the structure growth parameter and the amplitude of
the non-linear intrinsic alignment model to
$\Sigma_{8}\equiv\sigma_{8}\left[\Omega_{\rm
m}/0.3\right]^{0.60}=0.765^{+0.030}_{-0.030}$ and $A_{\rm
IA}=0.71^{+0.42}_{-0.42}$, respectively, in agreement with previous KiDS-1000
results based on two-point shear statistics. These results are robust against
modelling of the non-linear physics, different scale cuts and selections of
tomographic bins. The posterior is also consistent with that from the Dark
Energy Survey Year-1 peak count analysis presented in Harnois-Déraps et al
(2021), and hence we jointly analyse both surveys with a common pipeline. We
obtain $\Sigma_{8}^{\rm joint}\equiv\sigma_{8}\left[\Omega_{\rm
m}/0.3\right]^{0.57}=0.759^{+0.020}_{-0.017}$, in agreement with the Planck
$w$CDM results. The shear-CMB tension on this parameter increases to
$3.1\sigma$ when forcing $w=-1.0$, and to $4.1\sigma$ if comparing instead
with $S_{8,\Lambda{\rm CDM}}^{\rm joint}=0.736^{+0.016}_{-0.018}$, one of the
tightest constraints to date on this quantity. Residual biases in the
photometric redshifts of the DES-Y1 data and in the modelling of small scales
physics could lower this tension, however it is robust against other
systematics. Limits in the accuracy of our emulator prevent us from
constraining $\Omega_{\rm m}$.

###### keywords:

Gravitational lensing: weak – Methods: data analysis, numerical – Cosmology:
dark matter, dark energy & cosmological parameters

††pubyear: 2023††pagerange: KiDS-1000 and DES-Y1 combined: Cosmology from peak
count statistics–20

##  1 Introduction

Cosmic shear cosmology has entered an era of high precision, with recent
measurements from the Kilo Degree Survey111KiDS:kids.strw.leidenuniv.nl
(KiDS), the Dark Energy Survey 222DES:www.darkenergysurvey.org (DES) and the
Hyper Suprime Camera Survey333HSC:www.naoj.org/Projects/HSC (HSC) reaching a
precision of a few percent on parameters central to the standard model of
cosmology (e.g. Asgari et al., 2021; van den Busch et al., 2022; Li et al.,
2023a; Secco et al., 2022b; Amon et al., 2022; Dalal et al., 2023; Li et al.,
2023b). Based on the detection of weak correlations between the observed
shapes of galaxies imparted by the foreground large scale structure, cosmic
shear is mostly sensitive to the structure growth parameter
$S_{8}\equiv\sigma_{8}\sqrt{\Omega_{\rm m}/0.3}$, a combination of the matter
density parameter $\Omega_{\rm m}$ and of the amplitude of the linear matter
power spectrum smoothed on spheres of $8h^{-1}$Mpc, labelled as $\sigma_{8}$
(for lensing reviews, see e.g. Kilbinger, 2015; Mandelbaum, 2018). These
Stage-III lensing surveys have been steadily improving the data quality and
the analysis methods, in preparation for the next generation of cosmic shear
experiments such as the Rubin observatory444LSST: www.lsst.org (Ivezić et al.,
2019), Euclid555Euclid: www.euclid-ec.org (Laureijs et al., 2011) and the
Nancy Grace Roman space telescope666WFIRST: roman.gsfc.nasa.gov (Akeson et
al., 2019).

Despite the large effort that is being invested by international
collaborations in constructing accurate lensing catalogues of hundreds of
millions of galaxies, it is not entirely clear how to best analyse these vast
data, striking an optimal compromise between accuracy and precision. To date
the shear two-point (2pt) functions are still regarded as the baseline summary
statistics, having been tested for over a decade and achieving an unmatched
level of understanding and control in all aspects of the analysis, including
measurements tools (e.g. TreeCorr and NaMaster, see Jarvis et al., 2004;
Alonso et al., 2019), theoretical predictions (e.g. Kilbinger et al., 2017)
and the impact of systematics (see e.g. Mandelbaum, 2018). The main drawback
from these statistics is that they completely disregard the non-Gaussian
information that is stored in the non-linear matter field, more precisely in
the coupling between the phases of distinct Fourier modes, without which the
cosmic web would look like a Gaussian random field. This is obviously sub-
optimal, and this waste of information will be aggravated in the upcoming
cosmic shear experiments. Accessing this non-Gaussian information is an active
field of research: an array of novel weak lensing statistics are being
developed specifically to utilise this complementary small-scale information.
These new methods are reaching a level of maturity that makes them competitive
at analysing existing cosmic shear data, carefully balancing the precision vs
accuracy metric. Recent progress is largely due to the radically improved
modelling of the signal, thanks to the increased accuracy of cosmological
$N$-body codes and the availability of super-computers (see Angulo & Hahn,
2022, for a recent review on $N$-body codes). Recent examples of these
‘beyond-2pt’ cosmic shear data analyses include the three-point function (Fu
et al., 2014; Secco et al., 2022a; Burger et al., 2023), peak count statistics
(Kacprzak et al., 2016; Martinet et al., 2018; Shan et al., 2018; Harnois-
Déraps et al., 2021a; Zürcher et al., 2022; Liu et al., 2023; Gatti et al.,
2023a; Marques et al., 2023), density split statistics (Brouwer et al., 2018;
Gruen et al., 2018; Burger et al., 2022), shear clipping (Giblin et al.,
2018), persistent homology (Heydenreich et al., 2022), moments of convergence
maps (van Waerbeke et al., 2013; Gatti et al., 2020), cumulative distribution
functions (Anbajagane et al., 2023), likelihood-free inference (Jeffrey et
al., 2021; Lin et al., 2023; Gatti et al., 2023b) or convolutional neural
network inference (Fluri et al., 2019, 2022).

At the moment, these alternative methods exhibit a constraining power that is
similar to that of two-point functions, which is not surprising given the
noise levels of current lensing data, which make difficult the extraction of
information stored in the noisy higher-order moments. The situation will
change drastically with the upcoming surveys, where the cosmic web itself will
be detectable with lensing, at which point the non-Gaussian information will
take on a larger proportion of the signal.

All forecasts are clear about this: joint cosmic shear analyses that combine
two-point functions and any complementary probe improve the constraints on
cosmological parameters even in presence of systematic uncertainties (e.g. Li
et al., 2019; Schneider et al., 2019; Zürcher et al., 2020; Pyne & Joachimi,
2021; Harnois-Déraps et al., 2021b; Giblin et al., 2023; Euclid Collaboration:
Ajani et al., 2023). The main difficulty in many of these methods comes from
their accrued dependence on numerical simulations, which adds a significant
computational overhead to the data analysis compared to those for which an
analytical model exists. Typically, simulations are needed for modelling the
cosmological signal, for modelling some of the systematics such as baryonic
feedback or intrinsic alignments of galaxies, and for the estimation of the
covariance matrix (although this is not always necessary, as demonstrated by
the recent likelihood-free inference analyses mentioned above).

This paper contributes an important step to this effort: we carry out a
cosmological analysis based on lensing peak statistics measured from the
fourth data release of the Kilo Degree Survey (KiDS-1000 hereafter). We use
the exact same data as those used in the two-point function analyses of Asgari
et al. (2021, A21 hereafter), while ignoring for now other non-lensing KiDS
galaxy catalogues designed for galaxy clustering analyses (Vakili et al.,
2020; Bilicki et al., 2021). Our method finds peaks in aperture mass maps with
an aperture filter designed for the extraction of small-scale structure, with
maximal sensitivity to scales of less than 4 arcmin, as in Martinet et al.
(2018, hereafter M18) and Harnois-Déraps et al. (2021a, HD21). This contrasts
with the recent peak count analysis of Zürcher et al. (2022), in which peaks
are extracted from convergence maps with pixel resolution of about 7 arcmin.
Both methods have their advantages and downsides, ours strongly focuses on
small, non-linear scales, which, as demonstrated in HD21 and Martinet et al.
(2020), have a higher potential for complementarity with two-point functions.
Finding a posterior that is statistically consistent with that from HD21, we
combine both likelihoods and carry out a joint KiDS-1000 + DES DR1 data
(DES-Y1 hereafter) peak count analysis, finding the tightest constraints on
$S_{8}$ to date from peaks alone.

After describing the data and simulations in Sec. 2, we detail our measurement
techniques and analysis pipeline in Sec. 3, and we present our mitigation
strategy for the key systematic uncertainties in Sec. 4. We finally show our
results in Sec. 5 and discuss our findings afterwards. Supplementary material
is provided in the Appendices, including a thorough discussion of $B$-modes
(in Appendix A), supplementary pipeline validation tests (in Appendix B) and a
detailed discussion on goodness-of-fit for noisy covariance matrices (in
Appendix C).

##  2 Data and Simulations

We present in this section the survey data and the various simulation suites
that are used for the cosmological analysis.

###  2.1 KiDS-1000 data

The Kilo Degree Survey (Kuijken et al., 2015) is a multi-band photometric
galaxy survey explicitly designed for weak lensing cosmology. Carried out at
the European Southern Observatory by the VST-OmegaCAM, we analyse here the
public777KiDS-1000 data:http://kids.strw.leidenuniv.nl/DR4 fourth data release
(Kuijken et al., 2019). The observation conditions are of exceptional quality,
with a mean seeing of 0.7 arcsec in the $r$-band, used for shape measurements.
The photometric redshifts are obtained from a combination of nine optical and
infrared bands ($ugriZYJHK_{\rm s}$, see Wright et al., 2020), thanks to the
observations of the companion VIKING survey (VISTA Kilo-degree INfrared
Galaxy, Edge et al., 2013). The galaxies selected in this analysis exactly
match those used in the cosmic shear two-point function analyses of A21 and
van den Busch et al. (2022), covering an effective area of 777.4 deg2.

The KiDS DR4 data are reduced with the theli (Erben et al., 2013) and Astro-
WISE (Begeman et al., 2013) pipelines, following which the shear is inferred
from lensfit (Miller et al., 2013; Fenech Conti et al., 2017). Shear additive
and multiplicative biases ($c$\- and $m$-corrections) are measured to a high
accuracy (Giblin et al., 2021), where it is shown via a series of null tests
that known residual systematics in the shear measurement could lead to no more
than a 0.1$\sigma$ shift in the structure growth parameter
$S_{8}\equiv\sigma_{8}\sqrt{\Omega_{\rm m}/0.3}$, the composite quantity that
is best measured by cosmic shear. Note that strictly speaking, the results
from the tests carried out in Giblin et al. (2021) are only shown to hold for
two-point cosmic shear statistics.

Following A21, we split the full DR4 galaxies in five tomographic bins
according to their individual best-fit redshift $z_{\rm B}$ as measured by bpz
(Benítez, 2000), with bin edges set to [0.1, 0.3, 0.5, 0.7, 0.9 and 1.2]. The
tomographic redshift distributions, $n^{a}(z)$, are estimated via self-
organising maps (SOM, see Wright et al., 2020), which group galaxies based on
their nine-band photometric properties and assign redshifts based on similar
studies made on spectroscopic samples; galaxies for which no match is found
are rejected. We further reject galaxies for which the SOM redshift
catastrophically differs from the initial $z_{\rm B}$, resulting in the so-
called ‘Gold Sample’ introduced in Hildebrandt et al. (2021) and used in the
subsequent KiDS-1000 cosmic shear analyses mentioned above. As detailed in
A21, the means and the error of the SOM redshift distributions are calibrated
on KiDS-like mock data constructed from the MICE2 simulations (Fosalba et al.,
2015; van den Busch et al., 2020) and accounted for during the inference stage
of our analysis. The redshift accuracy is excellent due to the nine-band
photometry, which helps breaking degeneracies in the galaxy spectral energy
distributions: at worst, the difference between the mean redshift and that
estimated from the matched spectroscopic sample is $z_{\rm est}-z_{\rm
true}=0.013\pm 0.0118$, making this a sub-dominant source of uncertainty in
our measurement. Note that Hildebrandt et al. (2021) further show that the SOM
redshift distributions are fully consistent with independent estimates based
on clustering cross-correlations with spectroscopic reference samples,
providing extra robustness to the method. Fig. 1 shows the redshift
distributions estimated in the five tomographic bins, along with the
variations on these distributions allowed within our photometric uncertainty.

Figure 1: Tomographic redshift distribution of the KiDS-1000 data. The thinner lines represent the effect of photometric uncertainty on these distributions, characterised by $n^{a}(z)\rightarrow n^{a}(z+\Delta z_{a})$, with $\Delta z_{a}$ sampled 10 times from Gaussian distributions with widths listed in Table 1. All shifted $n^{a}(z)$ are then rebinned with the same $z$ bins. Table 1: Main properties of the KiDS-1000 data used in the current work. The Gold Sample redshift selection based on $z_{\rm B}$ is identical to that presented in Hildebrandt et al. (2021). The effective number densities are listed in the second column, in gal/arcmin2. The shape noise (per component) listed in the third column reflects the dispersion measured in the observed galaxy shapes, as documented in Giblin et al. (2021), while the fifth column shows the mean shape calibration coefficients. The redshift bias and errors listed in the fourth column are estimated from the SOM method in Hildebrandt et al. (2021), while the last column shows the additive $c_{1/2}$ terms, which has an uncertainty of $0.23\times 10^{-3}$ (Giblin et al., 2021).  tomo | $n_{\rm eff}$ | $\sigma_{\epsilon}$ | $z_{\rm est}-z_{\rm true}$ | $m$ | $(c_{1},c_{2})\times 10^{3}$
---|---|---|---|---|---
bin1 | 0.62 | 0.27 | $0.000\pm 0.0106$ | $-0.009\pm 0.019$ | (0.295, 0.156)
bin2 | 1.18 | 0.26 | $0.002\pm 0.0113$ | $-0.011\pm 0.020$ | (0.004, 0.621)
bin3 | 1.85 | 0.27 | $0.013\pm 0.0118$ | $-0.015\pm 0.017$ | (0.052, 0.728)
bin4 | 1.26 | 0.25 | $0.011\pm 0.0087$ | $\,\,\,\,0.002\pm 0.012$ | (-0.360, 0.948)
bin5 | 1.31 | 0.27 | $-0.006\pm 0.0097$ | $\,\,\,\,0.007\pm 0.010$ | (-1.363, 1.155)

The SOM selection and the shear inference pipelines are both repeated on KiDS-
like image simulations (Kannawadi et al., 2019), from which a relation between
apparent size, magnitude and the observed galaxy shape is used to calibrate
the inferred lensfit shear888A KiDS-1000 re-analysis has been presented in Li
et al. (2023a) after correcting for an anisotropic error in the lensfit
likelihood sampler. This error has not been corrected here, but their study
shows the correction has a negligible impact on the inferred cosmology..
Whereas previous cosmological analyses use a single $m$-calibration factor per
tomographic bin, the aperture mass map statistics exploited in this paper are
subject to local variations in the noise levels and seeing conditions, and we
therefore use the above-mentioned relation to extract a shear calibration per
object, $m_{a}$. This is not necessary, but allows to capture possible
correlations between the $m$-correction and the lensfit weights. These are
inevitably noisier than the average over the full tomographic bins, but a
large fraction of this noise cancels within our aperture mass map calculations
as well, while providing optimal estimates of the local noise contribution
(M18). Let us recall that this calibration corrects for known residual biases
such as shape detection biases (Fenech Conti et al., 2017; Kannawadi et al.,
2019) or blending of the images of galaxies (Hoekstra et al., 2015). While we
apply the $m$-correction per object, the averaged multiplicative biases per
redshift bin used in A21 enter our analysis at the inference level in the form
of nuisance parameters over which we marginalise. Table 1 summarises the
survey properties relevant to our analysis.

###  2.2 DES-Y1 data

The DES-Y1 measurements is based on the public year-1 data release from the
Dark Energy Survey Collaboration (Abbott et al., 2018), with source galaxy
selections that exactly follow the main cosmic shear results described in
Troxel et al. (2018). The lensing catalogue consists of 26 million galaxies
covering a footprint of 1320 deg2 with a galaxy density of 5.07 gal arcmin-2.
The per-galaxy shear signal is inferred with the Metacalibration method
(Sheldon & Huff, 2017). Every galaxy is assigned to one of the four
tomographic bins based on the photometric redshift posteriors estimated from
the the griz flux measurements, as detailed in Hoyle et al. (2018). Following
Troxel et al. (2018), the uncertainty on the shear multiplicative calibration
is given by $m_{a}=0.012\pm 0.023$.

Whereas the original DES-Y1 results estimated the tomographic $n^{a}(z)$ from
a Bayesian photometric redshift analysis calibrated on the COSMOS2015 field
(Laigle et al., 2016), the HD21 reanalysis instead opted for $n^{a}(z)$
estimates based on a direct reweighted calibration of matched spectroscopic
data (Lima et al., 2008, DIR hereafter), following the DES-Y1 reanalyses of
Joudaki et al. (2020) and Asgari et al. (2020). The uncertainty on the DIR
mean redshift distributions is $\Delta z_{a}$ = [0.008, 0.014, 0.011 and
0.009] for redshift bins $a=1...4$, respectively. Both methods have their pros
and cons. The calibration with COSMOS is by design based on a complete sample
but suffers from imperfect redshifts (see e.g. Alarcon et al., 2021). In
contrast, the spec-$z$ samples used for the DIR method have (close to) perfect
redshifts but are incomplete and not representative of the source sample,
which is alleviated by the reweighting, but often cannot be fully eliminated
(see Gruen & Brimioulle, 2017). Importantly, the DIR $n^{a}(z)$ favours
$S_{8}$ values that are smaller by $\Delta S_{8}=0.03$ compared to the COSMOS-
calibrated $n(z)$, which is a 0.8$\sigma$ shift (Joudaki et al., 2020).

###  2.3 Simulations

As mentioned in the introduction, the accuracy of simulation-based inference
pipelines fully depends on the quality of the numerical simulations it is
calibrated on. The same way 2pt analyses must carefully understand the scales,
cosmologies and redshifts that are well captured by their model, it is
critical for our peak count analysis to identify the range of validity of our
training simulations. The additional complexity here is that no simulation
suite serves all purposes, and therefore we must carefully investigate, for
all of them separately, the accuracy and limits of the measurements and how
these impact the peak count statistics. The simulations used in this work are
in many aspects identical to those presented in HD21, which we refer to for
further details. Specifically:

  1. 1.

the cosmological dependence of the peak count statistics is calibrated on the
$w$CDM cosmo-SLICS $N$-body simulations introduced in Harnois-Déraps et al.
(2019). They sample a wide volume in $S_{8}$, $\Omega_{\rm m}$, $w_{0}$ and
$h$ with 25 points arranged in a Latin hypercube (plus one $\Lambda$CDM
point), each evolved with a pair of $N$-body simulations designed to suppress
sample variance in 2pt functions, then ray-traced in ten light-cones of 100
deg2 ($10000$ deg2 in total area). These form our Cosmology Training Set, and
resolve the non-linear physics to better than 2% up to $k$-modes of 2.0
$h^{-1}$Mpc, when compared to the Cosmic Emulator (Heitmann et al., 2014).
Smaller scales gradually lose precision, affecting mainly their ability to
resolve substructure in most massive objects. The exact impact of this loss on
weak lensing peak counts is investigated in HD21 with a separate set of
simulations ran with a much higher force resolution, where it is found that
this leads to at most a 1% loss of the highest peaks, which is largely sub-
dominant compared to both baryonic physics and statistical errors. We revisit
this in Sec. 4 (see also point iv);

  2. 2.

the covariance matrix that captures the sample variance is estimated from 124
fully independent SLICS $N$-body simulations described in Harnois-Déraps & van
Waerbeke (2015). These are evolved from independent initial conditions at a
fixed cosmology, and make our Covariance Training Set. They resolve the same
non-linear physics as the cosmo-SLICS, and are shown in Harnois-Déraps et al.
(2019) and HD21 to produce marginalised errors on cosmological parameters that
are fully consistent with those obtained with an analytical calculation, when
analysing 2pt statistics. Burger et al. (2022) further show in the context of
density-split statistics that a covariance matrix estimated from the SLICS or
from a much larger number of log-normal Flask mocks (Xavier et al., 2016)
produce fully consistent results, as expected for these mildly non-linear
statistics. We further increase the effective number of covariance mocks by
randomly rotating 10 times the shape noise components. This works particularly
well given that the peak statistics is currently shape-noise dominated: while
the expectation value of standard 2pt statistics does not depend on the noise
(only their variance does), shape noise affects both the signal and covariance
of map-based statistics (see Appendix D of Heydenreich et al., 2021);

  3. 3.

for the KiDS-1000 analysis, the impact of galaxy intrinsic alignments is
measured from the IA-infused lensing simulations described in Harnois-Déraps
et al. (2021b). These are also constructed from the cosmo-SLICS and therefore
resolve the same physical scales. This IA Training Set assumes a linear
coupling between the projected non-linear tidal field and the intrinsic
ellipticity of every galaxy, and is therefore physically modelling the non-
linear linear alignment model of Bridle & King (2007) without explicit
redshift nor luminosity dependence. It is expected that this effective IA
model does not fully capture the alignment signal, and that a more physical
model such as the Tidal Alignment And Torquing model (Blazek et al., 2019) or
the halo-model of Fortuna et al. (2021) would provide a more accurate
description, however current cosmic shear surveys do not have the statistical
power to constrain parameters beyond the simpler NLA model (Secco et al.,
2022b), which is therefore deemed sufficient for the current analysis. The IA
infusion process has been shown in Harnois-Déraps et al. (2021b) to accurately
reproduce the NLA predictions for the 2pt correlation function down to scales
of a few arcmin, beyond which the NLA is expected to fail in a manner that is
undetectable in the current data. Burger et al. (2023) further show that these
same simulations agree with the IA modelling of three-point shear statistics.
The model fails at scales that correspond to high over-densities in our
simulations, which contribute to lensing peaks that are excluded from our
analysis. We infuse different levels of IA and marginalise over these choices
in the end, as described in Sec. 4. for the DES-Y1 analysis, IA are included
with a non-linear halo-based model, see HD21 for details;

  4. 4.

limits in the force resolution of the cosmo-SLICS are bound to impact the weak
lensing statistics in a manner that is not always predictable. We assess this
with the SLICS-HR suite (Harnois-Déraps & van Waerbeke, 2015), a high-
resolution version of the SLICS light-cones recently used in a combined
lensing-clustering cosmological analysis (Duncan et al., 2022). The SLICS-HR
consist of ten independent $10\times 10$ deg2 catalogues that are run at the
same cosmology and with the same particle count and volume as the SLICS, but
the $N$-body force accuracy has been increased such as to resolve $k$-modes up
to 10 $h^{-1}$Mpc. We use these to validate the full inference pipeline in
Sec. 4.7, acting as our Validation Set;

  5. 5.

the impact of baryon feedback is estimated with the Magneticum hydrodynamical
simulations999Magneticum simulations: www.magneticum.org, forming our Baryons
Training Set. These have been shown to reproduce a number of key observations
relevant to weak lensing studies (Castro et al., 2018), and notably the
feedback on the matter distribution closely matches that of the BAHAMAS
(McCarthy et al., 2017), another suite of hydro simulations with independent
prescriptions for their sub-grid physics. The training set consists of ten
$10\times 10$ deg2 pseudo-independent light-cones extracted101010The
Magneticum light-cones were built with the public slicer code:
https://github.com/TiagoBsCastro/SLICER from full hydrodynamical simulations,
and another 10 light-cones extracted from dark matter-only sister simulations,
evolved from the same initial conditions (more details on the used simulations
can be found in Martinet et al., 2021). There is a large uncertainty on the
exact impact of baryonic physics on the matter distribution (and therefore on
our lensing statistics), which we account for by linearly scaling the relative
baryonic bias with a nuisance parameter, $b_{\rm bary}$, which we marginalise
over at the inference stage111111Note that this parameter is not to be
confused with $A_{\rm bary}$ used in A21 (see their Table 2), which
specifically relates to one of the free parameters entering their HMCode halo
model..

  6. 6.

different $N$-body codes and ray-tracing methods, even at fixed cosmology,
will have a residual impact on the peaks statistics (Hilbert et al., 2020). We
explore these numerical systematics with the public full-sky weak lensing
simulations from Takahashi et al. (2017, T17 hereafter)121212T17:
http://cosmo.phys.hirosaki-u.ac.jp/takahasi/allsky_raytracing/, post-processed
into KiDS-1000 mock data (North and South patches) as in Burger et al. (2023).
The T17 simulations follow the non-linear evolution of 20483 particles in a
series of nested cosmological volumes with side length starting at $L$ =
450$h^{-1}$Mpc at low redshift, then increasing at higher redshifts. These
result in 108 pseudo-independent full-sky lensing maps, seven of which are
used in this work, with flat $\Lambda$CDM cosmological parameters set to
$(\Omega_{\rm m},\Omega_{\rm b},\sigma_{8},h,n_{\rm s},w_{0})=$(0.279, 0.046,
0.82, 0.7, 0.97, -1.0).

####  2.3.1 Assembling mock surveys

Most of these simulations have been introduced in HD21, in Heydenreich et al.
(2022), in Burger et al. (2022) and in the references listed in the previous
section; we encourage the interested reader to consult these for a more
complete technical description. To summarise some of the key properties, all
of the above-mentioned simulations are organised in light-cones of 100 deg2
each, populated with galaxy samples that match the tomographic $n^{a}(z)$
distributions, number densities and shape noise levels of the KiDS-1000 Gold
Sample and DES-Y1 data. Except for the IA-infused simulations, the galaxy
positions, the amplitude of their ellipticities $|\boldsymbol{\epsilon}_{\rm
data}|$ and their multiplicative shear calibration factors $m_{a}$ are exactly
reproduced in each of the mock survey realisations (i.e. in the Cosmology,
Covariance, Validation and Baryon Training Sets). To achieve this, the
KiDS-1000 data are split into 18 tiles that each fit within 100 deg2 regions,
as depicted in Fig. 2, and the shear and convergence from every simulation is
repeated across them, interpolated at the local galaxy positions. These tiles
are analysed separately and combined only at the level of the summary
statistics, ensuring that cross-tile correlations that exist in the data but
not in the simulations are explicitly ignored. This effect is minor for
localised non-Gaussian probes such as peak statistics, but is critical for
e.g. shear 2pt functions. The shear and convergence are interpolated from the
underlying simulations at the position of every galaxy, infused with the
$m_{a}$ from the data, then combined with the (randomly rotated) observed
ellipticity following:

$\displaystyle\boldsymbol{\epsilon}_{\rm
mock}=\frac{{\boldsymbol{\epsilon}}_{\rm data}^{\rm
rand}+{\boldsymbol{g}}}{1+{\boldsymbol{\epsilon}}_{\rm data}^{\rm
rand}{\boldsymbol{g}}^{*}}\,.$ (1)

In the above expressions, bold-font symbols are spin-2 complex quantities and
$\boldsymbol{g}$ is the $m$-biased simulated reduced shear. As described in
HD21, this involves rotating each tile at the equator, which preserves the
relative positions of galaxies but modifies their ellipticities, defined with
respect to the North pole.

Figure 2: Tiling strategy adopted to pave the full KiDS-1000 data (black) with
flat-sky $10\times 10$ deg2 simulations (red squares). Some of the tiles
slightly overlap due to the sky curvature, in which case the data is split at
the mean Dec in overlapping regions.

We repeat this construction for all light-cones of the Cosmology Training Set,
the Covariance Training Set, the Baryon Training Set and the Validation Set.
Additionally, the uncertainty in the photometric redshifts is forward-modelled
with a further 10 full survey realisations computed at the fiducial cosmology,
in which the $n(z)$ is shifted by small amounts (details provided in Sec. 4).
In total, this results in 414 simulated mosaic surveys that we analyse in
preparation for the inference stage, with the majority (260) contributing to
the Cosmology Training Set. Each mock further contains 10 random rotations of
$\boldsymbol{\epsilon}_{\rm data}$, to improve convergence of the signal.

The IA Training Set are treated slightly differently, since for these the
positions of the mock galaxies must be sampled from the simulated over-density
maps or halo catalogues, which do not correlate with the positions in the data
(see Harnois-Déraps et al., 2021b, for more details). The mosaic survey tiling
is therefore not possible, so we use instead 5 light-cones per IA model and
explore four alignment strengths in KiDS, and one model in DES. Although these
represent a lesser total area than the real data, their sole purpose is to
capture the relative impact of IA on the signal, computed from ratios in which
the sample variance cancels by design. We use $4\times 5\times 100$ deg2 of
training data, which is enough to capture this.

The simulations are free of additive biases by construction, however Giblin et
al. (2021) measures residual additive terms $c_{1/2}$ in the KiDS-1000 cosmic
shear catalogues, caused by the shape measurement method itself. These are
reported in Table 1 and subtracted from the observed ellipticities when
analysing the real data. We follow Troxel et al. (2018) by not accounting for
any low-level additive terms in the DES Y1 catalogue. The multiplicative
biases are not easily removed from the data, hence we instead infuse the mocks
with the $m_{a}$ terms per object, and treat thereafter data and simulations
on equal footings.

##  3 Methods

###  3.1 Aperture mass map statistics

There exists a number of methods for identifying and counting lensing peaks,
including finding maxima on convergence maps (Li et al., 2019), on wavelet-
transformed maps (Ajani et al., 2020) or on aperture mass maps (Schneider,
1996). We here opted for the aperture mass maps for the following reasons: as
argued in M18, this statistics is immune to masking-induced biases and strong
$B$-mode leakage common to methods based on reconstruction of convergence
maps, plus it benefits from a local estimation of the shape and Poisson noise,
yielding more accurate signal-to-noise maps.

Specifically, we cover each of the 18 tiles with a two-dimensional grid with a
pixel size of 0.59 arcmin. We next reconstruct the mass inside an aperture
filter $Q$ centred on each pixel, at position $\boldsymbol{\theta}$ on the
sky, from the sum of all tangential ellipticities $\epsilon_{a,{\rm t}}$
contained therein as:

$\displaystyle M_{\rm ap}(\boldsymbol{\theta})=\frac{1}{n_{\rm
gal}(\boldsymbol{\theta})\sum_{a}w_{a}(1+m_{a})}\sum_{a}w_{a}\epsilon_{a,{\rm
t}}({\boldsymbol{\theta}},{\boldsymbol{\theta}}_{a})Q(|{\boldsymbol{\theta}}-{\boldsymbol{\theta}}_{a}|,\theta_{\rm
ap},x_{c})\,.$ (2)

The tangential ellipticity about $\boldsymbol{\theta}$ is computed as
$\epsilon_{a,{\rm
t}}({\boldsymbol{\theta}},{\boldsymbol{\theta}}_{a})=-[\epsilon_{1}({\boldsymbol{\theta}}_{a})\
{\rm
cos}(2\phi({\boldsymbol{\theta}},{\boldsymbol{\theta}}_{a}))+\epsilon_{2}({\boldsymbol{\theta}}_{a})\
{\rm sin}(2\phi({\boldsymbol{\theta}},{\boldsymbol{\theta}}_{a}))]$, where
$\boldsymbol{\theta}_{a}$ is the position of galaxy $a$ and
$\phi({\boldsymbol{\theta}},{\boldsymbol{\theta}}_{a})$ is the angle between
both coordinates. The sum runs over all galaxies in the aperture, and $n_{\rm
gal}(\boldsymbol{\theta})$ is the local galaxy density in the filter when
centred at $\boldsymbol{\theta}$. As in M18 and HD21, our filter
$Q(\theta,\theta_{\rm ap},x_{c})$, abridged to $Q(\theta)$, matches that of
Schirmer et al. (2007), which is optimised for efficiently detecting NFW
haloes:

$\displaystyle Q(x)=\frac{{\rm tanh}(x/x_{c})}{x/x_{c}}\big{[}1+{\rm
exp}(6-150x)+{\rm exp}(-47+50x)\big{]}^{-1}.$ (3)

In the above expression, we use the standard value of $x_{c}=0.15$, while
$x=\theta/\theta_{\rm ap}$, with $\theta$ the distance to the filter centre.
We additionally use the same filter size, set to $\theta_{\rm ap}=12.5$
arcmin, which is shown in M18 to better detect the cosmological signal over
other filter sizes in KiDS data. We compute Eq. (2) at every pixel location to
construct our signal map. The variance about this map is computed at every
pixel location with:

$\displaystyle\sigma^{2}_{\rm ap}(\boldsymbol{\theta})=\frac{1}{2n_{\rm
gal}^{2}(\boldsymbol{\theta})\left[\sum_{a}w_{a}\right]^{2}}\sum_{a}w_{a}^{2}|{\boldsymbol{\epsilon}}_{a}|^{2}Q^{2}(|{\boldsymbol{\theta}}-{\boldsymbol{\theta}}_{a}|)\,,$
(4)

where again the sum runs over all galaxies in the filter. The $m$-calibration
estimated from the image simulations of Kannawadi et al. (2019) is meant to
correct the inferred shear, not the ellipticity, which explains why it appears
in the denominator of Eq. (2) but not in that of Eq. (4), which describes the
noise map. Finally, we take the ratio between Eq. (2) and the square root of
Eq. (4) at every pixel location to construct our signal-to-noise maps,
$\mathcal{S}/\mathcal{N}(\boldsymbol{\theta})\equiv M_{\rm
ap}(\boldsymbol{\theta})/\sqrt{\sigma_{\rm ap}^{2}(\boldsymbol{\theta})}$,
from which we identify peaks as local maxima and record their
$\mathcal{S}/\mathcal{N}$-values. We repeat this process for the 10
realisations of random rotations and report the average, except for the
Covariance Training Set, for which we do not take the average; instead, each
noise realisation leads to an estimate of the covariance matrix, of which we
take the average in the end.

As detailed in HD21, masking is dealt with naturally in aperture mass
statistics, and no special treatment needs to be enforced as long as data and
simulations are masked and analysed the same way. This is achieved by fixing
galaxy positions in the simulations to that of the observed data, which
ensures the impact of the mask is identical. In our case, we decided
nevertheless to act upon masked pixels. These are identified from the galaxy
catalogues as regions with an aperture galaxy density that is either
critically low or null, then removed from the final
$\mathcal{S}/\mathcal{N}(\boldsymbol{\theta})$ maps.

It has been shown that some additional information can be extracted by
combining the peak count statistics measured from multiple filter sizes (e.g.
Zürcher et al., 2022; Giblin et al., 2023), however M18 shows that this gain
is mild for Stage III surveys. We therefore opted for a single-scale analysis
here, but intend to revisit this in the future.

###  3.2 Tomography and selection

Tomographic decomposition of the lensing data allows us to probe the redshift
evolution of the large scale structures, which is largely driven by
$\Omega_{\rm m}$ and $w_{0}$ via their impact on the growth of perturbations.
A direct consequence of the improved sensitivity to these is a gain in
precision in $S_{8}$, arising from degeneracy breaking. This decomposition is
different for the KiDS and DES surveys, which we detail here.

####  3.2.1 KiDS-1000

From the five KiDS tomographic bins, we include both the auto- and the cross-
redshift measurements, as first defined in Martinet et al. (2020). To be
specific, peaks are identified from the individual tomographic galaxy
catalogues (the ‘auto’ redshift bins 1, 2, 3, 4, 5), from every possible
combination of bin pairs (1$\cup$2, 1$\cup$3, 1$\cup$4, 1$\cup$5… 4$\cup$5),
triplets (1$\cup$2$\cup$3, 1$\cup$2$\cup$4, 1$\cup$2$\cup$5…), quadruplets
(1$\cup$2$\cup$3$\cup$4, 1$\cup$2$\cup$4$\cup$5…) and quintets (i.e. no
tomography). As shown in HD21, Zürcher et al. (2022) and Heydenreich et al.
(2022), these ‘cross-tomographic’ catalogues contain a significant amount of
additional information that is not contained within the ‘auto’ case. The
tomographic peak function is presented in Fig. 3, showing in the different
panels the 30 different redshift bin combinations. For each case we overlay
the predictions from the Cosmology Training Set in colour with the data
measurements in black; the error bars are obtained from the Covariance
Training Set. A similar measurement is presented in Fig 15, where the data is
replaced by the mean over our Baryons Training Set; it was shown in HD21 that
using instead the mean over the Covariance Training Set effectively validates
that the ensemble average agrees with the predictions, a result we verified
still holds with the KiDS-1000 mocks. In these figures, we have subtracted the
peak function measured from pure shape noise fields, $N^{\rm noise}_{\rm
peaks}$, to better highlight the cosmological dependence of the signal.

In all cases, we measure the peak function in $\mathcal{S}/\mathcal{N}$ bins
of width 0.5 in the range [-2.5, 4.0], for a total of 13 bins per sub-panel
and 390 elements in total. The motivation behind this initial choice of range
is driven by a number of requirements, notably that of having a large number
of peaks per bin to ensure the data is Gaussian-distributed (with our
selection, every bin has at least 200 objects, while bins outside this range
have far fewer objects). Additionally, our analysis has strict requirements on
the modelling precision and on the level of contamination by residual
systematic effects, resulting in this bin selection being in fact
“aggressive”. We expand on this in Sec. 4, where we argue that instead the
range $[-1.0<\mathcal{S}/\mathcal{N}<3.0]$ is a better choice with lower
modelling errors, forming a ‘clean’ data vector of $7\times 30=210$ elements
in total that is used for the main cosmological analysis.

####  3.2.2 DES-Y1

Following HD21, our DES peak count analysis includes the auto- and cross-
redshift measurements up to pairs of tomographic bins, for a total of 10 bin
combinations. The peak function is measured in 12 $\mathcal{S}/\mathcal{N}$
bins in the range $[0.0<\mathcal{S}/\mathcal{N}<4.0]$, forming a data vector
with 120 elements. Although these details differ compared to the KiDS-1000
case described above, it is shown in HD21 to be accurate and competitive.

Figure 3: Tomographic weak lensing peak function $N_{\rm
peaks}^{\kappa}(\mathcal{S}/\mathcal{N})$ measured in the KiDS-1000 data
(black squares) and in the Cosmology Training Set simulations, colour-coded by
their $S_{8}$ value. The pure noise signal $N_{\rm peaks}^{\rm noise}$ has
been removed to better highlight the variations with respect to cosmology. The
panels show the results from different combinations of tomographic bins, in
which the red dashed lines represents the best-fit model inferred from our
fiducial analysis, see Sec. 5. Figure 4: This figure highlights the
correlations between the different elements of the KiDS-1000 data vector. From
left to right, the 30 blocks show the correlation coefficients for the
different redshift bin combinations, starting with singlets (i.e. auto-bins),
pairs, triplets, quadruplet and the no-tomographic case, with redshift
increasing towards the right and the top of the figure.

###  3.3 Analysis pipelines

Table 2: Priors used in the KiDS likelihood sampling. The ranges for the four cosmological parameters are determined by the cosmo-SLICS simulations, while the prescription from sampling the nuisance parameters describing the photometric redshifts $\Delta z_{a}$ and intrinsic alignments $A_{\rm IA}$ are taken from Joachimi et al. (2021). In particular, the redshift parameters are correlated and drawn from a multivariate Gaussian distribution with means $\boldsymbol{\mu}$ taken from Table 1 (fourth column) and a covariance matrix $\mathcal{C}_{z}$ described in Sec. 4.3. The shear calibration parameters $\Delta m_{a}$ are sampled from Gaussian priors centred on zero with a standard deviation $(\mu,\sigma)$ estimated in Giblin et al. (2021). The baryonic feedback parameter $b_{\rm bary}$ is used to scale the effect measured in the Baryon Training Set. Parameter | range | prior
---|---|---
Cosmology |  |
$\Omega_{\rm m}$ | [0.1, 0.55] | Flat
$S_{8}$ | [0.6, 0.9] | Flat
$h$ | [0.6, 0.82] | Flat
$w_{0}$ | [-2.0, -0.5] | Flat
Nuisance |  |
$\Delta z_{a}\times 10^{2}$ | [-10, 10] | $\mathcal{G}(\boldsymbol{\mu},\mathcal{C}_{z})$
$\Delta m_{1}\times 10^{2}$ | [-10, 10] | $\mathcal{G}(0.0,1.9)$
$\Delta m_{2}\times 10^{2}$ | [-10, 10] | $\mathcal{G}(0.0,2.0)$
$\Delta m_{3}\times 10^{2}$ | [-10, 10] | $\mathcal{G}(0.0,1.7)$
$\Delta m_{4}\times 10^{2}$ | [-10, 10] | $\mathcal{G}(0.0,1.2)$
$\Delta m_{5}\times 10^{2}$ | [-10, 10] | $\mathcal{G}(0.0,1.0)$
Astrophysics |  |
$A_{\rm IA}$ | [-5, 5] | Flat
$b_{\rm bary}$ | [0, 2] | Flat
Table 3: Priors used for sampling the nuisance parameters in the DES-Y1 peak statistics analysis. The sampling of photometric redshifts $\Delta z_{a}$ and shear bias $\Delta m_{a}$ nuisance parameters follows the original cosmic shear paper by Troxel et al. (2018). The baryonic feedback parameter $b_{\rm bary}$ is the same as in the KiDS-1000 likelihood, however there is no IA parameter here. Parameter | range | prior
---|---|---
$\Delta z_{1}\times 10^{2}$ | [-10, 10] | $\mathcal{G}(0.1,1.6)$
$\Delta z_{2}\times 10^{2}$ | [-10, 10] | $\mathcal{G}(1.9,1.3)$
$\Delta z_{3}\times 10^{2}$ | [-10, 10] | $\mathcal{G}(0.9,1.1)$
$\Delta z_{4}\times 10^{2}$ | [-10, 10] | $\mathcal{G}(1.8,2.2)$
$\Delta m_{a}\times 10^{2}$ | [-10, 10] | $\mathcal{G}(1.2,2.3)$
Astrophysics |  |
$b_{\rm bary}$ | [0, 2] | Flat

Our cosmological inference pipeline heavily builds on the methods presented in
HD21 and Heydenreich et al. (2022), which we briefly overview here. First, we
model the peak function by training a Gaussian Process Regression131313We use
the GPR toolkit provided by SciKitLearn (Pedregosa et al., 2011). emulator
(GPR) on the measurements from the Cosmology Training Set. The GPR can
subsequently produce $N^{\kappa}_{\rm peaks}$ predictions within a fraction of
a second everywhere inside the parameter volume covered by the cosmo-SLICS.
This therefore determines the prior ranges over $\Omega_{\rm m},S_{8},w_{0}$
and $h$, which we report in Table 2.

Second, we must estimate the covariance matrix, which captures the correlation
between the elements of our data vector, central to the error propagation. As
mentioned before, the Covariance Training Set consists of 124 full survey
realisations, each replicated with 10 shape noise realisation, producing 1240
pseudo-independent data vectors from which our covariance matrix C is
extracted. We show in Fig. 4 the cross-correlation coefficient matrix, defined
as $C_{ij}/\sqrt{C_{ii}C_{jj}}$, which better highlights the correlations
between the negative and positive peaks in each of the tomographic block. Also
visible is the significant amount of correlation (and anti-correlation)
present in the off-diagonal component. This matrix contains at most $390^{2}$
elements and is thus invertible (since 390 $<1240$, see Hartlap et al., 2007),
a criteria that is also naturally satisfied by the ‘clean’ KiDS-1000 data
vector, which contains only 210 entries, and by the DES-Y1 data vector, which
contains 120.

Having our model and covariance matrix, we are now in a position to evaluate
the likelihood $\mathcal{L}$ of the model $\boldsymbol{x}(\boldsymbol{\pi})$
with parameters $\boldsymbol{\pi}$, given the data vector $\boldsymbol{d}$. We
use the Sellentin & Heavens (2016) $t$-distribution likelihood, which is well
suited for nearly Gaussian data vectors with simulation-based covariance
matrices. It is constructed as:

$\displaystyle{\rm
ln}\mathcal{L}(\boldsymbol{\pi}|\boldsymbol{d})=\frac{N_{\rm sims}}{2}{\rm
ln}\bigg{[}1+\chi^{2}/(N_{\rm sims}-1)\bigg{]}+{\rm const},\mbox{\hskip
14.22636pt with}$ (5)
$\displaystyle\chi^{2}=[\boldsymbol{x}(\boldsymbol{\pi})-\boldsymbol{d}]^{\rm
T}{\rm C}^{-1}[\boldsymbol{x}(\boldsymbol{\pi})-\boldsymbol{d}].$ (6)

In the above, $N_{\rm sims}=1240$ is the number of realisations used to
evaluate the covariance matrix C. The model depends on the four cosmological
parameters $\Omega_{\rm m},S_{8},w_{0}$ and $h$, and on a set of 12 (9)
astrophysical and nuisance parameters for KiDS (DES), which characterise the
dependence of our signal on the systematic effects mentioned previously. This
is an excellent approximation to the more general likelihood suggested by
Percival et al. (2022) in our case. Finally, the posteriors are sampled both
by the nested sampling algorithm MultiNest (Feroz et al., 2009) and by
Nautilus (Lange, 2023), implemented within CosmoSIS (Zuntz et al., 2015).
While the latter sampler is more robust (Lange, 2023), the former has been
more widely used in the literature and is therefore useful to make fair
comparisons with previous analyses. We report from these the mean and 68%
credible intervals computed from the one-dimensional projected
posteriors141414We refer to these intervals as $1\sigma$ regions, even though
strictly speaking this notation should only apply to Gaussian posteriors., as
well as the maximum a posteriori for some of our key results.

Since our likelihood function differs from the widely used multi-variate
Gaussian, the goodness-of-fit evaluation must be adapted accordingly. For
Gaussian likelihoods, the $\chi^{2}_{\rm best-fit}$, estimated at the best-fit
parameters, is to a very good approximation sampling an underlying
$\chi^{2}_{\nu}$ probability distribution, which depends only on the number of
degrees of freedom $\nu$ – this is only an approximation however, because of
informative priors, non-linear modelling and correlated error bars (see e.g.
Joachimi et al., 2021). A good fit will have a $\chi^{2}_{\rm best-fit}$ close
to the maximum of the $\chi^{2}_{\nu}$ probability distribution, while a bad
fit will land far in the tail, leading to a probability to exceed (PTE) that
is smaller than our acceptance threshold, set to 0.01. For our Student-$t$
distribution likelihood, we still assess the goodness-of-fit with PTE values,
however the $\chi^{2}_{\nu}$ curve needs to be modified (see Appendix C for
details).

A few differences exist between the KiDS-1000 and DES-Y1 likelihoods which are
worth highlighting here, as these influence the construction of our joint
pipeline. First, the original DES-Y1 peak count analysis samples $\sigma_{8}$
instead of $S_{8}$; the latter is a better option as it exactly covers the
training volume and is therefore adopted for our DES-Y1 re-analysis. Second,
the treatment of intrinsic alignments are simpler in the DES-Y1 analysis: the
IA contribution is estimated from the alignments of dark matter haloes, which
are assumed to fully correlate with the alignment of central galaxies. This
non-linear prescription provides a single IA model that is then added to the
predictions, without marginalisation. As discussed in Sec. 5.2, not
marginalising over the IA in the KiDS analysis slightly underestimates the
total error. This is likely less important in the DES-Y1 likelihood since the
statistical error is larger.

This also connects with the third difference, which is that in the baseline
DES-Y1 measurement, only the auto-tomographic redshift bins are included, in
an attempt to avoid possible residual contamination from unmodelled IA in the
cross-redshift bins. This turns out to be an over-conservative data cut.
Indeed, the recent DES-Y1 persistent homology cosmic shear analysis from
Heydenreich et al. (2022) reveals that the constraints on $S_{8}$ are
negligibly affected by these IA terms: they show that a full tomographic
analysis including all cross-tomographic combinations shift the parameter by
at most $0.3\sigma$ towards higher $S_{8}$ values, even when the inferred
$A_{\rm IA}$ is as large as unity. Although their analysis is based on the
different statistics (they use persistent homology instead of peak count),
their results should hold here too, given that peaks are a subset of their
data vectors. Therefore residual IA cannot play an important role in the
DES-Y1 peak count analysis, justifying our choice to include the cross-
redshift bins here (up-to-pairs, but not the triplets nor the quadruplets
since these are not fully modelled yet for the DES-Y1 data).

A fourth difference in the likelihood concerns the treatment of the baryonic
feedback: in HD21 the peak statistics are measured in the Magneticum
simulations to ensure that the selected elements from the data vectors are
immune to unmodelled baryonic mechanisms, but no marginalisation is included.
This can potentially lead to a slightly over-optimistic precision on the DES
likelihood compared to the KiDS-1000 likelihood, which includes
marginalisation over the $b_{\rm bary}$ parameter. We therefore decided to
include in our joint analysis the same marginalisation machinery for both the
KiDS-1000 and the DES-Y1 pipelines. Moreover, we use a unique $b_{\rm bary}$
parameter to infuse baryonic feedback into both surveys, since these physical
processes describe physics that affect the foreground matter distribution
independently of survey-specific source selection. In total, the combined-
survey analysis marginalises over nine redshift bias parameters, nine shear
bias parameters, one IA and one baryon parameter. The sampling strategy of the
DES-related parameters are listed in Table 3. Finally, given the absence of
overlap between two survey footprints and the compatibility of the priors, the
two likelihoods can be directly added at each evaluation point, without
needing to consider cross-survey covariance.

##  4 Systematic Uncertainties

As the amount of high quality lensing data keeps increasing, the statistical
precision reaches unprecedented high levels, and consequently understanding
and controlling the residual systematics in every segment of the data analyses
has become one of the primary objectives and focus of development in the field
of beyond-2pt statistics. We investigate here a number of such systematic
effects that have been identified in the literature and mentioned earlier,
including residual uncertainties related to interpolation in the modelling,
shear calibration, photometric redshifts, astrophysics (intrinsic alignments
and baryonic feedback), simulation-based covariance matrix, non-linear
physics, source-lens coupling and likelihood sampling strategies. Some of the
systematics that are ignored in the current work are those related to the
effect of source blending, depth variations, PSF leakage, or the cosmology
dependence of the IA signal. These will likely become important in the future,
but can be safely omitted in current Stage-III lensing surveys (see HD21,
Zürcher et al., 2022). Amongst those that we investigate here, many are shown
to be sub-dominant or heavily suppressed by our range of
$\mathcal{S}/\mathcal{N}$, while others are forward-modelled with nuisance
parameters that are marginalised over in the likelihood analysis.

###  4.1 Modelling

Figure 5: Accuracy of the KiDS-1000 GPR emulator, computed with a leave-one-
out cross-validation test. The results are colour-coded with the $S_{8}$ value
of the removed training point, and compared with the statistical precision on
the measurement of the peak function (shown with the black dashed lines). The
black solid line indicates the $\Lambda$CDM node, and the different panels
show the auto- and cross-redshift (up to pairs) measurements; the other 15
tomographic combinations show a similar precision and hence are not shown. The
outliers seen in a few panels are of extreme $S_{8}$ values and as such
required the test emulator to extrapolate; this does not occur with the full
emulator, and should therefore not be considered when estimating the
interpolation error.

The accuracy of the peak function modelling has two aspects to be considered:
we must understand both how well the cosmology scaling is captured by the
emulator, and whether any elements of our data vector are affected by either
resolution limits of our simulations or the choice of gravity solver.
Regarding the first aspect, the cosmo-SLICS have been shown to match in
precision the commonly used CosmicEmu (Heitmann et al., 2014), and to even
outperform the latter in terms of range, benefiting from more training nodes
(Harnois-Déraps et al., 2019). The GPR interpolation uncertainty is fully
propagated in the likelihood, and is quantified with a leave-one-out cross-
validation test, in which the emulator is trained on all but one node,
producing predictions at that node that are then compared with the
measurement. As shown in Fig. 5, we cycle over all 26 nodes in this way and
estimate an upper limit on the error (since the actual GPR has all the nodes).
The accuracy degrades towards large positive or negative
$\mathcal{S}/\mathcal{N}$ values, but most of the cross-validation lines lie
well within the statistical precision on the data, shown with the dashed black
lines. There are a handful of exceptions with poorer accuracy, attributed to
removing extreme values of $S_{8}$ from the training set and therefore
effectively demanding the GPR to extrapolate. With these edge nodes included,
the full emulator has no such outliers. In fact, as argued in HD21, the most
reliable estimate of the emulator’s precision is evaluated by removing the
fiducial cosmology and training on the others, which is shown as the thick
black lines in the figure and always well within the statistical precision.
For the KiDS modelling, the interpolation error is mostly at the 1% error over
the range $-1.0<\mathcal{S}/\mathcal{N}<3.0$ (our ‘clean’ range), and is
otherwise always under 10%. Similarly, the DES interpolation error is
everywhere under 2% (see HD21), the difference coming from the choice of
$\mathcal{S}/\mathcal{N}$ cuts. This will likely be a limiting factor for
future data analysis with sub-percent accuracy requirement on the modelling,
and will be addressed by increasing the number of nodes in the next generation
of the Cosmology Training Set. At the moment however, the interpolation error
is low enough for our analysis. We nevertheless include it in our error budget
by averaging over the square of the residuals (after the outliers have been
removed):

$\displaystyle{\rm Cov}_{\rm interp}={\rm diag}\left<\left(N_{\rm peaks}^{\rm
GPR}-N_{\rm peaks}^{\rm sim}\right)^{2}\right>\,,$ (7)

and adding this to our statistical covariance. We could have instead used the
errors directly provided by the Gaussian Process Emulator, however Heydenreich
et al. (2021) have shown that the two methods yield posteriors with negligible
differences. This contribution, although small, helps with the goodness-of-fit
in the data analysis.

The second aspect, concerning the training sets themselves, is discussed below
in the section on $N$-body resolution.

###  4.2 Shape calibration

Table 1 shows the average multiplicative correction factors $m_{a}$ that must
be applied to the observed galaxy shapes in order to correct for a combination
of residual PSF leakage, blending and measurement noise, as assessed from
Giblin et al. (2021). While in A21 the uncertainty on the shape calibration is
absorbed directly in the analytical covariance matrix, our simulation-based
method works instead at the level of the data vector, as for all other
nuisance parameters. The $M_{\rm ap}$ estimator itself is unbiased (see Eq.
2), however we must propagate forward the uncertainty on the $m_{a}$
calibration. The impact of potentially mis-calibrated shape measurements is
estimated by infusing a non-corrected global term $m_{a}\rightarrow
m_{a}+\Delta m_{a}$ directly in the simulations and measuring the effect on
the different elements of the peak function $N_{\rm peaks}^{\kappa}$. As we
show later, this systematic effect is completely subdominant compared to the
others due to the tight priors on $\Delta m_{a}$ (reported in Table 2), and
hence it is sufficient to model its impact with a reduced accuracy. In HD21
the estimation is based on a linear regression (i.e. $\partial N_{\rm
peaks}/\partial\Delta m$ per data element) that is fit through 10 values of
$\Delta m_{a}$. We use here only two points, at $\pm 1\sigma$, which is
sufficient given the small values of $\Delta m_{a}$. The measured $\partial
N_{\rm peaks}/\partial\Delta m$ is further discussed in Appendix B, and is
used to modify the data vector for any value of $\Delta m_{a}$ (see Eq. 8)
sampled in the likelihood. For cross-redshift tomographic bins, we use the
mean shift, e.g. $\Delta m_{a}^{1\cup 2}=(\Delta m_{a}^{1}+\Delta
m_{a}^{2})/2$, which is consistent with what is currently done for all shear
two-point function analyses. We could instead use an $n_{\rm gal}$-weighted
mean to compute the $\Delta m_{a}$ shift in cross-redshift tomographic bin,
however this should have a negligible effect given the tight priors on these
parameters, and we therefore leave this for the future.

###  4.3 Photometric redshifts

The KiDS-1000 uncertainty on the redshift distributions has been fully
quantified in Hildebrandt et al. (2021), where it is shown that the mean of
the $n(z)$ is captured to a high accuracy, varying by no more than 0.014 at
the $1\sigma$ level151515The full shape of the $n(z)$ is less accurate than
its mean, and which consequences we leave for future work. . The posteriors on
the mean of the redshift distributions are used as priors on nuisance
parameters in the current work, summarised in Table 2. In this case however,
the five redshift bias parameters $\Delta z_{a}$ must be drawn from a
correlated distribution. This is achieved in a two-step operation where we
first draw five uncorrelated numbers from the priors, then rotate into the
correlated space using a Cholesky decomposition of the redshift covariance
matrix:

$\displaystyle C_{z}\times
10^{5}=\begin{bmatrix}11.20&2.600&1.562&0.056&0.622\\\
2.600&12.78&4.081&-1.692&-0.2140\\\ 1.562&4.081&13.81&-1.139&0.525\\\
0.056&-1.692&-1.139&7.551&3.054\\\
0.622&-0.2140&0.525&3.054&9.496\end{bmatrix}\,,$

which results in a correlated sampling of these five nuisance parameters (see
A21, Hildebrandt et al., 2021, for more details). We produced a dedicated set
of Redshift Training Set simulations in which the $n(z)$ are shifted, but
which are otherwise identical to the Cosmology Training Set at the fiducial
cosmology. Following HD21, we measure the peak function on full mock surveys
with 10 shifts, each with a slightly different value of $\Delta z_{a}$ sampled
from the prior, then extract a linear fit per data element and estimate
$\partial N_{\rm peaks}/\partial\Delta z_{a}$. This derivative is used to
forward model redshift uncertainties on our data vector for arbitrary $\Delta
z_{a}$ values. Again, we use the mean shift when considering cross-redshift
bins, and the DES-Y1 $\partial N_{\rm peaks}/\partial\Delta z_{a}$
measurements from HD21.

###  4.4 Astrophysics

Cosmic shear measurements are strongly affected by IA and baryon feedback.
Using the IA and the Baryons Training Sets described in Sec. 2.3, we estimate
in a similar way $\partial N_{\rm peaks}/\partial A_{\rm IA}$ and $\partial
N_{\rm peaks}/\partial b_{\rm bary}$, where $A_{\rm IA}$ and $b_{\rm bary}$
are free parameters that control the levels of IA and baryon contamination,
respectively. The IA derivative is obtained by linear fitting the peak
function’s response to changes in $A_{\rm IA}$, measured from the IA Training
Set infused with $A_{\rm IA}$ = 2.0,1.0, 0.0, $-$1.0 and $-$2.0. Since IA is
currently not well constrained and the NLA parameterisation is an effective
model, we adopt a wide top-hat prior over the range [$-$5.0 ; 5.0], as argued
in Joachimi et al. (2021). This extrapolates our fit to larger $A_{\rm IA}$
values, which can in principle become inaccurate, however in the end the
$3\sigma$ region of our posterior is fully contained within the training range
(see Sec. 5). Similarly, the baryon derivative is measured from the Baryons
Training Set, which we use to infuse a baryonic correction whose strength is
controlled by the parameter $b_{\rm bary}$. The case $b_{\rm bary}=0.0$
corresponds to a dark matter-only universe, while $b_{\rm bary}=1.0$
corresponds to the case where the feedback processes is exactly described by
the Magneticum physics. There is a large uncertainty on the amplitude of this
baryon correction, hence we scale the measured baryonic correction with a free
parameter $b_{\rm bary}$. Since the Magneticum suites are already a strong
model (see Martinet et al., 2021, for a comparison with other hydrodynamical
simulations), we sample the range $b_{\rm bary}\in[0.0,2.0]$, thereby spanning
a variety of realistic models (albeit imposing a fixed shape for the relative
signal). As seen later, low $b_{\rm bary}$ values are not well constrained by
the data while larger values are strongly disfavoured, hence we do not extend
the prior limit beyond 2.0.

###  4.5 Implementation of forward-modelled systematics

Four sources of systematics are forward-modelled in our pipeline. Following
Heydenreich et al. (2022), we construct systematics-infused data vector as:

$\displaystyle N_{\rm peaks}^{\rm syst}({\boldsymbol{\pi}},\Delta m_{a},\Delta
z_{a},A_{\rm IA},b_{\rm bary})=$ $\displaystyle\hskip 36.98866ptN_{\rm
peaks}^{\rm GPR}({\boldsymbol{\pi}})+\left[\partial N_{\rm
peaks}/\partial\Delta m_{a}\right]\Delta m_{a}+\left[\partial N_{\rm
peaks}/\partial\Delta z_{a}\right]\Delta z_{a}\,\,...$ $\displaystyle\hskip
56.9055pt+\left[\partial N_{\rm peaks}/\partial A_{\rm IA}\right]A_{\rm
IA}+\left[\partial N_{\rm peaks}/\partial b_{\rm bary}\right]b_{\rm bary}\,,$
(8)

where the twelve parameters ($\Delta m_{a},\Delta z_{a},A_{\rm IA},b_{\rm
bary}$) are sampled from the priors described in Table 2. We marginalise over
these nuisance parameters when inferring the values of the cosmological
parameters. Equation (8) assumes that these different systematics are
independent of cosmology and from each other, which we know is not entirely
true. It has been shown that the cosmology dependence of the baryon feedback
is a second order effect (McCarthy et al., 2017), supporting our simplified
approach, however the intrinsic alignments couple to the tidal field that is
in itself cosmology dependent. The shear calibration and redshift errors are
independent of cosmology a priori, however the derivatives of the peak
function with respect to $\Delta m_{a}$ and $\Delta z_{a}$ are not (see HD21),
a secondary effect we neglect here. Moreover, it has been shown that the
photometric and shape calibration errors are sometimes correlated (MacCrann et
al., 2022). Although these approximation will become important in Stage-IV
surveys, the current level of statistical precision allows us to relax the
modelling of these effects without hurting our results. We illustrate this
point in Section 5.2 by running inference MCMC chains in which the modelling
of some or all of these systematic effects are switched off: the minor impact
this has on the inference validates this approach. We also assume here that
these systematic effects have a linear dependence on the nuisance parameter,
which is probably not entirely true, but has been shown to be good enough for
Stage-III lensing data in Heydenreich et al. (2022, see their figure 7).

###  4.6 Other sources of systematics

In addition to the main systematic effects described in the last section, we
consider here other known sources of errors that could potentially impact our
results.

####  $N$-body resolution

Being completely simulation-based, our analysis relies on the quality of the
underlying training samples. As mentioned already in Sec. 2.3, the Cosmology
Training Set has been shown to closely reproduce the non-linear clustering of
the Cosmic Emulator (Heitmann et al., 2014), which is based on a completely
independent $N$-body code. This agreement between different gravity solvers is
key to assert the accuracy of the non-linear solution to structure formation
(see, e.g. Euclid Collaboration: Knabenhans et al., 2019, for a comparison
between different $N$-body solvers), and the convergence of the solution must
be assessed via a comparison with calculations carried out with a higher
force/mass resolution simulations. As shown in HD21, known limits in the mass
resolution of the cosmo-SLICS used for the peak function emulation mainly
affect high peaks. More precisely, $N^{\kappa}_{\rm
peaks}(\mathcal{S}/\mathcal{N}>4.0)$ is systematically under-predicted by tens
of percent, while the $\mathcal{S}/\mathcal{N}=4.0$ count is affected by no
more than 5%. This is in fact one of the main justification for our initial
choice of upper $\mathcal{S}/\mathcal{N}$ limit.

The KiDS-1000 data are deeper than DES-Y1, and hence the sensitivity to such
non-linear effects could be accrued here. We verify this by running our
cosmological inference on the peak count statistics measured from the SLICS-
HR, in which the increased force resolution results in a slightly larger
number of large positive and negative peaks. Details are presented in Appendix
B, but in short our data selection and marginalisation scheme almost
completely protects us against this, yielding no noticeable shifts on
$\Omega_{\rm m}$ nor $S_{8}$. As in HD21, we nevertheless compute a
multiplicative factor from the ratio between the SLICS-HR and the mean of the
SLICS and optionally apply it on our model predictions during the likelihood
sampling. The overall effect is smaller than the baryon and IA corrections,
hence marginalising over these latter two significantly washes out the impact
of inaccurately-modelled non-linear physics under question here. In the future
we intend to look into multi-fidelity emulators as in Ho et al. (2022). The
T17 and Magneticum simulations were produced with a different $N$-body solver,
and we show later that their reduced spatial resolution can affect quite
significantly the peak count statistics. We treat this as a further
uncertainty on the small scale physics and optionally add their scatter to the
theoretical error in the covariance matrix, similar to Eq. (7). In two-point
statistics analysis, this would be equivalent to including a theoretical error
in the covariance matrix to account for difference between the $P(k)$
predictions provided by HaloFIT (Takahashi et al., 2012), HMCode (Mead et al.,
2016) or the baccoEMU (Angulo et al., 2021), which can have a significant
impact on the results (Aricò et al., 2023).

#### Ray-tracing approximations

Our ray-tracing method in itself contains approximations and algorithmic
components that are bound to affect to some level the lensing statistics. Most
importantly, the finite thickness of the mass sheets and the randomisation
process between them destroys correlations along the line of sight; in
particular it can slice large galaxy clusters in two, and no structures larger
than 257.5 $h^{-1}$Mpc can exist along the line-of-sight in our light-cones
(except for the T17 mocks, which we discuss below). This suppresses some of
the large-scale power, as documented in Takahashi et al. (2017, see their
Appendix B). However, smaller structures, such as those probed by the peak
statistics, are left completely unaffected by this, which is why no forward
modelling is needed here. This has been measured specifically for peak
statistics in Zorrilla Matilla et al. (2020) where it was found to play a
subdominant role even for Stage-IV surveys. Of course, full on-the-fly light-
cones such as the ‘Onion Universe’ methods (Fosalba et al., 2008) avoid these
problems, which we will consider for future analyses. The T17 simulations have
thinner mass shells of 150.0 $h^{-1}$Mpc, but they are constructed such that
the structures are preserved in groups of three shells, thus yielding a
coherence length of 450.0 $h^{-1}$Mpc, further suppressing this residual
systematic effect.

Another source of error comes from the fact that our simulations assume the
Born approximation in the flat-sky limit, which introduces small inaccuracies
at high-$\ell$ and low redshift, respectively (Hilbert et al., 2020). However,
these are affecting the signal at a level much smaller than the statistical
accuracy of our lensing data, and are not expected to matter here.

#### Covariance matrix

Estimation of the covariance matrix is one of the main computational
challenges for non-Gaussian weak lensing probes, as it requires a large number
of simulations with a resolution that is high enough to capture the non-linear
physics being measured. Resorting to approximate methods such a Flask (Xavier
et al., 2016) and ICE-COLA (Izard et al., 2018) can significantly lower the
computational cost of creating such mocks, but at the price of a reduced
precision on the physics under investigation. We instead opted for mocks
produced by a full $N$-body suite, our Covariance Training Set, and are
therefore only limited by the number of mocks and their box size. To test the
convergence of our covariance matrix with respect to $N_{\rm sims}$, we run an
inference analysis in which we increase the number of pseudo-independent
realisations to 2120 (and adjusted the likelihood $N_{\rm sims}$ parameter
accordingly), and find an excellent match to the posterior, with only the tail
of the distribution being slightly modified. We could also have opted for a
data compression such as in Zürcher et al. (2022) but that is not necessary
given our results have converged, and our choice of likelihood accounts for
the noise in the covariance matrix.

The simulation box size could also affect our results, however it has been
shown in Harnois-Déraps et al. (2019) that the SLICS contains about 75% of the
“super-sample covariance” term (SSC), when applied to 2pt statistics, yielding
constraints on cosmological parameters that are highly accurate. Although this
has not been demonstrated to date, peak count statistics are thought to be
even less affected by the SSC, given that the covariance is close to being
Poissonian, not Gaussian. As such, it scales with the number of peaks
measured, which is independent of the survey window. In addition, as mentioned
earlier, Burger et al. (2022) finds for the density-split statistics an
excellent agreement between the SLICS covariance and that from full sky log-
normal Flask mocks (which contain an incomplete contribution from the
trispectrum term but the full SSC), supporting our claim that the partly
missing SSC must have a minimal influence on our error budget. This is also
consistent with the recent findings from Linke et al. (2023) according to
which the SSC term affects only the Fourier space estimators, whereas
covariance matrix measured from intra-survey real-space statistics such as the
$M_{\rm ap}$ are unbiased.

#### Source-lens coupling and blending

An important difference between real and mock galaxies is that those in the
data are clustered, which leads to a number of effects that are systematically
absent from the calibration sample. For example, the quality of the shape
measurements is lowered in regions of high density due to blending and
obscuration. More importantly, the uncertainty in photometric redshifts is
particularly severe in such areas, which often results in cluster members
being wrongly assigned a higher redshift. This subsequently creates a small
population of apparently high-redshift outliers that carry an unexpectedly
weak shear component, thus diluting the overall lensing signal. Correcting for
this can be partially achieved with ‘boost factors’, however it was shown in
HD21 and Zürcher et al. (2022) that even though the excess clustering around
high peaks is indeed measured in the data, the impact this has on the inferred
cosmology can be safely ignored. It was also shown in Gatti et al. (2023a)
that source clustering had a minimal effect on the peak count statistics,
supporting our choice to neglect this here.

#### Sampling the likelihood

Our likelihood sampling strategy, described in Sec. 3.3, assumes a flat prior
for the four main cosmological parameters and the two astrophysical parameters
($A_{\rm IA}$ and $b_{\rm bary}$), and Gaussian priors for the parameters
associated with photometric redshifts and shape calibration. This is not
strictly speaking a non-informative approach, however the prior edges about
the key measured parameters are sufficiently broad to have negligible impact
on the posterior. Since it is found in Lemos et al. (2022) that MultiNest
tends to yield slightly over-precise constraints, we use the Nautilus sampler
for our fiducial results, but report both.

We also note that our cosmology sampling strategy is different from the other
KiDS-1000 cosmic shear analyses, mainly due to the volume where our emulator
is valid. For example, A21 sample uniformly the parameters $S_{8}$,
$\omega_{\rm c}\equiv\Omega_{\rm c}h^{2}$, $\omega_{\rm b}\equiv\Omega_{\rm
b}h^{2}$, $h$ and $n_{\rm s}$. This choice is designed to avoid regions of
parameter space that are strongly disfavoured by external data, and it was
shown in Joachimi et al. (2021) that while it disfavoured high $\Omega_{\rm
m}$ values already in prior space, the resulting $S_{8}$ prior space is highly
uninformative. We could have taken a similar approach, however our emulator is
much quicker, and hence it is more natural to sample the full training space,
ensuring a wide sampling of $\Omega_{\rm m}$, $S_{8}$ and $w_{0}$.

Another aspect that currently limits our sampling strategy is the fact that we
hold the value of many parameters fixed, notably $\Omega_{\rm b}$ and $n_{\rm
s}$. In contrast, the DES-Y3 peak count analysis of Zürcher et al. (2022) use
derivatives to marginalise over variation in these parameters, following the
approach we adopt for IA and baryons. Neglecting to account for these has a
small effect on current data sets (the DES-Y3 joint peaks+power spectrum
analysis finds to be of about 0.13$\sigma$), which are thus ignored here.

####  $M_{\times}$ modes

The observed weak lensing signal can generally be decomposed into a
combination of $E$\- and $B$-modes, the latter of which can be estimated for
any measurement by rotating all galaxies by 45 degrees; therefore, for the
aperture mass map statistics, it is often referred to as
$M_{\times}({\boldsymbol{\theta}})$. The cosmic shear signal being a pure
$E$-mode generator to first order, measurements of $B$-modes are therefore
routinely used to assess the presence of residual systematics in lensing data
(see, e.g. Zürcher et al., 2022, for a recent application to peak statistics).
Whereas the two-point function $B$-mode signal is zero in absence of
systematics, the construction of aperture mass maps on a grid inevitably
injects non-zero $M_{\times}$-modes due to the missing contribution from sub-
pixel scales (Kilbinger et al., 2006). This can be important: for a small-
angle cut-off scale of 10 arcsec and an aperture of $\theta$=2.0 arcmin,
$B$-modes measured this way can reach about 10 percent the size of the
$E$-mode $M_{\rm ap}^{2}$ signal. This effect is accentuated for larger cut-
off scales and smaller opening angles $\theta$. Given our pixel scale of 35
arcsec, we do expect non-zero $M_{\times}$-modes to be introduced by our
aperture map making, which we fully quantify in Appendix A. We show therein
that the level of contamination is consistent with noise, that there is no
evidence for residual systematics in the data from this measurement, hence
that our cosmological analysis is clean of $B$-modes.

###  4.7 Peak count to cosmology pipeline validation

Figure 6: Full inference analysis on the Validation Set (SLICS-HR) peak count
data (blue) with Multinest, optionally infused with intrinsic alignments (red)
or baryon feedback (grey). We marginalise over these two effects plus shape
calibration and photometric uncertainty. The priors are shown by the dashed
grey lines at the edge of the panels, while the cross-hairs show the input
truths.

We test our KiDS-1000 cosmology inference pipeline by analysing simulated data
of known cosmology, infused with a controlled amount of residual systematics.
In order to avoid confirmation bias, these tests are carried out with the
Validation Set, which have not been used in the cosmology training nor for the
covariance estimation, with an $N$-body force resolution that is higher than
the other simulations used in this work161616We have further verified that the
$w$CDM cosmology is correctly inferred when analysing data from the Cosmology
Training Set but these tests are easier to satisfy since the data is used for
training the emulator. We discuss these in greater details in Appendix B.. In
addition, we use the forward-modelling approach presented in Sec. 4 to infuse
the simulated data vectors with either intrinsic alignments (assuming $A_{\rm
IA}=1.0$) or baryonic feedback (with $b_{\rm bary}=1.0$). Fig. 6 shows the
results for these three analysis cases. The maxima of the projected posterior
distributions are all centred on the input truth, except for the $b_{\rm
bary}$ parameter, which are away from zero even in the no-baryon cases. This
is a projection effect similar to those discussed in Joachimi et al. (2021),
Chintalapati et al. (2022) and Dark Energy Survey and Kilo-Degree Survey
Collaboration (2023), and we have verified that reducing the lower prior limit
to $b_{\rm bary}=-2.0$ pushed both the red and blue maxima towards the ground
truth.

In this test, there is a secondary solution for $\Omega_{\rm m}\sim 0.4$ that
is unexpected, and not observed in other peak count analyses (Martinet et al.,
2018; Zürcher et al., 2022; Marques et al., 2023). As detailed in Appendix B,
this feature persists when analysing data from the Cosmology Training Set at
the fiducial cosmology, from the Baryons Training Set and from the T17 mocks,
but can vanish at other cosmologies. This is caused first by the poor
sensitivity of the current lensing data to $\Omega_{\rm m}$, as also seen in
the large $\Omega_{\rm m}$ scatter reported in A21 between different two-point
functions, but more importantly by limits in our GPR emulator, whose residual
inaccuracy mostly affects this parameter. Inferring $\Omega_{\rm m}$ from
single mock survey realisations (as opposed to a mean over several light-cones
or shape noise realisations) yields posteriors drawn from either one of these
peaks, resulting in occasional strong biases on this cosmological parameter.
We thoroughly verify that only $\Omega_{\rm m}$ is affected by this, and
therefore do not report its value in our main analyses. See Appendix B for
full details.

However, the secondary $\Omega_{\rm m}$ solution corresponds to an $S_{8}$
posterior that is slightly lower than the main solution, which means that if a
particular realisation of the data prefers this region, it will on average
have an $S_{8}$ value about 0.03 lower, which is of the size of our
statistical precision. Conversely, realisations that prefer lower $\Omega_{\rm
m}$ tend to have $S_{8}$ values that are 0.02 higher than the input truth. We
further observe that this is not always the case: some individual mock survey
realisations from the Covariance Training Set have a best fit $\Omega_{\rm
m}\sim 0.45$, yet their $S_{8}$ is unbiased compared to the input truth. Given
that this $0.02-0.03$ shift is about a $1\sigma$ shift, this potentially
dominates the systematic error budget on $S_{8}$, which we therefore must
report as $\sigma$(syst)$=^{+0.03}_{-0.02}$.

This additional systematics error take its roots from the tilt in the
$[S_{8}-\Omega_{\rm m}]$ posterior, which indicates residual correlation
between these two parameters. We can suppress this tilt, and hence the
additional error, by replacing $S_{8}$ with
$\Sigma_{8}^{\alpha}\equiv\sigma_{8}[\Omega_{\rm m}/0.3]^{\alpha}$, where
$\alpha$ is the parameter that best fits the $[\Omega_{\rm m}-\sigma_{8}]$
degeneracy. According to this metric, $\Sigma_{8}^{\alpha}$ is the most
robustly measured quantity from peak statistics, with no need for a standalone
$\sigma$(syst) term, in this case a significant advantage. With the validation
data, we find $\Sigma^{\alpha}_{8}=0.824^{+0.033}_{-0.033}$, with
$\alpha=0.582$, in excellent agreement with the input truth of 0.811 with the
same $\alpha$. We report the measurements of both $S_{8}$ and
$\Sigma^{\alpha}_{8}$ in this paper, but while emphasise is on the former to
better compare with previous measurements from the literature, the latter is
more robust and has interesting properties which we highlight as well, notably
on increasing the agreement with previous KiDS-1000 measurements and lowering
the tension with external probes.

Back to Fig. 6, we observe that the posteriors on $w_{0}$ and $b_{\rm bary}$
are wide and significantly overlap with the prior limits, and we thus expect
to be unable to place meaningful constraints on these parameters with our main
KiDS-1000 analysis alone. We observe a degeneracy in the $[S_{8}-w_{0}]$ plane
here, however we show in Appendix B that it is not always seen when analysing
other cosmologies, making it impossible to draw physically meaningful
conclusions about this. Only the $[S_{8}-A_{\rm IA}]$ plane is well
constrained with the current KiDS-1000 peak count analysis: we achieve a 4.4%
precision measurement on $S_{8}$, with $S_{8}^{\rm SLICS-
HR}=0.816^{+0.039}_{-0.033}$ (truth is 0.813), and a precision of
$\sigma_{A_{\rm IA}}=0.45$ on $A_{\rm IA}$, sampling the likelihood with
Multinest.

The DES-Y1 pipeline validation is presented in HD21, while that for the joint
KiDS-DES is presented in Appendix B, showing again an excellent agreement
between the inferred cosmology and the input truth.

##  5 Results: KiDS-1000

Figure 7: KiDS-1000 constraints on the two best-measured parameters from peaks
count statistics, for different selection of redshift bins. Tomographic
analyses all break the $S_{8}-A_{\rm IA}$ degeneracy.

We present in this section the results from our cosmological inference
analyses, beginning with the fiducial KiDS-1000 pipeline, then reporting on
the importance of various selection cuts and systematic effects. For reasons
explained in Sec. 4.7, we report only the constraints on $S_{8}$ and $A_{\rm
IA}$; results are summarised in Table 4 and further condensed in Fig. 10.

From our fiducial full tomographic KiDS-1000 analysis of the measurements
presented in Fig. 3, we obtain:

$\displaystyle S_{8}^{\rm KiDS}=0.733^{+0.032+0.020\,(\rm
syst)}_{-0.032-0.030\,(\rm syst)}\,,\hskip 28.45274ptA_{\rm
IA}=0.71^{+0.49}_{-0.49}\,,$ (9) $\displaystyle\Sigma_{8}^{\rm
KiDS}=0.765^{+0.030}_{-0.030}\,,\hskip 28.45274pt\alpha=0.600\,,$ (10)

after marginalising over three cosmological parameters ($\Omega_{\rm m}$,
$w_{0}$ and $h$) and 11 nuisance parameters ($5\times\Delta m_{a}$,
$5\times\Delta z_{a}$, $b_{\rm bary}$). Unless explicitly mentioned, all
quoted parameter constraints correspond to the mean $\pm 1\sigma$ region of
the marginalised posterior, not to be confused with the point of maximum
likelihood in the higher-dimensional space. This is therefore a 3.9%
measurement of the structure growth parameter $\Sigma_{8}$. The best-fit model
is shown with the red line in Fig. 3. The joint constraints on two of these
parameters are shown in Fig. 7, along with results from different selections
of tomographic bins. Importantly, the strong $S_{8}-A_{\rm IA}$ degeneracy
seen in the no-tomographic case (the tilted dark purple contour) is lifted by
tomographic decompositions, which capture the different redshift dependence of
the cosmological and IA signals. This is an important verification of our IA
modelling, which, under the NLA framework, compensates large $S_{8}$ values
with large tidal alignments in the no-tomographic case only. Notably, all
tomographic additions contribute to further tightening the constraints, once
again demonstrating the power of using cross-redshift bins in non-Gaussian
statistics. We also observe that all cases shown in Fig. 7 are consistent,
providing statistical robustness to our measurement.

At our best-fit parameters the measurement yields a $\chi^{2}$ of 250, which
reduces to $\chi^{2}_{\rm red}$ = 1.22 after dividing by $\nu=(220-4.5)=205.5$
degrees of freedom. Note that although we use six unconstrained
parameters171717We do not count as free parameters those nuisance parameters
for which we impose a tight prior. in our likelihood evaluation (the four
cosmological parameters plus $A_{\rm IA}$ and $b_{\rm bary}$), it was shown in
Joachimi et al. (2021) that an effective number of $\nu=4.5$ free parameters
better describes the weak lensing data given the existing correlations and
degeneracies, results which we have used here181818It is not guaranteed that
the exact same effective number of degrees of freedom applies here, given that
the likelihood is not sampled over the same volume. We have checked that our
goodness-of-fit is robust over choices for this quantity, with PTE varying
between 0.48 and 0.37 over the range $2<\nu<7$.. Our PTE for this measurement
is 0.43, which is well above our threshold of 0.01, using the non-$\chi^{2}$
distribution described in Appendix C. It is worth noting that the KiDS-1000
shear two-point correlation functions and band power analyses had a lower
goodness-of-fit, with PTE=0.034 and 0.013, respectively.

Figure 8: Effect of $\mathcal{S}/\mathcal{N}$ cuts on the KiDS-1000
constraints. Large peaks ($\mathcal{S}/\mathcal{N}$ $>$ 3) slightly increase
the statistical precision on $S_{8}$, as seen by comparing the red and blue
contours (see values in Table 4). Negative peaks ($\mathcal{S}/\mathcal{N}$
$<$ -1), included in the red but excluded from the grey contours, help in
breaking the $[S_{8}-A_{\rm IA}]$ degeneracy. The grey and blue contours
correspond to the ‘clean’ and ‘aggressive’ cases, respectively.

In many previous analyses, sampling and marginalisation over $w_{0}$ is often
excluded, being considered an extension to the vanilla $\Lambda$CDM scenario.
In the present case, fixing $w_{0}$ to $-1.0$ when sampling the
likelihood191919This still uses the same $w$CDM emulator, but only varying the
other three cosmological parameters. results in minor changes to the reported
$(S_{8},\Sigma_{8},A_{\rm IA})$ constraints, leading to
$0.732^{+0.029}_{-0.029}$, $0.767^{+0.026}_{-0.026}$ and
$0.73^{+0.48}_{-0.48}$. Interestingly, we find that the impact of opening up
the $w_{0}$ dimension is far lower than for the two-point statistics, where
Tröster et al. (2021) finds a degradation by a factor of a few on the $S_{8}$
constraints (compare their figures 1 and 6). Different degeneracy-breaking
directions are likely causing this difference, which is promising for upcoming
measurements of $w_{0}$ with alternative statistics (see Martinet et al.,
2021, for a Stage-IV lensing forecast on the dark energy parameter with peak
statistics).

One of the key questions to be explored by beyond-2pt statistics concerns the
exact origin of the non-Gaussian cosmological information. Large peaks are
often associated with massive galaxy clusters, which are known to be highly
sensitive to the dark energy equation-of-state parameters for instance,
however the wide projection effect and the fact that baryons, IA and non-
linear physics maximally affect these large $\mathcal{S}/\mathcal{N}$ peaks
(Martinet et al., 2021; Harnois-Déraps et al., 2021b) complicate the picture.
To (partly) answer this question, we investigate the constraining power
contained in the highest ($\mathcal{S}/\mathcal{N}$$>$3) and lowest
($\mathcal{S}/\mathcal{N}$$<$0) bins by removing these sequentially from the
‘aggressive’ data vector ($-2.5\leq\mathcal{S}/\mathcal{N}\leq 4.0$) in the
likelihood. The results are shown in Fig. 8, where we observe that the
negative $\mathcal{S}/\mathcal{N}$ peaks significantly help break the
$[S_{8}-A_{\rm IA}]$ degeneracy, while the highest peaks help in tightening
the $S_{8}$ constraints. In an analysis that ignored the role of IA, M20 found
that the amount of information about $S_{8}$ that is contained in negative
peaks is quite small, however here we find that they actually play a key role
once IA are forward-modelled.

Table 4: Summary of our cosmological inference analyses. Posteriors on $\Omega_{\rm m}$, $h$ and $w_{0}$ are prior-limited, so their constraints are not reported here. Unless explicitly specified in the first column, the KiDS-1000 measurements are based on the ‘clean’ data vector, i.e. $-1.0<\mathcal{S}/\mathcal{N}$$\leq 3.0$. The last column presents the maximum a posteriori (MAP) values. Validation of the inference pipelines on mock data are presented in Appendix B.  Peak count analysis of KiDS-1000 |
---|---
| Nautilus | Multinest | MAP
| $S_{8}$ | $A_{\rm IA}$ | $S_{8}$ | $A_{\rm IA}$ | $S_{8}$
Fiducial | ${0.733^{+0.032}_{-0.032}}$ | ${0.71^{+0.49}_{-0.49}}$ | ${0.733^{+0.021}_{-0.027}}$ | ${0.74^{+0.43}_{-0.43}}$ | 0.726
$\Lambda$CDM | $0.732^{+0.029}_{-0.029}$ | $0.73^{+0.48}_{-0.48}$ | ${0.729^{+0.026}_{-0.026}}$ | ${0.73^{+0.43}_{-0.43}}$ | 0.718
Auto-only | $0.734^{+0.050}_{-0.050}$ | $-0.2^{+1.0}_{-1.0}$ | ${0.732^{+0.039}_{-0.048}}$ | ${-0.2^{+1.0}_{-1.0}}$ | 0.725
Up to pairs | $0.750^{+0.036}_{-0.050}$ | $0.10^{+0.69}_{-0.69}$ | ${0.742^{+0.032}_{-0.043}}$ | ${0.11^{+0.63}_{-0.63}}$ | 0.742
Up to triplets | $0.740^{+0.035}_{-0.035}$ | $0.35^{+0.53}_{-0.53}$ | ${0.740^{+0.029}_{-0.029}}$ | ${0.37^{+0.48}_{-0.48}}$ | 0.738
No tomo | $0.695^{+0.033}_{-0.087}$ | $-0.6^{+1.7}_{-2.2}$ | ${0.690^{+0.038}_{-0.068}}$ | ${-0.7^{+1.5}_{-2.0}}$ | 0.682
-2.5$<\\!\mathcal{S}/\mathcal{N}\\!\leq$4.0 | $0.720^{+0.036}_{-0.026}$ | $0.08^{+0.30}_{-0.30}$ | ${0.717^{+0.031}_{-0.022}}$ | ${0.07^{+0.27}_{-0.27}}$ | 0.728
-2.5$<\\!\mathcal{S}/\mathcal{N}\\!\leq$3.0 | $0.717^{+0.031}_{-0.031}$ | $0.14^{+0.31}_{-0.31}$ | $0.713^{+0.023}_{-0.023}$ | $0.22^{+0.27}_{-0.27}$ | 0.712
0.0$<\\!\mathcal{S}/\mathcal{N}\\!\leq$4.0 | $0.739^{+0.031}_{-0.026}$ | $0.77^{+0.47}_{-0.47}$ | $0.744^{+0.023}_{-0.023}$ | $0.80^{+0.38}_{-0.38}$ | 0.734
No IA | $0.726^{+0.024}_{-0.042}$ | $-$ | ${0.720^{+0.021}_{-0.031}}$ | $-$ | 0.725
No baryons | $0.732^{+0.032}_{-0.032}$ | $0.71^{+0.49}_{-0.49}$ | ${0.725^{+0.022}_{-0.027}}$ | ${0.70^{+0.43}_{-0.43}}$ | 0.725
No syst | $0.729^{+0.024}_{-0.056}$ | $-$ | ${0.723^{+0.022}_{-0.048}}$ |  $-$ $-$ | 0.709
No GPR error | $-$ | $-$ | $0.732^{+0.019}_{-0.025}$ | $0.71^{+0.40}_{-0.40}$ | 0.708
$N$-body error | $-$ | $-$ | $0.725^{+0.027}_{-0.027}$ | $0.70^{+0.43}_{-0.43}$ | 0.717
No bin1 | $0.734^{+0.043}_{-0.037}$ | $0.10^{+0.74}_{-0.74}$ | ${0.735^{+0.035}_{-0.035}}$ | ${0.10^{+0.70}_{-0.10}}$ | 0.727
No bin2 | $0.740^{+0.042}_{-0.048}$ | $-0.78^{+0.72}_{-0.72}$ | ${0.740^{+0.040}_{-0.040}}$ | ${-0.76^{+0.70}_{-0.70}}$ | 0.738
No bin3 | $0.775^{+0.049}_{-0.055}$ | $0.97^{+0.63}_{-0.63}$ | ${0.777^{+0.048}_{-0.048}}$ | ${0.95^{+0.60}_{-0.60}}$ | 0.836
No bin4 | $0.701^{+0.037}_{-0.037}$ | $0.41^{+0.68}_{-0.68}$ | ${0.702^{+0.029}_{-0.034}}$ | ${0.45^{+0.59}_{-0.59}}$ | 0.659
No bin5 | $0.720^{+0.036}_{-0.028}$ | $0.53^{+0.61}_{-0.61}$ | ${0.723^{+0.030}_{-0.025}}$ | ${0.53^{+0.57}_{-0.57}}$ | 0.716
Peak count analysis of DES-Y1 |
DH21 | – | – | $0.737^{+0.027}_{-0.031}$ | $-$ | $-$
This work | $0.743^{+0.036}_{-0.036}$ | $-$ | $0.742^{+0.030}_{-0.034}$ | $-$ | 0.712
Joint peak count analysis |
| Nautilus | Multinest |
Fiducial | ${0.732^{+0.020}_{-0.020}}$ | ${0.82^{+0.47}_{-0.47}}$ | ${0.732^{+0.012}_{-0.010}}$ | ${0.82^{+0.33}_{-0.33}}$ | 0.745
$\Lambda$CDM | $0.736^{+0.016}_{-0.018}$ | $0.81^{+0.46}_{-0.46}$ | ${0.736^{+0.012}_{-0.015}}$ | ${0.79^{+0.40}_{-0.40}}$ | 0.732
No baryons | $0.728^{+0.020}_{-0.016}$ | $0.82^{+0.46}_{-0.46}$ | ${0.725^{+0.018}_{-0.014}}$ | ${0.83^{+0.39}_{-0.39}}$ | 0.726
No IA | $0.726^{+0.020}_{-0.016}$ | – | ${0.729^{+0.015}_{-0.015}}$ | – | 0.721

###  5.1 Internal consistency

It has been found in previous cosmic shear analyses (e.g. A21, Hamana et al.,
2020; Amon et al., 2022) that internal consistency tests can help
differentiate residual systematics from statistical fluctuations. We therefore
stress-test our results by removing data from tomographic bins one at a time
before proceeding to the inference. For example, we consider results obtained
from an analysis where exactly no data from bin1 (i.e. 1, 1$\cup$2…1$\cup$2,
1$\cup$2$\cup$3 … 1$\cup$2$\cup$3$\cup$4$\cup$5) is used, then no data from
bin2, and so on. The results are shown in Fig. 9, where we observe that all
cases are self-consistent, in agreement with the full selection. Note that the
$S_{8}$ shifts per-bin are not expected to match exactly those measured with
other lensing probes due to different responses of the cosmic shear estimators
to noise in the data. For example, A21 found that removing the fifth
tomographic bin maximally degrades the precision on $S_{8}$, confirming the
large amount of information on this parameter carried by high redshift bins in
shear two-point functions. In contrast, we find here that removing the third
redshift bin has the worst impact on the precision. The third bin has the
greatest number density of galaxies, hence better captures the information in
peak statistics, whose mean value is affected by the noise level. The
constraints on $A_{\rm IA}$ fluctuate about the fiducial results by less than
$2\sigma$, while those on $S_{8}$ agree within $1\sigma$, as expected.

Figure 9: Internal consistency: effect of removing tomographic data from the
KiDS-1000 analysis. Figure 10: Summary of $S_{8}$ constraints from this work,
from recent cosmic shear data analyses and from Planck. This figure shows the
projected 1$\sigma$ errors.

###  5.2 Impact of systematics

We present in this section additional variations with respect to the fiducial
analysis, designed to better understand our results and assess their
robustness to residual systematics. We first investigate the impact of IA on
the uncertainty by fixing $A_{\rm IA}$ to 0.72, the best-fit value in the
fiducial analysis. Doing so, the error bars on $S_{8}$ shrink by less than
10%, while the mean value is not affected, by construction. Setting instead
the IA parameter to 0.0, we can estimate the bias on the inferred cosmology if
IA are completely neglected. We measure in this case $S_{8}^{\rm no-
IA}=0.725^{+0.024}_{-0.042}$, a 0.22$\sigma$ shift from the fiducial results.
Intrinsic alignments are therefore a modest part of the error budget,
suggesting that peak count analyses where IA are not modelled or held fixed
(e.g. M18, HD21, Marques et al., 2023) likely yield both biased low and
slightly optimistic constraints for $S_{8}$.

We next carry out a similar study this time removing the modelling of baryons,
fixing the associated nuisance parameter to $b_{\rm bary}=0.0$. As reported in
Table 4, the measurements are mostly unchanged. As shown in M21, any non-zero
residual feedback tends to lower the number of high $\mathcal{S}/\mathcal{N}$
peaks in all tomographic bins, which, when confronted to fixed data, must be
compensated with an increased value of inferred $S_{8}$. Therefore, removing
the baryon modelling goes the other way and reduces the inferred $S_{8}$. This
is not clearly seen with the Nautilus chains, but the Multinest runs shows
this shift with $0.2\sigma$ significance.

Then, removing modelling of all systematics (photo-$z$, shape calibration, IA
and baryons) results in $S_{8}$ values half way between the no-baryon and no-
IA cases, but the error bars are the larger. This suggests that
marginalisation over these systematics helps in finding the true maximal
likelihood. This is supported by the fact that the error on $S_{8}$ becomes
smaller than the fiducial case if $A_{\rm IA}$ and $b_{\rm bary}$ are fixed to
their best-fit value (0.72 and 0.5, respectively) instead of zero, leading to
$S_{8}^{\rm syst-fixed}=0.728\pm 0.030$, a slight reduction in error compared
to the fiducial error.

The contribution to the error budget coming from the GPR interpolation
uncertainty (Eq. 7) can be estimated from an MCMC run where the covariance
matrix excludes this term, and we observe that the error on $S_{8}$ is reduced
by just under 10%. Similarly, adding an error on small scale non-linear
physics estimated from the scatter between the cosmo-SLICS, Magneticum dark
matter-only and the T17 simulations (see Sec. 4.6), can degrade the error on
$S_{8}$ by 12%. This is an upper limit on the degradation, given that the
Cosmology Training Set has better resolution than these, hence the real
uncertainty is certainly smaller. We do not include this latter error in the
fiducial analysis here, because it is not accurately estimated, and instead
report an upper bound on the effect.

Finally, we compared our fiducial Nautilus results with those from the
Multinest nested sampler and recover negligible biases in the inferred
parameters, but with smaller error bars ($S_{8}=0.733\pm 0.032$ versus
${0.733^{+0.021}_{-0.027}}$ for Multinest). This is consistent with previous
findings (Lemos et al., 2022) and justifies our choice of Nautilus as our main
sampler. We nevertheless report results from both samplers to ease comparison
with previous results.

###  5.3 Comparison with previous KiDS-1000 results

The $S_{8}$ measurement presented here is not the first carried out from
KiDS-1000. Previous analyses include the measurements of A21 and van den Busch
et al. (2022), the latter of which used an upgraded photometric calibration
compared to the former, followed by that of Li et al. (2023a) based on
upgraded shear measurements. Loureiro et al. (2022) carried out a
pseudo-$C_{\ell}$ analysis, Fluri et al. (2022) used instead a convolutional
neural network, while Longley et al. (2023) re-analysed the data within the
LSST-DESC pipeline. We report these results as the purple symbols in Fig. 10,
where we see that they all seem to prefer slightly higher values of $S_{8}$
compared to our own measurements, albeit not by a significant amount. Given
the important differences in the analysis pipelines between these efforts, it
is reassuring to recover $<1\sigma$ agreements. The constraints from the
$w$CDM band power analysis from Tröster et al. (2021) are reported in Fig. 11
and are broadly consistent with our peak statistics constraints, even though
peaks are clearly more constraining on $S_{8}$ ($0.732\pm 0.032$ for peaks vs
$0.742\pm 0.047$ for band power), due to the reduced degeneracy in the
[$S_{8}-w_{0}$] plane. It is worth noting that both statistics provide similar
constraints on the $A_{\rm IA}$ parameter ($\sigma_{A_{\rm IA}}=0.42$ for
peaks, compared to $\sigma_{A_{\rm IA}}=0.36$ for band power), which is
reassuring given that both use the same NLA approach. This error is
significantly reduced ($\sigma_{A_{\rm IA}}=0.30$) when considering the more
aggressive data selection ( $-2.5\leq\mathcal{S}/\mathcal{N}\leq 4.0$), but
since the associated goodness-of-fit is poor, the results are not straight-
forward to interpret. We nevertheless expect tighter constraints on $A_{\rm
IA}$ to be achievable coming from non-Gaussian probes.

We finally remark that our constraints on $\Sigma_{8}^{\alpha}$ aligns
remarkably well with the Band Power measurements presented in A21 (they found
$\Sigma_{8}^{\alpha}=0.765^{+0.018}_{-0.024}$ with $\alpha=0.58$, compared to
our measurement of $\Sigma_{8}^{\alpha}=0.765^{+0.030}_{-0.030}$ with
$\alpha=0.60$).

##  6 Joint analysis with DES-Y1

Figure 11: Comparison with the previous cosmic shear results. The top part
shows a comparison with the KiDS-1000 band power analysis from Tröster et al.
(2021, based on Multinest), while the bottom part compares the KiDS results
with the DES-Y1 peak count analysis from HD21 (black) along with current re-
analysis (green) and detailed in Sec. 6. Note that the posteriors obtained
from the Nautilus sampler are typically wider and more accurate than those
from Multinest used in HD21. The excellent agreement seen in this figure
warrants the joint survey analysis.

The posterior obtained from the KiDS-1000 peak count analysis is fully
consistent with that from the peak count analysis of the Dark Energy Year 1
(DES-Y1) presented in HD21. In particular, the latter finds $S_{8}^{\rm
HD21}=0.737^{+0.027}_{-0.031}$, which significantly overlaps with our
$S_{8}^{\rm KiDS}$ $1\sigma$ results. Other parameters less well measured such
as $\Omega_{\rm m}$ and $w_{0}$ are also largely overlapping at the $1\sigma$
level (see the lower part of Fig. 11), which means the intersection between
the two likelihood hyper-volumes must be large enough to safely combine the
two data sets. Furthermore, both measurements are based on the similar
analysis pipeline and, in particular, exploit the same simulations to model
the cosmology dependence, thereby suppressing the risk of mis-interpreting the
joint data due to non-uniform modelling of the signal.

###  6.1 Results: DES-Y1 re-analysis

As detailed in Sec. 3.3, there are differences between our DES-Y1 pipeline and
that presented in HD21, including the $S_{8}$ sampling, the treatment of
baryons, the inclusion of the emulator uncertainty on the covariance and the
choice of sampler. The results from these re-analyses are presented in the
lower panel of Fig. 11 (in grey and blue). The difference induced on these
contours are small, but the goodness-of-fit improvement is important, with a
PTE of 0.53 (using the same PTE estimator as HD21, we obtain 0.25, which is
still a massive improvement compared to their PTE=0.005.

We remark that our joint pipeline contains a slight inconsistency: we include
IA with the NLA model in the KiDS-1000 data (with marginalisation over $A_{\rm
IA}$) and with the non-linear halo-based IA model for the DES-Y1 data (without
marginalisation, but with an on/off switch instead). We verify the impact of
this feature by analysing the likelihood with the DES IA model turned on and
off and report on the difference, which is sub-dominant ($\Delta
S_{8}=0.002$). We also compare the results from turning off the modelling of
baryons, and from replacing the $w$CDM by a $\Lambda$CDM analysis, finding in
all cases results consistent with the fiducial analysis. The re-analysis
presented in this work has slightly larger error bars compared to that of
HD21, due to the marginalisation over baryons, and to the fact that Nautilus
yields constraints slightly larger compared to Multinest, as summarised in
Table 4. Notably, we infer:

$\displaystyle S_{8}^{\rm DES}=0.743^{+0.036}_{-0.036}\,,\,\,\,\Sigma_{8}^{\rm
DES}=0.762^{+0.036}_{-0.036}\,,\,\,\,{\rm with}\,\alpha=0.559\,.$ (11)

Figure 12: Joint peak count analysis of the KiDS-1000 and DES-Y1 data. In the
upper panel, the green bands indicate the 1$\sigma$ and 2$\sigma$ confidence
intervals from the DES-Y1 re-analysis presented in this paper. The $A_{\rm
IA}$ parameter shown here describes only the IA signal in the KiDS likelihood,
since the IA is fixed to a non-linear model in the DES likelihood (see main
text for details).

###  6.2 Results: joint KiDS + DES

We present in this section the results from our joint KiDS-1000 + DES-Y1 peak
count analysis. Sampling the joint likelihood with our fiducial setup, we
achieve improved constraints on $S_{8}$ and $\Sigma_{8}^{\alpha}$ with:

$\displaystyle S_{8}^{\rm joint,{\it
w}CDM}=0.732^{+0.020}_{-0.020}\,,\,\,\,\Sigma_{8}^{\rm joint,{\it
w}CDM}=0.759^{+0.020}_{-0.017}$ (12)

and

$\displaystyle{S_{8}^{\rm joint,\Lambda
CDM}=0.735^{+0.016}_{-0.018}}\,,\,\,\,\Sigma_{8}^{\rm joint,\Lambda
CDM}=0.762^{+0.017}_{-0.017}$ (13)

computed with $\alpha=0.572$ in both cases. These are the tightest results
obtained from non-Gaussian cosmic shear statistics to date, comparable to the
recent joint $\Lambda$CDM analysis of the KiDS-1000 and DES-Y3 data (Dark
Energy Survey and Kilo-Degree Survey Collaboration, 2023), which measured
$S_{8}^{\rm NLA}=0.792^{+0.016}_{-0.013}$. The two-dimensional posterior is
shown in Fig. 12 (in blue) and compared to the fiducial KiDS-1000 (red) and
DES-Y1 (green) peak statistics constraints. Recall that the $A_{\rm IA}$
parameter affects only the KiDS likelihood since, as explained in the previous
section, the DES likelihood assumes instead a fixed halo-based IA model with
no free parameter. We should therefore use caution when interpreting this
parameter. The reported value is close to the point of maximum likelihood
(${S_{8}^{\rm ML}=0.728}$), and the size of the error bars on $S_{8}$ is
consistent with our expectation: for example, we read from Table 4 that the
$\Lambda$CDM KiDS-1000 analysis has a mean error of $\sigma_{S_{8}}=0.029$.
Scaling this precision by the square root of the area, we naively predict a
joint survey error of around 0.018, and obtain 0.017. The error would be
slightly larger had we included as well a marginalisation of the IA in the
DES-Y1 part of this analysis, possibly explaining this slight difference. At
the joint best-fit cosmology, the PTE values for the KiDS and DES pipelines
are basically unchanged, while the joint analysis has a $\chi^{2}_{\rm
red}=1.15$ and a PTE of 0.96, all satisfying our goodness-of-fit criteria.

If we restrict the joint analysis to $w_{0}=-1.0$, the $S_{8}$ values are
minimally affected while the uncertainty is reduced, as expected from lowering
the dimensionality of the likelihood. Alternatively, turning on the IA
modelling in the DES likelihood only yields a 0.2$\sigma$ downward shift, also
expected whenever IA modelling is added. The smallness of this shift is once
again showing that the intrinsic alignment do not significantly impact the
peak count statistics as measured in the DES-Y1 data. In comparison, setting
to zero the IA model in both KiDS and DES results in ${S_{8}^{\rm joint,no-
IA}=0.725^{+0.020}_{-0.016}}$. Holding fixed the baryonic feedback parameter
to $b_{\rm bary}=0.0$ has similar consequences on this joint analysis,
shifting the best fit value to $S_{8}^{\rm joint,no-
bary}=0.727^{+0.020}_{-0.016}$, a 0.25$\sigma$ shift compared to the fiducial
case. All these values are summarised in Table 4 and in Fig. 10 (with the
brown symbols).

The dark energy equation-of-state is constrained from this joint analysis,
with

$\displaystyle w_{0}^{\rm joint}=-1.12^{+0.42}_{-0.31}\,,$ (14)

which is the first measurement of this quantity from peak statistics, and
arguably one of the best from cosmic shear-only data analyses. The upper limit
is close to the prior edge on $w_{0}$, which might lead to a slight under-
estimation of the error on this side. However, this measurement is robust
against the choice of sampler ($w_{0}=-1.09^{+0.29}_{-0.29}$ for Multinest),
against baryon modelling ($w_{0}=-1.05^{+0.51}_{-0.22}$ setting $b_{\rm
bary}=0.0$), IA ($w_{0}=-1.13^{+0.44}_{-0.33}$ setting $A_{\rm IA}=0.0$) and
scale cuts ($w_{0}=-0.958^{+0.45}_{-0.093}$ when including the aggressive
$\mathcal{S}/\mathcal{N}$ cut in the KiDS-1000 data vector). As shown in
Martinet et al. (2020), aperture-mass maps statistics are highly sensitive to
dark energy and these results seem to be showing exactly that. Previously, the
shear two-point function measurement from Troxel et al. (2018) on DES-Y1
achieved $w_{0}=-0.77^{+0.30}_{-0.37}$ when varying the baryonic feedback
model, using the Multinest sampler. The GCNN analysis of Fluri et al. (2022)
was also able to set constraints on dark energy, with
$w_{0}=-0.93^{+0.32}_{-0.29}$, although they recognise that their results are
affected by the prior boundary on the low side, just like ours is on the high
side202020We review the A21 definition that constraints are uninformed by the
prior when the posterior drops below 0.135 of its maximum at the edges of a
uniform prior volume. In the case of $w_{0}$ we find that the posterior is
slightly above this threshold (0.156) at the upper edge. As this is at the
borderline of the A21 criteria, we therefore caution that the error on this
side might be slightly under-estimated.. Similarly, HD21 found $w_{0}>-1.5$,
also prior-dominated on one side. Other cosmic shear measurements of $w_{0}$
involve additional data (Abbott et al., 2023; Tröster et al., 2021), making
this an unfair comparison.

It is worth mentioning that all peak count analyses based on the cosmo-SLICS
yield $S_{8}$ constraints that are lower than the 2PCFs fiducial analyses.
This could be pointing to limitations in the training set, but is quite
speculative at this stage given that the $\Sigma_{8}$ values align well.
Further investigations and novel simulation suites would be required to
ascertain this, which we post-pone for future work.

###  6.3 Tension with Planck

Figure 13:  Joint constraint on $S_{8}$, $w_{0}$ and $\Sigma_{8}^{\alpha}$,
where $\alpha=0.572$, comparing here the combined cosmic shear surveys with
the CMB results. The tension is lower on $\Sigma_{8}$ than on $S_{8}$.

The $S_{8}$ tension between recent CMB anisotropy and weak lensing data
analyses is drawing a lot of attention, as it could point towards new physics
or hidden systematics (see, e.g. Abdalla et al., 2022a, for a review). The
Planck mission reports $S_{8}^{\it Planck}=0.830\pm 0.013$ (Planck
Collaboration et al., 2020), which is higher than many lensing results (see
Amon et al., 2023, and references therein). The tension $\tau$ can be
evaluated with a number of metrics, and we use here a relatively simple one
used in A21, which compares the difference in the mean with the combined
variances, ${\rm var}[S_{8}]^{\it Planck},{\rm var}[S_{8}]^{\rm peaks}$ :

$\displaystyle\tau=\frac{S_{8}^{\it Planck}-S_{8}^{\rm peaks}}{\sqrt{{\rm
var}[S_{8}]^{\it Planck}+{\rm var}[S_{8}]^{\rm peaks}}}$ (15)

With the fiducial setup shown in Table 4, and using the above definition, our
results from the KiDS-1000 peak count analysis are in $\tau=2.0\sigma$ tension
with the Planck nominal constraints on $S_{8}$ from their $w$CDM analysis.
Tröster et al. (2021) finds a similarly low tension with the same KiDS-1000
data in a $w$CDM analysis, either using this simple tension metric or a more
sophisticated method based on the full shape of the posteriors. Similarly, we
evaluate our joint KiDS-DES analysis to be in $\tau=2.7\sigma$ tension with
Planck, an increase that is driven by the decrease in error bars. Using
instead the MAP values (listed in Table 4) lowers the KiDS tension to
1.7$\sigma$, and the joint-survey to 1.65.

The tension reaches $\tau=4.1\sigma$ in our joint-survey $\Lambda$CDM
analysis, which could be pointing to a resolution of the $S_{8}$ tension that
includes modification to the dark energy equation of state, as suggested by
Tröster et al. (2021) and by the recent results from the Dark Energy
Spectroscopic Instrument (DESI Collaboration et al., 2024), however such
statement cannot be definitive until better $w$ measurements are obtained from
cosmic shear data.

Note that this tension is not only seen in weak lensing, but also in other
late-time probes, e.g. data involving galaxy clustering, as recently
summarised in Alonso et al. (2023, see their figure 7) and reviewed with
greater details in Abdalla et al. (2022b). The current work aligns with the
existing trend, without providing an obvious solution. Again, large
unaccounted contributions from IA and baryons could push the inferred $S_{8}$
value towards Planck, but our analysis prefers lower values: in particular, we
measure $b_{\rm bary}<1.05$ at 95% CL in the KiDS-1000 analysis, and $<0.82$
in the joint analysis, excluding baryonic feedback models that are stronger
than the Magneticum. Also, the redshift estimation methods used to analyse the
DES-Y1 data are suboptimal compared to recent developments (see e.g.
Hildebrandt et al., 2021), potentially causing biases of up to 0.03 in the
inferred $S_{8}$ value. The observed tension with Planck would be different if
that bias was real.

Interestingly, the tension on $\Sigma_{8}$ is reduced to $\tau=0.72\sigma$
with the KiDS-1000 data, and to $\tau=1.33\sigma$ with the joint data, in both
cases evaluated at the $\alpha$ value preferred by the cosmic shear
measurements212121The Planck $w$CDM measurements of $\Sigma_{8}$ are of
$0.793^{+0.019}_{-0.031}$ and $0.797^{+0.018}_{-0.028}$ for the KiDS-1000 and
joint-survey $\alpha$ parameters, respectively.. This is in part caused by a
degradation of the CMB constraints along this quantity, which completely
relaxes the tension, and in part by a different projection angle of the high-
dimensional posterior, as seen in Fig. 13. At the same time, when holding $w$
fixed, the tension on $\Sigma_{8}$ is again increased reaching $2.3\sigma$ for
KiDS, slightly lower than the $3\sigma$ reported in A21, and $3.1\sigma$ for
our joint analysis. To summarise, the tension with Planck on $S_{8}$ is
lowered in $w$CDM compared to $\Lambda$CDM, and is further lowered when
considering the more robust $\Sigma_{8}$ parameter instead of $S_{8}$. This is
discussed in more details in Appendix B.

A more in-depth analysis of this tension requires a robust determination of
the inferred $\Omega_{\rm m}$ parameter. We defer this improvement to future
work, which will benefit from a denser training set, yielding a more accurate
and robust emulator.

##  7 Conclusions

We report in this paper a 4.4% measurement of $S_{8}$ from the tomographic
peak count statistics measured from the KiDS-1000 data. Our simulation-based
inference method exploits the non-linear features extracted from aperture mass
statistics, sensitive to scales as small as 2.0 arcmin. We model the
cosmological dependence with simulated $w$CDM weak lensing light-cones, we
estimate the covariance matrix numerically, and forward model the effect of
intrinsic alignments, baryonic feedback, photometric redshift error and galaxy
shape mis-calibrations. We find a value of $S_{8}^{\rm
KiDS}=0.733^{+0.032}_{-0.032}$, which aligns well with previous KiDS-1000
measurements. We show that our results are robust to residual systematics and
that, of these, intrinsic alignment of galaxies plays the most important role,
shifting the best-fit $S_{8}$ value by 0.22$\sigma$ if left unmodelled. We
also show that the most robustly measured parameter in our analysis is
$\Sigma_{8}^{\rm KiDS}\equiv\sigma_{8}\left(\Omega_{\rm
m}/0.3\right)^{\alpha}=0.765^{+0.030}_{-0.030}$, with $\alpha=0.60$, in
excellent agreement with previous KiDS-1000 analysis.

The inferred posterior distribution is consistent with the peak count
measurement carried out on the DES-Y1 data using a similar analysis pipeline
(HD21), allowing us to jointly analyse the two data sets, which yields
$S_{8}^{\rm joint}=0.732^{+0.020}_{-0.020}$, one of the tightest constraints
on this parameter from lensing data alone. The combined data sets have enough
statistical precision to allow the first measurement of the dark energy
equation-of-state parameter from non-Gaussian statistics: $w_{0}^{\rm
joint}=-1.12^{+0.42}_{-0.31}$, in agreement with the $\Lambda$CDM scenario,
and robust to variations in the analysis choices.

Our best-fit $S_{8}^{\rm joint}$ is also in statistical agreement with all
previous KiDS-1000 analyses and with the HSC-Y3 and DES-Y3 $\gamma$-2PCF
results, but lower than the DES-Y3 measurements from peaks and moments, and in
$2.7\sigma$ tension with Planck. This joint-survey tension increases to
4$\sigma$ in our $\Lambda$CDM analysis, but lowers to 2.3 and 3.1$\sigma$ when
considering instead the $\Sigma_{8}$ parameter, for the KiDS-only and joint
survey analyses, respectively.

Our pipeline has been thoroughly tested, however we recognise it is
incomplete. As detailed in Sec. 4, we hold fixed a number of cosmological
parameters, which likely affect our results, including $\Omega_{\rm b},n_{\rm
s}$ and $m_{\nu}$. We also consider a single baryonic feedback model (although
we allow its amplitude to vary), knowing that other hydrodynamical simulations
would provide slightly different responses. Furthermore, we model IA with the
redshift-independent NLA model, which we know is an incomplete effective
model, and we neglect source clustering, as it was shown to be completely
subdominant. Finally, we have identified a systematic effect in our Gaussian
process emulator that is caused by the relatively small number of training
nodes, preventing us from extracting meaningful information about $\Omega_{\rm
m}$, however all other parameters are unaffected by this. This limits our
ability to further study in higher-dimensional space the potential $S_{8}$
tension with the CMB. Addressing these will therefore be the object of future
work. It will also be informative to compare our results to other non-Gaussian
probes of cosmic shear, and possibly combine the methods to further reduce the
uncertainty on $S_{8}$ and $w_{0}$.

## Acknowledgements

We would like to thank Dr. Eric Tittley from the Royal Observatory of
Edinburgh for his continuous help and support with the Cuillin cluster on
which the bulk of the computations were carried on, as well as Elena Sellentin
for discussions on the use of Hotelling’s $T^{2}$ statistics.

JHD acknowledges support from an STFC Ernest Rutherford Fellowship (project
reference ST/S004858/1). SH is supported by the U.D Department of Energy,
Office of Science, Office of High Energy Physics under Award Number DE-
SC0019301. BG acknowledges support from a UK STFC Consolidated Grant and from
the Royal Society of Edinburgh through the Saltire Early Career Fellowship
(ref. number 1914). NM acknowledges support from the French Agence Nationale
de la Recherche (PISCO project / ANR-22-CE31-0004-01). TT acknowledges funding
from the Swiss National Science Foundation under the Ambizione project
PZ00P2_193352. KD is supported by the grant agreements ANR-21-CE31-0019 /
490702358 from the French Agence Nationale de la Recherche / Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) for the LOCALIZATION
project. TC is supported by the INFN INDARK PD51 grant, by the Fondazione ICSC
National Recovery and Resilience Plan (PNRR) Project ID CN-00000013 “Italian
Research Center on High-Performance Computing, Big Data and Quantum Computing”
funded by MUR Missione 4 Componente 2 Investimento 1.4: “Potenziamento
strutture di ricerca e creazione di “campioni nazionali” di R$\&$S (M4C2-19)”
– Next Generation EU (NGEU), and by the FARE MIUR grant ‘ClustersXEuclid’. HH
is supported by a DFG Heisenberg grant (Hi 1495/5-1), the DFG Collaborative
Research Center SFB1491, as well as an ERC Consolidator Grant (No. 770935). BJ
acknowledges support by the ERC-selected UKRI Frontier Research Grant
EP/Y03015X/1 and by STFC Consolidated Grant ST/V000780/1. AHW is supported by
the Deutsches Zentrum für Luft- und Raumfahrt (DLR), made possible by the
Bundesministerium für Wirtschaft und Klimaschutz, and acknowledges funding
from the German Science Foundation DFG, via the Collaborative Research Center
SFB1491 ”Cosmic Interacting Matters - From Source to Signal”.

The KiDS-1000 results in this paper are based on data products from
observations made with ESO Telescopes at the La Silla Paranal Observatory
under programme IDs 177.A-3016, 177.A-3017 and 177.A-3018, and on data
products produced by Target/OmegaCEN, INAF-OACN, INAF-OAPD and the KiDS
production team, on behalf of the KiDS consortium. This project also used
public archival data from the Dark Energy Survey (DES). Funding for the DES
Projects has been provided by the U.S. Department of Energy, the U.S. National
Science Foundation, the Ministry of Science and Education of Spain, the
Science and Technology FacilitiesCouncil of the United Kingdom, the Higher
Education Funding Council for England, the National Center for Supercomputing
Applications at the University of Illinois at Urbana-Champaign, the Kavli
Institute of Cosmological Physics at the University of Chicago, the Center for
Cosmology and Astro-Particle Physics at the Ohio State University, the
Mitchell Institute for Fundamental Physics and Astronomy at Texas A&M
University, Financiadora de Estudos e Projetos, Fundação Carlos Chagas Filho
de Amparo à Pesquisa do Estado do Rio de Janeiro, Conselho Nacional de
Desenvolvimento Científico e Tecnológico and the Ministério da Ciência,
Tecnologia e Inovação, the Deutsche Forschungsgemeinschaft, and the
Collaborating Institutions in the Dark Energy Survey. The Collaborating
Institutions are Argonne National Laboratory, the University of California at
Santa Cruz, the University of Cambridge, Centro de Investigaciones
Energéticas, Medioambientales y Tecnológicas-Madrid, the University of
Chicago, University College London, the DES-Brazil Consortium, the University
of Edinburgh, the Eidgenössische Technische Hochschule (ETH) Zürich, Fermi
National Accelerator Laboratory, the University of Illinois at Urbana-
Champaign, the Institut de Ciències de l’Espai (IEEC/CSIC), the Institut de
Física d’Altes Energies, Lawrence Berkeley National Laboratory, the Ludwig-
Maximilians Universität München and the associated Excellence Cluster
Universe, the University of Michigan, the National Optical Astronomy
Observatory, the University of Nottingham, The Ohio State University, the
OzDES Membership Consortium, the University of Pennsylvania, the University of
Portsmouth, SLAC National Accelerator Laboratory, Stanford University, the
University of Sussex, and Texas A&M University. Based in part on observations
at Cerro Tololo Inter-American Observatory, National Optical Astronomy
Observatory, which is operated by the Association of Universities for Research
in Astronomy (AURA) under a cooperative agreement with the National Science
Foundation.

The chains presented in this paper are plotted with GetDist (Lewis, 2019),
available at https://getdist.readthedocs.io.

All authors contributed to the development and writing of this paper. The
authorship list is given in three groups: the lead authors (JHD, SH), followed
by two alphabetical groups. The first alphabetical group includes those who
are key contributors to both the scientific analysis and the data products
(BG, NM, TT). The second group covers those who have either made a significant
contribution to the data products, or to the scientific analysis.

## Data Availability

The SLICS numerical simulations can be found at http://slics.roe.ac.uk/, while
the SLICS-HR, the cosmo-SLICS and the Magneticum can be made available upon
request. This work also uses public KiDS-1000 and DES-Y1 data, which can be
found at https://kids.strw.leidenuniv.nl/DR4/index.php and
https://des.ncsa.illinois.edu/releases/y1a1, respectively.

## References

  * Abbott et al. (2018) Abbott, T. M. C. et al. 2018, ApJS, 239, 18, 1801.03181
  * Abbott et al. (2022) ——. 2022, Phys. Rev. D, 105, 023520, 2105.13549
  * Abbott et al. (2023) ——. 2023, Phys. Rev. D, 107, 083504, 2207.05766
  * Abdalla et al. (2022a) Abdalla, E. et al. 2022a, Journal of High Energy Astrophysics, 34, 49, 2203.06142
  * Abdalla et al. (2022b) ——. 2022b, Journal of High Energy Astrophysics, 34, 49, 2203.06142
  * Ajani et al. (2020) Ajani, V., Peel, A., Pettorino, V., Starck, J.-L., Li, Z., & Liu, J. 2020, Phys. Rev. D, 102, 103531
  * Akeson et al. (2019) Akeson, R. et al. 2019, arXiv e-prints, arXiv:1902.05569, 1902.05569
  * Alarcon et al. (2021) Alarcon, A. et al. 2021, MNRAS, 501, 6103, 2007.11132
  * Alonso et al. (2023) Alonso, D., Fabbian, G., Storey-Fisher, K., Eilers, A.-C., García-García, C., Hogg, D. W., & Rix, H.-W. 2023, arXiv e-prints, arXiv:2306.17748, 2306.17748
  * Alonso et al. (2019) Alonso, D., Sanchez, J., Slosar, A., & LSST Dark Energy Science Collaboration. 2019, MNRAS, 484, 4127
  * Amon et al. (2022) Amon, A. et al. 2022, Phys. Rev. D, 105, 023514
  * Amon et al. (2023) ——. 2023, MNRAS, 518, 477