© The Authors 2023

1. Introduction

The Λ cold dark matter (ΛCDM) model has been considered the standard model of cosmology for the past few decades. This relatively simple, six-parameter model describes a wide range of observations, from the cosmic microwave background (CMB) to the observed large-scale structure (LSS) of galaxies, with remarkable accuracy. As the reported uncertainties on cosmological parameters are approaching the per-cent level, a few tensions arise between CMB observations of the early Universe and observations of the local Universe that quantify the LSS (for example, in the Hubble parameter H₀, see Di Valentino et al. 2021a, and references therein). In the past few years, the matter clustering parameter $S_8={\sigma }_8\sqrt{{\mathrm{\Omega }}_{\text{m}}/0.3}$ has also become subject to tension (Hildebrandt et al. 2017; Planck Collaboration VI 2020; Joudaki et al. 2020; Heymans et al. 2021; Abbott et al. 2022; Di Valentino et al. 2021b, and references therein): The local Universe seems less clustered than observations of the CMB suggest. Assuming that these tensions are not due to unknown systematic effects, extensions to the ΛCDM model need to be explored. One of the most popular extensions is the wCDM model, where the equation of state of dark energy differs from w = −1.

The dark energy task force has established that the weak gravitational lensing effect from the LSS, also called cosmic shear, poses one of the most promising methods to constrain the equation of state of dark energy (Albrecht et al. 2006). 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) will be able to constrain potential extensions to the ΛCDM model and may help to decipher the nature of dark energy.

Tight constraints on cosmological parameters are essential to discriminate between the different modifications of the ΛCDM model. So far, two-point statistics have been established as the main analysis tool for cosmic shear (Schneider et al. 1998; Troxel et al. 2018; Hildebrandt et al. 2017, 2020; Hikage et al. 2019; Asgari et al. 2020). These statistics capture the entire information content of a Gaussian random field. Since the initial density field of the Universe is believed to be Gaussian, two-point statistics capture a large amount of cosmological information. However, in late times non-linear structure formation has introduced non-Gaussian features at the smaller scales of the matter distribution, whose information content cannot be captured by two-point statistics. To use this information, a variety of higher-order statistics has been introduced in recent years, including 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. 2022), and many others. In this work, we consider third-order shear statistics, which measure the skewness of the LSS at various scales.

In contrast to most of these higher-order statistics, three-point statistics in configuration space can be directly modelled from a matter bispectrum, allowing for a broad range of consistency checks that can be performed. Furthermore, their modelling does not require simulations that are adjusted to specific survey properties such as the number density and redshift distribution of sources, which allows us to apply them to any data set easily. However, these natural extensions to two-point statistics have yet received surprisingly little attention. Several papers have reported a massive potential information gain when combining two- and three-point statistics (Takada & Jain 2004; Kilbinger & Schneider 2005; Sato & Nishimichi 2013; Kayo et al. 2013). Fu et al. (2014) performed a combined analysis of two- and three-point statistics on 154 deg² CFHTLenS data, reporting a rather moderate gain in information content. Secco et al. (2022) measured the shear three-point correlation functions and third-order aperture statistics in the third-year data release of the Dark Energy Survey (Flaugher 2005; Sevilla-Noarbe et al. 2021), showing that they can be detected with high signal-to-noise. Huterer et al. (2006), Pyne & Joachimi (2021) showed that a combined analysis of two- and three-point statistics has an additional advantage: These two statistics react differently to observational and astrophysical systematics, meaning that a combined analysis allows us to constrain nuisance parameters internally without the need for any additional observations or simulations, yielding an additional null-test and much tighter bounds on cosmological parameters (in an optimistic case, a factor of 20 in the figure of merit of dark energy can be achieved).

This article aims to prepare a cosmological parameter analysis with cosmic shear three-point statistics by developing a pipeline to measure and model both the three-point correlation functions Γ⁽ⁱ⁾ and third-order aperture mass statistics $\langle {\mathcal{M}}_{\text{ap}}^3\rangle$. We compare different estimators, computational costs and information content of both statistics. The covariance calculation will be investigated in the second publication of this series (Linke et al., in prep.). We further show that these statistics can be measured with relative ease in a Stage-III survey, which constitutes a significant advantage over their Fourier-space counterpart, the convergence bispectrum. The aperture mass statistics have several additional advantages: They offer good data compression, are not subject to the mass-sheet degeneracy, and, most importantly, they decompose the signal into E- and B-modes, where to leading order only E-modes can be created by gravitational lensing.

In total, our modelling and validation pipeline can be summarised in the diagram in Fig. 1. Our modelling algorithm is publicly available¹. The paper is structured as follows: In Sect. 2 we introduce the N-body simulations we use to validate and test our modelling pipeline. We then present the convergence bispectrum in Sect. 3, the shear three-point correlation functions in Sect. 4 and the third-order aperture statistics in Sect. 5. For each of these statistics, we describe their theoretical background, how we chose to model them, how they are measured in simulations, and the validation tests we performed. We then compare their information content to one of the second-order shear statistics in a mock-MCMC in Sect. 6 and discuss our findings in Sect. 7.

Fig. 1. Schematic diagram of the modelling and validation pipeline introduced in this paper. The numbers on the arrows correspond to the respective section where this part of the pipeline is discussed.

2. Model validation and covariance estimation with N-body simulations

We use N-body simulations containing only dark matter to validate our model and estimate covariance matrices. One of the main advantages of third-order shear statistics is that they can be easily adapted to different survey specifications. To highlight this, we use different simulation suites with varying source redshift distributions, galaxy number densities and cosmologies for the validation and parameter estimation. In this paper, we use the full-sky gravitational lensing simulations described in Takahashi et al. (2017, hereafter T17), the Millennium Simulations (Springel et al. 2005; Hilbert 2008, hereafter MS), and the Scinet-LIghtCone Simulations (Harnois-Déraps et al. 2018, hereafter SLICS).

2.1. T17 simulations

The T17 are used in this work to perform a realistic analysis of a survey that mimics the KiDS-1000 data and are constructed from a series of nested cubic boxes with side lengths of L, 2L, 3L... placed around a fixed vertex representing the observer’s position, with L = 450 Mpc/h. Each box is replicated eight times and placed around the observer using periodic boundary conditions. With the N-body code GADGET2 (Springel et al. 2001) the gravitational evolution of 2048³ dark matter particles is simulated. Within each box, three spherical lens shells are constructed, each with a width of 150 Mpc/h, which are then used by the public code GRAYTRIX² to trace the light-ray trajectories from the observer to the last scattering surface³. The cosmological parameters of the simulation are Ω_m = 1 − Ω_Λ = 0.279, Ω_b = 0.046, h = 0.7, σ₈ = 0.82, and n_s = 0.97. The matter power spectrum agrees with theoretical predictions of the revised HALOFIT (Takahashi et al. 2012) within 5%(10%) for k < 5(6) h Mpc⁻¹ at z < 1.

For each of the 108 realisations, we build a realistic convergence map by taking a weighted average of all 38 convergence shells at different redshifts, where the weights were determined by the fiducial KiDS-1000 n(z) – see Fig. 2.

Fig. 2. Redshift distribution constructed from the T17 simulation given the fiducial n(z) of the KiDS-1000 data.

We then transform the pure convergence maps into realistic convergence maps by adding to each pixel a Gaussian random variable with a vanishing mean and standard deviation of

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

where A_pix is the pixel area of the convergence grid, and the effective number density n_gal = 6.17 arcmin⁻² and σ_ϵ = 0.265 are chosen such that they are consistent with the combined 1–5 tomographic bins of the KiDS-1000 data.

2.2. Millennium simulations

The MS were run with 2160³ particles in a 500 h⁻¹ Mpc box in a flat ΛCDM cosmology with h = 0.73, σ₈ = 0.9, Ω_m = 0.25, Ω_b = 0.045 and n_s = 1. Subsequently, shear- and convergence-maps of 64 independent lines of sight with an area of 4 × 4 deg² each were created at 36 different redshifts between z = 0 and z = 3 (Hilbert 2008; Hilbert et al. 2009). Each map is calculated on a grid of 4096 × 4096 pixel. For each line of sight, we use the full shear- and convergence maps at redshift z = 1. As we use the MS solely to validate our model, we do not add any noise to the maps.

2.3. Scinet-LIghtCone Simulations

The SLICS were run with 1536³ particles in a 505 h⁻¹ Mpc box, filling up 10 × 10 deg² light-cones up to z = 3. All SLICS were run in a flat ΛCDM cosmology with h = 0.69, σ₈ = 0.83, Ω_m = 0.29, Ω_b = 0.047 and n_s = 0.969. From these simulations, we use convergence maps as well as galaxy catalogues with a redshift distribution of

$$\begin{array}{c}\hfill n(z)\propto z^{2}exp[-{\left(\frac{z}{z_0}\right)}^{\beta }],\end{array}$$(2)

with z₀ = 0.637, β = 1.5 and the overall proportionality constant given by normalising the distribution to 30 gal/arcmin². We use mock galaxy catalogues provided for 924 (pseudo-)independent lines of sight. The SLICS are useful to estimate the constraining power of third-order statistics since the three-point correlation function can be calculated relatively quickly on the 100 deg² fields, and the 924 lines of sight enable the determination of a stable covariance matrix.

3. Convergence power spectrum and bispectrum

In this section, we briefly recap the basics of the weak gravitational lensing formalism, focusing on the shear statistics in Fourier space. More detailed reviews can be found in Bartelmann & Schneider (2001), Hoekstra & Jain (2008), Munshi et al. (2008), Bartelmann (2010).

