# On constraining Cosmology and the Halo Mass Function with Weak Gravitational
Lensing.

Shiming Gu (顾时铭)1 Marc-Antoine Dor,1,2 Ludovic van Waerbeke,1 Marika
Asgari,3,4,5 Alexander Mead,6,8 Tilman Tröster,5,7 and Ziang Yan (颜子昂)8
1Department of Physics and Astronomy, University of British Columbia,
Vancouver, BC, V6T 1Z1, Canada
2ENSTA Paris, 828, Boulevard des Maréchaux, 91762 Palaiseau Cedex, France
3E.A Milne Centre for Astrophysics, University of Hull, Cottingham Road,
Kingston upon Hull, HU6 7RX, UK
4Centre of Excellence for Data Science, AI, and Modelling (DAIM), University
of Hull, Cottingham Road, Kinston upon Hull, HU6 7RX
5Institute for Astronomy, University of Edinburgh, Royal Observatory,
Blackford Hill, Edinburgh, EH9 3HJ, UK
6Department of Computer Science, University of British Columbia, Vancouver, BC
V6T 1Z1, Canada
7Eidgenoessische Technische Hochschule Zürich, Rämistrasse 101, 8092 Zurich,
Switzerland
8Ruhr University Bochum, Faculty of Physics and Astronomy, Astronomical
Institute (AIRUB), German Centre for Cosmological Lensing, 44780 Bochum,
Germany
E-mail<EMAIL_ADDRESS>(SG)

