© The Authors 2024

1. Introduction

The most widely used set of statistics for extracting cosmological information from galaxy or weak gravitational lensing surveys is related to N-point correlation functions (NPCF s). Estimating those quantities from the position or shape of a set of observed tracers, one reveals parts of the statistical properties of the underlying large-scale cosmic structure. In particular, for weak gravitational lensing surveys, the analysis of shear two-point correlation functions (2PCFs) and their harmonic analog, the power spectrum, have made it possible to put tight constraints on the matter clustering parameter, $S_8\equiv {\sigma }_8\sqrt{{\mathrm{\Omega }}_{\mathrm{m}}/0.3}$, where σ₈ is the normalization of the matter power spectrum and Ω_m is the matter density parameter (Heymans et al. 2013; Hildebrandt et al. 2017; Troxel et al. 2018; Hikage et al. 2019; Hamana et al. 2020; Asgari et al. 2021; Amon et al. 2022; Secco et al. 2022; Dalal et al. 2023; Li et al. 2023).

While most of the theoretical and observational work has been done at the second-order level, it is well known that through this choice one is sensitive only to the Gaussian information of the density field. To extract additional information including contributions from nonlinear structure formation, one needs to resort to other higher-order statistics, some of which are related to higher-order moments (Jain & Seljak 1997; Schneider et al. 1998; Vicinanza et al. 2018; Gatti et al. 2020), peaks (Kruse & Schneider 1999; Hamana et al. 2004; Harnois-Déraps et al. 2021), the one-point probability distribution function (Patton et al. 2017; Barthelemy et al. 2020; Boyle et al. 2021; Giblin et al. 2023), or the topological properties (Kratochvil et al. 2012; Parroni et al. 2020; Heydenreich et al. 2022) of the density field. Besides the gain in cosmological information, it has been shown that when performing a joint analysis together with second-order statistics, higher-order statistics can also break degeneracies between parameters of cosmological, astrophysical, or systematic origin, and therefore have the potential to yield much tighter cosmological constraints (Kilbinger & Schneider 2005; Takada & Jain 2003; Semboloni et al. 2008; Pyne & Joachimi 2021).

Higher-order statistics from the correlation function hierarchy consist of the (N >  2)-point correlation functions, of which the three-point correlation function (3PCF) is the most prominent example (Schneider & Lombardi 2003; Zaldarriaga & Scoccimarro 2003). Besides the challenges in accurately modeling cosmological and observational effects on the shape of the 3PCF (Takada & Jain 2003; Heydenreich et al. 2023), it is also computationally much more expensive to estimate them, compared to those of the 2PCF. A naive estimator would require all triplets of tracers to be put into a predefined set of bins describing the triplets’ configuration. The resulting computational complexity of a survey of N_gal tracers is 𝒪(N_gal³), rendering this brute-force approach computationally prohibitive. To facilitate an estimation of the 3PCF of scalar (spin-0) and polar (spin-2) observables for current surveys, kd-tree based methods have been developed (Jarvis et al. 2004; Zhang & Pen 2005; Kilbinger et al. 2014) and applied to a wide range of datasets (Bernardeau et al. 2002; Semboloni et al. 2011; Fu et al. 2014; Secco et al. 2022). Even though the computational speedup compared to the brute-force estimator is remarkable, it is most likely still too time-consuming to apply it to stage IV cosmic shear surveys like Euclid (Laureijs et al. 2011) or the Vera Rubin Observatory Legacy Survey of Space and Time (Ivezić et al. 2019). A different class of estimators has been tailored toward the multipole components of the NPCF s of scalar fields (Chen & Szapudi 2005; Slepian & Eisenstein 2015; Philcox et al. 2022). As the multipoles can be estimated in quadratic time complexity for a discrete galaxy catalog, or with an N_pixlog(N_pix) scaling for data distributed on a grid of N_pix pixels, this estimation strategy makes it feasible to efficiently and accurately measure higher-order correlation functions in large datasets. For example, the first measurement of the 4PCF of the scalar galaxy position field was presented in Hou et al. (2023) and Philcox (2022).

In this work, we extend the multipole decomposition of the scalar 3PCF to the natural components of the 3PCF associated with spin-2 tracer fields, such as weak lensing shear. Additionally, we construct an efficient estimator that adaptively switches between a discrete and grid-based prescription for triangle sides of different lengths. With those modifications, we can measure the full 3PCF of any triangle configuration in a stage-III survey in a few central processing unit (CPU) hours without significantly sacrificing accuracy compared to the brute-force estimator.

This work belongs to a series of papers that aim to undertake a cosmological parameter analysis using third-order shear statistics. Heydenreich et al. (2023) introduces an analytical model for the third-order aperture statistics based on the Bihalofit formula (Takahashi et al. 2020) for the matter bispectrum. In Linke et al. (2023), an analytical covariance matrix for the third-order aperture statistics is derived in the presence of a finite survey extent. Burger et al. (2024) extends the analytical model to a tomographic setting, includes methods to account for the impact of intrinsic alignments and baryonic feedback, and presents a cosmological analysis using a joint data vector of the COSEBIs and the third-order aperture mass statistics measured in the KiDS-1000 data.

This paper is structured as follows. In Sect. 2, we introduce the shear 3PCF and their natural components, Γ_μ. In Sect. 3, we derive our estimator of the Γ_μ and discuss some approximations. In Sect. 4, we validate our estimator on a large suite of N-body simulations, the SLICS ensemble, and compare it with measurements obtained with TREECORR (Jarvis et al. 2004). In Sect. 5, we introduce the third-order aperture measures, show how they can be obtained from our shear 3PCF estimator, and validate our implementation. In Sect. 6, we apply our estimator to the fourth data release of the Kilo Degree Survey (Kuijken et al. 2019, hereafter KiDS-1000). After obtaining an empirical estimate of the expected correlation matrix for the KiDS-1000 data from a large suite of mock catalogs we present our measurement of the third-order aperture statistics in the KiDS-1000 data. We conclude in Sect. 7.

2. Third-order measures of cosmic shear

This section aims to introduce the main concepts related to weak gravitational lensing with an emphasis on the third-order shear correlation functions (for extensive reviews of the weak gravitational lensing, see Bartelmann & Schneider 2001; Kilbinger 2015; Dodelson 2017; Mandelbaum 2018). In order to keep the notation concise, we do not consider tomographic setups in this and the following section but collect the corresponding expressions in Sect. 3.5. Furthermore, we do not consider higher-order effects like reduced shear, source redshift clustering, or astrophysical effects, such as intrinsic alignments.

2.1. Basic notions

The two fundamental quantities in weak gravitational lensing are the convergence, κ, and the shear, γ, which describe the isotropic stretching and distortion experienced by a light bundle propagating through spacetime. The convergence, κ, can be written as a weighted projection of the density contrast, δ, along the line of sight:

$$\begin{array}{c}\hfill \kappa (\mathit{\theta })={\displaystyle {\int }_0^{{\chi }_{\max }}}\mathrm{d}\chi \prime W(\chi \prime )\delta [f_K(\chi \prime )\mathit{\theta };\chi \prime ],\end{array}$$(1)

where the associated projection kernel, W, is defined as

$$\begin{array}{c}\hfill W(\chi )\equiv \frac{3{\mathrm{\Omega }}_{\mathrm{m}}{H}_0^2}{2c^2}\frac{f_K(\chi )}{a(\chi )}{\displaystyle {\int }_{z(\chi )}^{\infty }}\mathrm{d}z\prime n(z\prime )\frac{f_K[\chi (z^{\prime })-\chi ]}{f_K[\chi (z\prime )]},\end{array}$$(2)

and f_K[χ(z)] denotes the comoving angular diameter distance of the comoving radial distance, χ, at redshift, z. We introduced the dimensionless matter density parameter, Ω_m, and the Hubble constant, H₀, as well as the scale factor, a, and the source redshift distribution, n(z). For the remainder of this work, we assume a flat universe for which f_K(χ)=χ.

The components of the complex shear field, γ, can be characterized with respect to a Cartesian coordinate frame, γ_c ≡ γ₁ + iγ₂. To evaluate the shear at position θ in a reference frame, which is rotated by an angle, ζ, with respect to the Cartesian basis, one has

$$\begin{array}{c}\hfill \gamma (\mathit{\theta };\zeta )\equiv -[{\gamma }_{\mathrm{t}}(\mathit{\theta };\zeta )+\mathrm{i}{\gamma }_{\times }(\mathit{\theta };\zeta )]\equiv -{\gamma }_{\mathrm{c}}(\mathit{\theta }){\mathrm{e}}^{-2\mathrm{i}\zeta },\end{array}$$(3)

where we introduced the tangential (γ_t) and cross-components (γ_×) of the shear with respect to the projection direction, ζ.

2.2. The shear three-point correlation function and its natural components

Due to the polar nature of the shear, there are eight different real-valued components for its three-point correlator. Schneider & Lombardi (2003) regroup those components into four complex-valued natural components, Γ_μ, that do not mix and only transform by a particular phase factor under rotations of the corresponding triangle. Parameterizing a triangle as depicted in Fig. 1 and choosing some arbitrary projection, [BOLD]𝒫, in which the individual shear components are rotated by angles, ζ_i, we defined the natural components of the shear 3PCF as¹

$$\begin{array}{cc}\hfill {\mathrm{\Gamma }}_0^{\mathcal{P}}({\vartheta }_1,{\vartheta }_2,\varphi )& =\langle \gamma ({\mathbf{X}}_1;{\zeta }_1)\gamma ({\mathbf{X}}_2;{\zeta }_2)\gamma ({\mathbf{X}}_3;{\zeta }_3)\rangle ,\hfill \end{array}$$(4)

$$\begin{array}{cc}\hfill {\mathrm{\Gamma }}_1^{\mathcal{P}}({\vartheta }_1,{\vartheta }_2,\varphi )& =\langle {\gamma }^{\ast }({\mathbf{X}}_1;{\zeta }_1)\gamma ({\mathbf{X}}_2;{\zeta }_2)\gamma ({\mathbf{X}}_3;{\zeta }_3)\rangle ,\hfill \end{array}$$(5)