We start by defining the density contrast at comoving position x and redshift z, $\delta (\mathit{x},z)=\frac{\rho (\mathit{x},z)}{\overline{\rho }(z)}-1$, where ρ(x, z) is the matter density at position x and redshift z and $\overline{\rho }(z)$ the average density at redshift z. From this density contrast, we define the convergence κ for sources at redshift z as a line-of-sight integration, weighted by the lensing efficiency

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

where f_K(χ) is the comoving angular diameter distance. We note that throughout this paper, we work in a flat Universe with Ω_m + Ω_Λ = 1, meaning that f_K[χ(z)] = χ(z), where χ(z) is the comoving distance at redshift z. However, everything we present in this section also works for open or closed universes. The convergence can not be directly observed, but it can be recovered from an observed shear field (Kaiser et al. 1995; Seitz et al. 1998; Jeffrey et al. 2020). The relations between shear, convergence and the matter density contrast allow us to relate all second- and third-order shear statistics to the well-understood matter power spectrum P_δ(k, z) and bispectrum B_δ(k₁, k₂, k₃, z).

3.1. Definition of power- and bispectrum

The matter power spectrum P_δ(k, z) and bispectrum B_δ(k₁, k₂, k₃, z) can be defined as

$$\begin{array}{cc}& ⟨\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 }(k_1,z),\hfill \end{array}$$(4)

$$\begin{array}{cc}& ⟨\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 }(k_1,k_2,k_3,z),\hfill \end{array}$$(5)

where $\widehat{\delta }$ describes the Fourier transform of δ and δ_D is the Dirac-delta distribution. The fact that the power- and bispectrum only depend on the moduli of the k-vectors can be easily derived from the statistical isotropy of the Universe.

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

$$\begin{array}{cc}& P_{\kappa }(\ell )=\frac{9{\mathrm{\Omega }}_{\mathrm{m}}^2{H}_0^4}{4c^4}{\displaystyle {\int }_0^{{\chi }_{\max }}\mathrm{d}\chi \frac{g^2(\chi )}{a^2(\chi )}}\hfill \\ \hfill & \times P[\frac{\ell }{f_K(\chi )},z(\chi )],\hfill \end{array}$$(6)

$$\begin{array}{cc}\hfill & B_{\kappa }({\ell }_1,{\ell }_2,{\ell }_3)=\frac{27H_0^6{\mathrm{\Omega }}_{\mathrm{m}}^3}{8c^6}{\displaystyle {\int }_0^{{\chi }_{\max }}\mathrm{d}\chi \frac{g^3(\chi )}{a^3(\chi )f_K(\chi )}}\hfill \\ \hfill & \times B_{\delta }[\frac{{\ell }_1}{f_K(\chi )},\frac{{\ell }_2}{f_K(\chi )},\frac{{\ell }_3}{f_K(\chi )},z(\chi )].\hfill \end{array}$$(7)

Here,

$$\begin{array}{c}\hfill g(\chi )={\displaystyle {\int }_{\chi }^{{\chi }_{\max }}\mathrm{d}{\chi }^{\prime }p({\chi }^{\prime })\frac{f_K({\chi }^{\prime }-\chi )}{f_K({\chi }^{\prime })}}\end{array}$$(8)

describes the lensing efficiency, where p(χ) is the (normalised) comoving distance probability distribution of sources. We note that usually one instead measures a redshift probability distribution n(z); in our modelling pipeline we thus instead write Eqs. (6)–(8) as integrals over the redshift z.

The Limber approximation breaks down for small values of ℓ. For example, the bispectrum from the Limber approximation overestimates the truth by up to an order of magnitude for ℓ ≪ 60 (corresponding to angular scales of roughly 6°), depending on the source redshift distribution (Deshpande & Kitching 2020). However, at these scales, the impact of non-linear structure formation is small, meaning that the matter distribution is well-described by a Gaussian and higher-order statistics like the bispectrum are small. In this work, we only consider shear statistics up to ∼4°; at these scales, we expect the Limber approximation to hold.

Instead of the three ℓ-values ℓ₁, ℓ₂ and ℓ₃, we can also describe the bispectrum as a function of the two vectors ℓ₁, ℓ₂ or their moduli ℓ₁, ℓ₂ and the angle φ between them. We define

$$\begin{array}{c}\hfill b({\ell }_1,{\ell }_2,\phi )=B_{\kappa }({\ell }_1,{\ell }_2,\sqrt{{\ell }_1^2+{\ell }_2^2+2{\ell }_1{\ell }_{2}cos\phi }).\end{array}$$(9)

3.2. Modelling the bispectrum

We use the state-of-the-art B IHALOFIT algorithm (Takahashi et al. 2020) to model the dark matter bispectrum on non-linear scales. In comparison to older bispectrum models (e.g., Gil-Marín et al. 2012; Scoccimarro & Couchman 2001), B IHALOFIT appears to trace the non-equilateral triangles much better: In comparison with N-body simulations, B IHALOFIT retains an accuracy of 10% or better, whereas the other two fitting formulae are subject to errors of more than 200% for squeezed triangle configurations (compare Figs. 5, 12, and 13 of Takahashi et al. 2020). The effects on higher-order shear statistics are substantial, as can be seen in Halder et al. (2021), who modelled a different third-order shear statistic from using both B IHALOFIT and the bispectrum model of Gil-Marín et al. (2012). In a direct comparison with N-body simulations, B IHALOFIT is accurate on all tested scales, whereas the other fitting formula breaks down on scales of ≲30′.

Another advantage of B IHALOFIT is that it only requires a linear power spectrum as its input, compared to a non-linear spectrum in Gil-Marín et al. (2012) or Scoccimarro & Couchman (2001). We use the fitting formula developed by Eisenstein & Hu (1999) to model the linear power spectrum.

One of the main advantages of third-order shear statistics is that one can rigorously test each stage of the modelling pipeline. We perform such a test on our bispectrum model in Appendix A and conclude that the Limber integrated B IHALOFIT bispectrum is consistent with the MS up to ℓ ≲ 10⁴ and deviates by up to 40% for larger values ℓ.

4. Shear three-point correlation functions

4.1. Definition of the shear three-point correlation functions

Shear three-point correlation functions (3pcf) are the natural extension to the widely used two-point correlation functions. Let γ_c = γ₁ + iγ₂ denote the complex shear in Cartesian coordinates. Considering a triangle of galaxies, as a first step, we project the shear γⁱ of each galaxy i to its tangential- and cross-components with respect to a point fixed with respect to the triangle, for example one of its centres,

$$\begin{array}{c}\hfill \gamma \equiv {\gamma }_{\mathrm{t}}+\mathrm{i}{\gamma }_{\times }=-{\gamma }_{\mathrm{c}}{\mathrm{e}}^{-2\mathrm{i}\zeta },\end{array}$$(10)

where ζ is the angle of the projection direction. Schneider & Lombardi (2003) established four natural components of the shear 3pcf, which remain invariant under rotations of the triangle. They are defined as

$$\begin{array}{cc}& {\mathrm{\Gamma }}^{(0)}=⟨\gamma \gamma \gamma ⟩,{\mathrm{\Gamma }}^{(1)}=⟨{\gamma }^{\ast }\gamma \gamma ⟩,\hfill \\ \hfill & {\mathrm{\Gamma }}^{(2)}=⟨\gamma {\gamma }^{\ast }\gamma ⟩,{\mathrm{\Gamma }}^{(3)}=⟨\gamma \gamma {\gamma }^{\ast }⟩,\hfill \end{array}$$(11)

where the ‘^*’ denotes complex conjugation. The choice of the reference point of the projection is to some degree arbitrary. Usually one of the triangles cenres is chosen, most often the orthocenter (the intersection of its three altitudes) or the centroid (the intersection of its three medians), as shown in Fig. 3. The natural components have the advantage that they are invariant under the choice of triangle centre up to multiplication with a complex phase factor, meaning that their moduli are invariant under the choice of triangle centre.

Fig. 3. Example triangle with our used notations. The sidelengths are denoted by x1 = X3 − X2 and cyclic permutations. The interior angles of the triangle are denoted by ϕi, whereas φi are the angles that the sidelengths take with respect to the x-axis. We also show two possible definitions of its centre: The orthocenter O is at the intersection of the three altitudes (blue dashed); the centroid C is at the intersection of its three medians (red dashed). Figure adapted from Schneider et al. (2005).

Parametrising the shear 3pcf by the triangle side-lengths, x₁, x₂ and x₃, where the indices 1, 2 and 3 are ordered in a counter-clockwise direction, the natural components exhibit a nice behaviour concerning cyclic permutation of arguments. While the first natural component Γ⁽⁰⁾ is invariant under cyclic permutations, the other three components transform into each other,

$$\begin{array}{cc}& {\mathrm{\Gamma }}^{(0)}(x_1,x_2,x_3)={\mathrm{\Gamma }}^{(0)}(x_2,x_3,x_1)={\mathrm{\Gamma }}^{(0)}(x_3,x_1,x_2),\hfill \\ \hfill & {\mathrm{\Gamma }}^{(1)}(x_1,x_2,x_3)={\mathrm{\Gamma }}^{(3)}(x_2,x_3,x_1)={\mathrm{\Gamma }}^{(2)}(x_3,x_1,x_2).\hfill \end{array}$$(12)

Similar behaviour can be observed for parity transformations:

