© The Authors 2025

1. Introduction

The Lambda cold dark matter (ΛCDM) model is the prevailing framework for understanding the evolution of the Universe, from the post-inflationary epoch to the development of the large-scale structure (LSS) observed today. It is characterised by six key parameters that must be empirically determined. Observational efforts to constrain these parameters fall into two broad categories: those that examine the early Universe through measurements of the cosmic microwave background (CMB), such as those provided by Planck (Planck Collaboration VI 2020), and those that probe the late-time Universe, notably via the LSS.

A longstanding area of interest in recent years has been the apparent tension between CMB-based and LSS-based determinations of the clustering amplitude parameter $S_8={\sigma }_8\sqrt{{\mathrm{\Omega }}_{\mathrm{m}}/0.3}$, where σ₈ quantifies the normalisation of the matter power spectrum and Ω_m denotes the total matter density. Several LSS surveys have reported lower values of S₈ compared to those inferred from the CMB, with tensions reaching 2 − 3σ significance (Di Valentino et al. 2021; Joudaki et al. 2020; Asgari et al. 2021). However, recent results from the KiDS-Legacy analysis (Wright et al. 2025) have shown that this tension can be resolved through a more robust statistical treatment of the data, suggesting that the discrepancy does not require new physics, but rather a refined methodological framework.

Nevertheless, the potential of weak gravitational lensing as a tool for probing cosmology remains undiminished. Through the measurement of background galaxy distortions induced by the intervening matter distribution, weak lensing offers a direct and unbiased probe of the total matter content, encompassing both baryonic and dark components, without assumptions about galaxy bias. This makes it uniquely suited for constraining key cosmological parameters such as Ω_m and S₈ (σ₈).

In the past decade, weak lensing analyses have primarily employed second-order statistics, such as the two-point shear correlation function, which have substantially improved parameter constraints (Hildebrandt et al. 2017; DES Collaboration 2022; Heymans et al. 2021; Dalal et al. 2023). However, second-order measures alone are limited in their ability to fully characterise the matter distribution, particularly at late times. The Gaussian approximation valid in the early Universe breaks down due to non-linear structure formation, which induces non-Gaussian features such as high-density peaks corresponding to galaxy clusters. These complexities are not captured by second-order statistics alone.

Moreover, degeneracies, especially between Ω_m and σ₈, remain a significant challenge when relying solely on two-point statistics (Takada & Jain 2004; Kilbinger & Schneider 2005; Euclid Collaboration: Ajani et al. 2023). To break these degeneracies and access the richer information encoded in the LSS, higher-order statistics and non-Gaussian probes have become an increasingly active area of development. These methods hold the potential not only to tighten parameter constraints, but also to offer a more complete statistical description of the cosmic matter distribution.

In the following, we focus on third-order statistics, which provide a measure of the skewness of the matter distribution. Since second- and third-order statistics have different dependences on Ω_m and S₈ (Takada & Jain 2004; Kilbinger & Schneider 2005; Kayo et al. 2013; Pyne & Joachimi 2021), combining them in a joint parameter analysis helps to break degeneracies, as demonstrated by Heydenreich et al. (2023). Various higher-order statistics have been introduced to address the limitations of second-order statistics in weak gravitational lensing. Examples include the shear three-point correlation function (Schneider & Lombardi 2003; Schneider et al. 2005), peak-count statistics (Martinet et al. 2018; Harnois-Déraps et al. 2021), persistent homology (Heydenreich et al. 2021), and density split statistics (Gruen et al. 2018; Burger et al. 2022). For our analysis, we adopted second- and third-order aperture mass statistics, ⟨ℳ_ap²⟩ and ⟨ℳ_ap³⟩, due to several advantages. Firstly, they are scalars, which effectively compress data and simplify handling (Heydenreich et al. 2023). Secondly, they are unaffected by the mass-sheet degeneracy, a transformation that adds a constant 1 − λ to the dimensionless surface mass density, while linearly scaling it by λ, leaving lensing observables unchanged (Falco et al. 1985; Schneider & Seitz 1995). Thirdly, they naturally decompose E- and B-modes, where B-modes should not be present in gravitational lensing to leading order, providing a valuable check for systematic errors. Lastly, they can be analytically modelled in a comparatively simple manner.

This paper is part of a series that presents a comprehensive approach to cosmological parameter estimation using third-order shear statistics. The series begins with an in-depth analysis of the analytical model for ⟨ℳ_ap³⟩ (Heydenreich et al. 2023). The second paper derives the analytical auto-covariance for ⟨ℳ_ap³⟩ (Linke et al. 2023, L+23 hereafter). The third paper Porth et al. (2024) introduces an efficient method for computing aperture mass statistics by measuring the shear three-point correlation function (3PCF). Finally, the last paper (Burger et al. 2024) presents the first numerical cosmological parameter analysis of the fourth data release of the Kilo-Degree Survey (KiDS-1000), incorporating second- and third-order shear statistics, intrinsic alignments, and the effects of baryonic feedback. The current paper establishes the foundation for an analytical joint analysis of ⟨ℳ_ap²⟩ and ⟨ℳ_ap³⟩ by deriving the final missing component: the analytical cross-covariance. Once derived, we validated the cross-covariance through a joint parameter analysis on simulated data. To achieve this, we used the analytical third-order auto-covariance from L+23 and the analytical second-order auto-covariance computed in Linke et al. (2024). By combining these covariances with the cross-covariance, we were able to construct posterior distributions using Markov chain Monte Carlo (MCMC) sampling.

Using an analytical covariance instead of a numerical one provides several advantages. A numerical covariance can be obtained either from simulations or directly from observational data. In the case of simulations, the number of realisations required must significantly exceed the number of components in the data vector, making the approach computationally expensive. Running these simulations and computing the covariance is particularly demanding, as data vectors are often large. In contrast, an analytical model significantly reduces computational cost while maintaining accuracy. Alternatively, the covariance matrix can be estimated using jackknife resampling or bootstrapping. These methods divide the data field into multiple patches that are small enough to ensure there are more patches than data vector components, and compute the sample covariance over these patches. However, this approach has limitations. It neglects correlations on scales larger than the patch size and assumes that the patches are mutually independent, which is not strictly true since adjacent patches share boundaries. As a result, this method introduces a bias in the covariance estimate, whereas an analytical covariance remains unaffected. Moreover, an analytical approach enables a straightforward evaluation of the covariance under different cosmological models. This is crucial for parameter inference and forecasting, as observational data by definition only samples a single cosmology, and generating large suites of simulations for multiple cosmologies would be prohibitively expensive.

The paper is organised as follows. Section 2 provides a concise review of aperture mass statistics. Section 3 introduces real-space estimators for aperture mass statistics and derives the analytical cross-covariance from these estimators. Additionally, we present a numerical method for validating the analytical results. Section 4 compares the analytical and numerical cross-covariances, where the analytical model is constructed by adapting the publicly available third-order modelling code threepoint¹. Section 5 outlines the joint cosmological parameter estimation and its validation against numerical contours. Finally, Sect. 6 discusses our findings and presents concluding remarks.

2. Theoretical background

