© The Authors 2023

1. Introduction

Using weak gravitational lensing of the cosmic large-scale structure (LSS) has become a precise and accurate method to infer cosmological parameters (Heymans et al. 2021; Hikage et al. 2019; Abbott et al. 2022). Weak lensing analyses have mainly concentrated on second-order statistics, such as shear two-point correlation functions or galaxy-galaxy-lensing. However, since the cosmic LSS is non-Gaussian, second-order statistics do not contain its total information content. To access the additional information, higher-order statistics (HOS) are needed. These HOS depend on cosmological parameters differently than second-order statistics, so they can be used to break degeneracies (Takada & Jain 2004; Kilbinger & Schneider 2005; Kayo et al. 2013) or constrain nuisance parameters such as galaxy bias or intrinsic alignment (Pyne & Joachimi 2021; Troxel & Ishak 2012).

Many HOS have been suggested to complement second-order analyses, for example, peak statistics (Kacprzak et al. 2016; Zürcher et al. 2021; Martinet et al. 2021), density split statistics (Friedrich et al. 2018; Gruen et al. 2018; Burger et al. 2023), or persistent homology (Heydenreich et al. 2021, 2022). Most of these statistics, though, cannot be modelled from analytical theories and instead require time-consuming realistic N-body simulations. Additionally, for many HOS, the measurement requires converting weak lensing shear estimates to mass convergence maps. This conversion becomes complicated in the presence of masks and for finite survey areas, which can lead to biased estimators (Seitz & Schneider 1996).

In this series of papers, we are considering the third-order aperture mass $\langle M_{\text{ap}}^3\rangle$, which is not affected by these problems. It is linearly related to the third-order shear correlation functions Γ_i (Jarvis et al. 2004; Schneider et al. 2005), which can be directly measured in shear catalogues and converted into $\langle M_{\text{ap}}^3\rangle$. This was demonstrated recently by Secco et al. (2022), who measured $\langle M_{\text{ap}}^3\rangle$ with high significance from the Dark Energy Survey Year 3 shear catalogues. Moreover, in contrast to many other HOS, $\langle M_{\text{ap}}^3\rangle$ and Γ_i can be modelled analytically without directly relying on simulations.

There are several advantages to using $\langle M_{\text{ap}}^3\rangle$ instead of Γ_i, among them that $\langle M_{\text{ap}}^3\rangle$ compress the information in the Γ_i from a data vector with thousands of entries to one containing only tens of entries in a non-tomographic setup, and that the $\langle M_{\text{ap}}^3\rangle$ have a clear separation into E- and B-modes, allowing for the detection of untreated systematic effects. Our goal is to lay the groundwork for a cosmological analysis with $\langle M_{\text{ap}}^3\rangle$. Such an analysis requires two ingredients: a model of $\langle M_{\text{ap}}^3\rangle$ and an estimate of its covariance matrix. The first ingredient is the subject of Heydenreich et al. (2023, Paper I hereafter), where we present and test an analytic model of $\langle M_{\text{ap}}^3\rangle$ based on the BiHalofit bispectrum model (Takahashi et al. 2020). Here, we are concerned with the second ingredient, which is the covariance.

The covariance of the estimator of a statistic can be obtained in three ways. The first possibility is applying the estimator to a large set of (quasi-)independent cosmological N-body simulations and using the sample covariance of the estimates. However, an unbiased and precise estimate requires many more realisations than entries in the data vector. For example, for a tomographic analysis of $\langle M_{\text{ap}}^3\rangle$ with four redshift bins and four different scale radii, the data vector contains 400 entries, so hundreds of simulations are needed for a simulated covariance estimate. Such an estimate would be very time-consuming, as both the creation of this many independent simulations and the measurement of $\langle M_{\text{ap}}^3\rangle$ in them is computationally demanding. The second possibility is to estimate the covariance directly from the data using jackknife resampling or bootstrap methods. For this, the survey area is divided into smaller patches, the statistics are estimated separately on each of these patches, and the sample covariance of these individual estimates is taken as the covariance estimate. However, this method has two main disadvantages. First, information on correlations larger than the patch size is lost. Since the number of patches must be at least as large as the length of the data vector, the required number of patches becomes large and their size small. Second, the method implicitly assumes that the individual patches are independent of each other. This assumption is not true, in particular, if small, neighbouring patches are considered. This causes biases in the covariance estimate. Third, triplets covering two or three patches need to be neglected for the estimation of the data vector. Otherwise, the covariance from jackknife or bootstrap resampling would overestimate the true covariance of the data.

Instead, cosmic shear covariances can also be calculated analytically for a given estimator. This was done, for example, for two-point statistics (Joachimi et al. 2008) or the matter bispectrum (Joachimi et al. 2009). While the derivation of an expression for a covariance can be tedious, once it is done, evaluating them for the characteristics of a specific survey is usually accurate and quick, at least compared to running hundreds of cosmological simulations. It is also easier to infer the dependence of the covariance on survey area and geometry from an analytic expression. Consequently, we follow this approach and derive an analytical expression for the covariance of $\langle M_{\text{ap}}^3\rangle$. We test our result by comparing it with the covariance measured in simulated mock data, both simple Gaussian shear fields and full cosmological simulations, and investigate the influence of different terms in the covariance by performing mock cosmological analyses.

In our derivation, we make five crucial findings. First, the covariance consists of six terms, two of which belong to the Gaussian part and four to the non-Gaussian part. These terms all depend on the survey area and window function (Sects. 3.1 and 3.2). Second, the individual covariance terms can be estimated from correlation functions of the aperture mass map. These correlation functions need to be known only on scales inside the survey area (Sect. 3.3). Third, under the ‘large-field approximation’, which assumes a broad survey window function, two covariance terms vanish, while the others scale inversely with the survey area. The vanishing terms decrease faster than 1/A with survey area A. We refer to them as finite-field terms. One of these terms is already present for Gaussian fields (Sect. 3.4). Fourth, the covariance of $\langle M_{\text{ap}}^3\rangle$ cannot be obtained from the bispectrum covariance unless the large-field approximation is valid. Using a linear transformation of the covariance of a Fourier space quantity to obtain the covariance for a real-space observable has been done for second-order statistics (Joachimi et al. 2021; Friedrich et al. 2021). However, for $\langle M_{\text{ap}}^3\rangle$ this approach already fails for Gaussian fields if the large-field approximation is not used (Sect. 4). Finally, our covariance model, based on an estimator using convergence maps, can be used to obtain upper and lower limits on the covariance for an estimator based on third-order shear correlation functions (Sect. 5.3.2).

