© The Authors 2024

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 matter 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{{\mathrm{\Omega }}_{\mathrm{m}}/0.3}$, where σ₈ is the normalisation of the power spectrum and Ω_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 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 {ℳ}_{\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 {ℳ}_{\text{ap}}^3\rangle$. We found that $\langle {ℳ}_{\text{ap}}^3\rangle$, even though they are combined from the shear 3PCF at different scales, contain a similar amount of information on Ω_m and S₈ as the 3PCF itself. The fact that $\langle {ℳ}_{\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 {ℳ}_{\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 {ℳ}_{\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[χ(z)] = χ(z), where χ(z) is the comoving distance at redshift z. Given the matter density, ρ(x, z), at comoving position, x, and redshift, z, the density contrast is $\delta (\mathit{x},z)=\frac{\rho (\mathit{x},z)}{\overline{\rho }(z)}-1$, where $\overline{\rho }(z)$ is the average matter density at redshift, z. The dimensionless surface mass density or convergence, κ, for sources at redshift, z, is determined by the line-of-sight integration as

$$\begin{array}{c}\hfill \kappa (\mathit{\vartheta },z)=\frac{3{\mathrm{\Omega }}_{\mathrm{m}}{H}_0^2}{2c^2}{\displaystyle {\int }_0^{\chi (z)}\mathrm{d}{\chi }^{\prime }\frac{{\chi }^{\prime }[\chi (z)-{\chi }^{\prime }]}{\chi (z)}\frac{\delta ({\chi }^{\prime }\mathit{\vartheta },z)}{a({\chi }^{\prime })},}\end{array}$$(1)

where ϑ is the angular position on the sky, H₀ is the Hubble constant, and a is the scale factor. The second argument of δ simultaneously describes the radial direction and the cosmological epoch, related through the light-cone condition |c dt|=a(z) dχ.

2.1. Limber projections of power- and bispectrum

Given the Fourier transform $\widehat{\delta }$ of the matter density contrast, the matter power spectrum, P_δδ(k, z), and bispectrum, B_δδδ(k₁, k₂, k₃, z), are:

$$\begin{array}{cc}\hfill ⟨\widehat{\delta }({\mathit{k}}_1,z)\widehat{\delta }({\mathit{k}}_2,z)⟩& ={(2\pi )}^3{\delta }_{\mathrm{D}}({\mathit{k}}_1+{\mathit{k}}_2)P_{\delta \delta }(k_1,z),\hfill \end{array}$$(2)

$$\begin{array}{cc}\hfill ⟨\widehat{\delta }({\mathit{k}}_1,z)\widehat{\delta }({\mathit{k}}_2,z)\widehat{\delta }({\mathit{k}}_3,z)⟩& ={(2\pi )}^3{\delta }_{\mathrm{D}}({\mathit{k}}_1+{\mathit{k}}_2+{\mathit{k}}_3)\hfill \\ \hfill & \times B_{\delta \delta \delta }(k_1,k_2,k_3,z),\hfill \end{array}$$(3)

where δ_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),

$$\begin{array}{cc}\hfill P_{\kappa \kappa }^{(ij)}(\ell )& ={\displaystyle {\int }_0^{{\chi }_{\max }}\mathrm{d}\chi \frac{g^{(i)}(\chi )g^{(j)}(\chi )}{a^2(\chi )}{P}_{\delta \delta }[\frac{\ell +1/2}{\chi },z(\chi )],}\hfill \end{array}$$(4)

$$\begin{array}{cc}\hfill B_{\kappa \kappa \kappa }^{(ijk)}({\ell }_1,{\ell }_2,{\ell }_3)& ={\displaystyle {\int }_0^{{\chi }_{\max }}\mathrm{d}\chi \frac{g^{(i)}(\chi )g^{(j)}(\chi )g^{(k)}(\chi )}{a^3(\chi )\chi }}\hfill \\ \hfill & \times B_{\delta \delta \delta }[\frac{{\ell }_1}{\chi },\frac{{\ell }_2}{\chi },\frac{{\ell }_3}{\chi },z(\chi )],\hfill \end{array}$$(5)

where ℓ₃ = |ℓ₁ + ℓ₂|, and g⁽ⁱ⁾(χ) denotes the lensing efficiency and is defined as

$$\begin{array}{c}\hfill g^{(i)}(\chi )=\frac{3{\mathrm{\Omega }}_{\mathrm{m}}{H}_0^2}{2c^2}{\displaystyle {\int }_{\chi }^{{\chi }_{\max }}\mathrm{d}{\chi }^{\prime }{n}^{(i)}({\chi }^{\prime })\frac{{\chi }^{\prime }-\chi }{{\chi }^{\prime }},}\end{array}$$(6)

with n⁽ⁱ⁾(χ) being the redshift probability distribution of the ith tomographic z_ph-bin. Since the Limber approximation breaks down for small values of ℓ (Kilbinger et al. 2017), we consider only scales below 4° to model $\langle {ℳ}_{\text{ap}}^3\rangle$. In principle, we could have also used the modified Limber approximation and shift ℓ_1, 2, 3 → ℓ_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_δδ, we use the revised HMCODE2020 model of Mead et al. (2021) and for the matter bispectrum B_δδδ 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

$$\begin{array}{c}\hfill f_{\mathrm{IA}}=-\frac{A_{\mathrm{IA}}\overline{C_1}\overline{\rho }(z)}{D_{+}(z)},\end{array}$$(7)

where D₊ is the linear growth factor at redshift z, $\overline{C_1}=5\times {10}^{-14}{M}_{\odot }^{-1}{h}^{-2}$ Mpc³, as calibrated in Brown et al. (2002), and A_IA captures the coupling strength between the matter density and the tidal field.

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

$$\begin{array}{c}\hfill {\delta }_{\mathrm{I}}(\mathit{x},z)=f_{\mathrm{IA}}\delta (\mathit{x},z).\end{array}$$(8)

With this, we find that:

$$\begin{array}{cc}\hfill P_{\delta {\delta }_{\mathrm{I}}}(k,z)& =f_{\mathrm{IA}}{P}_{\delta \delta }(k,z),\hfill \end{array}$$(9)

$$\begin{array}{cc}\hfill P_{{\delta }_{\mathrm{I}}{\delta }_{\mathrm{I}}}(k,z)& =f_{\mathrm{IA}}^2{P}_{\delta \delta }(k,z),\hfill \end{array}$$(10)

where P_δδ(k, z) is the non-linear matter power spectrum. The projected power spectra then follow to

$$\begin{array}{cc}& P_{\mathrm{GI}}^{(ij)}(\ell )={\displaystyle {\int }_0^{\infty }\mathrm{d}\chi \frac{g^{(i)}(\chi )n^{(j)}(\chi )}{a(\chi )\chi }{P}_{\delta {\delta }_{\mathrm{I}}}(\ell /\chi ,\chi ),}\hfill \end{array}$$(11)

$$\begin{array}{cc}& P_{\mathrm{II}}^{(ij)}(\ell )={\displaystyle {\int }_0^{\infty }\mathrm{d}\chi \frac{n^{(i)}({\chi }^{\prime })n^{(j)}(\chi )}{{\chi }^2}{P}_{{\delta }_{\mathrm{I}}{\delta }_{\mathrm{I}}}(\ell /\chi ,\chi ),}\hfill \end{array}$$(12)

and the total projected power spectrum becomes

$$\begin{array}{c}\hfill P^{(ij)}(\ell )=P_{\mathrm{GG}}^{(ij)}(\ell )+P_{\mathrm{GI}}^{(ij)}(\ell )+P_{\mathrm{GI}}^{(ji)}(\ell )+P_{\mathrm{II}}^{(ij)}(\ell ).\end{array}$$(13)

The term $P_{\mathrm{GG}}^{(ij)}(\ell )=P_{\kappa \kappa }^{(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 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

$$\begin{array}{cc}\hfill B_{\delta \delta {\delta }_{\mathrm{I}}}(k_1,k_2,k_3,z)& =f_{\mathrm{IA}}{B}_{\delta \delta \delta }(k_1,k_2,k_3,z),\hfill \end{array}$$(14)

$$\begin{array}{cc}\hfill B_{\delta {\delta }_{\mathrm{I}}{\delta }_{\mathrm{I}}}(k_1,k_2,k_3,z)& =f_{\mathrm{IA}}^2{B}_{\delta \delta \delta }(k_1,k_2,k_3,z),\hfill \end{array}$$(15)

$$\begin{array}{cc}\hfill B_{{\delta }_{\mathrm{I}}{\delta }_{\mathrm{I}}{\delta }_{\mathrm{I}}}(k_1,k_2,k_3,z)& =f_{\mathrm{IA}}^3{B}_{\delta \delta \delta }(k_1,k_2,k_3,z),\hfill \end{array}$$(16)

where B_δδδ is the non-linear matter bispectrum, which we calculated with BIHALOFIT. The projected bispectra are

$$\begin{array}{cc}& B_{\mathrm{GGI}}^{(ijk)}({\ell }_1,{\ell }_2,{\ell }_3)=\hfill \\ \hfill & {\displaystyle {\int }_0^{\infty }\mathrm{d}\chi \frac{g^{(i)}(\chi )g^{(j)}(\chi )n^{(k)}(\chi )}{a^2(\chi ){\chi }^2}{B}_{\delta \delta {\delta }_{\mathrm{I}}}({\ell }_1/\chi ,{\ell }_2/\chi ,{\ell }_3/\chi ,\chi ),}\hfill \end{array}$$(17)

$$\begin{array}{cc}& B_{\mathrm{GII}}^{(ijk)}({\ell }_1,{\ell }_2,{\ell }_3)=\hfill \\ \hfill & {\displaystyle {\int }_0^{\infty }\mathrm{d}\chi \frac{g^{(i)}(\chi )n^{(j)}(\chi )n^{(k)}(\chi )}{a(\chi ){\chi }^3}{B}_{\delta {\delta }_{\mathrm{I}}{\delta }_{\mathrm{I}}}({\ell }_1/\chi ,{\ell }_2/\chi ,{\ell }_3/\chi ,\chi ),}\hfill \end{array}$$(18)

$$\begin{array}{cc}& B_{\mathrm{III}}^{(ijk)}({\ell }_1,{\ell }_2,{\ell }_3)=\hfill \\ \hfill & {\displaystyle {\int }_0^{\infty }\mathrm{d}\chi \frac{n^{(i)}(\chi )n^{(j)}(\chi )n^{(k)}(\chi )}{{\chi }^4}{B}_{{\delta }_{\mathrm{I}}{\delta }_{\mathrm{I}}{\delta }_{\mathrm{I}}}({\ell }_1/\chi ,{\ell }_2/\chi ,{\ell }_3/\chi ,\chi ).}\hfill \end{array}$$(19)

The total projected bispectrum is:

$$\begin{array}{cc}\hfill B^{(ijk)}(\widehat{\ell })& =B_{\mathrm{GGG}}^{(ijk)}(\widehat{\ell })+B_{\mathrm{GGI}}^{(ijk)}(\widehat{\ell })+B_{\mathrm{GGI}}^{(kij)}(\widehat{\ell })+B_{\mathrm{GGI}}^{(jki)}(\widehat{\ell })\hfill \\ \hfill & +B_{\mathrm{GII}}^{(ijk)}(\widehat{\ell })+B_{\mathrm{GII}}^{(kij)}(\widehat{\ell })+B_{\mathrm{GII}}^{(jki)}(\widehat{\ell })+B_{\mathrm{III}}^{(ijk)}(\widehat{\ell }),\hfill \end{array}$$(20)

where the tuple is $\widehat{\ell }=({\ell }_1,{\ell }_2,{\ell }_3)$. The interpretation of all these terms is analogous to the ones from the power spectrum, where $B_{\mathrm{GGG}}^{(ijk)}(\widehat{\ell })=B_{\kappa \kappa \kappa }^{(ijk)}(\widehat{\ell })$ 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 {ℳ}_{\text{ap}}^3\rangle$ bins.

The aperture mass ℳ_ap at position ϑ with filter radius θ_ap is defined through the convergence κ, as follows:

$$\begin{array}{c}\hfill {\mathcal{M}}_{\mathrm{ap}}(\mathit{\vartheta };{\theta }_{\mathrm{ap}})={\displaystyle \int }{\mathrm{d}}^2{\vartheta }^{\prime }{U}_{{\theta }_{\mathrm{ap}}}(|{\mathit{\vartheta }}^{\mathbf{\prime }}|)\kappa (\mathit{\vartheta }+{\mathit{\vartheta }}^{\mathbf{\prime }}),\end{array}$$(21)

where U_θap(ϑ′) is a compensated filter such that ∫dϑ′ ϑ′ U_θap(ϑ′) = 0. The tangential shear component γ_t of the complex shear in Cartesian coordinates γ = γ₁ + iγ₂, is defined as ${\gamma }_{\mathrm{t}}=-\mathfrak{Re}(\gamma {\mathrm{e}}^{-2i\varphi })$, where ϕ is the polar angle of ϑ′. Given the tangential shear γ_t, the aperture mass ℳ_ap can also be calculated as

$$\begin{array}{c}\hfill {\mathcal{M}}_{\mathrm{ap}}(\mathit{\vartheta };{\theta }_{\mathrm{ap}})={\displaystyle \int }{\mathrm{d}}^2{\vartheta }^{\prime }{Q}_{{\theta }_{\mathrm{ap}}}(|{\mathit{\vartheta }}^{\mathbf{\prime }}|){\gamma }_{\mathrm{t}}(\mathit{\vartheta }+{\mathit{\vartheta }}^{\mathbf{\prime }}),\end{array}$$(22)

where Q_θap is related to U_θap via

$$\begin{array}{c}\hfill Q_{{\theta }_{\mathrm{ap}}}(\vartheta )=\frac{2}{{\vartheta }^2}{\displaystyle {\int }_0^{\vartheta }\mathrm{d}{\vartheta }^{\prime }{\vartheta }^{\prime }{U}_{{\theta }_{\mathrm{ap}}}({\vartheta }^{\prime })-U_{{\theta }_{\mathrm{ap}}}(\vartheta ).}\end{array}$$(23)

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

$$\begin{array}{cc}\hfill u(x)& =\frac{1}{2\pi }(1-\frac{x^2}{2}){\mathrm{e}}^{-x^2/2},\widehat{u}(\eta )=\frac{{\eta }^2}{2}{\mathrm{e}}^{-{\eta }^2/2},\hfill \\ \hfill Q_{{\theta }_{\mathrm{ap}}}(\vartheta )& =\frac{{\vartheta }^2}{4\pi {\theta }_{\mathrm{ap}}^4}exp(-\frac{{\vartheta }^2}{2{\theta }_{\mathrm{ap}}^2}).\hfill \end{array}$$(24)

2.4. Modelling aperture mass moments

The expectation value of the aperture mass ⟨ℳ_ap⟩(θ_ap), which approximates the ensemble average over all positions ϑ, vanishes by construction. However, the second-order (variance) of the aperture mass is nonzero and can be calculated as

$$\begin{array}{c}\hfill {⟨{\mathcal{M}}_{\mathrm{ap}}^2⟩}^{(ij)}({\theta }_{\mathrm{ap}})={\displaystyle \int }\frac{\mathrm{d}\ell \ell }{2\pi }{P}_{\kappa \kappa }^{(ij)}(\ell ){\widehat{u}}^2({\theta }_{\mathrm{ap}}\ell ).\end{array}$$(25)

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

$$\begin{array}{cc}\hfill {⟨{\mathcal{M}}_{\mathrm{ap}}^3⟩}^{(ijk)}({\theta }_{\mathrm{ap},1},{\theta }_{\mathrm{ap},2},{\theta }_{\mathrm{ap},3})& =\int \frac{{\mathrm{d}}^2{\ell }_1}{{(2\pi )}^2}\int \frac{{\mathrm{d}}^2{\ell }_2}{{(2\pi )}^2}\hfill \\ \hfill & \times B_{\kappa \kappa \kappa }^{(ijk)}({\ell }_1,{\ell }_2,{\ell }_3)\widehat{u}({\theta }_{\mathrm{ap},1}{\ell }_1)\widehat{u}({\theta }_{\mathrm{ap},2}{\ell }_2)\widehat{u}({\theta }_{\mathrm{ap},3}{\ell }_3),\hfill \end{array}$$(26)

where ℓ₃ = |ℓ₁ + ℓ₂|. Later, we differentiate between equal filter radii, for which θ_ap, 1 = θ_ap, 2 = θ_ap, 3 and non-equal filter radii, where the values of θ_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 ξ_±(ϑ) = ⟨γ_tγ_t⟩(ϑ)±⟨γ_×γ_×⟩(ϑ) on a finite angular range

$$\begin{array}{cc}& E_n^{(ij)}=\frac{1}{2}{\displaystyle {\int }_{{\vartheta }_{\min }}^{{\vartheta }_{\max }}\mathrm{d}\vartheta \vartheta [T_{+n}(\vartheta ){\xi }_{+}^{(ij)}(\vartheta )+T_{-n}(\vartheta ){\xi }_{-}^{(ij)}(\vartheta )],}\hfill \end{array}$$(27)

$$\begin{array}{cc}& B_n^{(ij)}=\frac{1}{2}{\displaystyle {\int }_{{\vartheta }_{\min }}^{{\vartheta }_{\max }}\mathrm{d}\vartheta \vartheta [T_{+n}(\vartheta ){\xi }_{+}^{(ij)}(\vartheta )-T_{-n}(\vartheta ){\xi }_{-}^{(ij)}(\vartheta )],}\hfill \end{array}$$(28)

where T_±n(ϑ) are filter functions with support in [ϑ_min, ϑ_max] (Schneider et al. 2010). If

$$\begin{array}{cc}\hfill {\displaystyle {\int }_{{\vartheta }_{\min }}^{{\vartheta }_{\max }}}\mathrm{d}\vartheta \vartheta T_{+n}(\vartheta )J_0(\ell \vartheta )& ={\displaystyle {\int }_{{\vartheta }_{\min }}^{{\vartheta }_{\max }}\mathrm{d}\vartheta \vartheta T_{-n}(\vartheta )J_4(\ell \vartheta )}\hfill \\ \hfill & :=W_n(\ell ),\hfill \end{array}$$(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 [ϑ_min, ϑ_max]. This is not given, for instance, for second-order aperture mass statistics, as this would require information of ξ_± 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:

$$\begin{array}{cc}& E_n^{(ij)}={\displaystyle {\int }_0^{\infty }\frac{\mathrm{d}\ell \ell }{2\pi }{P}_{\kappa \kappa }^{(ij)}(\ell )W_n(\ell ),}\hfill \end{array}$$(30)

$$\begin{array}{cc}& B_n^{(ij)}={\displaystyle {\int }_0^{\infty }\frac{\mathrm{d}\ell \ell }{2\pi }{P}_{\kappa \kappa ,B}^{(ij)}(\ell )W_n(\ell ).}\hfill \end{array}$$(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 {ℳ}_{\text{ap}}^3\rangle$ can be estimated from data in three ways, which we review below. The first uses the convergence field, κ, 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, κ, is available (e.g., from simulations), the easiest way to measure ℳ_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¹ (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_θap filter.

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

$$\begin{array}{c}\hfill E_n^{(ij)}=⟨{\kappa }^{(i)}(\mathit{\vartheta }){\kappa }_s^{(j)}(\mathit{\vartheta };n)⟩,\end{array}$$(32)

where

$$\begin{array}{c}\hfill {\kappa }_s^{(i)}(\mathit{\vartheta };n)=\frac{1}{2\pi }{\displaystyle \int }{\mathrm{d}}^2{\theta }_{\mathrm{ap}}\prime T_{+n}(|\mathit{\vartheta }-{\mathit{\vartheta }}^{\mathbf{\prime }}|){\kappa }^{(i)}({\mathit{\theta }}_{\mathrm{ap}}^{\prime })\end{array}$$(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 ℳ_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

$$\begin{array}{c}\hfill {{\widehat{\mathcal{M}}}_{\mathrm{ap}}}^{(i)}(\mathit{\vartheta };{\theta }_{\mathrm{ap}})={{\displaystyle [\sum _i{Q}_{{\theta }_{\mathrm{ap}}}(|\mathit{\vartheta }-{\mathit{\vartheta }}_i|)]}}^{-1}{\displaystyle \sum _j}{Q}_{{\theta }_{\mathrm{ap}}}(|\mathit{\vartheta }-{\mathit{\vartheta }}_j|){\epsilon }_{\mathrm{t},j}^{(i)},\end{array}$$(34)

where the sum over the filter functions serves as a normalisation, ε_t are the observed galaxy ellipticities converted into their tangential component, and ϑ_i are their respective positions. We note that we sampled the galaxies on a grid using a cloud-in-cell method², 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 }_{\mathrm{ap}}^{\max }$ to all aperture mass maps, where ${\theta }_{\mathrm{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 ξ_±(θ_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 {ℳ}_{\text{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³. It is collectively called ‘KiDS-1000’ as it covers ∼1000 deg² of images, which is then reduced to an effective area of 777.4 deg² 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, K_s, 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⁴ described in Wright et al. (2020). The redshift distributions of all five tomographic bins 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₈ parameter by (at most) 0.1σ 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.

Fig. 1. Redshift distribution of the five tomographic zph-bins of the KiDS-1000 data.

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³ 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⁵. The cosmological parameters of the matter and vacuum energy density are fixed to Ω_m = 1 − Ω_Λ = 0.279, the baryon density parameter to Ω_b = 0.046, the dimensionless Hubble constant to h = 0.7, the normalisation of the power spectrum to σ₈ = 0.82, and the spectral index to n_s = 0.97. The shear information of the T17 simulations is given in terms of 108 γ and κ full-sky realisations, where 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 γ and κ 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_±, n(θ_ap) or U(θ_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². Since the κ 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:

$$\begin{array}{c}\hfill \sigma =\frac{{\sigma }_{ϵ}}{\sqrt{n_{\mathrm{eff}}{A}_{\mathrm{pix}}}},\end{array}$$(35)

with the pixel area as A_pix = 0.74 arcmin² (NSIDE = 4096), n_eff as the effective galaxy number density that included the lensfit weights, and σ_ϵ 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², we rescaled the covariance by 860 deg²/777.4 deg² ≈ 1.1. The reference data vector is measured from one T17 realisation with a resolution of A_pix = 0.18 arcmin² (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 γ and κ 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 arcmin² (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 = γ/(1 − κ) of each object with a shape noise contribution, ϵ^s, to create observed ellipticities ϵ^obs (Seitz & Schneider 1997), as follows:

$$\begin{array}{c}\hfill {ϵ}^{\mathrm{obs}}=\frac{{ϵ}^{\mathrm{s}}+g}{1+{ϵ}^{\mathrm{s}}{g}^{\ast }}.\end{array}$$(36)

The quantities here are all complex numbers and the asterisk ‘’ indicates complex conjugation. To mock the ϵ^s, we use the observed ellipticities ${ϵ}_{\mathrm{KiDS}}^{\mathrm{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 ϵ^s matches the distribution of ${ϵ}_{\mathrm{KiDS}}^{\mathrm{obs}}$. The estimated mean and dispersion σ_ϵ 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⁶. We only use the first five E_n-moments determined from two-point correlation functions, measured from $0\stackrel{\prime }{.}5$ to 300′ in 400 radial bins. 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 γ⁷ and g is the mean of all 1944 T17 mock data; and for the κ data vector, we used one realisation with resolution A_pix = 0.18 arcmin² (NSIDE = 8192). The analytical model was computed using OSMOSIS (Zuntz et al. 2015a).

Fig. 2. Measured and modelled En 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 Apix = 0.18 arcmin2 and shape noise. The red dots represent the data vector measured from one full-sky T17 realisation with a resolution Apix = 0.18 arcmin2 and no shape noise. The grey band indicates the expected uncertainty from the KiDS-1000 survey.

The $\langle {ℳ}_{\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 {ℳ}_{\text{ap}}^3\rangle$ for aperture filter radii θ_ap ∈ {4′,6′,8′,10′,14′,18′,22′,26′,32′,36′}, where we used only equal scales. We show them here for some specific z_ph-bin combinations. In Fig. B.2, we show $\langle {ℳ}_{\text{ap}}^3\rangle$ also for non-equal aperture filter radii θ_ap ∈ {4′,8′,16′,32′}.

Fig. 3. Same as Fig. 2, but here the measured data and model are the $\langle {ℳ}_{\text{ap}}^3\rangle$ vector for equal-scale aperture filter radii θap ∈ {4′,6′,8′,10′,14′,18′,22′,26′,32′,36′} in the T17 mock data for some selected zph-bin combinations. The full data set is shown in Fig. B.1.

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 Ω_m = 0.2905, Ω_Λ = 0.7095, Ω_b = 0.0473, h = 0.6898, σ₈ = 0.836, w₀ = −1.0 and n_s = 0.969. We use the full set of 50 simulated galaxy catalogues, each covering 100 deg² and reproducing the KiDS-1000 n(z) (see Fig. 1) and galaxy number densities n_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 ${ϵ}_{1/2}^{\mathrm{IA}}$ of these galaxies are computed as

$$\begin{array}{c}\hfill {ϵ}_1^{\mathrm{IA}}=-f_{\mathrm{IA}}({\partial }_{\mathit{xx}}-{\partial }_{\mathit{yy}})\varphi ,{ϵ}_2^{\mathrm{IA}}=-2f_{\mathrm{IA}}{\partial }_{\mathit{xy}}\varphi ,\end{array}$$(37)

where f_IA is defined in Eq. (7), and ϕ is the gravitational potential. The partial derivatives of the gravitational potential describe the Cartesian components of the projected tidal field tensors. The ${ϵ}_{1,2}^{\mathrm{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 × 10⁸ h⁻¹ M_⊙ and gas particle mass 1.4 × 10⁸ h⁻¹ M_⊙. The underlying matter fields were constructed from the Magneticum Pathfinder simulations⁸ in the Run-2 with a comoving volume of side 352 h⁻¹ Mpc and Run-2b with a comoving volume of side 640 h⁻¹ 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², 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^DM, using dark matter only and the data vector, d^DM + BA, where dark matter and baryons are included. With them we modified the model vector, m, to

$$\begin{array}{c}\hfill {\mathit{m}}^{\prime }=[1+A_{\mathrm{ba}}(\frac{{\mathit{d}}^{\mathrm{DM}+\mathrm{BA}}}{{\mathit{d}}^{\mathrm{DM}}}-1)]\mathit{m},\end{array}$$(38)

where we introduce the parameter A_ba. A value of zero for A_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_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_ba = 5. Although extrapolating to these large values is probably incorrect, we introduced A_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 Ω_m and S₈ while fixing the remaining parameters to the T17 cosmology. Additionally, we varied the intrinsic alignment amplitude, A_IA, and the parameter A_ba that accounts for the strengths of the response to baryonic feedback⁹ and was described in Sect. 5.3. Lastly, we always marginalise over the shifts of the redshift distributions δ⟨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 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, $\stackrel{\sim }{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, m, that depends on n_p parameters, p, if the covariance matrix $\stackrel{\sim }{C}$ is measured from n_r mock data, is given by:

$$\begin{array}{c}\hfill \mathit{P}\left[\mathit{m}(\mathit{p})|\mathit{d},\stackrel{\sim }{C}\right]\propto {|\stackrel{\sim }{C}|}^{-\frac{1}{2}}{(1+\frac{{\chi }^2}{n_{\mathrm{r}}-1})}^{-m/2},\end{array}$$(39)

where d is the measured data vector and

$$\begin{array}{c}\hfill {\chi }^2={[\mathit{m}(\mathit{p})-\mathit{d}]}^{\mathrm{T}}{\stackrel{\sim }{C}}^{-1}[\mathit{m}(\mathit{p})-\mathit{d}].\end{array}$$(40)

The power-law index m is

$$\begin{array}{c}\hfill m=n_{\mathrm{p}}+2+\frac{n_{\mathrm{r}}-1+B(n_{\mathrm{d}}-n_{\mathrm{p}})}{1+B(n_{\mathrm{d}}-n_{\mathrm{p}})},\end{array}$$(41)

with n_d being the number of data points and

$$\begin{array}{c}\hfill B=\frac{n_{\mathrm{r}}-n_{\mathrm{d}}-2}{(n_{\mathrm{r}}-n_{\mathrm{d}}-1)(n_{\mathrm{r}}-n_{\mathrm{d}}-4)}.\end{array}$$(42)

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

$$\begin{array}{c}\hfill \mathrm{FoM}=\frac{1}{\sqrt{detC}},\end{array}$$(43)

where C is the parameter covariance matrix of the S₈ − Ω_m plane resulting from the MCMC process.

Depending on the used filter combinations, the modelling of $\langle {ℳ}_{\text{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 {Ω_m, S₈, A_IA} and 5000 model points for the real data analysis in the parameter space {Ω_m, S₈, A_IA, δ⟨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 δ⟨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 {ℳ}_{\text{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¹⁰, 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 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 {ℳ}_{\text{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 arcmin², and that we use the κ 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 {ℳ}_{\text{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.

Fig. 4. Posterior distribution for data vectors and covariance measured from the T17 convergence maps catalogues. Here only specific parts of the $\langle {ℳ}_{\text{ap}}^3\rangle$ data vector are used. Here ‘only equal zph-bins’ means that only the auto tomographic bins are used. Lastly, all non-equal filter radii are discarded for ‘only equal filter radii’.

Next, we investigated the choice of aperture filter radii. The aperture filter radii under consideration are θ_ap ∈ {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 θ_ap ∈ {4′,8′,14′,32′}.

Fig. 5. Posterior distribution for different filter radii combinations if the data vector and covariances are measured from the T17 convergence maps.

Next, we considered the full data vector, meaning the $\langle {ℳ}_{\text{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₈. We start with the element with the highest S/N, then consecutively add those vector elements that maximise the S₈ Fisher information content. This Fisher information is calculated as (Tegmark et al. 1997)

$$\begin{array}{c}\hfill F_{S_8}={\left(\frac{\partial \mathit{m}(S_8)}{\partial S_8}\right)}^T{C}^{-1}\left(\frac{\partial \mathit{m}(S_8)}{\partial S_8}\right),\end{array}$$(44)

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

$$\begin{array}{c}\hfill \frac{\partial \mathit{m}(S_8)}{\partial S_8}\approx \frac{-\mathit{m}(S_8^{++})+8\mathit{m}(S_8^{+})-8\mathit{m}(S_8^{-})+\mathit{m}(S_8^{-})}{12\mathrm{\Delta }{S}_8},\end{array}$$(45)

where $S_8^{\pm }=S_8^{\text{T17}}\pm 0.02S_8^{\text{T17}}$ and $S_8^{\pm \pm }=S_8^{\text{T17}}\pm 0.04S_8^{\text{T17}}$. For this analysis, we used only the equal-scale filter radii θ_ap ∈ {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 {ℳ}_{\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 θ_ap ∈ {4′,8′,14′,32′} results in a better FoM then using the best 215 elements. This is likely because we optimised only the S₈ 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 Ω_m–S₈ 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 {ℳ}_{\text{ap}}^3\rangle$ elements with the equal-scale filter radii θ_ap ∈ {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 = γ/(1 − κ), 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_pix = 0.05 arcmin². 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 γ 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₈ and Ω_m. Compared to the E_n-only, the S₈-Ω_m FoM increases by 93% and for the $\langle {ℳ}_{\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 Ω_m, meaning that the improvements are lower bounds. Compared to the analysis based on κ mock data (see Fig. 5), the posteriors are broader because the covariance based on κ 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 κ analysis is not subject to masks 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 ΔΩ_m < 0.02 σ_Ωm and in the clustering amplitude ΔS₈ < 0.05 σ_S8 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.

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.

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 {ℳ}_{\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 {ℳ}_{\text{ap}}^3\rangle$ with all equal filter radii θ_ap ∈ {4′,8′,14′,32′}. To control whether the model accurately fits the data, we minimised the ${\chi }_{\text{real}}^2$ from the real data and the ${\chi }_{\text{T17}}^2$ from each of the 1944 T17 mock data used to compute the covariance matrix. To estimate the probability of measuring a ${\chi }^2>{\chi }_{\text{real}}^2$ (p-value), we counted the number of ${\chi }_{\text{T17}}^2$ that are greater than ${\chi }_{\text{real}}^2$ divided by 1944. The resulting p-value for the E_n-only, the $\langle {ℳ}_{\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 χ² 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 numerical 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_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.

Fig. 7. Measured and modelled En 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 {ℳ}_{\text{ap}}^3\rangle$ or only En is used.

Fig. 8. Measured and modelled $\langle {ℳ}_{\text{ap}}^3\rangle$ with filter radii θap ∈ {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 {ℳ}_{\text{ap}}^3\rangle$ or only $\langle {ℳ}_{\text{ap}}^3\rangle$ is used.

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₈ by 23% and on Ω_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_IA and A_ba are basically untouched, showing that $\langle {ℳ}_{\text{ap}}^3\rangle$ is not helpful for constraining these nuisance parameters.

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 En are compared to $\langle {ℳ}_{\text{ap}}^3\rangle$ using all available combinations for the aperture filter radii θap ∈ {4′,8′,14′,32′} and redshift bin combinations. We note that the prior (see Table 2) range on Ωm is enlarged compared to the previous validation plots. The FoM values of the En-only or $\langle {ℳ}_{\text{ap}}^3\rangle$-only case should not be compared to Fig. 6, as priors of Ωm bind their contours.

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 ${\mathrm{\Omega }}_{\mathrm{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 Ω_cdmh², Ω_bh² and h, while we fixed h and Ω_b and varied only Ω_m. We also find that if we use the same pipeline as A21 but allow larger Ω_cdmh², the posteriors increase towards larger Ω_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 σ in Ω_m and within 0.14 σ in S₈. Similar to A21, we also perform an internal consistency check, removing one z_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_ph-bin is most important to constrain Ω_m and S₈. 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₈. The IA, however, is important, especially if lower z_ph-bins are included and must be accurately modelled.

Finally, we notice that all statistics are consistent with A_IA = 0, although both statistics alone seem to favour positive A_IA. Interestingly, the joint analysis is shifted to lower A_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_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_ba < 1.4 at 68% confidence with our statistic. Here, we should note that $\langle {ℳ}_{\text{ap}}^3\rangle$ does not help in constraining A_ba, which is probably due to the fact that changes in S₈ absorb all A_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 {ℳ}_{\text{ap}}^3\rangle$ (Schneider et al. 2005). In particular, we incorporated intrinsic alignment 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 {ℳ}_{\text{ap}}^3\rangle$ for θ_ap ∈ {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 θ_ap ∈ {4′,8′,14′,32′}. We also tested whether selecting elements that maximise the Fisher information matrix on S₈ 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₈ = 0.772 ± 0.022 and an ${\mathrm{\Omega }}_{\mathrm{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 Ω_m is also interesting for combined weak lensing and galaxy clustering analysis because the constraining power in Ω_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 Ω_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.

Acknowledgments

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-N]. 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

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 Ω_m and S₈ drastically decreases. For the fifth z_ph-bin, this also affects the constraining power on A_IA.

Fig. A.1. Same as black contours in Fig. 9 but while removing one zph-bin at a time.

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.

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 ∈ { − 1, 1} and divided it by the data vector at A_IA = 0. This ratio is multiplied with the averaged T17 γ mock data vector to infuse IA. For this analysis, we used only equal-scale filter radii θ_ap ∈ {4′,8′,14′,32′}. As shown in Fig. , our model nicely recovers the input IA amplitude without shifting the other parameter posteriors. Interestingly, the larger A_IA the larger the FoM, which is due to the increased degeneracy breaking seen in the Ω_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₈ = 0.78, and Ω_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₈ = 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₈ or IA. For the higher z_ph-bins, the S₈ 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 {ℳ}_{\text{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.

Fig. A.3. Illustration of several effects that change the En model vector. The fiducial model $E_n^{\text{fid}}$ is without reduced shear, AIA = 0, Aba = 0, S8 = 0.78 and Ω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.

Fig. A.4. Same as Fig. A.3 but for the $\langle {ℳ}_{\text{ap}}^3\rangle$ model vector.

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 {ℳ}_{\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 {ℳ}_{\text{ap}}^3\rangle$ data vector if non-equal filter radii are used and measured from the convergence T17 maps.

Fig. B.1. Same as Fig. 3, but here the measured data and model are the $\langle {ℳ}_{\text{ap}}^3\rangle$ vector for equal-scale aperture filter radii θap ∈ {4′,6′,8′,10′,14′,18′,22′,26′,32′,36′} in the T17 mock data for zph-bin combinations.

Fig. B.2. Measured data and model $\langle {ℳ}_{\text{ap}}^3\rangle$ vector for all combinations of aperture filter radii θap ∈ {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 {ℳ}_{\text{ap}}^3\rangle$ is scaled by the third root of the product of the corresponding filter radii.

All Tables

All Figures

	Fig. 1. Redshift distribution of the five tomographic zph-bins of the KiDS-1000 data.
In the text

Fig. 2. Measured and modelled En 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 Apix = 0.18 arcmin2 and shape noise. The red dots represent the data vector measured from one full-sky T17 realisation with a resolution Apix = 0.18 arcmin2 and no shape noise. The grey band indicates the expected uncertainty from the KiDS-1000 survey.

In the text

	Fig. 3. Same as Fig. 2, but here the measured data and model are the $\langle {ℳ}_{\text{ap}}^3\rangle$ vector for equal-scale aperture filter radii θap ∈ {4′,6′,8′,10′,14′,18′,22′,26′,32′,36′} in the T17 mock data for some selected zph-bin combinations. The full data set is shown in Fig. B.1.
In the text

	Fig. 4. Posterior distribution for data vectors and covariance measured from the T17 convergence maps catalogues. Here only specific parts of the $\langle {ℳ}_{\text{ap}}^3\rangle$ data vector are used. Here ‘only equal zph-bins’ means that only the auto tomographic bins are used. Lastly, all non-equal filter radii are discarded for ‘only equal filter radii’.
In the text

	Fig. 5. Posterior distribution for different filter radii combinations if the data vector and covariances are measured from the T17 convergence maps.
In the text

	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.
In the text

	Fig. 7. Measured and modelled En 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 {ℳ}_{\text{ap}}^3\rangle$ or only En is used.
In the text

	Fig. 8. Measured and modelled $\langle {ℳ}_{\text{ap}}^3\rangle$ with filter radii θap ∈ {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 {ℳ}_{\text{ap}}^3\rangle$ or only $\langle {ℳ}_{\text{ap}}^3\rangle$ is used.
In the text

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 En are compared to $\langle {ℳ}_{\text{ap}}^3\rangle$ using all available combinations for the aperture filter radii θap ∈ {4′,8′,14′,32′} and redshift bin combinations. We note that the prior (see Table 2) range on Ωm is enlarged compared to the previous validation plots. The FoM values of the En-only or $\langle {ℳ}_{\text{ap}}^3\rangle$-only case should not be compared to Fig. 6, as priors of Ωm bind their contours.

In the text

	Fig. A.1. Same as black contours in Fig. 9 but while removing one zph-bin at a time.
In the text

	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.
In the text

	Fig. A.3. Illustration of several effects that change the En model vector. The fiducial model $E_n^{\text{fid}}$ is without reduced shear, AIA = 0, Aba = 0, S8 = 0.78 and Ω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.
In the text

	Fig. A.4. Same as Fig. A.3 but for the $\langle {ℳ}_{\text{ap}}^3\rangle$ model vector.
In the text

	Fig. B.1. Same as Fig. 3, but here the measured data and model are the $\langle {ℳ}_{\text{ap}}^3\rangle$ vector for equal-scale aperture filter radii θap ∈ {4′,6′,8′,10′,14′,18′,22′,26′,32′,36′} in the T17 mock data for zph-bin combinations.
In the text

Fig. B.2. Measured data and model $\langle {ℳ}_{\text{ap}}^3\rangle$ vector for all combinations of aperture filter radii θap ∈ {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 {ℳ}_{\text{ap}}^3\rangle$ is scaled by the third root of the product of the corresponding filter radii.

In the text