This section briefly overviews the key quantities of weak gravitational lensing and aperture masses. A more detailed discussion can be found in Schneider (1996), Bartelmann & Schneider (2001), Hoekstra & Jain (2008) or Bartelmann (2010). This paper assumes a flat universe, such that the comoving angular-diameter distance equals the comoving radial distance f_K(χ) = χ. Furthermore, we use the flat-sky approximation. The matter density contrast is $\delta (\chi \mathit{\vartheta },\chi )=\rho (\chi \mathit{\vartheta },\chi )/\overline{\rho }(\chi )-1$, where $\rho (\chi \mathit{\vartheta },\chi )$ is the matter density at angular position ϑ and comoving distance χ on the backward light cone, and $\overline{\rho }(\chi )$ is the average density at χ. Defining the Fourier transform of variables with a tilde, the polyspectra of the density contrast 𝒫_n^(3d) are defined by

$$\begin{array}{c}\hfill \langle \stackrel{\sim }{\delta }({\ell }_1/\chi ,\chi )\cdots \stackrel{\sim }{\delta }({\ell }_n/\chi ,\chi ){\rangle }_{\mathrm{c}}={(2\pi )}^3{\delta }_{\mathrm{D}}({\ell }_1/\chi +\cdots +{\ell }_n/\chi ){\mathcal{P}}_n^{(3\mathrm{d})}({\ell }_1/\chi ,\cdots ,{\ell }_n/\chi ;\chi ),\end{array}$$(1)

where ⟨…⟩_c are cumulants or connected correlation functions, and δ_D is a Dirac delta function.

2.1. The convergence

The convergence, or the dimensionless surface mass density, can be calculated from δ by a weighted line-of-sight integration

$$\begin{array}{c}\hfill \kappa (\vartheta )=\frac{3H_0^2{\mathrm{\Omega }}_{\mathrm{m}}}{2c^2}{\displaystyle {\int }_0^{\infty }}\mathrm{d}\chi q(\chi )\chi \frac{\delta (\chi \vartheta ,\chi )}{a(\chi )},\text{with the weight function}q(\chi )={\displaystyle {\int }_{\chi }^{\infty }}\mathrm{d}\chi \prime p(\chi \prime )\frac{\chi \prime -\chi }{\chi \prime },\end{array}$$(2)

the Hubble constant H₀, the probability distribution of source galaxies in comoving distance p(χ) dχ, the speed of light c, and the scale factor a(χ), normalised to unity today. The nth-order polyspectra 𝒫_n of the convergence are defined by

$$\begin{array}{c}\hfill \langle \kappa ({\vartheta }_1)\cdots \kappa ({\vartheta }_n){\rangle }_{\mathrm{c}}={\displaystyle \int \frac{{\mathrm{d}}^2{\ell }_1}{{(2\pi )}^2}}\cdots \frac{{\mathrm{d}}^2{\ell }_n}{{(2\pi )}^2}{\mathcal{P}}_n({\ell }_1,\cdots ,{\ell }_n){(2\pi )}^2{\delta }_{\mathrm{D}}({\ell }_1+\cdots +{\ell }_n){\mathrm{e}}^{-\mathrm{i}({\ell }_1·{\vartheta }_1+\cdots +{\ell }_n·{\vartheta }_n)}.\end{array}$$(3)

Here and in the following, integrals without explicit boundaries are integrals over ℝ². Particularly, the power spectrum P, the bispectrum B, and the tetraspectrum P₅ are the convergence polyspectra important for our present discussion and are defined as

$$\begin{array}{c}\hfill P(\ell )={\mathcal{P}}_2(\ell ,-\ell ),B({\ell }_1,{\ell }_2)={\mathcal{P}}_3({\ell }_1,{\ell }_2,-{\ell }_1-{\ell }_2),\text{and}{P}_5({\ell }_1,{\ell }_2,{\ell }_3,{\ell }_4)={\mathcal{P}}_5({\ell }_1,{\ell }_2,{\ell }_3,{\ell }_4,-{\ell }_1-{\ell }_2-{\ell }_3-{\ell }_4),\end{array}$$(4)

where we made use of the homogeneity and isotropy of the Universe. Under the Limber approximation, they are related to the polyspectra of δ (𝒫_n^(3d)) by a weighted integration along the line of sight (Limber 1954; Kaiser & Jaffe 1997; Kayo et al. 2013):

$$\begin{array}{c}\hfill {\mathcal{P}}_n({\ell }_1,\cdots ,{\ell }_n)={\left(\frac{3H_0^2{\mathrm{\Omega }}_{\mathrm{m}}}{2c^2}\right)}^n{\displaystyle {\int }_0^{\infty }}\mathrm{d}\chi \frac{q^n(\chi )}{{\chi }^{n-2}{a}^n(\chi )}{\mathcal{P}}_n^{(3\mathrm{d})}({\ell }_1/\chi ,\cdots ,{\ell }_n/\chi ;\chi ).\end{array}$$(5)

The convergence defined in Eq. (2) and its nth-order connected correlation function (Eq. (3)) will be used to calculate the aperture mass in the following subsection, from which our statistics are derived. If shape noise from the intrinsic ellipticity of galaxies, ϵ_i, is present, it must be accounted for by adding σ_ϵi²/(2N) to the power spectrum P(ℓ), where σ_ϵi is the two-component intrinsic ellipticity dispersion and N is the galaxy number density. A proof of this correction can be found in Appendix B of L+23. Since shape noise is Gaussian, it does not affect the bispectrum or tetraspectrum.

2.2. Definition of the aperture mass

The aperture mass ℳ_ap can be calculated by either a weighted two-dimensional integration over the convergence or the tangential shear γ_t. It was first introduced in Schneider (1996), and can be defined by either

$$\begin{array}{c}\hfill {\mathcal{M}}_{\mathrm{ap}}(\varphi ,\theta ):={\displaystyle \int {\mathrm{d}}^2}\vartheta U_{\theta }(|\varphi -\vartheta |)\kappa (\vartheta )\text{or}{\mathcal{M}}_{\mathrm{ap}}(\varphi ,\theta ):={\displaystyle \int {\mathrm{d}}^2}\vartheta Q_{\theta }(|\varphi -\vartheta |){\gamma }_{\mathrm{t}}(\vartheta ),\end{array}$$(6)

where ϕ is the aperture centre coordinate, θ the filter radius, and U_θ and Q_θ are radially symmetric filter functions. U_θ is compensated, obeying the relation

$$\begin{array}{c}\hfill {\displaystyle {\int }_0^{\infty }}\mathrm{d}\vartheta \vartheta U_{\theta }(\vartheta )=0,\text{and is connected to}{Q}_{\theta }\text{through}{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}$$(7)

For our purpose, we use the exponential filter function proposed in Crittenden et al. (2002),

$$\begin{array}{c}\hfill U_{\theta }(\vartheta )=\frac{1}{2\pi {\theta }^2}(1-\frac{{\vartheta }^2}{2{\theta }^2}){\mathrm{e}}^{-{\vartheta }^2/(2{\theta }^2)},\text{such that}{Q}_{\theta }(\vartheta )=\frac{{\vartheta }^2}{4\pi {\theta }^4}{\mathrm{e}}^{-{\vartheta }^2/(2{\theta }^2)}.\end{array}$$(8)