This paper is structured as follows. In Sect. 2 we give a short overview of third-order shear statistics, and we introduce our notation, the third-order aperture mass $\langle M_{\text{ap}}^3\rangle$, and its estimator. We derive the covariance of this estimator in Sect. 3. In Sect. 4 we show that an alternative derivation of the $\langle M_{\text{ap}}^3\rangle$ covariance from the bispectrum covariance is not correct. We validate the model by comparing its predictions to simulated data in Sect. 5. Finally, in Sect. 6, we perform mock cosmological analyses to compare the results from the validation data to our model and investigate the impact of the individual covariance terms. We conclude with a summary and discussion in Sect. 7. Throughout this paper, we assume a spatially flat universe. Our covariance modelling code is publically available¹.

2. Theoretical background

In this section, we briefly review some fundamental weak lensing quantities relevant to third-order shear statistics and the covariance calculation in the remainder of this work. We refer to Bartelmann & Schneider (2001), Hoekstra & Jain (2008) or Bartelmann (2010) for in-depth reviews on weak lensing.

One of the fundamental quantities in weak gravitational lensing is the convergence κ(ϑ), which is the normalised surface mass density at angular position ϑ. In a flat universe, it is related to the density contrast δ(χϑ, χ) at angular position ϑ and comoving distance χ via

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

with the Hubble constant H₀, the matter density parameter Ω_m, the speed-of-light c, the cosmic scale factor a(χ) at χ, normalised to unity today, and the probability distribution p(χ) dχ of source galaxies in comoving distance.

We treat κ(ϑ) as a homogenous and isotropic random field, which is characterised by the full set of its polyspectra 𝒫_n, defined by

$$\begin{array}{c}\hfill {\langle \stackrel{\sim }{\kappa }({\mathit{\ell }}_1)\cdots \stackrel{\sim }{\kappa }({\mathit{\ell }}_n)\rangle }_{\mathrm{c}}={\mathcal{P}}_n({\mathit{\ell }}_1,\cdots ,{\mathit{\ell }}_n){(2\pi )}^2{\delta }_{\mathrm{D}}({\mathit{\ell }}_1+\cdots +{\mathit{\ell }}_n),\end{array}$$(2)

where $\stackrel{\sim }{\kappa }(\mathit{\ell })$ is the Fourier transform of κ and the ⟨…⟩_c are connected correlation functions. Interesting for us in this work is the powerspectrum P, bispectrum B, trispectrum T, and the pentaspectrum P₆, defined as

$$\begin{array}{c}\hfill P(\ell )=P_2(\mathit{\ell },-\mathit{\ell }),B({\mathit{\ell }}_1,{\mathit{\ell }}_2)={\mathcal{P}}_n({\mathit{\ell }}_1,{\mathit{\ell }}_2,-{\mathit{\ell }}_1-{\mathit{\ell }}_2),T({\mathit{\ell }}_1,{\mathit{\ell }}_2,{\mathit{\ell }}_3)={\mathcal{P}}_n({\mathit{\ell }}_1,{\mathit{\ell }}_2,{\mathit{\ell }}_3,-{\mathit{\ell }}_1-{\mathit{\ell }}_2-{\mathit{\ell }}_3),\end{array}$$(3)

and

$$\begin{array}{cc}& P_6({\mathit{\ell }}_1,{\mathit{\ell }}_2,{\mathit{\ell }}_3,{\mathit{\ell }}_4,{\mathit{\ell }}_5)={\mathcal{P}}_n({\mathit{\ell }}_1,{\mathit{\ell }}_2,{\mathit{\ell }}_3,{\mathit{\ell }}_4,{\mathit{\ell }}_5,-{\mathit{\ell }}_1-{\mathit{\ell }}_2-{\mathit{\ell }}_3-{\mathit{\ell }}_4-{\mathit{\ell }}_5).\hfill \end{array}$$(4)

The κ-polyspectra of are related to the polyspectra ${\mathcal{P}}_n^{(\text{3d})}$ of the three-dimensional density contrast δ, using the Limber-approximation (Kaiser & Jaffe 1997; Kayo et al. 2013),

$$\begin{array}{cc}\hfill {\mathcal{P}}_n({\mathit{\ell }}_1,\cdots ,{\mathit{\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{(\chi )}^n}{\mathcal{P}}_n^{(3\mathrm{d})}({\mathit{\ell }}_1/\chi ,\cdots ,{\mathit{\ell }}_n/\chi ,\chi ).}\hfill \end{array}$$(5)

We model the three-dimensional power spectrum with the revised Halofit prescription by Takahashi et al. (2012) and the three-dimensional bispectrum with BiHalofit (Takahashi et al. 2020). For the tri- and pentaspectrum, we use the halo model (Cooray & Sheth 2002) with a Sheth–Tormen halo mass function and halo bias (Sheth & Tormen 1999) and Navarro–Frenk–White-halo profiles. To simplify our calculation, we use only the 1-halo term for the trispectrum, which we expect to be dominant at the considered ℓ-scales. For the pentaspectrum, we use the 1- and part of the 2-halo term (see Sect. 3.4 and Appendix C).

The convergence cannot be directly observed, but it is related to the weak lensing shear γ, which describes the change in the observed ellipticity of source galaxies due to the lensing effect. For third-order shear statistics, we are interested in the three-point correlation function of the shear, whose so-called natural components Γ_i were derived by Schneider & Lombardi (2003). However, these components are challenging to handle in practice since they require a large number of bins (typically of the order of 10³) and are complex to model because they relate to the matter bispectrum via multi-dimensional integrals involving oscillating Bessel functions. As shown in Paper I, more practical quantities with similar information content are the third-order aperture statistics $\langle M_{\text{ap}}^3\rangle$.

Aperture statistics are moments of the aperture mass M_ap, a smoothed convergence map, which depends on its position α and a characteristic smoothing scale θ,

$$\begin{array}{c}\hfill M_{\mathrm{ap}}(\mathit{\alpha },\theta )={\displaystyle \int }{\mathrm{d}}^2\vartheta U_{\theta }(|\mathit{\alpha }-\mathit{\vartheta }|)\kappa (\mathit{\vartheta }),\mathrm{with}{U}_{\theta }(\vartheta )=\frac{1}{{\theta }^2}u\left(\frac{\vartheta }{\theta }\right),\end{array}$$(6)

where U_θ is a filter function with aperture scale radius θ. As long as this is a compensated filter function, meaning ∫dϑϑ U_θ(ϑ) = 0, the aperture mass can be estimated from the tangential shear γ_t instead of the convergence (Schneider 1996),

$$\begin{array}{c}\hfill M_{\mathrm{ap}}(\mathit{\alpha },\theta )={\displaystyle \int {\mathrm{d}}^2\vartheta Q_{\theta }(|\mathit{\alpha }-\mathit{\vartheta }|){\gamma }_{\mathrm{t}}(\mathit{\alpha };\mathit{\vartheta }),\mathrm{with}{Q}_{\theta }(\vartheta )=\frac{2}{{\vartheta }^2}{\int }_0^{\vartheta }\mathrm{d}{\vartheta }^{\prime }{\vartheta }^{\prime }{U}_{\theta }({\vartheta }^{\prime })-U_{\theta }(\vartheta ).}\end{array}$$(7)

Here, we are interested in the covariance of the third-order aperture statistics $\langle M_{\text{ap}}^3\rangle$, given by

$$\begin{array}{c}\hfill \langle M_{\mathrm{ap}}^3\rangle ({\theta }_1,{\theta }_2,{\theta }_3)=\langle M_{\mathrm{ap}}(\mathit{\alpha },{\theta }_1)M_{\mathrm{ap}}(\mathit{\alpha },{\theta }_2)M_{\mathrm{ap}}(\mathit{\alpha },{\theta }_3)\rangle ,\end{array}$$(8)

which is our main observable for third-order shear statistics. $\langle M_{\text{ap}}^3\rangle$ is related to the convergence bispectrum by

$$\begin{array}{cc}\hfill \langle M_{\mathrm{ap}}^3\rangle ({\theta }_1,{\theta }_2,{\theta }_3)& ={\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)\stackrel{\sim }{u}(|{\mathit{\ell }}_1+{\mathit{\ell }}_2|{\theta }_3)B({\mathit{\ell }}_1,{\mathit{\ell }}_2).}\hfill \end{array}$$(9)