$$\begin{array}{cc}\hfill {\mathrm{\Gamma }}^{(0)}(x_1,x_2,x_3)=& {\mathrm{\Gamma }}^{(0)\ast }(x_2,x_1,x_3),\hfill \\ \hfill {\mathrm{\Gamma }}^{(1)}(x_1,x_2,x_3)=& {\mathrm{\Gamma }}^{(1)\ast }(x_1,x_3,x_2),\hfill \\ \hfill {\mathrm{\Gamma }}^{(2)}(x_1,x_2,x_3)=& {\mathrm{\Gamma }}^{(2)\ast }(x_3,x_2,x_1),\hfill \\ \hfill {\mathrm{\Gamma }}^{(3)}(x_1,x_2,x_3)=& {\mathrm{\Gamma }}^{(3)\ast }(x_2,x_1,x_3).\hfill \end{array}$$(13)

4.2. Modelling the shear three-point correlation functions

We model the natural components ${\left\{{\mathrm{\Gamma }}^{(i)}\right\}}_{i=0,1,2,3}$ of the three-point correlation functions of cosmic shear using the methods described in Schneider et al. (2005, hereafter S+05): When we project the shear of all three galaxies to the orthocenter of the triangle, the projection direction at the triangle vertex X_i is orthogonal to the orientation φ_i of the triangle side x_i. Thus, the shear transforms as

$$\begin{array}{c}\hfill {\gamma }^{(\mathrm{o})}({\mathit{X}}_i)={\gamma }_{\mathrm{c}}({\mathit{X}}_i){\mathrm{e}}^{-2\mathrm{i}{\phi }_i}.\end{array}$$(14)

We now utilise the relation between convergence and shear in Fourier space (Kaiser & Squires 1993),

$$\begin{array}{c}\hfill {\widehat{\gamma }}_{\mathrm{c}}(\mathit{\ell })={\mathrm{e}}^{2\mathrm{i}{\beta }_{\ell }}\widehat{\kappa }(\mathit{\ell }),\end{array}$$(15)

where β_ℓ is the polar angle of ℓ, to write

$$\begin{array}{cc}& {\mathrm{\Gamma }}^{(0)}(x_1,x_2,x_3)=⟨{\gamma }^{(\mathrm{o})}({\mathit{X}}_1){\gamma }^{(\mathrm{o})}({\mathit{X}}_2){\gamma }^{(\mathrm{o})}({\mathit{X}}_3)⟩\hfill \\ \hfill & ={\displaystyle \int \frac{{\mathrm{d}}^2{\ell }_1}{{(2\pi )}^2}\int \frac{{\mathrm{d}}^2{\ell }_2}{{(2\pi )}^2}\int \frac{{\mathrm{d}}^2{\ell }_3}{{(2\pi )}^2}⟨\stackrel{\sim }{\kappa }({\mathit{\ell }}_1)\stackrel{\sim }{\kappa }({\mathit{\ell }}_2)\stackrel{\sim }{\kappa }({\mathit{\ell }}_3)⟩}\hfill \\ \hfill & \times exp[-\mathrm{i}({\mathit{\ell }}_1·{\mathit{X}}_1+{\mathit{\ell }}_2·{\mathit{X}}_2+{\mathit{\ell }}_3·{\mathit{X}}_3)]\hfill \\ \hfill & \times exp\left[2\mathrm{i}{\displaystyle \sum _i}({\beta }_i-{\phi }_i)\right].\hfill \end{array}$$(16)

This can then be transformed into equation (15) of S+05⁴:

$$\begin{array}{cc}& {\mathrm{\Gamma }}^{(0)}(x_1,x_2,x_3)=\frac{2\pi }{3}{\displaystyle {\int }_0^{\infty }\frac{\mathrm{d}{\ell }_1{\ell }_1}{{(2\pi )}^2}{\int }_0^{\infty }}\frac{\mathrm{d}{\ell }_2{\ell }_2}{{(2\pi )}^2}{\int }_0^{2\pi }\mathrm{d}\phi \hfill \\ \hfill & \times b({\ell }_1,{\ell }_2,\phi ){\mathrm{e}}^{2\mathrm{i}\overline{\beta }}[{\mathrm{e}}^{\mathrm{i}({\varphi }_1-{\varphi }_2-6{\alpha }_3)}{J}_6(A_3)\hfill \\ \hfill & +{\mathrm{e}}^{\mathrm{i}({\varphi }_3-{\varphi }_2-6{\alpha }_1)}{J}_6(A_1)+{\mathrm{e}}^{i({\varphi }_3-{\varphi }_1-6{\alpha }_2)}{J}_6(A_2)],\hfill \end{array}$$(17)

with

$$\begin{array}{cc}\hfill A_3& ={[{({\ell }_1{x}_2)}^2+{({\ell }_2{x}_1)}^2+x_1{x}_2{\ell }_1{\ell }_{2}cos(\phi +{\varphi }_3)]}^{\frac{1}{2}},\hfill \\ \hfill |{\mathit{\ell }}_1+{\mathit{\ell }}_2{|}^{2}cos2\overline{\beta }& =({\ell }_1^2+{\ell }_2^2)cos\phi +2{\ell }_1{\ell }_2,\hfill \\ \hfill |{\mathit{\ell }}_1+{\mathit{\ell }}_2{|}^{2}sin2\overline{\beta }& =({\ell }_1^2-{\ell }_2^2)sin\phi ,\hfill \\ \hfill A_{3}cos{\alpha }_3& =({\ell }_1{x}_2+{\ell }_2{x}_1)cos\left(\frac{\phi +{\varphi }_3}{2}\right),\hfill \\ \hfill A_{3}sin{\alpha }_3& =({\ell }_1{x}_2-{\ell }_2{x}_1)sin\left(\frac{\phi +{\varphi }_3}{2}\right).\hfill \end{array}$$(18)

The quantities A_1, 2 and α_1, 2 are obtained by cyclic permutation of indices. The angles ϕ_i are the interior angles of the triangle, as shown in Fig. 3. By introducing polar coordinates $R=\sqrt{{\ell }_1^2+{\ell }_2^2}$, ψ = arctan(ℓ₂/ℓ₁), we get

$$\begin{array}{cc}& {\mathrm{\Gamma }}^{(0)}(x_1,x_2,x_3)=\frac{1}{3{(2\pi )}^3}{\displaystyle {\int }_0^{2\pi }\mathrm{d}\phi {\int }_0^{\pi /2}\mathrm{d}\psi {\int }_0^{\infty }\mathrm{d}R}\hfill \\ \hfill & \times R^{3}sin\psi cos\psi b(Rcos\psi ,Rsin\psi ,\phi ){\mathrm{e}}^{2\mathrm{i}\overline{\beta }}\hfill \\ \hfill & \times [{\mathrm{e}}^{\mathrm{i}({\varphi }_1-{\varphi }_2-6{\alpha }_3)}{J}_6(R{A}_3^{\prime })+{\mathrm{e}}^{\mathrm{i}({\varphi }_3-{\varphi }_2-6{\alpha }_1)}{J}_6(R{A}_1^{\prime })\hfill \\ \hfill & +{\mathrm{e}}^{i({\varphi }_3-{\varphi }_1-6{\alpha }_2)}{J}_6(R{A}_2^{\prime })],\hfill \end{array}$$(19)

with

$$\begin{array}{cc}\hfill A_3^{\prime }=\frac{A_3}{R}& =[{(x_{2}cos\psi )}^2+{(x_{1}sin\psi )}^2\hfill \\ \hfill & {+x_1{x}_{2}sin2\psi cos(\phi +{\varphi }_3)]}^{\frac{1}{2}},\hfill \\ \hfill cos2\overline{\beta }& =(cos\phi +2cos\psi sin\psi ),\hfill \\ \hfill sin2\overline{\beta }& =({cos}^2\psi -{sin}^2\psi )sin\phi ,\hfill \\ \hfill A_3^{\prime }cos{\alpha }_3& =(x_{2}cos\psi +x_{1}sin\psi )cos\left(\frac{\phi +{\varphi }_3}{2}\right),\hfill \\ \hfill A_3^{\prime }sin{\alpha }_3& =(x_{2}cos\psi -x_{1}sin\psi )sin\left(\frac{\phi +{\varphi }_3}{2}\right).\hfill \end{array}$$(20)

Defining

$$\begin{array}{c}\hfill E_3={\mathrm{e}}^{\mathrm{i}({\varphi }_1-{\varphi }_2-6{\alpha }_3)},\end{array}$$(21)

and E₁ and E₂ via cyclic permutations of indices, we can write

$$\begin{array}{cc}& {\mathrm{\Gamma }}^{(0)}(x_1,x_2,x_3)=\frac{1}{6{(2\pi )}^5}{\displaystyle {\int }_0^{\pi /2}\mathrm{d}\psi sin2\psi {\int }_0^{2\pi }\mathrm{d}\phi }\hfill \\ \hfill & \times {\mathrm{e}}^{2\mathrm{i}\overline{\beta }}{\displaystyle \sum _{i=1}^3\frac{E_i}{A_i^{\prime 4}}{\int }_0^{\infty }\mathrm{d}R{R}^3}\hfill \\ \hfill & \times b(\frac{R}{A_i^{\prime }}cos(\psi ),\frac{R}{A_i^{\prime }}sin(\psi ),\phi )J_6(R).\hfill \end{array}$$(22)