$$\begin{array}{cc}\hfill {\mathrm{\Gamma }}_2^{\mathcal{P}}({\vartheta }_1,{\vartheta }_2,\varphi )& =\langle \gamma ({\mathbf{X}}_1;{\zeta }_1){\gamma }^{\ast }({\mathbf{X}}_2;{\zeta }_2)\gamma ({\mathbf{X}}_3;{\zeta }_3)\rangle ,\hfill \end{array}$$(6)

$$\begin{array}{cc}\hfill {\mathrm{\Gamma }}_3^{\mathcal{P}}({\vartheta }_1,{\vartheta }_2,\varphi )& =\langle \gamma ({\mathbf{X}}_1;{\zeta }_1)\gamma ({\mathbf{X}}_2;{\zeta }_2){\gamma }^{\ast }({\mathbf{X}}_3;{\zeta }_3)\rangle ,\hfill \end{array}$$(7)

Fig. 1. Upper panel: Parameterization of a triplet of shears used in this work. For some shears at position X1, we denote the connecting lines to the two other shears as ϑi and the enclosing angle as ϕ. The dashed red lines show the directions of the × projection Eq. (9) for which the three projection axes intersect in X1. Lower panel: Definition of the geometric quantities involved in the centroid projection Eq. (10). We see that the q vectors connect the centroid of the triangle with the triangles’ vertices.

where ϑ_i ≡ |ϑ_i| and we assumed that the projection directions, ζ_i, are defined in terms of the vertices of the triangle; this implies that the natural components are invariant under rotations of the triangle such that a parameterization in terms of three variables is sufficient. We note that any of the natural components can be written in terms of the elements

$$\begin{array}{c}\hfill \{{\gamma }_{\mathrm{ttt}},{\gamma }_{\times \mathrm{tt}},{\gamma }_{\mathrm{t}\times \mathrm{t}},{\gamma }_{\mathrm{tt}\times },{\gamma }_{\mathrm{t}\times \times },{\gamma }_{\times \mathrm{t}\times },{\gamma }_{\times \times \mathrm{t}},{\gamma }_{\times \times \times }\},\end{array}$$(8)

where we introduced shorthand notation like γ_ttt ≡ ⟨γ_tγ_tγ_t⟩ for the various triple products of shears (Schneider & Lombardi 2003). Inverting the relations (4)–(7), one can relate the components in Eq. (8) to the natural components.

2.3. Choice of projections for the shear three-point correlation function

In this work, we employed two particular projections. The first of them, which we dubbed the × projection, is defined by the rotation angles

$$\begin{array}{c}\hfill {\zeta }_1^{\times }=\frac{1}{2}({\phi }_1+{\phi }_2),{\zeta }_2^{\times }={\phi }_1,{\zeta }_3^{\times }={\phi }_2,\end{array}$$(9)

and the corresponding projection axes are depicted in the upper panel of Fig. 1. As we will see later, this projection turns out to be useful for the estimator derived in Sect. 3. For the second projection shown in the lower panel of Fig. 1, we first defined the centroid of the triangle as the intersection of its three medians and then projected the shear at position X_i along the direction, q_i, of the median passing through X_i:

$$\begin{array}{c}\hfill {\mathrm{e}}^{2\mathrm{i}{\zeta }_1^{\mathrm{cent}}}=\frac{{\breve{\mathit{q}}}_1}{{\breve{\mathit{q}}}_1^{\ast }},{\mathrm{e}}^{2\mathrm{i}{\zeta }_2^{\mathrm{cent}}}=\frac{{\breve{\mathit{q}}}_2}{{\breve{\mathit{q}}}_2^{\ast }},{\mathrm{e}}^{2\mathrm{i}{\zeta }_3^{\mathrm{cent}}}=\frac{{\breve{\mathit{q}}}_3}{{\breve{\mathit{q}}}_3^{\ast }},\end{array}$$(10)

where the q_i are given as

$$\begin{array}{c}\hfill {\mathit{q}}_1=-\frac{{\mathit{\vartheta }}_1+{\mathit{\vartheta }}_2}{3},{\mathit{q}}_2=\frac{2{\mathit{\vartheta }}_1-{\mathit{\vartheta }}_2}{3},{\mathit{q}}_3=\frac{2{\mathit{\vartheta }}_2-{\mathit{\vartheta }}_1}{3},\end{array}$$(11)

and where we defined ${\breve{\mathit{q}}}_i\equiv q_{i,1}+\mathrm{i}{q}_{i,2}$. The corresponding ${\mathrm{\Gamma }}_{\mu }^{\mathrm{cent}}$ are useful for connecting the shear 3PCF to the third-order aperture statistics and are being used by other estimators for the 3PCF such as TREECORR.

We can convert between the two projections using Eq. (3). For example, the zeroth natural component in the × projection can be written as

$$\begin{array}{cc}\hfill {\mathrm{\Gamma }}_0^{\times }({\vartheta }_1,{\vartheta }_2,\varphi )& =-⟨{\gamma }_{\mathrm{c}}(\theta ){\gamma }_{\mathrm{c}}(\theta +{\vartheta }_1){\gamma }_{\mathrm{c}}(\theta +{\vartheta }_2)⟩{\mathrm{e}}^{-3\mathrm{i}({\phi }_1+{\phi }_2)}\hfill \\ \hfill & ={\mathrm{\Gamma }}_0^{\mathrm{cent}}({\vartheta }_1,{\vartheta }_2,\varphi )\frac{{\left|{\breve{\mathit{q}}}_1{\breve{\mathit{q}}}_2{\breve{\mathit{q}}}_3\right|}^2}{{\left({\breve{\mathit{q}}}_1^{\ast }{\breve{\mathit{q}}}_2^{\ast }{\breve{\mathit{q}}}_3^{\ast }\right)}^2}{\mathrm{e}}^{-3\mathrm{i}({\phi }_1+{\phi }_2)},\hfill \end{array}$$

from which one can immediately read off Γ₀^cent. Repeating similar computations for the other components and using φ₂ = φ₁ + ϕ, we arrived at

$$\begin{array}{cc}\hfill {\mathrm{\Gamma }}_0^{\mathrm{cent}}({\vartheta }_1,{\vartheta }_2,\varphi )& ={\mathrm{\Gamma }}_0^{\times }({\vartheta }_1,{\vartheta }_2,\varphi )\frac{{\breve{\mathit{q}}}_1^{\ast }}{{\breve{\mathit{q}}}_1}\frac{{\breve{\mathit{q}}}_2^{\ast }}{{\breve{\mathit{q}}}_2}\frac{{\breve{\mathit{q}}}_3^{\ast }}{{\breve{\mathit{q}}}_3}{\mathrm{e}}^{6\mathrm{i}{\phi }_1}{\mathrm{e}}^{3\mathrm{i}\varphi },\hfill \end{array}$$(12)

$$\begin{array}{cc}\hfill {\mathrm{\Gamma }}_1^{\mathrm{cent}}({\vartheta }_1,{\vartheta }_2,\varphi )& ={\mathrm{\Gamma }}_1^{\times }({\vartheta }_1,{\vartheta }_2,\varphi )\frac{{\breve{\mathit{q}}}_1}{{\breve{\mathit{q}}}_1^{\ast }}\frac{{\breve{\mathit{q}}}_2^{\ast }}{{\breve{\mathit{q}}}_2}\frac{{\breve{\mathit{q}}}_3^{\ast }}{{\breve{\mathit{q}}}_3}{\mathrm{e}}^{2\mathrm{i}{\phi }_1}{\mathrm{e}}^{\mathrm{i}\varphi },\hfill \end{array}$$(13)

$$\begin{array}{cc}\hfill {\mathrm{\Gamma }}_2^{\mathrm{cent}}({\vartheta }_1,{\vartheta }_2,\varphi )& ={\mathrm{\Gamma }}_2^{\times }({\vartheta }_1,{\vartheta }_2,\varphi )\frac{{\breve{\mathit{q}}}_1^{\ast }}{{\breve{\mathit{q}}}_1}\frac{{\breve{\mathit{q}}}_2}{{\breve{\mathit{q}}}_2^{\ast }}\frac{{\breve{\mathit{q}}}_3^{\ast }}{{\breve{\mathit{q}}}_3}{\mathrm{e}}^{2\mathrm{i}{\phi }_1}{\mathrm{e}}^{3\mathrm{i}\varphi },\hfill \end{array}$$(14)

$$\begin{array}{cc}\hfill {\mathrm{\Gamma }}_3^{\mathrm{cent}}({\vartheta }_1,{\vartheta }_2,\varphi )& ={\mathrm{\Gamma }}_3^{\times }({\vartheta }_1,{\vartheta }_2,\varphi )\frac{{\breve{\mathit{q}}}_1^{\ast }}{{\breve{\mathit{q}}}_1}\frac{{\breve{\mathit{q}}}_2^{\ast }}{{\breve{\mathit{q}}}_2}\frac{{\breve{\mathit{q}}}_3}{{\breve{\mathit{q}}}_3^{\ast }}{\mathrm{e}}^{2\mathrm{i}{\phi }_1}{\mathrm{e}}^{-\mathrm{i}\varphi },\hfill \end{array}$$(15)

where the rotation angles can be expressed as

$$\begin{array}{cc}\hfill \frac{{\breve{\mathit{q}}}_1^{\ast }}{{\breve{\mathit{q}}}_1}{\mathrm{e}}^{2\mathrm{i}{\phi }_1}& =\frac{{\vartheta }_1+{\vartheta }_2{\mathrm{e}}^{-\mathrm{i}\varphi }}{{\vartheta }_1+{\vartheta }_2{\mathrm{e}}^{\mathrm{i}\varphi }},\frac{{\breve{\mathit{q}}}_2^{\ast }}{{\breve{\mathit{q}}}_2}{\mathrm{e}}^{2\mathrm{i}{\phi }_1}=\frac{{\vartheta }_2{\mathrm{e}}^{-\mathrm{i}\varphi }-2{\vartheta }_1}{{\vartheta }_2{\mathrm{e}}^{\mathrm{i}\varphi }-2{\vartheta }_1},\hfill \\ \hfill \frac{{\breve{\mathit{q}}}_3^{\ast }}{{\breve{\mathit{q}}}_3}{\mathrm{e}}^{2\mathrm{i}{\phi }_1}& =\frac{{\vartheta }_1-2{\vartheta }_2{\mathrm{e}}^{-\mathrm{i}\varphi }}{{\vartheta }_1-2{\vartheta }_2{\mathrm{e}}^{\mathrm{i}\varphi }}.\hfill \end{array}$$