For the comparison of the analytical covariance to mock data in Sect. 5 and the estimate of the cosmological parameter analysis in Sect. 6 we choose the exponential filter function from Crittenden et al. (2002),

$$\begin{array}{c}\hfill u(x)=\frac{1}{2\pi }(1-\frac{x^2}{2})exp(-\frac{x^2}{2}).\end{array}$$(10)

The normalisation of the aperture filter functions is chosen such that ∫d²ϑQ(ϑ) = 1.

Based on Eq. (6), M_ap(α, θ) can be estimated from convergence fields using

$$\begin{array}{c}\hfill {\widehat{M}}_{\mathrm{ap}}(\mathit{\alpha },\theta )={\displaystyle {\int }_{A^{\prime }}{\mathrm{d}}^2\vartheta U_{\theta }(|\mathit{\alpha }-\mathit{\vartheta }|)\kappa (\mathit{\vartheta }),}\end{array}$$(11)

where A′ is the full area of the convergence field κ. This estimator can be visualised as placing an aperture of radius θ centred at position α and averaging the convergence values within the aperture weighted by the filter function U_θ (see Fig. 1). The $\langle M_{\text{ap}}^3\rangle ({\theta }_1,{\theta }_2,{\theta }_3)$ can then be estimated by averaging the product of the ${\widehat{M}}_{\text{ap}}(\mathit{\alpha },\theta )$ for the three filter radii over the aperture positions α. However, for α close to the border of the survey area, ${\widehat{M}}_{\text{ap}}(\mathit{\alpha },\theta )$ is biased because the aperture includes regions outside the survey where κ is not known. Consequently, for an unbiased estimate of $\langle M_{\text{ap}}^3\rangle$, we only average over an area A smaller than A′, which leads to the estimator

$$\begin{array}{cc}\hfill {\widehat{M}}_{\mathrm{ap}}^3({\theta }_1,{\theta }_2,{\theta }_3)& =\frac{1}{A}{\displaystyle {\int }_A{\mathrm{d}}^2\alpha \prod _{i=1}^3{\int }_{A^{\prime }}{\mathrm{d}}^2{\mathit{\vartheta }}_i{U}_{{\theta }_i}(|\mathit{\alpha }-{\mathit{\vartheta }}_i|)\kappa ({\mathit{\vartheta }}_i).}\hfill \end{array}$$(12)

Fig. 1. Illustration of aperture mass estimation. The area A′ is the size of the full convergence field, on which we place apertures with scale radius θ, illustrated by the circles, to obtain Map(α, θ). Apertures centred on positions α within the smaller area A lie completely within A′, while apertures centred on positions α′ outside of A extend outside of A′, so Map(α′,θ) is biased.

We have assumed here that apertures are placed densely so that the separation between aperture centres is much smaller than the aperture radius θ, to use a continuous integral over α instead of a sum. If A′ is large enough that all U_θ centred on points within A vanish outside of A′, we can replace A′ with the full ℝ². We also introduce the survey window function W_A(α), which is unity for α inside A and zero otherwise. With W, we convert integrals over A into integrals over ℝ², which we denote as two-dimensional integrals without integration borders in the following. The estimator for $\langle M_{\text{ap}}^3\rangle$ becomes

$$\begin{array}{cc}\hfill {\widehat{M}}_{\mathrm{ap}}^3({\theta }_1,{\theta }_2,{\theta }_3)& =\frac{1}{A}{\displaystyle \int {\mathrm{d}}^2\alpha W_A(\mathit{\alpha })[\prod _{i=1}^3\int {\mathrm{d}}^2{\mathit{\vartheta }}_i{U}_{{\theta }_i}(|\mathit{\alpha }-{\mathit{\vartheta }}_i|)\kappa ({\mathit{\vartheta }}_i)].}\hfill \end{array}$$(13)

This estimator is easily applied to simulations, where convergence maps without masks are available (see Sect. 5.2.1). However, this estimator should not be used for survey data that includes masked areas. In that case, one would instead measure the correlation functions Γ_i and convert these into $\langle M_{\text{ap}}^3\rangle$ (see the discussion in Sect. 5.3 of Heydenreich et al. 2023). In this approach, one does not need to decrease the survey area from A′ to A because no apertures are laid down. Therefore, the covariance for this estimator is smaller than the covariance for the ${{\widehat{M}}_{\text{ap}}}^3$ in Eq. (13). However, the unbiased construction of the aperture mass on the whole area A′ requires shear estimates outside of A′ as well. Therefore, even Γ_i measured from shear estimates on all of A′ cannot include the full information that would be contained in an unbiased estimate of ${{\widehat{M}}_{\text{ap}}}^3$ on A′. Consequently, in the absence of masked regions, we expect the magnitude of the covariance for the estimator based on the Γ_i to be between the covariance for ${{\widehat{M}}_{\text{ap}}}^3$ for a survey area of A and a survey area of A′. We verify this expectation in Sect. 5.

3. Derivation of aperture statistics covariance from the real-space estimator

We now derive the covariance $C_{{{\widehat{M}}_{\text{ap}}}^3}$ of ${{\widehat{M}}_{\text{ap}}}^3$. An overview of our calculation is given in Fig. 2, with references to the main equations in our derivation. In the derivation, we find that the covariance comprises several terms, including various permutations of scale radii. For simplicity, we write here ‘Perm’. to indicate permutations; the complete list of permutations for all terms is given in Appendix A. We are not explicitly addressing the effect of shape noise in this section, but shape noise can be easily included in the derived expressions by replacing the power spectrum P with $P+{\sigma }_{ϵ}^2/2n$, where n is the galaxy number density and ${\sigma }_{ϵ}^2$ the two-component ellipticity dispersion. We show the validity of this treatment in Appendix B.