The R-integration filters the bispectrum with a 6-th order Bessel function, making the integration routine difficult to solve numerically, as the functional form of the bispectrum prevents the application of fast Hankel transform algorithms like FFTLog (Hamilton 2000)⁵. We thus use the method developed in Ogata (2005) to solve the R-integration and integrate the remaining dimensions using the CUBATURE library⁶. To model Γ⁽¹⁾, we apply the same transformations to Eq. (18) of S+05; for Γ⁽²⁾ and Γ⁽³⁾ we perform a cyclic permutation of the input variables as outlined in Eq. (12).

While a triangle of galaxy positions for which we evaluate the three-point correlation function can be described by its three side-lengths x₁, x₂ and x₃, it is certainly not a good idea to use these variables for a binning scheme; for example, for x₁ > x₂ + x₃ a triangle can not be defined, which means that the 3pcf would not be defined for many bins. A better way to bin the triangles was introduced by Jarvis et al. (2004). Assuming x₁ > x₂ > x₃ they defined a triangle via the values r ∈ [0, ∞], u ∈ [0, 1] and v ∈ [ − 1, 1] by

$$\begin{array}{c}\hfill r=x_2,u=\frac{x_3}{x_2},v=\pm \frac{x_1-x_2}{x_3}.\end{array}$$(23)

Here, v is positive for triangles where x₁, x₂ and x₃ are oriented clockwise and negative for a counter-clockwise orientation. This binning choice allows us to bin the triangle size r in logarithmic steps without having bins where the 3pcf is not defined. In all cases, we bin the shear 3pcf logarithmic in r and linear in u and v.

We note that Eq. (12) implies that the four shear 3pcf for x₁ > x₂ > x₃ (as in the binning scheme of Jarvis et al. 2004) already contain the entire information content of the third-order shear signal, as does knowledge of Γ⁽⁰⁾ and Γ⁽¹⁾ for all combinations of x₁, x₂ and x₃. In a similar manner, Eq. (12) implies that Γ⁽ⁱ⁾(r, u, v) = Γ^(i)*(r, u, −v) holds, where the ‘^*’ denotes complex conjugation. To ensure compatibility with results from the measured 3pcf, we transform the modelled functions from the centroid (as used in S+05) to the orthocenter (as used in Jarvis et al. 2004, compare Sect. 4.3).

For a potential cosmological parameter analysis, the three-point correlation functions face a few hurdles: Assuming we bin the three-point correlation functions in 10 bins for each r, u and v, then our data vector consists of 8000 entries. This makes estimating a covariance matrix using simulations practically impossible and leads to a modelling time that is unfeasible even for a non-tomographic analysis.

4.3. Measuring the shear three-point correlation functions

We use the public tree-code TREECORR (Jarvis et al. 2004) to measure the three-point correlation functions Γ_i. This algorithm estimates the quantity

$$\begin{array}{c}\hfill {\widehat{\mathrm{\Gamma }}}^{(0)}=\frac{{\sum }_{\mathit{ijk}}{w}_i{\epsilon }_i{w}_j{\epsilon }_j{w}_k{\epsilon }_k}{{\sum }_{\mathit{ijk}}{w}_i{w}_j{w}_k},\end{array}$$(24)

where the ε_i = ε_t, i + iε_×, i are the observed ellipticities of galaxies and w_i the associated weights. The other natural components ${\widehat{\mathrm{\Gamma }}}^{(1)},\cdots$ are estimated in the same manner. This estimator has the advantage that it is not impacted by the survey geometry: As long as at least one galaxy triplet falls into each bin, it remains unbiased (Simon et al. 2008)⁷. The disadvantage of the estimator is that its computational complexity scales with $\mathcal{O}(N_{\text{gal}}^3)$, which is not feasible to execute even for moderate values of N_gal ≳ 10⁶. That is why TREECORR constructs a hierarchical ball tree out of the galaxy sample and calculates the correlation functions from this tree. This results in a remarkable speed-up and allows us to calculate the shear 3pcf for an ensemble of about 10⁷ source galaxies distributed over a 10 × 10 deg² field in about 1500 CPUh. A disadvantage of the tree-code is that its execution time scales massively with the number of bins: If b is the logarithmic bin size, then the run-time scales roughly with b⁻⁴.

The TREECORR algorithm also has a BINSLOP parameter, which allows balls of the KD-tree to overlap the edges of a bin. This parameter heavily affects computation time, and while the expectation value is relatively stable between different values of BINSLOP, the covariance is subject to change (Secco et al. 2022).

4.4. Validation

We test our developed integration routine described in Sect. 4.2 with a lensing potential for which we can derive analytic expressions for both the convergence bispectrum and the shear 3pcf. As discussed in more detail in Appendix B we found that the analytic 3pcf and the one obtained by integration from the bispectrum agree to the sub-percent level.

To validate our model for the shear 3pcf, we measure the shear signal at a redshift of z = 1 in the MS. We choose to measure the signal in 10³ bins (10 bins in each r, u and v), with logarithmic r-bins from 0′​.1 to 120′. To speed up computation time, we randomly select every tenth pixel of the 4096² pixel grid. As we do not include shape noise, we expect the loss of signal to be small. The results can be seen in Fig. 4. We conclude that we can model the shear 3pcf reliably down to sub-arcminute scales for almost all triangle configurations. Only for almost degenerate, flattened triangle configurations (v > 0.9) do we see that the model and simulations differ significantly. This might signify that B IHALOFIT breaks down at the corresponding triangle configurations in Fourier-space. Alternatively, this might point towards a break-down of the tree-code’s accuracy at these very degenerate triangles. As we can see in Sect. 5.4.2, these points play a negligible role in the conversion to aperture mass statistics. Overall we see that the agreement between shear 3pcf is better than the one at the bispectrum level (compare Appendix A.2). We observe the same effects for the other natural components of the 3pcf (compare Fig. C.1).

Fig. 4. First natural component Γ(0) of the shear three-point correlation functions modelled from B IHALOFIT compared to measurements from the MS that were extracted using TREECORR. Each panel corresponds to one fixed triangle configuration; the u and v values are listed at the top, and the corresponding shape of the triangle is shown on the bottom right. The x-axis corresponds to the length of the red side of the triangle. We show both the real part (red) and the imaginary part (blue) for the model (dashed line) and the simulations (crosses). Under each panel we show the ratio between measurement and model. The error bars denote the error on the mean of the 64 lines of sight of the MS.

5. Aperture mass statistics

5.1. Definition of aperture mass statistics

An alternative way to analyse cosmic shear is via aperture mass maps (Schneider 1996; Bartelmann & Schneider 2001). Their advantage is that they can separate the signal into so-called E- and B-modes (Schneider et al. 2002), where B-modes can, to leading order, not be created by the weak gravitational lensing effect. Thus, the absence of B-modes provides a crucial null-test in a cosmic shear analysis (see, e.g., Hildebrandt et al. 2017; Asgari et al. 2019). Compared to convergence maps (Kaiser & Squires 1993; Seitz & Schneider 2001; Gatti et al. 2022), aperture mass maps are constructed in a way that they do not suffer from the well-known mass-sheet degeneracy (Falco et al. 1985; Schneider & Seitz 1995).

The aperture mass ℳ_ap at position ϑ and filter radius θ are defined as

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

here, U_θ(ϑ) is a compensated filter (i.e. ∫dϑ ϑ U(ϑ) = 0). Given a shear field γ, the aperture mass ℳ_ap and its respective B-mode counterpart ℳ_⊥ can be calculated as

$$\begin{array}{cc}& {\mathcal{M}}_{\mathrm{ap}}(\mathit{\vartheta };\theta )+\mathrm{i}{\mathcal{M}}_{\perp }(\mathit{\vartheta };\theta )={\displaystyle \int }{\mathrm{d}}^2{\vartheta }^{\prime }{Q}_{\theta }(|\mathit{\vartheta }-{\mathit{\vartheta }}^{\mathbf{\prime }}|)\hfill \\ \hfill & \times [{\gamma }_{\mathrm{t}}({\mathit{\vartheta }}^{\mathbf{\prime }})+\mathrm{i}{\gamma }_{\times }({\mathit{\vartheta }}^{\mathbf{\prime }})],\hfill \end{array}$$(26)

where γ_t and γ_× are projected along the vector ϑ′−ϑ (see Eq. (10)) and Q_θ is related to U_θ via

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

For simplicity of notation, we define U_θ(ϑ) = θ⁻²u(ϑ/θ) and denote by $\widehat{u}(\eta )$ the Fourier transform of u. In this work, we opt for the filter function introduced by 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 }(\vartheta )& =\frac{{\vartheta }^2}{4\pi {\theta }^4}exp(-\frac{{\vartheta }^2}{2{\theta }^2}).\hfill \end{array}$$(28)

While the construction of aperture mass maps has its uses (see for example Harnois-Déraps et al. 2021; Heydenreich et al. 2021), we do not care about the structure of an aperture mass map but rather about its statistical properties. We define for arbitrary combinations of E- and B-mode aperture mass statistics

$$\begin{array}{cc}& ⟨{\mathcal{M}}_{\mathrm{ap}}^m{\mathcal{M}}_{\perp }^n⟩({\theta }_1,\dots ,{\theta }_n)=⟨{\mathcal{M}}_{\mathrm{ap}}(\mathit{\vartheta };{\theta }_1)\cdots {\mathcal{M}}_{\mathrm{ap}}(\mathit{\vartheta };{\theta }_m)\hfill \\ \hfill & \times {\mathcal{M}}_{\perp }(\mathit{\vartheta };{\theta }_m+1)\cdots {\mathcal{M}}_{\perp }(\mathit{\vartheta };{\theta }_m+n){⟩}_{\mathit{\vartheta }}.\hfill \end{array}$$(29)