The Fourier transform of U_θ is

$$\begin{array}{c}\hfill {\stackrel{\sim }{U}}_{\theta }(\ell )=\frac{{\ell }^2{\theta }^2}{2}{\mathrm{e}}^{-{\ell }^2{\theta }^2/2}=:\stackrel{\sim }{u}(\ell \theta ).\end{array}$$(9)

The second- and third-order aperture mass statistics are essential for our present discussion and are defined, respectively, by

$$\begin{array}{c}\hfill ⟨{\mathcal{M}}_{\mathrm{ap}}^2⟩({\theta }_1,{\theta }_2):=\langle {\mathcal{M}}_{\mathrm{ap}}(\varphi ,{\theta }_1){\mathcal{M}}_{\mathrm{ap}}(\varphi ,{\theta }_2)\rangle \text{and}⟨{\mathcal{M}}_{\mathrm{ap}}^3⟩({\theta }_1,{\theta }_2,{\theta }_3):=\langle {\mathcal{M}}_{\mathrm{ap}}(\varphi ,{\theta }_1){\mathcal{M}}_{\mathrm{ap}}(\varphi ,{\theta }_2){\mathcal{M}}_{\mathrm{ap}}(\varphi ,{\theta }_3)\rangle ,\end{array}$$(10)

where ⟨…⟩ denotes the spatial average over ϕ. We note that since ⟨ℳ_ap⟩ = 0, as the expectation value of the density contrast δ is by definition zero, the second- and third-order cumulants of the aperture mass are equal to the spatial averages: ⟨ℳ_ap³⟩_c = ⟨ℳ_ap³⟩ and ⟨ℳ_ap²⟩_c = ⟨ℳ_ap²⟩. Inserting the definition of the aperture mass calculated from the convergence (Eq. (6)) into the Eq. (10), using that the spatial averages are equal to the cumulants and substituting the cumulants of the convergence with the polyspectra as defined in Eq. (3), one finds

$$\begin{array}{cc}\hfill ⟨{\mathcal{M}}_{\mathrm{ap}}^2⟩({\theta }_1,{\theta }_2)& ={\displaystyle \int {\mathrm{d}}^2{\vartheta }_1{U}_{{\theta }_1}}\left(|\varphi -{\vartheta }_1|\right){\displaystyle \int {\mathrm{d}}^2}{\vartheta }_2{U}_{{\theta }_2}\left(|\varphi -{\vartheta }_2|\right){\displaystyle \int \frac{{\mathrm{d}}^2{\ell }_1}{{(2\pi )}^2}\int \frac{{\mathrm{d}}^2{\ell }_2}{{(2\pi )}^2}}{\mathcal{P}}_2({\ell }_1,{\ell }_2){(2\pi )}^2{\delta }_{\mathrm{D}}({\ell }_1+{\ell }_2){\mathrm{e}}^{-\mathrm{i}({\ell }_1·{\vartheta }_1+{\ell }_2·{\vartheta }_2)}\hfill \\ \hfill & ={\displaystyle \int \frac{{\mathrm{d}}^2{\ell }_1}{{(2\pi )}^2}\int \frac{{\mathrm{d}}^2{\ell }_2}{{(2\pi )}^2}}\stackrel{\sim }{u}({\ell }_1{\theta }_1)\stackrel{\sim }{u}({\ell }_2{\theta }_2){\mathcal{P}}_2({\ell }_1,{\ell }_2){(2\pi )}^2{\delta }_{\mathrm{D}}({\ell }_1+{\ell }_2){\mathrm{e}}^{-\mathrm{i}({\ell }_1+{\ell }_2)·\varphi }\hfill \\ \hfill & ={\displaystyle \int \frac{{\mathrm{d}}^2\ell }{{(2\pi )}^2}}\stackrel{\sim }{u}(\ell {\theta }_1)\stackrel{\sim }{u}(\ell {\theta }_2)P(\ell ),\hfill \end{array}$$(11)

for the second-order aperture mass statistics, where going from the first to the second line, we performed the Fourier transformation on the U_θ filters. In the last step, we integrated the Dirac delta function. It is important to note that taking equal aperture radii for the second-order aperture mass statistics (i.e., θ₁ = θ₂) captures all the information in this order. With the chosen $\stackrel{\sim }{u}$ filter (Eq. (9)), the dispersion at two distinct aperture radii can be expressed as a dispersion at an average aperture radius (Schneider et al. 2005)

$$\begin{array}{c}\hfill \langle {\mathcal{M}}_{\mathrm{ap}}({\theta }_1){\mathcal{M}}_{\mathrm{ap}}({\theta }_2)\rangle ={\displaystyle \int \frac{{\mathrm{d}}^2\ell }{{(2\pi )}^2}}P(\ell )\frac{l^4{\theta }_1^2{\theta }_2^2}{4}{\mathrm{e}}^{-l^2({\theta }_1^2+{\theta }_2^2)/2}=\frac{4{\theta }_1^2{\theta }_2^2}{{({\theta }_1^2+{\theta }_2^2)}^2}\langle {\mathcal{M}}_{\mathrm{ap}}^2\left(\sqrt{\frac{{\theta }_1^2+{\theta }_2^2}{2}}\right)\rangle .\end{array}$$(12)

Therefore, using equal aperture radii in the dispersion preserves all relevant information, and it is sufficient to evaluate ⟨ℳ_ap²⟩(θ):=⟨ℳ_ap²⟩(θ, θ) with a single aperture radius as input. Equivalently, for the third-order

$$\begin{array}{cc}\hfill ⟨{\mathcal{M}}_{\mathrm{ap}}^3⟩({\theta }_1,{\theta }_2,{\theta }_3)& ={\displaystyle \int \frac{{\mathrm{d}}^2{\ell }_1}{{(2\pi )}^2}}{\displaystyle \int \frac{{\mathrm{d}}^2{\ell }_2}{{(2\pi )}^2}}{\displaystyle \int \frac{{\mathrm{d}}^2{\ell }_3}{{(2\pi )}^2}}\stackrel{\sim }{u}({\ell }_1{\theta }_1)\stackrel{\sim }{u}({\ell }_2{\theta }_2)\stackrel{\sim }{u}({\ell }_3{\theta }_3){\mathcal{P}}_3({\ell }_1,{\ell }_2,{\ell }_3){(2\pi )}^2{\delta }_{\mathrm{D}}({\ell }_1+{\ell }_2+{\ell }_3){\mathrm{e}}^{-\mathrm{i}({\ell }_1+{\ell }_2+{\ell }_3)·\varphi }\hfill \\ \hfill & ={\displaystyle \int \frac{{\mathrm{d}}^2{\ell }_1}{{(2\pi )}^2}}{\displaystyle \int \frac{{\mathrm{d}}^2{\ell }_2}{{(2\pi )}^2}}\stackrel{\sim }{u}({\ell }_1{\theta }_1)\stackrel{\sim }{u}({\ell }_2{\theta }_2)\stackrel{\sim }{u}(\mid {\ell }_1+{\ell }_2\mid {\theta }_3)B({\ell }_1,{\ell }_2).\hfill \end{array}$$(13)