(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

The discrepancy between the weak lensing (WL) and the Planck measurements of
${S_{8}}$ has been a subject of several studies. Assuming that residual
systematics are not the cause, these studies tend to show that a strong
suppression of the amplitude of the mass power spectrum $P(k)$ in the late
universe at high $k$ could resolve it. The WL signal at the small scale is
sensitive to various effects not related to lensing, such as baryonic effects
and intrinsic alignment. These effects are still poorly understood, therefore,
the accuracy of $P(k)$ depends on the modelling precision of these effects. A
common approach for calculating $P(k)$ relies on a halo model. Amongst the
various components necessary for the construction of $P(k)$ in the halo model
framework, the halo mass function (HMF) is an important one. Traditionally,
the HMF has been assumed to follow a fixed model, motivated by dark matter-
only numerical simulations. Recent literature shows that baryonic physics,
amongst several other factors, could affect the HMF. In this study, we
investigate the impact of allowing the HMF to vary. This provides a way of
testing the validity of the halo model-HMF calibration using data. In the
context of the aforementioned ${S_{8}}$ discrepancy, we find that the Planck
cosmology is not compatible with the vanilla HMF for both the DES-y3 and the
KiDS-1000 data. Moreover, when the cosmology and the HMF parameters are
allowed to vary, the Planck cosmology is no longer in tension. The modified
HMF predicts a matter power spectrum with a $\sim 25\%$ power loss at $k\sim
1~{}{\rm h/Mpc}$, in agreement with the recent studies that try to mitigate
the ${S_{8}}$ tension with modifications in $P(k)$. We show that Stage IV
surveys will be able to measure the HMF parameters with a few percent
accuracy.

###### keywords:

gravitational lensing: weak – cosmology: cosmological parameters – dark matter
– large-scale structure of Universe

††pubyear: 2023††pagerange: On constraining Cosmology and the Halo Mass
Function with Weak Gravitational Lensing.–A.2

##  1 Introduction

The presence of dark matter in our Universe has been firmly established, with
baryonic matter accounting for approximately 20% of the total mass, while dark
matter (DM) comprises around 80% (e.g. Planck Collaboration et al., 2020;
Aiola et al., 2020; eBOSS Collaboration et al., 2020). Weak gravitational
lensing is a well-known technique employed to investigate the mass
distribution by analyzing the distorted shapes of background galaxies (Gunn,
1967; Blandford et al., 1991). Over the past decade, several cosmic shear
surveys have been conducted, including the Canada-France-Hawaii-Telescope
Lensing Survey (CFHTLens) (Heymans et al., 2012), the Kilo Degree Survey
(KiDS) (Kuijken et al., 2015), the Dark Energy Survey (DES) (Drlica-Wagner et
al., 2018), and the Hyper Suprime-Cam (HSC) (Aihara et al., 2018).

The cosmological interpretation of weak lensing measurements necessitates the
computation of the matter power spectrum $P(k)$, which is a function of
redshift and cosmological parameters. However, calculating $P(k)$ accurately
is challenging, particularly at the small scale, where some contributions are
not constrained well enough. These contributions arise from factors such as
intrinsic alignment, baryonic feedback, and neutrino streaming, among others.
Consequently, semi-analytical models of $P(k)$ are calibrated with numerical
simulations (e.g. Peacock & Smith, 2000; Smith et al., 2003; Takahashi et al.,
2012; Bird et al., 2012). The issue is that differences in the calibration
recipes, employed by different models, lead to different $P(k)$, which
consequently impact the cosmological parameter inference at the level of a few
percents. While this problem has been recognized for quite some time, several
modifications to the calibration procedure have been implemented. However, as
future lensing surveys strive for higher accuracy, the calibration of $P(k)$
may once again become the main limiting factor, thereby preventing lensing
surveys from reaching their full potential. The ${S_{8}}$ tension observed
between low and high
redshifts111${S_{8}}={\sigma_{8}}\sqrt{{\Omega_{\mathrm{m}}}/0.3}$ could be an
indication of this issue. If this tension is genuine, its underlying cause
(residual systematics, modelling of $P(k)$, or new physics), remains unknown.

Recently, the joint analysis of the Dark Energy Survey (DES-y3) and the Kilo-
Degree Survey (KiDS-1000) was conducted by DES and KiDS collaboration et al.
(2023). They discovered that by swapping the analysis pipeline of each survey,
the ${S_{8}}$ tension between the surveys increased significantly, hence
exacerbating the tension between KiDS and Planck while decreasing it between
DES and Planck. However, this issue was resolved by adopting a hybrid analysis
pipeline that incorporates elements from both surveys. The solution emerged
from a tuning of the intrinsic alignment and baryonic effects models, tested
on numerical simulations. Consequently, the ${S_{8}}$ tension with the cosmic
microwave background diminished, while KiDS-1000 and DES-y3 exhibited good
agreement. It should be noted that this does not provide a definitive solution
for modelling $P(k)$, but rather emphasizes the challenges associated with
developing a reliable and predictive high-precision model for $P(k)$. In DES
and KiDS collaboration et al. (2023) and other papers aimed at mitigating the
${S_{8}}$ tension, the focus consistently revolves around the modelling of
$P(k)$. For instance, Amon & Efstathiou (2022) proposed a one-parameter
modification of the matter power spectrum on a non-linear scale.222However,
they demonstrated that achieving a reduction in the ${S_{8}}$ tension would
necessitate an unrealistic suppression of power at the non-linear scale.

Not all options impacting the modelling of $P(k)$ have been explored. For
instance, semi-analytical models, such as the halo model based HMcode (Mead et
al., 2015), are widely used to calculate $P(k)$. At the heart of any halo
model, there is a Halo Mass Function (HMF) (e.g., Press & Schechter, 1974;
Sheth & Tormen, 1999) a fundamental building block from which $P(k)$ is
constructed. The parameters describing the HMF are set to fixed values,
usually determined by DM only $N$-body simulations, and other modifications
are supposed to take care of all the complexities of non-linear physics. The
HMF is not an exact description of the haloes distribution in the universe,
but it is a proxy that serves this purpose. As summarized below, there are
good reasons for wanting the parameters of the HMF to be adjustable, which is
the main motivation of our study.

A direct measurement of the HMF would significantly improve the modelling of
$P(k)$. Unfortunately, this is not possible, and, at best, the HMF can be
partially constrained using baryonic probes as tracers of DM haloes, which
involves various baryon-based quantities. Castro et al. (2016) discussed the
constraints of Sheth-Tormen HMF models (Sheth & Tormen, 1999; Despali et al.,
2016) by analyzing observations of galaxy cluster counts, cluster power
spectra, and Type Ia supernova lensing. One possibility is to use the size and
velocity dispersion relations of galaxy clusters and/or groups (e.g. Bahcall &
Cen, 1993; Böhringer et al., 2017), but this method relies on a simulation-
calibrated relationship between mass and X-ray luminosities (Hoekstra et al.,
2011). Previous work using $N$-body simulations has shown that modifications
of the Halo Mass Function (HMF) can arise from distinct initial conditions
and/or dark matter physics. For example, Bagla et al. (2009) demonstrated that
the HMF changes depending on the spectral indices of the primordial power
spectrum, while Costanzi et al. (2013) showed that the HMF’s shape depends on
the neutrino mass and matter density. Barreira et al. (2014) showed that the
spherical collapse is changed under an alternate theory of gravity, thereby
causing the HMF calibration from DM-only N-body simulations to be inaccurate.
Using hydro simulations, Bocquet et al. (2016) showed that baryonic physics
impacts the shape of the HMF, and Baugh et al. (2019) showed that the HMF
varies by $25\%$ between the WMAP7-based and Planck-based mock catalogues due
to diverse gas cooling rates triggered by different $\omega_{\mathrm{b}}$
values (White & Frenk, 1991). The HMF also depends on the specifics of the
halo-finding algorithm (Knebe et al., 2011). In summary, it was shown in
various papers that the HMF is not universal that it depends on underlying
physical processes, cosmology and selection effects (Ondaro-Mallea et al.,
2022; Diemer, 2021; Castro et al., 2021). Therefore, we believe that it is
important to incorporate the HMF parameters as degrees of freedom in the
modelling of $P(k)$ in order to better represent our ignorance and see if the
fixed values of the HMF parameters are recovered by the fitting process. If
they are not, this would indicate a problem with the data and/or the
modelling. Our approach is inspired by the shape measurement process to test
residual systematics, where a shear calibration parameter $m$ was introduced
at the model fitting stage, which should be zero if no residual systematics
are present in the galaxy shapes measurement; all lensing surveys generally
find a small but non-zero value. Similarly, the free HMF parameters are
introduced at the model fitting level in order to verify if the theoretically
inferred values are acceptable.

This paper is structured as follows. In section 2, we provide an introduction
to the halo mass function and explain how it affects weak lensing observables.
In section 1, we discuss the modelling process and the data used in this
study. Parameter interference from the observables is discussed in section
4.1. In section 4, we present the parameter constraint on the HMF and
cosmology based on the new HMF. Finally, we summarize our findings and
conclude in section 5.

##  2 Theory

###  2.1 The Shear Correlation Function

This study is based on the comparison between the measured shear correlation
function and the theory, using the halo model approach and the DES-y3 and KiDS
surveys. We decided to base our analysis on cosmic shear only and not include
galaxy clustering nor galaxy-galaxy lensing because we want an estimator
completely free of any assumption of galaxy bias, assembly bias, or any
consideration regarding on how galaxies populate haloes.

The shear two-point correlation functions is $\xi_{\pm}^{\rm ij}(\theta)$,
where $i$ and $j$ refer to redshift bins. The estimator $\hat{\xi}_{\pm}^{\rm
ij}(\theta)$ is given by:

$\hat{\xi}_{\pm}^{\mathrm{ij}}(\theta)=\frac{\sum_{a,b}w_{\mathrm{a}}w_{\mathrm{b}}\left[\epsilon^{i}_{t}({{\hbox{$\bf\theta$}}}_{a})\epsilon^{j}_{\mathrm{t}}({{\hbox{$\bf\theta$}}}_{\mathrm{b}})\pm\epsilon^{i}_{\times}({{\hbox{$\bf\theta$}}}_{\mathrm{a}})\epsilon^{j}_{\times}({{\hbox{$\bf\theta$}}}_{\mathrm{b}})\right]}{\sum_{a,b}w_{\mathrm{a}}w_{\mathrm{b}}},$
(1)

where ${{\hbox{$\bf\theta$}}}_{\mathrm{a}}$ and
${{\hbox{$\bf\theta$}}}_{\mathrm{b}}$ are galaxies’ position vectors and the
summation is taken over galaxy pairs $(a,b)$. The galaxies have an angular
separation
$|{{\hbox{$\bf\theta$}}}_{\mathrm{a}}-{{\hbox{$\bf\theta$}}}_{\mathrm{b}}|$
within an interval $\Delta\theta$ centered on $\theta$. The galaxy weights are
represented as $w_{\mathrm{a}}$ and $w_{\mathrm{b}}$, while
$\epsilon_{\mathrm{t},\times}^{i,j}$ refers to the tangential and cross
ellipticities of galaxies in redshift bins $i$ and $j$.

The estimator $\hat{\xi}_{\pm}^{\rm ij}(\theta)$ can be expressed in terms of
the convergence power spectrum $C_{\kappa\kappa}^{ij}(\ell)$ (Hamilton, 2000):

$\xi_{\pm}^{ij}(\theta)=\frac{1}{2\pi}\int_{0}^{\infty}\text{d}\ell\ \ell
J_{0,4}(\ell\theta)C_{\kappa\kappa}^{ij}(\ell),$ (2)

Here, $J_{0,4}$ represents the 0- and 4-th orders of the Bessel functions of
the first kind. For scales $\ell>10$, it is possible to approximate the sky as
flat, and the angular cross-spectrum $C_{\kappa\kappa}^{ij}(\ell)$ is related
to the 3D power spectrum $P_{\delta}$ using the Limber approximation (Limber,
1953):

$C_{\kappa\kappa}^{ij}(\ell)=\int_{0}^{\chi_{\rm
H}}\text{d}\chi\frac{q_{i}(\chi)q_{j}(\chi)}{[f_{K}(\chi)]^{2}}P\left(k=\frac{\ell+\dfrac{1}{2}}{f_{K}(\chi)},\chi\right),$
(3)

In the above equation, $\chi$ is the radial comoving distance, $f_{K}(\chi)$
is the angular diameter distance, $\chi_{\rm H}$ is the horizon distance, and
$P(k,\chi)$ is the 3D matter power spectrum. The lensing efficiency function,
for the redshift bin $i$, $q_{i}(\chi)$, is given by:

$q_{\mathrm{i}}(\chi)=\frac{3H_{0}^{2}{\Omega_{\mathrm{m}}}}{2c^{2}}\int_{\chi}^{\chi_{\rm
H}}\frac{\text{d}\chi^{\prime}}{a(\chi)}n_{\mathrm{i}}(\chi^{\prime})\frac{f_{K}(\chi^{\prime}-\chi)f_{K}(\chi)}{f_{K}(\chi^{\prime})},$
(4)

where, $H_{0}$ is the Hubble Constant (Hubble, 1936), ${\Omega_{\mathrm{m}}}$
the matter density fraction of the universe, $c$ the speed of light, $a(\chi)$
the dimensionless scale factor, and $n_{\mathrm{i}}(\chi)$ is the effective
number of galaxies in bin $i$ normalized to $\int_{0}^{\chi_{\rm
H}}\text{d}\chi n(\chi)=1$.

The power spectrum $P(k)$ is the contribution of two parts. The first part is
due to matter clustering within a halo and is called the one-halo term,
$P_{\text{1H}}(k)$. The second part arises from clustering between haloes and
is called the two-halo term, $P_{\text{2H}}(k)$:

$P(k)=P_{\text{1H}}(k)+P_{\text{2H}}(k).$ (5)

In the following two Sections, we introduce the Halo Model in general terms
and then describe the specific HMcode (Mead et al., 2015) implementation that
we use to do the calculations.

###  2.2 The Halo Profile and the Halo Mass Function

A recent review of the halo model can be found in Asgari et al. (2023). The
theoretical ingredients for the halo model are the halo mass density profile
$\rho(r,M)$ and the halo mass function (HMF) $f(\nu)$, where $\nu$ is the peak
threshold of a halo of mass $M$, defined as:

$\nu=\frac{\delta_{c}(z)}{\sigma(M,z)},$ (6)

where $\delta_{c}(z)$ is the linear-theory-based critical threshold for
spherical collapse, and the variance $\sigma^{2}(M,z)$ is computed from the
linear matter power spectrum $\Delta^{2}_{\mathrm{lin}}(k,z)$:

$\sigma^{2}(M,z)=\int_{0}^{\infty}\text{d}k\
k^{-1}\Delta^{2}_{\mathrm{lin}}(k,z)W^{2}(k,r),$ (7)

where $W(k,r)$ is a Fourier transformed spherical top-hat function. The
function $f(\nu)$ can be expressed in terms of the number density $n(M)$ of
haloes of mass $M$:

$f(\nu)=n(M)\frac{\bar{\rho}}{M}\frac{\text{d}\nu}{\text{d}M},$ (8)

where $\bar{\rho}$ is the mean matter density. The shape of $f(\nu)$ is
determined by two parameters, $p$ and $q$, as described in Sheth & Tormen
(1999). The parameter $p$ is associated with the tidal field between
overlapping DM haloes, which determines the abundance of low-mass haloes. The
parameter $q$ is associated with the boundary crossing between two haloes,
which modifies the spherical collapse of the DM halo by correcting the peak
threshold $\nu^{\prime 2}=q\nu^{2}$ (Sheth et al., 2001). The function
$f(\nu)$ and the normalisation constant $A$ are given by:

$f(\nu)=A\sqrt{\frac{2q}{\pi}}\left(1+\left(\frac{1}{\nu^{\prime
2}}\right)^{-p}\right)e^{-\frac{\nu^{\prime 2}}{2}},$ (9)
$A(p)=\frac{1}{1+\pi^{-1/2}2^{-p}\Gamma(\frac{1}{2}-p)}.$ (10)

Although $f(\nu)$ was originally constructed as a universal function that has
no cosmological evolution or power spectrum shape at all, its form can still
change drastically when baryonic feedback or exotic dark matter models are
considered, as has been discussed in the introduction. For instance, the
parameter $q$, which scales with the peak threshold $\nu$, is relevant to
baryonic feedback (van Daalen et al., 2020) and neutrinos (Costanzi et al.,
2013), as both effects change the peak threshold for spherical collapse. On
the other hand, $p$ determines the abundance gradient from low-mass to high-
mass haloes. As shown by simulations of exotic dark matter models, this
gradient is sensitive to specific dark matter properties (e.g., Lovell, 2020;
Kulkarni & Ostriker, 2022). The best-fit values of $(p,q)$ given by Sheth &
Tormen (1999) are $p=0.3$ and $q=0.707$, which yields $A=0.322$. Another
approach to the Sheth-Tormen mass function is applied by Despali et al.
(2016), where $A$ is set to be a free parameter, which can affect the vanilla
two-halo term at large $k$ (Mead et al., 2020). In this paper, we will
consider several mathematically normalized ST-like HMF models given by
Equation 10 (labelled ST-X), as well as a few models that are not normalized
(labelled DG-X).

Finally, it is important to clarify the concept of the halo definition. The
original ST model uses the spherical overdensity group finder, based on the
virial mass definition (Tormen, 1998), which takes into account the evolution
of halos over cosmic time (Jenkins et al., 2001). This virial mass definition
was also used in the calibration of HMcode (Mead et al., 2015). Therefore, as
our primary objective is to check the consistency of the WL analysis
pipelines, we will adhere to the virial definition that has been used in
conjunction with the ST model and HMcode calibration.

###  2.3 The Halo Model and HMcode

In this paper, we propose a modification of the matter power spectrum by
allowing the Sheth-Tormen HMF parameters to vary. Our work has a similar
motivation to Amon & Efstathiou (2022), which uses a had-oc modification of
$P(k)$ at the small scale. Amon & Efstathiou (2022) introduced the free
parameter $A_{\rm mod}$ such that:

$P(k)=P^{\rm L}(k)+A_{\rm mod}[P^{\rm NL}-P^{\rm L}]$ (11)

where $P^{\rm L}(k)$ and $P^{\rm NL}$ are the linear and non-linear power
spectra respectively. Using the KiDS-1000 data, they found that $A_{\rm
mod}=0.69\pm 0.1$ best fit the data when forcing the Planck18 cosmology. This
corresponds to a power loss of the order of $\sim 30\%$ in the non-linear
scales, in addition to the fiducial KiDS-1000 fiducial cosmology. The
underlying causes for this power loss could be the superposition of many
different effects. In our work, the proposed modification of $P(k)$ consists
of an alteration of the HMF, which is justified by the fact that there are
several clearly identified reasons for the HMF to be modified compared to its
vanilla functional form, as discussed in the Introduction.

For our calculations, we will use HMcode (Mead et al., 2016). In order to
solve the problem of missing power at intermediate scales with the halo model,
HMcode implements a modified one-halo $P^{\prime}_{\text{1H}}$ and two-halo
$P^{\prime}_{\text{2H}}$ terms such that:

$P(k)=\left[(P^{\prime}_{\text{1H}}(k))^{\alpha}+(P^{\prime}_{\text{2H}}(k))^{\alpha}\right]^{\frac{1}{\alpha}},$
(12)

where $\alpha$ is a smoothing parameter that depends on the scale where the
halo collapse occurs. The modified one-halo term takes into account the
damping from the halo exclusion effect (Smith et al., 2007). Specifically,
$P^{\prime}_{\mathrm{1H}}=\left[1-e^{-(k/k_{*})^{2}}\right]P_{\mathrm{1H}}$,
with $k_{*}$ being the one-halo damping wavenumber given by Mead et al.
(2016). The two-halo term is damped at quasi-linear scales following the
calculation of Crocce & Scoccimarro (2006), which means that only
$P^{\prime}_{\text{1H}}$ in HMcode depends on $p$ and $q$, not the 2-halo term
(unlike the vanilla halo model).

Within the halo model, the one-halo term is given by:

$P_{\text{1H}}(k)=\frac{1}{\bar{\rho}}\int_{0}^{\infty}\text{d}\nu\
M(\nu)\tilde{\rho}^{2}(\nu^{\eta}k,M)f(\nu),$ (13)

where, $M$ represents the halo mass and $\eta$ is the halo bloating parameter.
The normalized Fourier transform of the halo density profile is given by
$\tilde{\rho}^{2}(k,M)$, which was calibrated to the Navarro-Frenk-White (NFW)
profile. The NFW profile is defined by the equation:

$\rho(r)\propto\frac{r_{\mathrm{s}}}{r}\left(1+\frac{r}{r_{\mathrm{s}}}\right)^{-1},$
(14)

where $r_{\mathrm{s}}=r_{\mathrm{virial}}/c(M,z)$ is the scale radius and
$c(M,z)$ is the halo concentration. The HMcode uses the halo concentration
recipe developed by Bullock et al. (2001). It is worth noting that the halo
concentration relation can have a significant impact on small scales and may
have an internal degeneracy with HMF parameters. However, in the HMcode
framework, changing the halo concentration recipe can result in a worse $P(k)$
fit to simulations compared to just changing the HMF recipe, something which
has already been tested in the HMcode. Therefore, we keep the original Bullock
et al. (2001) halo concentration relation in the HMcode.

$c(M,z)=B\frac{1+z_{\mathrm{HF}}(M)}{1+z},$ (15)

Here, the halo-forming redshift $z_{\mathrm{HF}}(M)$ can be found in Bullock
et al. (2001). The parameter $B$ represents the minimum halo concentration and
is a useful indicator for baryonic feedback. A smaller value of $B$
corresponds to a stronger cooling or a larger AGN feedback (van Daalen et al.,
2020; Mead et al., 2021). This parameter was added to the KiDS-VIKING-450
(Hildebrandt et al., 2020) and KiDS-1000 pipelines (Asgari et al., 2021) with
a flat prior and a very limited range.

##  3 Data and Likelihood

Two lensing surveys are used in this work: the Kilo-Degree Survey fourth data
release (KiDS-1000 Kuijken et al., 2019; Giblin et al., 2021) and the Dark
Energy Survey year 3 data release (DES-y3 DES Collaboration et al., 2022). In
addition to those two datasets, we also use the LSST mock lightcone from the
Scinet Light Cone Simulations (SLICS Harnois-Déraps et al., 2018) to provide
comparisons between the current state-III results and the future Stage-IV
surveys.

###  3.1 KiDS

The Kilo-Degree Survey fourth data release dataset (KiDS-1000) covers an
unmasked field of view of $1006$ deg2. This survey is based on data obtained
from two telescopes, both of which are part of the Very-Large Telescope (VLT).
The main survey telescope, the $2.65$m VLT Survey Telescope (VST), provides
$4$ optical bands ($ugri$) with the VST-OmegaCAM (Kuijken, 2011), which has a
mean seeing of $0.7^{\prime\prime}$ in the $r$-band. Additionally, the $4.1$m
Visible and Infrared Survey Telescope for Astronomy (VISTA) provides another
$5$ near-infrared bands ($ZYJHK_{s}$) through the VISTA Kilo-degree Infrared
Galaxy survey (VIKING, Edge et al., 2013), which were introduced in the
updated sample KiDS-VIKING-450 of the third release of the same survey,
KiDS-450 (de Jong et al., 2017; Hildebrandt et al., 2017). The VIKING infrared
bands reach depths of $r\leq 25$ and cover redshifts up to $z\leq 6$.

The KiDS galaxy catalogue is divided into five photometric-redshift bins. The
photometric-redshifts used in KiDS-1000 were obtained by applying the Bayesian
photometric-redshift (BPZ) technique (Benítez, 2000; Raichoor et al., 2014),
and were refined with the self-organising map (SOM) method (SOM, Kohonen,
2001; Geach, 2012; Wright et al., 2020), using spectroscopic data for
calibration. The redshift distribution for each bin was determined based on
the spectroscopic redshifts of similar galaxies. Any galaxy that could not be
grouped with spectroscopic galaxies were removed from the ’Gold’ catalogue.
The accuracy of these redshift distributions was confirmed with the mock data
from the MICE2 simulation (Fosalba et al., 2015b, a; Crocce et al., 2015) (van
den Busch et al., 2022). The fiducial data products and corresponding
covariance matrices used in our study are taken from Asgari et al. (2021);
Giblin et al. (2021) and Joachimi et al. (2020), respectively, as official
KiDS-1000 samples.

###  3.2 DES

The Dark Energy Survey Year 3 data release (DES-y3) is shallower and wider
than KiDS-1000, covering $5000$ deg2 of the sky. DES employs the Dark Energy
Camera (DECam) (DECam Flaugher et al., 2015), which is mounted on the 4m
Victor M. Blanco Telescope located at Cerro-Tololo Observatory. The survey
comprises of five bands: four optical bands ($griz$) and one infrared band
($Y$), although, effectively only three bands (riz) are used in the shapes and
photometric redshift estimates. The survey’s depth is $i\leq 24$, which
roughly corresponds to $r\leq 24.4$. The average seeing is approximately
$1^{\prime}$ in the $r$-band (Sevilla-Noarbe et al., 2021).

For our tomographic re-analysis of DES-y3, we use the same four redshift bins
and data vector as in DES-y3 cosmic shear (Amon et al., 2022b; Secco et al.,
2022). The DES photometric redshifts are determined using the calibration
procedure outlined in Myles et al. (2021). This procedure involves three
estimates for $n(z)$: Bayesian probability distributions obtained from the
self-organizing map (SOMPZ) (Hartley et al., 2022), whose Gaussian prior is
derived from the combination of distance constraints from clustering redshifts
(Newman, 2008; Gatti et al., 2022) and shear ratios (Jain & Taylor, 2003;
Sánchez et al., 2022). All three methods are independent of each other. We use
the covariance matrices from the fiducial DES-y3 Fang et al. (2020) analyses
based on the CosmoLike framework (Krause & Eifler, 2017).

###  3.3 SLICS LSST Mock Lightcone

Scinet Light Cone Simulations (SLICS) (Harnois-Déraps & van Waerbeke, 2015) is
a large set of simulations designed for data quality assessment and covariance
estimation of multiple weak lensing data analyses, including the KiDS-1000 and
LSST. It is based on a series of $1025$ N-body simulation boxes using the
best-fit of WMAP9 + BAO + SN parameters in Hinshaw et al. (2013)
($\Omega_{m}=0.2905$, $\Omega_{b}=0.0473$, $h=0.6898$, $\sigma_{8}=0.826$,
$n_{s}=0.969$). In order to perform the mock analysis of future Stage-IV
surveys of LSST, we use all $819$ LSST-like Mock lightcones from SLICS
(Harnois-Déraps et al., 2018) to construct our covariance matrices for all 10
LSST-like tomographic bins, as well as the corresponding 10-bin tomographic
redshift distributions.

###  3.4 Likelihood

Figure 1: Left Panels: The halo mass functions at $z=0$ are shown with
different HMF input parameter. The dotted line in the top-left sub-panel shows
the canonical ST model with $q=0.707$, $p=0.3$, ${\sigma_{8}}=0.81$, and
$A(p)=0.322$. The solid and dashed lines on this panel show the impact
changing $\sigma_{8}$ (corresponding to a 1$\sigma$ uncertainty in Planck
Collaboration et al. (2020)), while keeping everything else unchanged. The
other sub-panels on the left show the same lines as in the top-left sub-panel
with two additional lines (a solid and a dashed) where some HMF parameters are
changed. The modified values for $p$, $q$ and $A$ correspond to the
uncertainty on these parameters from Despali et al. (2016).
Right Panels: The HMcode 3-D matter power spectra at $z=0$ are shown with
different HMF input parameters. The layout is the same as in the left panels,
except that the right panels show the 3-D matter power spectra.

Figure 2: Impact of modified HMF parameters on the shear correlation function
$\xi_{+}$ (left panels) and $\xi_{-}$ (right panels). The layout and the
curves is the same as in Figure 1. The Limber projection is based on the
redshift distribution of KiDS-1000 bin $5$.

We use the Code for Anisotropies in the Microwave Background (CAMB; Lewis &
Challinor 2011) to compute the matter power spectrum $P_{\delta}$. The non-
linear power spectra is calculated with our modified version of HMcode-2016
(Mead et al., 2016) in order to include the HMF parameters $p$ and $q$ as free
parameters.

The parameter space is sampled with the multinest algorithm, which is
integrated in the Cosmosis framework (Zuntz et al., 2015). multinest is a
method that integrates over the posterior probability space, providing good
precision and efficiency when dealing with high-dimensional parameter spaces.
In order to characterize our chain outputs, we use the multinest-weighted mean
instead of the maximum posterior, as it better depicts the Bayesian nature of
the posterior distribution. It is worth noting that multinest may
underestimate parameter errors, as shown by Lemos et al. (2022) where the
error on $S_{8}$ can be underestimated by 5-10%.

Each survey has its own specific redshift calibration, source selection, and
corresponding scale cuts. For this reason, it is particularly interesting to
compare how these two data sets will constrain the modified HMF. We use the
$\Lambda$CDM optimised scale cut for DES-Y3. For KiDS-1000, we adopt the
fiducial scale cut
($\theta(\xi_{+})\in[0.5^{\prime},300^{\prime}],\theta(\xi_{-})\in[4.0^{\prime},300^{\prime}]$)
(i.e. the one used by Asgari et al. (2021)). In Appendix A, we also provide
additional KiDS-1000 scale cuts comparisons. To sample the HMF and profile
parameters, we use the following flat priors:

$\displaystyle 2.0\leq$ $\displaystyle B$ $\displaystyle\leq 3.13\text{\
\cite[citep]{(follows \@@bibref{AuthorsPhrase1Year}{Asgari+21}{\@@citephrase{,
}}{})}},$ (16) $\displaystyle-0.5<$ $\displaystyle p$ $\displaystyle<0.5,$
(17) $\displaystyle 0.0<$ $\displaystyle q$ $\displaystyle<2.0,$ (18)
$\displaystyle 0.0<$ $\displaystyle A$ $\displaystyle\leq 0.55,\text{\
unnormalised\ HMF\ only}.$ (19)

For KiDS-1000, we use the same setup as Asgari et al. (2021), except for the
HMcode and HMF parameters. For DES-y3, we made two modifications: we use
HMcode and Multinest instead of Halofit and Polychord (Amon et al., 2022a;
Secco et al., 2022). In Figure 9 of Appendix A, we compare the DES-Y3
posterior contours between these two setups. The difference is small, but it
is worth to mention, even if it is irrelevant to our discussion in this paper.
DES and KiDS collaboration et al. (2023) also found a similar difference
between the two sampling codes. All other modelling choices, such as intrinsic
alignments333In KiDS-1000 (Asgari et al., 2021), the intrinsic alignment model
is the redshift-independent NLA model, while DES-y3 (Amon et al., 2022b; Secco
et al., 2022) uses the TATT model (Blazek et al., 2019)., are kept unchanged
compared to the original KiDS-1000 and DES-Y3 analyses.

Note that the amplitude $A$ of the HMF is occasionally considered as a free
parameter, as initially proposed by Despali et al. (2016). The question arises
whether we should also set $A$ free in addition to ($p$, $q$). In the appendix
of Mead et al. (2021), it was demonstrated that abandoning the normalization
requirement of the HMF does not affect the normalization of the one-halo term,
implying that $A$ could be estimated if a sufficiently large range of scales
is observed in addition to ($p$, $q$). We attempted to constrain ($p$, $q$,
$A$) simultaneously, but discovered that $A$ and $q$ are strongly degenerate,
and $p$ cannot be constrained. As a result, we employ A(p) from equation 10 to
normalize the HMF in our main analysis.

##  4 Results

###  4.1 Impact of the halo mass function parameters on the halo abundance,
$P(k)$ and shear correlation functions.

Cosmology | ${\Omega_{\mathrm{m}}}$ | ${S_{8}}$ | $h$ | $n_{s}$ | $M_{\nu}$ | $\Omega_{b}h^{2}$ | Reference
---|---|---|---|---|---|---|---
KiDS-1000 | 0.253 | 0.760 | 0.729 | 0.938 | 0.06 | 0.022 |  Asgari et al. (2021) $\xi_{\pm}$ posterior mean
DES-y3 | 0.289 | 0.772 | 0.677 | 0.966 | 0.06 | 0.022 |  Amon et al. (2022b); Secco et al. (2022) $\Lambda$CDM-Optimised $\xi_{\pm}$ posterior mean
Planck 2018 | $0.316$ | $0.832$ | 0.677 | 0.966 | 0.06 | 0.022 |  Planck Collaboration et al. (2020) TT,TE,EE+lowE+lensing+BAO posterior mean
Table 1: Cosmological parameters set-up of our study. All other independent cosmological parameters such as $N_{\nu,\mathrm{eff}}$ are fixed to the fiducial analysis and dependent cosmological parameters such $\Omega_{\Lambda}$ are varied with these parameters in order to keep a flat cosmology. In the case of DES, they did not report certain cosmological parameters, we used the Planck values for $n_{s}$ and $h$, since DES did not report what values they used. # | $p$ | $q$ | $A$ | Note
---|---|---|---|---
ST1 | $0.300$ | $0.707$ | A(p) |  Sheth & Tormen (1999)
ST2 | $0.360$ | $0.880$ | A(p) | A HMF parameter fit of the Model-ST1 HMF under the $200\rho_{c}$ Halo Mass Definition
ST3 | $0.402$ | $0.903$ | A(p) | Combination of ST4 and ST5
ST4 | $0.402$ | $0.707$ | A(p) | Best Planck Collaboration et al. (2020) ${S_{8}}$ recovery when fixing $q$ in KV450
ST5 | $0.300$ | $0.903$ | A(p) | Best Planck Collaboration et al. (2020) ${S_{8}}$ recovery when fixing $p$ in KV450
ST6 | $0.000$ | $1.700$ | A(p) | Best Planck Collaboration et al. (2020) ${S_{8}}$ recovery when fixing $p=0$ in KV450
ST7 | $0.237$ | $0.614$ | A(p) | A random-chosen high amplitude HMF
DG1 | $0.258$ | $0.766$ | 0.330 | Original Despali et al. (2016) of virial haloes
DG2 | $0.258$ | $0.766$ | 0.250 | Parameters giving a $25\%$ Lower HMF Amplitude
DG3 | $0.258$ | $0.766$ | 0.410 | Parameters giving a $25\%$ Higher HMF Amplitude
ST-Free | Free | Free | A(p) | A flat HMF parameter prior based on the Sheth & Tormen (1999) Model
DG-Free | Free | Free | Free | A flat HMF parameter prior based on the Despali et al. (2016) Model
Table 2: List of the various HMF models used in this paper. Models beginning
with ST are based on the Sheth & Tormen (1999) halo mass function models,
where A=A(p) indicates that the HMF is normalized using equation 10. Models
beginning with DG are based on the Despali et al. (2016) HMF models, where the
HMFs are not necessarily normalized. The reasons why we chose these particular
HMF models are provided in the ‘Note’ column of this table. Figure 3: Plot of
the $q-p$ plane posterior contours for KiDS-1000 and DES-y3 data, with fixed
cosmology either to their fiducial cosmic shear constraints (salmon pink for
KiDS-1000 and cyan for DES-y3, shown as dashed line filled contours) or to the
Planck 2018 (dark orange for KiDS and bright azure for DES, shown as solid
line unfilled contours). All posteriors are sampled using multinest. The
canonical ST parameters are indicated with a marker. Figure 4: Posterior
contours ($68\%$ and $95\%$) in the ${\Omega_{\mathrm{m}}}-{\sigma_{8}}$,
${\Omega_{\mathrm{m}}}-{S_{8}}$ and ${\sigma_{8}}-{S_{8}}$ planes for
KiDS-1000, DES-y3 and Planck. The salmon-pink filled contours (Sheth & Tormen
(1999) halo mass function, labelled Model-ST1) corresponds to the fiducial
analysis of Asgari et al. (2021). The solid dark orange contour (Model-ST5)
correspond to $(p,q)=(0.3,0.903)$ (See table 2), the best Planck Collaboration
et al. (2020)-matching model in our table of halo mass function models. The
filled cyan (DES-y3 with Model-ST1) and solid bright azure (DES-y3 with Model-
ST5) contours are posteriors obtained under the same logic from the DES-y3
pipeline with HMcode and multinest. The slate grey filled contours indicate
the $95\%$ posterior of Planck 2018 TTTEEE+LowE+Lensing+BAO result (Planck
Collaboration et al., 2020) as the comparison.

As discussed in Section 2.2, the HMF in Sheth & Tormen (1999) is constructed
using two parameters, $p$ and $q$, and an additional parameter $A$ if the
normalisation is not fixed. In this section, we investigate the impact of a
change of these parameters on the abundance of halos $n(M)$, the matter power
spectrum $P(k)$ and the shear correlation functions $\xi_{\pm}$. This will
allow us to visualise how the HMF degrees of freedom impact the quantities
important for weak lensing calculations, and to what extent they might be
degenerate with $\sigma_{8}$ 444For weak lensing, ${S_{8}}$ would be the
parameter to look at, but when $\Omega_{m}$ is fixed, $\sigma_{8}$ and
${S_{8}}$ are the same.

Figure 1 (left panels) shows the change in the HMF relative to the fiducial
Sheth-Torman HMF (q = 0.707, p=0.3, A(p)=0.322) and Planck2018 cosmology when
either $\sigma_{8}$ or Sheth-Torman parameters values are varied. The top-left
sub-panel serves as a reference, where we show the matter power spectrum for
three different values of $\sigma_{8}$ by $\pm 1\sigma$ using the Planck
Collaboration et al. (2020) as a fiducial value of 0.81 (horizontal dotted
line). The corresponding dotted, solid and dashed curves are shown in all
other left sub-panels for comparison to changes when the HMF parameters are
varied. In these other left sub-panels, additional curves are shown with the
HMF parameters varied within the uncertainty reported in Despali et al.
(2016). The second left sub-panel shows that varying $q$ by $\pm 0.05$ has a
similar effect to $\sigma_{8}$, indicating that $q$ is partially degenerate
with $\sigma_{8}$ as expected from Eqs. 9 and 12, since, in our formalism,
only the 1-halo term depends on the HMF parameters. One can see that $q$ and
$\sigma_{8}$ have a larger impact on the abundance of haloes with
$M>10^{13}M_{\odot}$. In the third left sub-panel, $p$ is varied by $0.009$
along with $A$ to keep the HMF normalised, while in the fourth left sub-panel,
$A$ is fixed to its fiducial value while $p$ is varied by the same amount.
Comparing these last two sub-panels, we observe that varying or fixing $A$ has
a significant effect on the HMF.

When $A$ is fixed, varying $p$ primarily impacts the abundance of low-mass
haloes. When $A$ is used to normalise the HMF, haloes of all masses are
impacted, with more sensitivity towards higher masses. In this case a lower
value of $p$ increases the abundance of haloes. Therefore, $p$ mostly impacts
the abundance of halos with virial masses lower than $10^{13}$ solar masses
(when $A$ is fixed), while $q$ governs the abundance of higher-mass-halos
($M_{\mathrm{Halo}}>10^{13}M_{\odot}$).

The right panels of Figure 1 show the variation in the HMcode matter power
spectrum, using a similar plotting scheme as in the left panels. Changing the
HMF parameters $p$ and/or $q$ results in a notable shape change of the matter
power spectrum, particularly around the $k\sim 1$ scale. Generally, in HMcode,
the alteration of the halo models only affects the power spectrum at $k>0.01$
because the two-halo term remains unaffected (see Eq. 12). The second right
sub-panel demonstrates that a higher $q$-value produces suppression of the
power spectrum at the cluster scale ($k\sim 1$) and a rise of power at small
scales, resembling the baryonic feedback models (Schneider & Teyssier, 2015).
It is interesting to note that $N$-body simulations favour a higher $q$-value
compared to hydrodynamical simulations (Bocquet et al., 2016; Schaye et al.,
2023). On the other hand, increasing the parameter $p$ produces a power
spectrum change up to $k=10\mathrm{h\ Mpc}^{-1}$ that more closely resembles a
higher AGN feedback (Chisari et al., 2018). This is due to the fact that the
dark matter clustering at this scale is governed by the gas temperature of
baryons (Beltz-Mohrmann & Berlind, 2021). Additionally, alternative dark
matter models such as Axions have similar effects on this scale (e.g., Marsh,
2015).

Figure 2 shows the shear correlation functions $\xi_{\pm}$ (equation 2),
following a similar layout as the panels in Figure 1 (The left panels of
Figure 2 are for $\xi_{+}$ and the right panels for $\xi_{-}$). It is observed
that for angular scales $\theta<30^{\prime}$, changes in $p$ or $q$ have a
similar effect as $\sigma_{8}$, but different scale dependence between
$\xi_{+}$ and $\xi_{-}$. Here, all angular scales are affected, and changing
$q$ has a more significant impact. This is the consequence of the fact that
the shear correlation function is mixing all physical scales in the
projection, where the clear separation between the 1-halo and the 2-halo terms
in Fourier space is lost in the configuration space.

The comparison between the power spectra, the shear auto-correlation functions
and the HMF (from figure 2 and 1), shows that modifying the HMF parameters
leads to changes in these quantities which $\sigma_{8}$ cannot reproduce. This
means that the HMF parameters and $\sigma_{8}$ are not degenerate, so we
expect that the HMF parameters can be constrained in addition to $\sigma_{8}$.
In the next section, we will explore quantitatively the cosmological
parameters and the HMF parameters constraints with the DES and KiDS surveys.

We now present the results of our re-analysis of KiDS-1000 and DES-y3 data
when both the halo mass function and cosmological parameters are permitted to
vary either independently or jointly.

Figure 5: Marginalised posterior contours for a flat $\Lambda$CDM model with the normalised Sheth & Tormen (1999) halo mass function model with the normalisation from the Equation 10. The solid dark red and indigo blue contours represent the results obtained by the KiDS-1000 and HMcode-based DES-y3 pipeline, respectively. The slate grey contours indicate the $95\%$ posterior of the Planck 2018 TTTEEE+LowE+Lensing+BAO result (Planck Collaboration et al., 2020) for comparison. We have also added the cosmological contours with fixed HMF Model-ST1 from Figure 4 and the HMF constraints with fixed Planck 2018 cosmology from Figure 3 into this figure with their original colours as a comparison. The markers represent the fitted HMF parameters $p$ and $q$ from Sheth & Tormen (1999), with the normalisation A(p) (equation 10). Model | KiDS $p$ | KiDS $q$ | KiDS $A$ | DES $p$ | DES $q$ | DES $A$ | KiDS ${S_{8}}$ | KiDS $\chi^{2}$ | DES ${S_{8}}$ | DES $\chi^{2}$
---|---|---|---|---|---|---|---|---|---|---
ST1 | 0.300 | 0.707 | A(p) | 0.300 | 0.707 | A(p) | $0.765_{-0.020}^{+0.019}$ | 259.1 | $0.799_{-0.016}^{+0.016}$ | 282.7
ST2 | 0.360 | 0.880 | A(p) | 0.360 | 0.880 | A(p) | $0.863_{-0.021}^{+0.024}$ | 257.8 | $0.867_{-0.023}^{+0.027}$ | 285.5
ST3 | 0.402 | 0.903 | A(p) | 0.402 | 0.903 | A(p) | $0.920_{-0.031}^{+0.030}$ | 269.6 | $0.877_{-0.027}^{+0.027}$ | 288.2
ST4 | 0.402 | 0.707 | A(p) | 0.402 | 0.707 | A(p) | $0.856_{-0.022}^{+0.024}$ | 257.6 | $0.852_{-0.025}^{+0.024}$ | 284.5
ST5 | 0.300 | 0.903 | A(p) | 0.300 | 0.903 | A(p) | $0.826_{-0.019}^{+0.021}$ | 256.5 | $0.851_{-0.021}^{+0.020}$ | 284.0
ST6 | 0.000 | 1.700 | A(p) | 0.000 | 1.700 | A(p) | $0.874_{-0.027}^{+0.027}$ | 273.2 | $0.887_{-0.022}^{+0.027}$ | 286.1
ST7 | 0.257 | 0.614 | A(p) | 0.257 | 0.614 | A(p) | $0.702_{-0.018}^{+0.018}$ | 260.3 | $0.738_{-0.012}^{+0.011}$ | 282.4
DG1 | 0.258 | 0.766 | 0.330 | 0.258 | 0.766 | 0.330 | $0.778_{-0.020}^{+0.019}$ | 258.9 | $0.812_{-0.017}^{+0.017}$ | 282.8
DG2 | 0.258 | 0.766 | 0.250 | 0.258 | 0.766 | 0.250 | $0.823_{-0.021}^{+0.021}$ | 256.9 | $0.844_{-0.020}^{+0.020}$ | 284.0
DG3 | 0.258 | 0.766 | 0.410 | 0.258 | 0.766 | 0.410 | $0.748_{-0.022}^{+0.022}$ | 264.4 | $0.783_{-0.015}^{+0.014}$ | 282.7
Lensing | $0.11^{+0.24}_{-0.09}$ | $0.91^{+0.16}_{-0.24}$ | A(p) | $-0.01^{+0.28}_{-0.18}$ | $0.91^{+0.16}_{-0.17}$ | A(p) | 0.760 | 258.9 | 0.772 | 283.3
| $-0.03^{+0.28}_{-0.28}$ | $0.83^{+0.13}_{-0.18}$ | $0.34^{+0.07}_{-0.10}$ | $-0.01^{+0.26}_{-0.30}$ | $0.77^{+0.14}_{-0.14}$ | $0.35^{+0.10}_{-0.09}$ | 0.760 | 258.6 | 0.772 | 283.0
Planck18 | $0.07^{+0.27}_{-0.12}$ | $1.23^{+0.24}_{-0.33}$ | A(p) | $-0.15^{+0.17}_{-0.26}$ | $1.25^{+0.20}_{-0.16}$ | A(p) | 0.832 | 263.4 | 0.832 | 284.8
| $-0.02^{+0.27}_{-0.34}$ | $1.10^{+0.18}_{-0.23}$ | $0.35^{+0.08}_{-0.10}$ | $-0.04^{+0.21}_{-0.36}$ | $1.06^{+0.16}_{-0.16}$ | $0.42^{+0.11}_{-0.06}$ | 0.832 | 263.3 | 0.832 | 285.0
ST-Free | $0.10^{+0.26}_{-0.10}$ | $0.93^{+0.19}_{-0.27}$ | A(p) | $-0.03^{+0.27}_{-0.22}$ | $0.93^{+0.21}_{-0.21}$ | A(p) | $0.771_{-0.046}^{+0.045}$ | 258.3 | $0.767_{-0.033}^{+0.029}$ | 282.2
DG-Free | $0.01^{+0.27}_{-0.28}$ | $0.87^{+0.18}_{-0.22}$ | $0.37^{+0.08}_{-0.10}$ | $0.00_{-0.27}^{+0.25}$ | $0.76_{-0.17}^{+0.17}$ | $0.36^{+0.11}_{-0.01}$ | $0.782_{-0.053}^{+0.047}$ | 259.7 | $0.761_{-0.033}^{+0.031}$ | 282.6
Table 3: Table of posterior weighted-mean values of ${S_{8}}$ for each HMF
model, along with the HMF parameters for each cosmological model, along with
their corresponding marginal errors. The values are reported using weighted
quantiles of $16\%$ and $84\%$ percentile of each parameter via the Getdist
package (Lewis, 2019). For the HMF models, A=A(p) indicates that the HMF is
normalized via equation 10. Fixed parameters in the chain are reported without
error bars. The $\chi^{2}$ value in the table corresponds to the minimum
$\chi^{2}$ from the chain and may be strongly influenced by the randomness in
the multinest sampler.

###  4.2 Halo mass function constraints for fixed cosmology

Here, we show the constraints on ($p$, $q$), for the fiducial and Planck 2018
cosmologies. For the KiDS and DES fiducial cosmologies, the mean posterior
cosmology is assumed, which are $\xi_{\pm}$ posterior mean of Asgari et al.
(2021) for KiDS-1000, Fiducial $\xi_{\pm}$ posterior mean from Amon et al.
(2022b); Secco et al. (2022) for DES-y3 data. For the Planck 2018 cosmology,
we use the TT,TE,EE+lowP+lensing+BAO posterior mean. Although the cosmology is
fixed for each data set, the nuisance parameters are still marginalized over,
using the procedures described in Asgari et al. (2021), Amon et al. (2022a),
and Secco et al. (2022). Table 1 presents the cosmological parameters that
were used.

The results are presented in Figure 3. Despite the degeneracy between $p$ and
$q$, the standard values of the Sheth-Tormen (ST) model, denoted as Model-ST1,
are still consistent with the lensing fiducial cosmologies of each survey.
This is expected, it is shown here as a self-consistency test that the
fiducial cosmology is not in tension with the ST1 values of $(p,q)$ even if
they are allowed to vary. However, when the Planck cosmology is assumed, the
ST1 values are rejected at more than $2\sigma$. The Planck cosmology generally
prefers a higher $q$ value. In the previous section, we showed that a higher
$q$ corresponds to a lower abundance of high-mass haloes compared to ST1. One
can show that, under the Planck cosmology, and for halo mass
$M>10^{14}M_{\odot}$, the posterior HMF has $29.5\%_{-8.5\%}^{+8.5\%}$ less
mass in DES and $48.8\%_{-9.4\%}^{+8.3\%}$ less mass in KiDS compared to a
standard Sheth & Tormen (1999) HMF. The difference arises because the Planck
cosmology predicts more high-mass haloes than what the weak lensing fiducial
cosmologies suggest, a tension that could be addressed by changing the HMF
parameters. While the amplitude of the matter power spectrum ${S_{8}}$ affects
all scales of $P(k)$, the HMF parameters $p$ and $q$ only affect non-linear
scales $k>1$. This result is consistent with the DES-y1 cluster abundance
analysis (Costanzi et al., 2019; Abbott et al., 2020), which suggests that a
richer mass abundance prediction requires a lower value of ${S_{8}}$ to
compensate for the observation. One also notes that, while for the KiDS and
DES-Y3 fiducial cosmologies the $(p,q)$ contours encompass the ST1 values, the
same contours shift to higher $q$ and tend to split from each other with the
Planck cosmology. It means that, in order to be consistent with the Planck18
cosmology, the best fit $P(k)$ model must be different between KiDS-1000 and
DES-y3. This is of course an undesirable feature, and it is an illustration of
how our approach can also be used as a self-consistency check between
different surveys. It will be interesting to see how stage IV surveys (e.g.
LSST versus Euclid) can perform on this test. In section 4.4 we will show how
precisely a LSST-type survey can measure $(p,q)$.

Figure 6: The (${S_{8}}$,$\Omega_{m}$) forecast for a LSST-like survey with
corresponding Field-of-View and galaxy number density. The darkred contours
correspond to a joint HMF and cosmology likelihood using an LSST-like setup.
The limegreen and gold contours correspond to the LSST-like cosmological
analysis with fixed $(p,q)=(0.3,0.707)$ and HMF-only analysis with a fixed
cosmology the same as the SLICS simulation (Harnois-Déraps et al., 2018)
respectively. These three sets of LSST-like contours use the Pyccl-generated
theoretical mock data vectors and LSST-like covariance matrices computed from
the SLICS simulation (Harnois-Déraps et al., 2018). We also add both KiDS-1000
and DES-y3 joint HMF-cosmology contours here only in order to compare the
constraining power between the Stage-III and the Stage-IV surveys.

###  4.3 Cosmological constraints for fixed halo mass function parameters

In order to test if a unique $n(M)$, common to KiDs and DES, is able to give
consistent cosmological results (and whether the posterior cosmology might be
consistent with Planck’s or not), it would be useful to perform a series of
tests with a set of fixed HMF parameters different from ST1. The list of HMF
parameters, labelled with numbers, that we will use are given in table 2.
These models were determined during the exploration phase of our work using
the KiDS-VIKING-450 data. Most of these models have higher $p$ and/or $q$
values than the original ST1 model. They were chosen so that $P(k,z=0)$ at
$k\geq 1$ h Mpc-1 is lowered by approximately $20\%$ to $30\%$ compared to the
ST1 values, similar to some of the baryonic feedback models, such as gas
outflows (Schneider et al., 2020), counterbalancing the effect of a higher
${\sigma_{8}}$. This also corresponds to a deficit of roughly $30\%$ to $50\%$
of the high-mass end of the halo mass function
($M_{\mathrm{Halo}}>10^{14}M_{\odot}$). The marginal errors on $(p,q)$
reported in Table 3 correspond to $16\%$ and $84\%$ percentiles of each
parameter calculated with the Getdist package (Lewis, 2019).

Figure 4 shows the contours for ST1 and ST5 (ST5 is chosen as it best matches
the Planck ${S_{8}}$ from our list in Table 2). The solid blue and red
contours show the cosmological constraints using ST1. The open contours show
how the cosmological parameters posteriors are displaced with model ST5. As
expected, with ST5, the change of HMF parameters lowers the abundance of high-
mass haloes ($M_{\mathrm{Halo}}>10^{13}M_{\odot}$) and therefore increases the
inferred amplitude of $P(k)$]. Figure 4 shows that looking at the shifts from
the solid to the line contours, the $(p,q)$ degrees of freedom can easily lead
to a few percent changes in the cosmological and HMF parameters. It is
important to note that the new degrees of freedom do not make the model ill-
constrained. This can be seen by looking at the $\chi^{2}$ values in Table 3,
they do not dramatically change between the different models. Models ST7 (see
Table 3), which uses $(p,q)$ values lower than the ST1 model, also provides a
good fit to both the KiDS-1000 and DES-y3 data, with a lower ${S_{8}}$, while
the reduced $\chi^{2}$ remains very close to the of the original cosmological
analysis with the HMF parameters fixed to ST1. Table 3 shows that, for the
models we have chosen, the spectral features introduced by the different
$(p,q)$ values are all acceptable by the data, which explains why, a change in
${S_{8}}$ can still absorb the changes caused by the different HMF parameters.
We will see in the next subsection that this degeneracy is broken with stage
IV surveys 555Note that adding external, low redshift, cosmological prior from
supernovae SNIa and Baryon Acoustic Oscillation would help in constraining
$(p,q)$. Our conclusion is that the KiDS-1000 and DES-y3 data have some
spectral feature in the $P(k)$, essentially a power loss at the small scale,
which a non-standard $(p,q)$ can absorb, but the constraining power of these
surveys is not strong enough to reject low $(p,q)$ values. By comparing
Figures 4 and 3 we clearly understand why the Planck18 cosmology would
definitely not work with the ST1 values: they are too low compared to the
level of power loss that exists in the data.