By construction, ⟨ℳ_ap⟩(θ) vanishes. In a parity-symmetric field, all odd powers of B-mode components vanish (Schneider 2003), meaning that the relevant B-mode counterparts to $\langle {\mathcal{M}}_{\text{ap}}^2\rangle$ and $\langle {\mathcal{M}}_{\text{ap}}^3\rangle$ are $\langle {\mathcal{M}}_{\perp }^2\rangle$ and $\langle {\mathcal{M}}_{\text{ap}}{\mathcal{M}}_{\perp }^2\rangle$, respectively.

5.2. Modelling aperture mass statistics

Given a model for the convergence power spectrum, the second-order aperture statistics can be calculated as

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

As a model for the non-linear power spectrum, we use the revised HALOFIT model of Takahashi et al. (2012). Equivalently, the third-order aperture statistics $\langle {\mathcal{M}}_{\text{ap}}^3\rangle$ can be derived from a bispectrum model via (compare Jarvis et al. 2004; Schneider et al. 2005)⁸

$$\begin{array}{cc}& ⟨{\mathcal{M}}_{\mathrm{ap}}^3⟩({\theta }_1,{\theta }_2,{\theta }_3)=\frac{2}{{(2\pi )}^3}{\displaystyle {\int }_0^{\infty }\mathrm{d}{\ell }_1{\ell }_1{\int }_0^{\infty }\mathrm{d}{\ell }_2{\ell }_2{\int }_0^{\pi }\mathrm{d}\phi }\hfill \\ \hfill & \times \widehat{u}({\theta }_1{\ell }_1)\widehat{u}({\theta }_2{\ell }_2)\widehat{u}\left({\theta }_3\sqrt{{\ell }_1^2+{\ell }_2^2+2{\ell }_1{\ell }_{2}cos\phi }\right)\hfill \\ \hfill & \times b({\ell }_1,{\ell }_2,\phi ).\hfill \end{array}$$(31)

We use the public CUBATURE library to solve this integration. In our implementation, the integration kernel is executed on a graphics processing unit (GPU), yielding a significant speed-up over parallelisation on central processing units (CPUs).

5.3. Measuring aperture mass statistics

There are three main methods to estimate the third-order aperture mass statistics $\langle {\mathcal{M}}_{\text{ap}}^3\rangle$, first, via the convergence field κ, second, via the shear field or, in practice, from the observed galaxy ellipticities, and third, via the third-order correlation functions Γ⁽ⁱ⁾.

5.3.1. Measuring aperture mass statistics from convergence maps

The most straightforward way is to measure aperture mass maps directly on a convergence field using Eq. (25). In a real survey, this is difficult, as the convergence is not directly observable. In principle, one could compute the aperture mass statistics of a reconstructed convergence field, but this is not a good way to estimate aperture statistics, as the convergence reconstruction yields a convergence map that is necessarily smoothed and potentially also contains other systematic effects caused by masks or boundaries of the survey (Pires et al. 2020). While not really applicable to real data, this method yields a quick and unbiased way to estimate aperture statistics in lightcones from simulations, as for them, convergence maps are readily available. However, one faces the issue of boundary effects when the integral in Eq. (25) extends past the simulation boundary. To avoid this issue, we cut off a slice of width 4θ from the computed aperture mass maps⁹.

To extract the data vectors from the (full-sky) convergence maps of the T17 simulations, we smooth the maps with HEALPY function SMOOTHING, with a given beam window function created by the function BEAM2BL, which in turn is determined by the corresponding U_θ filter. For each filter radius θ, this yields a full-sky aperture mass map ${\widehat{\mathcal{M}}}_{\text{ap}}(\mathit{\vartheta };\theta )$, without the need to cut off boundaries.

5.3.2. Measuring aperture mass statistics from galaxy catalogues

Another way to estimate aperture statistics is from an ensemble of observed galaxy ellipticities, using

$$\begin{array}{c}\hfill {\widehat{\mathcal{M}}}_{\mathrm{ap}}(\mathit{\vartheta };\theta )+\mathrm{i}{\widehat{\mathcal{M}}}_{\perp }(\mathit{\vartheta };\theta )=\frac{1}{n_{\mathrm{gal}}}{\displaystyle \sum _i{Q}_{\theta }(|\mathit{\vartheta }-{\mathit{\vartheta }}_i|)({\epsilon }_{\mathrm{t},i}+\mathrm{i}{\epsilon }_{\times ,i}),}\end{array}$$(32)

where ε_t/× are the observed galaxy ellipticities converted into their tangential and croos components according to Eq. (10); ϑ_i are their respective positions. Here, n_gal can be the global number density of galaxies (Bartelmann & Schneider 2001) or the number density of galaxies within the aperture radius (Martinet et al. 2018). For this work, we define n_gal as the number of galaxies weighted by the Q-filter function:

$$\begin{array}{c}\hfill n_{\mathrm{gal}}={\displaystyle \sum _i{Q}_{\theta }(|\mathit{\vartheta }-{\mathit{\vartheta }}_i|).}\end{array}$$(33)

We tested all three definitions of n_gal using the SLICS and found that, for randomly distributed galaxies, setting n_gal as the number density within the aperture radius or the one weighted by the Q-filter function induces sub-percent differences on the third-order aperture masses $\langle {\mathcal{M}}_{\text{ap}}^3\rangle$. However, setting n_gal as the global galaxy density can induce differences of about 5% in $\langle {\mathcal{M}}_{\text{ap}}^3\rangle$.

In the following, we adopt n_gal to be the number of galaxies weighted by the Q-filter function. Rewriting Eq. (10) as

$$\begin{array}{c}\hfill {\gamma }_t+\mathrm{i}{\gamma }_{\times }=-({\gamma }_1+\mathrm{i}{\gamma }_2)\frac{{(\mathit{\vartheta }-{\mathit{\vartheta }}^{\prime })}^{\ast }}{(\mathit{\vartheta }-{\mathit{\vartheta }}^{\prime })},\end{array}$$(34)

where the vector ϑ − ϑ′ denotes the projection direction of the tangential shear in complex notation, we can rewrite Eqs. (26) and (32) as a convolution. To calculate aperture mass fields, we distribute galaxies on a grid using a cloud-in-cell method. From this, we compute the aperture masses using a Fast Fourier Transform (FFT), allowing us to compute an aperture mass map in 𝒪(N_pix log N_pix) operations. To compute second- and third-order aperture statistics, we apply the same cut-off of 4θ_ap to the aperture mass maps. From these aperture mass maps, we obtain estimates for the second- and third-order aperture statistics by multiplication of the respective aperture mass maps on each pixel, and then taking the average of all pixel values.

5.3.3. Measuring aperture statistics from three-point correlation functions

While the abovementioned method to estimate aperture statistics is extremely fast, it can not be applied to realistic survey data. Assuming a relatively large aperture radius of θ_ap = 30′, we would have to cut off a 2°-strip around every edge or mask in the survey footprint, meaning that we would disregard most of the data. While active research is being conducted to circumvent these problems (Porth et al. 2020; Porth & Smith 2021), the arguably best method to estimate third-order aperture statistics from real data is to derive them from the measured 3pcf, as has been introduced in Jarvis et al. (2004), generalised in Schneider et al. (2005) and applied to survey data in Fu et al. (2014) and Secco et al. (2022). The shear 3pcf can be estimated straightforwardly from a survey with arbitrarily complex geometry, meaning that the converted aperture statistics are not biased by boundary effects. One caveat is that this conversion requires the knowledge of the 3pcf for all triangle configurations, particularly for infinitesimally small or extremely large ones, both of which can not be measured. The incomplete knowledge of the correlation functions can lead to a mixing of E- and B-modes for the aperture statistics (Kilbinger et al. 2006). However, for third-order aperture statistics, this effect appears to be not as severe as for their second-order counterpart (at least for the diagonal part of the aperture statistics, this has been demonstrated in Shi et al. 2014, we are testing this assumption for non-diagonal aperture mass statistics in Sect. 5.4.1). Despite its advantages, this method comes at the price of computation time: Calculating the shear 3pcf for a realistic number of source galaxies takes orders of magnitude longer¹⁰ than the direct estimation of aperture statistics, so calculating the shear 3pcf of an ensemble of simulations (as would be necessary to estimate a covariance matrix) comes at a prohibitively high computational cost.

The method to transform shear three-point correlation functions into aperture statistics is already implemented in TREECORR. To compute the third-order aperture statistics, the quantities

$$\begin{array}{cc}& ⟨\mathit{MMM}⟩({\theta }_1,{\theta }_2,{\theta }_3)=A_1{\displaystyle \int }\mathrm{d}{y}_1{\displaystyle \int }\mathrm{d}{y}_2{\displaystyle {\int }_0^{2\pi }\mathrm{d}\psi }\hfill \\ \hfill & \times {\mathrm{\Gamma }}_{\mathrm{cen}}^{(0)}(y_1,y_2,\psi )F_1(y_1,y_2,\psi ),\hfill \\ \hfill & ⟨MM{M}^{\ast }⟩({\theta }_1,{\theta }_2;{\theta }_3)=A_2{\displaystyle \int \mathrm{d}{y}_1\int \mathrm{d}{y}_2{\int }_0^{2\pi }\mathrm{d}\psi }\hfill \\ \hfill & \times {\mathrm{\Gamma }}_{\mathrm{cen}}^{(3)}(y_1,y_2,\psi )F_2(y_1,y_2,\psi )\hfill \end{array}$$(35)