Fig. 2. Schematic representation of the calculation of the covariance $C_{{{\widehat{M}}_{\text{ap}}}^3}$. The covariance is given by the difference between $\langle {{\widehat{M}}_{\text{ap}}}^3{{\widehat{M}}_{\text{ap}}}^3\rangle$ and $\langle {{\widehat{M}}_{\text{ap}}}^3\rangle \langle {{\widehat{M}}_{\text{ap}}}^3\rangle$, the first of which can be decomposed (indicated by solid arrows) into one Gaussian (denoted by G) and three non-Gaussian parts (denoted by BB, PT, and P6). We discuss the Gaussian part in Sect. 3.1 and the non-Gaussian parts in Sect. 3.2. These parts can be further decomposed into terms depending on different permutations of the aperture scale radii, called TPPP, 1 to TP6. For large survey areas, the Ti can be approximated, indicated by dashed arrows, as shown in Sect. 3.4. Under this approximation, two terms vanish, which is why we term them ‘finite-field terms’. Also shown are equation numbers for the relevant expressions.

The covariance is

$$\begin{array}{cc}\hfill C_{{\widehat{M}}_{\mathrm{ap}}^3}({\mathrm{\Theta }}_1,{\mathrm{\Theta }}_2)& =\langle {\widehat{M}}_{\mathrm{ap}}^3{\widehat{M}}_{\mathrm{ap}}^3\rangle ({\mathrm{\Theta }}_1,{\mathrm{\Theta }}_2)-\langle {\widehat{M}}_{\mathrm{ap}}^3({\mathrm{\Theta }}_1)\rangle \langle {\widehat{M}}_{\mathrm{ap}}^3({\mathrm{\Theta }}_2)\rangle ,\hfill \end{array}$$(14)

where Θ₁ = (θ₁, θ₂, θ₃) and Θ₂ = (θ₄, θ₅, θ₆). With Eq. (13) and κ_i := κ(ϑ_i),

$$\begin{array}{cc}\hfill \langle {\widehat{M}}_{\mathrm{ap}}^3{\widehat{M}}_{\mathrm{ap}}^3\rangle ({\mathrm{\Theta }}_1,{\mathrm{\Theta }}_2)& =\frac{1}{A^2}{\displaystyle \int {\mathrm{d}}^2{\alpha }_1\int {\mathrm{d}}^2{\alpha }_2{W}_A({\mathit{\alpha }}_1)W_A({\mathit{\alpha }}_2)[\prod _{i=1}^3\int {\mathrm{d}}^2{\vartheta }_i{U}_{{\theta }_i}(|{\mathit{\alpha }}_1-{\mathit{\vartheta }}_i|)][\prod _{j=4}^6\int {\mathrm{d}}^2{\vartheta }_j{U}_{{\theta }_j}(|{\mathit{\alpha }}_2-{\mathit{\vartheta }}_j|)]}\hfill \\ \hfill & \times \langle {\kappa }_1{\kappa }_2{\kappa }_3{\kappa }_4{\kappa }_5{\kappa }_6\rangle .\hfill \end{array}$$(15)

The six-point correlation can be written in terms of connected correlation functions, denoted by ⟨⟩_c. This, together with ⟨κ(ϑ)⟩ = 0, leads to

$$\begin{array}{cc}\hfill \langle {\widehat{M}}_{\mathrm{ap}}^3{\widehat{M}}_{\mathrm{ap}}^3\rangle ({\mathrm{\Theta }}_1,{\mathrm{\Theta }}_2)& =\frac{1}{A^2}{\displaystyle \int {\mathrm{d}}^2{\alpha }_1\int {\mathrm{d}}^2{\alpha }_2{W}_A({\mathit{\alpha }}_1)W_A({\mathit{\alpha }}_2)[\prod _{i=1}^3\int {\mathrm{d}}^2{\vartheta }_i{U}_{{\theta }_i}(|{\mathit{\alpha }}_1-{\mathit{\vartheta }}_i|)][\prod _{j=4}^6\int {\mathrm{d}}^2{\vartheta }_j{U}_{{\theta }_j}(|{\mathit{\alpha }}_2-{\mathit{\vartheta }}_j|)]}\hfill \\ \hfill & \times [({\langle {\kappa }_1{\kappa }_2\rangle }_{\mathrm{c}}{\langle {\kappa }_3{\kappa }_4\rangle }_{\mathrm{c}}{\langle {\kappa }_5{\kappa }_6\rangle }_{\mathrm{c}}+14\mathrm{Perm}.)+({\langle {\kappa }_1{\kappa }_2{\kappa }_3\rangle }_{\mathrm{c}}{\langle {\kappa }_4{\kappa }_5{\kappa }_6\rangle }_{\mathrm{c}}+9\mathrm{Perm}.)\hfill \\ \hfill & +({\langle {\kappa }_1{\kappa }_2\rangle }_{\mathrm{c}}{\langle {\kappa }_3{\kappa }_4{\kappa }_5{\kappa }_6\rangle }_{\mathrm{c}}+14\mathrm{Perm}.)+{\langle {\kappa }_1{\kappa }_2{\kappa }_3{\kappa }_4{\kappa }_5{\kappa }_6\rangle }_{\mathrm{c}}]\hfill \\ \hfill & ={\langle {\widehat{M}}_{\mathrm{ap}}^3{\widehat{M}}_{\mathrm{ap}}^3\rangle }_{\mathrm{G}}({\mathrm{\Theta }}_1,{\mathrm{\Theta }}_2)+{\langle {\widehat{M}}_{\mathrm{ap}}^3{\widehat{M}}_{\mathrm{ap}}^3\rangle }_{\mathit{BB}}({\mathrm{\Theta }}_1,{\mathrm{\Theta }}_2)+{\langle {\widehat{M}}_{\mathrm{ap}}^3{\widehat{M}}_{\mathrm{ap}}^3\rangle }_{\mathit{PT}}({\mathrm{\Theta }}_1,{\mathrm{\Theta }}_2)\hfill \\ \hfill & +{\langle {\widehat{M}}_{\mathrm{ap}}^3{\widehat{M}}_{\mathrm{ap}}^3\rangle }_{P_6}({\mathrm{\Theta }}_1,{\mathrm{\Theta }}_2).\hfill \end{array}$$(16)

The first term ${\langle {{\widehat{M}}_{\text{ap}}}^3{{\widehat{M}}_{\text{ap}}}^3\rangle }_{\mathrm{G}}$ is the Gaussian part of the covariance – the only part that is non-zero for Gaussian random fields. It depends only on the two-point correlation of κ, or, equivalently, the matter power spectrum P(ℓ). The other terms comprise the non-Gaussian part of the covariance. They depend on higher-order polyspectra, namely the bi-, tri- and the pentaspectrum.