2.4. Bin-averaged shear three-point correlation function

When estimating the Γ_μ^𝒫 from a finite set of galaxy ellipticities, one does not have direct access to every possible triangle shape but instead collects the point triplets into radial and angular bins. Such a measurement then provides an estimator of the bin-averaged shear 3PCF, ${\overline{\mathrm{\Gamma }}}_{\mu }^{\mathcal{P}}$:

$$\begin{array}{cc}\hfill {\overline{\mathrm{\Gamma }}}_{\mu }^{\mathcal{P}}(& {\mathrm{\Theta }}_i,{\mathrm{\Theta }}_j,{\mathrm{\Phi }}_k)\hfill \\ \hfill & \equiv {\displaystyle {\int }_{{\mathrm{\Theta }}_i}}\frac{\mathrm{d}{\vartheta }_i{\vartheta }_i}{A_i}{\displaystyle {\int }_{{\mathrm{\Theta }}_j}}\frac{\mathrm{d}{\vartheta }_j{\vartheta }_j}{A_j}{\displaystyle {\int }_{{\mathrm{\Phi }}_k}}\frac{\mathrm{d}\varphi }{|\mathrm{\Delta }{\mathrm{\Phi }}_k|}{\mathrm{\Gamma }}_{\mu }^{\mathcal{P}}({\vartheta }_i,{\vartheta }_j,\varphi ),\hfill \end{array}$$(16)

where for each radial bin, Θ_i, we have ϑ_i ∈ [Θ_low, i, Θ_up, i], and where we defined A_i ≡ (Θ_up, i² − Θ_low, i²)/2. Similarly, for the angular bins we have ϕ ∈ [Φ_low, k, Φ_up, k] and |ΔΦ_k|≡Φ_up, k − Φ_low, k.

3. An efficient estimator of the shear three-point correlation function

Given a weak lensing catalog consisting of N_gal galaxies at positions, θ_i, with ellipticities², γ_c, i, and weights, w_i, the standard estimator of the bin-averaged natural components of the shear 3PCF consists of assigning all galaxy triplets to their corresponding triangle configuration bin and then averaging over those. In particular, the estimator of the zeroth natural component using the × projection Eq. (9) becomes

$$\begin{array}{c}\hfill {\widehat{\overline{\mathrm{\Gamma }}}}_0^{\times }({\mathrm{\Theta }}_1,{\mathrm{\Theta }}_2,\mathrm{\Phi })\equiv \frac{{\mathrm{{\rm Y}}}_0^{\times }({\mathrm{\Theta }}_1,{\mathrm{\Theta }}_2,\mathrm{\Phi })}{\mathcal{N}({\mathrm{\Theta }}_1,{\mathrm{\Theta }}_2,\mathrm{\Phi })},\end{array}$$(17)

where

$$\begin{array}{cc}\hfill {\mathrm{{\rm Y}}}_0^{\times }({\mathrm{\Theta }}_1,{\mathrm{\Theta }}_2,\mathrm{\Phi })& \equiv -{\displaystyle \sum _{i,j,k=1}^{N_{\mathrm{gal}}}}{w}_i{\gamma }_{\mathrm{c},i}{w}_j{\gamma }_{\mathrm{c},j}{w}_k{\gamma }_{\mathrm{c},k}{\mathrm{e}}^{-3\mathrm{i}({\phi }_{\mathit{ij}}+{\phi }_{\mathit{ik}})}\hfill \\ \hfill & \times \mathcal{B}({\theta }_{\mathit{ij}}\in {\mathrm{\Theta }}_1)\mathcal{B}({\theta }_{\mathit{ik}}\in {\mathrm{\Theta }}_2)\mathcal{B}({\varphi }_{\mathit{ijk}}\in \mathrm{\Phi }),\hfill \end{array}$$(18)

$$\begin{array}{cc}\hfill \mathcal{N}({\mathrm{\Theta }}_1,{\mathrm{\Theta }}_2,\mathrm{\Phi })& \equiv {\displaystyle \sum _{i,j,k=1}^{N_{\mathrm{gal}}}}{w}_i{w}_j{w}_k\hfill \\ \hfill & \times \mathcal{B}({\theta }_{\mathit{ij}}\in {\mathrm{\Theta }}_1)\mathcal{B}({\theta }_{\mathit{ik}}\in {\mathrm{\Theta }}_2)\mathcal{B}({\varphi }_{\mathit{ijk}}\in \mathrm{\Phi }).\hfill \end{array}$$(19)

In the preceding equations, we defined θ_ij ≡ |θ_ij|≡|θ_j − θ_i| and ϕ_ijk ≡ φ_ik − φ_ij, with φ_ij denoting the polar angle of θ_ij. We furthermore introduced the bin selection function, ℬ(x ∈ X), which is at unity if x ∈ X and zero otherwise; the different combinations of the bin selection functions then specify the various triangle configurations.

3.1. Multipole decomposition

Let us first consider the numerator in Eq. (17). Instead of a discrete angular binning, one can instead perform a multipole decomposition of Υ₀^×. For an infinitely dense angular binning, one can decompose

$$\begin{array}{c}\hfill {\mathrm{{\rm Y}}}_0^{\times }({\mathrm{\Theta }}_1,{\mathrm{\Theta }}_2,\varphi )\equiv \frac{1}{2\pi }{\displaystyle \sum _{n=-\infty }^{\infty }}{\mathrm{{\rm Y}}}_{0,n}^{\times }({\mathrm{\Theta }}_1,{\mathrm{\Theta }}_2){\mathrm{e}}^{\mathrm{i}n\varphi },\end{array}$$(20)

in which the nth multipole can be successively manipulated to yield

$$\begin{array}{cc}\hfill {\mathrm{{\rm Y}}}_{0,n}^{\times }& ({\mathrm{\Theta }}_1,{\mathrm{\Theta }}_2)\hfill \\ \hfill & \equiv {\displaystyle {\int }_0^{2\pi }}\mathrm{d}\varphi {\mathrm{{\rm Y}}}_0^{\times }({\mathrm{\Theta }}_1,{\mathrm{\Theta }}_2,\varphi ){\mathrm{e}}^{-\mathrm{i}n\varphi }\hfill \\ \hfill & =-{\displaystyle {\int }_0^{2\pi }}\mathrm{d}\varphi {\displaystyle \sum _{i,j,k=1}^{N_{\mathrm{gal}}}}{w}_i{\gamma }_{\mathrm{c},i}{w}_j{\gamma }_{\mathrm{c},j}{w}_k{\gamma }_{\mathrm{c},k}{\mathrm{e}}^{-3\mathrm{i}({\phi }_{\mathit{ij}}+{\phi }_{\mathit{ik}})}\hfill \\ \hfill & \times \mathcal{B}({\theta }_{\mathit{ij}}\in {\mathrm{\Theta }}_1)\mathcal{B}({\theta }_{\mathit{ik}}\in {\mathrm{\Theta }}_2)\mathcal{B}({\varphi }_{\mathit{ijk}}\in \{\varphi \}){\mathrm{e}}^{-\mathrm{i}n\varphi }\hfill \\ \hfill & =-{\displaystyle \sum _{i,j,k=1}^{N_{\mathrm{gal}}}}{w}_i{\gamma }_{\mathrm{c},i}{w}_j{\gamma }_{\mathrm{c},j}{w}_k{\gamma }_{\mathrm{c},k}{\mathrm{e}}^{-3\mathrm{i}({\phi }_{\mathit{ij}}+{\phi }_{\mathit{ik}})}{\mathrm{e}}^{-\mathrm{i}n({\phi }_{\mathit{ik}}-{\phi }_{\mathit{ij}})}\hfill \\ \hfill & \times \mathcal{B}({\theta }_{\mathit{ij}}\in {\mathrm{\Theta }}_1)\mathcal{B}({\theta }_{\mathit{ik}}\in {\mathrm{\Theta }}_2)\hfill \\ \hfill & =-{\displaystyle \sum _{i=1}^{N_{\mathrm{gal}}}}{w}_i{\gamma }_{\mathrm{c},i}{\displaystyle \sum _{j=1}^{N_{\mathrm{gal}}}}{w}_j{\gamma }_{\mathrm{c},j}{\mathrm{e}}^{\mathrm{i}(n-3){\phi }_{\mathit{ij}}}\mathcal{B}({\theta }_{\mathit{ij}}\in {\mathrm{\Theta }}_1)\hfill \\ \hfill & \times {\displaystyle \sum _{k=1}^{N_{\mathrm{gal}}}}{w}_k{\gamma }_{\mathrm{c},k}{\mathrm{e}}^{-\mathrm{i}(n+3){\phi }_{\mathit{ik}}}\mathcal{B}({\theta }_{\mathit{ik}}\in {\mathrm{\Theta }}_2)\hfill \\ \hfill & \equiv -{\displaystyle \sum _{i=1}^{N_{\mathrm{gal}}}}{w}_i{\gamma }_{\mathrm{c},i}{G}_{n-3}^{\mathrm{disc}}({\mathit{\theta }}_i;{\mathrm{\Theta }}_1)G_{-n-3}^{\mathrm{disc}}({\mathit{\theta }}_i;{\mathrm{\Theta }}_2),\hfill \end{array}$$(21)

where in the second step we inserted Eq. (18) and in the final step defined the quantity

$$\begin{array}{cc}\hfill G_n^{\mathrm{disc}}({\mathit{\theta }}_i;\mathrm{\Theta })& ={\displaystyle \sum _{k=1}^{N_{\mathrm{gal}}}}{w}_k{\gamma }_{\mathrm{c},k}{\mathrm{e}}^{\mathrm{i}n{\phi }_{\mathit{ik}}}\mathcal{B}({\theta }_{\mathit{ik}}\in \mathrm{\Theta })\hfill \\ \hfill & \equiv {\displaystyle \sum _{k=1}^{N_{\mathrm{gal}}}}{w}_k{\gamma }_{\mathrm{c},k}{g}_{n;\mathrm{\Theta }}({\mathit{\theta }}_{\mathit{ik}}).\hfill \end{array}$$(22)

Repeating similar calculations for the numerators of the other three natural components, we find