need to be computed. Here, A_1, 2 and F_1, 2 are the prefactors and filter functions, which are specified in S+05 (see Eqs. (62) and (71)). The aperture statistics $\langle {\mathcal{M}}_{\text{ap}}^3\rangle$ and $\langle {\mathcal{M}}_{\text{ap}}{\mathcal{M}}_{\perp }^2\rangle$ are linear combinations of these quantities:

$$\begin{array}{cc}& ⟨{\mathcal{M}}_{\mathrm{ap}}^3⟩({\theta }_1,{\theta }_2,{\theta }_3)\hfill \\ \hfill & =\mathfrak{R}[⟨M^2{M}^{\ast }⟩({\theta }_1,{\theta }_2;{\theta }_3)+⟨M^2{M}^{\ast }⟩({\theta }_1,{\theta }_3;{\theta }_2)\hfill \\ \hfill & +⟨M^2{M}^{\ast }⟩({\theta }_2,{\theta }_3;{\theta }_1)+⟨M^3⟩({\theta }_1,{\theta }_2,{\theta }_3)]/4,\hfill \\ \hfill & ⟨{\mathcal{M}}_{\mathrm{ap}}{\mathcal{M}}_{\perp }^2⟩({\theta }_1;{\theta }_2,{\theta }_3)\hfill \\ \hfill & =\mathfrak{R}[⟨M^2{M}^{\ast }⟩({\theta }_1,{\theta }_2;{\theta }_3)+⟨M^2{M}^{\ast }⟩({\theta }_1,{\theta }_3;{\theta }_2)\hfill \\ \hfill & -⟨M^2{M}^{\ast }⟩({\theta }_2,{\theta }_3;{\theta }_1)-⟨M^3⟩({\theta }_1,{\theta }_2,{\theta }_3)]/4.\hfill \end{array}$$(36)

We do not consider the quantities $\langle {\mathcal{M}}_{\text{ap}}^2{\mathcal{M}}_{\perp }\rangle$ and $\langle {\mathcal{M}}_{\perp }^3\rangle$, as they vanish for any parity-symmetric field (Schneider 2003). In TREECORR, this method is implemented in the following way: First the Γ⁽¹⁾ and Γ⁽²⁾ are transformed into Γ⁽³⁾ via Eq. (12). For each bin (r, u, v), the transformation matrix $\frac{\text{d}\{r,u,v\}}{\text{d}\{{\mathit{y}}_1,{\mathit{y}}_2,\psi \}}$ is then calculated. The value of the respective integral is computed as the sum of the values of Γ⁽³⁾ multiplied by the determinant of the transformation matrix and the value of the filter functions F_1, 2 at the bin centre. While numerically very cheap, this is probably not the most efficient way to compute that integral. In case the filter functions vary significantly over a bin (compare Figs. 2 and 3 of Schneider et al. 2005), it might be more appropriate to calculate the average of the filter function in a bin, for example. We tried to improve the integration results by interpolating the measured shear 3pcf and performing the same integral, achieving rather moderate improvements. We leave an optimisation of the conversion from shear 3pcf to aperture mass statistics for future work.

5.4. Validation

5.4.1. Binning choice of three-point correlation functions

As a first step, we want to validate the conversion Γ⁽ⁱ⁾⇢$\langle {\mathcal{M}}_{\text{ap}}^3\rangle$ performed by the TREECORR algorithm. In particular, we investigate the number of bins necessary to achieve an unbiased estimate for third-order aperture statistics and quantify the leakage of E- and B-modes. While the latter has already been investigated by Shi et al. (2014), we extend upon these results by using a realistic convergence bispectrum model and by taking into account non-diagonal aperture mass statistics.

Our modelling pipeline provides us with the ability to test this conversion. As we start from the same convergence bispectrum B_κ, the aperture mass statistics achieved by the conversion $B_{\kappa }\stackrel{(22)}{\to }{\mathrm{\Gamma }}^{(i)}\stackrel{(35)}{⤏}⟨{{\mathcal{M}}_{\text{ap}}}^3⟩$ and by direct modelling $B_{\kappa }\stackrel{(31)}{\to }⟨{{\mathcal{M}}_{\text{ap}}}^3⟩$ have to be consistent. Furthermore, the modelled Γ⁽ⁱ⁾ are pure E-mode functions, so any B-modes $\langle {\mathcal{M}}_{\text{ap}}{\mathcal{M}}_{\perp }^2\rangle$ that we observe have to be created by the transformation ${\mathrm{\Gamma }}^{(i)}\stackrel{(35)}{⤏}⟨{{\mathcal{M}}_{\text{ap}}}^3⟩$. These tests would be unfeasible to perform with simulations due to the prohibitively high computational cost of extracting the shear 3pcf from an extensive simulation set for different bin sizes.

The results of our tests can be seen in Fig. 5. We see that as long as the three filter radii θ₁, θ₂ and θ₃ are similar, the conversion appears to work reasonably well. Only when we bin Γ⁽ⁱ⁾ in 7³ bins do we get significant deviations, meaning that this is certainly not a sufficient number of bins. The conversion becomes less accurate when two filter radii are much smaller than the third one, as shown in the top-right corner of Fig. 5. There also the results for 10³ bins show significant deviations, whereas the results for 15³ bins seem overall consistent with the ones from 20³ bins. We also observe a non-negligible amount of B-mode leakage for these cases, even for 15³ and 20³ bins. We are planning to exclude all combinations of filter radii from cosmological parameter analyses where we observe a B-mode leakage of more than 10% for the 3pcf in 15³ bins.

Fig. 5. Third-order aperture mass statistics $\langle {\mathcal{M}}_{\text{ap}}^3\rangle$ when modelled directly from the bispectrum (Eq. (31), blue line) compared to the ones constructed from modelled shear 3pcf (Eq. (35), coloured crosses). The dotted lines denote the respective B-modes. We plot $\langle {\mathcal{M}}_{\text{ap}}^3\rangle$(θ1, θ2, θ3), where θ1 is constant in each row, θ2 is constant in each column and θ3 varies along the x-axis in each panel. The yellow, green, red and purple dots denote the $\langle {\mathcal{M}}_{\text{ap}}^3\rangle$ we get from the shear 3pcf when we bin r, u and v in 7 × 7 × 7, 10 × 10 × 10, 15 × 15 × 15 and 20 × 20 × 20 bins, respectively.

Inspecting Fig. 6 we further note that the function Γ⁽¹⁾(r, u, v) is relatively smooth and well-behaved with respect to r, but strongly varies as a function of u and v, especially when u ≈ 0 and v ≈ ±1. This implies that r can be binned rather coarsely, as long as u and v are finely binned. This is in contrast to binning choices in other studies, for example Secco et al. (2022), who preferred a fine binning in r (55 bins) and a coarser binning in u and v (10 bins).

Fig. 6. Three-point correlation function Γ(1) as a function of its three parameters r (left), u (middle) and v (right). We note that the jump in the imaginary part of Γ(1) around v = 0 is due to the fact that Γ(i)(r, u, v) = Γ(i)*(r, u, −v) holds.

5.4.2. Comparison to N-body simulations

We compare the modelled aperture mass statistics with the ones we measure in the MS in Fig. 7. As expected from our discussions in Sect. 5.4.1, we observe that the conversion from shear 3pcf to third-order aperture masses fails when two aperture radii are small, and the third one is large (as observed in the top-right panel). In these cases, we also register significant B-modes. We also note that in most cases, the uncertainties on the direct measurements are significantly larger than those from the shear 3pcf. Cutting off a stripe around the boundary in order to avoid edge effect in the estimate of $\langle {\mathcal{M}}_{\text{ap}}^3\rangle$ (as discussed in Sect. 5.3), leads to a loss of information compared to measuring the 3pcfs, for which all triplets of available points¹¹ in the field are used.

Fig. 7. Third-order aperture masses from our modelling pipeline using B IHALOFIT (Eq. (31), blue) compared to the ones measured directly from the MS. The direct measurements on the convergence maps (Eq. (25)) are shown in green, and the measurements converted from shear 3pcf (Eq. (35)) are shown in orange. Both statistics were computed from 32 lines of sight of the MS, the error bars denote the error on the mean. We note that for the largest angle of 32′, a direct measurement could not be obtained due to the limited size of the individual light cones.

6. Cosmological parameter estimation

To perform a cosmological parameter inference, we use our described pipeline to create a model vector. For the covariance, we rely on N-body simulations, where we use the method of Percival et al. (2022) to debias the estimated covariance matrix $\stackrel{\sim }{C}$. Given a data vector d and a covariance matrix $\stackrel{\sim }{C}$ measured from n_r simulated survey realisations, the posterior distribution of a model vector m(Θ) depending on n_Θ parameters, is

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

where

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

The power-law index m is

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

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}$$(40)

If m = n_r the formalism of Sellentin & Heavens (2016) is recovered.