The simplification to equal aperture radii used above is not applicable to the third-order statistics. This can be observed from Eq. (9), which sharply peaks around $\ell \theta =\sqrt{2}$. For ⟨ℳ_ap³⟩, if one were to take θ₁ = θ₂ = θ₃ in Eq. (13), only equilateral configurations of the bispectrum with ${\ell }_1\sim {\ell }_2\sim \mid {\mathit{\ell }}_1+{\mathit{\ell }}_2\mid$ would be probed (Schneider et al. 2005). However, the non-equilateral aperture radii configurations depend on different bispectrum configurations and provide additional information in this order. Therefore, distinct aperture radii should be allowed for third-order statistics.

3. Derivation of the analytical cross-covariance

In a real survey or simulation, convolving with the weight function U_θ over an infinite area is not feasible. Therefore, following L+23, we define a real-space estimator of the aperture mass as

$$\begin{array}{c}\hfill {\widehat{\mathcal{M}}}_{\mathrm{ap}}(\varphi ,\theta ):={\int }_{A\prime }{\mathrm{d}}^2\vartheta U_{\theta }(|\varphi -\vartheta |)\kappa (\vartheta ),\end{array}$$(14)

where A′ represents the finite survey area. By spatially averaging products of this real-space estimator over the aperture centre ϕ, we estimate the second- and third-order aperture mass statistics. However, since we are working with a finite survey area, this spatial average is itself an approximation. Moreover, the averaging region must be smaller than the full survey area, as the aperture mass involves integrating over a region surrounding the aperture centre. Placing aperture centres too close to the survey boundary would introduce biases due to boundary effects. To mitigate these effects, we estimate the second- and third-order aperture mass statistics by averaging over an area A, which is entirely contained within the survey region A′ as illustrated in Fig. 1. The area A is chosen such that 99.9% of the effective support of Q_θ is within A′ when the aperture centre ϕ is placed on the boundary of A and the largest aperture radius θ_max is used. This condition is satisfied when the boundary of A is positioned 4 θ_max inside the boundary of A′ (Heydenreich et al. 2023). The field boundary of A is determined using Q_θ rather than U_θ because U_θ is a compensated filter function, meaning its support is not well-defined in the same way. Given this set-up, the integrals over the finite survey area A′ can be approximated as integrals over the entire real space ℝ², since the convolution with U_θ contributes negligibly outside of A′. The real-space statistics of the aperture mass are thus given by

$$\begin{array}{c}\hfill {\widehat{\mathcal{M}}}_{\mathrm{ap}}^2(\theta )=\frac{1}{A}{\int }_A{\mathrm{d}}^2\varphi \left[{\displaystyle \prod _{i=1}^2}{\displaystyle \int {\mathrm{d}}^2}{\vartheta }_i{U}_{\theta }\left(|\varphi -{\vartheta }_i|\right)\kappa ({\vartheta }_i)\right]\text{and}{\widehat{\mathcal{M}}}_{\mathrm{ap}}^3({\theta }_1,{\theta }_2,{\theta }_3)=\frac{1}{A}{\int }_A{\mathrm{d}}^2\varphi \left[{\displaystyle \prod _{i=1}^3}{\displaystyle \int {\mathrm{d}}^2}{\vartheta }_i{U}_{{\theta }_i}\left(|\varphi -{\vartheta }_i|\right)\kappa ({\vartheta }_i)\right].\end{array}$$(15)

Fig. 1. Aperture centre placement in a data field (from L+23). Given a data field of area A′ (not necessarily a square field), apertures are only placed within a smaller field A, which lies well within A′. The boundary of A is chosen such that 99.9% of the effective support of Qθ lies within A (which corresponds to ∼4 θ) if the aperture centre lies on this boundary. To illustrate this, two aperture centres are depicted: ϕ is within A and ϕ″ is outside of A. It can be seen that the circle around ϕ″, illustrating the effective support of Qθ, extends beyond A′, thus introducing bias into the calculated aperture mass at this point. To exclude this effect, we only average the aperture centre over the red region A.

The cross-covariance of the real-space estimators is then defined as

$$\begin{array}{c}\hfill C_{{\widehat{\mathcal{M}}}_{\mathrm{ap}}^{3,2}}({\theta }_1,{\theta }_2,{\theta }_3;{\theta }_4)=\langle {\widehat{\mathcal{M}}}_{\mathrm{ap}}^3{\widehat{\mathcal{M}}}_{\mathrm{ap}}^2\rangle ({\theta }_1,{\theta }_2,{\theta }_3;{\theta }_4)-\langle {\widehat{\mathcal{M}}}_{\mathrm{ap}}^3\rangle ({\theta }_1,{\theta }_2,{\theta }_3)\langle {\widehat{\mathcal{M}}}_{\mathrm{ap}}^2\rangle ({\theta }_4).\end{array}$$(16)

Here and in the following, when a comma ‘,’ is used, the function is invariant under permutations of its arguments, whereas when a semi-colon ‘;’ is used, they are not invariant. The first term in Eq. (16) is given as

$$\begin{array}{cc}\hfill \langle {\widehat{\mathcal{M}}}_{\mathrm{ap}}^3{\widehat{\mathcal{M}}}_{\mathrm{ap}}^2\rangle ({\theta }_1,{\theta }_2,{\theta }_3;{\theta }_4)=& \frac{1}{A^2}{\int }_A{\mathrm{d}}^2{\varphi }_1{\int }_A{\mathrm{d}}^2{\varphi }_2\left[{\displaystyle \prod _{i=1}^3}{\displaystyle \int {\mathrm{d}}^2}{\vartheta }_i{U}_{{\theta }_i}(|{\varphi }_1-{\vartheta }_i|)\right]\left[{\displaystyle \prod _{j=4}^5}{\displaystyle \int {\mathrm{d}}^2}{\vartheta }_j{U}_{{\theta }_4}(|{\varphi }_2-{\vartheta }_j|)\right]\hfill \\ \hfill & \times \langle \kappa ({\vartheta }_1)\kappa ({\vartheta }_2)\kappa ({\vartheta }_3)\kappa ({\vartheta }_4)\kappa ({\vartheta }_5)\rangle .\hfill \end{array}$$(17)

Here, the fifth-order correlation function can be separated into cumulants,