$$\begin{array}{cc}\hfill {\mathrm{{\rm Y}}}_{1,n}^{\times }({\mathrm{\Theta }}_1,{\mathrm{\Theta }}_2)& =-{\displaystyle \sum _{i=1}^{N_{\mathrm{gal}}}}{w}_i{\gamma }_{\mathrm{c},i}^{\ast }{G}_{n-1}^{\mathrm{disc}}({\mathit{\theta }}_i;{\mathrm{\Theta }}_1)G_{-n-1}^{\mathrm{disc}}({\mathit{\theta }}_i;{\mathrm{\Theta }}_2),\hfill \end{array}$$(23)

$$\begin{array}{cc}\hfill {\mathrm{{\rm Y}}}_{2,n}^{\times }({\mathrm{\Theta }}_1,{\mathrm{\Theta }}_2)& =-{\displaystyle \sum _{i=1}^{N_{\mathrm{gal}}}}{w}_i{\gamma }_{\mathrm{c},i}{\left(G_{-n-1}^{\mathrm{disc}}({\mathit{\theta }}_i;{\mathrm{\Theta }}_1)\right)}^{\ast }{G}_{-n-3}^{\mathrm{disc}}({\mathit{\theta }}_i;{\mathrm{\Theta }}_2),\hfill \end{array}$$(24)

$$\begin{array}{cc}\hfill {\mathrm{{\rm Y}}}_{3,n}^{\times }({\mathrm{\Theta }}_1,{\mathrm{\Theta }}_2)& =-{\displaystyle \sum _{i=1}^{N_{\mathrm{gal}}}}{w}_i{\gamma }_{\mathrm{c},i}{G}_{n-3}^{\mathrm{disc}}({\mathit{\theta }}_i;{\mathrm{\Theta }}_1){\left(G_{n-1}^{\mathrm{disc}}({\mathit{\theta }}_i;{\mathrm{\Theta }}_2)\right)}^{\ast },\hfill \end{array}$$(25)

where for the final two equations we used

$$\begin{array}{c}\hfill {\left(G_{-n}^{\mathrm{disc}}({\mathit{\theta }}_i;\mathrm{\Theta })\right)}^{\ast }={\displaystyle \sum _{k=1}^{N_{\mathrm{gal}}}}{w}_k{\gamma }_{\mathrm{c},k}^{\ast }{g}_{n;\mathrm{\Theta }}({\mathit{\theta }}_{\mathit{ik}}).\end{array}$$

For the remainder of this work, we only write down the μ = 0 component of any equation derived from the Υ_μ; the expressions for the other components can be inferred from Eqs. (23)–(25).

Now switching to the denominator, Eq. (19), we can perform a derivation similar to the one leading to Eq. (21) to find the multipole expansion,

$$\begin{array}{c}\hfill \mathcal{N}({\mathrm{\Theta }}_1,{\mathrm{\Theta }}_2,\varphi )\equiv \frac{1}{2\pi }{\displaystyle \sum _{n=-\infty }^{\infty }}{\mathcal{N}}_n({\mathrm{\Theta }}_1,{\mathrm{\Theta }}_2){\mathrm{e}}^{\mathrm{i}n\varphi },\end{array}$$(26)

with multipole components

$$\begin{array}{cc}\hfill {\mathcal{N}}_n({\mathrm{\Theta }}_1,{\mathrm{\Theta }}_2)& ={\displaystyle \sum _{i=1}^{N_{\mathrm{gal}}}}{w}_i{W}_n^{\mathrm{disc}}({\mathit{\theta }}_i;{\mathrm{\Theta }}_1){\left(W_n^{\mathrm{disc}}({\mathit{\theta }}_i;{\mathrm{\Theta }}_2)\right)}^{\ast },\hfill \end{array}$$(27)

$$\begin{array}{cc}\hfill W_n^{\mathrm{disc}}({\mathit{\theta }}_i;\mathrm{\Theta })& \equiv {\displaystyle \sum _{k=1}^{N_{\mathrm{gal}}}}{w}_k{g}_{n;\mathrm{\Theta }}({\mathit{\theta }}_{\mathit{ik}}).\hfill \end{array}$$(28)

The decomposition of the denominator has already been explored in the context of galaxy clustering in which the corresponding 3PCF is a scalar (Slepian & Eisenstein 2015; Philcox et al. 2022). Here we have shown that using the × projection, Eq. (9), gives structurally equivalent expressions for the multipole components of the 3PCF of polar fields.

Regarding the computational complexity, we see that each multipole in the decompositions of the Υ_μ^× and 𝒩 involves a double sum over the whole ellipticity catalog, making the estimation process of the shear 3PCF formally scale as 𝒪(N_gal²). In practice, the largest separation, θ_max, for which we want to compute the 3PCF might be much smaller than the survey extent; in this case, one can employ an efficient search routine that readily finds the various galaxies located in the rings around the θ_i. Assuming a constant runtime³ for this search algorithm, the scaling is brought down to $\mathcal{O}\left(N_{\mathrm{gal}}\overline{n}{\theta }_{\max }^2\right)$.

3.2. Grid-based approximation

While the time complexity of the multipole estimator derived in Sect. 3.1 provides a significant speedup compared to the naive estimator and is competitive with tree-based implementations, the quadratic scaling is still prohibitive for applications in stage-IV surveys. However, we note that most of the complexity stems from the largest separations for which neighboring galaxies make similar contributions to the G_n^disc and W_n^{disc 4}. If one approximates the discrete ellipticity catalog by distributing the w_i and the w_iγ_c, i on a regular spatial grid with pixel size Δ and elements w_i^(Δ) and (wγ_c)_i^(Δ) then for any radial bin with lower bound Θ_low ≫ Δ this approximation will give very similar results to the discrete catalog. The corresponding expression for the G_n and W_n can now be reformulated in terms of a convolution:

$$\begin{array}{cc}\hfill G_n^{\mathrm{\Delta }}({\mathit{\theta }}_i^{(\mathrm{\Delta })};\mathrm{\Theta })& ={\displaystyle \sum _{k=1}^{N_{\mathrm{pix}}^{\mathrm{\Delta }}}}{(w{\gamma }_{\mathrm{c}})}_k^{(\mathrm{\Delta })}{\mathrm{e}}^{\mathrm{i}n{\phi }_{\mathit{ik}}^{(\mathrm{\Delta })}}\mathcal{B}({\theta }_{\mathit{ik}}^{(\mathrm{\Delta })}\in \mathrm{\Theta })\hfill \\ \hfill & =[{(w{\gamma }_{\mathrm{c}})}^{(\mathrm{\Delta })}\star g_{n;\mathrm{\Theta }}^{(\mathrm{\Delta })}]\left({\mathit{\theta }}_i^{(\mathrm{\Delta })}\right),\hfill \end{array}$$(29)

$$\begin{array}{cc}\hfill W_n^{\mathrm{\Delta }}({\mathit{\theta }}_i^{(\mathrm{\Delta })};\mathrm{\Theta })& =[w^{(\mathrm{\Delta })}\star g_{n;\mathrm{\Theta }}^{(\mathrm{\Delta })}]\left({\mathit{\theta }}_i^{\left(\mathrm{\Delta }\right)}\right),\hfill \end{array}$$(30)

where all quantities with the superscript (Δ) are evaluated at the pixel centers and where “⋆” denotes the convolution operator. The transformation to multipole components, Υ_μ, n^×, is equivalent to the discrete case (21), but the outer sum now runs over the number of pixels, $N_{\mathrm{pix}}^{(\mathrm{\Delta })}$, instead of the number of galaxies:

$$\begin{array}{cc}\hfill {\mathrm{{\rm Y}}}_{0,n}^{\times }(& {\mathrm{\Theta }}_1,{\mathrm{\Theta }}_2)\hfill \\ \hfill & =-{\displaystyle \sum _{i=1}^{N_{\mathrm{pix}}^{(\mathrm{\Delta })}}}{\left(w{\gamma }_{\mathrm{c}}\right)}_i^{(\mathrm{\Delta })}{G}_{n-3}^{(\mathrm{\Delta })}({\mathit{\theta }}_i^{(\mathrm{\Delta })};{\mathrm{\Theta }}_1)G_{-n-3}^{(\mathrm{\Delta })}({\mathit{\theta }}_i^{(\mathrm{\Delta })};{\mathrm{\Theta }}_2).\hfill \end{array}$$(31)

The possibility of using fast Fourier transform (FFT) methods for computing the convolution gives a time complexity of the grid-based estimator of $\mathcal{O}(N_{\mathrm{pix}}^{(\mathrm{\Delta })}log\left(N_{\mathrm{pix}}^{(\mathrm{\Delta })}\right))$, significantly outperforming all of the previously mentioned estimators whenever applicable.

3.3. Combined estimator

For an ellipticity catalogue from a cosmic shear survey, we can now devise an efficient estimation procedure of the shear 3PCF, as is depicted in Fig. 2. In the first step, we chose a scale, θ_disc, up to which the multipoles were computed using the discrete prescription outlined in Sect. 3.1. We then picked a base resolution scale, Δ₁ ≡ Δ, and some coarser resolutions, Δ_d ≡ 2^dΔ (d = 2, 3, …), and distributed the ellipticity catalog on the corresponding grids. For each of the scales, θ >  θ_disc, we now selected a resolution, Δ_i(θ), and computed the G_n using the grid-based method from Sect. 3.2. To allocate the multipoles of two scales for which the G_n were computed with different methods, the G_n with higher resolutions needed to be mapped to the coarser grid before the Υ_μ, n^× could be updated. In particular, for the case where one of the G_n was computed via the discrete method and the other one via the grid-based one, we first defined the quantity, $G_n^{\prime (\mathrm{\Delta })}({\mathit{\theta }}_{p_i}^{(\mathrm{\Delta })};{\mathrm{\Theta }}_2)$, as the value of G_n^Δ at the pixel, p, in which the ith galaxy resides, normalized by the number of galaxies in the pth pixel. We now have

Fig. 2. General strategy for computing the three-point correlators used in this work. For small angular scales, we used the prescription for discrete data to compute the Gn and the Υi, n×, while for larger scales, we used the grid-based approximation with varying resolution scales, Δ. For the off-diagonal elements, we first computed the Gn of each radial bin using the corresponding method and then mapped them to the grid of the larger resolution scale from which the Υi, n× were obtained.