Figure 7: Left panel: This panel shows the reconstructed halo mass function
$n(M)$ obtained from the posteriors of this work using KiDS-1000 data,
relative to the fiducial $n(M)$ calculated from the fiducial cosmological
constraints of KiDS-1000 $\xi_{\pm}$ (labelled as ‘Fiducial’). The dashed line
labelled as ST1+Planck18 corresponds to the Planck18 cosmology with standard
Sheth-Tormen (ST) parameters $(p,q)=(0.3,0.707)$. The red contours correspond
to the posterior obtained under a fixed Planck18 cosmology where only the ST
parameters $(p,q)$ are allowed to vary, which corresponds to the dark orange
contour from figure 3. The bright azure contour labelled as "KiDS-1000 Cosmo
w. Fixed HMF ST1" corresponds to the fiducial cosmology chain where the
cosmological parameters are allowed to vary but with a fixed Sheth-Tormen HMF,
and this is represented by the salmon pink contour from figure 4. Finally, the
gold contour labelled as "KiDS-1000 Cosmo+HMF Flat Priors" corresponds to the
dark red contour of figure 5, where both HMF and cosmological parameters are
sampled under a flat prior. This plot is calculated for redshift $z=0$. Right
panel: The right panel shows the same information as the left panel, but
obtained from the DES-y3 data.

