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

Pierre A Burger, Lucas Porth, Sven Heydenreich, Laila Linke, Niek Wielders, Peter Schneider, Marika Asgari, Tiago Castro, Klaus Dolag, Joachim Harnois-Deraps, et al.

► To cite this version:

Pierre A Burger, Lucas Porth, Sven Heydenreich, Laila Linke, Niek Wielders, et al. KiDS-1000 cosmology: Combined second- and third-order shear statistics. Astron.Astrophys., 2024, 683, pp.A103. 10.1051/0004-6361/202347986 . hal-04230048

HAL Id: hal-04230048

https://hal.science/hal-04230048

Submitted on 19 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_1.png" /></figure>

HAL
open science

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

Pierre A Burger, Lucas Porth, Sven Heydenreich, Laila Linke, Niek Wielders, Peter Schneider, Marika Asgari, Tiago Castro, Klaus Dolag, Joachim Harnois-Deraps, et al.

► To cite this version:

Pierre A Burger, Lucas Porth, Sven Heydenreich, Laila Linke, Niek Wielders, et al. KiDS-1000 cosmology: Combined second- and third-order shear statistics. Astronomy & Astrophysics - A&A, 2024, 683, pp.A103. 10.1051/0004-6361/202347986. hal-04230048

HAL Id: hal-04230048

https://hal.science/hal-04230048v1

Submitted on 19 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# KiDS-1000 cosmology: Combined second- and third-order shear statistics

Pierre A. Burger<fnref n="1" />, Lucas Porth<fnref n="1" />, Sven Heydenreich<fnref n="1" />,<fnref n="3" />, Laila Linke<fnref n="1" />,<fnref n="2" />, Niek Wielders<fnref n="1" />, Peter Schneider<fnref n="1" />, Marika Asgari<fnref n="4" />, Tiago Castro<fnref n="5" />,<fnref n="6" />,<fnref n="7" />, Klaus Dolag<fnref n="8" />,<fnref n="9" />, Joachim Harnois-Déraps<fnref n="10" />, Hendrik Hildebrandt<fnref n="11" />, Konrad Kuijken<fnref n="12" />, and Nicolas Martinet<fnref n="13" />

Received 15 September 2023 / Accepted 6 December 2023

## ABSTRACT

**Aims.** In this work, we perform the first cosmological parameter analysis of the fourth release of Kilo Degree Survey (KiDS-1000) data with second- and third-order shear statistics. This paper builds on a series of studies aimed at describing the roadmap to third-order shear statistics.

**Methods.** We derived and tested a combined model of the second-order shear statistic, namely, the COSEBIs and the third-order aperture mass statistics ($\langle M_{ap}^2 \rangle$) in a tomographic set-up. We validated our pipeline with N-body mock simulations of the KiDS-1000 data release. To model the second- and third-order statistics, we used the latest version of HMCODE2020 for the power spectrum and BiHALOFIT for the bispectrum. Furthermore, we used an analytic description to model intrinsic alignments and hydro-dynamical simulations to model the effect of baryonic feedback processes. Lastly, we decreased the dimension of the data vector significantly by considering only equal smoothing radii for the $\langle M_{ap}^2 \rangle$ part of the data vector. This makes it possible to carry out a data analysis of the KiDS-1000 data release using a combined analysis of COSEBIs and third-order shear statistics.

**Results.** We first validated the accuracy of our modelling by analysing a noise-free mock data vector, assuming the KiDS-1000 error budget, finding a shift in the maximum of the posterior distribution of the matter density parameter, $\Delta\Omega_m < 0.02 \sigma_{\Omega_m}$, and of the structure growth parameter, $\Delta S_8 < 0.05 \sigma_{S_8}$. Lastly, we performed the first KiDS-1000 cosmological analysis using a combined analysis of second- and third-order shear statistics, where we constrained $\Omega_m = 0.248_{-0.055}^{+0.062}$ and $S_8 = \sigma_8 \sqrt{\Omega_m/0.3} = 0.772 \pm 0.022$. The geometric average on the errors of $\Omega_m$ and $S_8$ of the combined statistics decreases, compared to the second-order statistic, by a factor of 2.2.

**Key words.** gravitation – gravitational lensing: weak – methods: analytical – methods: numerical – cosmological parameters – large-scale structure of Universe

## 1. Introduction

Gravitational lensing describes the deflection of light by massive objects. It is sensitive to baryonic and dark matter and, therefore, ideal for probing the total matter distribution in the Universe. Since the distribution of matter is highly sensitive to cosmological parameters, it is excellent to test and probe the standard model of cosmology, called the Λ cold dark matter model (ΛCDM). Although the ΛCDM model can describe observations of the early Universe, such as the cosmic microwave background (CMB; e.g., Planck Collaboration VI 2020), or the Local Universe, such as the observed large-scale structure (LSS) of mat-

ter and galaxies (Sánchez et al. 2017), with remarkable accuracy, it is being put under stress due to tension observed between early and local probes. A ~2σ tension is in the structure growth parameter $S_8 = \sigma_8 \sqrt{\Omega_m/0.3}$, where $\sigma_8$ is the normalisation of the power spectrum and $\Omega_m$ is the total matter density parameter (Hildebrandt et al. 2017; Planck Collaboration VI 2020; Joudaki et al. 2020; Heymans et al. 2021; DES Collaboration 2022; Di Valentino et al. 2021; Dalal et al. 2023), suggests that the Local Universe is less clustered than what is expected from early-time measurements when extrapolated under the ΛCDM model. If these tensions are not due to systematics, extensions to the ΛCDM model are necessary, which will be tested with

<footnote n="1">¹ University of Bonn, Argelander-Institut für Astronomie, Auf dem Hügel 71, 53121 Bonn, Germany
e-mail: pburger@astro.uni-bonn.de</footnote>

<footnote n="2">² Universität Innsbruck, Institut für Astro- und Teilchenphysik, Technikerstr. 25/8, 6020 Innsbruck, Austria</footnote>

<footnote n="3">³ Department of Astronomy and Astrophysics, University of California, Santa Cruz, 1156 High Street, Santa Cruz, CA 95064, USA</footnote>

<footnote n="4">⁴ E.A. Milne Centre, University of Hull, Cottingham Road, Hull HU6 7RX, UK</footnote>

<footnote n="5">⁵ INAF – Osservatorio Astronomico di Trieste, via Tiepolo 11, 34131 Trieste, Italy</footnote>

<footnote n="6">⁶ INFN – Sezione di Trieste, 34100 Trieste, Italy</footnote>

<footnote n="7">⁷ IFPU – Institute for Fundamental Physics of the Universe, via Beirut 2, 34151 Trieste, Italy</footnote>

<footnote n="8">⁸ Universitäts-Sternwarte, Fakultät für Physik, Ludwig-Maximilians-Universität München, Scheinerstr. 1, 81679 München, Germany</footnote>

<footnote n="9">⁹ Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Straße 1, 85741 Garching, Germany</footnote>

<footnote n="10">¹⁰ School of Mathematics, Statistics and Physics, Newcastle University, Newcastle upon Tyne, NE1 7RU, UK</footnote>

<footnote n="11">¹¹ Ruhr University Bochum, Faculty of Physics and Astronomy, Astronomical Institute (AIRUB), German Centre for Cosmological
Lensing, 44780 Bochum, Germany</footnote>

<footnote n="12">¹² Leiden Observatory, Leiden University, PO Box 9513, 2300 RA Leiden, The Netherlands</footnote>

<footnote n="13">¹³ Aix-Marseille Université, CNRS, CNES, LAM, Marseille, France</footnote>
# KiDS-1000 cosmology: Combined second- and third-order shear statistics

Pierre A. Burger<fnref n="1" />, Lucas Porth<fnref n="1" />, Sven Heydenreich<fnref n="1" />,<fnref n="3" />, Laila Linke<fnref n="1" />,<fnref n="2" />, Niek Wielders<fnref n="1" />, Peter Schneider<fnref n="1" />, Marika Asgari<fnref n="4" />, Tiago Castro<fnref n="5" />,<fnref n="6" />,<fnref n="7" />, Klaus Dolag<fnref n="8" />,<fnref n="9" />, Joachim Harnois-Déraps<fnref n="10" />, Hendrik Hildebrandt<fnref n="11" />, Konrad Kuijken<fnref n="12" />, and Nicolas Martinet<fnref n="13" />

Received 15 September 2023 / Accepted 6 December 2023

## ABSTRACT

**Aims.** In this work, we perform the first cosmological parameter analysis of the fourth release of Kilo Degree Survey (KiDS-1000) data with second- and third-order shear statistics. This paper builds on a series of studies aimed at describing the roadmap to third-order shear statistics.

**Methods.** We derived and tested a combined model of the second-order shear statistic, namely, the COSEBIs and the third-order aperture mass statistics ($\langle M_{ap}^2 \rangle$) in a tomographic set-up. We validated our pipeline with N-body mock simulations of the KiDS-1000 data release. To model the second- and third-order statistics, we used the latest version of HMCODE2020 for the power spectrum and BiHALOFIT for the bispectrum. Furthermore, we used an analytic description to model intrinsic alignments and hydro-dynamical simulations to model the effect of baryonic feedback processes. Lastly, we decreased the dimension of the data vector significantly by considering only equal smoothing radii for the $\langle M_{ap}^2 \rangle$ part of the data vector. This makes it possible to carry out a data analysis of the KiDS-1000 data release using a combined analysis of COSEBIs and third-order shear statistics.

**Results.** We first validated the accuracy of our modelling by analysing a noise-free mock data vector, assuming the KiDS-1000 error budget, finding a shift in the maximum of the posterior distribution of the matter density parameter, $\Delta\Omega_m < 0.02 \sigma_{\Omega_m}$, and of the structure growth parameter, $\Delta S_8 < 0.05 \sigma_{S_8}$. Lastly, we performed the first KiDS-1000 cosmological analysis using a combined analysis of second- and third-order shear statistics, where we constrained $\Omega_m = 0.248_{-0.055}^{+0.062}$ and $S_8 = \sigma_8 \sqrt{\Omega_m/0.3} = 0.772 \pm 0.022$. The geometric average on the errors of $\Omega_m$ and $S_8$ of the combined statistics decreases, compared to the second-order statistic, by a factor of 2.2.

**Key words.** gravitation – gravitational lensing: weak – methods: analytical – methods: numerical – cosmological parameters – large-scale structure of Universe

## 1. Introduction

Gravitational lensing describes the deflection of light by massive objects. It is sensitive to baryonic and dark matter and, therefore, ideal for probing the total matter distribution in the Universe. Since the distribution of matter is highly sensitive to cosmological parameters, it is excellent to test and probe the standard model of cosmology, called the Λ cold dark matter model (ΛCDM). Although the ΛCDM model can describe observations of the early Universe, such as the cosmic microwave background (CMB; e.g., Planck Collaboration VI 2020), or the Local Universe, such as the observed large-scale structure (LSS) of mat-

ter and galaxies (Sánchez et al. 2017), with remarkable accuracy, it is being put under stress due to tension observed between early and local probes. A ~2σ tension is in the structure growth parameter $S_8 = \sigma_8 \sqrt{\Omega_m/0.3}$, where $\sigma_8$ is the normalisation of the power spectrum and $\Omega_m$ is the total matter density parameter (Hildebrandt et al. 2017; Planck Collaboration VI 2020; Joudaki et al. 2020; Heymans et al. 2021; DES Collaboration 2022; Di Valentino et al. 2021; Dalal et al. 2023), suggests that the Local Universe is less clustered than what is expected from early-time measurements when extrapolated under the ΛCDM model. If these tensions are not due to systematics, extensions to the ΛCDM model are necessary, which will be tested with

<footnote n="1">¹ University of Bonn, Argelander-Institut für Astronomie, Auf dem Hügel 71, 53121 Bonn, Germany
e-mail: pburger@astro.uni-bonn.de</footnote>

<footnote n="2">² Universität Innsbruck, Institut für Astro- und Teilchenphysik, Technikerstr. 25/8, 6020 Innsbruck, Austria</footnote>

<footnote n="3">³ Department of Astronomy and Astrophysics, University of California, Santa Cruz, 1156 High Street, Santa Cruz, CA 95064, USA</footnote>

<footnote n="4">⁴ E.A. Milne Centre, University of Hull, Cottingham Road, Hull HU6 7RX, UK</footnote>

<footnote n="5">⁵ INAF – Osservatorio Astronomico di Trieste, via Tiepolo 11, 34131 Trieste, Italy</footnote>

<footnote n="6">⁶ INFN – Sezione di Trieste, 34100 Trieste, Italy</footnote>

<footnote n="7">⁷ IFPU – Institute for Fundamental Physics of the Universe, via Beirut 2, 34151 Trieste, Italy</footnote>

<footnote n="8">⁸ Universitäts-Sternwarte, Fakultät für Physik, Ludwig-Maximilians-Universität München, Scheinerstr. 1, 81679 München, Germany</footnote>

<footnote n="9">⁹ Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Straße 1, 85741 Garching, Germany</footnote>

<footnote n="10">¹⁰ School of Mathematics, Statistics and Physics, Newcastle University, Newcastle upon Tyne, NE1 7RU, UK</footnote>

<footnote n="11">¹¹ Ruhr University Bochum, Faculty of Physics and Astronomy, Astronomical Institute (AIRUB), German Centre for Cosmological
Lensing, 44780 Bochum, Germany</footnote>

<footnote n="12">¹² Leiden Observatory, Leiden University, PO Box 9513, 2300 RA Leiden, The Netherlands</footnote>

<footnote n="13">¹³ Aix-Marseille Université, CNRS, CNES, LAM, Marseille, France</footnote>
the next generation of cosmic shear surveys such as *Euclid* (Laureijs et al. 2011) or the *Vera Rubin Observatory Legacy Survey* of Space and Time (LSST, Ivezic et al. 2008).

Commonly, two-point statistics for weak lensing and galaxy positions are used to infer cosmological parameters since they can be modelled accurately and systematic inaccuracies are well understood (Schneider et al. 1998; Troxel et al. 2018; Hildebrandt et al. 2017, 2020; Hikage et al. 2019; Asgari et al. 2019, 2020). Two-point statistics are excellent for capturing the entire information content of a Gaussian random field. However, non-linear gravitational instabilities created a significant amount of non-Gaussian features during the evolution of the Universe, such that the local matter distribution departed strongly from a Gaussian field. Therefore, higher order statistics are needed to extract all the available information in the local LSS of matter and galaxies. Furthermore, such higher order statistics usually depend differently on cosmological parameters and systematic effects such as intrinsic alignment (IA), meaning that a joint investigation of second- and higher order statistics tightens cosmological parameter constraints (see, e.g., Kilbinger & Schneider 2005; Bergé et al. 2010; Pires et al. 2012; Fu et al. 2014; Pyne & Joachimi 2021). Recently used examples of higher order statistics are the peak count statistics (Martinet et al. 2018; Harnois-Déraps et al. 2021), persistent homology (Heydenreich et al. 2021), density split statistics (Gruen et al. 2018; Burger et al. 2023), and the integrated three-point correlation function used in (Halder et al. 2021; Halder & Barreira 2022), along with a second- and third-order convergence moment analysis (Gatti et al. 2022).

This work considers second- and third-order shear statistics, where the former probes the variance and the latter the skewness of the LSS at various scales. Our chosen second-order statistic is the $E_n$-modes of the COSEBIs (Schneider et al. 2010), and the third-order statistic $\langle M_{\text{ap}}^3 \rangle$ is described in Schneider et al. (2005) and recently measured in the Dark Energy Survey Year (DES) 3 Results (Secco et al. 2022). Furthermore, this work belongs to a series of papers that aim for cosmological parameter analyses using third-order shear statistics. In the first paper, Heydenreich et al. (2023, hereafter H23), we validated the analytical fitting formulae for a non-tomographic analysis and the conversion from three-point correlation function (3PCF) of cosmic shear to $\langle M_{\text{ap}}^3 \rangle$. We found that $\langle M_{\text{ap}}^3 \rangle$, even though they are combined from the shear 3PCF at different scales, contain a similar amount of information on $\Omega_m$ and $S_8$ as the 3PCF itself. The fact that $\langle M_{\text{ap}}^3 \rangle$ can be measured from 3PCF is very convenient since this allows unbiased estimates for any survey geometry. A fast computation method of the aperture mass statistics by measuring the shear 3PCF is tackled in Porth et al. (2023). Lastly, adding third-order statistics increases the dimension of the data vector significantly, and an analytical expression for the covariance is preferred, which is derived and validated for $\langle M_{\text{ap}}^3 \rangle$ in a non-tomographic setup in Linke et al. (2023). However, as this analysis considers combining second-order statistics with $\langle M_{\text{ap}}^3 \rangle$ in a tomographic setup, and a joint covariance matrix has not been derived yet (Wielders et al., in prep.), we must still rely on numerical simulations to determine the covariance matrix.

This article presents the first cosmological parameter analysis using the fourth data release of KiDS (KiDS-1000), combining second- and third-order shear statistics. We show the cosmological results, preceded by validation of several extensions of the analytical model to allow for a tomographic analysis and to include astrophysical effects such as IA (e.g., Joachimi et al. 2015) or baryonic feedback processes

(Chisari et al. 2015). We aim to find the smallest set of scales that retain most of the cosmological information.

The paper is structured as follows: In Sect. 2, we review the basics of the second- and third-order shear statistics, extend the modelling to a tomographic analysis and describe our method to model the IA analytically. In Sect. 4, we describe the KiDS-1000 data and in Sect. 5, we introduce the simulation data used to validate our model and to correct it for baryonic feedback processes. In Sect. 6, we briefly review our method to perform cosmological parameter interference. In Sects. 7 and 8, we validate our analysis pipeline against several systematics. The final cosmological results are presented in Sect. 9, and we present our conclusions in Sect. 10.

## 2. Theoretical background

This section offers a review of the basics of weak gravitational lensing formalism and aperture statistics. For more detailed reviews, we refer to Bartelmann & Schneider (2001), Hoekstra & Jain (2008), Munshi et al. (2008), Bartelmann (2010), and Kilbinger (2015). In this work, we assume a spatially flat universe, such that the comoving angular-diameter distance is expressed as $f_K[\chi(z)] = \chi(z)$, where $\chi(z)$ is the comoving distance at redshift $z$. Given the matter density, $\rho(x, z)$, at comoving position, $x$, and redshift, $z$, the density contrast is $\delta(x, z) = \frac{\rho(x, z) - \rho(z)}{\bar{\rho}(z)} - 1$, where $\bar{\rho}(z)$ is the average matter density at redshift, $z$. The dimensionless surface mass density or convergence, $\kappa$, for sources at redshift, $z$, is determined by the line-of-sight integration as