$$\begin{array}{cc}\hfill {\mathrm{{\rm Y}}}_{0,n}^{\times }(& {\mathrm{\Theta }}_1,{\mathrm{\Theta }}_2)\hfill \\ \hfill & ={\displaystyle \sum _{i=1}^{N_{\mathrm{gal}}}}{w}_i{\gamma }_{\mathrm{c},i}{G}_{n-3}^{\mathrm{disc}}({\mathit{\theta }}_i;{\mathrm{\Theta }}_1)G_{-n-3}^{\prime (\mathrm{\Delta })}({\mathit{\theta }}_{p_i}^{(\mathrm{\Delta })};{\mathrm{\Theta }}_2)\hfill \\ \hfill & ={\displaystyle \sum _{p=1}^{N_{\mathrm{pix}}^{(\mathrm{\Delta })}}}{G}_{-n-3}^{(\mathrm{\Delta })}({\mathit{\theta }}_p^{(\mathrm{\Delta })};{\mathrm{\Theta }}_2){\displaystyle \sum _{i=1}^{N_{\mathrm{gal},p}^{(\mathrm{\Delta })}}}{w}_i{\gamma }_{\mathrm{c},i}{G}_{n-3}^{\mathrm{disc}}({\mathit{\theta }}_i;{\mathrm{\Theta }}_1)\hfill \\ \hfill & \equiv {\displaystyle \sum _{p=1}^{N_{\mathrm{pix}}^{(\mathrm{\Delta })}}}{\left(w\gamma G_{n-3}^{\mathrm{disc}}\right)}^{(\mathrm{\Delta })}({\mathit{\theta }}_p^{(\mathrm{\Delta })};{\mathrm{\Theta }}_1)G_{-n-3}^{(\mathrm{\Delta })}({\mathit{\theta }}_p^{(\mathrm{\Delta })};{\mathrm{\Theta }}_2),\hfill \end{array}$$(32)

where the second sum in the second line runs over all $N_{\mathrm{gal},p}^{(\mathrm{\Delta })}$ galaxies residing in the pth pixel. This implies that the products, $\mathit{w}\gamma G_{n-3}^{\mathrm{disc}}$, need to be evaluated at the true galaxy positions before they are mapped to the corresponding pixel. In Appendix A, we give more explicit details on the time and space complexity of the combined estimator.

3.4. Edge corrections

Looking back at the general form of the estimator for the shear 3PCF (Eq. 17), we see that for a brute-force-like estimator the shear 3PCF, Γ_μ, can be evaluated by a simple division of the two correlators, Υ_μ and 𝒩, for each angular and polar bin combination. On the other hand, when working on a multipole basis, we have

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_{\mu }^{\times }({\mathrm{\Theta }}_1,{\mathrm{\Theta }}_2,\varphi )=\frac{1}{2\pi }{\displaystyle \sum _{n=-\infty }^{\infty }}{\mathrm{\Gamma }}_{\mu ,n}^{\times }({\mathrm{\Theta }}_1,{\mathrm{\Theta }}_2){\mathrm{e}}^{\mathrm{i}n\varphi }.\end{array}$$(33)

To work out the analog of the division on a multipole basis, we first multiplied Eq. (17) by 𝒩 and then inserted the multipole expansion of all quantities:

$$\begin{array}{c}\hfill {\displaystyle \sum _{n,n\prime =-\infty }^{\infty }}{\overline{\mathrm{\Gamma }}}_{\mu ,n}^{\times }\left({\mathrm{\Theta }}_1,{\mathrm{\Theta }}_2\right)\&{\mathcal{N}}_{n\prime }({\mathrm{\Theta }}_1,{\mathrm{\Theta }}_2){\mathrm{e}}^{\mathrm{i}(n+n\prime )\varphi }\\ \hfill & ={\displaystyle \sum _{n=-\infty }^{\infty }}{\mathrm{{\rm Y}}}_{\mu ,n}^{\times }({\mathrm{\Theta }}_1,{\mathrm{\Theta }}_2){\mathrm{e}}^{\mathrm{i}n\varphi }.\hfill \end{array}$$(34)

Multiplying Eq. (34) by e^−iℓϕ and integrating over ϕ, we find that

$$\begin{array}{c}\hfill {\displaystyle \sum _{n=-\infty }^{\infty }}{\mathcal{N}}_{\ell -n}({\mathrm{\Theta }}_1,{\mathrm{\Theta }}_2){\overline{\mathrm{\Gamma }}}_{\mu ,n}^{\times }\left({\mathrm{\Theta }}_1,{\mathrm{\Theta }}_2\right)={\mathrm{{\rm Y}}}_{\mu ,\ell }^{\times }({\mathrm{\Theta }}_1,{\mathrm{\Theta }}_2).\end{array}$$(35)

Thus, we see that in general there exists a coupling between different multipoles. While Eq. (35) formally describes an infinite system of equations, one expects 𝒩_{ℓ − n} to be dominated by its diagonal, such that truncating the expansion at some n_max will provide a sufficient approximation⁵. Inverting the resulting (2n_max + 1) – dimensional system equips us with an expression for the edge-corrected multipoles of the shear 3PCF:

$$\begin{array}{c}\hfill {\overline{\mathrm{\Gamma }}}_{\mu ,n}^{\times }\left({\mathrm{\Theta }}_1,{\mathrm{\Theta }}_2\right)={\displaystyle \sum _{\ell =-n_{\max }}^{n_{\max }}}{\left({\mathcal{C}}^{-1}\right)}_{n\ell }({\mathrm{\Theta }}_1,{\mathrm{\Theta }}_2)\frac{{\mathrm{{\rm Y}}}_{\mu ,\ell }^{\times }({\mathrm{\Theta }}_1,{\mathrm{\Theta }}_2)}{{\mathcal{N}}_0({\mathrm{\Theta }}_1,{\mathrm{\Theta }}_2)},\end{array}$$(36)

where we defined the components of the mode coupling matrix, 𝒞, as

$$\begin{array}{c}\hfill {\mathcal{C}}_{\ell n}({\mathrm{\Theta }}_1,{\mathrm{\Theta }}_2)=\frac{{\mathcal{N}}_{\ell -n}({\mathrm{\Theta }}_1,{\mathrm{\Theta }}_2)}{{\mathcal{N}}_0({\mathrm{\Theta }}_1,{\mathrm{\Theta }}_2)}.\end{array}$$(37)

For a uniform distribution of tracers, 𝒩 does not carry any angular dependence such that 𝒞 becomes the identity matrix. We note that in the cosmological context, a nonuniform distribution of tracers can either be realized by a clustering of the tracers or a nontrivial survey geometry. We further note that this mode coupling effect has already been explored in Slepian & Eisenstein (2015) in the context of edge-correcting the 3PCF of scalar fields.

3.5. Including tomography

Our estimator can straightforwardly be generalized to a tomographic setting with n_z tomographic redshift bins {Z_k} (k ∈ {1, ⋯, n_z}), where we defined the galaxy at location X_i in each galaxy triplet to lie in the Z_ith redshift bin. In the case of gridded data, we simply need to evaluate the G_n^Δ for each tomographic bin and then allocate the Υ_μ, n^× for all relevant combinations:

$$\begin{array}{cc}& {\mathrm{{\rm Y}}}_{0,n}^{\times }({\mathrm{\Theta }}_1,{\mathrm{\Theta }}_2;Z_1,Z_2,Z_3)=-{\displaystyle \sum _{i=1}^{N_{\mathrm{pix}}^{(\mathrm{\Delta })}}}{\left(w_i{\gamma }_{\mathrm{c},i}\right)}^{(Z_1,\mathrm{\Delta })}\hfill \\ \hfill & \times G_{n-3}^{(Z_2,\mathrm{\Delta })}({\mathit{\theta }}_i;{\mathrm{\Theta }}_1)G_{-n-3}^{(Z_3,\mathrm{\Delta })}({\mathit{\theta }}_i;{\mathrm{\Theta }}_2),\hfill \end{array}$$(38)

$$\begin{array}{cc}& G_n^{(Z_2,\mathrm{\Delta })}({\mathit{\theta }}_i;{\mathrm{\Theta }}_1)=[{(w{\gamma }_{\mathrm{c}})}^{(Z_2,\mathrm{\Delta })}\star g_{n;\mathrm{\Theta }}^{(\mathrm{\Delta })}]({\mathit{\theta }}_i).\hfill \end{array}$$(39)

For discrete data, the modifications are similar, but one additionally needs to keep track of the redshift bins of both galaxies contributing to the G_n:

$$\begin{array}{cc}\hfill {\mathrm{{\rm Y}}}_{0,n}^{\times }({\mathrm{\Theta }}_1,{\mathrm{\Theta }}_2;Z_1,Z_2,Z_3)& =-{\displaystyle \sum _{i=1}^{N_{\mathrm{gal}}^{(Z_1)}}}{\left(w_i{\gamma }_{\mathrm{c},i}\right)}^{(Z_1,\mathrm{\Delta })}\hfill \\ \hfill \times G_{n-3}^{(Z_2,\mathrm{disc})}& ({\mathit{\theta }}_i;{\mathrm{\Theta }}_1)G_{-n-3}^{(Z_3,\mathrm{disc})}({\mathit{\theta }}_i;{\mathrm{\Theta }}_2),\hfill \end{array}$$(40)

$$\begin{array}{cc}\hfill G_n^{(Z_2,\mathrm{disc})}({\mathit{\theta }}_i;\mathrm{\Theta })& \equiv {\displaystyle \sum _{k=1}^{N_{\mathrm{gal}}^{(Z_2)}}}{w}_k{\gamma }_{\mathrm{c},k}{g}_{n;\mathrm{\Theta }}({\mathit{\theta }}_{\mathit{ik}}).\hfill \end{array}$$(41)

The transformation to the natural components works analogously to Eq. (36).

4. Validation of the estimator