$$\begin{array}{cc}\hfill \langle {\widehat{\mathcal{M}}}_{\mathrm{ap}}^3{\widehat{\mathcal{M}}}_{\mathrm{ap}}^2\rangle ({\theta }_1,{\theta }_2,{\theta }_3;{\theta }_4)=& \frac{1}{A^2}{\int }_A{\mathrm{d}}^2{\varphi }_1{\int }_A{\mathrm{d}}^2{\varphi }_2\left[{\displaystyle \prod _{i=1}^3}{\displaystyle \int {\mathrm{d}}^2}{\vartheta }_i{U}_{{\theta }_i}(|{\varphi }_1-{\vartheta }_i|)\right]\left[{\displaystyle \prod _{j=4}^5}{\displaystyle \int {\mathrm{d}}^2}{\vartheta }_j{U}_{{\theta }_4}(|{\varphi }_2-{\vartheta }_j|)\right]\hfill \\ \hfill & \times \{\langle \kappa ({\vartheta }_1)\kappa ({\vartheta }_2)\kappa ({\vartheta }_3){\rangle }_{\mathrm{c}}\langle \kappa ({\vartheta }_4)\kappa ({\vartheta }_5){\rangle }_{\mathrm{c}}\hfill \\ \hfill & +\langle \kappa ({\vartheta }_1)\kappa ({\vartheta }_4)\kappa ({\vartheta }_5){\rangle }_{\mathrm{c}}\langle \kappa ({\vartheta }_2)\kappa ({\vartheta }_3){\rangle }_{\mathrm{c}}+2\mathrm{Perm}.\mathrm{of}({\theta }_1,{\theta }_2,{\theta }_3)\hfill \\ \hfill & +2[\langle \kappa ({\vartheta }_1)\kappa ({\vartheta }_2)\kappa ({\vartheta }_4){\rangle }_{\mathrm{c}}\langle \kappa ({\vartheta }_3)\kappa ({\vartheta }_5){\rangle }_{\mathrm{c}}+2\mathrm{Perm}.\mathrm{of}({\theta }_1,{\theta }_2,{\theta }_3)]\hfill \\ \hfill & +\langle \kappa ({\vartheta }_1)\kappa ({\vartheta }_2)\kappa ({\vartheta }_3)\kappa ({\vartheta }_4)\kappa ({\vartheta }_5){\rangle }_{\mathrm{c}}\}\hfill \\ \hfill =:& T_{\mathrm{PB},1}({\theta }_1,{\theta }_2,{\theta }_3;{\theta }_4)+T_{\mathrm{PB},2}({\theta }_1,{\theta }_2,{\theta }_3;{\theta }_4)+T_{\mathrm{PB},3}({\theta }_1,{\theta }_2,{\theta }_3;{\theta }_4)+T_{\mathrm{P}5}({\theta }_1,{\theta }_2,{\theta }_3;{\theta }_4),\hfill \end{array}$$(18)

where we defined the T-expressions such that T_PB, 1 corresponds to the term in the second line, T_PB, 2 to the third, T_PB, 3 to the fourth, and T_P5 to the fifth. The permutations of T_PB, 2 and T_PB, 3 can be found in Table 1. The third term (T_PB, 3) includes a factor of two because swapping $\kappa ({\mathit{\vartheta }}_4)$ with $\kappa ({\mathit{\vartheta }}_5)$ results in equivalent terms, as both are connected to θ₄. Remember ⟨κ⟩ = 0 because the expectation value of the density contrast δ is zero; thus, the second- and third-order cumulants of the convergence are equal to the respective moments: ⟨κ³⟩_c = ⟨κ³⟩ and ⟨κ²⟩_c = ⟨κ²⟩ (but ⟨κ⁵⟩_c ≠ ⟨κ⁵⟩). As a result, T_PB, 1 cancels the second term on the r.h.s. of Eq. (16), observing it can be written as the expectation value of ${\widehat{\mathcal{M}}}_{\mathrm{ap}}^2$ times the expectation value of ${\widehat{\mathcal{M}}}_{\mathrm{ap}}^3$, such that

$$\begin{array}{c}\hfill C_{{\widehat{\mathcal{M}}}_{\mathrm{ap}}^{3,2}}({\theta }_1,{\theta }_2,{\theta }_3;{\theta }_4)=T_{\mathrm{PB},2}({\theta }_1,{\theta }_2,{\theta }_3;{\theta }_4)+T_{\mathrm{PB},3}({\theta }_1,{\theta }_2,{\theta }_3;{\theta }_4)+T_{\mathrm{P}5}({\theta }_1,{\theta }_2,{\theta }_3;{\theta }_4).\end{array}$$(19)

The second term T_PB, 2 can be rewritten by substituting the definition for the convergence polyspectra (Eq. (3)) into Eq. (18), cancelling ℓ₃ and ℓ₅ by integrating over the Dirac delta functions and performing the Fourier transformations of the U_θ functions

$$\begin{array}{cc}\hfill T_{\mathrm{PB},2}({\theta }_1,{\theta }_2,{\theta }_3;{\theta }_4)=& \frac{1}{A^2}{\int }_A{\mathrm{d}}^2{\varphi }_1{\int }_A{\mathrm{d}}^2{\varphi }_2\left[{\displaystyle \prod _{i=1}^3}{\displaystyle \int {\mathrm{d}}^2}{\vartheta }_i{U}_{{\theta }_i}(|{\varphi }_1-{\vartheta }_i|)\right]\left[{\displaystyle \prod _{j=4}^5}{\displaystyle \int {\mathrm{d}}^2}{\vartheta }_j{U}_{{\theta }_4}(|{\varphi }_2-{\vartheta }_j|)\right]{\displaystyle \prod _{k=1}^5}{\displaystyle \int \frac{{\mathrm{d}}^2{\ell }_k}{{(2\pi )}^6}}\hfill \\ \hfill & \times {\delta }_{\mathrm{D}}({\ell }_2+{\ell }_3){\delta }_{\mathrm{D}}({\ell }_1+{\ell }_4+{\ell }_5){\mathcal{P}}_3({\ell }_1,{\ell }_4,{\ell }_5){\mathcal{P}}_2({\ell }_2,{\ell }_3){\mathrm{e}}^{-\mathrm{i}({\ell }_1·{\vartheta }_1+{\ell }_2·{\vartheta }_2+{\ell }_3·{\vartheta }_3+{\ell }_4·{\vartheta }_4+{\ell }_5·{\vartheta }_5)}\hfill \\ \hfill & +2\mathrm{Perm}.\mathrm{of}({\theta }_1,{\theta }_2,{\theta }_3)\hfill \\ \hfill =& \frac{1}{A^2}{\int }_A{\mathrm{d}}^2{\varphi }_1{\int }_A{\mathrm{d}}^2{\varphi }_2\left[{\displaystyle \prod _{i=1}^3}{\displaystyle \int {\mathrm{d}}^2}{\vartheta }_i{U}_{{\theta }_i}(|{\varphi }_1-{\vartheta }_i|)\right]\left[{\displaystyle \prod _{j=4}^5}{\displaystyle \int {\mathrm{d}}^2}{\vartheta }_j{U}_{{\theta }_4}(|{\varphi }_2-{\vartheta }_j|)\right]{\displaystyle \int \frac{{\mathrm{d}}^2{\ell }_1}{{(2\pi )}^2}}{\displaystyle \int \frac{{\mathrm{d}}^2{\ell }_2}{{(2\pi )}^2}}\hfill \\ \hfill & \times {\displaystyle \int \frac{{\mathrm{d}}^2{\ell }_4}{{(2\pi )}^2}}B({\ell }_1,{\ell }_4)P({\ell }_2){\mathrm{e}}^{-\mathrm{i}({\ell }_1·({\vartheta }_1-{\vartheta }_5)+{\ell }_2·({\vartheta }_2-{\vartheta }_3)+{\ell }_4·({\vartheta }_4-{\vartheta }_5))}+2\mathrm{Perm}.\mathrm{of}({\theta }_1,{\theta }_2,{\theta }_3)\hfill \\ \hfill =& \frac{1}{A^2}{\int }_A{\mathrm{d}}^2{\varphi }_1{\int }_A{\mathrm{d}}^2{\varphi }_2{\displaystyle \int \frac{{\mathrm{d}}^2{\ell }_1}{{(2\pi )}^2}}{\displaystyle \int \frac{{\mathrm{d}}^2{\ell }_2}{{(2\pi )}^2}}{\displaystyle \int \frac{{\mathrm{d}}^2{\ell }_4}{{(2\pi )}^2}}B({\ell }_1,{\ell }_4)P({\ell }_2)\hfill \\ \hfill & \times \stackrel{\sim }{u}({\ell }_1{\theta }_1)\stackrel{\sim }{u}({\ell }_2{\theta }_2)\stackrel{\sim }{u}({\ell }_2{\theta }_3)\stackrel{\sim }{u}({\ell }_4{\theta }_4)\stackrel{\sim }{u}(|{\ell }_1+{\ell }_4|{\theta }_4){\mathrm{e}}^{-\mathrm{i}{\ell }_1·({\varphi }_1-{\varphi }_2)}+2\mathrm{Perm}.\mathrm{of}({\theta }_1,{\theta }_2,{\theta }_3).\hfill \end{array}$$(20)