Normally, one needs to evaluate this likelihood function in high-dimensional parameter space to perform a cosmological parameter analysis with third-order statistics. The creation of the model vector for third-order aperture statistics with 35 combinations of aperture radii takes about one minute on an NVIDIA A40 GPU, and needs to be evaluated at about 10⁴ points even for sampling methods like POLYCHORD (Handley et al. 2015). This means that a complete likelihood analysis with third-order aperture statistics is possible but takes a long time, whereas a complete likelihood analysis with the 3pcf, where the modelling of the 3pcf in 10³ bins takes two to three hours, is not feasible. To solve this issue, we use a neural network emulator called COSMOPOWER (Spurio Mancini et al. 2022), which was first developed to emulate power spectra but can easily be adapted for arbitrary vectors. Since a neural network emulator needs as many as possible evaluation points, we calculate our model at 7500 points¹² in a four-dimensional Latin hypercube describing a flat wCDM cosmological model, varying the parameters Ω_m, S₈, w₀ and h with flat priors of Ω_m ∈ [0.1, 0.55], S₈ ∈ [0.6, 0.9],w₀ ∈ [ − 2, −0.5] and h ∈ [0.6, 0.82]. We leave all other parameters fixed at the values corresponding to the SLICS, which we introduced in Sect. 2. We then train the emulator with 6500 points and use the remaining 1000 as a validation test, as shown in Fig. C.2. This neural network-based emulator performs extraordinarily well, modelling all statistics with a sub-per-cent accuracy. As our bispectrum model is only accurate to about 10% (Takahashi et al. 2020), the emulator uncertainty plays a negligible role in the modelling process. This emulator yields extremely fast predictions for the third-order statistics, so that we can perform a full likelihood analysis using the Metropolis-Hastings sampling method.

Even though we trained our emulator on the cosmological parameters Ω_m, S₈, w₀ and h, we notice that in the absence of tomography we can not constrain w₀ or h. Therefore we fix these parameters in all subsequent MCMC analyses to their fiducial values w₀ = −1 and h = 0.7. For Ω_m and S₈ we use the same priors that we took for the creation of the training sample.

6.1. Shear three-point correlation functions vs. third-order aperture masses

One aspect of this work is investigating the information loss using aperture mass statistics instead of the 3pcf itself. Although the aperture mass statistics are well suited for a cosmological parameter inference due to their E-/B-mode decomposition and fast modelling times, they should not be used if the loss of information is too severe.

Unfortunately, we cannot quantify the full information content of the shear 3pcf, as the data vector contains about 10⁴ entries, and therefore a reliable covariance matrix is not accessible. To circumvent this problem, we model the shear 3pcf in 10³ bins in r, u and v, where we bin r logarithmically from 0.​′1 to 100′, at 500 of the 7500 training nodes. We then perform a principal component analysis (PCA) to decide on the 20 most relevant principal components of the 3pcf data vector¹³.

Using this PCA, we determine the covariance matrix for the shear 3pcf, which we measure from 200 10 × 10 deg² tiles of the SLICS in the same configuration as the training data. For the model vector, we again use COSMOPOWER, which is trained on the PCA components of the 500 models used in determining the principal components.

In Fig. 8 we compare the constraining power of the PCA analysis of 3pcf to the $\langle {\mathcal{M}}_{\text{ap}}^3\rangle$ analysis, where the covariance for PCA analysis of 3pcf is measured from 200 SLICS realisations; for $\langle {\mathcal{M}}_{\text{ap}}^3\rangle$ we use all 927 available SLICS realisations and take into consideration all combinations of the filter radii 0′​.5, 1′,2′,4′,8′,16′, and 32′. The different number of realisations are considered by Eq. (39). In both cases, the data vector was created by the COSMOPOWER Emulator. It is clearly seen that the constraining power from the PCA analysis of 3pcf is only slightly better than the one from $\langle {\mathcal{M}}_{\text{ap}}^3\rangle$, and this slight difference may very well be explained by the use of different scales between the 3pcf and aperture mass statistics, and the fact that the aperture mass maps are smaller due to the cutoff at the boundaries. Overall, the advantages of the $\langle {\mathcal{M}}_{\text{ap}}^3\rangle$ justify their use, even considering their potentially slightly lower constraining power.

Fig. 8. Comparison of the posterior distributions of different third-order statistics using the SLICS simulations to determine the covariance. In red we show the posteriors of the PCA from the 3pcf and in blue the one resulting from the $\langle {\mathcal{M}}_{\text{ap}}^3\rangle$ statistic.

6.2. Combination of second- and third-order aperture masses

To assess the constraining power for third-order aperture statistics, especially when combined with second-order shear statistics, we perform a mock analysis for a non-tomographic KiDS-1000-like setup. To estimate the covariance matrix, we made use of all 108 realisations of the T17 simulations with a resolution NSIDE = 4096 (corresponding to a pixel size of 0′​.74²). From each realisation, we extract 18 HEALPix squares of size ≈860 deg² that do not share common borders. This results in 1944 independent realisations from which the covariance matrix is estimated, such that with a data vector size of ∼50, the statistical noise of the covariance matrix can be neglected. Since the area of the square patches is slightly larger than the one from KiDS-1000 with ≈777.4 deg², the covariance needs to be re-scaled by a factor of 1.11.

Furthermore, in order to have a data vector as unbiased and noise-free as possible, we estimate it with one full-sky realisation with a resolution NSIDE = 8192 (corresponding to a pixel size of 0.18 arcmin). Lastly, we use the filter scales of (4, 8, 16, 32) arcmin, as the model and simulations are inconsistent for smaller scales.

The resulting posterior distribution is shown in the left panel of Fig. 9, where we first notice that the combination of second- and third-order statistics significantly increases the constraining power, especially in the Ω_m-σ₈ panel to the different degeneracy directions of the individual statistics (compare Table 1). Indeed, a joint analysis increases the constraints on S₈ by 42% with respect to second-order statistics; the constraints on Ω_m and σ₈ increase by at least 68% and 54%, respectively¹⁴. The figure of merit (Albrecht et al. 2006) in the Ω_m-σ₈ plane increases by a factor of 5.88. Finally, we note that the true cosmology is well within 1 σ of the expected KiDS-1000 uncertainty for all three statistics.

Fig. 9. Left panel: Comparison of the posterior distributions of second- and third-order aperture statistics with a joint analysis within a KiDS-1000-like setup. The data vector and covariance matrix are estimated by the T17 simulations, and we use filter scales of (4, 8, 16, 32) arcmin. The Hubble parameter h and the dark energy equation-of-state are fixed to the T17 values. Right panel: Comparison of the combined second- and third-order aperture statistics with the combination of second-order aperture statistics with the equal-scale third-order aperture statistics. The corresponding data vector and model vector can be found in Fig. C.3.

Additionally, we investigate the constraining power if only equal-scale aperture masses $\langle {\mathcal{M}}_{\text{ap}}^3\rangle$(θ, θ, θ) are used. As displayed in the right panel of Fig. 9, the loss of constraining power by the limitation to equal-scale aperture masses is small, although not zero. This is in stark contrast to Kilbinger & Schneider (2005), who found a strong difference in constraining power using a Fisher forecast. However, their analysis was conducted using a covariance matrix from a significantly smaller set of simulations. Our results are roughly in line with the findings of Fu et al. (2014), who found rather marginal differences in an MCMC. We note that our analysis is quite simplified: We only vary two cosmological parameters and do not include any systematic effects. The introduction of additional parameters or tomography might lead to some degeneracies in equal-scale aperture mass statistics that could be broken by combinations of different filter radii. Furthermore, future analyses might be able to use a larger dynamic range of scales. In this case the non-equal filter radii contain information about more flattened and squeezed bispectrum configurations, which can yield additional constraining power.

7. Discussion

In this work, which is the first of a series on cosmological analysis with third-order shear statistics, we have shown that a cosmological parameter analysis with third-order aperture mass statistics is feasible and beneficial for Stage-III surveys. Both the shear 3pcf and the third-order aperture statistics can be modelled from the matter bispectrum, and we found that our models based on the B IHALOFIT bispectrum model are accurate enough for Stage-III surveys. In particular, the flat-sky and Limber approximations are valid for our selected range of scales, so the accuracy of our model is mainly limited by the accuracy of the bispectrum model. We note that we have not yet tested the impact of astrophysical or observational systematics.

We developed a test for binning strategies of the three-point correlation functions in order to obtain unbiased estimates for aperture mass statistics and found that, at our selected scales, a measurement of the three-point correlation functions in 15³ bins yields good results. In particular, we found that the leakage between E- and B-modes is at the percent level, and the bias in the aperture mass statistics is well below the sample variance for a Stage-III survey. This extends upon the findings of Shi et al. (2014), who found a percent-level leakage for diagonal aperture mass statistics utilising a simplified bispectrum model. We emphasise that, in addition to the effect of a minimum scale in the shear 3pcf investigated by Shi et al. (2014), our approach also quantifies the leakage that stems from the binning choices in the shear 3pcf, resulting in an inaccurate evaluation of the integral in Eq. (35).

We have tested the information loss when converting the shear three-point correlation functions to third-order aperture mass statistics by performing a principal component analysis of the three-point correlation functions and comparing the constraining power of the principal components to the one of third-order aperture statistics. We found comparable information content, suggesting that third-order aperture statistics constitute a good data compression method for the shear three-point correlation functions. In addition to an easier modelling, the aperture statistics have the added advantage that they cleanly separate E- and B-modes.

We demonstrate that the computational load of a cosmological parameter analysis with third-order aperture statistics is manageable, particularly when utilising an emulator to speed up the MCMC. We make a COSMOSIS-module of our modelling algorithm publicly available¹⁵.