We validated our estimator on the Scinet-LIghtCone Simulations (SLICS) simulation suite (Harnois-Déraps et al. 2018). It consists of a suite of dark-matter-only N-body simulations that track the evolution of 1536³ particles within a box of side length 505 h⁻¹Mpc. From each simulation, convergence maps are constructed up until z = 3 within a field of view of 100 deg² on a grid of 7745² pixels. For this work, we used the KiDS-450-like ensemble, in which 819 fully independent ellipticity catalogs mimicking the redshift distribution and number density of the KiDS-450 data (de Jong et al. 2017) were constructed from past light cones. We further used a maximum multipole of n_max ≡ 20 for all tests presented in this and in the following section. For the remainder of this work, when transforming the shear 3PCF from its multipole basis to the real space basis, we use 100 linearly spaced angular bins between [0, 2π], independent of the chosen value of n_max.

4.1. Accuracy of the combined estimator

To test the applicability of the grid-based approximation and the associated combined estimator, we applied all of those schemes to a single noiseless mock catalog of the SLICS ensemble. In particular, we evaluated the 3PCF multipoles in 33 log-spaced radial bins between $0\stackrel{\prime }{.}25$ and 80′, which were then transformed to the real-space shear 3PCF. In Fig. 3, we compared the runs using the discrete estimator, the grid-based approximation for resolution scales $\mathrm{\Delta }\in \{0\stackrel{\prime }{.}25,0\stackrel{\prime }{.}5,1\stackrel{\prime }{.}0\}$, as well as the combined estimator that uses a grid-based approximation for scales 20Δ ≤ Θ_low ≤ Θ_up ≤ 40Δ and afterward switches to the next resolution scale up until Δ = 1′⁶. For this binning scheme, the combined estimator performed ≈20 times faster than the discrete estimator, even outperforming the pure grid-based approximation of $\mathrm{\Delta }=0\stackrel{\prime }{.}25$. Focusing our attention on the equilateral configurations, we see that the grid-based approximation is in good agreement with the discrete estimator once the separation, Θ, is much larger than the pixelization scale, Δ, but it starts to bias the statistics once the separation is less than around ten times the pixelization scale. For the equilateral configurations, the combined estimator is, by construction, equal to the discrete estimator for Θ <  5′, and it then follows the measurements of the grid-based implementations. For our example shown in Fig. 3, the combined estimator is in agreement with the discrete estimator at a sub-percent level.

Fig. 3. Performance of the combined estimator on a noiseless mock catalog with ≈3.1 × 106 ellipticities on a 100 deg2 area. In the upper panels, solid lines indicate the real part of the 3PCF, while dashed lines represent the imaginary part. Left-hand side: Comparison of various discrete, grid-based estimators, and a combined implementation of the multipole estimator for equilateral configurations of the zeroth natural component of the shear 3PCF. The timings correspond to the runtime on eight CPU cores and the dashed vertical lines indicate the regions of constant grid resolution for the combined estimator. In the bottom panel, we show the ratios of the real parts between all approximate implementations and the discrete estimator. The gray-shaded region displays a one-percent interval. We note that, by construction, the curve corresponding to the combined estimator always lies on top of the curve corresponding to the resolution, Δd, in the interval 20Δd ≤ Θlow ≤ Θup ≤ 40Δd. Right-hand side: Comparison of the discrete and the combined estimator for different triangle shapes. We fix two triangle sides and vary the enclosing angle, ϕ. In the bottom panel, the thick red line shows the difference between the real part of both estimators, normalized by the largest value of the discrete estimator for the given radial bin combination, maxD. The gray lines display the corresponding ratios of all other radial bin combinations, where again the normalization was done with respect to the largest value of the discrete estimator for the corresponding radial bin combination.

To assess the accuracy of the blocks in the (Θ₁, Θ₂) plane shown in Fig. 2, in which the first scale was evaluated using the discrete prescription, while the second scale uses the grid-based approximation, we picked one such combination of radial bins and then varied the enclosed angle, ϕ, between them. For the example shown in the right-hand panel in Fig. 3, we chose a scale for which the G_n^(Δ) are evaluated on a grid with Δ = 1′ and we find an excellent agreement between the two methods. We also show the ratios between the two schemes for all other bin combinations. Although for the majority of the curves the ratio is very close to unity, there do appear to be some spikes at which the relative difference between both estimators is at the five percent level. However, those bin configurations correspond to cases in which are only a few triplets present, and therefore they will not carry significant signal-to-noise (S/N) in realistic data with non-vanishing shape noise.

4.2. Comparison with TREECORR

The state-of-the-art method of computing shear correlation functions is based on tree codes – the most well-known implementation is the publicly available code TREECORR (Jarvis et al. 2004)⁷. In it, the ${\overline{\mathrm{\Gamma }}}_{\mu }$ are estimated by explicitly accounting for all triplets. For an efficient evaluation, TREECORR first organizes the data into a hierarchical ball-tree structure from which the 3PCF is then computed. A further speedup can be achieved when allowing the bin edges to be fuzzy. For reasonable values of the associated bin_slop parameter, this approximation does not bias the mean, but does increase the measurement uncertainty (Secco et al. 2022). For an efficient binning of triangle shapes, TREECORR uses the quantities (r, u, v), which are related to the three triangle side lengths, d₁ ≥ d₂ ≥ d₃, as

$$\begin{array}{c}\hfill r=d_2;u=\frac{d_3}{d_2};v=\pm \frac{d_1-d_2}{d_3},\end{array}$$(42)

such that with u ∈ [0, 1], v ∈ [ − 1, 1] all triangle configurations with d₁ ≥ d₂ ≥ d₃ are covered. As is shown in Jarvis et al. (2004), this coverage is sufficient to obtain integrated measures of the 3PCF, such as the third-order aperture mass.

In Fig. 4, we compare the estimated components of the first two natural components, ${\mathrm{\Gamma }}_{0,1}^{\mathrm{cent}}$, in a single noiseless mock catalog of the SLICS ensemble. For this, we first transformed the bin centers of the (Θ₁, Θ₂, ϕ) basis in the (r, u, v) basis and chose the value of the TREECORR measurement within the bin that the transformed bin center fell into. In order to facilitate an accurate matching between the two binning parameterizations, one needs to use very fine bins for both estimators. For our test, we computed the 3PCF between 5′ and 50′, where for the multipole-based estimator we used 40 logarithmically spaced radial bins without making use of the grid-based approximation. For the TREECORR measurement, we chose the same r binning and used 20 (40) linearly spaced bins covering the maximum range of the u (v) parameters. To achieve a reasonable run time, we only selected every tenth galaxy, leaving a total of ≈3.1 × 10⁵ ellipticities. The associated measurements were carried out on 16 CPU cores and took ≈1.5 minutes for the discrete multipole-based estimator and ≈164 minutes for TREECORR when using a value of 0.8 for the bin_slop parameter. In Fig. 4, we find good agreement overall between both estimators across the different triangle configurations; the only significant differences are visible for strongly squeezed and folded triangles. We note that the small-scale wiggles visible in the TREECORR measurements are due to the different binning parameterizations and do not reflect an additional noise component in the estimator (see Appendix B for further details).

Fig. 4. First two natural components of the shear 3PCF on a catalog with ≈3.1 × 105 shapes on a 100 deg2 area. We fixed two triangle sides and varied the enclosing angle, ϕ. Solid lines indicate the real part of the 3PCF, while dashed lines stand for the imaginary part. The regions between the dashed black lines correspond to the intervals in which the ordering of the lengths of the triangle is constant. We note that due to the different binning parameterizations of TREECORR and our estimator, we do not expect perfect agreement between the curves.

Besides the significantly lower computational cost compared to the tree code, the multipole-based estimator also automatically allows for a fine ϕ binning of the shear 3PCF, and thus provides a good approximation even for strongly degenerate triangles⁸. Finally, we note that the relative speedup with respect to TREECORR increases when the combined estimator is employed and when the number density of galaxies increases. For a thorough scaling analysis of the combined estimator, we again refer to Appendix A.

5. Application to third-order aperture mass measures

Due to the large size of a data vector built from an estimate of the shear 3PCF, it is unfeasible to perform a cosmological analysis using this probe. The preferred statistic is the third-order aperture mass, which provides both a compression of the data vector and a separation of the third-order shear signal in its E- and B-modes.

5.1. Aperture mass

The aperture mass introduced in Kaiser (1995), Schneider (1996) is a measure for the projected overdensity within a circular region of radius θ around a particular location, ϑ:

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

where the filter function, U_θ, is compensated, meaning that aperture mass measures evade the mass-sheet degeneracy (Falco et al. 1985). Within the flat-sky approximation, the aperture mass can also be written as the real part of a complex aperture measure, being defined in terms of the shear,

$$\begin{array}{cc}\hfill \mathcal{M}(\mathit{\vartheta };\theta )& ={\mathcal{M}}_{\mathrm{ap}}(\mathit{\vartheta };\theta )+\mathrm{i}{\mathcal{M}}_{\times }(\mathit{\vartheta };\theta )\hfill \\ \hfill & ={\displaystyle \int {\mathrm{d}}^2}\mathit{\vartheta }\prime Q_{\theta }(|\mathit{\vartheta }\prime |)\gamma (\mathit{\vartheta }+\mathit{\vartheta }\prime ;\phi \prime ),\hfill \end{array}$$(44)

in which Q_θ is another filter function uniquely related to the choice of the U_θ filter, and φ′ denotes the polar angle of the separation vector, ϑ′. With the aperture measure being able to separate the shear signal into its E- and B-modes, it can be used for both gaining cosmological information and pinning down potential systematics. For this work, we used the exponential filter introduced in Crittenden et al. (2002),

$$\begin{array}{c}\hfill Q_{\theta }(\vartheta )=\frac{{\vartheta }^2}{4\pi {\theta }^4}exp\left[(-\frac{{\vartheta }^2}{2{\theta }^2})\right].\end{array}$$(45)

5.2. Third-order aperture mass measures

The E/B decomposition also extends to the third-order aperture mass measures. For example, combining three aperture measures as

$$\begin{array}{c}\hfill \langle {\mathcal{M}}_{\mathrm{ap}}^3\rangle ({\theta }_1,{\theta }_2,{\theta }_3)\equiv \langle {\mathcal{M}}_{\mathrm{ap}}({\theta }_1){\mathcal{M}}_{\mathrm{ap}}({\theta }_2){\mathcal{M}}_{\mathrm{ap}}({\theta }_3)\rangle ,\end{array}$$(46)