Collecting all terms containing information on the survey area and geometry in the geometry factor $G_A(\mathit{\ell })$ defined as (L+23)

$$\begin{array}{c}\hfill G_A(\ell )=\frac{1}{A^2}{\int }_A{\mathrm{d}}^2{\varphi }_1{\int }_A{\mathrm{d}}^2{\varphi }_2{\mathrm{e}}^{-\mathrm{i}\ell ·({\varphi }_1-{\varphi }_2)},\end{array}$$(21)

we find the final result

$$\begin{array}{cc}\hfill T_{\mathrm{PB},2}({\theta }_1,{\theta }_2,{\theta }_3;{\theta }_4)=& {\displaystyle \int \frac{{\mathrm{d}}^2{\ell }_1}{{(2\pi )}^2}}{\displaystyle \int \frac{{\mathrm{d}}^2{\ell }_2}{{(2\pi )}^2}}{\displaystyle \int \frac{{\mathrm{d}}^2{\ell }_4}{{(2\pi )}^2}}B({\ell }_1,{\ell }_4)P({\ell }_2)\stackrel{\sim }{u}({\ell }_1{\theta }_1)\stackrel{\sim }{u}({\ell }_2{\theta }_2)\stackrel{\sim }{u}({\ell }_2{\theta }_3)\stackrel{\sim }{u}({\ell }_4{\theta }_4)\stackrel{\sim }{u}(|{\ell }_1+{\ell }_4|{\theta }_4)\hfill \\ \hfill & \times G_A({\ell }_1)+2\mathrm{Perm}.\mathrm{of}({\theta }_1,{\theta }_2,{\theta }_3).\hfill \end{array}$$(22)

By equivalent calculations, one finds

$$\begin{array}{cc}\hfill T_{\mathrm{PB},3}({\theta }_1,{\theta }_2,{\theta }_3;{\theta }_4)=& 2{\displaystyle \int \frac{{\mathrm{d}}^2{\ell }_1}{{(2\pi )}^2}}{\displaystyle \int \frac{{\mathrm{d}}^2{\ell }_2}{{(2\pi )}^2}}{\displaystyle \int \frac{{\mathrm{d}}^2{\ell }_3}{{(2\pi )}^2}}B({\ell }_2,{\ell }_3)P({\ell }_1)\stackrel{\sim }{u}({\ell }_1{\theta }_1)\stackrel{\sim }{u}({\ell }_2{\theta }_2)\stackrel{\sim }{u}({\ell }_3{\theta }_3)\stackrel{\sim }{u}(|{\ell }_2+{\ell }_3|{\theta }_4)\stackrel{\sim }{u}({\ell }_1{\theta }_4)\hfill \\ \hfill & \times G_A({\ell }_1+{\ell }_2+{\ell }_3)+2\mathrm{Perm}.\mathrm{of}({\theta }_1,{\theta }_2,{\theta }_3),\hfill \end{array}$$(23)

and

$$\begin{array}{cc}\hfill T_{\mathrm{P}5}({\theta }_1,{\theta }_2,{\theta }_3;{\theta }_4)=& {\displaystyle \int \frac{{\mathrm{d}}^2{\ell }_1}{{(2\pi )}^2}}{\displaystyle \int \frac{{\mathrm{d}}^2{\ell }_2}{{(2\pi )}^2}}{\displaystyle \int \frac{{\mathrm{d}}^2{\ell }_3}{{(2\pi )}^2}}{\displaystyle \int \frac{{\mathrm{d}}^2{\ell }_4}{{(2\pi )}^2}}{P}_5({\ell }_1,{\ell }_2,{\ell }_3,{\ell }_4)\stackrel{\sim }{u}({\ell }_1{\theta }_1)\stackrel{\sim }{u}({\ell }_2{\theta }_2)\stackrel{\sim }{u}({\ell }_3{\theta }_3)\stackrel{\sim }{u}({\ell }_4{\theta }_4)\hfill \\ \hfill & \times \stackrel{\sim }{u}(|{\ell }_1+{\ell }_2+{\ell }_3+{\ell }_4|{\theta }_4)G_A({\ell }_1+{\ell }_2+{\ell }_3).\hfill \end{array}$$(24)

Equations ((22)–(24)) can subsequently be substituted into Eq. (19) to obtain the analytical covariance of the real-space estimator. Next, we investigate how these terms scale with the survey field size.

3.1. Large survey area approximation

This subsection illustrates how the cross-covariance terms behave as the field size increases, by taking the limit A → ∞ in the geometry factor (Eq. (21)). In this limit, one finds (L+23)

$$\begin{array}{c}\hfill {\displaystyle \underset{A\to \infty }{lim}}{G}_A(\ell )=\frac{{(2\pi )}^2}{A}{\delta }_{\mathrm{D}}(\ell ).\end{array}$$(25)

Substituting and eliminating ℓ₁ by integrating over the Dirac delta function, the large-field approximation of T_PB, 2 becomes

$$\begin{array}{cc}\hfill T_{\mathrm{PB},2}^{\infty }({\theta }_1,{\theta }_2,{\theta }_3;{\theta }_4)=& \frac{1}{A}{\displaystyle \int \frac{{\mathrm{d}}^2{\ell }_2}{{(2\pi )}^2}}{\displaystyle \int \frac{{\mathrm{d}}^2{\ell }_4}{{(2\pi )}^2}}B(0,{\ell }_4)P({\ell }_2)\stackrel{\sim }{u}(0)\stackrel{\sim }{u}({\ell }_2{\theta }_2)\stackrel{\sim }{u}({\ell }_2{\theta }_3)\stackrel{\sim }{u}({\ell }_4{\theta }_4)\stackrel{\sim }{u}({\ell }_4{\theta }_4)\hfill \\ \hfill & +2\mathrm{Permutations}\mathrm{of}({\theta }_1,{\theta }_2,{\theta }_3)=0,\hfill \end{array}$$(26)