Finally, we compare the constraining power between second-order, third-order, and joint aperture statistics analysis. While second-order aperture statistics are not being used in modern cosmological parameter analyses, we assume that all second-order shear statistics exhibit similar constraining power and parameter degeneracies (compare Asgari et al. 2021). We find that third-order aperture statistics alone have a lower constraining power than their second-order counterpart, but they exhibit a different degeneracy direction in the Ω_m–σ₈-plane so that a joint analysis almost doubles the constraining power on the structure growth parameter S₈ and increases the figure of merit in the Ω_m–σ₈-plane by a factor of 5.9. However, the information gain predicted by Fisher analyses, especially for the difference between diagonal and full third-order aperture statistics (compare Kilbinger & Schneider 2005), appears overly optimistic. This suggests that a Fisher forecast might not be an optimal tool to forecast parameter constraints and mock sampling methods give more realistic constraints.

While we have demonstrated here that cosmological analyses with third-order shear statistics are feasible and promising, there are steps left to do before applying our methods to a concrete cosmological survey. The first of these is the development of a model for the covariance of third-order statistics, which is essential for a tomographic analysis and will be addressed in the following paper of this series (Linke et al., in prep). Additionally, as mentioned above, our model does not yet incorporate systematic and astrophysical effects, like baryonic feedback or intrinsic alignments of source galaxies (Semboloni et al. 2013; Halder & Barreira 2022; Pyne et al. 2022). We will develop and test strategies for treating these effects in future works of this series.

Acknowledgments

We thank Benjamin Joachimi and Mike Jarvis for providing valuable insights to this project. We would like to thank Joachim Harnois-Déraps for making public the SLICS mock data, which can be found at http://slics.roe.ac.uk/. This work was funded by the TRA Matter (University of Bonn) as part of the Excellence Strategy of the federal and state governments. This work has been supported by the Deutsche Forschungsgemeinschaft through the project SCHN 342/15-1. SH acknowledges support from the German Research Foundation (DFG SCHN 342/13), the International Max-Planck Research School (IMPRS) and the German Academic Scholarship Foundation. Author contributions. All authors contributed to the development and writing of this paper. SH wrote the pipeline to model the bispectrum and shear 3pcf and the methods to measure the third-order statistics. LL implemented the modelling algorithm for aperture mass statistics, the GPU-integration, and the cosmosis module, aside from making various improvements to the codes. PB was responsible for everything regarding the T17 simulations and the MCMC runs, including the COSMOPOWER emulator. PS gave countless valuable insights into third-order shear statistics.

References

Appendix A: Testing the BIHALOFIT bispectrum model

A.1. Measuring the bispectrum

To measure the convergence bispectrum B_κ from simulations, we adapt the estimator developed by Watkinson et al. (2017). While their algorithm has been presented for three-dimensional density fields, it can be adapted to two-dimensional convergence fields.

Given a convergence field κ(θ) and its Fourier transform $\widehat{\kappa }(\mathit{\ell })$, for an ℓ-bin ${\overline{\ell }}_i=[{\ell }_{\min },{\ell }_{\max }]$ we define $\widehat{\kappa }(\mathit{\ell };{\overline{\ell }}_i)$ as

$$\begin{array}{c}\hfill \widehat{\kappa }(\mathit{\ell };{\overline{\ell }}_i)=\{\begin{array}{cc}\widehat{\kappa }(\mathit{\ell })\hfill & {\ell }_{\min }\le |\mathit{\ell }|<{\ell }_{\max }\hfill \\ 0\hfill & \mathrm{otherwise}\hfill \end{array},\end{array}$$(A.1)

and $\kappa (\mathit{\theta };{\overline{\ell }}_i)$ as its inverse Fourier transform. We also define $I(\mathit{\theta };{\overline{\ell }}_i)$ as the inverse Fourier transform of $\widehat{I}(\mathit{\ell };{\overline{\ell }}_i)$ with $\widehat{I}$ defined as in Eq. (A.1):

$$\begin{array}{c}\hfill \widehat{I}(\mathit{\ell };{\overline{\ell }}_i)=\{\begin{array}{cc}1\hfill & {\ell }_{\min }\le |\mathit{\ell }|<{\ell }_{\max }\hfill \\ 0\hfill & \mathrm{otherwise}\hfill \end{array}.\end{array}$$(A.2)

The estimator for the convergence bispectrum is then defined as

$$\begin{array}{c}\hfill B_{\kappa }({\overline{\ell }}_1,{\overline{\ell }}_2,{\overline{\ell }}_3)=\frac{{\mathrm{\Omega }}^2}{N_{\mathrm{pix}}^3}\frac{{\sum }_i^{N_{\mathrm{pix}}}\kappa ({\mathit{\theta }}_i;{\overline{\ell }}_1)\kappa ({\mathit{\theta }}_i;{\overline{\ell }}_2)\kappa ({\mathit{\theta }}_i;{\overline{\ell }}_3)}{{\sum }_i^{N_{\mathrm{pix}}}I({\mathit{\theta }}_i;{\overline{\ell }}_1)I({\mathit{\theta }}_i;{\overline{\ell }}_2)I({\mathit{\theta }}_i;{\overline{\ell }}_3)},\end{array}$$(A.3)

where Ω is the solid angle of the respective field and N_pix is the number of pixels. The advantage of this estimator is its speed: With seven Fourier transforms, we can extract the complete averaged bispectrum of a field, where the three Fourier transforms required for the computation of $I(\mathit{\theta };\overline{{\ell }_i})$ only need to be performed once, even when computing the bispectra of multiple fields. Furthermore, $\kappa (\mathit{\theta };{\overline{\ell }}_i)$ can be stored for computing bispectra of different triangle configurations containing the same $\overline{\ell }$-bin. Nevertheless, the estimator suffers from one significant drawback: The Fourier transform assumes periodicity of the field κ, which is normally not given for convergence maps (in contrast to the three-dimensional N-body simulation cubes, which usually exhibit periodic boundary conditions). Therefore the estimator can be biased for ℓ-scales approaching the scales of either the field or individual pixels. We, therefore, discard scales smaller than 5 pixels or larger than a third of the field size.

A.2. Validation

To validate our implementation of the BIHALOFIT algorithm and the Limber integration, we compare our model bispectrum with one extracted from the MS. For all 64 lines-of-sight, we take a convergence map at redshift z = 1 and use the estimator described in Sect. A.1 to extract the bispectrum for a set of triangle configurations. We compare these to our model predictions in Fig. A.1. We recover the bispectrum signal quite well, although we fall short of the 10-20% accuracy reported in Takahashi et al. (2020). However, that may very well be due to the limited sample size provided in the MS or the smoothing inside the MS induced by the ray-tracing (Hilbert et al. 2009).

Fig. A.1. Top panels: convergence bispectrum at z = 1 extracted from the MS (orange) and modelled from BIHALOFIT (blue). Bottom panels: ratio of measurement over the model. The error bars denote the error on the mean of all 64 lines of sight. For each point, we compute the average in a bin of width Δℓ = 0.13ℓ, both in the measurements and the model. For the model, we use equation (8) of Joachimi et al. (2009) as a weight for different triangle configurations within one bin.

Appendix B: Testing the integration routine for Γ_i

When we define the three-point correlation function of the deflection potential as

$$\begin{array}{c}\hfill ⟨\psi (\mathit{X})\psi (\mathit{Y})\psi (\mathit{Z})⟩=\frac{1}{8{\alpha }^3}{\mathrm{e}}^{-\alpha [{(\mathit{X}-\mathit{Y})}^2+{(\mathit{Y}-\mathit{Z})}^2+{(\mathit{X}-\mathit{Y})}^2]},\end{array}$$(B.1)

we can analytically compute both the three-point correlation functions and the bispectrum. For this, we define

$$\begin{array}{c}\hfill {\partial }_X={\partial }_{X_1}+\mathrm{i}{\partial }_{X_2},{\mathrm{\nabla }}_X^2={\partial }_X{\partial }_X^{\ast },\end{array}$$(B.2)

and use the relations

$$\begin{array}{cc}\hfill ⟨\kappa (\mathit{X})\kappa (\mathit{Y})\kappa (\mathit{Z})⟩=& \left(\frac{1}{2}{\mathrm{\nabla }}_X^2\right)\left(\frac{1}{2}{\mathrm{\nabla }}_Y^2\right)\left(\frac{1}{2}{\mathrm{\nabla }}_Z^2\right)⟨\psi (\mathit{X})\psi (\mathit{Y})\psi (\mathit{Z})⟩,\hfill \end{array}$$(B.3)

$$\begin{array}{cc}\hfill ⟨\gamma (\mathit{X})\gamma (\mathit{Y})\gamma (\mathit{Z})⟩=& \left(\frac{1}{2}{\partial }_X^2\right)\left(\frac{1}{2}{\partial }_Y^2\right)\left(\frac{1}{2}{\partial }_Z^2\right)⟨\psi (\mathit{X})\psi (\mathit{Y})\psi (\mathit{Z})⟩,\hfill \end{array}$$(B.4)

$$\begin{array}{cc}\hfill ⟨\gamma (\mathit{X})\gamma (\mathit{Y}){\gamma }^{\ast }(\mathit{Z})⟩=& \left(\frac{1}{2}{\partial }_X^2\right)\left(\frac{1}{2}{\partial }_Y^2\right)\left(\frac{1}{2}{\partial }_Z^{\ast 2}\right)⟨\psi (\mathit{X})\psi (\mathit{Y})\psi (\mathit{Z})⟩.\hfill \end{array}$$(B.5)