one can see from Eq. (43) that the resulting quantity depends on the third-order correlation function of the convergence field. By writing down the nonredundant products that appear when combining three complex measures, ℳ, one obtains four real-valued quantities ⟨ℳ_ap³⟩, ⟨ℳ_ap²ℳ_×⟩, ⟨ℳ_apℳ_ײ⟩, and ⟨ℳ_׳⟩, which are all related to a filtered version of the third-order shear correlator. Of those, only ⟨ℳ_ap³⟩ carries cosmological information, while the two measures, ⟨ℳ_ap²ℳ_×⟩ and ⟨ℳ_׳⟩, vanish for parity symmetric shear fields (Schneider 2003), and a statistically significant signal of ⟨ℳ_apℳ_ײ⟩ indicates the presence of B-modes (Schneider et al. 2005, hereafter SKL05). Using the centroid projection Eq. (10) for the natural components, SKL05 derived an explicit form of the transformation between the Γ_μ^cen and the third-order aperture mass measures that are defined in terms of two third-order correlators of the complex aperture measure, ℳ:

$$\begin{array}{cc}\hfill \langle {\mathcal{M}}^3\rangle ({\theta }_1,{\theta }_2,{\theta }_3)& ={\displaystyle {\int }_0^{\infty }}\mathrm{d}{\vartheta }_1{\displaystyle {\int }_0^{\infty }}\mathrm{d}{\vartheta }_2{\displaystyle {\int }_0^{2\pi }}\mathrm{d}\varphi \hfill \\ \hfill & \times {\mathrm{\Gamma }}_0^{\mathrm{cen}}({\vartheta }_1,{\vartheta }_2,\varphi )F_0({\theta }_1,{\theta }_2,{\theta }_3;{\vartheta }_1,{\vartheta }_2,\varphi ),\hfill \\ \hfill \langle {\mathcal{M}}^{\ast }{\mathcal{M}}^2\rangle ({\theta }_1,{\theta }_2,{\theta }_3)& ={\displaystyle {\int }_0^{\infty }}\mathrm{d}{\vartheta }_1{\displaystyle {\int }_0^{\infty }}\mathrm{d}{\vartheta }_2{\displaystyle {\int }_0^{2\pi }}\mathrm{d}\varphi \hfill \\ \hfill & \times {\mathrm{\Gamma }}_1^{\mathrm{cen}}({\vartheta }_1,{\vartheta }_2,\varphi )F_1({\theta }_1,{\theta }_2,{\theta }_3;{\vartheta }_1,{\vartheta }_2,\varphi ),\hfill \end{array}$$(47)

where the filter functions, F_0, 1, in our convention only differ from the ones in SKL05 by cyclic permutations of the indices⁹. We give the explicit expressions of the filters in our convention in Appendix D.

While these two correlators are sufficient for the theoretical conversion, this might not be the case when the shear 3PCF is estimated from noisy data that is distributed over several tomographic redshift bins. In this case, one needs to include two additional correlators if the full S/N has to be retrieved:

$$\begin{array}{cc}\hfill \langle \mathcal{M}{\mathcal{M}}^{\ast }\mathcal{M}\rangle ({\theta }_1,{\theta }_2,{\theta }_3)& ={\displaystyle {\int }_0^{\infty }}\mathrm{d}{\vartheta }_1{\displaystyle {\int }_0^{\infty }}\mathrm{d}{\vartheta }_2{\displaystyle {\int }_0^{2\pi }}\mathrm{d}\varphi \hfill \\ \hfill & \times {\mathrm{\Gamma }}_2^{\mathrm{cen}}({\vartheta }_1,{\vartheta }_2,\varphi )F_2({\theta }_1,{\theta }_2,{\theta }_3;{\vartheta }_1,{\vartheta }_2,\varphi ),\hfill \\ \hfill \langle {\mathcal{M}}^2{\mathcal{M}}^{\ast }\rangle ({\theta }_1,{\theta }_2,{\theta }_3)& ={\displaystyle {\int }_0^{\infty }}\mathrm{d}{\vartheta }_1{\displaystyle {\int }_0^{\infty }}\mathrm{d}{\vartheta }_2{\displaystyle {\int }_0^{2\pi }}\mathrm{d}\varphi \hfill \\ \hfill & \times {\mathrm{\Gamma }}_3^{\mathrm{cen}}({\vartheta }_1,{\vartheta }_2,\varphi )F_3({\theta }_1,{\theta }_2,{\theta }_3;{\vartheta }_1,{\vartheta }_2,\varphi ).\hfill \end{array}$$(48)

The individual third-order aperture mass measures then follow by a linear combination of the previous quantities; in particular, for ⟨ℳ_ap³⟩, one has

$$\begin{array}{cc}& \langle {\mathcal{M}}_{\mathrm{ap}}^3\rangle ({\theta }_1,{\theta }_2,{\theta }_3)\hfill \\ \hfill & =\mathfrak{R}[\langle {\mathcal{M}}^2{\mathcal{M}}^{\ast }\rangle ({\theta }_1,{\theta }_2,{\theta }_3)+\langle \mathcal{M}{\mathcal{M}}^{\ast }\mathcal{M}\rangle ({\theta }_1,{\theta }_2,{\theta }_3)\hfill \\ \hfill & +\langle {\mathcal{M}}^{\ast }{\mathcal{M}}^2\rangle ({\theta }_1,{\theta }_2,{\theta }_3)+\langle {\mathcal{M}}^3\rangle ({\theta }_1,{\theta }_2,{\theta }_3)]/4,\hfill \end{array}$$(49)

and the linear combinations of the remaining aperture measures can be constructed in a similar fashion.

For a tomographic setup, the structure of the equations above remains equal, but we make the substitution

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_{\mu }^{\mathrm{cent}}({\vartheta }_1,{\vartheta }_2,\varphi )\to {\mathrm{\Gamma }}_{\mu }^{\mathrm{cent}}({\vartheta }_1,{\vartheta }_2,\varphi ;Z_1,Z_2,Z_3),\end{array}$$(50)

where the Γ_μ^cent(ϑ₁, ϑ₂, ϕ; Z₁, Z₂, Z₃) are estimated using the method described in Sect. 3.5. This in turn means that all the other quantities, like ⟨ℳ³⟩ or ⟨ℳ_ap³⟩, also become functions of the tomographic redshift bin combination (Z₁, Z₂, Z₃). We note that while the third-order aperture mass is symmetric with respect to permutations in Z_i, this is not the case for the shear 3PCF, and therefore one needs to average over all possible Z_i permutations to obtain the full S/N. We again refer to Appendix D for all explicit expressions.

5.3. Estimators of third-order aperture measures

5.3.1. Estimation via shear 3PCF

As was established in the previous subsection, one can estimate the third-order aperture mass measures by integrating over the 3PCF on a centroid basis. To do this in practice, we first used the multipole-based estimator (Eq. 36) introduced in Sect. 3 to obtain the shear 3PCF in the × projection, then transformed them to a centroid basis via Eqs. (12–15), and finally performed the numerical integrals contained in Eq. (49). In Appendix F.2, we show the level of accuracy to which the 3PCF needs to be computed in order to guarantee a stable integration. As the computational bottleneck of this pipeline is the estimation of the shear 3PCF, the same computational speedup applies to the estimation of third-order aperture measures when using the multipole-based estimator instead of traditional estimation methods of the shear 3PCF.

5.3.2. Estimation via a direct estimator

The original estimator of the third-order aperture mass statistics in simulated data on an area, A, was introduced in Schneider et al. (1998). Instead of evaluating the statistics in terms of the shear 3PCF, it directly samples a third-order generalization of Eq. (44),

$$\begin{array}{cc}\hfill \langle {\mathcal{M}}_{\mathrm{ap}}^3\rangle ({\theta }_1,{\theta }_2,{\theta }_3)& =\int \frac{{\mathrm{d}}^2\mathit{\vartheta }}{A}\hfill \\ \hfill & \times {\displaystyle \prod _{i=1}^3}\hfill & \hfill \int {\mathrm{d}}^2{\mathit{\vartheta }}_i{Q}_{{\theta }_i}\left(|{\mathit{\vartheta }}_i|\right){\gamma }_{\mathrm{t}}(\mathit{\vartheta }+{\mathit{\vartheta }}_i;{\varphi }_i),\end{array}$$(51)

from the ellipticities at the galaxy positions. In other words, at a particular angular position, ϑ, one has

$$\begin{array}{cc}\hfill {\widehat{\mathcal{M}}}_{\mathrm{ap}}^3(& \mathit{\vartheta };{\theta }_1,{\theta }_2,{\theta }_3)\hfill \\ \hfill & =\frac{{\sum }_{i\ne j\ne k}^{N_{\mathrm{gal}}}{w}_i{w}_j{w}_k{Q}_{{\theta }_1,i}{Q}_{{\theta }_2,j}{Q}_{{\theta }_3,k}{\gamma }_{\mathrm{t},i}{\gamma }_{\mathrm{t},j}{\gamma }_{\mathrm{t},k}}{{\displaystyle {\sum }_{i\ne j\ne k}^{N_{\mathrm{gal}}}}{w}_i{w}_j{w}_k{Q}_{{\theta }_1,i}{Q}_{{\theta }_2,j}{Q}_{{\theta }_3,k}},\hfill \end{array}$$(52)

where we abbreviated Q_θl, i ≡ Q_θl(|ϑ_i − ϑ|) and γ_t, i ≡ γ_t(ϑ_i; ζ_i), with ϑ_i being the angular position of the ith galaxy and ζ_i denoting the direction, ϑ_i − ϑ. The corresponding estimate for ⟨ℳ_ap³⟩ is then obtained by covering the survey footprint with apertures and averaging over the individual estimates. As is shown in Porth & Smith (2021), the nested sums appearing in Eq. (52) can be decomposed such that the estimation procedure scales linearly with the number of galaxies in the survey¹⁰. Although the employed Q filter (45) does formally have infinite support, one can show that by only including galaxies separated by at most 4θ from the aperture center, one recovers a practically unbiased result (Heydenreich et al. 2023). Thus, in order to avoid border effects, we do not sample apertures in regions that are separated by less than 4max(θ₁, θ₂, θ₃) from the fields’ boundary. We recall that while the direct estimator can be used on unmasked simulated data, like the SLICS ensemble, this is not the case for real data with a nontrivial angular survey mask.

5.4. Comparison of both estimators