$$ \kappa(\boldsymbol{\theta}, z) = \frac{3\Omega_m H_0^2}{2c^2} \int_0^{\chi(z)} d\chi' \frac{\chi'[\chi(z) - \chi']}{\chi(z)} \frac{\delta(\chi'\boldsymbol{\theta}, z)}{a(\chi')}, \quad (1) $$

where $\boldsymbol{\theta}$ is the angular position on the sky, $H_0$ is the Hubble constant, and $a$ is the scale factor. The second argument of $\delta$ simultaneously describes the radial direction and the cosmological epoch, related through the light-cone condition $|c dt| = a(z) d\chi$.

### 2.1. Limber projections of power- and bispectrum

Given the Fourier transform $\hat{\delta}$ of the matter density contrast, the matter power spectrum, $P_{\delta\delta}(k, z)$, and bispectrum, $B_{\delta\delta\delta}(k_1, k_2, k_3, z)$, are:

$$ \langle \hat{\delta}(\mathbf{k}_1, z) \hat{\delta}(\mathbf{k}_2, z) \rangle = (2\pi)^3 \delta_D(\mathbf{k}_1 + \mathbf{k}_2) P_{\delta\delta}(k_1, z), \quad (2) $$

$$ \langle \hat{\delta}(\mathbf{k}_1, z) \hat{\delta}(\mathbf{k}_2, z) \hat{\delta}(\mathbf{k}_3, z) \rangle = (2\pi)^3 \delta_D(\mathbf{k}_1 + \mathbf{k}_2 + \mathbf{k}_3) \times B_{\delta\delta\delta}(k_1, k_2, k_3, z), \quad (3) $$

where $\delta_D$ is the Dirac-delta distribution. The statistical isotropy of the Universe implies that the power- and bispectrum only depend on the moduli of the $k$-vectors.

The projected power- and bispectrum can then be computed using the Limber approximation (Limber 1954; Kaiser & Jaffe 1997; Bernardeau et al. 1997; Schneider et al. 1998; LoVerde & Afshordi 2008),

$$ P_{\kappa\kappa}^{(ij)}(\ell) = \int_{0}^{\chi_{\text{max}}} d\chi \frac{g^{(i)}(\chi) g^{(j)}(\chi)}{a^2(\chi)} P_{\delta\delta} \left[ \frac{\ell + 1/2}{\chi}, z(\chi) \right], \quad (4) $$

$$ B_{\kappa\kappa\kappa}^{(ijk)}(\ell_1, \ell_2, \ell_3) = \int_0^{\chi_{\text{max}}} d\chi \frac{g^{(i)}(\chi) g^{(j)}(\chi) g^{(k)}(\chi)}{a^3(\chi)\chi} \times B_{\delta\delta\delta} \left[ \frac{\ell_1}{\chi}, \frac{\ell_2}{\chi}, \frac{\ell_3}{\chi}, z(\chi) \right], \quad (5) $$
the next generation of cosmic shear surveys such as *Euclid* (Laureijs et al. 2011) or the *Vera Rubin* Observatory Legacy Survey of Space and Time (LSST, Ivezic et al. 2008).

Commonly, two-point statistics for weak lensing and galaxy positions are used to infer cosmological parameters since they can be modelled accurately and systematic inaccuracies are well understood (Schneider et al. 1998; Troxel et al. 2018; Hildebrandt et al. 2017, 2020; Hikage et al. 2019; Asgari et al. 2019, 2020). Two-point statistics are excellent for capturing the entire information content of a Gaussian random field. However, non-linear gravitational instabilities created a significant amount of non-Gaussian features during the evolution of the Universe, such that the local matter distribution departed strongly from a Gaussian field. Therefore, higher order statistics are needed to extract all the available information in the local LSS of matter and galaxies. Furthermore, such higher order statistics usually depend differently on cosmological parameters and systematic effects such as intrinsic alignment (IA), meaning that a joint investigation of second- and higher order statistics tightens cosmological parameter constraints (see, e.g., Kilbinger & Schneider 2005; Bergé et al. 2010; Pires et al. 2012; Fu et al. 2014; Pyne & Joachimi 2021). Recently used examples of higher order statistics are the peak count statistics (Martinet et al. 2018; Harnois-Déraps et al. 2021), persistent homology (Heydenreich et al. 2021), density split statistics (Gruen et al. 2018; Burger et al. 2023), and the integrated three-point correlation function used in (Halder et al. 2021; Halder & Barreira 2022), along with a second- and third-order convergence moment analysis (Gatti et al. 2022).

This work considers second- and third-order shear statistics, where the former probes the variance and the latter the skewness of the LSS at various scales. Our chosen second-order statistic is the $E_n$-modes of the COSEBIs (Schneider et al. 2010), and the third-order statistic $\langle M_{\text{ap}}^3 \rangle$ is described in Schneider et al. (2005) and recently measured in the Dark Energy Survey Year (DES) 3 Results (Secco et al. 2022). Furthermore, this work belongs to a series of papers that aim for cosmological parameter analyses using third-order shear statistics. In the first paper, Heydenreich et al. (2023, hereafter H23), we validated the analytical fitting formulae for a non-tomographic analysis and the conversion from three-point correlation function (3PCF) of cosmic shear to $\langle M_{\text{ap}}^3 \rangle$. We found that $\langle M_{\text{ap}}^3 \rangle$, even though they are combined from the shear 3PCF at different scales, contain a similar amount of information on $\Omega_m$ and $S_8$ as the 3PCF itself. The fact that $\langle M_{\text{ap}}^3 \rangle$ can be measured from 3PCF is very convenient since this allows unbiased estimates for any survey geometry. A fast computation method of the aperture mass statistics by measuring the shear 3PCF is tackled in Porth et al. (2023). Lastly, adding third-order statistics increases the dimension of the data vector significantly, and an analytical expression for the covariance is preferred, which is derived and validated for $\langle M_{\text{ap}}^3 \rangle$ in a non-tomographic setup in Linke et al. (2023). However, as this analysis considers combining second-order statistics with $\langle M_{\text{ap}}^3 \rangle$ in a tomographic setup, and a joint covariance matrix has not been derived yet (Wielders et al., in prep.), we must still rely on numerical simulations to determine the covariance matrix.

This article presents the first cosmological parameter analysis using the fourth data release of KiDS (KiDS-1000), combining second- and third-order shear statistics. We show the cosmological results, preceded by validation of several extensions of the analytical model to allow for a tomographic analysis and to include astrophysical effects such as IA (e.g., Joachimi et al. 2015) or baryonic feedback processes

(Chisari et al. 2015). We aim to find the smallest set of scales that retain most of the cosmological information.

The paper is structured as follows: In Sect. 2, we review the basics of the second- and third-order shear statistics, extend the modelling to a tomographic analysis and describe our method to model the IA analytically. In Sect. 4, we describe the KiDS-1000 data and in Sect. 5, we introduce the simulation data used to validate our model and to correct it for baryonic feedback processes. In Sect. 6, we briefly review our method to perform cosmological parameter interference. In Sects. 7 and 8, we validate our analysis pipeline against several systematics. The final cosmological results are presented in Sect. 9, and we present our conclusions in Sect. 10.

## 2. Theoretical background

This section offers a review of the basics of weak gravitational lensing formalism and aperture statistics. For more detailed reviews, we refer to Bartelmann & Schneider (2001), Hoekstra & Jain (2008), Munshi et al. (2008), Bartelmann (2010), and Kilbinger (2015). In this work, we assume a spatially flat universe, such that the comoving angular-diameter distance is expressed as $f_K[\chi(z)] = \chi(z)$, where $\chi(z)$ is the comoving distance at redshift $z$. Given the matter density, $\rho(x, z)$, at comoving position, $x$, and redshift, $z$, the density contrast is $\delta(x, z) = \frac{\rho(x, z) - \rho(z)}{\bar{\rho}(z)} - 1$, where $\bar{\rho}(z)$ is the average matter density at redshift, $z$. The dimensionless surface mass density or convergence, $\kappa$, for sources at redshift, $z$, is determined by the line-of-sight integration as

$$ \kappa(\boldsymbol{\theta}, z) = \frac{3\Omega_m H_0^2}{2c^2} \int_0^{\chi(z)} d\chi' \frac{\chi'[\chi(z) - \chi']}{\chi(z)} \frac{\delta(\chi'\boldsymbol{\theta}, z)}{a(\chi')}, \quad (1) $$

where $\boldsymbol{\theta}$ is the angular position on the sky, $H_0$ is the Hubble constant, and $a$ is the scale factor. The second argument of $\delta$ simultaneously describes the radial direction and the cosmological epoch, related through the light-cone condition $|c dt| = a(z) d\chi$.

### 2.1. Limber projections of power- and bispectrum

Given the Fourier transform $\hat{\delta}$ of the matter density contrast, the matter power spectrum, $P_{\delta\delta}(k, z)$, and bispectrum, $B_{\delta\delta\delta}(k_1, k_2, k_3, z)$, are:

$$ \langle \hat{\delta}(\mathbf{k}_1, z) \hat{\delta}(\mathbf{k}_2, z) \rangle = (2\pi)^3 \delta_D(\mathbf{k}_1 + \mathbf{k}_2) P_{\delta\delta}(k_1, z), \quad (2) $$

$$ \langle \hat{\delta}(\mathbf{k}_1, z) \hat{\delta}(\mathbf{k}_2, z) \hat{\delta}(\mathbf{k}_3, z) \rangle = (2\pi)^3 \delta_D(\mathbf{k}_1 + \mathbf{k}_2 + \mathbf{k}_3) \times B_{\delta\delta\delta}(k_1, k_2, k_3, z), \quad (3) $$

where $\delta_D$ is the Dirac-delta distribution. The statistical isotropy of the Universe implies that the power- and bispectrum only depend on the moduli of the $k$-vectors.

The projected power- and bispectrum can then be computed using the Limber approximation (Limber 1954; Kaiser & Jaffe 1997; Bernardeau et al. 1997; Schneider et al. 1998; LoVerde & Afshordi 2008),

$$ P_{\kappa\kappa}^{(ij)}(\ell) = \int_{0}^{\chi_{\text{max}}} d\chi \frac{g^{(i)}(\chi) g^{(j)}(\chi)}{a^2(\chi)} P_{\delta\delta} \left[ \frac{\ell + 1/2}{\chi}, z(\chi) \right], \quad (4) $$

$$ B_{\kappa\kappa\kappa}^{(ijk)}(\ell_1, \ell_2, \ell_3) = \int_0^{\chi_{\text{max}}} d\chi \frac{g^{(i)}(\chi) g^{(j)}(\chi) g^{(k)}(\chi)}{a^3(\chi)\chi} \times B_{\delta\delta\delta} \left[ \frac{\ell_1}{\chi}, \frac{\ell_2}{\chi}, \frac{\ell_3}{\chi}, z(\chi) \right], \quad (5) $$
where $\ell_3 = |\ell_1 + \ell_2|$, and $g^{(i)}(\chi)$ denotes the lensing efficiency and is defined as

$$g^{(i)}(\chi) = \frac{3\Omega_m H_0^2}{2c^2} \int_{\chi}^{\chi_{\text{max}}} d\chi' n^{(i)}(\chi') \frac{\chi' - \chi}{\chi'}, \quad (6)$$

with $n^{(i)}(\chi)$ being the redshift probability distribution of the ith tomographic zph-bin. Since the Limber approximation breaks down for small values of $\ell$ (Kilbinger et al. 2017), we consider only scales below 4° to model $\langle M_{\text{ap}}^3 \rangle$. In principle, we could have also used the modified Limber approximation and shift $\ell_{1,2,3} \rightarrow \ell_{1,2,3} + 1/2$, but since we found a difference at maximum for the largest filter radii of 1.5%, we neglected it here. To model the non-linear matter power spectrum $P_{\delta\delta}$, we use the revised HMCODE2020 model of Mead et al. (2021) and for the matter bispectrum $B_{\delta\delta\delta}$ the BIHALOFIT of Takahashi et al. (2020).

## 2.2. Non-linear alignment model

The impact of galaxy IA is a known contaminating signal to the cosmic shear measurements that must be accounted for in all weak lensing studies. To model the effects of intrinsic alignment, we use the non-linear alignment (NLA) model (Bridle & King 2007), which is a one-parameter model described as

$$f_{\text{IA}} = -\frac{A_{\text{IA}} \bar{C}_1 \bar{\rho}(z)}{D_+(z)}, \quad (7)$$

where $D_+$ is the linear growth factor at redshift $z$, $\bar{C}_1 = 5 \times 10^{-14} M_{\odot}^{-1} h^{-2} \text{Mpc}^3$, as calibrated in Brown et al. (2002), and $A_{\text{IA}}$ captures the coupling strength between the matter density and the tidal field.

By considering the tidal alignment field, $\delta_I$, as a biased tracer of the matter density contrast field $\delta$, neglecting all higher order bias terms, we get:

$$\delta_I(\mathbf{x}, z) = f_{\text{IA}}\delta(\mathbf{x}, z). \quad (8)$$

With this, we find that:

$$P_{\delta\delta_1}(k, z) = f_{\text{IA}} P_{\delta\delta}(k, z), \quad (9)$$

$$P_{\delta_1\delta_1}(k, z) = f_{\text{IA}}^2 P_{\delta\delta}(k, z), \quad (10)$$

where $P_{\delta\delta}(k, z)$ is the non-linear matter power spectrum. The projected power spectra then follow to

$$P_{\text{GI}}^{(ij)}(\ell) = \int_{0}^{\infty} d\chi \frac{g^{(i)}(\chi) n^{(j)}(\chi)}{a(\chi)\chi} P_{\delta\delta_1}(\ell/\chi, \chi), \quad (11)$$

$$P_{\text{II}}^{(ij)}(\ell) = \int_{0}^{\infty} d\chi \frac{n^{(i)}(\chi') n^{(j)}(\chi)}{\chi^2} P_{\delta_1\delta_1}(\ell/\chi, \chi), \quad (12)$$

and the total projected power spectrum becomes

$$P^{(ij)}(\ell) = P_{\text{GG}}^{(ij)}(\ell) + P_{\text{GI}}^{(ij)}(\ell) + P_{\text{GI}}^{(ji)}(\ell) + P_{\text{II}}^{(ij)}(\ell). \quad (13)$$

The term $P_{\text{GG}}^{(ij)}(\ell) = P_{\text{KK}}^{(ij)}(\ell)$ describes the actual lensing signal, which is given in Eq. (4) with the weighting kernel in Eq. (6). The II contribution describes how two galaxies spatially close together tend to be aligned. The term $\text{GI}(z_i > z_j)$ describes the fact that high matter density regions align the lower redshift galaxies but also affect the shear of the background galaxies. While the II term is dominant if galaxies of the same tomographic bin are considered, GI and IG start to dominate if galaxies of separated tomographic bins are considered. With $z_i < z_j$,

the term GI is expected to vanish for two tomographic bins with no significant redshift overlap.

Following the ansatz for modelling IA in the power spectrum, we get

$$B_{\delta\delta\delta_1}(k_1, k_2, k_3, z) = f_{\text{IA}} B_{\delta\delta\delta}(k_1, k_2, k_3, z), \quad (14)$$

$$B_{\delta\delta\delta_1}(k_1, k_2, k_3, z) = f_{\text{IA}}^2 B_{\delta\delta\delta}(k_1, k_2, k_3, z), \quad (15)$$

$$B_{\delta_1\delta_1\delta_1}(k_1, k_2, k_3, z) = f_{\text{IA}}^3 B_{\delta\delta\delta}(k_1, k_2, k_3, z), \quad (16)$$

where $B_{\delta\delta\delta}$ is the non-linear matter bispectrum, which we calculated with BIHALOFIT. The projected bispectra are

$$B_{\text{GGI}}^{(ijk)}(\ell_1, \ell_2, \ell_3) = \int_{0}^{\infty} d\chi \frac{g^{(i)}(\chi) g^{(j)}(\chi) n^{(k)}(\chi)}{a^2(\chi)\chi^2} B_{\delta\delta\delta_1}(\ell_1/\chi, \ell_2/\chi, \ell_3/\chi, \chi), \quad (17)$$

$$B_{\text{GII}}^{(ijk)}(\ell_1, \ell_2, \ell_3) = \int_{0}^{\infty} d\chi \frac{g^{(i)}(\chi) n^{(j)}(\chi) n^{(k)}(\chi)}{a(\chi)\chi^3} B_{\delta\delta_1\delta_1}(\ell_1/\chi, \ell_2/\chi, \ell_3/\chi, \chi), \quad (18)$$

$$B_{\text{III}}^{(ijk)}(\ell_1, \ell_2, \ell_3) = \int_{0}^{\infty} d\chi \frac{n^{(i)}(\chi) n^{(j)}(\chi) n^{(k)}(\chi)}{\chi^4} B_{\delta_1\delta_1\delta_1}(\ell_1/\chi, \ell_2/\chi, \ell_3/\chi, \chi). \quad (19)$$

The total projected bispectrum is:

$$B^{(ijk)}(\hat{l}) = B_{\text{GGG}}^{(ijk)}(\hat{l}) + B_{\text{GGI}}^{(ijk)}(\hat{l}) + B_{\text{GGI}}^{(kij)}(\hat{l}) + B_{\text{GGI}}^{(jki)}(\hat{l}) + B_{\text{GII}}^{(ijk)}(\hat{l}) + B_{\text{GII}}^{(kij)}(\hat{l}) + B_{\text{GII}}^{(jki)}(\hat{l}) + B_{\text{III}}^{(ijk)}(\hat{l}), \quad (20)$$

where the tuple is $\hat{l} = (\ell_1, \ell_2, \ell_3)$. The interpretation of all these terms is analogous to the ones from the power spectrum, where $B_{\text{GGG}}^{(ijk)}(\hat{l}) = B_{\text{KKK}}^{(ijk)}(\hat{l})$ is the actual pure lensing signal given by Eq. (5) with the weighting kernel in Eq. (6).

## 2.3. Aperture mass statistics

One of the major problems of weak lensing mass reconstruction techniques is the mass-sheet degeneracy (hereafter MSD; Falco et al. 1985; Schneider & Seitz 1995), which corresponds to adding a uniform surface mass density without affecting lensing observables such as the shear. However, it is possible to define quantities invariant under the MSD, one example being the aperture mass statistics (Schneider 1996; Bartelmann & Schneider 2001). Another advantage of aperture mass statistics is that they separate the signal into so-called E- and B-modes (Schneider et al. 2002), where, to leading order, the weak gravitational lensing effect cannot create B-modes. Lastly, as shown in H23, aperture statistics are an excellent strategy to compress thousands of bins of 3PCF into a few hundred $\langle M_{\text{ap}}^3 \rangle$ bins.

The aperture mass $M_{\text{ap}}$ at position $\mathbf{v}$ with filter radius $\theta_{\text{ap}}$ is defined through the convergence $\kappa$, as follows:

$$M_{\text{ap}}(\mathbf{v}; \theta_{\text{ap}}) = \int d^2 \mathbf{v}' U_{\theta_{\text{ap}}}(|\mathbf{v}'|) \kappa(\mathbf{v} + \mathbf{v}',) \quad (21)$$

where $U_{\theta_{\text{ap}}}(\mathbf{v}')$ is a compensated filter such that $\int d\mathbf{v}' \mathbf{v}' U_{\theta_{\text{ap}}}(\mathbf{v}') = 0$. The tangential shear component $\gamma_t$ of the complex shear in Cartesian coordinates $\gamma = \gamma_1 + i\gamma_2$, is defined as $\gamma_t = -\Re(e^{2i\phi})$, where $\phi$ is the polar angle of $\mathbf{v}'$. Given the tangential shear $\gamma_t$, the aperture mass $M_{\text{ap}}$ can also be calculated as

$$M_{\text{ap}}(\mathbf{v}; \theta_{\text{ap}}) = \int d^2 \mathbf{v}' Q_{\theta_{\text{ap}}}(|\mathbf{v}'|) \gamma_t(\mathbf{v} + \mathbf{v}',) \quad (22)$$
where $\ell_3 = |\ell_1 + \ell_2|$, and $g^{(i)}(\chi)$ denotes the lensing efficiency and is defined as

$$g^{(i)}(\chi) = \frac{3\Omega_m H_0^2}{2c^2} \int_{\chi}^{\chi_{\text{max}}} d\chi' n^{(i)}(\chi') \frac{\chi' - \chi}{\chi'}, \quad (6)$$

with $n^{(i)}(\chi)$ being the redshift probability distribution of the ith tomographic zph-bin. Since the Limber approximation breaks down for small values of $\ell$ (Kilbinger et al. 2017), we consider only scales below 4° to model $\langle M_{\text{ap}}^3 \rangle$. In principle, we could have also used the modified Limber approximation and shift $\ell_{1,2,3} \rightarrow \ell_{1,2,3} + 1/2$, but since we found a difference at maximum for the largest filter radii of 1.5%, we neglected it here. To model the non-linear matter power spectrum $P_{\delta\delta}$, we use the revised HMCODE2020 model of Mead et al. (2021) and for the matter bispectrum $B_{\delta\delta\delta}$ the BIHALOFIT of Takahashi et al. (2020).

## 2.2. Non-linear alignment model

The impact of galaxy IA is a known contaminating signal to the cosmic shear measurements that must be accounted for in all weak lensing studies. To model the effects of intrinsic alignment, we use the non-linear alignment (NLA) model (Bridle & King 2007), which is a one-parameter model described as

$$f_{\text{IA}} = -\frac{A_{\text{IA}} \bar{C}_1 \bar{\rho}(z)}{D_+(z)}, \quad (7)$$

where $D_+$ is the linear growth factor at redshift $z$, $\bar{C}_1 = 5 \times 10^{-14} M_{\odot}^{-1} h^{-2} \text{Mpc}^3$, as calibrated in Brown et al. (2002), and $A_{\text{IA}}$ captures the coupling strength between the matter density and the tidal field.

By considering the tidal alignment field, $\delta_I$, as a biased tracer of the matter density contrast field $\delta$, neglecting all higher order bias terms, we get:

$$\delta_I(\mathbf{x}, z) = f_{\text{IA}} \delta(\mathbf{x}, z). \quad (8)$$

With this, we find that:

$$P_{\delta\delta_1}(k, z) = f_{\text{IA}} P_{\delta\delta}(k, z), \quad (9)$$

$$P_{\delta_1 \delta_1}(k, z) = f_{\text{IA}}^2 P_{\delta\delta}(k, z), \quad (10)$$

where $P_{\delta\delta}(k, z)$ is the non-linear matter power spectrum. The projected power spectra then follow to

$$P_{\text{GI}}^{(ij)}(\ell) = \int_{0}^{\infty} d\chi \frac{g^{(i)}(\chi) n^{(j)}(\chi)}{a(\chi)\chi} P_{\delta\delta_1}(\ell/\chi, \chi), \quad (11)$$

$$P_{\text{II}}^{(ij)}(\ell) = \int_{0}^{\infty} d\chi \frac{n^{(i)}(\chi') n^{(j)}(\chi)}{\chi^2} P_{\delta_1 \delta_1}(\ell/\chi, \chi), \quad (12)$$

and the total projected power spectrum becomes

$$P^{(ij)}(\ell) = P_{\text{GG}}^{(ij)}(\ell) + P_{\text{GI}}^{(ij)}(\ell) + P_{\text{GI}}^{(ji)}(\ell) + P_{\text{II}}^{(ij)}(\ell). \quad (13)$$

The term $P_{\text{GG}}^{(ij)}(\ell) = P_{\text{KK}}^{(ij)}(\ell)$ describes the actual lensing signal, which is given in Eq. (4) with the weighting kernel in Eq. (6). The II contribution describes how two galaxies spatially close together tend to be aligned. The term $\text{GI}(z_i > z_j)$ describes the fact that high matter density regions align the lower redshift galaxies but also affect the shear of the background galaxies. While the II term is dominant if galaxies of the same tomographic bin are considered, GI and IG start to dominate if galaxies of separated tomographic bins are considered. With $z_i < z_j$,

the term GI is expected to vanish for two tomographic bins with no significant redshift overlap.

Following the ansatz for modelling IA in the power spectrum, we get

$$B_{\delta\delta\delta_1}(k_1, k_2, k_3, z) = f_{\text{IA}} B_{\delta\delta\delta}(k_1, k_2, k_3, z), \quad (14)$$

$$B_{\delta\delta_1\delta_1}(k_1, k_2, k_3, z) = f_{\text{IA}}^2 B_{\delta\delta\delta}(k_1, k_2, k_3, z), \quad (15)$$

$$B_{\delta_1\delta_1\delta_1}(k_1, k_2, k_3, z) = f_{\text{IA}}^3 B_{\delta\delta\delta}(k_1, k_2, k_3, z), \quad (16)$$

where $B_{\delta\delta\delta}$ is the non-linear matter bispectrum, which we calculated with BIHALOFIT. The projected bispectra are

$$B_{\text{GGI}}^{(ijk)}(\ell_1, \ell_2, \ell_3) = \int_{0}^{\infty} d\chi \frac{g^{(i)}(\chi) g^{(j)}(\chi) n^{(k)}(\chi)}{a^2(\chi)\chi^2} B_{\delta\delta\delta_1}(\ell_1/\chi, \ell_2/\chi, \ell_3/\chi, \chi), \quad (17)$$

$$B_{\text{GII}}^{(ijk)}(\ell_1, \ell_2, \ell_3) = \int_{0}^{\infty} d\chi \frac{g^{(i)}(\chi) n^{(j)}(\chi) n^{(k)}(\chi)}{a(\chi)\chi^3} B_{\delta\delta_1\delta_1}(\ell_1/\chi, \ell_2/\chi, \ell_3/\chi, \chi), \quad (18)$$

$$B_{\text{III}}^{(ijk)}(\ell_1, \ell_2, \ell_3) = \int_{0}^{\infty} d\chi \frac{n^{(i)}(\chi) n^{(j)}(\chi) n^{(k)}(\chi)}{\chi^4} B_{\delta_1\delta_1\delta_1}(\ell_1/\chi, \ell_2/\chi, \ell_3/\chi, \chi). \quad (19)$$

The total projected bispectrum is:

$$B^{(ijk)}(\hat{l}) = B_{\text{GGG}}^{(ijk)}(\hat{l}) + B_{\text{GGI}}^{(ijk)}(\hat{l}) + B_{\text{GGI}}^{(kij)}(\hat{l}) + B_{\text{GGI}}^{(jki)}(\hat{l}) \\ + B_{\text{GII}}^{(ijk)}(\hat{l}) + B_{\text{GII}}^{(kij)}(\hat{l}) + B_{\text{GII}}^{(jki)}(\hat{l}) + B_{\text{III}}^{(ijk)}(\hat{l}), \quad (20)$$

where the tuple is $\hat{l} = (\ell_1, \ell_2, \ell_3)$. The interpretation of all these terms is analogous to the ones from the power spectrum, where $B_{\text{GGG}}^{(ijk)}(\hat{l}) = B_{\text{KKK}}^{(ijk)}(\hat{l})$ is the actual pure lensing signal given by Eq. (5) with the weighting kernel in Eq. (6).

## 2.3. Aperture mass statistics

One of the major problems of weak lensing mass reconstruction techniques is the mass-sheet degeneracy (hereafter MSD; Falco et al. 1985; Schneider & Seitz 1995), which corresponds to adding a uniform surface mass density without affecting lensing observables such as the shear. However, it is possible to define quantities invariant under the MSD, one example being the aperture mass statistics (Schneider 1996; Bartelmann & Schneider 2001). Another advantage of aperture mass statistics is that they separate the signal into so-called E- and B-modes (Schneider et al. 2002), where, to leading order, the weak gravitational lensing effect cannot create B-modes. Lastly, as shown in H23, aperture statistics are an excellent strategy to compress thousands of bins of 3PCF into a few hundred $\langle M_{\text{ap}}^3 \rangle$ bins.

The aperture mass $M_{\text{ap}}$ at position $\mathbf{v}$ with filter radius $\theta_{\text{ap}}$ is defined through the convergence $\kappa$, as follows:

$$M_{\text{ap}}(\mathbf{v}; \theta_{\text{ap}}) = \int d^2 \mathbf{v}' U_{\theta_{\text{ap}}}(|\mathbf{v}'|) \kappa(\mathbf{v} + \mathbf{v}',) \quad (21)$$

where $U_{\theta_{\text{ap}}}(\mathbf{v}')$ is a compensated filter such that $\int d\mathbf{v}' \mathbf{v}' U_{\theta_{\text{ap}}}(\mathbf{v}') = 0$. The tangential shear component $\gamma_t$ of the complex shear in Cartesian coordinates $\gamma = \gamma_1 + i\gamma_2$, is defined as $\gamma_t = -\Re(e^{2i\phi})$, where $\phi$ is the polar angle of $\mathbf{v}'$. Given the tangential shear $\gamma_t$, the aperture mass $M_{\text{ap}}$ can also be calculated as

$$M_{\text{ap}}(\mathbf{v}; \theta_{\text{ap}}) = \int d^2 \mathbf{v}' Q_{\theta_{\text{ap}}}(|\mathbf{v}'|) \gamma_t(\mathbf{v} + \mathbf{v}',) \quad (22)$$
where $Q_{\theta_{ap}}$ is related to $U_{\theta_{ap}}$ via

$$Q_{\theta_{ap}}(\vartheta) = \frac{2}{\vartheta^2} \int_{0}^{\vartheta} d\vartheta' \vartheta' U_{\theta_{ap}}(\vartheta') - U_{\theta_{ap}}(\vartheta). \qquad (23)$$

We define $U_{\theta_{ap}}(\vartheta) = \theta_{ap}^{-2}u(\vartheta/\theta_{ap})$, denote by $\hat{u}(\eta)$ the Fourier transform of $u$ and use the filter function introduced in Crittenden et al. (2002),

$$\begin{aligned} u(x) &= \frac{1}{2\pi} \left(1 - \frac{x^2}{2}\right) e^{-x^2/2}, \quad \hat{u}(\eta) = \frac{\eta^2}{2} e^{-\eta^2/2}, \\ Q_{\theta_{ap}}(\vartheta) &= \frac{\vartheta^2}{4\pi\theta_{ap}^4} \exp\left(-\frac{\vartheta^2}{2\theta_{ap}^2}\right). \end{aligned} \qquad (24)$$

## 2.4. Modelling aperture mass moments

The expectation value of the aperture mass $\langle M_{ap} \rangle(\theta_{ap})$, which approximates the ensemble average over all positions $\vartheta$, vanishes by construction. However, the second-order (variance) of the aperture mass is nonzero and can be calculated as

$$\langle M_{ap}^2 \rangle^{(ij)}(\theta_{ap}) = \int \frac{d\ell \, \ell}{2\pi} P_{\kappa\kappa}^{(ij)}(\ell) \hat{u}^2(\theta_{ap}\ell). \qquad (25)$$

Equivalently, the third-order moment of the aperture statistics $\langle M_{ap}^3 \rangle$, can be computed from the convergence bispectrum via (Jarvis et al. 2004; Schneider et al. 2005)

$$\langle M_{ap}^3 \rangle^{(ijk)}(\theta_{ap,1}, \theta_{ap,2}, \theta_{ap,3}) = \int \frac{d^2 \ell_1}{(2\pi)^2} \int \frac{d^2 \ell_2}{(2\pi)^2} \times B_{\kappa\kappa}^{(ijk)}(\ell_1, \ell_2, \ell_3) \hat{u}(\theta_{ap,1} \ell_1) \hat{u}(\theta_{ap,2} \ell_2) \hat{u}(\theta_{ap,3} \ell_3), \qquad (26)$$

where $\ell_3 = |\ell_1 + \ell_2|$. Later, we differentiate between equal filter radii, for which $\theta_{ap,1} = \theta_{ap,2} = \theta_{ap,3}$ and non-equal filter radii, where the values of $\theta_{ap,i}$ are all allowed to vary.

## 2.5. Modelling COSEBIs

Schneider et al. (2010) introduced the complete orthogonal sets of E/B-integrals (COSEBIs), which are defined via the two-point shear statistics $\xi_{\pm}(\vartheta) = \langle \gamma_t \gamma_t \rangle(\vartheta) \pm \langle \gamma_{\times} \gamma_{\times} \rangle(\vartheta)$ on a finite angular range

$$E_n^{(ij)} = \frac{1}{2} \int_{\vartheta_{\min}}^{\vartheta_{\max}} d\vartheta \, \vartheta [T_{+n}(\vartheta)\xi_{+}^{(ij)}(\vartheta) + T_{-n}(\vartheta)\xi_{-}^{(ij)}(\vartheta)], \qquad (27)$$

$$B_n^{(ij)} = \frac{1}{2} \int_{\vartheta_{\min}}^{\vartheta_{\max}} d\vartheta \, \vartheta [T_{+n}(\vartheta)\xi_{+}^{(ij)}(\vartheta) - T_{-n}(\vartheta)\xi_{-}^{(ij)}(\vartheta)], \qquad (28)$$

where $T_{\pm n}(\vartheta)$ are filter functions with support in $[\vartheta_{\min}, \vartheta_{\max}]$ (Schneider et al. 2010). If

$$\begin{aligned} \int_{\vartheta_{\min}}^{\vartheta_{\max}} d\vartheta \, \vartheta T_{+n}(\vartheta) J_0(\ell \vartheta) &= \int_{\vartheta_{\min}}^{\vartheta_{\max}} d\vartheta \, \vartheta T_{-n}(\vartheta) J_4(\ell \vartheta) \\ &:= W_n(\ell), \end{aligned} \qquad (29)$$

where $J_{0,4}$ are the zeroth and fourth order Bessel function, then the COSEBIs offer the advantage of cleanly separating all well-defined E- and B-modes within the range $[\vartheta_{\min}, \vartheta_{\max}]$. This is not given, for instance, for second-order aperture mass statistics, as this would require information of $\xi_{\pm}$ over the full space. Analytically, the COSEBIs can be calculated from the E-mode power

spectrum $P_{\kappa\kappa}^{(ij)}(\ell)$ defined in Eq. (4) and a B-mode angular power spectra $P_{\kappa\kappa,B}^{(ij)}(\ell)$ as:

$$E_n^{(ij)} = \int_0^\infty \frac{d\ell \, \ell}{2\pi} P_{\kappa\kappa}^{(ij)}(\ell) W_n(\ell), \qquad (30)$$

$$B_n^{(ij)} = \int_0^\infty \frac{d\ell \, \ell}{2\pi} P_{\kappa\kappa,B}^{(ij)}(\ell) W_n(\ell). \qquad (31)$$

Since B-modes cannot be created by gravitational lensing directly, we neglected the modelling of $B_n$ for this work. Therefore, we focus, henceforth, only on the $E_n$ modes from the COSEBIs.

# 3. Measuring shear statistics

As discussed in H23 $\langle M_{ap}^3 \rangle$ can be estimated from data in three ways, which we review below. The first uses the convergence field, $\kappa$, the second uses the shear field, and the third uses correlation functions. The latter is used to measure the real data and compute the covariance matrix used for the real data analysis.

## 3.1. Measuring shear statistics from convergence maps

If the convergence field, $\kappa$, is available (e.g., from simulations), the easiest way to measure $M_{ap}$ is to use Eq. (21). If the convergence field does not contain masks, this estimator is also unbiased as long as a border of the size of the filter function is removed from the aperture mass field before calculating the spatial average. However, no boundaries need to be removed for unmasked full-sky convergence fields to get an unbiased estimator. To estimate the aperture mass from (full-sky) convergence maps, the maps are smoothed with the HEALPY<fnref n="1" /> (Zonca et al. 2019) function SMOOTHING. This function needs a beam window function created by the function BEAM2BL, which is determined from the corresponding $U_{\theta_{ap}}$ filter.

In the absence of B-modes the $E_n$ modes can also be calculated by using convolutions as (Schneider et al. 2010):

$$E_n^{(ij)} = \langle \kappa^{(i)}(\vartheta) \kappa_s^{(j)}(\vartheta; n) \rangle, \qquad (32)$$

where

$$\kappa_s^{(i)}(\vartheta; n) = \frac{1}{2\pi} \int d^2\theta'_{ap} T_{+n}(|\vartheta - \vartheta'|) \kappa^{(i)}(\theta'_{ap}) \qquad (33)$$

with $T_{+n}$ introduced in Eq. (29).

## 3.2. Measuring shear statistics from galaxy catalogues

In a real survey, the convergence is not observable and can be inferred only from the measured shear. However, estimating $M_{ap}$ from these reconstructed convergence maps is not accurate as the reconstructed convergence is necessarily smoothed and potentially also contains other systematic effects caused by masks or boundaries of the survey (Seitz & Schneider 1997).

However, as motivated by Eq. (22), the aperture mass can also be estimated from an unmasked shear field (e.g., for simulations). Given a galaxy catalogue, where galaxies are only found at specific positions such that the number of galaxies within an aperture varies on the sky, Eq. (22) needs to get modified to

$$\tilde{M}_{ap}^{(i)}(\vartheta; \theta_{ap}) = \left[ \sum_{i} Q_{\theta_{ap}}(|\vartheta - \vartheta_i|) \right]^{-1} \sum_{j} Q_{\theta_{ap}}(|\vartheta - \vartheta_j|) \varepsilon_{t,j}^{(i)}, \qquad (34)$$

<footnote n="1">¹ http://healpix.sourceforge.net</footnote>
where $Q_{\theta_{ap}}$ is related to $U_{\theta_{ap}}$ via

$$Q_{\theta_{ap}}(\vartheta) = \frac{2}{\vartheta^2} \int_0^\vartheta d\vartheta' \vartheta' U_{\theta_{ap}}(\vartheta') - U_{\theta_{ap}}(\vartheta). \quad (23)$$

We define $U_{\theta_{ap}}(\vartheta) = \theta_{ap}^{-2}u(\vartheta/\theta_{ap})$, denote by $\hat{u}(\eta)$ the Fourier transform of $u$ and use the filter function introduced in Crittenden et al. (2002),

$$\begin{align}
u(x) &= \frac{1}{2\pi} \left(1 - \frac{x^2}{2}\right) e^{-x^2/2}, \quad \hat{u}(\eta) = \frac{\eta^2}{2} e^{-\eta^2/2}, \nonumber \\
Q_{\theta_{ap}}(\vartheta) &= \frac{\vartheta^2}{4\pi\theta_{ap}^4} \exp\left(-\frac{\vartheta^2}{2\theta_{ap}^2}\right). \tag{24}
\end{align}$$

## 2.4. Modelling aperture mass moments

The expectation value of the aperture mass $\langle M_{ap} \rangle(\theta_{ap})$, which approximates the ensemble average over all positions $\vartheta$, vanishes by construction. However, the second-order (variance) of the aperture mass is nonzero and can be calculated as

$$\langle M_{ap}^2 \rangle^{(ij)}(\theta_{ap}) = \int \frac{d\ell \, \ell}{2\pi} P_{\kappa\kappa}^{(ij)}(\ell) \hat{u}^2(\theta_{ap}\ell). \quad (25)$$

Equivalently, the third-order moment of the aperture statistics $\langle M_{ap}^3 \rangle$, can be computed from the convergence bispectrum via (Jarvis et al. 2004; Schneider et al. 2005)

$$\langle M_{ap}^3 \rangle^{(ijk)}(\theta_{ap,1}, \theta_{ap,2}, \theta_{ap,3}) = \int \frac{d^2\ell_1}{(2\pi)^2} \int \frac{d^2\ell_2}{(2\pi)^2} \times B_{\kappa\kappa}^{(ijk)}(\ell_1, \ell_2, \ell_3) \hat{u}(\theta_{ap,1}\ell_1) \hat{u}(\theta_{ap,2}\ell_2) \hat{u}(\theta_{ap,3}\ell_3), \quad (26)$$

where $\ell_3 = |\ell_1 + \ell_2|$. Later, we differentiate between equal filter radii, for which $\theta_{ap,1} = \theta_{ap,2} = \theta_{ap,3}$ and non-equal filter radii, where the values of $\theta_{ap,i}$ are all allowed to vary.

## 2.5. Modelling COSEBIs

Schneider et al. (2010) introduced the complete orthogonal sets of E/B-integrals (COSEBIs), which are defined via the two-point shear statistics $\xi_{\pm}(\vartheta) = \langle \gamma_t \gamma_t'(\vartheta) \rangle \pm \langle \gamma_{\times} \gamma_{\times}' \rangle(\vartheta)$ on a finite angular range

$$E_n^{(ij)} = \frac{1}{2} \int_{\vartheta_{\min}}^{\vartheta_{\max}} d\vartheta \, \vartheta \left[ T_{+n}(\vartheta) \xi_{+}^{(ij)}(\vartheta) + T_{-n}(\vartheta) \xi_{-}^{(ij)}(\vartheta) \right], \quad (27)$$

$$B_n^{(ij)} = \frac{1}{2} \int_{\vartheta_{\min}}^{\vartheta_{\max}} d\vartheta \, \vartheta \left[ T_{+n}(\vartheta) \xi_{+}^{(ij)}(\vartheta) - T_{-n}(\vartheta) \xi_{-}^{(ij)}(\vartheta) \right], \quad (28)$$

where $T_{\pm n}(\vartheta)$ are filter functions with support in $[\vartheta_{\min}, \vartheta_{\max}]$ (Schneider et al. 2010). If

$$\begin{align}
\int_{\vartheta_{\min}}^{\vartheta_{\max}} d\vartheta \, \vartheta T_{+n}(\vartheta) J_0(\ell \vartheta) &= \int_{\vartheta_{\min}}^{\vartheta_{\max}} d\vartheta \, \vartheta T_{-n}(\vartheta) J_4(\ell \vartheta) \nonumber \\
&:= W_n(\ell), \tag{29}
\end{align}$$

where $J_{0,4}$ are the zeroth and fourth order Bessel function, then the COSEBIs offer the advantage of cleanly separating all well-defined E- and B-modes within the range $[\vartheta_{\min}, \vartheta_{\max}]$. This is not given, for instance, for second-order aperture mass statistics, as this would require information of $\xi_{\pm}$ over the full space. Analytically, the COSEBIs can be calculated from the E-mode power

spectrum $P_{\kappa\kappa}^{(ij)}(\ell)$ defined in Eq. (4) and a B-mode angular power spectra $P_{\kappa\kappa,B}^{(ij)}(\ell)$ as:

$$E_n^{(ij)} = \int_0^\infty \frac{d\ell \, \ell}{2\pi} P_{\kappa\kappa}^{(ij)}(\ell) W_n(\ell), \quad (30)$$

$$B_n^{(ij)} = \int_0^\infty \frac{d\ell \, \ell}{2\pi} P_{\kappa\kappa,B}^{(ij)}(\ell) W_n(\ell). \quad (31)$$

Since B-modes cannot be created by gravitational lensing directly, we neglected the modelling of $B_n$ for this work. Therefore, we focus, henceforth, only on the $E_n$ modes from the COSEBIs.

# 3. Measuring shear statistics

As discussed in H23 $\langle M_{ap}^3 \rangle$ can be estimated from data in three ways, which we review below. The first uses the convergence field, $\kappa$, the second uses the shear field, and the third uses correlation functions. The latter is used to measure the real data and compute the covariance matrix used for the real data analysis.

## 3.1. Measuring shear statistics from convergence maps

If the convergence field, $\kappa$, is available (e.g., from simulations), the easiest way to measure $M_{ap}$ is to use Eq. (21). If the convergence field does not contain masks, this estimator is also unbiased as long as a border of the size of the filter function is removed from the aperture mass field before calculating the spatial average. However, no boundaries need to be removed for unmasked full-sky convergence fields to get an unbiased estimator. To estimate the aperture mass from (full-sky) convergence maps, the maps are smoothed with the HEALPY<fnref n="1" /> (Zonca et al. 2019) function SMOOTHING. This function needs a beam window function created by the function BEAM2BL, which is determined from the corresponding $U_{\theta_{ap}}$ filter.

In the absence of B-modes the $E_n$ modes can also be calculated by using convolutions as (Schneider et al. 2010):

$$E_n^{(ij)} = \langle \kappa^{(i)}(\vartheta) \kappa_s^{(j)}(\vartheta; n) \rangle, \quad (32)$$

where

$$\kappa_s^{(i)}(\vartheta; n) = \frac{1}{2\pi} \int d^2\theta'_{ap} T_{+n}(|\vartheta - \vartheta'|) \kappa^{(i)}(\theta'_{ap}) \quad (33)$$

with $T_{+n}$ introduced in Eq. (29).

## 3.2. Measuring shear statistics from galaxy catalogues

In a real survey, the convergence is not observable and can be inferred only from the measured shear. However, estimating $M_{ap}$ from these reconstructed convergence maps is not accurate as the reconstructed convergence is necessarily smoothed and potentially also contains other systematic effects caused by masks or boundaries of the survey (Seitz & Schneider 1997).

However, as motivated by Eq. (22), the aperture mass can also be estimated from an unmasked shear field (e.g., for simulations). Given a galaxy catalogue, where galaxies are only found at specific positions such that the number of galaxies within an aperture varies on the sky, Eq. (22) needs to get modified to

$$\tilde{M}_{ap}^{(i)}(\vartheta; \theta_{ap}) = \left[ \sum_{i} Q_{\theta_{ap}}(|\vartheta - \vartheta_i|) \right]^{-1} \sum_{j} Q_{\theta_{ap}}(|\vartheta - \vartheta_j|) \varepsilon_{t,j}^{(i)}, \quad (34)$$

<footnote n="1">¹ http://healpix.sourceforge.net</footnote>
where the sum over the filter functions serves as a normalisation, $\epsilon_t$ are the observed galaxy ellipticities converted into their tangential component, and $\theta_i$ are their respective positions. We note that we sampled the galaxies on a grid using a cloud-in-cell method<fnref n="2" />, so that convolutions can determine the aperture masses. Since we used this approach solely for finite fields, we applied the same cut-off of $4\theta_{ap}^{max}$ to all aperture mass maps, where $\theta_{ap}^{max}$ is the largest filter radius. For both approaches, where the aperture mass is determined, the second- and third-order aperture statistics follow by multiplying the respective aperture mass maps pixel-wise and then taking the average of all pixel values.

### 3.3. Measuring shear statistics from correlation functions

The first two methods to estimate unbiased aperture statistics are not applicable to observed data with masks. Another disadvantage is the removal of the edges, which can be significant for a complex survey footprint, leading to decreased statistical certainty. The best method to estimate aperture shear statistics or the $E_n$ is by measuring them from two- and three-point shear correlation functions (Schneider et al. 2002, 2005; Jarvis et al. 2004). The advantage of correlation functions is that they can be estimated for any survey geometry. The measurement of the two-point shear correlation functions $\xi_{\pm}(\theta_{ap})$ is unbiased and easily performed by TREECORR (Jarvis et al. 2004) and the conversion to COSEBIs are given in Eqs. (27) and (28).

For the measurement of the aperture statistics, $\langle M_{ap}^3 \rangle$, we refer to our companion paper, namely, Porth et al. (2023). It describes an efficient estimation procedure of the natural components of the shear 3PCF (Schneider & Lombardi 2003), which is then transformed into aperture statistics. The corresponding equations for a tomographic set-up can be found in their Sect. 5.2.

## 4. Observed data

This analysis explores the fourth data release of KiDS (Kuijken et al. 2015, 2019; de Jong et al. 2015, 2017). The weak lensing data observed with the high-quality VST-OmegaCAM camera is public<fnref n="3" />. It is collectively called ‘KiDS-1000’ as it covers ~1000 deg<fnref n="2" /> of images, which is then reduced to an effective area of 777.4 deg<fnref n="2" /> after masking. The significant advantage, compared to previous weak lensing surveys and data releases, is its overlap with its partner survey VIKING (VISTA Kilo-degree Infrared Galaxy survey, Edge et al. 2013), which observes galaxy images at infrared wavelength. Therefore, galaxies were observed in nine optical and near-infrared bands, *u*, *g*, *r*, *i*, *Z*, *Y*, *J*, *H*, *Ks*, allowing for better control over redshift uncertainties (Hildebrandt et al. 2021, hereafter H21).

The KiDS-1000 cosmic shear catalogue is divided into five tomographic $z_{ph}$-bins, whose redshifts are calibrated using the self-organising map (SOM) method<fnref n="4" /> described in Wright et al. (2020). The redshift distributions of all five tomographic bins

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Redshift distribution of the five tomographic $z_{ph}$-bins of the KiDS-1000 data.</figcaption></figure>

are shown in Fig. 1, and were initially presented in H21. The residual systematic uncertainties on the redshift distributions are listed in Table 1 and are included in this work as nuisance parameters which we marginalise over. We note that they are correlated, which we account for by using their correlation matrix for the marginalisation. The galaxy shear ellipticities and their corresponding weights, $w$, are estimated by the lensfit tool (Miller et al. 2013; Fenech Conti et al. 2017; Kannawadi et al. 2019) and are described in more detail in Giblin et al. (2021). These shear-related systematic effects shift the $S_8$ parameter by (at most) $0.1\sigma$ when measured by cosmic shear two-point functions. The resulting systematics are stated in Table 1, where we marginalise over the shear multiplicative $m$-bias correction in the resulting posteriors.

## 5. Simulated data

We use several simulated data sets created to resemble the observed KiDS-1000 data, to validate our inference pipeline, study the impact of key systematic uncertainties, and forecast the expected KiDS-1000 analysis.

In particular, we used the full-sky gravitational lensing simulations described in Takahashi et al. (2017, hereafter T17) to generate data vectors and numerical covariance matrices. The cosmo-SLICS+IA simulations, described in Harnois-Déraps et al. (2022), were used to test the modelling of IA; while the Magneticum lensing simulations, first introduced in Hirschmann et al. (2014), were used to infuse Baryon feedback on the model whose strength we regulate with a free parameter in the posterior estimation.

### 5.1. Takahashi simulations

Since the T17 simulations are used in this series of previous works, we only reiterate the essential details here. The T17 simulations follow the non-linear evolution of $2048^3$ particles evolved in a large series of nested cosmological volumes with side length starting at $L = 450$ Mpc $h$ at low redshift, and increasing at higher redshift, resulting in 108 different full-sky realisations. These were produced by the Gadget-3 *N*-body code (Springel 2005) and are publicly available<fnref n="5" />. The cosmological parameters of the matter and vacuum energy density are fixed to $\Omega_m = 1 - \Omega_\Lambda = 0.279$, the baryon density parameter to $\Omega_b = 0.046$, the dimensionless Hubble constant to $h = 0.7$, the normalisation of the power spectrum to $\sigma_8 = 0.82$, and the spectral index to $n_s = 0.97$. The shear information of the T17 simulations is given in terms of 108 $\gamma$ and $\kappa$ full-sky realisations, where

<footnote n="2">² When placing the properties (shear) of a galaxy in a pixel grid, instead of shifting it to the pixel centre, we assume that the galaxy itself has the size of a pixel and distribute its properties to all neighbouring pixel centres weighted by their relative distance to the galaxy. This improves the accuracy on scales smaller than the pixel size.</footnote>

<footnote n="3">³ The KiDS data products are public and available through http://kids.strw.leidenuniv.nl/DR4</footnote>

<footnote n="4">⁴ By use of the nine-band photometry the SOM method allocates groups of galaxies to corresponding spectroscopic samples. If no matches are found, these galaxies are removed from the catalogue.</footnote>

<footnote n="5">⁵ T17 simulations: http://cosmo.phys.hirosaki-u.ac.jp/takahasi/allskyRaytracing/</footnote>
where the sum over the filter functions serves as a normalisation, $\epsilon_t$ are the observed galaxy ellipticities converted into their tangential component, and $\theta_i$ are their respective positions. We note that we sampled the galaxies on a grid using a cloud-in-cell method<fnref n="2" />, so that convolutions can determine the aperture masses. Since we used this approach solely for finite fields, we applied the same cut-off of $4\theta_{ap}^{max}$ to all aperture mass maps, where $\theta_{ap}^{max}$ is the largest filter radius. For both approaches, where the aperture mass is determined, the second- and third-order aperture statistics follow by multiplying the respective aperture mass maps pixel-wise and then taking the average of all pixel values.

### 3.3. Measuring shear statistics from correlation functions

The first two methods to estimate unbiased aperture statistics are not applicable to observed data with masks. Another disadvantage is the removal of the edges, which can be significant for a complex survey footprint, leading to decreased statistical certainty. The best method to estimate aperture shear statistics or the $E_n$ is by measuring them from two- and three-point shear correlation functions (Schneider et al. 2002, 2005; Jarvis et al. 2004). The advantage of correlation functions is that they can be estimated for any survey geometry. The measurement of the two-point shear correlation functions $\xi_{\pm}(\theta_{ap})$ is unbiased and easily performed by TREECORR (Jarvis et al. 2004) and the conversion to COSEBIs are given in Eqs. (27) and (28).

For the measurement of the aperture statistics, $\langle M_{ap}^3 \rangle$, we refer to our companion paper, namely, Porth et al. (2023). It describes an efficient estimation procedure of the natural components of the shear 3PCF (Schneider & Lombardi 2003), which is then transformed into aperture statistics. The corresponding equations for a tomographic set-up can be found in their Sect. 5.2.

## 4. Observed data

This analysis explores the fourth data release of KiDS (Kuijken et al. 2015, 2019; de Jong et al. 2015, 2017). The weak lensing data observed with the high-quality VST-OmegaCAM camera is public<fnref n="3" />. It is collectively called ‘KiDS-1000’ as it covers ~1000 deg<fnref n="2" /> of images, which is then reduced to an effective area of 777.4 deg<fnref n="2" /> after masking. The significant advantage, compared to previous weak lensing surveys and data releases, is its overlap with its partner survey VIKING (VISTA Kilo-degree Infrared Galaxy survey, Edge et al. 2013), which observes galaxy images at infrared wavelength. Therefore, galaxies were observed in nine optical and near-infrared bands, *u*, *g*, *r*, *i*, *Z*, *Y*, *J*, *H*, *Ks*, allowing for better control over redshift uncertainties (Hildebrandt et al. 2021, hereafter H21).

The KiDS-1000 cosmic shear catalogue is divided into five tomographic $z_{ph}$-bins, whose redshifts are calibrated using the self-organising map (SOM) method<fnref n="4" /> described in Wright et al. (2020). The redshift distributions of all five tomographic bins

<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Redshift distribution of the five tomographic $z_{ph}$-bins of the KiDS-1000 data.</figcaption></figure>

are shown in Fig. 1, and were initially presented in H21. The residual systematic uncertainties on the redshift distributions are listed in Table 1 and are included in this work as nuisance parameters which we marginalise over. We note that they are correlated, which we account for by using their correlation matrix for the marginalisation. The galaxy shear ellipticities and their corresponding weights, $w$, are estimated by the lensfit tool (Miller et al. 2013; Fenech Conti et al. 2017; Kannawadi et al. 2019) and are described in more detail in Giblin et al. (2021). These shear-related systematic effects shift the $S_8$ parameter by (at most) $0.1\sigma$ when measured by cosmic shear two-point functions. The resulting systematics are stated in Table 1, where we marginalise over the shear multiplicative $m$-bias correction in the resulting posteriors.

## 5. Simulated data

We use several simulated data sets created to resemble the observed KiDS-1000 data, to validate our inference pipeline, study the impact of key systematic uncertainties, and forecast the expected KiDS-1000 analysis.

In particular, we used the full-sky gravitational lensing simulations described in Takahashi et al. (2017, hereafter T17) to generate data vectors and numerical covariance matrices. The cosmo-SLICS+IA simulations, described in Harnois-Déraps et al. (2022), were used to test the modelling of IA; while the Magneticum lensing simulations, first introduced in Hirschmann et al. (2014), were used to infuse Baryon feedback on the model whose strength we regulate with a free parameter in the posterior estimation.

### 5.1. Takahashi simulations

Since the T17 simulations are used in this series of previous works, we only reiterate the essential details here. The T17 simulations follow the non-linear evolution of $2048^3$ particles evolved in a large series of nested cosmological volumes with side length starting at $L = 450$ Mpc $h$ at low redshift, and increasing at higher redshift, resulting in 108 different full-sky realisations. These were produced by the Gadget-3 *N*-body code (Springel 2005) and are publicly available<fnref n="5" />. The cosmological parameters of the matter and vacuum energy density are fixed to $\Omega_m = 1 - \Omega_\Lambda = 0.279$, the baryon density parameter to $\Omega_b = 0.046$, the dimensionless Hubble constant to $h = 0.7$, the normalisation of the power spectrum to $\sigma_8 = 0.82$, and the spectral index to $n_s = 0.97$. The shear information of the T17 simulations is given in terms of 108 $\gamma$ and $\kappa$ full-sky realisations, where

<footnote n="2">² When placing the properties (shear) of a galaxy in a pixel grid, instead of shifting it to the pixel centre, we assume that the galaxy itself has the size of a pixel and distribute its properties to all neighbouring pixel centres weighted by their relative distance to the galaxy. This improves the accuracy on scales smaller than the pixel size.</footnote>

<footnote n="3">³ The KiDS data products are public and available through http://kids.strw.leidenuniv.nl/DR4</footnote>

<footnote n="4">⁴ By use of the nine-band photometry the SOM method allocates groups of galaxies to corresponding spectroscopic samples. If no matches are found, these galaxies are removed from the catalogue.</footnote>

<footnote n="5">⁵ T17 simulations: http://cosmo.phys.hirosaki-u.ac.jp/takahasi/allsky raytracing/</footnote>
<table><caption>**Table 1.** Overview of the observational KiDS-1000 data.</caption><thead><tr><th>Name</th><th>n<sub>eff</sub>[arcmin<sup>-2</sup>]</th><th>δ(z)</th><th>σ<sub>ε</sub></th><th>m-bias ×10<sup>3</sup></th></tr></thead><tbody><tr><td>z<sub>ph</sub>-bin: 1</td><td>0.62</td><td>0.000 ± 0.0106</td><td>0.270</td><td>-10 ± 19</td></tr><tr><td>z<sub>ph</sub>-bin: 2</td><td>1.18</td><td>0.002 ± 0.0113</td><td>0.258</td><td>-9 ± 20</td></tr><tr><td>z<sub>ph</sub>-bin: 3</td><td>1.85</td><td>0.013 ± 0.0118</td><td>0.273</td><td>-11 ± 17</td></tr><tr><td>z<sub>ph</sub>-bin: 4</td><td>1.26</td><td>0.011 ± 0.0087</td><td>0.254</td><td>8 ± 12</td></tr><tr><td>z<sub>ph</sub>-bin: 5</td><td>1.31</td><td>-0.006 ± 0.0097</td><td>0.270</td><td>12 ± 10</td></tr></tbody></table>

**Notes.** The second column shows the effective number density, which also accounts for the lensfit weights. The third column shows the mean and uncertainty on the redshift bias taken from H21. The redshift bias uncertainties are correlated, which we accounted for by the correlation matrix given in H21. The fourth column displays the measured ellipticity dispersion per component, $\sigma_{\epsilon}$, as given in Giblin et al. (2021). The fifth column shows the shear multiplicative m-bias correction updated in van den Busch et al. (2022).

each realisation is divided into 38 ascending redshift slices. To reproduce the KiDS-1000 data for a given tomographic $z_{ph}$-bin shown in Fig. 1, we built the weighted average of the first 30 $\gamma$ and $\kappa$ redshift slices. The weight for each redshift slice is measured by integrating the $n(z)$ over the corresponding width of the redshift slice.

We decided on two approaches (convergence maps vs. galaxy ellipticity) to measure the data vectors and covariance matrices from the T17 simulations. While the first approach (Sect. 5.1.1) has the advantage that a large number of mock data is available to measure a reliable covariance matrix even for large data vectors, the second approach (Sect. 5.1.2) has the advantage that the exact galaxy positions are used, which implies that the holes (masks) match the data.

### 5.1.1. Convergence mock data

The first approach is to convolve the full sky convergence maps with the $T_{\pm,n}(\theta_{ap})$ or $U(\theta_{ap})$ filters, then multiply with corresponding other smoothed maps and then take the spatial average. For the covariance matrix, we divided the smoothed convergence maps into 48 sub-patches, where each patch has an area of 860 deg$^2$. Since the $\kappa$ maps are first convolved with the filter functions and the sub-patches have common edges, the individual patches are not fully independent from each other, which decreases the covariance matrix. However, we found in our testing that selecting only 18 sub-patches, that have no common edges, gives almost identical results. Shape noise is added to the convergence maps by drawing random numbers from a Gaussian distribution with a vanishing mean and a standard deviation, as follows:

$$ \sigma = \frac{\sigma_{\epsilon}}{\sqrt{n_{\text{eff}} A_{\text{pix}}}}, \qquad (35) $$

with the pixel area as $A_{pix} = 0.74 \text{ arcmin}^2$ (NSIDE = 4096), $n_{\text{eff}}$ as the effective galaxy number density that included the lensfit weights, and $\sigma_{\epsilon}$ as the shape noise contribution given in Table 1. With 48 sub-patches and 108 realisations, the covariance matrix is measured from 5184 mock data. As the actual KiDS-1000 area is roughly 777.4 deg$^2$, we rescaled the covariance by 860 deg$^2/777.4 \text{ deg}^2 \approx 1.1$. The reference data vector is measured from one T17 realisation with a resolution of $A_{pix} = 0.18 \text{ arcmin}^2$ (NSIDE = 8192) without shape noise.

### 5.1.2. Galaxy shear catalogue mock data

For real surveys with a complex topology, the convolved maps give biased results (Seitz & Schneider 1997). Therefore, we

decided on our second approach to measure our statistics from second- and third-order shear correlation functions. For this, we created galaxy shear catalogues from projected $\gamma$ and $\kappa$ fields by extracting the shear information only at the true positions of the observed galaxies in KiDS-1000. Since the correlation functions need to be measured to very small scales, we used realisations with a pixel resolution of $A_{pix} = 0.18 \text{ arcmin}^2$ (NSIDE = 8192). We checked that this resolution is sufficient by comparing the data vectors obtained from catalogues constructed from the same initial conditions on the reference resolution (NSIDE = 8192) and a higher resolution (NSIDE = 16 384), finding a maximum deviation of 1% for an aperture filter scale of 4'. To increase the number of mock data, we shifted the galaxy positions 18 times by 20° along the lines of constant declination and extracted the shear information at the new positions. Afterwards, the shifted galaxy positions and the corresponding lensing information are shifted back to the original footprint. The back shifting is done only for simplicity and has no physical reason. With this procedure, 1944 almost independent mock data were created, from which the covariance and the reference data vector were measured. To add shape noise, we combined the two-component reduced shear $g = \gamma/(1-\kappa)$ of each object with a shape noise contribution, $\epsilon^s$, to create observed ellipticities $\epsilon^{\text{obs}}$ (Seitz & Schneider 1997), as follows:

$$ \epsilon^{\text{obs}} = \frac{\epsilon^s + g}{1 + \epsilon^s g^*}. \qquad (36) $$

The quantities here are all complex numbers and the asterisk '*' indicates complex conjugation. To mock the $\epsilon^s$, we use the observed ellipticities $\epsilon_{\text{KiDS}}^{\text{obs}}$ of the KiDS-1000 data and randomly rotate them to erase the underlying correlated shear signal. This procedure offers the advantage that the resulting distribution of $\epsilon^s$ matches the distribution of $\epsilon_{\text{KiDS}}^{\text{obs}}$. The estimated mean and dispersion $\sigma_{\epsilon}$ are given in Table 1. Furthermore, when computing our statistics via correlation functions, we need to consider the corresponding weight $w$ for each shape measurement, which ensures that we use the correct effective number density. Since the lensing weights and the intrinsic ellipticities of source galaxies are correlated, adding the rotated ellipticities to the shear signal preserves this correlation.

### 5.1.3. Data vector measurement and modelling

We follow the analysis of Asgari et al. (2021, hereafter A21) as it is the fiducial cosmic shear analysis of the KiDS-1000$^6$. We only use the first five $E_n$-moments determined from two-point correlation functions, measured from 0.5 to 300' in 400 radial bins.

<footnote>$^6$ For updated KiDS-1000 cosmic shear analysis, we refer the reader to DES & KiDS Collaboration (2023) and Li et al. (2023).</footnote>
<table><caption>**Table 1.** Overview of the observational KiDS-1000 data.</caption><thead><tr><th>Name</th><th>n<sub>eff</sub>[arcmin<sup>-2</sup>]</th><th>δ(z)</th><th>σ<sub>ε</sub></th><th>m-bias ×10<sup>3</sup></th></tr></thead><tbody><tr><td>z<sub>ph</sub>-bin: 1</td><td>0.62</td><td>0.000 ± 0.0106</td><td>0.270</td><td>-10 ± 19</td></tr><tr><td>z<sub>ph</sub>-bin: 2</td><td>1.18</td><td>0.002 ± 0.0113</td><td>0.258</td><td>-9 ± 20</td></tr><tr><td>z<sub>ph</sub>-bin: 3</td><td>1.85</td><td>0.013 ± 0.0118</td><td>0.273</td><td>-11 ± 17</td></tr><tr><td>z<sub>ph</sub>-bin: 4</td><td>1.26</td><td>0.011 ± 0.0087</td><td>0.254</td><td>8 ± 12</td></tr><tr><td>z<sub>ph</sub>-bin: 5</td><td>1.31</td><td>-0.006 ± 0.0097</td><td>0.270</td><td>12 ± 10</td></tr></tbody></table>

**Notes.** The second column shows the effective number density, which also accounts for the lensfit weights. The third column shows the mean and uncertainty on the redshift bias taken from H21. The redshift bias uncertainties are correlated, which we accounted for by the correlation matrix given in H21. The fourth column displays the measured ellipticity dispersion per component, $\sigma_{\epsilon}$, as given in Giblin et al. (2021). The fifth column shows the shear multiplicative m-bias correction updated in van den Busch et al. (2022).

each realisation is divided into 38 ascending redshift slices. To reproduce the KiDS-1000 data for a given tomographic $z_{ph}$-bin shown in Fig. 1, we built the weighted average of the first 30 $\gamma$ and $\kappa$ redshift slices. The weight for each redshift slice is measured by integrating the $n(z)$ over the corresponding width of the redshift slice.

We decided on two approaches (convergence maps vs. galaxy ellipticity) to measure the data vectors and covariance matrices from the T17 simulations. While the first approach (Sect. 5.1.1) has the advantage that a large number of mock data is available to measure a reliable covariance matrix even for large data vectors, the second approach (Sect. 5.1.2) has the advantage that the exact galaxy positions are used, which implies that the holes (masks) match the data.

### 5.1.1. Convergence mock data

The first approach is to convolve the full sky convergence maps with the $T_{\pm,n}(\theta_{ap})$ or $U(\theta_{ap})$ filters, then multiply with corresponding other smoothed maps and then take the spatial average. For the covariance matrix, we divided the smoothed convergence maps into 48 sub-patches, where each patch has an area of 860 deg$^2$. Since the $\kappa$ maps are first convolved with the filter functions and the sub-patches have common edges, the individual patches are not fully independent from each other, which decreases the covariance matrix. However, we found in our testing that selecting only 18 sub-patches, that have no common edges, gives almost identical results. Shape noise is added to the convergence maps by drawing random numbers from a Gaussian distribution with a vanishing mean and a standard deviation, as follows:

$$ \sigma = \frac{\sigma_{\epsilon}}{\sqrt{n_{\text{eff}} A_{\text{pix}}}}, \qquad (35) $$

with the pixel area as $A_{pix} = 0.74 \text{ arcmin}^2$ (NSIDE = 4096), $n_{\text{eff}}$ as the effective galaxy number density that included the lensfit weights, and $\sigma_{\epsilon}$ as the shape noise contribution given in Table 1. With 48 sub-patches and 108 realisations, the covariance matrix is measured from 5184 mock data. As the actual KiDS-1000 area is roughly 777.4 deg$^2$, we rescaled the covariance by 860 deg$^2/777.4 \text{ deg}^2 \approx 1.1$. The reference data vector is measured from one T17 realisation with a resolution of $A_{pix} = 0.18 \text{ arcmin}^2$ (NSIDE = 8192) without shape noise.

### 5.1.2. Galaxy shear catalogue mock data

For real surveys with a complex topology, the convolved maps give biased results (Seitz & Schneider 1997). Therefore, we

decided on our second approach to measure our statistics from second- and third-order shear correlation functions. For this, we created galaxy shear catalogues from projected $\gamma$ and $\kappa$ fields by extracting the shear information only at the true positions of the observed galaxies in KiDS-1000. Since the correlation functions need to be measured to very small scales, we used realisations with a pixel resolution of $A_{pix} = 0.18 \text{ arcmin}^2$ (NSIDE = 8192). We checked that this resolution is sufficient by comparing the data vectors obtained from catalogues constructed from the same initial conditions on the reference resolution (NSIDE = 8192) and a higher resolution (NSIDE = 16384), finding a maximum deviation of 1% for an aperture filter scale of 4'. To increase the number of mock data, we shifted the galaxy positions 18 times by 20° along the lines of constant declination and extracted the shear information at the new positions. Afterwards, the shifted galaxy positions and the corresponding lensing information are shifted back to the original footprint. The back shifting is done only for simplicity and has no physical reason. With this procedure, 1944 almost independent mock data were created, from which the covariance and the reference data vector were measured. To add shape noise, we combined the two-component reduced shear $g = \gamma/(1-\kappa)$ of each object with a shape noise contribution, $\epsilon^s$, to create observed ellipticities $\epsilon^{\text{obs}}$ (Seitz & Schneider 1997), as follows:

$$ \epsilon^{\text{obs}} = \frac{\epsilon^s + g}{1 + \epsilon^s g^*}. \qquad (36) $$

The quantities here are all complex numbers and the asterisk '*' indicates complex conjugation. To mock the $\epsilon^s$, we use the observed ellipticities $\epsilon_{\text{KiDS}}^{\text{obs}}$ of the KiDS-1000 data and randomly rotate them to erase the underlying correlated shear signal. This procedure offers the advantage that the resulting distribution of $\epsilon^s$ matches the distribution of $\epsilon_{\text{KiDS}}^{\text{obs}}$. The estimated mean and dispersion $\sigma_{\epsilon}$ are given in Table 1. Furthermore, when computing our statistics via correlation functions, we need to consider the corresponding weight $w$ for each shape measurement, which ensures that we use the correct effective number density. Since the lensing weights and the intrinsic ellipticities of source galaxies are correlated, adding the rotated ellipticities to the shear signal preserves this correlation.

### 5.1.3. Data vector measurement and modelling

We follow the analysis of Asgari et al. (2021, hereafter A21) as it is the fiducial cosmic shear analysis of the KiDS-1000$^6$. We only use the first five $E_n$-moments determined from two-point correlation functions, measured from 0.5 to 300' in 400 radial bins.

<footnote>$^6$ For updated KiDS-1000 cosmic shear analysis, we refer the reader to DES & KiDS Collaboration (2023) and Li et al. (2023).</footnote>
The resulting $E_n$ are shown in Fig. 2, where we checked that the $B_n$ are consistent with zero. The $E_n$ data vector for $\gamma^7$ and $g$ is the mean of all 1944 T17 mock data; and for the $\kappa$ data vector, we used one realisation with resolution $A_{\text{pix}} = 0.18$ arcmin$^2$ (NSIDE = 8192). The analytical model was computed using OSMOSIS (Zuntz et al. 2015a).

The $\langle M_{\text{ap}}^3 \rangle$ data vectors are measured from the same mock data that are also used to compute the $E_n$. The analytical model pipeline is described in H23. In Fig. 3, we show $\langle M_{\text{ap}}^3 \rangle$ for aperture filter radii $\theta_{\text{ap}} \in \{4', 6', 8', 10', 14', 18', 22', 26', 32', 36'\}$, where we used only equal scales. We show them here for some specific $z_{\text{ph}}$-bin combinations. In Fig. B.2, we show $\langle M_{\text{ap}}^3 \rangle$ also for non-equal aperture filter radii $\theta_{\text{ap}} \in \{4', 8', 16', 32'\}$.

## 5.2. Cosmo-SLICS+IA

The cosmo-SLICS+IA simulations are cosmic shear mock galaxy catalogues infused with the non-linear alignment model of Bridle & King (2007), which is ideally suited for testing and validating our analytical IA model. They are based on N-body simulations of identical box size and particle density as the SLICS (Harnois-Déraps et al. 2018), which were already used in previous works of this series; we therefore only summarise the essential details.

The cosmology corresponds to the fiducial cosmo-SLICS model presented in Harnois-Déraps et al. (2019), with $\Omega_m = 0.2905$, $\Omega_\Lambda = 0.7095$, $\Omega_b = 0.0473$, $h = 0.6898$, $\sigma_8 = 0.836$, $w_0 = -1.0$ and $n_s = 0.969$. We use the full set of 50 simulated galaxy catalogues, each covering 100 deg$^2$ and reproducing the KiDS-1000 $n(z)$ (see Fig. 1) and galaxy number densities $n_{\text{eff}}$ specified in H21. As we use these mock data only to infuse IA effects into the T17 data vector, shape noise is not included. Following the methods described in Harnois-Déraps et al. (2022), the intrinsic ellipticity components $\epsilon_{1/2}^{\text{IA}}$ of these galaxies are computed as

$$ \epsilon_{1}^{\text{IA}} = -f_{\text{IA}}(\partial_{xx} - \partial_{yy})\phi, \quad \epsilon_{2}^{\text{IA}} = -2f_{\text{IA}}\partial_{xy}\phi, \qquad (37) $$

where $f_{\text{IA}}$ is defined in Eq. (7), and $\phi$ is the gravitational potential. The partial derivatives of the gravitational potential describe the Cartesian components of the projected tidal field tensors. The $\epsilon_{1,2}^{\text{IA}}$ terms were then combined with the noise-free cosmic shear signal using Eq. (36), resulting in an IA-contaminated weak lensing sample that is consistent with the NLA model. As these mock data cover a square patch without any masking, we used the methods described in Sect. 3.2 to estimate the aperture statistics. Furthermore, they were used only to validate the IA modelling, which we discuss in Appendix A.

## 5.3. Magneticum simulations

The feedback processes due to baryonic matter significantly affect the distribution of the LSS, such that the clustering of the matter is reduced on intra-cluster scales by up to tens of per cent (van Daalen et al. 2011). However, in a quantitative sense, this suppression is not well understood (Chisari et al. 2015). We use the Magneticum simulations (Hirschmann et al. 2014) to investigate the impact of baryonic feedback processes. Magneticum was run using GADGET3 code which is a more efficient version of Gadget 2 (Springel 2005) that includes modern smoothed particle hydrodynamics (Beck et al. 2016). The

dark matter particle mass is $6.9 \times 10^8 h^{-1} M_{\odot}$ and gas particle mass $1.4 \times 10^8 h^{-1} M_{\odot}$. The underlying matter fields were constructed from the Magneticum Pathfinder simulations$^8$ in the Run-2 with a comoving volume of side 352 $h^{-1}$ Mpc and Run-2b with a comoving volume of side 640 $h^{-1}$ Mpc, and are described in Hirschmann et al. (2014) and Ragagnin et al. (2017), respectively. The Magneticum simulations account for radiative cooling, star formation, supernovae and active galactic nuclei.

The Magneticum shear maps used in this work were first presented in Castro et al. (2018). Although the underlying creation of these shear maps is similar to the production of the cosmo-SLICS, meaning that the mock data follow the same redshift distribution and number density as the KiDS-1000 data, there are three main differences. First, only ten instead of 50 pseudo-independent light cones are available. Second, the galaxies are placed at the exact galaxy positions as in the observed data rather than randomly. Third, as these mock data are flat sky simulations and cover only an area of 100 deg$^2$, the KiDS-1000 footprint is divided into 18 regions. This results in 18 catalogues for each of the ten pseudo-independent light cones. Due to the non-trivial geometry of the masks, we used correlation functions to estimate the shear statistics.

To incorporate baryonic feedback processes into our modelling pipeline using the Magneticum simulations, we calculated the data vector, $d^{\text{DM}}$, using dark matter only and the data vector, $d^{\text{DM+BA}}$, where dark matter and baryons are included. With them we modified the model vector, $m$, to

$$ m' = \left[ 1 + A_{\text{ba}} \left( \frac{d^{\text{DM+BA}}}{d^{\text{DM}}} - 1 \right) \right] m, \qquad (38) $$

where we introduce the parameter $A_{\text{ba}}$. A value of zero for $A_{\text{ba}}$ means no baryonic feedback, and a value of one is exactly the strength of baryonic feedback as in the Magneticum. Furthermore, we interpolate between zero and one and extrapolate to five, which we chose arbitrarily. We note that for the joint analysis of second- and third-order statistics, we need at least extrapolation to $A_{\text{ba}} = 2$ (at 95% confidence), which approximately agrees with the baryonic strength of other hydrodynamical simulations (Martinet et al. 2021). The third-order analysis alone is bound by the prior $A_{\text{ba}} = 5$. Although extrapolating to these large values is probably incorrect, we introduced $A_{\text{ba}}$ only to give the model some flexibility to account for baryonic feedback effects without having a direct physical meaning such as the ejected gas.

# 6. Cosmological parameter inference methodology

In the following sections, we determine multiple posterior distributions using Markov chain Monte Carlo (MCMC) samplings for different model ingredients.

We varied the two cosmological parameters $\Omega_m$ and $S_8$ while fixing the remaining parameters to the T17 cosmology. Additionally, we varied the intrinsic alignment amplitude, $A_{\text{IA}}$, and the parameter $A_{\text{ba}}$ that accounts for the strengths of the response to baryonic feedback$^9$ and was described in Sect. 5.3. Lastly, we always marginalise over the shifts of the redshift distributions $\delta(z)$ given in H21 and the $m$-bias (Giblin et al. 2021) of all five tomographic bins shown in Table 1. For these nuisance parameters, we assume Gaussian priors and account for the fact that the

<footnote>$^8$ www.magneticum.org</footnote>

<footnote>$^9$ This parameter is used for both the $E_n$ and the $\langle M_{\text{ap}}^3 \rangle$ part of the data vector. We found that $A_{\text{ba}} = [0.0, 1.0, 3.0]$ roughly corressponds to $\log_{10} T_{\text{AGN}} = [7.1, 7.6, 8.5]$, which is the baryon parameter used in HMCODE2020.</footnote>

<footnote>$^7$ To compute the individual $\gamma$ data vectors with shape noise we have used Eq. (36) and replaced</footnote>
The resulting $E_n$ are shown in Fig. 2, where we checked that the $B_n$ are consistent with zero. The $E_n$ data vector for $\gamma^7$ and $g$ is the mean of all 1944 T17 mock data; and for the $\kappa$ data vector, we used one realisation with resolution $A_{\text{pix}} = 0.18 \, \text{arcmin}^2$ (NSIDE = 8192). The analytical model was computed using OSMOSIS (Zuntz et al. 2015a).

The $\langle M_{\text{ap}}^3 \rangle$ data vectors are measured from the same mock data that are also used to compute the $E_n$. The analytical model pipeline is described in H23. In Fig. 3, we show $\langle M_{\text{ap}}^3 \rangle$ for aperture filter radii $\theta_{\text{ap}} \in \{4', 6', 8', 10', 14', 18', 22', 26', 32', 36'\}$, where we used only equal scales. We show them here for some specific $z_{\text{ph}}$-bin combinations. In Fig. B.2, we show $\langle M_{\text{ap}}^3 \rangle$ also for non-equal aperture filter radii $\theta_{\text{ap}} \in \{4', 8', 16', 32'\}$.

## 5.2. Cosmo-SLICS+IA

The cosmo-SLICS+IA simulations are cosmic shear mock galaxy catalogues infused with the non-linear alignment model of Bridle & King (2007), which is ideally suited for testing and validating our analytical IA model. They are based on N-body simulations of identical box size and particle density as the SLICS (Harnois-Déraps et al. 2018), which were already used in previous works of this series; we therefore only summarise the essential details.

The cosmology corresponds to the fiducial cosmo-SLICS model presented in Harnois-Déraps et al. (2019), with $\Omega_m = 0.2905$, $\Omega_\Lambda = 0.7095$, $\Omega_b = 0.0473$, $h = 0.6898$, $\sigma_8 = 0.836$, $w_0 = -1.0$ and $n_s = 0.969$. We use the full set of 50 simulated galaxy catalogues, each covering 100 deg$^2$ and reproducing the KiDS-1000 $n(z)$ (see Fig. 1) and galaxy number densities $n_{\text{eff}}$ specified in H21. As we use these mock data only to infuse IA effects into the T17 data vector, shape noise is not included. Following the methods described in Harnois-Déraps et al. (2022), the intrinsic ellipticity components $\epsilon_{1/2}^{\text{IA}}$ of these galaxies are computed as

$$ \epsilon_{1}^{\text{IA}} = -f_{\text{IA}}(\partial_{xx} - \partial_{yy})\phi, \quad \epsilon_{2}^{\text{IA}} = -2f_{\text{IA}}\partial_{xy}\phi, \qquad (37) $$

where $f_{\text{IA}}$ is defined in Eq. (7), and $\phi$ is the gravitational potential. The partial derivatives of the gravitational potential describe the Cartesian components of the projected tidal field tensors. The $\epsilon_{1,2}^{\text{IA}}$ terms were then combined with the noise-free cosmic shear signal using Eq. (36), resulting in an IA-contaminated weak lensing sample that is consistent with the NLA model. As these mock data cover a square patch without any masking, we used the methods described in Sect. 3.2 to estimate the aperture statistics. Furthermore, they were used only to validate the IA modelling, which we discuss in Appendix A.

## 5.3. Magneticum simulations

The feedback processes due to baryonic matter significantly affect the distribution of the LSS, such that the clustering of the matter is reduced on intra-cluster scales by up to tens of per cent (van Daalen et al. 2011). However, in a quantitative sense, this suppression is not well understood (Chisari et al. 2015). We use the Magneticum simulations (Hirschmann et al. 2014) to investigate the impact of baryonic feedback processes. Magneticum was run using GADGET3 code which is a more efficient version of Gadget 2 (Springel 2005) that includes modern smoothed particle hydrodynamics (Beck et al. 2016). The

dark matter particle mass is $6.9 \times 10^8 h^{-1} M_{\odot}$ and gas particle mass $1.4 \times 10^8 h^{-1} M_{\odot}$. The underlying matter fields were constructed from the Magneticum Pathfinder simulations$^8$ in the Run-2 with a comoving volume of side 352 $h^{-1}$ Mpc and Run-2b with a comoving volume of side 640 $h^{-1}$ Mpc, and are described in Hirschmann et al. (2014) and Ragagnin et al. (2017), respectively. The Magneticum simulations account for radiative cooling, star formation, supernovae and active galactic nuclei.

The Magneticum shear maps used in this work were first presented in Castro et al. (2018). Although the underlying creation of these shear maps is similar to the production of the cosmo-SLICS, meaning that the mock data follow the same redshift distribution and number density as the KiDS-1000 data, there are three main differences. First, only ten instead of 50 pseudo-independent light cones are available. Second, the galaxies are placed at the exact galaxy positions as in the observed data rather than randomly. Third, as these mock data are flat sky simulations and cover only an area of 100 deg$^2$, the KiDS-1000 footprint is divided into 18 regions. This results in 18 catalogues for each of the ten pseudo-independent light cones. Due to the non-trivial geometry of the masks, we used correlation functions to estimate the shear statistics.

To incorporate baryonic feedback processes into our modelling pipeline using the Magneticum simulations, we calculated the data vector, $d^{\text{DM}}$, using dark matter only and the data vector, $d^{\text{DM+BA}}$, where dark matter and baryons are included. With them we modified the model vector, $m$, to

$$ m' = \left[ 1 + A_{\text{ba}} \left( \frac{d^{\text{DM+BA}}}{d^{\text{DM}}} - 1 \right) \right] m, \qquad (38) $$

where we introduce the parameter $A_{\text{ba}}$. A value of zero for $A_{\text{ba}}$ means no baryonic feedback, and a value of one is exactly the strength of baryonic feedback as in the Magneticum. Furthermore, we interpolate between zero and one and extrapolate to five, which we chose arbitrarily. We note that for the joint analysis of second- and third-order statistics, we need at least extrapolation to $A_{\text{ba}} = 2$ (at 95% confidence), which approximately agrees with the baryonic strength of other hydrodynamical simulations (Martinet et al. 2021). The third-order analysis alone is bound by the prior $A_{\text{ba}} = 5$. Although extrapolating to these large values is probably incorrect, we introduced $A_{\text{ba}}$ only to give the model some flexibility to account for baryonic feedback effects without having a direct physical meaning such as the ejected gas.

# 6. Cosmological parameter inference methodology

In the following sections, we determine multiple posterior distributions using Markov chain Monte Carlo (MCMC) samplings for different model ingredients.

We varied the two cosmological parameters $\Omega_m$ and $S_8$ while fixing the remaining parameters to the T17 cosmology. Additionally, we varied the intrinsic alignment amplitude, $A_{\text{IA}}$, and the parameter $A_{\text{ba}}$ that accounts for the strengths of the response to baryonic feedback$^9$ and was described in Sect. 5.3. Lastly, we always marginalise over the shifts of the redshift distributions $\delta(z)$ given in H21 and the $m$-bias (Giblin et al. 2021) of all five tomographic bins shown in Table 1. For these nuisance parameters, we assume Gaussian priors and account for the fact that the

<footnote>$^8$ www.magneticum.org</footnote>

<footnote>$^9$ This parameter is used for both the $E_n$ and the $\langle M_{\text{ap}}^3 \rangle$ part of the data vector. We found that $A_{\text{ba}} = [0.0, 1.0, 3.0]$ roughly corressponds to $\log_{10} T_{\text{AGN}} = [7.1, 7.6, 8.5]$, which is the baryon parameter used in HMCODE2020.</footnote>

<footnote>$^7$ To compute the individual $\gamma$ data vectors with shape noise we have used Eq. (36) and replaced</footnote>
<figure><img src="image_3.png" /><figcaption>Fig. 2. Measured and modelled $E_n$ vector for the first five moments in the T17 mock data. The green and blue dots are the mean of all 1944 KiDS-1000 mock data that are also used to compute the covariance matrix with a resolution $A_{pix} = 0.18 \text{ arcmin}^2$ and shape noise. The red dots represent the data vector measured from one full-sky T17 realisation with a resolution $A_{pix} = 0.18 \text{ arcmin}^2$ and no shape noise. The grey band indicates the expected uncertainty from the KiDS-1000 survey.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 3. Same as Fig. 2, but here the measured data and model are the $\langle M_{ap}^3 \rangle$ vector for equal-scale aperture filter radii $\theta_{ap} \in \{4', 6', 8', 10', 14', 18', 22', 26', 32', 36'\}$ in the T17 mock data for some selected $z_{ph}$-bin combinations. The full data set is shown in Fig. B.1.</figcaption></figure>

uncertainties of shifts of the redshift distributions are correlated using its correlation matrix. We give an overview of the priors in Table 2.

If the covariance matrix, $\tilde{C}$, is measured from simulations, it is a random variable. We followed the method from Percival et al. (2022), which leads to credible intervals that can also be interpreted as confidence intervals with approximately the same coverage probability. The posterior distribution of a model vector, $\mathbf{m}$, that depends on $n_p$ parameters, $\mathbf{p}$, if the covariance matrix $\tilde{C}$ is measured from $n_r$ mock data, is given by:

$$ \mathbf{P}[\mathbf{m}(\mathbf{p})|\mathbf{d}, \tilde{C}] \propto |\tilde{C}|^{-\frac{1}{2}} \left(1 + \frac{\chi^2}{n_r - 1}\right)^{-m/2}, \quad (39) $$

where $\mathbf{d}$ is the measured data vector and

$$ \chi^2 = [\mathbf{m}(\mathbf{p}) - \mathbf{d}]^\mathrm{T} \tilde{C}^{-1} [\mathbf{m}(\mathbf{p}) - \mathbf{d}]. \quad (40) $$

The power-law index $m$ is

$$ m = n_p + 2 + \frac{n_r - 1 + B(n_d - n_p)}{1 + B(n_d - n_p)}, \quad (41) $$

with $n_d$ being the number of data points and

$$ B = \frac{n_r - n_d - 2}{(n_r - n_d - 1)(n_r - n_d - 4)}. \quad (42) $$

To give a quantified value for the comparison of different modelling choices, we use the Figure of Merit (FoM), which we calculate as

$$ \text{FoM} = \frac{1}{\sqrt{\det \tilde{C}}}, \quad (43) $$
<figure><img src="image_4.png" /><figcaption>Fig. 2. Measured and modelled $E_n$ vector for the first five moments in the T17 mock data. The green and blue dots are the mean of all 1944 KiDS-1000 mock data that are also used to compute the covariance matrix with a resolution $A_{pix} = 0.18 \text{ arcmin}^2$ and shape noise. The red dots represent the data vector measured from one full-sky T17 realisation with a resolution $A_{pix} = 0.18 \text{ arcmin}^2$ and no shape noise. The grey band indicates the expected uncertainty from the KiDS-1000 survey.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 3. Same as Fig. 2, but here the measured data and model are the $\langle M_{ap}^3 \rangle$ vector for equal-scale aperture filter radii $\theta_{ap} \in \{4', 6', 8', 10', 14', 18', 22', 26', 32', 36'\}$ in the T17 mock data for some selected $z_{ph}$-bin combinations. The full data set is shown in Fig. B.1.</figcaption></figure>

uncertainties of shifts of the redshift distributions are correlated using its correlation matrix. We give an overview of the priors in Table 2.

If the covariance matrix, $\tilde{C}$, is measured from simulations, it is a random variable. We followed the method from Percival et al. (2022), which leads to credible intervals that can also be interpreted as confidence intervals with approximately the same coverage probability. The posterior distribution of a model vector, $\mathbf{m}$, that depends on $n_p$ parameters, $\mathbf{p}$, if the covariance matrix $\tilde{C}$ is measured from $n_r$ mock data, is given by:

$$ \mathbf{P}[\mathbf{m}(\mathbf{p})|\mathbf{d}, \tilde{C}] \propto |\tilde{C}|^{-\frac{1}{2}} \left(1 + \frac{\chi^2}{n_r - 1}\right)^{-m/2}, \quad (39) $$

where $\mathbf{d}$ is the measured data vector and

$$ \chi^2 = [\mathbf{m}(\mathbf{p}) - \mathbf{d}]^\mathrm{T} \tilde{C}^{-1} [\mathbf{m}(\mathbf{p}) - \mathbf{d}]. \quad (40) $$

The power-law index $m$ is

$$ m = n_p + 2 + \frac{n_r - 1 + B(n_d - n_p)}{1 + B(n_d - n_p)}, \quad (41) $$

with $n_d$ being the number of data points and

$$ B = \frac{n_r - n_d - 2}{(n_r - n_d - 1)(n_r - n_d - 4)}. \quad (42) $$

To give a quantified value for the comparison of different modelling choices, we use the Figure of Merit (FoM), which we calculate as

$$ \text{FoM} = \frac{1}{\sqrt{\det \tilde{C}}}, \quad (43) $$
<table><caption>**Table 2.** All varied parameters and their flat prior knowledge.</caption><thead><tr><th>Parameter</th><th>Validation</th><th>Real</th></tr></thead><tbody><tr><td>Ω<sub>m</sub></td><td>[0.10, 0.50]</td><td>[0.10, 0.75]</td></tr><tr><td>S<sub>8</sub></td><td>[0.6, 1.0]</td><td>[0.6, 1.0]</td></tr><tr><td>A<sub>IA</sub></td><td>[-1.8, 1.8]</td><td>[-1.5, 1.5]</td></tr><tr><td>A<sub>ba</sub></td><td>[0.0, 5.0]</td><td>[0.0, 5.0]</td></tr></tbody></table>

**Notes.** Uniformly distributed priors on the parameters used in our cosmological inferences. The priors given in Table 1 on the multiplicative shear *m*-bias and photometric redshift errors $\delta(z)$ are used both for the real data analysis and for the mock data analysis, where for the latter the expectation values are set to zero. The $\delta(z)$ for the sources follow a joint normal distribution with covariance matrix $C_{\delta(z)}$ shown in Fig. 6 of H21. The *m*-bias follows an uncorrelated normal distribution. We increased the upper prior range on $\Omega_m$ to 0.75 for the real data analysis to have $E_n$ posteriors that are not bound by the prior. The $A_{IA}$ of the real data analysis are narrower to improve the emulator accuracy. The upper bound of $A_{ba}$ is chosen arbitrarily and does play a role only for $\langle M_{ap}^3 \rangle$-only analysis.

where $C$ is the parameter covariance matrix of the $S_8-\Omega_m$ plane resulting from the MCMC process.

Depending on the used filter combinations, the modelling of $\langle M_{ap}^3 \rangle$ takes around 10–20 min, which stands as an obstacle to directly running an MCMC using the model. To circumvent this issue, we use the emulation tool contained in CosmoPower (Spurio Mancini et al. 2022). We trained the emulator on 3000 model points for the analysis based on mock data in the parameter space {$\Omega_m, S_8, A_{IA}$} and 5000 model points for the real data analysis in the parameter space {$\Omega_m, S_8, A_{IA}, \delta(z)$}, which we distributed in a Latin hypercube for the given prior range, and fixed all other parameters to the ones used in the T17 simulations. The *m*-bias and the $A_{ba}$ do not need to be emulated as they follow a simple correction formalism that is applied during the MCMC sampling. The $\delta(z)$ of each $z_{ph}$-bin for the training and testing are distributed uncorrelated between [-0.06, 0.06]. The accuracy of the emulator is tested by comparing the emulator prediction with the model at 500 independent points in the same parameter space. The fractional error of each vector element is smaller than 2% for $E_n$ and smaller than 5% for $\langle M_{ap}^3 \rangle$, which is well within the accuracy of the HMCODE2020 or the BIHALOFIT model itself. We used a Metropolis-Hastings sampler for MCMC<fnref n="10" />, where we used 1000 walkers running each 20 000 steps and cut the first 2000 steps away to ensure that the posteriors are not biased by the burn-in phase.

## 7. Reduction of data vector combinations and filter radii

Generally, as long as the covariance is converged, the more information is used, the higher the constraining power. However, for third-order statistics with four different filter radii and five tomographic bins, the data vector contains 700 elements. To obtain a reliable covariance matrix, roughly ~10 times the dimension of the data vector is needed. This implies that modelling and measuring the model/data vector and covariance matrix are time-consuming and require thousands of simulations. Therefore, it is inevitable for this and for future analyses to compress the data vector without substantial information loss. We note that this section is only concerned with reducing the dimension of the

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Posterior distribution for data vectors and covariance measured from the T17 convergence maps catalogues. Here only specific parts of the $\langle M_{ap}^3 \rangle$ data vector are used. Here 'only equal $z_{ph}$-bins' means that only the auto tomographic bins are used. Lastly, all non-equal filter radii are discarded for 'only equal filter radii'.</figcaption></figure>

data vector and is not meant as a proper forecast, which we do in the following sections. The spatial resolution of the mock data used to compute the covariance induces an error smaller than <1% in the $E_n$ covariance matrix. However, as we only reduce the elements of $\langle M_{ap}^3 \rangle$, while using all combinations of the $E_n$, this is unproblematic for the purpose of this section. We note again that the data vector was measured from one realisation with a resolution of $A_{pix} = 0.18 \text{ arcmin}^2$, and that we use the $\kappa$ based covariance for this section to ensure the covariance matrix is converged even for the data vector with the largest dimension. As a first check, we show in Fig. 4 multiple element choices that discard a significant part of the $\langle M_{ap}^3 \rangle$ data vector. Using only the five auto-tomographic bins is not a good choice for reduction, as seen in the fact that the FoM is reduced by 32%. However, using only equal-scale aperture radii is a better reduction, reducing the data vector to 28% while reducing the FoM by only ≈8%. Using only equal-scale aperture radii also has the advantage of being modelled faster. Therefore, we continue using only the equal-scale aperture radii for the rest of this work, which reduces the dimension of the data vector from 775 to 215.

Next, we investigated the choice of aperture filter radii. The aperture filter radii under consideration are $\theta_{ap} \in \{4', 6', 8', 10', 14', 18', 22', 26', 32', 36'\}$. The resulting posteriors are displayed in Fig. 5, which reveals that using only the large filter radii above >22' have the worst constraining power. However, using only the four smallest filter radii is also not the best choice, as they are largely correlated and miss the large-scale information. The best choice of filter function is to use a filter from each range such as $\theta_{ap} \in \{4', 8', 14', 32'\}$.

Next, we considered the full data vector, meaning the $\langle M_{ap}^3 \rangle$, and the $E_n$ for the next compression strategy. The idea is to decrease the number of elements by considering only those with the highest constraining power on $S_8$. We start with the element

<footnote n="10">¹⁰ The code we made use of can be found here: https://github.com/justinsling/affine</footnote>
<table><caption>**Table 2.** All varied parameters and their flat prior knowledge.</caption><thead><tr><th>Parameter</th><th>Validation</th><th>Real</th></tr></thead><tbody><tr><td>Ω<sub>m</sub></td><td>[0.10, 0.50]</td><td>[0.10, 0.75]</td></tr><tr><td>S<sub>8</sub></td><td>[0.6, 1.0]</td><td>[0.6, 1.0]</td></tr><tr><td>A<sub>IA</sub></td><td>[-1.8, 1.8]</td><td>[-1.5, 1.5]</td></tr><tr><td>A<sub>ba</sub></td><td>[0.0, 5.0]</td><td>[0.0, 5.0]</td></tr></tbody></table>

**Notes.** Uniformly distributed priors on the parameters used in our cosmological inferences. The priors given in Table 1 on the multiplicative shear *m*-bias and photometric redshift errors $\delta(z)$ are used both for the real data analysis and for the mock data analysis, where for the latter the expectation values are set to zero. The $\delta(z)$ for the sources follow a joint normal distribution with covariance matrix $C_{\delta(z)}$ shown in Fig. 6 of H21. The *m*-bias follows an uncorrelated normal distribution. We increased the upper prior range on $\Omega_m$ to 0.75 for the real data analysis to have $E_n$ posteriors that are not bound by the prior. The $A_{IA}$ of the real data analysis are narrower to improve the emulator accuracy. The upper bound of $A_{ba}$ is chosen arbitrarily and does play a role only for $\langle M_{ap}^3 \rangle$-only analysis.

where $C$ is the parameter covariance matrix of the $S_8-\Omega_m$ plane resulting from the MCMC process.

Depending on the used filter combinations, the modelling of $\langle M_{ap}^3 \rangle$ takes around 10–20 min, which stands as an obstacle to directly running an MCMC using the model. To circumvent this issue, we use the emulation tool contained in CosmoPower (Spurio Mancini et al. 2022). We trained the emulator on 3000 model points for the analysis based on mock data in the parameter space {$\Omega_m, S_8, A_{IA}$} and 5000 model points for the real data analysis in the parameter space {$\Omega_m, S_8, A_{IA}, \delta(z)$}, which we distributed in a Latin hypercube for the given prior range, and fixed all other parameters to the ones used in the T17 simulations. The *m*-bias and the $A_{ba}$ do not need to be emulated as they follow a simple correction formalism that is applied during the MCMC sampling. The $\delta(z)$ of each $z_{ph}$-bin for the training and testing are distributed uncorrelated between [-0.06, 0.06]. The accuracy of the emulator is tested by comparing the emulator prediction with the model at 500 independent points in the same parameter space. The fractional error of each vector element is smaller than 2% for $E_n$ and smaller than 5% for $\langle M_{ap}^3 \rangle$, which is well within the accuracy of the HMCODE2020 or the BIHALOFIT model itself. We used a Metropolis-Hastings sampler for MCMC<fnref n="10" />, where we used 1000 walkers running each 20 000 steps and cut the first 2000 steps away to ensure that the posteriors are not biased by the burn-in phase.

## 7. Reduction of data vector combinations and filter radii

Generally, as long as the covariance is converged, the more information is used, the higher the constraining power. However, for third-order statistics with four different filter radii and five tomographic bins, the data vector contains 700 elements. To obtain a reliable covariance matrix, roughly ~10 times the dimension of the data vector is needed. This implies that modelling and measuring the model/data vector and covariance matrix are time-consuming and require thousands of simulations. Therefore, it is inevitable for this and for future analyses to compress the data vector without substantial information loss. We note that this section is only concerned with reducing the dimension of the

<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Posterior distribution for data vectors and covariance measured from the T17 convergence maps catalogues. Here only specific parts of the $\langle M_{ap}^3 \rangle$ data vector are used. Here 'only equal $z_{ph}$-bins' means that only the auto tomographic bins are used. Lastly, all non-equal filter radii are discarded for 'only equal filter radii'.</figcaption></figure>

data vector and is not meant as a proper forecast, which we do in the following sections. The spatial resolution of the mock data used to compute the covariance induces an error smaller than <1% in the $E_n$ covariance matrix. However, as we only reduce the elements of $\langle M_{ap}^3 \rangle$, while using all combinations of the $E_n$, this is unproblematic for the purpose of this section. We note again that the data vector was measured from one realisation with a resolution of $A_{pix} = 0.18 \text{ arcmin}^2$, and that we use the $\kappa$ based covariance for this section to ensure the covariance matrix is converged even for the data vector with the largest dimension. As a first check, we show in Fig. 4 multiple element choices that discard a significant part of the $\langle M_{ap}^3 \rangle$ data vector. Using only the five auto-tomographic bins is not a good choice for reduction, as seen in the fact that the FoM is reduced by 32%. However, using only equal-scale aperture radii is a better reduction, reducing the data vector to 28% while reducing the FoM by only ≈8%. Using only equal-scale aperture radii also has the advantage of being modelled faster. Therefore, we continue using only the equal-scale aperture radii for the rest of this work, which reduces the dimension of the data vector from 775 to 215.

Next, we investigated the choice of aperture filter radii. The aperture filter radii under consideration are $\theta_{ap} \in \{4', 6', 8', 10', 14', 18', 22', 26', 32', 36'\}$. The resulting posteriors are displayed in Fig. 5, which reveals that using only the large filter radii above >22' have the worst constraining power. However, using only the four smallest filter radii is also not the best choice, as they are largely correlated and miss the large-scale information. The best choice of filter function is to use a filter from each range such as $\theta_{ap} \in \{4', 8', 14', 32'\}$.

Next, we considered the full data vector, meaning the $\langle M_{ap}^3 \rangle$, and the $E_n$ for the next compression strategy. The idea is to decrease the number of elements by considering only those with the highest constraining power on $S_8$. We start with the element

<footnote n="10">¹⁰ The code we made use of can be found here: https://github.com/justinsling/affine</footnote>
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Posterior distribution for different filter radii combinations if the data vector and covariances are measured from the T17 convergence maps.</figcaption></figure>

<table><caption>**Table 3.** Overview of the $\langle M_{\text{ap}}^3 \rangle$ data vector choices and their relative constraining power.</caption><thead><tr><th>Choice</th><th>Figure</th><th>rel. FoM</th><th>n<sub>d</sub></th></tr></thead><tbody><tr><td>All combinations</td><td>4</td><td>1.00</td><td>775</td></tr><tr><td>Only equal z<sub>ph</sub>-bins</td><td>4</td><td>0.68</td><td>175</td></tr><tr><td>θ<sub>ap</sub> ∈ {4', 8', 16', 32'}</td><td>4</td><td>0.92</td><td>215</td></tr><tr><td>θ<sub>ap</sub> ∈ {4', 8', 14', 32'}</td><td>5</td><td>0.92</td><td>215</td></tr><tr><td>θ<sub>ap</sub> ∈ {4', 14'}</td><td>–</td><td>0.87</td><td>145</td></tr><tr><td>θ<sub>ap</sub> ∈ {4', 14', 32'}</td><td>–</td><td>0.89</td><td>180</td></tr><tr><td>θ<sub>ap</sub> ∈ {4', 6', 8', 10'}</td><td>5</td><td>0.88</td><td>215</td></tr><tr><td>θ<sub>ap</sub> ∈ {10', 14', 18', 22'}</td><td>5</td><td>0.72</td><td>215</td></tr><tr><td>θ<sub>ap</sub> ∈ {22', 26', 32', 36'}</td><td>5</td><td>0.55</td><td>215</td></tr><tr><td>All equal-filter radii</td><td>–</td><td>0.92</td><td>425</td></tr><tr><td>Best 100 elements</td><td>–</td><td>0.80</td><td>100</td></tr><tr><td>Best 215 elements</td><td>–</td><td>0.89</td><td>200</td></tr></tbody></table>

**Notes.** For the measurement of the relative FoM the individual FoM is divided by the FoM of the $E_n + \langle M_{\text{ap}}^3 \rangle$ of the ‘all combinations’ case. From the second row onwards, only equal-scale filter radii are used.

with the highest S/N, then consecutively add those vector elements that maximise the $S_8$ Fisher information content. This Fisher information is calculated as (Tegmark et al. 1997)

$$F_{S_8} = \left( \frac{\partial \mathbf{m}(S_8)}{\partial S_8} \right)^T C^{-1} \left( \frac{\partial \mathbf{m}(S_8)}{\partial S_8} \right), \quad (44)$$

where the partial derivatives are computed with a five-point stencil beam (Fornberg 1988),

$$\frac{\partial \mathbf{m}(S_8)}{\partial S_8} \approx \frac{-\mathbf{m}(S_8^{++}) + 8\mathbf{m}(S_8^+) - 8\mathbf{m}(S_8^-) + \mathbf{m}(S_8^-)}{12\Delta S_8}, \quad (45)$$

where $S_8^{\pm} = S_8^{\text{T17}} \pm 0.02 S_8^{\text{T17}}$ and $S_8^{\pm\pm} = S_8^{\text{T17}} \pm 0.04 S_8^{\text{T17}}$. For this analysis, we used only the equal-scale filter radii

$\theta_{\text{ap}} \in \{4', 6', 8', 10', 14', 18', 22', 26', 32', 36'\}$. The first ~200 elements are sufficient to get converged posteriors. It is also interesting to see which elements help increase the constraining power. Unsurprisingly, the first elements are all COSEBIs, but among the first 100, approximately half are $\langle M_{\text{ap}}^3 \rangle$ elements. Furthermore, it is interesting to see that cross-tomographic bin elements are more likely to be selected by our method, which is expected because they have a higher S/N than auto-tomographic bin elements. Nevertheless, we also observe that using equal-scale filter radii $\theta_{\text{ap}} \in \{4', 8', 14', 32'\}$ results in a better FoM then using the best 215 elements. This is likely because we optimised only the $S_8$ parameter here while fixing all others. Therefore, we also checked for the equal scale filter case if a principal component analysis (PCA) applied to the covariance matrix performs better in compressing the data. However, it needs more elements to get the same constraining power as our FoM maximiser. This is probably because a PCA considers only the covariance matrix and ignores the derivatives, meaning that a PCA is not necessarily sensitive to cosmology.

Finally, we give in Table 3 an overview of the FoM for the $\Omega_m-S_8$ plane and the size of the data vector for some more element choices. As the covariance matrix is measured from 5184 mock data, we can assume that for all data vector dimensions under consideration, the covariance matrix is converged. Nevertheless, for the next sections, we use the covariance matrix measured from 1944 galaxy shear catalogues, limiting us to a maximum dimension of our data vector ~200 elements. Given the investigations in this section, we restrict the further sections to the case where we use all $E_n$ elements and all $\langle M_{\text{ap}}^3 \rangle$ elements with the equal-scale filter radii $\theta_{\text{ap}} \in \{4', 8', 14', 32'\}$, as this resulted in the best FoM for the restricted dimension of the data vector.

## 8. Validation of data vector estimator

In a real survey analysis, two further difficulties arise. The first issue is that the lensing information is not given in terms of convergence maps but by point estimates (galaxies). The finite area where galaxies are measured implies that no lensing information is available outside that area. It is, therefore, necessary to measure the statistic from correlation functions described in Sect. 3.3. Second, these point estimates are the reduced shear $g = \gamma/(1 - \kappa)$, which increases the signal and therefore needs to be accounted for. To correct these effects, we measured data vectors without shape noise and with the largest available resolution $A_{\text{pix}} = 0.05 \text{ arcmin}^2$. Since these data vectors result only from one realisation without shape noise, we can expect that the ratio is a good approximation for the reduced shear effect. Although the deviations are small (as seen in Figs. A.3 and A.4), we decided to scale the model vectors by the ratio of $g$ and $\gamma$ data vector.

The resulting posteriors are shown in Fig. 6. Our first observation is that similar to the finding of H23, combining second-order with third-order shear statistics significantly improves the constraining power on $S_8$ and $\Omega_m$. Compared to the $E_n$-only, the $S_8-\Omega_m$ FoM increases by 93% and for the $\langle M_{\text{ap}}^3 \rangle$-only by 233%. For these improvements, we have not considered that the posteriors of the individual statistics are bound by the priors on $\Omega_m$, meaning that the improvements are lower bounds. Compared to the analysis based on $\kappa$ mock data (see Fig. 5), the posteriors are broader because the covariance based on $\kappa$ maps also uses information outside the patches since the boundaries were not removed. This decreases the variance between the patches. A further difference is that the $\kappa$ analysis is not subject to masks
<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Posterior distribution for different filter radii combinations if the data vector and covariances are measured from the T17 convergence maps.</figcaption></figure>

<table><caption>**Table 3.** Overview of the $\langle M_{\text{ap}}^3 \rangle$ data vector choices and their relative constraining power.</caption><thead><tr><th>Choice</th><th>Figure</th><th>rel. FoM</th><th>n<sub>d</sub></th></tr></thead><tbody><tr><td>All combinations</td><td>4</td><td>1.00</td><td>775</td></tr><tr><td>Only equal z<sub>ph</sub>-bins</td><td>4</td><td>0.68</td><td>175</td></tr><tr><td>θ<sub>ap</sub> ∈ {4', 8', 16', 32'}</td><td>4</td><td>0.92</td><td>215</td></tr><tr><td>θ<sub>ap</sub> ∈ {4', 8', 14', 32'}</td><td>5</td><td>0.92</td><td>215</td></tr><tr><td>θ<sub>ap</sub> ∈ {4', 14'}</td><td>—</td><td>0.87</td><td>145</td></tr><tr><td>θ<sub>ap</sub> ∈ {4', 14', 32'}</td><td>—</td><td>0.89</td><td>180</td></tr><tr><td>θ<sub>ap</sub> ∈ {4', 6', 8', 10'}</td><td>5</td><td>0.88</td><td>215</td></tr><tr><td>θ<sub>ap</sub> ∈ {10', 14', 18', 22'}</td><td>5</td><td>0.72</td><td>215</td></tr><tr><td>θ<sub>ap</sub> ∈ {22', 26', 32', 36'}</td><td>5</td><td>0.55</td><td>215</td></tr><tr><td>All equal-filter radii</td><td>—</td><td>0.92</td><td>425</td></tr><tr><td>Best 100 elements</td><td>—</td><td>0.80</td><td>100</td></tr><tr><td>Best 215 elements</td><td>—</td><td>0.89</td><td>200</td></tr></tbody></table>

**Notes.** For the measurement of the relative FoM the individual FoM is divided by the FoM of the $E_n + \langle M_{\text{ap}}^3 \rangle$ of the ‘all combinations’ case. From the second row onwards, only equal-scale filter radii are used.

with the highest S/N, then consecutively add those vector elements that maximise the $S_8$ Fisher information content. This Fisher information is calculated as (Tegmark et al. 1997)

$$F_{S_8} = \left( \frac{\partial \mathbf{m}(S_8)}{\partial S_8} \right)^T C^{-1} \left( \frac{\partial \mathbf{m}(S_8)}{\partial S_8} \right), \quad (44)$$

where the partial derivatives are computed with a five-point stencil beam (Fornberg 1988),

$$\frac{\partial \mathbf{m}(S_8)}{\partial S_8} \approx \frac{-\mathbf{m}(S_8^{++}) + 8\mathbf{m}(S_8^+) - 8\mathbf{m}(S_8^-) + \mathbf{m}(S_8^-)}{12\Delta S_8}, \quad (45)$$

where $S_8^{\pm} = S_8^{\text{T17}} \pm 0.02 S_8^{\text{T17}}$ and $S_8^{\pm\pm} = S_8^{\text{T17}} \pm 0.04 S_8^{\text{T17}}$. For this analysis, we used only the equal-scale filter radii

$\theta_{\text{ap}} \in \{4', 6', 8', 10', 14', 18', 22', 26', 32', 36'\}$. The first ~200 elements are sufficient to get converged posteriors. It is also interesting to see which elements help increase the constraining power. Unsurprisingly, the first elements are all COSEBIs, but among the first 100, approximately half are $\langle M_{\text{ap}}^3 \rangle$ elements. Furthermore, it is interesting to see that cross-tomographic bin elements are more likely to be selected by our method, which is expected because they have a higher S/N than auto-tomographic bin elements. Nevertheless, we also observe that using equal-scale filter radii $\theta_{\text{ap}} \in \{4', 8', 14', 32'\}$ results in a better FoM then using the best 215 elements. This is likely because we optimised only the $S_8$ parameter here while fixing all others. Therefore, we also checked for the equal scale filter case if a principal component analysis (PCA) applied to the covariance matrix performs better in compressing the data. However, it needs more elements to get the same constraining power as our FoM maximiser. This is probably because a PCA considers only the covariance matrix and ignores the derivatives, meaning that a PCA is not necessarily sensitive to cosmology.

Finally, we give in Table 3 an overview of the FoM for the $\Omega_m-S_8$ plane and the size of the data vector for some more element choices. As the covariance matrix is measured from 5184 mock data, we can assume that for all data vector dimensions under consideration, the covariance matrix is converged. Nevertheless, for the next sections, we use the covariance matrix measured from 1944 galaxy shear catalogues, limiting us to a maximum dimension of our data vector ~200 elements. Given the investigations in this section, we restrict the further sections to the case where we use all $E_n$ elements and all $\langle M_{\text{ap}}^3 \rangle$ elements with the equal-scale filter radii $\theta_{\text{ap}} \in \{4', 8', 14', 32'\}$, as this resulted in the best FoM for the restricted dimension of the data vector.

## 8. Validation of data vector estimator

In a real survey analysis, two further difficulties arise. The first issue is that the lensing information is not given in terms of convergence maps but by point estimates (galaxies). The finite area where galaxies are measured implies that no lensing information is available outside that area. It is, therefore, necessary to measure the statistic from correlation functions described in Sect. 3.3. Second, these point estimates are the reduced shear $g = \gamma/(1 - \kappa)$, which increases the signal and therefore needs to be accounted for. To correct these effects, we measured data vectors without shape noise and with the largest available resolution $A_{\text{pix}} = 0.05 \text{ arcmin}^2$. Since these data vectors result only from one realisation without shape noise, we can expect that the ratio is a good approximation for the reduced shear effect. Although the deviations are small (as seen in Figs. A.3 and A.4), we decided to scale the model vectors by the ratio of $g$ and $\gamma$ data vector.

The resulting posteriors are shown in Fig. 6. Our first observation is that similar to the finding of H23, combining second-order with third-order shear statistics significantly improves the constraining power on $S_8$ and $\Omega_m$. Compared to the $E_n$-only, the $S_8-\Omega_m$ FoM increases by 93% and for the $\langle M_{\text{ap}}^3 \rangle$-only by 233%. For these improvements, we have not considered that the posteriors of the individual statistics are bound by the priors on $\Omega_m$, meaning that the improvements are lower bounds. Compared to the analysis based on $\kappa$ mock data (see Fig. 5), the posteriors are broader because the covariance based on $\kappa$ maps also uses information outside the patches since the boundaries were not removed. This decreases the variance between the patches. A further difference is that the $\kappa$ analysis is not subject to masks
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Posterior distribution for data vectors and covariance measured from the T17 galaxy shear catalogues. The covariance matrix and reference data vector were measured from 1944 noisy g mock data.</figcaption></figure>

and uses a slightly larger effective area. Although we rescaled the covariance to correct for the different effective areas, we found in Linke et al. (2023) that this rescaling is not necessarily accurate. Lastly, we notice that our modelling within the KiDS-1000 uncertainty is accurate, which we quantified by measuring the shift of the maximum of the posterior (MP) from the true values in the matter density parameter $\Delta\Omega_m < 0.02 \sigma_{\Omega_m}$ and in the clustering amplitude $\Delta S_8 < 0.05 \sigma_{S_8}$ with respect to the averaged noisy mock data vector and the KiDS-1000 uncertainty. We define the MP as the maximum of the one-dimensional marginal distributions.

## 9. Cosmological results

Finally, we are ready to present the first cosmological constraints from the KiDS-1000 data using second- and third-order shear statistics, displayed in Figs. 7 and 8. The data vectors are described in more detail in A21 for the $E_n$, and the $\langle M_{\text{ap}}^3 \rangle$ in Porth et al. (2023). Given the fact that the covariance matrix is measured only from 1944 realisations using the reduced shear g, we decide to build our data and model vector from all $E_n$ modes and $\langle M_{\text{ap}}^3 \rangle$ with all equal filter radii $\theta_{\text{ap}} \in \{4', 8', 14', 32'\}$. To control whether the model accurately fits the data, we minimised the $\chi^2_{\text{real}}$ from the real data and the $\chi^2_{\text{T17}}$ from each of the 1944 T17 mock data used to compute the covariance matrix. To estimate the probability of measuring a $\chi^2 > \chi^2_{\text{real}}$ (p-value), we counted the number of $\chi^2_{\text{T17}}$ that are greater than $\chi^2_{\text{real}}$ divided by 1944. The resulting p-value for the $E_n$-only, the $\langle M_{\text{ap}}^3 \rangle$-only, and the combination are given in Table 4. The resulting p-values for all combinations are better than 0.1, indicating that our covariance is well matched to the observed data and our model is accurate enough to describe the data. Our maximum posterior $\chi^2$ value increases to 81 if we swap our $E_n$ covariance matrix to the analytical expression in Joachimi et al. (2021) calculated at the MP parameter values in A21. This is unsurprising as the numeri-

cal covariance is computed at the T17 cosmology, giving a signal larger than that computed at the MP of A21. Furthermore, our covariance matrix is measured from reduced shear mock data that slightly increases the covariance matrix. We also notice that the model and the data seem inconsistent for some $z_{\text{ph}}$-bin combinations. However, adjacent COSEBI modes are highly correlated, so visually inspecting the model's goodness of fit to the data is misleading. We refer to A21 for further details on this discrepancy.

We show the resulting posteriors in Fig. 9, where we marginalised over the shift in the redshift distribution and multiplicative shear correction, both stated in Table 1. We improve the constraints on $S_8$ by 23% and on $\Omega_m$ by 47% compared to the $E_n$-only case. This shows how powerful a combined analysis of the second- and third-order shear statistics is. The constraints on $A_{\text{IA}}$ and $A_{\text{ba}}$ are basically untouched, showing that $\langle M_{\text{ap}}^3 \rangle$ is not helpful for constraining these nuisance parameters.

Compared to the maximum of the one-dimensional marginal distributions constraints of $E_n$ measurements given in A21 ($S_8 = 0.758_{-0.026}^{+0.017}$ and $\Omega_m = 0.253_{-0.074}^{+0.088}$), we have slightly larger constraints when using only the $E_n$. This is because we use a numerical covariance, which is larger than the analytical one due to the underlying cosmology and the fact that we model the reduced shear effect. Furthermore, we note that A21 varied $\Omega_{\text{cdm}}h^2$, $\Omega_b h^2$ and h, while we fixed h and $\Omega_b$ and varied only $\Omega_m$. We also find that if we use the same pipeline as A21 but allow larger $\Omega_{\text{cdm}}h^2$, the posteriors increase towards larger $\Omega_m$ and therefore get more consistent with our results. Furthermore, we use a different sampler compared to A21 and different baryon feedback process modelling. Nevertheless, our results from the $E_n$ analysis are consistent with A21 within $0.03 \sigma$ in $\Omega_m$ and within $0.14 \sigma$ in $S_8$. Similar to A21, we also perform an internal consistency check, removing one $z_{\text{ph}}$-bin at a time and finding consistent results. We discuss this in more detail in Appendix A, and we find that (as expected) the fifth $z_{\text{ph}}$-bin is most important to constrain $\Omega_m$ and $S_8$. We further discuss in Appendix A modelling checks regarding the infusion of IA, baryon feedback and the reduced shear correction. The baryon feedback and reduced shear corrections are always sub-dominant compared to a shift in $S_8$. The IA, however, is important, especially if lower $z_{\text{ph}}$-bins are included and must be accurately modelled.

Finally, we notice that all statistics are consistent with $A_{\text{IA}} = 0$, although both statistics alone seem to favour positive $A_{\text{IA}}$. Interestingly, the joint analysis is shifted to lower $A_{\text{IA}}$ and is slightly more constraining, which we do not observe in the validation in Sect. 8 where all posteriors accurately peak at the input value and all constraining power comes from the $E_n$. This might indicate that for real data with non-zero $A_{\text{IA}}$, third-order shear statistics can contribute (at least a bit) to constraining IA, in line with predictions by, for instance, Pyne & Joachimi (2021) and Troxel & Ishak (2012). Furthermore, the baryonic parameter can be confined to $A_{\text{ba}} < 1.4$ at 68% confidence with our statistic. Here, we should note that $\langle M_{\text{ap}}^3 \rangle$ does not help in constraining $A_{\text{ba}}$, which is probably due to the fact that changes in $S_8$ absorb all $A_{\text{ba}}$ effects.

## 10. Conclusions

This work validates the combined modelling of second- and third-order shear statistics, showing that its accuracy is well-suited for a KiDS-1000 analysis. Our second-order shear statistic of choice is the $E_n$-modes of the COSEBIs (Schneider et al. 2010) and the third-order statistic is $\langle M_{\text{ap}}^3 \rangle$ (Schneider et al. 2005). In particular, we incorporated intrinsic alignment
<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Posterior distribution for data vectors and covariance measured from the T17 galaxy shear catalogues. The covariance matrix and reference data vector were measured from 1944 noisy g mock data.</figcaption></figure>

and uses a slightly larger effective area. Although we rescaled the covariance to correct for the different effective areas, we found in Linke et al. (2023) that this rescaling is not necessarily accurate. Lastly, we notice that our modelling within the KiDS-1000 uncertainty is accurate, which we quantified by measuring the shift of the maximum of the posterior (MP) from the true values in the matter density parameter $\Delta\Omega_m < 0.02 \sigma_{\Omega_m}$ and in the clustering amplitude $\Delta S_8 < 0.05 \sigma_{S_8}$ with respect to the averaged noisy mock data vector and the KiDS-1000 uncertainty. We define the MP as the maximum of the one-dimensional marginal distributions.

## 9. Cosmological results

Finally, we are ready to present the first cosmological constraints from the KiDS-1000 data using second- and third-order shear statistics, displayed in Figs. 7 and 8. The data vectors are described in more detail in A21 for the $E_n$, and the $\langle M_{\text{ap}}^3 \rangle$ in Porth et al. (2023). Given the fact that the covariance matrix is measured only from 1944 realisations using the reduced shear g, we decide to build our data and model vector from all $E_n$ modes and $\langle M_{\text{ap}}^3 \rangle$ with all equal filter radii $\theta_{\text{ap}} \in \{4', 8', 14', 32'\}$. To control whether the model accurately fits the data, we minimised the $\chi^2_{\text{real}}$ from the real data and the $\chi^2_{\text{T17}}$ from each of the 1944 T17 mock data used to compute the covariance matrix. To estimate the probability of measuring a $\chi^2 > \chi^2_{\text{real}}$ (p-value), we counted the number of $\chi^2_{\text{T17}}$ that are greater than $\chi^2_{\text{real}}$ divided by 1944. The resulting p-value for the $E_n$-only, the $\langle M_{\text{ap}}^3 \rangle$-only, and the combination are given in Table 4. The resulting p-values for all combinations are better than 0.1, indicating that our covariance is well matched to the observed data and our model is accurate enough to describe the data. Our maximum posterior $\chi^2$ value increases to 81 if we swap our $E_n$ covariance matrix to the analytical expression in Joachimi et al. (2021) calculated at the MP parameter values in A21. This is unsurprising as the numeri-

cal covariance is computed at the T17 cosmology, giving a signal larger than that computed at the MP of A21. Furthermore, our covariance matrix is measured from reduced shear mock data that slightly increases the covariance matrix. We also notice that the model and the data seem inconsistent for some $z_{\text{ph}}$-bin combinations. However, adjacent COSEBI modes are highly correlated, so visually inspecting the model's goodness of fit to the data is misleading. We refer to A21 for further details on this discrepancy.

We show the resulting posteriors in Fig. 9, where we marginalised over the shift in the redshift distribution and multiplicative shear correction, both stated in Table 1. We improve the constraints on $S_8$ by 23% and on $\Omega_m$ by 47% compared to the $E_n$-only case. This shows how powerful a combined analysis of the second- and third-order shear statistics is. The constraints on $A_{\text{IA}}$ and $A_{\text{ba}}$ are basically untouched, showing that $\langle M_{\text{ap}}^3 \rangle$ is not helpful for constraining these nuisance parameters.

Compared to the maximum of the one-dimensional marginal distributions constraints of $E_n$ measurements given in A21 ($S_8 = 0.758_{-0.026}^{+0.017}$ and $\Omega_m = 0.253_{-0.074}^{+0.088}$), we have slightly larger constraints when using only the $E_n$. This is because we use a numerical covariance, which is larger than the analytical one due to the underlying cosmology and the fact that we model the reduced shear effect. Furthermore, we note that A21 varied $\Omega_{\text{cdm}}h^2$, $\Omega_b h^2$ and h, while we fixed h and $\Omega_b$ and varied only $\Omega_m$. We also find that if we use the same pipeline as A21 but allow larger $\Omega_{\text{cdm}}h^2$, the posteriors increase towards larger $\Omega_m$ and therefore get more consistent with our results. Furthermore, we use a different sampler compared to A21 and different baryon feedback process modelling. Nevertheless, our results from the $E_n$ analysis are consistent with A21 within $0.03 \sigma$ in $\Omega_m$ and within $0.14 \sigma$ in $S_8$. Similar to A21, we also perform an internal consistency check, removing one $z_{\text{ph}}$-bin at a time and finding consistent results. We discuss this in more detail in Appendix A, and we find that (as expected) the fifth $z_{\text{ph}}$-bin is most important to constrain $\Omega_m$ and $S_8$. We further discuss in Appendix A modelling checks regarding the infusion of IA, baryon feedback and the reduced shear correction. The baryon feedback and reduced shear corrections are always sub-dominant compared to a shift in $S_8$. The IA, however, is important, especially if lower $z_{\text{ph}}$-bins are included and must be accurately modelled.

Finally, we notice that all statistics are consistent with $A_{\text{IA}} = 0$, although both statistics alone seem to favour positive $A_{\text{IA}}$. Interestingly, the joint analysis is shifted to lower $A_{\text{IA}}$ and is slightly more constraining, which we do not observe in the validation in Sect. 8 where all posteriors accurately peak at the input value and all constraining power comes from the $E_n$. This might indicate that for real data with non-zero $A_{\text{IA}}$, third-order shear statistics can contribute (at least a bit) to constraining IA, in line with predictions by, for instance, Pyne & Joachimi (2021) and Troxel & Ishak (2012). Furthermore, the baryonic parameter can be confined to $A_{\text{ba}} < 1.4$ at 68% confidence with our statistic. Here, we should note that $\langle M_{\text{ap}}^3 \rangle$ does not help in constraining $A_{\text{ba}}$, which is probably due to the fact that changes in $S_8$ absorb all $A_{\text{ba}}$ effects.

## 10. Conclusions

This work validates the combined modelling of second- and third-order shear statistics, showing that its accuracy is well-suited for a KiDS-1000 analysis. Our second-order shear statistic of choice is the $E_n$-modes of the COSEBIs (Schneider et al. 2010) and the third-order statistic is $\langle M_{\text{ap}}^3 \rangle$ (Schneider et al. 2005). In particular, we incorporated intrinsic alignment
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Measured and modelled $E_n$ for the first five components. The blue points show the measurements from the KiDS-1000 data (A21). The blue error bars indicate the KiDS-1000 uncertainty. The different dashed lines show analytical descriptions at the MP if the combination of $E_n + \langle M_{\text{ap}}^3 \rangle$ or only $E_n$ is used.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Measured and modelled $\langle M_{\text{ap}}^3 \rangle$ with filter radii $\theta_{\text{ap}} \in \{4', 8', 14', 32'\}$. The blue error bars indicate the KiDS-1000 uncertainty. The different dashed lines show analytical descriptions at the MP if the combination of $E_n + \langle M_{\text{ap}}^3 \rangle$ or only $\langle M_{\text{ap}}^3 \rangle$ is used.</figcaption></figure>

<table><caption>**Table 4.** MP values with their marginal 68% credible intervals.</caption><thead><tr><th></th><th>Ω<sub>m</sub></th><th>S<sub>8</sub></th><th>σ<sub>8</sub></th><th>A<sub>IA</sub></th><th>A<sub>ba</sub></th><th>χ<sup>2</sup></th><th>p-value</th></tr></thead><tbody><tr><td>E<sub>n</sub>-only</td><td>0.246<sup>+0.137</sup><sub>-0.076</sub></td><td>0.765<sup>+0.025</sup><sub>-0.032</sub></td><td>0.77<sup>+0.17</sup><sub>-0.16</sub></td><td>0.31<sup>+0.39</sup><sub>-0.46</sub></td><td>0.19<sup>+1.29</sup><sub>-0.18</sub></td><td>72</td><td>0.42</td></tr><tr><td>⟨M<sub>ap</sub><sup>3</sup>⟩-only</td><td>0.248<sup>+0.253</sup><sub>-0.082</sub></td><td>0.791<sup>+0.060</sup><sub>-0.046</sub></td><td>0.61<sup>+0.22</sup><sub>-0.05</sub></td><td>0.49<sup>+0.83</sup><sub>-0.76</sub></td><td>0.39<sup>+2.83</sup><sub>-0.38</sub></td><td>147</td><td>0.25</td></tr><tr><td>E<sub>n</sub> + ⟨M<sub>ap</sub><sup>3</sup>⟩</td><td>0.248<sup>+0.062</sup><sub>-0.055</sub></td><td>0.772 ± 0.022</td><td>0.82<sup>+0.12</sup><sub>-0.09</sub></td><td>0.12<sup>+0.40</sup><sub>-0.46</sub></td><td>0.38<sup>+1.07</sup><sub>-0.37</sub></td><td>224</td><td>0.26</td></tr></tbody></table>

**Notes.** The projected MP and their 68% confidence intervals result from the MCMC chains shown in Fig. 9. The values corresponding to the best $\chi^2$ might differ slightly. We fixed $h = 0.6898$, $w_0 = -1$ and $n_s = 0.969$ but marginalised over the $\delta\langle\zeta\rangle$ and m-bias uncertainties. The p-values were estimated by counting the T17 mock data that give larger $\chi^2$ than the real $\chi^2$ relative to the total number of mock data.
<figure><img src="image_9.png" /><figcaption>**Fig. 7.** Measured and modelled $E_n$ for the first five components. The blue points show the measurements from the KiDS-1000 data (A21). The blue error bars indicate the KiDS-1000 uncertainty. The different dashed lines show analytical descriptions at the MP if the combination of $E_n + \langle M_{\text{ap}}^3 \rangle$ or only $E_n$ is used.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 8.** Measured and modelled $\langle M_{\text{ap}}^3 \rangle$ with filter radii $\theta_{\text{ap}} \in \{4', 8', 14', 32'\}$. The blue error bars indicate the KiDS-1000 uncertainty. The different dashed lines show analytical descriptions at the MP if the combination of $E_n + \langle M_{\text{ap}}^3 \rangle$ or only $\langle M_{\text{ap}}^3 \rangle$ is used.</figcaption></figure>

<table><caption>**Table 4.** MP values with their marginal 68% credible intervals.</caption><thead><tr><th></th><th>Ω<sub>m</sub></th><th>S<sub>8</sub></th><th>σ<sub>8</sub></th><th>A<sub>IA</sub></th><th>A<sub>ba</sub></th><th>χ<sup>2</sup></th><th>p-value</th></tr></thead><tbody><tr><td>E<sub>n</sub>-only</td><td>0.246<sup>+0.137</sup><sub>-0.076</sub></td><td>0.765<sup>+0.025</sup><sub>-0.032</sub></td><td>0.77<sup>+0.17</sup><sub>-0.16</sub></td><td>0.31<sup>+0.39</sup><sub>-0.46</sub></td><td>0.19<sup>+1.29</sup><sub>-0.18</sub></td><td>72</td><td>0.42</td></tr><tr><td>⟨M<sup>3</sup><sub>ap</sub>⟩-only</td><td>0.248<sup>+0.253</sup><sub>-0.082</sub></td><td>0.791<sup>+0.060</sup><sub>-0.046</sub></td><td>0.61<sup>+0.22</sup><sub>-0.05</sub></td><td>0.49<sup>+0.83</sup><sub>-0.76</sub></td><td>0.39<sup>+2.83</sup><sub>-0.38</sub></td><td>147</td><td>0.25</td></tr><tr><td>E<sub>n</sub> + ⟨M<sup>3</sup><sub>ap</sub>⟩</td><td>0.248<sup>+0.062</sup><sub>-0.055</sub></td><td>0.772 ± 0.022</td><td>0.82<sup>+0.12</sup><sub>-0.09</sub></td><td>0.12<sup>+0.40</sup><sub>-0.46</sub></td><td>0.38<sup>+1.07</sup><sub>-0.37</sub></td><td>224</td><td>0.26</td></tr></tbody></table>

**Notes.** The projected MP and their 68% confidence intervals result from the MCMC chains shown in Fig. 9. The values corresponding to the best $\chi^2$ might differ slightly. We fixed $h = 0.6898$, $w_0 = -1$ and $n_s = 0.969$ but marginalised over the $\delta\langle\zeta\rangle$ and m-bias uncertainties. The p-values were estimated by counting the T17 mock data that give larger $\chi^2$ than the real $\chi^2$ relative to the total number of mock data.
<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Posterior distribution for the real KiDS-1000 data vector while the covariance is measured from the 1944 T17 galaxy shear catalogues. Here the $E_n$ are compared to $\langle M_{\text{ap}}^3 \rangle$ using all available combinations for the aperture filter radii $\theta_{\text{ap}} \in \{4', 8', 14', 32'\}$ and redshift bin combinations. We note that the prior (see Table 2) range on $\Omega_m$ is enlarged compared to the previous validation plots. The FoM values of the $E_n$-only or $\langle M_{\text{ap}}^3 \rangle$-only case should not be compared to Fig. 6, as priors of $\Omega_m$ bind their contours.</figcaption></figure>

modelling based on the non-linear alignment model of Bridle & King (2007) and validated its accuracy against simulations infused with IA effects. This test is also interesting for other simulation-based analyses for which IA cannot be modelled analytically. We incorporated the impact of the baryonic feedback process by measuring a response function using the Magneticum simulations. Since the amplitude of this response function has no physical meaning, it is considered a nuisance parameter, which does not bias our cosmological parameter predictions.

We investigate which parts of the data vector can be neglected without losing too much cosmological information.

This is important because the data vector for third-order shear statistics, due to the possibility it offers of combining three different filters with three different redshift bins, inflating the data vector to several hundred elements easily. Therefore, both a numerical and an analytical covariance are difficult to compute. We find that cross-tomographic redshift bins contain a large amount of cosmological information. Using only equal-scale filter radii but all available tomographic bin combinations was the best data compression strategy, which comes with the advantage that equal-scale filters are faster to compute analytically as they require a lower integral accuracy. Next, we investigated the chosen filter radii. For this we measured and modelled the $\langle M_{\text{ap}}^3 \rangle$
<figure><img src="image_11.png" /><figcaption>**Fig. 9.** Posterior distribution for the real KiDS-1000 data vector while the covariance is measured from the 1944 T17 galaxy shear catalogues. Here the $E_n$ are compared to $\langle \mathcal{M}_{\text{ap}}^3 \rangle$ using all available combinations for the aperture filter radii $\theta_{\text{ap}} \in \{4', 8', 14', 32'\}$ and redshift bin combinations. We note that the prior (see Table 2) range on $\Omega_m$ is enlarged compared to the previous validation plots. The FoM values of the $E_n$-only or $\langle \mathcal{M}_{\text{ap}}^3 \rangle$-only case should not be compared to Fig. 6, as priors of $\Omega_m$ bind their contours.</figcaption></figure>

modelling based on the non-linear alignment model of Bridle & King (2007) and validated its accuracy against simulations infused with IA effects. This test is also interesting for other simulation-based analyses for which IA cannot be modelled analytically. We incorporated the impact of the baryonic feedback process by measuring a response function using the Magneticum simulations. Since the amplitude of this response function has no physical meaning, it is considered a nuisance parameter, which does not bias our cosmological parameter predictions.

We investigate which parts of the data vector can be neglected without losing too much cosmological information.

This is important because the data vector for third-order shear statistics, due to the possibility it offers of combining three different filters with three different redshift bins, inflating the data vector to several hundred elements easily. Therefore, both a numerical and an analytical covariance are difficult to compute. We find that cross-tomographic redshift bins contain a large amount of cosmological information. Using only equal-scale filter radii but all available tomographic bin combinations was the best data compression strategy, which comes with the advantage that equal-scale filters are faster to compute analytically as they require a lower integral accuracy. Next, we investigated the chosen filter radii. For this we measured and modelled the $\langle \mathcal{M}_{\text{ap}}^3 \rangle$.
for $\theta_{ap} \in \{4', 6', 8', 10', 14', 18', 22', 26', 32', 36'\}$. We find that filter radii above 30' do not contribute much constraining power and that having many small filter radii is unnecessary. The best option is to use filter radii $\theta_{ap} \in \{4', 8', 14', 32'\}$. We also tested whether selecting elements that maximise the Fisher information matrix on $S_8$ is more optimal but find no difference compared to selecting all equal-scale filter radii. However, this method can be used to speed up the modelling and measurement of future analyses by discarding irrelevant elements.

Next, we validated the assumption that the correlation function estimator gives accurate results. This is important since only correlation function estimators have the potential to give unbiased results for real data. We used realistic mock data created from the T17 simulations. In particular, we created several galaxy catalogues where the positions of the galaxies are exactly at the KiDS-1000 galaxy positions. We had to rely on correlation functions to measure the second- and third-order statistics, which give unbiased results also if the data has a complex topology. Our first finding is that our modelling and measurement result in unbiased cosmological parameters given the KiDS-1000 uncertainty. Second, we find that using the reduced shear, or the shear itself, does not change the results and can, therefore, be ignored for the KiDS-1000 data analysis.

We conclude this paper with an analysis of the real KiDS-1000 data. Overall, our $E_n$ constraints are less informative than the original KiDS-1000 analysis. This is mostly because we used a numerical covariance matrix from T17 simulations. However, the chosen sampler, the priors on the cosmological parameters, and the modelling strategy of the baryonic feedback processes impact our constraints, too. We find an $S_8 = 0.772 \pm 0.022$ and an $\Omega_m = 0.248_{-0.055}^{+0.062}$, which are improved compared to the $E_n$-only case by 23% and by 47%, respectively. With a $p$-value of 0.25, we also find a good agreement of model and data given the KiDS-1000 uncertainty. This demonstrates that combining second- and third-order statistics is powerful in constraining cosmological parameters. The gain in constraining power in $\Omega_m$ is also interesting for combined weak lensing and galaxy clustering analysis because the constraining power in $\Omega_m$ for clustering analysis comes with the issue of further nuisance parameters such as galaxy bias. However, since second- and third-order shear statistics constrain $\Omega_m$ quite well, combining it with clustering statistics might enable us to learn more about these nuisance parameters.

We leave the optimisation of the IA and baryonic feedback modelling for future analysis. An especially interesting improvement would be a more physically motivated description of the baryon feedback processes, which are identical for power and bispectrum. As all baryon feedback models rely on hydro simulations, we have to use the same simulations for power and bispectrum. Furthermore, although we found that the reduced shear and limber approximation is sufficient for a KiDS-1000 analysis, we probably have to model these effects for future Stage IV surveys. Lastly, we ignore the effect of source clustering for this work. Although Gatti et al. (2024) found it relevant for third-order weak lensing statistics based on convergence mass maps, we expect it to be less critical for our analysis, which uses shear catalogues and no mass map reconstruction. For future Stage IV surveys, this needs to be investigated.

**Acknowledgements.** This paper went through the KiDS review process, and we want to thank the KiDS internal reviewer for the fruitful comments to improve this paper. We thank Mike Jarvis for maintaining treecorr and his constructive and fruitful comments. Furthermore, we thank Alessio Mancini for developing the CosmoPower emulator, which made the analysis pipeline of this work significantly faster. Some of the results in this paper have been derived using the

HEALPY and HEALIX package (currently http://healpix.sourceforge.net; Górski et al. 2005; Zonca et al. 2019). The figures in this work were created with MATPLOTLIB (Hunter 2007) and CHAINCONSUMER (Hinton 2016). We further make use of COSMOSIS (Zuntz et al. 2015a), NUMPY (Harris et al. 2020) and SCIPY (Jones et al. 2001) software packages. L.P. acknowledges support from the DLR grant 50QE2002. SH is supported by the US Department of Energy, Office of Science, Office of High Energy Physics under Award Number DE-SC0019301. We acknowledge use of the lux supercomputer at UC Santa Cruz, funded by NSF MRI grant AST 1828315. L.L. is supported by the Austrian Science Fund (FWF) [ESP 357-NJ]. T.C. is supported by the INFN INDARK PD51 grant and by the FARE MIUR grant 'ClustersXEuclid' R165SBKTMA. K.D. acknowledges support by the COMPLEX project from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program grant agreement ERC-2019-AdG 882679 as well as support by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy - EXC-2094 - 390783311. J.H.D. acknowledges support from an STFC Ernest Rutherford Fellowship (project reference ST/S004858/1). H.H. 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). K.K. acknowledges support from the Royal Society and Imperial College, N.M. acknowledges support from the Centre National d'Etudes Spatiales (CNES) fellowship. Author contributions: All authors contributed to the development and writing of this paper. The authorship list is given in three groups: the lead authors (P.A.B., L.P., S.H., L.L., N.W., P.S.), followed by an alphabetical group. This alphabetical group includes those who have either made a significant contribution to the data products or to the scientific analysis. The KiDS results in this paper are based on observations made with ESO Telescopes at the La Silla Paranal Observatory under programme IDs 177.A-3016, 177.A-3017, 177.A-3018 and 179.A-2004, and on data products produced by the KiDS consortium. The KiDS production team acknowledges support from: Deutsche Forschungsgemeinschaft, ERC, NOVA and NWO-M grants; Target; the University of Padova, and the University Federico II (Naples). Data processing for VIKING has been contributed by the VISTA Data Flow System at CASU, Cambridge and WFAU, Edinburgh. The calculations for the hydrodynamical simulations were carried out at the Leibniz Supercomputer Center (LRZ) under the project pr83li (Magneticum). J.H.D. acknowledges support from an STFC Ernest Rutherford Fellowship (project reference ST/S004858/1).

## References

Asgari, M., Heymans, C., Hildebrandt, H., et al. 2019, A&A, 624, A134

Asgari, M., Tröster, T., Heymans, C., et al. 2020, A&A, 634, A127

Asgari, M., Lin, C.-A., Joachimi, B., et al. 2021, A&A, 645, A104

Bartelmann, M. 2010, Class. Quant. Grav., 27, 233001

Bartelmann, M., & Schneider, P. 2001, Phys. Rep., 340, 291

Beck, A. M., Murante, G., Arth, A., et al. 2016, MNRAS, 455, 2110

Bergé, J., Amara, A., & Réfrégier, A. 2010, ApJ, 712, 992

Bernardeau, F., van Waaerbeke, L., & Mellier, Y. 1997, A&A, 322, 1

Bridle, S., & King, L. 2007, New J. Phys., 9, 444

Brown, M. L., Taylor, A. N., Hambly, N. C., & Dye, S. 2002, MNRAS, 333, 501

Burger, P. A., Friedrich, O., Harnois-Déraps, J., et al. 2023, A&A, 669, A69

Castro, T., Quartin, M., Giocoli, C., Borgani, S., & Dolag, K. 2018, MNRAS, 478, 1305

Chisari, N., Codis, S., Laigle, C., et al. 2015, MNRAS, 454, 2736

Crittenden, R. G., Natarajan, P., Pen, U.-L., & Theuns, T. 2002, ApJ, 568, 20

Dalal, R., Li, X., Nicola, A., et al. 2023, arXiv e-prints [arXiv:2304.00701]

de Jong, J. T. A., Verdoes Kleijn, G. A., Boxhoorn, D. R., et al. 2015, A&A, 582, A62

de Jong, J. T. A., Verdoes Kleijn, G. A., Erben, T., et al. 2017, A&A, 604, A134

DES Collaboration (Abbott, T. M. C., et al.) 2022, Phys. Rev. D, 105, 023520

DES & KiDS Collaboration (Abbott, T. M. C., et al.) 2023, Open J. Astrophys., 6, 36

Di Valentino, E., Anchordoqui, L. A., Akarsu, Ö., et al. 2021, Astropart. Phys., 131, 102604

Edge, A., Sutherland, W., Kuijken, K., et al. 2013, The Messenger, 154, 32

Falco, E. E., Gorenstein, M. V., & Shapiro, I. I. 1985, ApJ, 289, L1

Fenech Conti, I., Herbonnet, R., Hoekstra, H., et al. 2017, MNRAS, 467, 1627

Fornberg, B. 1988, Math. Comp., 184, 699

Fu, L., Kilbinger, M., Erben, T., et al. 2014, MNRAS, 441, 2725

Gatti, M., Jain, B., Chang, C., et al. 2022, Phys. Rev. D, 106, 083509P

Gatti, M., Jeffrey, N., Whiteway, L., et al. 2024, MNRAS, 527, L115

Giblin, B., Heymans, C., Asgari, M., et al. 2021, A&A, 645, A105

Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759

Gruen, D., Friedrich, O., Krause, E., et al. 2018, Phys. Rev. D, 98, 023507

Halder, A., & Barreira, A. 2022, MNRAS, 515, 4639

Halder, A., Friedrich, O., Seitz, S., & Varga, T. N. 2021, MNRAS, 506, 2780

Harnois-Déraps, J., Amon, A., Choi, A., et al. 2018, MNRAS, 481, 1337
for $\theta_{ap} \in \{4', 6', 8', 10', 14', 18', 22', 26', 32', 36'\}$. We find that filter radii above 30' do not contribute much constraining power and that having many small filter radii is unnecessary. The best option is to use filter radii $\theta_{ap} \in \{4', 8', 14', 32'\}$. We also tested whether selecting elements that maximise the Fisher information matrix on $S_8$ is more optimal but find no difference compared to selecting all equal-scale filter radii. However, this method can be used to speed up the modelling and measurement of future analyses by discarding irrelevant elements.

Next, we validated the assumption that the correlation function estimator gives accurate results. This is important since only correlation function estimators have the potential to give unbiased results for real data. We used realistic mock data created from the T17 simulations. In particular, we created several galaxy catalogues where the positions of the galaxies are exactly at the KiDS-1000 galaxy positions. We had to rely on correlation functions to measure the second- and third-order statistics, which give unbiased results also if the data has a complex topology. Our first finding is that our modelling and measurement result in unbiased cosmological parameters given the KiDS-1000 uncertainty. Second, we find that using the reduced shear, or the shear itself, does not change the results and can, therefore, be ignored for the KiDS-1000 data analysis.

We conclude this paper with an analysis of the real KiDS-1000 data. Overall, our $E_n$ constraints are less informative than the original KiDS-1000 analysis. This is mostly because we used a numerical covariance matrix from T17 simulations. However, the chosen sampler, the priors on the cosmological parameters, and the modelling strategy of the baryonic feedback processes impact our constraints, too. We find an $S_8 = 0.772 \pm 0.022$ and an $\Omega_m = 0.248_{-0.055}^{+0.062}$, which are improved compared to the $E_n$-only case by 23% and by 47%, respectively. With a $p$-value of 0.25, we also find a good agreement of model and data given the KiDS-1000 uncertainty. This demonstrates that combining second- and third-order statistics is powerful in constraining cosmological parameters. The gain in constraining power in $\Omega_m$ is also interesting for combined weak lensing and galaxy clustering analysis because the constraining power in $\Omega_m$ for clustering analysis comes with the issue of further nuisance parameters such as galaxy bias. However, since second- and third-order shear statistics constrain $\Omega_m$ quite well, combining it with clustering statistics might enable us to learn more about these nuisance parameters.

We leave the optimisation of the IA and baryonic feedback modelling for future analysis. An especially interesting improvement would be a more physically motivated description of the baryon feedback processes, which are identical for power and bispectrum. As all baryon feedback models rely on hydro simulations, we have to use the same simulations for power and bispectrum. Furthermore, although we found that the reduced shear and limber approximation is sufficient for a KiDS-1000 analysis, we probably have to model these effects for future Stage IV surveys. Lastly, we ignore the effect of source clustering for this work. Although Gatti et al. (2024) found it relevant for third-order weak lensing statistics based on convergence mass maps, we expect it to be less critical for our analysis, which uses shear catalogues and no mass map reconstruction. For future Stage IV surveys, this needs to be investigated.

**Acknowledgements.** This paper went through the KiDS review process, and we want to thank the KiDS internal reviewer for the fruitful comments to improve this paper. We thank Mike Jarvis for maintaining treecorr and his constructive and fruitful comments. Furthermore, we thank Alessio Mancini for developing the CosmoPower emulator, which made the analysis pipeline of this work significantly faster. Some of the results in this paper have been derived using the

HEALPY and HEALIX package (currently http://healpix.sourceforge.net; Górski et al. 2005; Zonca et al. 2019). The figures in this work were created with MATPLOTLIB (Hunter 2007) and CHAINCONSUMER (Hinton 2016). We further make use of COSMOSIS (Zuntz et al. 2015a), NUMPY (Harris et al. 2020) and SCIPY (Jones et al. 2001) software packages. L.P. acknowledges support from the DLR grant 50QE2002. SH is supported by the US Department of Energy, Office of Science, Office of High Energy Physics under Award Number DE-SC0019301. We acknowledge use of the lux supercomputer at UC Santa Cruz, funded by NSF MRI grant AST 1828315. L.L. is supported by the Austrian Science Fund (FWF) [ESP 357-NJ]. T.C. is supported by the INFN INDARK PD51 grant and by the FARE MIUR grant 'ClustersXEuclid' R165SBKTMA. K.D. acknowledges support by the COMPLEX project from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program grant agreement ERC-2019-AdG 882679 as well as support by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy - EXC-2094 - 390783311. J.H.D. acknowledges support from an STFC Ernest Rutherford Fellowship (project reference ST/S004858/1). H.H. 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). K.K. acknowledges support from the Royal Society and Imperial College, N.M. acknowledges support from the Centre National d'Etudes Spatiales (CNES) fellowship. Author contributions: All authors contributed to the development and writing of this paper. The authorship list is given in three groups: the lead authors (P.A.B., L.P., S.H., L.L., N.W., P.S.), followed by an alphabetical group. This alphabetical group includes those who have either made a significant contribution to the data products or to the scientific analysis. The KiDS results in this paper are based on observations made with ESO Telescopes at the La Silla Paranal Observatory under programme IDs 177.A-3016, 177.A-3017, 177.A-3018 and 179.A-2004, and on data products produced by the KiDS consortium. The KiDS production team acknowledges support from: Deutsche Forschungsgemeinschaft, ERC, NOVA and NWO-M grants; Target; the University of Padova, and the University Federico II (Naples). Data processing for VIKING has been contributed by the VISTA Data Flow System at CASU, Cambridge and WFAU, Edinburgh. The calculations for the hydrodynamical simulations were carried out at the Leibniz Supercomputer Center (LRZ) under the project pr83li (Magneticum). J.H.D. acknowledges support from an STFC Ernest Rutherford Fellowship (project reference ST/S004858/1).

## References

* Asgari, M., Heymans, C., Hildebrandt, H., et al. 2019, A&A, 624, A134

* Asgari, M., Tröster, T., Heymans, C., et al. 2020, A&A, 634, A127

* Asgari, M., Lin, C.-A., Joachimi, B., et al. 2021, A&A, 645, A104

* Bartelmann, M. 2010, Class. Quant. Grav., 27, 233001

* Bartelmann, M., & Schneider, P. 2001, Phys. Rep., 340, 291

* Beck, A. M., Murante, G., Arth, A., et al. 2016, MNRAS, 455, 2110

* Bergé, J., Amara, A., & Réfrégier, A. 2010, ApJ, 712, 992

* Bernardeau, F., van Waerbeke, L., & Mellier, Y. 1997, A&A, 322, 1

* Bridle, S., & King, L. 2007, New J. Phys., 9, 444

* Brown, M. L., Taylor, A. N., Hambly, N. C., & Dye, S. 2002, MNRAS, 333, 501

* Burger, P. A., Friedrich, O., Harnois-Déraps, J., et al. 2023, A&A, 669, A69

* Castro, T., Quartin, M., Giocoli, C., Borgani, S., & Dolag, K. 2018, MNRAS, 478, 1305

* Chisari, N., Codis, S., Laigle, C., et al. 2015, MNRAS, 454, 2736

* Crittenden, R. G., Natarajan, P., Pen, U.-L., & Theuns, T. 2002, ApJ, 568, 20

* Dalal, R., Li, X., Nicola, A., et al. 2023, arXiv e-prints [arXiv:2304.00701]

* de Jong, J. T. A., Verdoes Kleijn, G. A., Boxhoorn, D. R., et al. 2015, A&A, 582, A62

* de Jong, J. T. A., Verdoes Kleijn, G. A., Erben, T., et al. 2017, A&A, 604, A134

* DES Collaboration (Abbott, T. M. C., et al.) 2022, Phys. Rev. D, 105, 023520

* DES & KiDS Collaboration (Abbott, T. M. C., et al.) 2023, Open J. Astrophys., 6, 36

* Di Valentino, E., Anchordoqui, L. A., Akarsu, Ö., et al. 2021, Astropart. Phys., 131, 102604

* Edge, A., Sutherland, W., Kuijken, K., et al. 2013, The Messenger, 154, 32

* Falco, E. E., Gorenstein, M. V., & Shapiro, I. I. 1985, ApJ, 289, L1

* Fenech Conti, I., Herbonnet, R., Hoekstra, H., et al. 2017, MNRAS, 467, 1627

* Fornberg, B. 1988, Math. Comp., 184, 699

* Fu, L., Kilbinger, M., Erben, T., et al. 2014, MNRAS, 441, 2725

* Gatti, M., Jain, B., Chang, C., et al. 2022, Phys. Rev. D, 106, 083509P

* Gatti, M., Jeffrey, N., Whiteway, L., et al. 2024, MNRAS, 527, L115

* Giblin, B., Heymans, C., Asgari, M., et al. 2021, A&A, 645, A105

* Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759

* Gruen, D., Friedrich, O., Krause, E., et al. 2018, Phys. Rev. D, 98, 023507

* Halder, A., & Barreira, A. 2022, MNRAS, 515, 4639

* Halder, A., Friedrich, O., Seitz, S., & Varga, T. N. 2021, MNRAS, 506, 2780

* Harnois-Déraps, J., Amon, A., Choi, A., et al. 2018, MNRAS, 481, 1337
Harnois-Déraps, J., Giblin, B., & Joachimi, B. 2019, A&A, 631, A160

Harnois-Déraps, J., Martinet, N., Castro, T., et al. 2021, MNRAS, 506, 1623

Harnois-Déraps, J., Martinet, N., & Reischke, R. 2022, MNRAS, 509, 3868

Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357

Heydenreich, S., Brück, B., & Harnois-Déraps, J. 2021, A&A, 648, A74

Heydenreich, S., Linke, L., Burger, P., & Schneider, P. 2023, A&A, 672, A44

Heymans, C., Tröster, T., Asgari, M., et al. 2021, A&A, 646, A140

Hikage, C., Oguri, M., Hamana, T., et al. 2019, PASJ, 71, 43

Hildebrandt, H., Viola, M., Heymans, C., et al. 2017, MNRAS, 465, 1454

Hildebrandt, H., Köhlinger, F., van den Busch, J. L., et al. 2020, A&A, 633, A69

Hildebrandt, H., van den Busch, J. L., Wright, A. H., et al. 2021, A&A, 647, A124

Hinton, S. R. 2016, J. Open Source Softw., 1, 00045

Hirschmann, M., Dolag, K., Saro, A., et al. 2014, MNRAS, 442, 2304

Hoekstra, H., & Jain, B. 2008, Ann. Rev. Nucl. Part. Sci., 58, 99

Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90

Ivezic, Z., Axelrod, T., Brandt, W. N., et al. 2008, Serb. Astron. J., 176, 1

Jarvis, M., Bernstein, G., & Jain, B. 2004, MNRAS, 352, 338

Joachimi, B., Cacciato, M., Kitching, T. D., et al. 2015, Space. Sci. Rev., 193, 1

Joachimi, B., Lin, C. A., Asgari, M., et al. 2021, A&A, 646, A129

Jones, E., Oliphant, T., Peterson, P., et al. 2001, SciPy: Open Source Scientific Tools for Python, http://www.scipy.org/

Joudaki, S., Hildebrandt, H., Traykova, D., et al. 2020, A&A, 638, L1

Kaiser, N., & Jaffe, A. 1997, ApJ, 484, 545

Kannawadi, A., Hoekstra, H., Miller, L., et al. 2019, A&A, 624, A92

Kilbinger, M. 2015, Rep. Progr. Phys., 78, 086901

Kilbinger, M., & Schneider, P. 2005, A&A, 442, 69

Kilbinger, M., Heymans, C., Asgari, M., et al. 2017, MNRAS, 472, 2126

Kuijken, K., Heymans, C., Hildebrandt, H., et al. 2015, MNRAS, 454, 3500

Kuijken, K., Heymans, C., Dvornik, A., et al. 2019, A&A, 625, A2

Laureijs, R., Amiaux, J., Arduini, S., et al. 2011, arXiv e-prints [arXiv:1110.3193]

Li, S. S., Hoekstra, H., Kuijken, K., et al. 2023, A&A, 679, A133

Limber, D. N. 1954, ApJ, 119, 655

Linke, L., Heydenreich, S., Burger, P. A., & Schneider, P. 2023, A&A, 672, A185

LoVerde, M., & Afshordi, N. 2008, Phys. Rev. D, 78, 123506

Martinet, N., Schneider, P., Hildebrandt, H., et al. 2018, MNRAS, 474, 712

Martinet, N., Castro, T., Harnois-Déraps, J., et al. 2021, A&A, 648, A115

Mead, A. J., Brieden, S., Tröster, T., & Heymans, C. 2021, MNRAS, 502, 1401

Miller, L., Heymans, C., Kitching, T. D., et al. 2013, MNRAS, 429, 2858

Munshi, D., Valageas, P., van Waerbeke, L., & Heavens, A. 2008, Phys. Rep., 462, 67

Percival, W. J., Friedrich, O., Sellentin, E., & Heavens, A. 2022, MNRAS, 510, 3207

Pires, S., Leonard, A., & Starck, J.-L. 2012, MNRAS, 423, 983

Planck Collaboration VI. 2020, A&A, 641, A6

Porth, L., Heydenreich, S., Burger, P., Linke, L., & Schneider, P. 2023, arXiv e-prints [arXiv:2309.08601]

Pyne, S., & Joachimi, B. 2021, MNRAS, 503, 2300

Ragagnin, A., Dolag, K., Biffi, V., et al. 2017, Astron. Comput., 20, 52

Sánchez, A. G., Scoccimarro, R., Crocce, M., et al. 2017, MNRAS, 464, 1640

Schneider, P. 1996, MNRAS, 283, 837

Schneider, P., & Lombardi, M. 2003, A&A, 397, 809

Schneider, P., & Seitz, C. 1995, A&A, 294, 411

Schneider, P., van Waerbeke, L., Jain, B., & Kruse, G. 1998, MNRAS, 296, 873

Schneider, P., van Waerbeke, L., & Mellier, Y. 2002, A&A, 389, 729

Schneider, P., Kilbinger, M., & Lombardi, M. 2005, A&A, 431, 9

Schneider, P., Eifler, T., & Krause, E. 2010, A&A, 520, A116

Secco, L. F., Jarvis, M., Jain, B., et al. 2022, Phys. Rev. D, 105, 103537

Seitz, C., & Schneider, P. 1997, A&A, 318, 687

Springel, V. 2005, MNRAS, 364, 1105

Spurio Mancini, A., Piras, D., Alsing, J., Joachimi, B., & Hobson, M. P. 2022, MNRAS, 511, 1771

Takahashi, R., Hamana, T., Shirasaki, M., et al. 2017, ApJ, 850, 24

Takahashi, R., Nishimichi, T., Namikawa, T., et al. 2020, ApJ, 895, 113

Tegmark, M., Taylor, A. N., & Heavens, A. F. 1997, ApJ, 480, 22

Troxel, M. A., & Ishak, M. 2012, MNRAS, 419, 1804

Troxel, M. A., MacCrann, N., Zuntz, J., et al. 2018, Phys. Rev. D, 98, 043528

van Daalen, M. P., Schaye, J., Booth, C. M., & Dalla Vecchia, C. 2011, MNRAS, 415, 3649

van den Busch, J. L., Wright, A. H., Hildebrandt, H., et al. 2022, A&A, 664, A170

Wright, A. H., Hildebrandt, H., van den Busch, J. L., et al. 2020, A&A, 640, L14

Zonca, A., Singer, L., Lenz, D., et al. 2019, J. Open Source Softw., 4, 1298

Zuntz, J., Paterno, M., Jennings, E., et al. 2015, Astron. Comput., 12, 45
Harnois-Déraps, J., Giblin, B., & Joachimi, B. 2019, A&A, 631, A160

Harnois-Déraps, J., Martinet, N., Castro, T., et al. 2021, MNRAS, 506, 1623

Harnois-Déraps, J., Martinet, N., & Reischke, R. 2022, MNRAS, 509, 3868

Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357

Heydenreich, S., Brück, B., & Harnois-Déraps, J. 2021, A&A, 648, A74

Heydenreich, S., Linke, L., Burger, P., & Schneider, P. 2023, A&A, 672, A44

Heymans, C., Tröster, T., Asgari, M., et al. 2021, A&A, 646, A140

Hikage, C., Oguri, M., Hamana, T., et al. 2019, PASJ, 71, 43

Hildebrandt, H., Viola, M., Heymans, C., et al. 2017, MNRAS, 465, 1454

Hildebrandt, H., Köhlinger, F., van den Busch, J. L., et al. 2020, A&A, 633, A69

Hildebrandt, H., van den Busch, J. L., Wright, A. H., et al. 2021, A&A, 647, A124

Hinton, S. R. 2016, J. Open Source Softw., 1, 00045

Hirschmann, M., Dolag, K., Saro, A., et al. 2014, MNRAS, 442, 2304

Hoekstra, H., & Jain, B. 2008, Ann. Rev. Nucl. Part. Sci., 58, 99

Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90

Ivezic, Z., Axelrod, T., Brandt, W. N., et al. 2008, Serb. Astron. J., 176, 1

Jarvis, M., Bernstein, G., & Jain, B. 2004, MNRAS, 352, 338

Joachimi, B., Cacciato, M., Kitching, T. D., et al. 2015, Space. Sci. Rev., 193, 1

Joachimi, B., Lin, C. A., Asgari, M., et al. 2021, A&A, 646, A129

Jones, E., Oliphant, T., Peterson, P., et al. 2001, SciPy: Open Source Scientific Tools for Python, http://www.scipy.org/

Joudaki, S., Hildebrandt, H., Traykova, D., et al. 2020, A&A, 638, L1

Kaiser, N., & Jaffe, A. 1997, ApJ, 484, 545

Kannawadi, A., Hoekstra, H., Miller, L., et al. 2019, A&A, 624, A92

Kilbinger, M. 2015, Rep. Progr. Phys., 78, 086901

Kilbinger, M., & Schneider, P. 2005, A&A, 442, 69

Kilbinger, M., Heymans, C., Asgari, M., et al. 2017, MNRAS, 472, 2126

Kuijken, K., Heymans, C., Hildebrandt, H., et al. 2015, MNRAS, 454, 3500

Kuijken, K., Heymans, C., Dvornik, A., et al. 2019, A&A, 625, A2

Laureijs, R., Amiaux, J., Arduini, S., et al. 2011, arXiv e-prints
[arXiv:1110.3193]

Li, S. S., Hoekstra, H., Kuijken, K., et al. 2023, A&A, 679, A133

Limber, D. N. 1954, ApJ, 119, 655

Linke, L., Heydenreich, S., Burger, P. A., & Schneider, P. 2023, A&A, 672, A185

LoVerde, M., & Afshordi, N. 2008, Phys. Rev. D, 78, 123506

Martinet, N., Schneider, P., Hildebrandt, H., et al. 2018, MNRAS, 474, 712

Martinet, N., Castro, T., Harnois-Déraps, J., et al. 2021, A&A, 648, A115

Mead, A. J., Brieden, S., Tröster, T., & Heymans, C. 2021, MNRAS, 502, 1401

Miller, L., Heymans, C., Kitching, T. D., et al. 2013, MNRAS, 429, 2858

Munshi, D., Valageas, P., van Waerbeke, L., & Heavens, A. 2008, Phys. Rep., 462, 67

Percival, W. J., Friedrich, O., Sellentin, E., & Heavens, A. 2022, MNRAS, 510, 3207

Pires, S., Leonard, A., & Starck, J.-L. 2012, MNRAS, 423, 983

Planck Collaboration VI. 2020, A&A, 641, A6

Porth, L., Heydenreich, S., Burger, P., Linke, L., & Schneider, P. 2023, arXiv e-prints [arXiv:2309.08601]

Pyne, S., & Joachimi, B. 2021, MNRAS, 503, 2300

Ragagnin, A., Dolag, K., Biffi, V., et al. 2017, Astron. Comput., 20, 52

Sánchez, A. G., Scoccimarro, R., Crocce, M., et al. 2017, MNRAS, 464, 1640

Schneider, P. 1996, MNRAS, 283, 837

Schneider, P., & Lombardi, M. 2003, A&A, 397, 809

Schneider, P., & Seitz, C. 1995, A&A, 294, 411

Schneider, P., van Waerbeke, L., Jain, B., & Kruse, G. 1998, MNRAS, 296, 873

Schneider, P., van Waerbeke, L., & Mellier, Y. 2002, A&A, 389, 729

Schneider, P., Kilbinger, M., & Lombardi, M. 2005, A&A, 431, 9

Schneider, P., Eifler, T., & Krause, E. 2010, A&A, 520, A116

Secco, L. F., Jarvis, M., Jain, B., et al. 2022, Phys. Rev. D, 105, 103537

Seitz, C., & Schneider, P. 1997, A&A, 318, 687

Springel, V. 2005, MNRAS, 364, 1105

Spurio Mancini, A., Piras, D., Alsing, J., Joachimi, B., & Hobson, M. P. 2022, MNRAS, 511, 1771

Takahashi, R., Hamana, T., Shirasaki, M., et al. 2017, ApJ, 850, 24

Takahashi, R., Nishimichi, T., Namikawa, T., et al. 2020, ApJ, 895, 113

Tegmark, M., Taylor, A. N., & Heavens, A. F. 1997, ApJ, 480, 22

Troxel, M. A., & Ishak, M. 2012, MNRAS, 419, 1804

Troxel, M. A., MacCrann, N., Zuntz, J., et al. 2018, Phys. Rev. D, 98, 043528

van Daalen, M. P., Schaye, J., Booth, C. M., & Dalla Vecchia, C. 2011, MNRAS, 415, 3649

van den Busch, J. L., Wright, A. H., Hildebrandt, H., et al. 2022, A&A, 664, A170

Wright, A. H., Hildebrandt, H., van den Busch, J. L., et al. 2020, A&A, 640, L14

Zonca, A., Singer, L., Lenz, D., et al. 2019, J. Open Source Softw., 4, 1298

Zuntz, J., Paterno, M., Jennings, E., et al. 2015, Astron. Comput., 12, 45
## Appendix A: Consistency checks and modelling choices

As a first consistency check, we again inferred the joint analysis posteriors but removed one $z_{ph}$-bin at a time. The results are shown in Fig. A.1 and reveal (most importantly) that all $z_{ph}$-bins are internally consistent with each other. However, we further observe that $z_{ph}$-bin one and $z_{ph}$-bin two are not important for inferring cosmological parameters, as expected, given their S/N as shown in Porth et al. (2023). For the $A_{IA}$ parameter, in turn, $z_{ph}$-bin one is very important, as the low redshift is most sensitive to IA. Next, we observe that the third $z_{ph}$-bin results are basically the same as for all five $z_{ph}$-bins, meaning that the third bin can eventually be discarded in future KiDS analyses if the dimension needs to be decreased further. Lastly, if either the fourth or fifth $z_{ph}$-bin is removed, the constraining power on $\Omega_m$ and $S_8$ drastically decreases. For the fifth $z_{ph}$-bin, this also affects the constraining power on $A_{IA}$.

<figure><img src="image_11.png" /><figcaption>**Fig. A.1.** Same as black contours in Fig. 9 but while removing one $z_{ph}$-bin at a time.</figcaption></figure>

Now we investigate the accuracy of the modelling strategy of IA described in Sect. 2.2. To validate our IA modelling, we measured the data vector from the cosmo-SLICS+IA for $A_{IA} = \{-1, 0, 1\}$ and then took the ratio of the data vector at $A_{IA} \in \{-1, 1\}$ and divided it by the data vector at $A_{IA} = 0$. This ratio is multiplied with the averaged T17 $\gamma$ mock data vector to infuse IA. For this analysis, we used only equal-scale filter radii $\theta_{ap} \in \{4', 8', 14', 32'\}$. As shown in Fig. A.2, our model nicely recovers the input IA amplitude without shifting the other parameter posteriors. Interestingly, the larger $A_{IA}$ the larger the

<figure><img src="image_12.png" /><figcaption>**Fig. A.2.** Posterior distribution for data vectors and covariance measured from the 1944 T17 noisy $g$ mock data. Here the T17 data vector is infused with IA measured from the cosmo-SLICS+IA. The modelling is described in Sect. 2.2.</figcaption></figure>

FoM, which is due to the increased degeneracy breaking seen in the $\Omega_m-A_{IA}$ panel.

As a further modelling check, we investigate in Fig. A.3 and Fig. A.4 the impact on our model vector of several ingredients, namely the reduced shear correction, the infusion of baryonic feedback processes and intrinsic alignments (IA). As a reference, we always use a model without reduced shear, $A_{IA} = 0$, $A_{ba} = 0$, $S_8 = 0.78$, and $\Omega_m = 0.25$. We then only change one of these effects while fixing the others. We also show the impact on the model vector when we change the cosmology to $S_8 = 0.8$. The change of the model vectors is divided by the square root of the covariance diagonal, indicating the relevant importance given the KiDS-1000 uncertainty. The reduced shear and baryon feedback correction are always subdominant compared to the impact of $S_8$ or IA. For the higher $z_{ph}$-bins, the $S_8$ change dominates, and the IA dominates if $z_{ph}$-bin one or two are included. This result is unsurprising given that the lower $z_{ph}$-bins are more affected by IA as the overall shear signal is much lower. However, we observe especially for the $\langle M_{ap}^3 \rangle$ that the largest IA S/N emerge if lower $z_{ph}$-bins are combined with large $z_{ph}$-bins, which is due to the fact that the large $z_{ph}$-bins drastically deplete the noise. The takeaway from this investigation is that for the current KiDS-1000 analysis, our model correction for the reduced shear and baryon feedback is sufficient. The biggest impact, especially if lower $z_{ph}$-bins are included, comes from IA, which is the biggest concern of current weak lensing analyses.
## Appendix A: Consistency checks and modelling choices

As a first consistency check, we again inferred the joint analysis posteriors but removed one $z_{ph}$-bin at a time. The results are shown in Fig. A.1 and reveal (most importantly) that all $z_{ph}$-bins are internally consistent with each other. However, we further observe that $z_{ph}$-bin one and $z_{ph}$-bin two are not important for inferring cosmological parameters, as expected, given their S/N as shown in Porth et al. (2023). For the $A_{IA}$ parameter, in turn, $z_{ph}$-bin one is very important, as the low redshift is most sensitive to IA. Next, we observe that the third $z_{ph}$-bin results are basically the same as for all five $z_{ph}$-bins, meaning that the third bin can eventually be discarded in future KiDS analyses if the dimension needs to be decreased further. Lastly, if either the fourth or fifth $z_{ph}$-bin is removed, the constraining power on $\Omega_m$ and $S_8$ drastically decreases. For the fifth $z_{ph}$-bin, this also affects the constraining power on $A_{IA}$.

<figure><img src="image_12.png" /><figcaption>**Fig. A.1.** Same as black contours in Fig. 9 but while removing one $z_{ph}$-bin at a time.</figcaption></figure>

Now we investigate the accuracy of the modelling strategy of IA described in Sect. 2.2. To validate our IA modelling, we measured the data vector from the cosmo-SLICS+IA for $A_{IA} = \{-1, 0, 1\}$ and then took the ratio of the data vector at $A_{IA} \in \{-1, 1\}$ and divided it by the data vector at $A_{IA} = 0$. This ratio is multiplied with the averaged T17 $\gamma$ mock data vector to infuse IA. For this analysis, we used only equal-scale filter radii $\theta_{ap} \in \{4', 8', 14', 32'\}$. As shown in Fig. A.2, our model nicely recovers the input IA amplitude without shifting the other parameter posteriors. Interestingly, the larger $A_{IA}$ the larger the

<figure><img src="image_13.png" /><figcaption>**Fig. A.2.** Posterior distribution for data vectors and covariance measured from the 1944 T17 noisy $g$ mock data. Here the T17 data vector is infused with IA measured from the cosmo-SLICS+IA. The modelling is described in Sect. 2.2.</figcaption></figure>

FoM, which is due to the increased degeneracy breaking seen in the $\Omega_m-A_{IA}$ panel.

As a further modelling check, we investigate in Fig. A.3 and Fig. A.4 the impact on our model vector of several ingredients, namely the reduced shear correction, the infusion of baryonic feedback processes and intrinsic alignments (IA). As a reference, we always use a model without reduced shear, $A_{IA} = 0$, $A_{ba} = 0$, $S_8 = 0.78$, and $\Omega_m = 0.25$. We then only change one of these effects while fixing the others. We also show the impact on the model vector when we change the cosmology to $S_8 = 0.8$. The change of the model vectors is divided by the square root of the covariance diagonal, indicating the relevant importance given the KiDS-1000 uncertainty. The reduced shear and baryon feedback correction are always subdominant compared to the impact of $S_8$ or IA. For the higher $z_{ph}$-bins, the $S_8$ change dominates, and the IA dominates if $z_{ph}$-bin one or two are included. This result is unsurprising given that the lower $z_{ph}$-bins are more affected by IA as the overall shear signal is much lower. However, we observe especially for the $\langle M_{ap}^3 \rangle$ that the largest IA S/N emerge if lower $z_{ph}$-bins are combined with large $z_{ph}$-bins, which is due to the fact that the large $z_{ph}$-bins drastically deplete the noise. The takeaway from this investigation is that for the current KiDS-1000 analysis, our model correction for the reduced shear and baryon feedback is sufficient. The biggest impact, especially if lower $z_{ph}$-bins are included, comes from IA, which is the biggest concern of current weak lensing analyses.
<figure><img src="image_13.png" /><figcaption>**Fig. A.3.** Illustration of several effects that change the $E_n$ model vector. The fiducial model $E_n^{\text{fid}}$ is without reduced shear, $A_{\text{IA}} = 0$, $A_{\text{ba}} = 0$, $S_8 = 0.78$ and $\Omega_m = 0.25$. For the different points, one of the effects is changed. We scaled the model differences by the KiDS-1000 uncertainty, indicating the relevance of the change in the context of this work.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Fig. A.4.** Same as Fig. A.3 but for the $\langle M_{\text{ap}}^3 \rangle$ model vector.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Fig. A.3.** Illustration of several effects that change the $E_n$ model vector. The fiducial model $E_n^{\text{fid}}$ is without reduced shear, $A_{\text{IA}} = 0$, $A_{\text{ba}} = 0$, $S_8 = 0.78$ and $\Omega_m = 0.25$. For the different points, one of the effects is changed. We scaled the model differences by the KiDS-1000 uncertainty, indicating the relevance of the change in the context of this work.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Fig. A.4.** Same as Fig. A.3 but for the $\langle M_{\text{ap}}^3 \rangle$ model vector.</figcaption></figure>
## Appendix B: Additional material

This section offers a collection of figures to support the analysis in the main text. In Fig. B.1, we show the full $\langle M_{\text{ap}}^3 \rangle$ data vector

for equal filter radii, from which we showed a subpart in Fig. 3. In Fig. B.2, we show the $\langle M_{\text{ap}}^3 \rangle$ data vector if non-equal filter radii are used and measured from the convergence T17 maps.

<figure><img src="image_15.png" /><figcaption>**Fig. B.1.** Same as Fig. 3, but here the measured data and model are the $\langle M_{\text{ap}}^3 \rangle$ vector for equal-scale aperture filter radii $\theta_{\text{ap}} \in \{4', 6', 8', 10', 14', 18', 22', 26', 32', 36'\}$ in the T17 mock data for $z_{\text{ph}}$-bin combinations.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Fig. B.2.** Measured data and model $\langle M_{\text{ap}}^3 \rangle$ vector for all combinations of aperture filter radii $\theta_{\text{ap}} \in \{4', 8', 16', 32'\}$. The data vector results from one full-sky convergence realisation without shape noise. The grey band shows the expected uncertainty of KiDS-1000 estimated with the convergence maps. The $\langle M_{\text{ap}}^3 \rangle$ is scaled by the third root of the product of the corresponding filter radii.</figcaption></figure>
## Appendix B: Additional material

This section offers a collection of figures to support the analysis in the main text. In Fig. B.1, we show the full $\langle M_{\text{ap}}^3 \rangle$ data vector

for equal filter radii, from which we showed a subpart in Fig. 3. In Fig. B.2, we show the $\langle M_{\text{ap}}^3 \rangle$ data vector if non-equal filter radii are used and measured from the convergence T17 maps.

<figure><img src="image_16.png" /><figcaption>**Fig. B.1.** Same as Fig. 3, but here the measured data and model are the $\langle M_{\text{ap}}^3 \rangle$ vector for equal-scale aperture filter radii $\theta_{\text{ap}} \in \{4', 6', 8', 10', 14', 18', 22', 26', 32', 36'\}$ in the T17 mock data for $z_{\text{ph}}$-bin combinations.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Fig. B.2.** Measured data and model $\langle M_{\text{ap}}^3 \rangle$ vector for all combinations of aperture filter radii $\theta_{\text{ap}} \in \{4', 8', 16', 32'\}$. The data vector results from one full-sky convergence realisation without shape noise. The grey band shows the expected uncertainty of KiDS-1000 estimated with the convergence maps. The $\langle M_{\text{ap}}^3 \rangle$ is scaled by the third root of the product of the corresponding filter radii.</figcaption></figure>