Figure 8: The same as Fig 7, but for the reconstructed matter power spectrum
$P(k)$ at redshift $z=0$ relative to the fiducial cosmologies and HMcode power
spectra (labelled ‘Fiducial’). The normalisation factor of
${\omega_{\mathrm{m}}}={\Omega_{\mathrm{m}}}h^{2}$ is used to represent the
effect of lensing projections.

###  4.4 Joint constraints on cosmological and Halo Mass Function parameters

Having learned how $(p,q)$ impacts $P(k)$, and the interplay with the
cosmological parameters, we can now explore the joint constraint on cosmology
and the HMF. Figure 5 shows the joint constraints on the HMF and cosmological
parameters, with the flat priors. The ${S_{8}}$ posterior of the two surveys
are in very good agreement, as well as the $(p,q)$ posteriors 666Remember that
it was not the case for the fixed Planck18 cosmology as the blue and orange
contours in Figure 5 show. The KiDS-1000 ${S_{8}}$ contour with flat HMF
parameters prior is similar to Figure 5 in Amon & Efstathiou (2022), which
shows that their $A_{\rm mod}$ parametrization and our modified HMF has a
similar effect on the power spectrum. On the other hand, the DES-y3 with flat
HMF prior tends to prefer a smaller $\sigma_{8}$ and larger $\Omega_{m}$
compared to their fiducial analysis, in a way which makes the ${S_{8}}$
posterior agrees with KiDS-1000. This is illustrated by the shift of
$(\Omega_{m},\sigma_{8})$ contours which are both still aligned with the
${S_{8}}$ degeneracy. How can we interpret this result? It is known that
moving along the ${S_{8}}$ degeneracy is particularly sensitive to the power
spectrum slope. By comparing the $\chi^{2}$ values of the HMF flat prior cases
to the fixed ST1 cases, for both KiDS-1000 and DES-y3 in Table 3, we clearly
see that the $(p,q)$ degrees of freedom do not create a tension in the model,
and, in fact, the $\chi^{2}$ are almost unchanged.