3.1. Gaussian part

We first concentrate on the Gaussian part ${\langle {{\widehat{M}}_{\text{ap}}}^3{{\widehat{M}}_{\text{ap}}}^3\rangle }_{\mathrm{G}}$, given by

$$\begin{array}{cc}\hfill {\langle {\widehat{M}}_{\mathrm{ap}}^3{\widehat{M}}_{\mathrm{ap}}^3\rangle }_{\mathrm{G}}({\mathrm{\Theta }}_1,{\mathrm{\Theta }}_2)& =\frac{1}{A^2}{\displaystyle \int {\mathrm{d}}^2{\alpha }_1\int {\mathrm{d}}^2{\alpha }_2{W}_A({\mathit{\alpha }}_1)W_A({\mathit{\alpha }}_2)[\prod _{i=1}^3\int {\mathrm{d}}^2{\vartheta }_i{U}_{{\theta }_i}(|{\mathit{\alpha }}_1-{\mathit{\vartheta }}_i|)][\prod _{j=4}^6\int {\mathrm{d}}^2{\vartheta }_j{U}_{{\theta }_j}(|{\mathit{\alpha }}_2-{\mathit{\vartheta }}_j|)]}\hfill \\ \hfill & \times [{\langle {\kappa }_1{\kappa }_2\rangle }_{\mathrm{c}}{\langle {\kappa }_3{\kappa }_4\rangle }_{\mathrm{c}}{\langle {\kappa }_5{\kappa }_6\rangle }_{\mathrm{c}}+14\mathrm{Perm}.].\hfill \end{array}$$(17)

We split the permutations in Eq. (17) into two groups, as

$$\begin{array}{cc}\hfill {\langle {\widehat{M}}_{\mathrm{ap}}^3{\widehat{M}}_{\mathrm{ap}}^3\rangle }_{\mathrm{G}}({\mathrm{\Theta }}_1,{\mathrm{\Theta }}_2)& =\frac{1}{A^2}{\displaystyle \int {\mathrm{d}}^2{\alpha }_1\int {\mathrm{d}}^2{\alpha }_2{W}_A({\mathit{\alpha }}_1)W_A({\mathit{\alpha }}_2)[\prod _{i=1}^3\int {\mathrm{d}}^2{\vartheta }_i{U}_{{\theta }_i}(|{\mathit{\alpha }}_1-{\mathit{\vartheta }}_i|)][\prod _{j=4}^6\int {\mathrm{d}}^2{\vartheta }_j{U}_{{\theta }_j}(|{\mathit{\alpha }}_2-{\mathit{\vartheta }}_j|)]}\hfill \\ \hfill & \times [{\langle {\kappa }_1{\kappa }_4\rangle }_{\mathrm{c}}{\langle {\kappa }_3{\kappa }_5\rangle }_{\mathrm{c}}{\langle {\kappa }_2{\kappa }_6\rangle }_{\mathrm{c}}+5\mathrm{Perm}.]\hfill \\ \hfill & +\frac{1}{A^2}{\displaystyle \int {\mathrm{d}}^2{\alpha }_1\int {\mathrm{d}}^2{\alpha }_2{W}_A({\mathit{\alpha }}_1)W_A({\mathit{\alpha }}_2)[\prod _{i=1}^3\int {\mathrm{d}}^2{\vartheta }_i{U}_{{\theta }_i}(|{\mathit{\alpha }}_1-{\mathit{\vartheta }}_i|)][\prod _{j=4}^6\int {\mathrm{d}}^2{\vartheta }_j{U}_{{\theta }_j}(|{\mathit{\alpha }}_2-{\mathit{\vartheta }}_j|)]}\hfill \\ \hfill & \times [{\langle {\kappa }_1{\kappa }_2\rangle }_{\mathrm{c}}{\langle {\kappa }_3{\kappa }_4\rangle }_{\mathrm{c}}{\langle {\kappa }_5{\kappa }_6\rangle }_{\mathrm{c}}+8\mathrm{Perm}.]\hfill \\ \hfill & =T_{PPP,1}({\mathrm{\Theta }}_1,{\mathrm{\Theta }}_2)+T_{PPP,2}({\mathrm{\Theta }}_1,{\mathrm{\Theta }}_2).\hfill \end{array}$$(18)

The first group consists of six terms, in which for all three ⟨κ_i κ_j⟩_c, the first index i ∈ {1, 2, 3}, and the second index j ∈ {4, 5, 6}. The second group consists of the nine other permutations. To evaluate T_PPP, 1 and T_PPP, 2, we use Eq. (2) for n = 2 and

$$\begin{array}{c}\hfill {\displaystyle \int }{\mathrm{d}}^2\vartheta U_{\theta }(|\mathit{\alpha }-\mathit{\vartheta }|){\mathrm{e}}^{-\mathrm{i}\mathit{\ell }·\mathit{\vartheta }}=\stackrel{\sim }{u}(\ell \theta ){\mathrm{e}}^{-\mathrm{i}\mathit{\ell }·\mathit{\alpha }},\end{array}$$(19)

and find

$$\begin{array}{cc}\hfill T_{PPP,1}({\mathrm{\Theta }}_1,{\mathrm{\Theta }}_2)& =\frac{1}{A^2}{\displaystyle [\prod _{i=1}^3\int \frac{{\mathrm{d}}^2{\ell }_i}{{(2\pi )}^2}P({\ell }_i)][\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 }_1{\theta }_4)\stackrel{\sim }{u}({\ell }_2{\theta }_5)\stackrel{\sim }{u}({\ell }_3{\theta }_6)+5\mathrm{Perm}.]}\hfill \\ \hfill & \times {\displaystyle \int {\mathrm{d}}^2{\alpha }_1\int {\mathrm{d}}^2{\alpha }_2{W}_A({\mathit{\alpha }}_1)W_A({\mathit{\alpha }}_2){\mathrm{e}}^{-\mathrm{i}({\mathit{\alpha }}_1-{\mathit{\alpha }}_2)({\mathit{\ell }}_1+{\mathit{\ell }}_2+{\mathit{\ell }}_3)}}\hfill \\ \hfill & ={\displaystyle [\prod _{i=1}^3\int \frac{{\mathrm{d}}^2{\ell }_i}{{(2\pi )}^2}P({\ell }_i)]G_A({\mathit{\ell }}_1+{\mathit{\ell }}_2+{\mathit{\ell }}_3)[\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 }_1{\theta }_4)\stackrel{\sim }{u}({\ell }_2{\theta }_5)\stackrel{\sim }{u}({\ell }_3{\theta }_6)+5\mathrm{Perm}.],}\hfill \end{array}$$(20)