because $\stackrel{\sim }{u}(0)=0$, as can be seen from Eq. (9). Thus, T_PB, 2 arises solely from the limited size of the survey field and will therefore be referred to as a finite-field term.

For T_PB, 3 and T_P5, as given in Eqs. (23) and (24) respectively, the foregoing limiting procedure leads to

$$\begin{array}{cc}\hfill T_{\mathrm{PB},3}^{\infty }({\theta }_1,{\theta }_2,{\theta }_3;{\theta }_4)=& \frac{2}{A}{\displaystyle \int \frac{{\mathrm{d}}^2{\ell }_1}{{(2\pi )}^2}}{\displaystyle \int \frac{{\mathrm{d}}^2{\ell }_2}{{(2\pi )}^2}}B({\ell }_2,-{\ell }_1-{\ell }_2)P({\ell }_1)\stackrel{\sim }{u}({\ell }_1{\theta }_1)\stackrel{\sim }{u}({\ell }_2{\theta }_2)\stackrel{\sim }{u}(|{\ell }_1+{\ell }_2|{\theta }_3)\stackrel{\sim }{u}({\ell }_1{\theta }_4)\stackrel{\sim }{u}({\ell }_1{\theta }_4)\hfill \\ \hfill & +2\mathrm{Permutations}\mathrm{of}({\theta }_1,{\theta }_2,{\theta }_3),\hfill \end{array}$$(27)

and

$$\begin{array}{c}\hfill T_{\mathrm{P}5}^{\infty }({\theta }_1,{\theta }_2,{\theta }_3;{\theta }_4)=\frac{1}{A}{\displaystyle \int \frac{{\mathrm{d}}^2{\ell }_1}{{(2\pi )}^2}}{\displaystyle \int \frac{{\mathrm{d}}^2{\ell }_2}{{(2\pi )}^2}}{\displaystyle \int \frac{{\mathrm{d}}^2{\ell }_4}{{(2\pi )}^2}}{P}_5({\ell }_1,{\ell }_2,-{\ell }_1-{\ell }_2,{\ell }_4)\stackrel{\sim }{u}({\ell }_1{\theta }_1)\stackrel{\sim }{u}({\ell }_2{\theta }_2)\stackrel{\sim }{u}(|{\ell }_1+{\ell }_2|{\theta }_3)\stackrel{\sim }{u}({\ell }_4{\theta }_4)\stackrel{\sim }{u}({\ell }_4{\theta }_4),\end{array}$$(28)

such that the cross-covariance in the large-field approximation simplifies to

$$\begin{array}{c}\hfill C_{{\widehat{\mathcal{M}}}_{\mathrm{ap}}^{3,2}}^{\infty }({\theta }_1,{\theta }_2,{\theta }_3;{\theta }_4)=T_{\mathrm{PB},3}^{\infty }({\theta }_1,{\theta }_2,{\theta }_3;{\theta }_4)+T_{\mathrm{P}5}^{\infty }({\theta }_1,{\theta }_2,{\theta }_3;{\theta }_4).\end{array}$$(29)

While T_PB, 3 and T_P5 scale approximately as 1/A, T_PB, 2 exhibits a steeper decline with increasing field size.

Since the geometry term (Eq. (21)) is a rapidly oscillating function, its integration presents numerical challenges. Because our computational resources are largely sufficient for this task (albeit with some limitations; see Sect. 6.1), we retain the full expressions. However, in scenarios with limited computing power, or when computational efficiency is prioritised, large-field approximations may serve as useful alternatives.

A concluding remark on the derived cross-covariance and its large-field approximation is that to compute their values, we estimate P₅ in T_P5 using a halo decomposition (Appendix A). T_P5 is then approximated by keeping the two dominant terms, the one-halo term and one of the two-halo terms:

$$\begin{array}{c}\hfill T_{\mathrm{P}5}({\theta }_1,{\theta }_2,{\theta }_3;{\theta }_4)\approx T_{\mathrm{P}5,1\mathrm{h}}({\theta }_1,{\theta }_2,{\theta }_3;{\theta }_4)+T_{\mathrm{P}5,2\mathrm{h}}({\theta }_1,{\theta }_2,{\theta }_3;{\theta }_4).\end{array}$$(30)

In the large-field approximation T_P5, 2h^∞ = 0, which means this is a finite field term, and thus we find T_P5^∞ ≈ T_P5, 1h^∞.

3.2. Aperture mass correlation functions

This section demonstrates that the covariance terms (Eqs. (22)–(24)) can alternatively be constructed with aperture mass correlation functions. These correlation functions can be numerically estimated from N-body simulations and therefore enable us to verify the individual components of our analytical cross-covariance.

The aperture mass correlation functions are defined as

$$\begin{array}{cc}& {\mathcal{M}}_{\mathrm{ap}}^{n,m}({\theta }_1,\cdots ,{\theta }_m;{\theta }_{m+1},\cdots ,{\theta }_n;\eta ):=\langle {\mathcal{M}}_{\mathrm{ap}}(\varphi ,{\theta }_1)\cdots {\mathcal{M}}_{\mathrm{ap}}(\varphi ,{\theta }_m){\mathcal{M}}_{\mathrm{ap}}(\varphi +\eta ,{\theta }_{m+1})\cdots {\mathcal{M}}_{\mathrm{ap}}(\varphi +\eta ,{\theta }_n)\rangle \hfill \\ \hfill & =[{\displaystyle \prod _{i=1}^m}{\displaystyle \int {\mathrm{d}}^2}{\vartheta }_i{U}_{{\theta }_1}({\vartheta }_1)\cdots U_{{\theta }_m}({\vartheta }_m)][{\displaystyle \prod _{j=m+1}^n}{\displaystyle \int {\mathrm{d}}^2}{\vartheta }_j{U}_{{\theta }_{m+1}}(\mid {\vartheta }_{m+1}+\eta \mid )\cdots U_{{\theta }_n}(\mid {\vartheta }_n+\eta \mid )]\langle \kappa ({\vartheta }_1)\cdots \kappa ({\vartheta }_n)\rangle ,\hfill \end{array}$$(31)

where, without loss of generality, we set ϕ = 0 after the equality sign by shifting the origin. Substituting the correlation functions into Eq. (18), T_PB, 2 is rewritten as

$$\begin{array}{c}\hfill T_{\mathrm{PB},2}({\theta }_1,{\theta }_2,{\theta }_3;{\theta }_4)=\frac{1}{A^2}{\mathcal{M}}_{\mathrm{ap}}^{2,0}({\theta }_2,{\theta }_3;0){\int }_A{\mathrm{d}}^2{\varphi }_1{\int }_A{\mathrm{d}}^2{\varphi }_2{\mathcal{M}}_{\mathrm{ap}}^{3,1}({\theta }_1;{\theta }_4,{\theta }_4;{\varphi }_1-{\varphi }_2)+2\mathrm{Permutations}\mathrm{of}({\theta }_1,{\theta }_2,{\theta }_3).\end{array}$$(32)