The meaning of Figure 5 is that with the HMF flat prior, the KiDS and DES
contours moved such that $(p,q)$ and ${S_{8}}$ posteriors are in better
agreement, but the distance between the $\Omega_{m}$ posteriors has increased.
This behaviour could suggest an internal inconsistency in the model, revealed
by the difficulty of finding a good agreement for the $(p,q),{S_{8}}$ and
$\Omega_{m}$ posteriors simultaneously. The DES and KiDS contours on Figure 5
remain large and, at this stage, it only suggests a trend about how the
contours want to move, rather than representing a definitive result. We can
however make the following comment: By introducing the HMF parameters, we have
given the possibility for the best fit to deviate from the conventional
theoretical wisdom $(p,q)=(0.3,0.707)$. If this model is right, then the fit
should recover the ST1 values within the errors, without creating tensions
elsewhere. Figure 5 only suggests that it might not be the case, but more
statistical power is necessary. The hybrid analysis in DES and KiDS
collaboration et al. (2023) leads to a reduction of the ${S_{8}}$ tension,
without completely eliminating it. It would be interesting to rerun their
analysis with our HMF flat prior setup and observe how the posteriors on
${S_{8}}$ and $\Omega_{m}$ want to move.

In order to see if future stage IV surveys could measure the HMF parameters,
we have performed a similar Multinest-based analysis on a Legacy Survey of
Space and Time, LSST-like, survey. By using unbiased and noiseless mock data
vectors generated by Pyccl (Chisari et al., 2019) for the fiducial cosmology
of SLICS, we set constraints on parameters by analysing this mock data vector
for fixed cosmological but variable HMF parameters, fixed HMF but variable
cosmological parameters, and free HMF and cosmological parameters. We show the
contours of these LSST-like contours, accompanied by current stage III
KiDS-1000 and DES-y3 joint constraints on cosmological and HMF parameters, in
Figure 6.