and

$$\begin{array}{cc}\hfill T_{PPP,2}({\mathrm{\Theta }}_1,{\mathrm{\Theta }}_2)& =\frac{1}{A^2}{\displaystyle [\prod _{i=1}^3\int \frac{{\mathrm{d}}^2{\ell }_i}{{(2\pi )}^2}P({\ell }_i)][\stackrel{\sim }{u}({\ell }_1{\theta }_1)\stackrel{\sim }{u}({\ell }_1{\theta }_2)\stackrel{\sim }{u}({\ell }_2{\theta }_3)\stackrel{\sim }{u}({\ell }_2{\theta }_4)\stackrel{\sim }{u}({\ell }_3{\theta }_5)\stackrel{\sim }{u}({\ell }_3{\theta }_6)+8\mathrm{Perm}.]}\hfill \\ \hfill & \times {\displaystyle \int {\mathrm{d}}^2{\alpha }_1\int {\mathrm{d}}^2{\alpha }_2{W}_A({\mathit{\alpha }}_1)W_A({\mathit{\alpha }}_2){\mathrm{e}}^{-\mathrm{i}({\mathit{\alpha }}_1-{\mathit{\alpha }}_2)·{\mathit{\ell }}_2}}\hfill \\ \hfill & ={\displaystyle [\prod _{i=1}^3\int \frac{{\mathrm{d}}^2{\ell }_i}{{(2\pi )}^2}P({\ell }_i)]G_A({\mathit{\ell }}_2)[\stackrel{\sim }{u}({\ell }_1{\theta }_1)\stackrel{\sim }{u}({\ell }_1{\theta }_2)\stackrel{\sim }{u}({\ell }_2{\theta }_3)\stackrel{\sim }{u}({\ell }_2{\theta }_4)\stackrel{\sim }{u}({\ell }_3{\theta }_5)\stackrel{\sim }{u}({\ell }_3{\theta }_6)+8\mathrm{Perm}.],}\hfill \end{array}$$(21)

where the geometry factor G_A(ℓ), defined by

$$\begin{array}{c}\hfill G_A(\mathit{\ell })=\frac{1}{A^2}{\displaystyle \int {\mathrm{d}}^2{\alpha }_1\int {\mathrm{d}}^2{\alpha }_2{W}_A({\mathit{\alpha }}_1)W_A({\mathit{\alpha }}_2){\mathrm{e}}^{-\mathrm{i}({\mathit{\alpha }}_1-{\mathit{\alpha }}_2)·\mathit{\ell }},}\end{array}$$(22)

contains the full dependence of the covariance on the survey area and geometry. For a square survey with side length ${\vartheta }_{\text{max}}=\sqrt{A}$, and ℓ = (ℓ_x, ℓ_y) the factor is

$$\begin{array}{c}\hfill G_{A,\mathrm{square}}(\mathit{\ell })=\frac{4{sin}^2({\ell }_x{\vartheta }_{\max }/2)}{{\ell }_x^2{\vartheta }_{\max }^2}\frac{4{sin}^2({\ell }_y{\vartheta }_{\max }/2)}{{\ell }_y^2{\vartheta }_{\max }^2}.\end{array}$$(23)

The geometry factor G_A is related to a function E_A(η), which, for a point α within A, gives the probability for a second point at α + η to lie within A as well. It is given by

$$\begin{array}{c}\hfill E_A(\mathit{\eta })=\frac{1}{A}{\displaystyle \int {\mathrm{d}}^2\alpha W_A(\mathit{\alpha })W_A(\mathit{\alpha }+\mathit{\eta }),\mathrm{so}\mathrm{that}{G}_A(\mathit{\ell })=\frac{1}{A}{\int }_A{\mathrm{d}}^2\eta E_A(\mathit{\eta }){\mathrm{e}}^{-\mathrm{i}\mathit{\eta }·\mathit{\ell }}.}\end{array}$$(24)

For a square area, E_A(η) was determined by Heydenreich et al. (2020). It is unity for vanishing separation |η| = 0 and smoothly declines to zero for η outside the square.

3.2. Non-Gaussian part

We now derive the non-Gaussian part of the covariance. First, we calculate the covariance part ${\langle {{\widehat{M}}_{\text{ap}}}^3{{\widehat{M}}_{\text{ap}}}^3\rangle }_{\mathit{BB}}$, which depends on third-order correlations of κ. We divide this term into two parts, as

$$\begin{array}{cc}\hfill {\langle {\widehat{M}}_{\mathrm{ap}}^3{\widehat{M}}_{\mathrm{ap}}^3\rangle }_{\mathit{BB}}({\mathrm{\Theta }}_1,{\mathrm{\Theta }}_2)& =\frac{1}{A^2}{\displaystyle \int {\mathrm{d}}^2{\alpha }_1\int {\mathrm{d}}^2{\alpha }_2[\prod _{i=1}^3\int {\mathrm{d}}^2{\vartheta }_i{U}_{{\theta }_i}(|{\mathit{\alpha }}_1-{\mathit{\vartheta }}_i|)][\prod _{j=4}^6\int {\mathrm{d}}^2{\vartheta }_j{U}_{{\theta }_j}(|{\mathit{\alpha }}_2-{\mathit{\vartheta }}_j|)]}\hfill \\ \hfill & \times W_A({\mathit{\alpha }}_1)W_A({\mathit{\alpha }}_2)\{{\langle {\kappa }_1{\kappa }_2{\kappa }_3\rangle }_{\mathrm{c}}{\langle {\kappa }_4{\kappa }_5{\kappa }_6\rangle }_{\mathrm{c}}+[{\langle {\kappa }_4{\kappa }_2{\kappa }_3\rangle }_{\mathrm{c}}{\langle {\kappa }_1{\kappa }_5{\kappa }_6\rangle }_{\mathrm{c}}+8\mathrm{Perm}.]\}\hfill \\ \hfill & =\langle {\widehat{M}}_{\mathrm{ap}}^3({\theta }_1,{\theta }_2,{\theta }_3)\rangle \langle {\widehat{M}}_{\mathrm{ap}}^3({\theta }_4,{\theta }_5,{\theta }_6)\rangle +T_{\mathit{BB}}({\mathrm{\Theta }}_1,{\mathrm{\Theta }}_2).\hfill \end{array}$$(25)

The first term cancels with the second term in the definition of the covariance of the estimator in Eq. (14). The term T_BB is