We compared the results of both estimation procedures on the SLICS ensemble, using all 819 mock catalogs. We first estimated the 3PCF in the interval $[0\stackrel{\prime }{.}3125,160\prime ]$ using a logarithmic bin width of $b\equiv \frac{{\mathrm{\Theta }}_{\mathrm{up},\mathrm{i}+1}}{{\mathrm{\Theta }}_{\mathrm{up},\mathrm{i}}}=1.05$. We then transformed the 3PCF to the equal-scale aperture statistics, where for the aperture radii we chose 30 logarithmically spaced values between [1′,25′]. Finally, we estimated ⟨ℳ_ap³⟩ using the direct estimator (52) for the same set of aperture radii. In Fig. 5, we show the measurements obtained from both estimators, together with their standard deviation across the SLICS ensemble. We find a percent-level agreement on aperture scales between 2′ and 10′. While for smaller aperture radii, the multipole-based estimator begins to yield systematically lower results, for larger aperture radii there appears to be a scatter around the zero axis. Both of these cases were expected. On the one hand, as the shear 3PCF has not been estimated down to zero separation, the integral transformation (47) cannot be fully sampled for small aperture radii. Shi et al. (2014) have demonstrated that this results in a systematic underprediction of the estimated ⟨ℳ_ap³⟩ for aperture radii of less than ∼5 times the smallest angular separation considered in the 3PCF. At large scales, on the other hand, one would not expect both estimators to perfectly agree with each other. Cutting off the survey boundary for the direct estimator implies that both estimators do not use the data in the same way. This can also be seen in the standard deviation of both estimators across the ensemble. While for small aperture radii the error bars are of a similar size, the multipole-based estimator gradually extracts more S/N compared to the direct estimator for increasing aperture radii, as not all triangles can be accounted for in the latter due to the restricted sampling domain of the apertures.

Fig. 5. Comparison of the third-order aperture mass statistics between the multipole-based estimator (blue line with error bars) and the direct estimator (dashed red line with error band) in the SLICS ensemble. In the lower panel, we plot the ratio between the two measurements (blue line) and show a one percent interval as the gray-shaded region.

6. Application to the KiDS-1000 data

6.1. The KiDS-1000 data

For the first application of our estimator to real data, we used the gold sample of weak lensing and photometric redshift measurements from the fourth data release of the Kilo-Degree Survey (Kuijken et al. 2019; Wright et al. 2020; Hildebrandt et al. 2021; Giblin et al. 2021)¹¹. The gold sample consists of around 21 million galaxies distributed over an area of 1006 square degrees (777 square degrees after masking) and is further divided into five tomographic redshift bins.

6.2. Covariance matrix

We obtain an estimate for the covariance matrix by using a suite of full-sky gravitational lensing simulations introduced in Takahashi et al. (2017) (hereafter the T17 ensemble). The underlying dark-matter-only N-body simulations were run on a ΛCDM cosmology with Ω_m = 0.279, Ω_Λ = 1 − Ω_m, h = 0.7, σ₈ = 0.82, and spectral index n_s = 0.97. We constructed 1944 mock ellipticity catalogs from the T17 ensemble. In particular, we matched the angular positions, as well as the weights, the shape noise level per tomographic redshift bin, and the correlation between the measured ellipticities and the weights to the public KiDS-1000 ellipticity catalog (Giblin et al. 2021). We note that those specifications give rise to a different correlation structure than would be given by the SLICS ensemble, in which the galaxies are randomly placed within an unmasked area and in which each ellipticity carries an equal weight. For further details about the creation of the mocks, we refer to Burger et al. (2024).

To measure the shear 3PCF via the multipole estimator we first decomposed the survey area into 16 smaller patches. We proceeded along the lines of Appendix C to measure the 3PCF multipoles in each patch and then joined them, after which we obtained an unbiased estimate for the third-order aperture statistics. We measured the shear 3PCF in 38 angular bins in the interval $[0\stackrel{\prime }{.}75,240\prime ]$, yielding a mean logarithmic angular bin width of b ≈ 1.15. The choice of the lower bound was due to the pixelization scale of underlying shear maps in the T17 ensemble, which are given in a Healpix format with nside=8192, for which the individual pixels have a side length of $\approx 0\stackrel{\prime }{.}4$. As we show in Appendix F.2, this range of angular scales yields a numerically stable transformation of the shear 3PCF to the aperture mass measures for aperture scales between 4′ and 40′. We do not expect much additional information to be contained on larger aperture scales, as the third-order aperture mass only filters out non-Gaussian contributions to the cosmic shear field. As we aim to obtain an accurate covariance matrix for the unequal-scale aperture statistics in the non-tomographic setup, we computed the multipoles up to the thirtieth order in this case. In the tomographic case, we only considered equal aperture scales in the measurement, for which we show in Appendix F.2 that choosing n_max = 10 is sufficient. Finally, we transformed the shear 3PCF to the aperture statistics on two sets of radii, the “fine-binning set” consisting of 30 logarithmically spaced radii in the interval [4′,40′] and the “coarse-binning set” for which θ ∈ {4′,6′,8′,10′,14′,18′,22′,26′,32′,36′}. Using these configurations, the estimation of the third-order aperture measures per mock takes around 14 hours on a single CPU in the non-tomographic case, while 4.5 hours are required on 13 cores when including all 125 tomographic bin combinations for the 3PCF.

We produced mock covariance matrices for both a tomographic and a non-tomographic setup and display the resulting correlation matrices in Fig. 6. In particular, we see that the third-order aperture mass is highly correlated for radial bins of a similar angular scale. Similarly, there is a significant correlation between aperture measures being computed from different tomographic bin combinations. We used these sample correlation matrices to assess the detection significance of the third-order aperture measures in the KiDS-1000 data.

Fig. 6. Correlation matrices of third-order aperture mass measures in the T17 ensemble. Shown are the results for the non-tomographic, equal-scale statistics (top left), the non-tomographic unequal-scale statistics (bottom left), and the tomographic equal-scale statistics (right-hand side). We order the statistics such that θi ≤ θj ≤ θk and Z1 ≤ Z2 ≤ Z3, meaning that each of the “small” 352 squares in the covariance matrix on the right-hand side corresponds to the cross-covariance of the statistics for fixed tomographic redshift bin combinations.

Besides the correlation structure, we have validated that there is no significant presence of B-modes in the mock catalogs and that the measurement of the E-mode matches the theoretical model introduced in Heydenreich et al. (2023). We refer to Appendix E for the corresponding figure.

6.3. Measurement results

In this subsection, we present our measurement of the third-order aperture measures in the KiDS-1000 data. Throughout this subsection, we make repeated use of the “chi-square”

$$\begin{array}{c}\hfill {\chi }^2={\left[\mathbf{m}(\mathrm{\Theta })-\mathbf{d}\right]}^T{\mathrm{Cov}}^{-1}\left[\mathbf{m}(\mathrm{\Theta })-\mathbf{d}\right],\end{array}$$(53)

in which d denotes the data vector, Cov stands for the covariance matrix of d, and m(Θ) is a model vector described by a set of parameters, Θ. As the covariance matrix is itself estimated from the T17 ensemble, we needed to include the Hartlap factor (Anderson 2003; Hartlap et al. 2007) to obtain an unbiased estimate for the precision matrix, Cov⁻¹. We quote the detection significance of d by virtue of its p-value associated with the null hypothesis of a vanishing signal, m(Θ) ≡ 0. Unless otherwise stated, we estimated the p-value of a statistic measured in the KiDS-1000 data by first computing the χ² for each of the 1944 mocks and then defining p as the fraction of mocks that have a χ² smaller than the χ² in the real data. Finally, we introduced the S/N of a data vector by following the definition of Secco et al. (2022):

$$\begin{array}{c}\hfill \mathrm{S}/\mathrm{N}\equiv \{\begin{array}{cc}\sqrt{{\chi }^2-N_{\mathrm{d}.\mathrm{o}.\mathrm{f}.}}\hfill & \mathrm{if}{\chi }^2>N_{\mathrm{d}.\mathrm{o}.\mathrm{f}.}+1\hfill \\ \mathrm{Null}\hfill & \mathrm{else}\hfill \end{array},\end{array}$$(54)

where N_d.o.f. denotes the degrees of freedom, which for our setup is equal to the dimension of the data vector, d. For a thorough cosmological analysis of the measurements, we refer to Burger et al. (2024).

6.3.1. Non-tomographic measurement

We estimated the non-tomographic third-order aperture measures in the KiDS-1000 data with the same angular binning that was used to obtain the corresponding covariance matrix in the T17 ensemble. While the figures were done using the fine-binning set, all statistical tests were performed on the coarsely binned set of aperture radii. We further repeated the measurement using various different setups like a finer angular binning with the combined estimator, using only the discrete estimator, excluding the patch-overlap and a run using TREECORR. In all cases, we find that the measurements scatter around each other at a level of ∼10% of the surveys’ error budget at scales ≲20′, while for larger scales the scatter increases a bit between the estimates, utilizing the overlap between footprints and the estimators that do not. All conclusions drawn from each of the setups are the same as the ones presented below.

In Fig. 7, we show our measurement of the non-tomographic third-order aperture measures, out of which only ⟨ℳ_ap³⟩ seems to be significantly different from zero. To test the null hypothesis “no cosmological ⟨ℳ_ap³⟩ signal”, we used the p-value introduced above but calculated the χ² using a covariance matrix that is consistent with a noise-only signal (i.e., the ⟨ℳ_׳⟩ covariance matrix). In this case, we computed the p-value using the percentile-point function of the χ² distribution, which yields a value of p = 5.7 × 10⁻²⁸, strongly rejecting the null hypothesis. Repeating this hypothesis test (but computing the p-value as was described before Eq. 54) for the remaining three modes, we find the two parity modes and the cross mode have p-values that are consistent with a non-detection; that is, p = 0.90 for ⟨ℳ_ap²ℳ_×⟩, p = 0.26 for ⟨ℳ_׳⟩, and p = 0.09 for ⟨ℳ_apℳ_ײ⟩. Besides the detection significance in the noise-only case, we can also compute the p-value and the S/N of ⟨ℳ_ap³⟩, for which we find p = 0.05 and S/N = 3.19. We list these values, as well as the values for all other setups described in this and the following subsection, in Table 1.