Figure 6 shows that a stage IV survey can easily probe the HMF parameters in
addition to probing cosmology. For instance, in the LSST-like set-up with flat
prior of HMF parameters, the ${S_{8}}$ constraint reads
$0.814_{-0.013}^{+0.012}$, while constraints of HMF parameters are
$q=0.741^{+0.049}_{-0.055}$ and $p=0.279^{+0.033}_{-0.032}$. Meanwhile,
KiDS-1000 gives ${S_{8}}=0.771_{-0.046}^{0.045}$, $q=0.93^{+0.19}_{-0.27}$,
and $p=0.10^{+0.26}_{-0.10}$; DES-y3 also gives a very similar constraint to
KiDS-1000 with ${S_{8}}=0.767_{-0.033}^{+0.029}$, $q=0.93^{+0.21}_{-0.21}$,
and $p=-0.03^{+0.27}_{-0.22}$ – stage III Errorbars on these parameters are
roughly $\sim 4\times$ larger than our prediction of stage IV LSST-like
surveys in both cosmological and HMF parameters. Note that we only used the
Sheth and Tormen model, but our approach can easily be generalized to any halo
model. We also notice that the widening of the LSST contours is significant
when the HMF parameters are not fixed, e.g., the $S_{8}$ constraint widened
from $S_{8}=0.813_{-0.004}^{+0.004}$, where the HMF parameters are fixed, to
$S_{8}=0.814_{-0.013}^{+0.012}$, where the HMF parameters are not fixed. This
is an important factor to consider for any kind of halo model being used for
cosmological inference with stage IV surveys in the future.

###  4.5 Full Posterior on the Halo Mass Function and $P(k)$

In order to conclude our study, we have to quantify the range of changes in
$n(M)$ and $P(k)$ corresponding to the cosmology contours in Figure 5. This
comparison is performed by taking the contours from Figures 3, 4 and 5 and
calculating the corresponding $1-\sigma$ regions for $n(M)$ and $P(k)$. The
results are shown on Figure 7 for $n(M)$ and Figure 8 for $P(k)$. Note that
$P(k)$ is scaled by ${\omega_{\mathrm{m}}}={\Omega_{\mathrm{m}}}h^{2}$ to
account for projection effects. For each figure, the left panel is for
KiDS-1000 and the right panel for DES-y3. In all panels, the horizontal line
represents the normalization to the fiducial cosmology with the standard ST
parameters. The black dashed line, on the other hand, shows the calculation
for the Planck18 cosmology. The comparison between the dashed line and the
horizontal dotted line gives another view of how the ${S_{8}}$ tension
manifests itself in $n(M)$ and $P(k)$. On all panels, the red shaded region
shows the $1\sigma$ posterior of the ($p$, $q$) parameters with Planck18
cosmology from Figure 3. The bright azure region corresponds to the $1\sigma$
cosmological contours of both KiDS and DES under a canonical Sheth & Tormen
(1999) (Figure 4). Finally, the gold contours represent the corresponding
$1\sigma$ region of Figure 5, where both cosmological and HMF parameters are
allowed to vary. All plots show calculations at redshift $z=0$.

Our analysis suggests that the DES $n(M)$ prefers a $29.5\%_{-8.5\%}^{+8.5\%}$
lower total mass and the KiDS $n(M)$ prefers a $48.8\%_{-9.4\%}^{+8.3\%}$
lower total mass for all halos with masses
$M_{\mathrm{Halo}}>10^{14}M_{\odot}$, compared to Planck Collaboration et al.
(2020) cosmological model with the standard Sheth & Tormen (1999) HMF
parameters. This agrees with recent studies on halo mass abundance (e.g., Li
et al. 2019; Abbott et al. 2020), which both suggest astrophysical
explanations for this low abundance. Regarding $P(k)$, the red region in
Figure 8 shows a power spectrum suppression by $\sim 25\%$ at $k\simeq
1\mathrm{h\ Mpc}^{-1}$ for the Planck18 cosmology. This is a very strong
amplitude suppression because it comes in addition to the internal
modifications of the halo model in HMcode, which accounts for the baryonic
effect. Interestingly, our result is similar to the value $A_{\rm mod}\sim
0.7$ obtained from Equation 11 found by Amon & Efstathiou (2022). In other
words, the unrealistic high Active Galactic Nuclei (AGN) temperature reported
by Amon & Efstathiou (2022) is also found with the non-AGN version of HMcode
when the HMF parameters are free to vary. Galaxy-galaxy lensing (GGL)
encounters a similar issue, known as the ’Lensing-is-Low’ problem. Leauthaud
et al. (2017) found that the measured surface mass density contrast
$\Delta\Sigma_{\mathrm{data}}$ at small scale is approximately 20% smaller
than the theoretical predictions if we assume the Planck18 cosmology. In
comparison, for KiDS-1000, (Amon et al., 2022a) showed that, assuming the
fiducial lensing cosmology from (Heymans et al., 2021), the data–to-theory
ratio $\Delta\Sigma_{\mathrm{data}}/\Delta\Sigma_{\mathrm{theory}}$ is close
to one. This result, again, outlines a fundamental incompatibility between the
Planck18 cosmology and the shape of $P(k)$.