$$\begin{array}{cc}\hfill T_{\mathit{BB}}({\mathrm{\Theta }}_1,{\mathrm{\Theta }}_2)& =\frac{1}{A^2}{\displaystyle \int {\mathrm{d}}^2{\alpha }_1\int {\mathrm{d}}^2{\alpha }_2[\prod _{i=1}^3\int {\mathrm{d}}^2{\vartheta }_i{U}_{{\theta }_i}(|{\mathit{\alpha }}_1-{\mathit{\vartheta }}_i|)][\prod _{j=4}^6\int {\mathrm{d}}^2{\vartheta }_j{U}_{{\theta }_j}(|{\mathit{\alpha }}_2-{\mathit{\vartheta }}_j|)]W_A({\mathit{\alpha }}_1)W_A({\mathit{\alpha }}_2)}\hfill \\ \hfill & \times [{\langle {\kappa }_4{\kappa }_2{\kappa }_3\rangle }_{\mathrm{c}}{\langle {\kappa }_1{\kappa }_5{\kappa }_6\rangle }_{\mathrm{c}}+8\mathrm{Perm}.].\hfill \end{array}$$(26)

Using Eq. (2) for n = 3 to introduce the bispectrum and the definition of G_A in Eq. (22) leads to

$$\begin{array}{cc}\hfill T_{\mathit{BB}}({\mathrm{\Theta }}_1,{\mathrm{\Theta }}_2)& ={\displaystyle \int \frac{{\mathrm{d}}^2{\ell }_1}{{(2\pi )}^2}\int \frac{{\mathrm{d}}^2{\ell }_2}{{(2\pi )}^2}\int \frac{{\mathrm{d}}^2{\ell }_4}{{(2\pi )}^2}\int \frac{{\mathrm{d}}^2{\ell }_5}{{(2\pi )}^2}B({\mathit{\ell }}_1,{\mathit{\ell }}_2)B({\mathit{\ell }}_4,{\mathit{\ell }}_5)G_A({\mathit{\ell }}_1-{\mathit{\ell }}_4)}\hfill \\ \hfill & \times [\stackrel{\sim }{u}({\ell }_4{\theta }_1)\stackrel{\sim }{u}({\ell }_2{\theta }_2)\stackrel{\sim }{u}(|{\mathit{\ell }}_1+{\mathit{\ell }}_2|{\theta }_3)\stackrel{\sim }{u}({\ell }_1{\theta }_4)\stackrel{\sim }{u}({\ell }_5{\theta }_5)\stackrel{\sim }{u}(|{\mathit{\ell }}_4+{\mathit{\ell }}_5|{\theta }_6)+8\mathrm{Perm}.].\hfill \end{array}$$(27)

Second, we calculate the covariance part ${\langle {{\widehat{M}}_{\text{ap}}}^3{{\widehat{M}}_{\text{ap}}}^3\rangle }_{\mathit{PT}}$. We divide this term into two parts as

$$\begin{array}{cc}\hfill {\langle {\widehat{M}}_{\mathrm{ap}}^3{\widehat{M}}_{\mathrm{ap}}^3\rangle }_{\mathit{PT}}({\mathrm{\Theta }}_1,{\mathrm{\Theta }}_2)& =\frac{1}{A^2}{\displaystyle \int {\mathrm{d}}^2{\alpha }_1\int {\mathrm{d}}^2{\alpha }_2[\prod _{i=1}^3\int {\mathrm{d}}^2{\vartheta }_i{U}_{{\theta }_i}(|{\mathit{\alpha }}_1-{\mathit{\vartheta }}_i|)][\prod _{j=4}^6\int {\mathrm{d}}^2{\vartheta }_j{U}_{{\theta }_j}(|{\mathit{\alpha }}_2-{\mathit{\vartheta }}_j|)]}\hfill \\ \hfill & \times W_A({\mathit{\alpha }}_1)W_A({\mathit{\alpha }}_2)[{\langle {\kappa }_1{\kappa }_4\rangle }_{\mathrm{c}}{\langle {\kappa }_2{\kappa }_3{\kappa }_5{\kappa }_6\rangle }_{\mathrm{c}}+8\mathrm{Perm}.+{\langle {\kappa }_1{\kappa }_2\rangle }_{\mathrm{c}}{\langle {\kappa }_3{\kappa }_4{\kappa }_5{\kappa }_6\rangle }_{\mathrm{c}}+5\mathrm{Perm}.]\hfill \\ \hfill & =T_{PT,1}({\mathrm{\Theta }}_1,{\mathrm{\Theta }}_2)+T_{PT,2}({\mathrm{\Theta }}_1,{\mathrm{\Theta }}_2).\hfill \end{array}$$(28)

Here, the first permutations contain all terms where the two-point correlation ⟨κ_i κ_j⟩ has i ∈ {1, 2, 3} and j ∈ {4, 5, 6}. The second permutations contain all other terms. We introduce the power- and trispectrum with Eq. (2), and find

$$\begin{array}{cc}\hfill T_{PT,1}({\mathrm{\Theta }}_1,{\mathrm{\Theta }}_2)& ={\displaystyle \int \frac{{\mathrm{d}}^2{\ell }_1}{{(2\pi )}^2}\int \frac{{\mathrm{d}}^2{\ell }_2}{{(2\pi )}^2}\int \frac{{\mathrm{d}}^2{\ell }_3}{{(2\pi )}^2}\int \frac{{\mathrm{d}}^2{\ell }_4}{{(2\pi )}^2}P({\ell }_1)T({\mathit{\ell }}_2,{\mathit{\ell }}_3,{\mathit{\ell }}_4)G_A({\mathit{\ell }}_1+{\mathit{\ell }}_2+{\mathit{\ell }}_3)}\hfill \\ \hfill & \times [\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 }_1{\theta }_4)\stackrel{\sim }{u}({\ell }_4{\theta }_5)\stackrel{\sim }{u}(|{\mathit{\ell }}_2+{\mathit{\ell }}_3+{\mathit{\ell }}_4|{\theta }_6)+8\mathrm{Perm}.],\hfill \end{array}$$(29)

and

$$\begin{array}{cc}\hfill T_{PT,2}({\mathrm{\Theta }}_1,{\mathrm{\Theta }}_2)& ={\displaystyle \int \frac{{\mathrm{d}}^2{\ell }_1}{{(2\pi )}^2}\int \frac{{\mathrm{d}}^2{\ell }_2}{{(2\pi )}^2}\int \frac{{\mathrm{d}}^2{\ell }_3}{{(2\pi )}^2}\int \frac{{\mathrm{d}}^2{\ell }_4}{{(2\pi )}^2}P({\ell }_1)T({\mathit{\ell }}_2,{\mathit{\ell }}_3,{\mathit{\ell }}_4)G_A({\mathit{\ell }}_2)}\hfill \\ \hfill & \times [\stackrel{\sim }{u}({\ell }_1{\theta }_1)\stackrel{\sim }{u}({\ell }_1{\theta }_2)\stackrel{\sim }{u}({\ell }_2{\theta }_3)\stackrel{\sim }{u}({\ell }_3{\theta }_4)\stackrel{\sim }{u}({\ell }_4{\theta }_5)\stackrel{\sim }{u}(|{\mathit{\ell }}_2+{\mathit{\ell }}_3+{\mathit{\ell }}_4|{\theta }_6)+5\mathrm{Perm}.].\hfill \end{array}$$(30)