By replacing the ϕ₂ integral by an η = ϕ₁−ϕ₂ integral and integrating over d[2]ϕ₁, we obtain $A{E}_A(\mathit{\eta })$, where $E_A(\mathit{\eta })$ represents the probability that for a point ϕ₁ within A, the point ϕ₁ + η also lies within A. For a square survey area, $E_A(\mathit{\eta })$ is given in Appendix A.2 of Heydenreich et al. (2020). By expanding these algebraic steps, we find

$$\begin{array}{c}\hfill T_{\mathrm{PB},2}({\theta }_1,{\theta }_2,{\theta }_3;{\theta }_4)=\frac{1}{A}{\mathcal{M}}_{\mathrm{ap}}^{2,0}({\theta }_2,{\theta }_3;0){\displaystyle \int {\mathrm{d}}^2}\eta E_A(\eta ){\mathcal{M}}_{\mathrm{ap}}^{3,1}({\theta }_1;{\theta }_4,{\theta }_4;\eta )+2\mathrm{Permutations}\mathrm{of}({\theta }_1,{\theta }_2,{\theta }_3).\end{array}$$(33)

After generating the aperture mass correlation function maps from simulations, we azimuthally average them over the polar angle of η due to the statistical isotropy of the Universe. Defining the aperture mass correlations as a function of distance η by ${\overline{\mathcal{M}}}_{\mathrm{ap}}^{n,m}$ (Eq. (31) cannot be used in this context, as it requires a vector-valued input η) and integrating $E_A(\mathit{\eta })$ over the polar angle of η, one finds

$$\begin{array}{c}\hfill T_{\mathrm{PB},2}({\theta }_1,{\theta }_2,{\theta }_3;{\theta }_4)=\frac{2\pi }{A}{\overline{\mathcal{M}}}_{\mathrm{ap}}^{2,0}({\theta }_2,{\theta }_3;0){\displaystyle {\int }_0^{\infty }}\mathrm{d}\eta E_A(\eta ){\overline{\mathcal{M}}}_{\mathrm{ap}}^{3,1}({\theta }_1;{\theta }_4,{\theta }_4;\eta )+2\mathrm{Permutations}\mathrm{of}({\theta }_1,{\theta }_2,{\theta }_3),\end{array}$$(34)

where for a survey of square area, E_A(η) is given in Appendix A.2 of Heydenreich et al. (2020)² as

$$\begin{array}{c}\hfill E_A(\eta )=\{\begin{array}{cc}\frac{1}{A\pi }[A\pi -(4\sqrt{A}-\eta )\eta ]\hfill & \text{for}\eta \le \sqrt{A},\hfill \\ \frac{2}{\pi }[2\sqrt{\frac{{\eta }^2}{A}-1}-1-\frac{{\eta }^2}{2A}-arccos\left(\frac{\sqrt{A}}{\eta }\right)+arcsin\left(\frac{\sqrt{A}}{\eta }\right)]\hfill & \text{for}\sqrt{A}\le \eta \le \sqrt{2A},\hfill \\ 0\hfill & \text{for}\sqrt{2A}\le \eta .\hfill \end{array}\end{array}$$(35)

Following the same procedure for T_PB, 1 and T_PB, 3 leads to

$$\begin{array}{cc}\hfill T_{\mathrm{PB},1}({\theta }_1,{\theta }_2,{\theta }_3;{\theta }_4)=& {\mathcal{M}}_{\mathrm{ap}}^{3,0}({\theta }_1,{\theta }_2,{\theta }_3;0){\mathcal{M}}_{\mathrm{ap}}^{2,0}({\theta }_4,{\theta }_4;0)\hfill \\ \hfill =& {\overline{\mathcal{M}}}_{\mathrm{ap}}^{3,0}({\theta }_1,{\theta }_2,{\theta }_3;0){\overline{\mathcal{M}}}_{\mathrm{ap}}^{2,0}({\theta }_4,{\theta }_4;0),\hfill \end{array}$$(36)

$$\begin{array}{cc}\hfill T_{\mathrm{PB},3}({\theta }_1,{\theta }_2,{\theta }_3;{\theta }_4)=& \frac{2}{A}\int {\mathrm{d}}^2\eta E_A(\eta ){\mathcal{M}}_{\mathrm{ap}}^{3,2}({\theta }_2,{\theta }_3;{\theta }_4;\eta ){\mathcal{M}}_{\mathrm{ap}}^{2,1}({\theta }_1;{\theta }_4;\eta )+2\mathrm{Permutations}\mathrm{of}({\theta }_1,{\theta }_2,{\theta }_3)\hfill \\ \hfill =& \frac{4\pi }{A}{\int }_0^{\infty }\mathrm{d}\eta E_A(\eta ){\overline{\mathcal{M}}}_{\mathrm{ap}}^{3,2}({\theta }_2,{\theta }_3;{\theta }_4;\eta ){\overline{\mathcal{M}}}_{\mathrm{ap}}^{2,1}({\theta }_1;{\theta }_4;\eta )+2\mathrm{Permutations}\mathrm{of}({\theta }_1,{\theta }_2,{\theta }_3).\hfill \end{array}$$(37)

We note that because Eq. (18) is written in cumulants, T_P5 cannot be simply connected to the fifth-order correlation function. Therefore, to calculate T_P5 in terms of correlation functions, we subtract T_PB, 1, T_PB, 2, and T_PB, 3 as defined above, such that

$$\begin{array}{cc}\hfill T_{\mathrm{P}5}({\theta }_1,{\theta }_2,{\theta }_3;{\theta }_4)=& \frac{1}{A}{\displaystyle \int {\mathrm{d}}^2}\eta E_A(\eta ){\mathcal{M}}_{\mathrm{ap}}^{5,3}({\theta }_1,{\theta }_2,{\theta }_3;{\theta }_4,{\theta }_4;\eta )\hfill \\ \hfill & -T_{\mathrm{PB},1}({\theta }_1,{\theta }_2,{\theta }_3;{\theta }_4)-T_{\mathrm{PB},2}({\theta }_1,{\theta }_2,{\theta }_3;{\theta }_4)-T_{\mathrm{PB},3}({\theta }_1,{\theta }_2,{\theta }_3;{\theta }_4)\hfill \\ \hfill =& \frac{2\pi }{A}{\displaystyle {\int }_0^{\infty }}\mathrm{d}\eta E_A(\eta ){\overline{\mathcal{M}}}_{\mathrm{ap}}^{5,3}({\theta }_1,{\theta }_2,{\theta }_3;{\theta }_4,{\theta }_4;\eta )\hfill \\ \hfill & -T_{\mathrm{PB},1}({\theta }_1,{\theta }_2,{\theta }_3;{\theta }_4)-T_{\mathrm{PB},2}({\theta }_1,{\theta }_2,{\theta }_3;{\theta }_4)-T_{\mathrm{PB},3}({\theta }_1,{\theta }_2,{\theta }_3;{\theta }_4).\hfill \end{array}$$(38)