##  5 Conclusion

In this study, we used weak gravitational lensing data from KiDS-1000 and
DES-Y3 to constrain the cosmological parameters and the HMF. The $(p,q)$
parameters add two new degrees of freedom to the calculation of the matter
power spectrum $P(k)$. We examined various HMF configurations, where the
cosmological parameters were either fixed to specific values (KiDS, DES, or
Planck cosmologies) or allowed to vary. To maintain the maximum consistency of
our analysis with the original KiDS-1000 and DES-Y3 studies, we used their
original pipelines and only made minimal changes to the samplers and modelling
choices, as outlined in section 1.

We modified HMcode (Mead et al., 2015, 2016) to account for these changes. We
first assessed how different parameter choices impact the Sheth-Tormen HMF
(Sheth & Tormen, 1999), the matter power spectrum $P(k)$ and the shear
correlation function $\xi_{\pm}$. We found that the amplitude ${S_{8}}$ is
generally higher when both the cosmology and $(p,q)$ are free. The new
${S_{8}}$ posterior is in better agreement with the Planck18 value, and we do
recover the ST1 $(p,q)$ when allowing for both cosmological and HMF parameters
to vary, but at the expense of reducing consistency along the
$(\sigma_{8},\Omega_{m})$ degeneracy. Moreover, $P(k)$ must be reduced by 25%
at scale $k\sim 1~{}{\rm h/Mpc}$, in agreement with other lensing studies. We
should be cautious about the interpretation of this result: this is not a
measurement of the mass function, nor a proof that it should be modified: if
we found enough evidence that p and q are not equal to ST1 values, then that
would indicate a missing component in the modelling, but since the modelling
is tuned to the simulations, it also indicates a missing factor in the
simulations, which may or may not be related to HMF. This is what is happening
for the Planck18 cosmology, where the $(p,q)$ contours strongly exclude the
ST1 values. A useful analogy is the shear calibration parameter $m$, which
should be zero if no residual systematics are present in the shear
measurement, but all lensing surveys generally find a small but non-zero
value, indicating a bias in the raw shear measurement.

It is interesting that DES and KiDS collaboration et al. (2023) found that
their hybrid pipeline reconciles DES-y3, KiDS-1000 and Planck18 by a unified
set of choices for the intrinsic alignment, baryonic feedback, the choice of
angular scales and halo model code. It would be very interesting to see if,
with this setup, the $(p,q)$ parameters would converge to $(0.3,0.707)$, and
this is left for future work.

The ${S_{8}}$ tension has triggered a lot of excitement about new physics
because the residual systematics coming from shape measurement and photometric
redshift calibration seem to be well under control. Unfortunately, we have to
admit that the understanding of the theory (and specifically the calculation
of $P(k)$ and $C_{\ell}$ with all the complexities) is now lagging behind the
constraining power of existing lensing surveys. This will be even worse for
stage IV surveys, therefore there is an urge to make progress in this area.
Fortunately, the potential of weak lensing science is still vastly under-
utilized on data: there is a lot of constraining power on intrinsic alignment,
baryonic physics and the non-linear growth of structures which can come from
combining high order statistics, peak statistics, cross-correlations, and what
we could call "global statistics" such as void- lensing and density split
statistics. Almost no effort has been done so far to use these tools together
to constrain the complexities of $P(k)$, this is the next frontier.

## Acknowledgements

We thank Gary Hinshaw and Robert Reischke for fruitful discussions. SG and LVW
acknowledge support by the University of British Columbia, Canada’s NSERC. TT
acknowledges funding from the Leverhulme Trust and the Swiss National Science
Foundation under the Ambizione project PZ00P2_193352.

Author contributions: All authors have contributed to the development and
writing of this paper. The authorship list is divided into three groups: the
lead authors (SG, MAD, & LW) followed by two alphabetical groups. The first
alphabetical group (MA, AM, TT) includes key contributors to both the
scientific analysis and the data products. The second group (ZY) includes
those who have either made significant contributions to the data products or
the scientific analysis.

## Data Availability

In this paper the data is analysed with these open-source python packages:
astropy (Astropy Collaboration et al., 2013, 2018), Camb (Lewis & Challinor,
2011), Cosmosis (Zuntz et al., 2015), Getdist (Lewis, 2019), KCAP (Tröster et
al., 2021), matplotlib (Hunter, 2007), Multinest (Feroz et al., 2009), numpy
(Harris et al., 2020), pyccl (Chisari et al., 2019), pyHMcode (Tröster et al.,
2022), and scipy (Virtanen et al., 2020). The data analysis pipeline and
generated posterior chain files used in this article will be shared on
reasonable request to the corresponding author (SG).