Finally, we consider ${\langle {{\widehat{M}}_{\text{ap}}}^3{{\widehat{M}}_{\text{ap}}}^3\rangle }_{P_6}$, which is

$$\begin{array}{cc}\hfill {\langle {\widehat{M}}_{\mathrm{ap}}^3{\widehat{M}}_{\mathrm{ap}}^3\rangle }_{P_6}({\mathrm{\Theta }}_1,{\mathrm{\Theta }}_2)& =\frac{1}{A^2}{\displaystyle \int {\mathrm{d}}^2{\alpha }_1\int {\mathrm{d}}^2{\alpha }_2{W}_A({\mathit{\alpha }}_1)W_A({\mathit{\alpha }}_2)}\hfill \\ \hfill & \times [{\displaystyle \prod _{i=1}^3\int {\mathrm{d}}^2{\vartheta }_i{U}_{{\theta }_i}(|{\mathit{\alpha }}_1-{\mathit{\vartheta }}_i|)][\prod _{j=4}^6\int {\mathrm{d}}^2{\vartheta }_j{U}_{{\theta }_j}(|{\mathit{\alpha }}_2-{\mathit{\vartheta }}_j|)]{\langle {\kappa }_1{\kappa }_2{\kappa }_3{\kappa }_4{\kappa }_5{\kappa }_6\rangle }_{\mathrm{c}}}\hfill \\ \hfill & =T_{P_6}({\mathrm{\Theta }}_1,{\mathrm{\Theta }}_2).\hfill \end{array}$$(31)

The term T_P6 is, with Eq. (2) for n = 6,

$$\begin{array}{cc}\hfill T_{P_6}({\mathrm{\Theta }}_1,{\mathrm{\Theta }}_2)& ={\displaystyle \int \frac{{\mathrm{d}}^2{\ell }_1}{{(2\pi )}^2}\int \frac{{\mathrm{d}}^2{\ell }_2}{{(2\pi )}^2}\int \frac{{\mathrm{d}}^2{\ell }_3}{{(2\pi )}^2}\int \frac{{\mathrm{d}}^2{\ell }_4}{{(2\pi )}^2}\int \frac{{\mathrm{d}}^2{\ell }_5}{{(2\pi )}^2}{P}_6({\mathit{\ell }}_1,{\mathit{\ell }}_2,{\mathit{\ell }}_3,{\mathit{\ell }}_4,{\mathit{\ell }}_5)G_A({\mathit{\ell }}_1+{\mathit{\ell }}_2+{\mathit{\ell }}_3)}\hfill \\ \hfill & \times \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)\stackrel{\sim }{u}({\ell }_5{\theta }_5)\stackrel{\sim }{u}(|{\mathit{\ell }}_1+{\mathit{\ell }}_2+{\mathit{\ell }}_3+{\mathit{\ell }}_4+{\mathit{\ell }}_5|{\theta }_6).\hfill \end{array}$$(32)

3.3. Connection of covariance terms to aperture mass correlation functions

We now show that the individual covariance terms are related to correlation functions of the aperture mass. While this might appear like a meaningless algebraic exercise, it has an important corollary: Since the correlation functions can be measured, we can estimate all individual terms of the covariance in the validation data. Thereby, the analytic expressions can be validated.

We define the nth order aperture mass correlation functions $M_{\text{ap}}^{n,m}$, as

$$\begin{array}{c}\hfill M_{\mathrm{ap}}^{n,m}({\theta }_1,\cdots ,{\theta }_m;{\theta }_{m+1},\cdots ,{\theta }_n;\mathit{\eta })=\langle M_{\mathrm{ap}}(\mathit{\vartheta },{\theta }_1)\cdots M_{\mathrm{ap}}(\mathit{\vartheta },{\theta }_m)M_{\mathrm{ap}}(\mathit{\vartheta }+\mathit{\eta },{\theta }_{m+1})M_{\mathrm{ap}}(\mathit{\vartheta }+\mathit{\eta },{\theta }_n)\rangle .\end{array}$$(33)

This function is symmetric in its first m aperture radii and in its second n − m aperture radii. The correlation functions can be expressed in terms of the convergence as

$$\begin{array}{c}\hfill M_{\mathrm{ap}}^{n,m}({\theta }_1,\cdots ,{\theta }_m;{\theta }_{m+1},\cdots ,{\theta }_n;\mathit{\eta })=[{\displaystyle \prod _{i=1}^m\int {\mathrm{d}}^2{\vartheta }_i{U}_{{\theta }_i}({\vartheta }_i)][\prod _{j=m+1}^n\int {\mathrm{d}}^2{\vartheta }_j{U}_{{\theta }_j}(|{\mathit{\vartheta }}_j+\mathit{\eta }|)]\langle {\kappa }_1\cdots {\kappa }_n\rangle ,}\end{array}$$(34)

where κ_i = κ(ϑ_i) and the average over the convergence contains both connected and unconnected terms. Using n = 2, Eq. (18) and the definition of E(η) from Eq. (24), the first Gaussian covariance term can be written as

$$\begin{array}{cc}\hfill T_{PPP,1}({\mathrm{\Theta }}_1,{\mathrm{\Theta }}_2)& =\frac{1}{A^2}{\displaystyle \int {\mathrm{d}}^2{\alpha }_1\int {\mathrm{d}}^2{\alpha }_2{W}_A({\mathit{\alpha }}_1)W_A({\mathit{\alpha }}_2)M_{\mathrm{ap}}^{2,1}({\theta }_1;{\theta }_4;|{\mathit{\alpha }}_1-{\mathit{\alpha }}_2|)M_{\mathrm{ap}}^{2,1}({\theta }_2;{\theta }_5;|{\mathit{\alpha }}_1-{\mathit{\alpha }}_2|)}\hfill \\ \hfill & \times M_{\mathrm{ap}}^{2,1}({\theta }_3;{\theta }_6;|{\mathit{\alpha }}_1-{\mathit{\alpha }}_2|)+5\mathrm{Perm}.\hfill \\ \hfill & =\frac{1}{A}{\displaystyle {\int }_A{\mathrm{d}}^2\eta E_A(\mathit{\eta })M_{\mathrm{ap}}^{2,1}({\theta }_1;{\theta }_4;\mathit{\eta })M_{\mathrm{ap}}^{2,1}({\theta }_2;{\theta }_5;\mathit{\eta })M_{\mathrm{ap}}^{2,1}({\theta }_3;{\theta }_6;\mathit{\eta })+5\mathrm{Perm}.}\hfill \end{array}$$(35)