## References

  * Abbott et al. (2020) Abbott T. M. C., et al., 2020, Phys. Rev. D, 102, 023509
  * Aihara et al. (2018) Aihara H., et al., 2018, PASJ, 70, S4
  * Aiola et al. (2020) Aiola S., et al., 2020, arXiv e-prints, p. arXiv:2007.07288
  * Amon & Efstathiou (2022) Amon A., Efstathiou G., 2022, arXiv e-prints, p. arXiv:2206.11794
  * Amon et al. (2022a) Amon A., et al., 2022a, arXiv e-prints, p. arXiv:2202.07440
  * Amon et al. (2022b) Amon A., et al., 2022b, Phys. Rev. D, 105, 023514
  * Asgari et al. (2021) Asgari M., et al., 2021, A&A, 645, A104
  * Asgari et al. (2023) Asgari M., Mead A. J., Heymans C., 2023, arXiv e-prints, p. arXiv:2303.08752
  * Astropy Collaboration et al. (2013) Astropy Collaboration et al., 2013, A&A, 558, A33
  * Astropy Collaboration et al. (2018) Astropy Collaboration et al., 2018, AJ, 156, 123
  * Bagla et al. (2009) Bagla J. S., Khandai N., Kulkarni G., 2009, arXiv e-prints, p. arXiv:0908.2702
  * Bahcall & Cen (1993) Bahcall N. A., Cen R., 1993, ApJ, 407, L49
  * Barreira et al. (2014) Barreira A., Li B., Hellwing W. A., Lombriser L., Baugh C. M., Pascoli S., 2014, J. Cosmology Astropart. Phys., 2014, 029
  * Baugh et al. (2019) Baugh C. M., et al., 2019, MNRAS, 483, 4922
  * Beltz-Mohrmann & Berlind (2021) Beltz-Mohrmann G. D., Berlind A. A., 2021, arXiv e-prints, p. arXiv:2103.05076
  * Benítez (2000) Benítez N., 2000, ApJ, 536, 571
  * Bird et al. (2012) Bird S., Viel M., Haehnelt M. G., 2012, MNRAS, 420, 2551
  * Blandford et al. (1991) Blandford R. D., Saust A. B., Brainerd T. G., Villumsen J. V., 1991, MNRAS, 251, 600
  * Blazek et al. (2019) Blazek J. A., MacCrann N., Troxel M. A., Fang X., 2019, Phys. Rev. D, 100, 103506
  * Bocquet et al. (2016) Bocquet S., Saro A., Dolag K., Mohr J. J., 2016, MNRAS, 456, 2361
  * Böhringer et al. (2017) Böhringer H., Chon G., Fukugita M., 2017, A&A, 608, A65
  * Bullock et al. (2001) Bullock J. S., Kolatt T. S., Sigad Y., Somerville R. S., Kravtsov A. V., Klypin A. A., Primack J. R., Dekel A., 2001, MNRAS, 321, 559
  * Castro et al. (2016) Castro T., Marra V., Quartin M., 2016, MNRAS, 463, 1666
  * Castro et al. (2021) Castro T., Borgani S., Dolag K., Marra V., Quartin M., Saro A., Sefusatti E., 2021, MNRAS, 500, 2316
  * Chisari et al. (2018) Chisari N. E., et al., 2018, MNRAS, 480, 3962
  * Chisari et al. (2019) Chisari N. E., et al., 2019, ApJS, 242, 2
  * Costanzi et al. (2013) Costanzi M., Villaescusa-Navarro F., Viel M., Xia J.-Q., Borgani S., Castorina E., Sefusatti E., 2013, J. Cosmology Astropart. Phys., 2013, 012
  * Costanzi et al. (2019) Costanzi M., et al., 2019, MNRAS, 488, 4779
  * Crocce & Scoccimarro (2006) Crocce M., Scoccimarro R., 2006, Phys. Rev. D, 73, 063520
  * Crocce et al. (2015) Crocce M., Castander F. J., Gaztañaga E., Fosalba P., Carretero J., 2015, MNRAS, 453, 1513
  * DES Collaboration et al. (2022) DES Collaboration et al., 2022, Phys. Rev. D, 105, 023520
  * DES and KiDS collaboration et al. (2023) DES and KiDS collaboration D., et al., 2023, arXiv e-prints, p. arXiv:2305.17173
  * Despali et al. (2016) Despali G., Giocoli C., Angulo R. E., Tormen G., Sheth R. K., Baso G., Moscardini L., 2016, MNRAS, 456, 2486
  * Diemer (2021) Diemer B., 2021, ApJ, 909, 112
  * Drlica-Wagner et al. (2018) Drlica-Wagner A., et al., 2018, ApJS, 235, 33
  * Edge et al. (2013) Edge A., Sutherland W., Kuijken K., Driver S., McMahon R., Eales S., Emerson J. P., 2013, The Messenger, 154, 32
  * Fang et al. (2020) Fang X., Eifler T., Krause E., 2020, MNRAS, 497, 2699
  * Feroz et al. (2009) Feroz F., Hobson M. P., Bridges M., 2009, MNRAS, 398, 1601
  * Flaugher et al. (2015) Flaugher B., et al., 2015, AJ, 150, 150
  * Fosalba et al. (2015a) Fosalba P., Gaztañaga E., Castander F. J., Crocce M., 2015a, MNRAS, 447, 1319
  * Fosalba et al. (2015b) Fosalba P., Crocce M., Gaztañaga E., Castander F. J., 2015b, MNRAS, 448, 2987
  * Gatti et al. (2022) Gatti M., et al., 2022, MNRAS, 510, 1223
  * Geach (2012) Geach J. E., 2012, MNRAS, 419, 2633
  * Giblin et al. (2021) Giblin B., et al., 2021, A&A, 645, A105
  * Gunn (1967) Gunn J. E., 1967, ApJ, 150, 737
  * Hamilton (2000) Hamilton A. J. S., 2000, MNRAS, 312, 257
  * Harnois-Déraps & van Waerbeke (2015) Harnois-Déraps J., van Waerbeke L., 2015, MNRAS, 450, 2857
  * Harnois-Déraps et al. (2018) Harnois-Déraps J., et al., 2018, MNRAS, 481, 1337
  * Harris et al. (2020) Harris C. R., et al., 2020, Nature, 585, 357
  * Hartley et al. (2022) Hartley W. G., et al., 2022, MNRAS, 509, 3547
  * Heymans et al. (2012) Heymans C., et al., 2012, MNRAS, 427, 146
  * Heymans et al. (2021) Heymans C., et al., 2021, A&A, 646, A140
  * Hildebrandt et al. (2017) Hildebrandt H., et al., 2017, MNRAS, 465, 1454
  * Hildebrandt et al. (2020) Hildebrandt H., et al., 2020, A&A, 633, A69
  * Hinshaw et al. (2013) Hinshaw G., et al., 2013, ApJS, 208, 19
  * Hoekstra et al. (2011) Hoekstra H., Donahue M., Conselice C. J., McNamara B. R., Voit G. M., 2011, ApJ, 726, 48
  * Hubble (1936) Hubble E. P., 1936, Realm of the Nebulae
  * Hunter (2007) Hunter J. D., 2007, Computing in Science & Engineering, 9, 90
  * Jain & Taylor (2003) Jain B., Taylor A., 2003, Phys. Rev. Lett., 91, 141302
  * Jenkins et al. (2001) Jenkins A., Frenk C. S., White S. D. M., Colberg J. M., Cole S., Evrard A. E., Couchman H. M. P., Yoshida N., 2001, MNRAS, 321, 372
  * Joachimi et al. (2020) Joachimi B., et al., 2020, arXiv e-prints, p. arXiv:2007.01844
  * Knebe et al. (2011) Knebe A., et al., 2011, MNRAS, 415, 2293
  * Kohonen (2001) Kohonen T., 2001, Self-Organizing Maps. Springer series in information sciences, 2001, xx, 501
  * Krause & Eifler (2017) Krause E., Eifler T., 2017, MNRAS, 470, 2100
  * Kuijken (2011) Kuijken K., 2011, The Messenger, 146, 8
  * Kuijken et al. (2015) Kuijken K., et al., 2015, MNRAS, 454, 3500
  * Kuijken et al. (2019) Kuijken K., et al., 2019, A&A, 625, A2
  * Kulkarni & Ostriker (2022) Kulkarni M., Ostriker J. P., 2022, MNRAS, 510, 1425
  * Leauthaud et al. (2017) Leauthaud A., et al., 2017, MNRAS, 467, 3024
  * Lemos et al. (2022) Lemos P., et al., 2022, MNRAS,
  * Lewis (2019) Lewis A., 2019, arXiv e-prints, p. arXiv:1910.13970
  * Lewis & Challinor (2011) Lewis A., Challinor A., 2011, CAMB: Code for Anisotropies in the Microwave Background, Astrophysics Source Code Library, record ascl:1102.026 (ascl:1102.026)
  * Li et al. (2019) Li P., Lelli F., McGaugh S., Pawlowski M. S., Zwaan M. A., Schombert J., 2019, ApJ, 886, L11
  * Limber (1953) Limber D. N., 1953, ApJ, 117, 134
  * Lovell (2020) Lovell M. R., 2020, MNRAS, 493, L11
  * Marsh (2015) Marsh D. J. E., 2015, Phys. Rev. D, 91, 123520
  * Mead et al. (2015) Mead A. J., Peacock J. A., Heymans C., Joudaki S., Heavens A. F., 2015, MNRAS, 454, 1958
  * Mead et al. (2016) Mead A. J., Heymans C., Lombriser L., Peacock J. A., Steele O. I., Winther H. A., 2016, MNRAS, 459, 1468
  * Mead et al. (2020) Mead A. J., Tröster T., Heymans C., Van Waerbeke L., McCarthy I. G., 2020, A&A, 641, A130
  * Mead et al. (2021) Mead A. J., Brieden S., Tröster T., Heymans C., 2021, MNRAS, 502, 1401
  * Myles et al. (2021) Myles J., et al., 2021, MNRAS, 505, 4249
  * Newman (2008) Newman J. A., 2008, ApJ, 684, 88
  * Ondaro-Mallea et al. (2022) Ondaro-Mallea L., Angulo R. E., Zennaro M., Contreras S., Aricò G., 2022, MNRAS, 509, 6077
  * Peacock & Smith (2000) Peacock J. A., Smith R. E., 2000, MNRAS, 318, 1144
  * Planck Collaboration et al. (2020) Planck Collaboration et al., 2020, A&A, 641, A6
  * Press & Schechter (1974) Press W. H., Schechter P., 1974, ApJ, 187, 425
  * Raichoor et al. (2014) Raichoor A., et al., 2014, ApJ, 797, 102
  * Sánchez et al. (2022) Sánchez C., et al., 2022, Phys. Rev. D, 105, 083529
  * Schaye et al. (2023) Schaye J., et al., 2023, arXiv e-prints, p. arXiv:2306.04024
  * Schneider & Teyssier (2015) Schneider A., Teyssier R., 2015, J. Cosmology Astropart. Phys., 2015, 049
  * Schneider et al. (2020) Schneider A., Stoira N., Refregier A., Weiss A. J., Knabenhans M., Stadel J., Teyssier R., 2020, J. Cosmology Astropart. Phys., 2020, 019
  * Secco et al. (2022) Secco L. F., et al., 2022, Phys. Rev. D, 105, 023515
  * Sevilla-Noarbe et al. (2021) Sevilla-Noarbe I., et al., 2021, ApJS, 254, 24
  * Sheth & Tormen (1999) Sheth R. K., Tormen G., 1999, MNRAS, 308, 119
  * Sheth et al. (2001) Sheth R. K., Mo H. J., Tormen G., 2001, MNRAS, 323, 1
  * Smith et al. (2003) Smith R. E., et al., 2003, MNRAS, 341, 1311
  * Smith et al. (2007) Smith R. E., Scoccimarro R., Sheth R. K., 2007, Phys. Rev. D, 75, 063512
  * Takahashi et al. (2012) Takahashi R., Sato M., Nishimichi T., Taruya A., Oguri M., 2012, ApJ, 761, 152
  * Tormen (1998) Tormen G., 1998, MNRAS, 297, 648
  * Tröster et al. (2021) Tröster T., et al., 2021, A&A, 649, A88
  * Tröster et al. (2022) Tröster T., et al., 2022, A&A, 660, A27
  * Virtanen et al. (2020) Virtanen P., et al., 2020, Nature Methods, 17, 261
  * White & Frenk (1991) White S. D. M., Frenk C. S., 1991, ApJ, 379, 52
  * Wright et al. (2020) Wright A. H., Hildebrandt H., van den Busch J. L., Heymans C., Joachimi B., Kannawadi A., Kuijken K., 2020, A&A, 640, L14
  * Zuntz et al. (2015) Zuntz J., et al., 2015, Astronomy and Computing, 12, 45
  * de Jong et al. (2017) de Jong J. T. A., et al., 2017, A&A, 604, A134
  * eBOSS Collaboration et al. (2020) eBOSS Collaboration et al., 2020, arXiv e-prints, p. arXiv:2007.08991
  * van Daalen et al. (2020) van Daalen M. P., McCarthy I. G., Schaye J., 2020, MNRAS, 491, 2424
  * van den Busch et al. (2022) van den Busch J. L., et al., 2022, A&A, 664, A170

##  Appendix A Sampling and Scale Cuts

###  A.1 Impact of Sampling

Figure 9: Posterior contours ($68\%$ and $95\%$) in the
${\Omega_{\mathrm{m}}}-{\sigma_{8}}$, ${\Omega_{\mathrm{m}}}-{S_{8}}$ and
${\sigma_{8}}-{S_{8}}$ planes for KiDS-1000, DES-y3, and Planck 2018\. The
filled salmon pink contours (KiDS-1000 Fiducial) are the $\xi_{\pm}$ analysis
using the fiducial scale cut from Asgari et al. (2021). The dashed rosybrown
contours (KiDS-1000 Full) use all data points from KiDS-1000 $\xi_{\pm}$. The
solid indigo blue and the dashed bright azure contours are the $\Lambda$CDM-
optimised DES-Y3 cosmic shear results with Halofit, using polychord (fiducial)
and multinest respectively. The filled cyan contours are the HMcode-based
DES-y3 results under KiDS-1000 prior of $B$ parameters. The slate grey
contours are the TTTEEE+low$\ell$+low E+lensing+BAO result from Planck
Collaboration et al. (2020) as the comparison to lensing results.

Here, we show how a particular choice of sampling may affect the posterior
distribution in cosmic shear calculations. Figure 9 shows different
combinations of sampling (multinest versus polychord) and non-linear power
spectrum calculations (Halofit versus HMcode). The original DES-y3 setup,
$\Lambda$CDM-optimised cosmic shear results, which uses polychord with
Halofit, is shown by the indigo contours. It shows a $2.8\sigma$ tension in
${S_{8}}$ compared to the Planck Collaboration et al. (2020) cosmology. The
multinest sampling is shown by the cyan dashed contours, and finally, the
multinest using HMcode is shown by the cyan solid contours. The latter has
only a $1.6\sigma$ tension with the Planck value. When we use HMcode, we
always set the halo concentration parameter $B$ from Bullock et al. (2001) as
a free parameter ($B\in[2,3.13]$) – like the set-up of the KiDS-1000 fiducial
analysis from Asgari et al. (2021).

###  A.2 Impact of scale cuts

Figure 10: Posterior contours ($68\%$ and $95\%$) in the
${\Omega_{\mathrm{m}}}-{\sigma_{8}}-{S_{8}}-B$ parameter space for KiDS-1000
given the Standard Sheth-Tormen HMF for different scale cuts. The rosybrown
solid contour is for no scale cut ($\theta\in(0.5^{\prime},300^{\prime})$).
The olive solid line, lime dashed line and orchid dotted line contours are for
the scale cuts bins $\theta\in(0.5^{\prime},4.2^{\prime})$,
$\theta\in(4.2^{\prime},35.5^{\prime})$ and
$\theta\in(35.5^{\prime},300^{\prime})$, respectively.

Different choices of scale cuts may also impact the posterior distributions by
altering the sensitivity of the statistical estimator to various physical
effects. For instance, it is known that small scales are more sensitive to
baryonic physics. The pink solid contours and the brown dashed contours in
Figure 9 respectively show the fiducial KiDS-1000 analysis (also appeared as
ST-1 cosmological analysis in our main text) and the analysis with no scale
cuts at all ($\theta\in[0.5^{\prime},300^{\prime}]$ for both $\xi_{\pm}$). The
KiDS-1000 fiducial analysis applied a $\theta>4^{\prime}$ cut. The difference
between the two is marginal.

Figure 10 shows the result of various scale cuts on KiDS-1000 data, assuming
the HMF parameters are fixed to ST1. We use three bins of different scale-cut,
separated by the geometric progressed values of $4.2^{\prime}$ and
$35.5^{\prime}$. We can see that the different scale bins do not show a
particular trend for a low or high ${S_{8}}$, which means that different
choices of scale cuts do not impact the cosmological inference significantly,
at least for this data set and survey area.