Issue 
A&A
Volume 672, April 2023



Article Number  A44  
Number of page(s)  19  
Section  Cosmology (including clusters of galaxies)  
DOI  https://doi.org/10.1051/00046361/202244820  
Published online  27 March 2023 
A roadmap to cosmological parameter analysis with thirdorder shear statistics
I. Modelling and validation^{⋆}
ArgelanderInstitut für Astronomie, Auf dem Hügel 71, 53121 Bonn, Germany
email: sven@astro.unibonn.de
Received:
27
August
2022
Accepted:
11
November
2022
In this work, which is the first of a series to prepare a cosmological parameter analysis with thirdorder cosmic shear statistics, we model both the shear threepoint correlation functions Γ^{(i)} and the thirdorder aperture statistics from the B IHALOFIT bispectrum model and validate these statistics with a series of Nbody simulations. We then investigate how to bin the shear threepoint correlation functions to achieve an unbiased estimate for thirdorder aperture statistics in real data. Finally, we perform a cosmological parameter analysis on KiDS1000like mock data with second and thirdorder statistics. In the absence of systematic effects, we recover all cosmological parameters with very little bias. Furthermore, we find that a joint analysis almost doubles the constraining power on S_{8} and increases the figure of merit in the Ω_{m}σ_{8} plane by a factor of 5.9 with respect to an analysis with only secondorder shear statistics.
Key words: gravitational lensing: weak / cosmological parameters / largescale structure of Universe
Our modelling pipeline is publicly available at https://github.com/sheydenreich/threepoint/releases/
© The Authors 2023
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.
1. Introduction
The Λ cold dark matter (ΛCDM) model has been considered the standard model of cosmology for the past few decades. This relatively simple, sixparameter model describes a wide range of observations, from the cosmic microwave background (CMB) to the observed largescale structure (LSS) of galaxies, with remarkable accuracy. As the reported uncertainties on cosmological parameters are approaching the percent level, a few tensions arise between CMB observations of the early Universe and observations of the local Universe that quantify the LSS (for example, in the Hubble parameter H_{0}, see Di Valentino et al. 2021a, and references therein). In the past few years, the matter clustering parameter has also become subject to tension (Hildebrandt et al. 2017; Planck Collaboration VI 2020; Joudaki et al. 2020; Heymans et al. 2021; Abbott et al. 2022; Di Valentino et al. 2021b, and references therein): The local Universe seems less clustered than observations of the CMB suggest. Assuming that these tensions are not due to unknown systematic effects, extensions to the ΛCDM model need to be explored. One of the most popular extensions is the wCDM model, where the equation of state of dark energy differs from w = −1.
The dark energy task force has established that the weak gravitational lensing effect from the LSS, also called cosmic shear, poses one of the most promising methods to constrain the equation of state of dark energy (Albrecht et al. 2006). The next generation of cosmic shear surveys such as Euclid (Laureijs et al. 2011) or the Vera Rubin Observatory Legacy Survey of Space and Time (LSST, Ivezic et al. 2008) will be able to constrain potential extensions to the ΛCDM model and may help to decipher the nature of dark energy.
Tight constraints on cosmological parameters are essential to discriminate between the different modifications of the ΛCDM model. So far, twopoint statistics have been established as the main analysis tool for cosmic shear (Schneider et al. 1998; Troxel et al. 2018; Hildebrandt et al. 2017, 2020; Hikage et al. 2019; Asgari et al. 2020). These statistics capture the entire information content of a Gaussian random field. Since the initial density field of the Universe is believed to be Gaussian, twopoint statistics capture a large amount of cosmological information. However, in late times nonlinear structure formation has introduced nonGaussian features at the smaller scales of the matter distribution, whose information content cannot be captured by twopoint statistics. To use this information, a variety of higherorder statistics has been introduced in recent years, including peak count statistics (Martinet et al. 2018; HarnoisDéraps et al. 2021), persistent homology (Heydenreich et al. 2021), density split statistics (Gruen et al. 2018; Burger et al. 2022), and many others. In this work, we consider thirdorder shear statistics, which measure the skewness of the LSS at various scales.
In contrast to most of these higherorder statistics, threepoint statistics in configuration space can be directly modelled from a matter bispectrum, allowing for a broad range of consistency checks that can be performed. Furthermore, their modelling does not require simulations that are adjusted to specific survey properties such as the number density and redshift distribution of sources, which allows us to apply them to any data set easily. However, these natural extensions to twopoint statistics have yet received surprisingly little attention. Several papers have reported a massive potential information gain when combining two and threepoint statistics (Takada & Jain 2004; Kilbinger & Schneider 2005; Sato & Nishimichi 2013; Kayo et al. 2013). Fu et al. (2014) performed a combined analysis of two and threepoint statistics on 154 deg^{2} CFHTLenS data, reporting a rather moderate gain in information content. Secco et al. (2022) measured the shear threepoint correlation functions and thirdorder aperture statistics in the thirdyear data release of the Dark Energy Survey (Flaugher 2005; SevillaNoarbe et al. 2021), showing that they can be detected with high signaltonoise. Huterer et al. (2006), Pyne & Joachimi (2021) showed that a combined analysis of two and threepoint statistics has an additional advantage: These two statistics react differently to observational and astrophysical systematics, meaning that a combined analysis allows us to constrain nuisance parameters internally without the need for any additional observations or simulations, yielding an additional nulltest and much tighter bounds on cosmological parameters (in an optimistic case, a factor of 20 in the figure of merit of dark energy can be achieved).
This article aims to prepare a cosmological parameter analysis with cosmic shear threepoint statistics by developing a pipeline to measure and model both the threepoint correlation functions Γ^{(i)} and thirdorder aperture mass statistics . We compare different estimators, computational costs and information content of both statistics. The covariance calculation will be investigated in the second publication of this series (Linke et al., in prep.). We further show that these statistics can be measured with relative ease in a StageIII survey, which constitutes a significant advantage over their Fourierspace counterpart, the convergence bispectrum. The aperture mass statistics have several additional advantages: They offer good data compression, are not subject to the masssheet degeneracy, and, most importantly, they decompose the signal into E and Bmodes, where to leading order only Emodes can be created by gravitational lensing.
In total, our modelling and validation pipeline can be summarised in the diagram in Fig. 1. Our modelling algorithm is publicly available^{1}. The paper is structured as follows: In Sect. 2 we introduce the Nbody simulations we use to validate and test our modelling pipeline. We then present the convergence bispectrum in Sect. 3, the shear threepoint correlation functions in Sect. 4 and the thirdorder aperture statistics in Sect. 5. For each of these statistics, we describe their theoretical background, how we chose to model them, how they are measured in simulations, and the validation tests we performed. We then compare their information content to one of the secondorder shear statistics in a mockMCMC in Sect. 6 and discuss our findings in Sect. 7.
Fig. 1. Schematic diagram of the modelling and validation pipeline introduced in this paper. The numbers on the arrows correspond to the respective section where this part of the pipeline is discussed. 
2. Model validation and covariance estimation with Nbody simulations
We use Nbody simulations containing only dark matter to validate our model and estimate covariance matrices. One of the main advantages of thirdorder shear statistics is that they can be easily adapted to different survey specifications. To highlight this, we use different simulation suites with varying source redshift distributions, galaxy number densities and cosmologies for the validation and parameter estimation. In this paper, we use the fullsky gravitational lensing simulations described in Takahashi et al. (2017, hereafter T17), the Millennium Simulations (Springel et al. 2005; Hilbert 2008, hereafter MS), and the ScinetLIghtCone Simulations (HarnoisDéraps et al. 2018, hereafter SLICS).
2.1. T17 simulations
The T17 are used in this work to perform a realistic analysis of a survey that mimics the KiDS1000 data and are constructed from a series of nested cubic boxes with side lengths of L, 2L, 3L... placed around a fixed vertex representing the observer’s position, with L = 450 Mpc/h. Each box is replicated eight times and placed around the observer using periodic boundary conditions. With the Nbody code GADGET2 (Springel et al. 2001) the gravitational evolution of 2048^{3} dark matter particles is simulated. Within each box, three spherical lens shells are constructed, each with a width of 150 Mpc/h, which are then used by the public code GRAYTRIX^{2} to trace the lightray trajectories from the observer to the last scattering surface^{3}. The cosmological parameters of the simulation are Ω_{m} = 1 − Ω_{Λ} = 0.279, Ω_{b} = 0.046, h = 0.7, σ_{8} = 0.82, and n_{s} = 0.97. The matter power spectrum agrees with theoretical predictions of the revised HALOFIT (Takahashi et al. 2012) within 5%(10%) for k < 5(6) h Mpc^{−1} at z < 1.
For each of the 108 realisations, we build a realistic convergence map by taking a weighted average of all 38 convergence shells at different redshifts, where the weights were determined by the fiducial KiDS1000 n(z) – see Fig. 2.
Fig. 2. Redshift distribution constructed from the T17 simulation given the fiducial n(z) of the KiDS1000 data. 
We then transform the pure convergence maps into realistic convergence maps by adding to each pixel a Gaussian random variable with a vanishing mean and standard deviation of
where A_{pix} is the pixel area of the convergence grid, and the effective number density n_{gal} = 6.17 arcmin^{−2} and σ_{ϵ} = 0.265 are chosen such that they are consistent with the combined 1–5 tomographic bins of the KiDS1000 data.
2.2. Millennium simulations
The MS were run with 2160^{3} particles in a 500 h^{−1} Mpc box in a flat ΛCDM cosmology with h = 0.73, σ_{8} = 0.9, Ω_{m} = 0.25, Ω_{b} = 0.045 and n_{s} = 1. Subsequently, shear and convergencemaps of 64 independent lines of sight with an area of 4 × 4 deg^{2} each were created at 36 different redshifts between z = 0 and z = 3 (Hilbert 2008; Hilbert et al. 2009). Each map is calculated on a grid of 4096 × 4096 pixel. For each line of sight, we use the full shear and convergence maps at redshift z = 1. As we use the MS solely to validate our model, we do not add any noise to the maps.
2.3. ScinetLIghtCone Simulations
The SLICS were run with 1536^{3} particles in a 505 h^{−1} Mpc box, filling up 10 × 10 deg^{2} lightcones up to z = 3. All SLICS were run in a flat ΛCDM cosmology with h = 0.69, σ_{8} = 0.83, Ω_{m} = 0.29, Ω_{b} = 0.047 and n_{s} = 0.969. From these simulations, we use convergence maps as well as galaxy catalogues with a redshift distribution of
with z_{0} = 0.637, β = 1.5 and the overall proportionality constant given by normalising the distribution to 30 gal/arcmin^{2}. We use mock galaxy catalogues provided for 924 (pseudo)independent lines of sight. The SLICS are useful to estimate the constraining power of thirdorder statistics since the threepoint correlation function can be calculated relatively quickly on the 100 deg^{2} fields, and the 924 lines of sight enable the determination of a stable covariance matrix.
3. Convergence power spectrum and bispectrum
In this section, we briefly recap the basics of the weak gravitational lensing formalism, focusing on the shear statistics in Fourier space. More detailed reviews can be found in Bartelmann & Schneider (2001), Hoekstra & Jain (2008), Munshi et al. (2008), Bartelmann (2010).
We start by defining the density contrast at comoving position x and redshift z, , where ρ(x, z) is the matter density at position x and redshift z and the average density at redshift z. From this density contrast, we define the convergence κ for sources at redshift z as a lineofsight integration, weighted by the lensing efficiency
where f_{K}(χ) is the comoving angular diameter distance. We note that throughout this paper, we work in a flat Universe with Ω_{m} + Ω_{Λ} = 1, meaning that f_{K}[χ(z)] = χ(z), where χ(z) is the comoving distance at redshift z. However, everything we present in this section also works for open or closed universes. The convergence can not be directly observed, but it can be recovered from an observed shear field (Kaiser et al. 1995; Seitz et al. 1998; Jeffrey et al. 2020). The relations between shear, convergence and the matter density contrast allow us to relate all second and thirdorder shear statistics to the wellunderstood matter power spectrum P_{δ}(k, z) and bispectrum B_{δ}(k_{1}, k_{2}, k_{3}, z).
3.1. Definition of power and bispectrum
The matter power spectrum P_{δ}(k, z) and bispectrum B_{δ}(k_{1}, k_{2}, k_{3}, z) can be defined as
where describes the Fourier transform of δ and δ_{D} is the Diracdelta distribution. The fact that the power and bispectrum only depend on the moduli of the kvectors can be easily derived from the statistical isotropy of the Universe.
The convergence power and bispectrum can then be computed using the Limber approximation (Limber 1954; Peebles 1980; Kaiser & Jaffe 1997; Bernardeau et al. 1997; Schneider et al. 1998),
Here,
describes the lensing efficiency, where p(χ) is the (normalised) comoving distance probability distribution of sources. We note that usually one instead measures a redshift probability distribution n(z); in our modelling pipeline we thus instead write Eqs. (6)–(8) as integrals over the redshift z.
The Limber approximation breaks down for small values of ℓ. For example, the bispectrum from the Limber approximation overestimates the truth by up to an order of magnitude for ℓ ≪ 60 (corresponding to angular scales of roughly 6°), depending on the source redshift distribution (Deshpande & Kitching 2020). However, at these scales, the impact of nonlinear structure formation is small, meaning that the matter distribution is welldescribed by a Gaussian and higherorder statistics like the bispectrum are small. In this work, we only consider shear statistics up to ∼4°; at these scales, we expect the Limber approximation to hold.
Instead of the three ℓvalues ℓ_{1}, ℓ_{2} and ℓ_{3}, we can also describe the bispectrum as a function of the two vectors ℓ_{1}, ℓ_{2} or their moduli ℓ_{1}, ℓ_{2} and the angle φ between them. We define
3.2. Modelling the bispectrum
We use the stateoftheart B IHALOFIT algorithm (Takahashi et al. 2020) to model the dark matter bispectrum on nonlinear scales. In comparison to older bispectrum models (e.g., GilMarín et al. 2012; Scoccimarro & Couchman 2001), B IHALOFIT appears to trace the nonequilateral triangles much better: In comparison with Nbody simulations, B IHALOFIT retains an accuracy of 10% or better, whereas the other two fitting formulae are subject to errors of more than 200% for squeezed triangle configurations (compare Figs. 5, 12, and 13 of Takahashi et al. 2020). The effects on higherorder shear statistics are substantial, as can be seen in Halder et al. (2021), who modelled a different thirdorder shear statistic from using both B IHALOFIT and the bispectrum model of GilMarín et al. (2012). In a direct comparison with Nbody simulations, B IHALOFIT is accurate on all tested scales, whereas the other fitting formula breaks down on scales of ≲30′.
Another advantage of B IHALOFIT is that it only requires a linear power spectrum as its input, compared to a nonlinear spectrum in GilMarín et al. (2012) or Scoccimarro & Couchman (2001). We use the fitting formula developed by Eisenstein & Hu (1999) to model the linear power spectrum.
One of the main advantages of thirdorder shear statistics is that one can rigorously test each stage of the modelling pipeline. We perform such a test on our bispectrum model in Appendix A and conclude that the Limber integrated B IHALOFIT bispectrum is consistent with the MS up to ℓ ≲ 10^{4} and deviates by up to 40% for larger values ℓ.
4. Shear threepoint correlation functions
4.1. Definition of the shear threepoint correlation functions
Shear threepoint correlation functions (3pcf) are the natural extension to the widely used twopoint correlation functions. Let γ_{c} = γ_{1} + iγ_{2} denote the complex shear in Cartesian coordinates. Considering a triangle of galaxies, as a first step, we project the shear γ^{i} of each galaxy i to its tangential and crosscomponents with respect to a point fixed with respect to the triangle, for example one of its centres,
where ζ is the angle of the projection direction. Schneider & Lombardi (2003) established four natural components of the shear 3pcf, which remain invariant under rotations of the triangle. They are defined as
where the ‘^{*}’ denotes complex conjugation. The choice of the reference point of the projection is to some degree arbitrary. Usually one of the triangles cenres is chosen, most often the orthocenter (the intersection of its three altitudes) or the centroid (the intersection of its three medians), as shown in Fig. 3. The natural components have the advantage that they are invariant under the choice of triangle centre up to multiplication with a complex phase factor, meaning that their moduli are invariant under the choice of triangle centre.
Fig. 3. Example triangle with our used notations. The sidelengths are denoted by x_{1} = X_{3} − X_{2} and cyclic permutations. The interior angles of the triangle are denoted by ϕ_{i}, whereas φ_{i} are the angles that the sidelengths take with respect to the xaxis. We also show two possible definitions of its centre: The orthocenter O is at the intersection of the three altitudes (blue dashed); the centroid C is at the intersection of its three medians (red dashed). Figure adapted from Schneider et al. (2005). 
Parametrising the shear 3pcf by the triangle sidelengths, x_{1}, x_{2} and x_{3}, where the indices 1, 2 and 3 are ordered in a counterclockwise direction, the natural components exhibit a nice behaviour concerning cyclic permutation of arguments. While the first natural component Γ^{(0)} is invariant under cyclic permutations, the other three components transform into each other,
Similar behaviour can be observed for parity transformations:
4.2. Modelling the shear threepoint correlation functions
We model the natural components of the threepoint correlation functions of cosmic shear using the methods described in Schneider et al. (2005, hereafter S+05): When we project the shear of all three galaxies to the orthocenter of the triangle, the projection direction at the triangle vertex X_{i} is orthogonal to the orientation φ_{i} of the triangle side x_{i}. Thus, the shear transforms as
We now utilise the relation between convergence and shear in Fourier space (Kaiser & Squires 1993),
where β_{ℓ} is the polar angle of ℓ, to write
This can then be transformed into equation (15) of S+05^{4}:
with
The quantities A_{1, 2} and α_{1, 2} are obtained by cyclic permutation of indices. The angles ϕ_{i} are the interior angles of the triangle, as shown in Fig. 3. By introducing polar coordinates , ψ = arctan(ℓ_{2}/ℓ_{1}), we get
with
Defining
and E_{1} and E_{2} via cyclic permutations of indices, we can write
The Rintegration filters the bispectrum with a 6th order Bessel function, making the integration routine difficult to solve numerically, as the functional form of the bispectrum prevents the application of fast Hankel transform algorithms like FFTLog (Hamilton 2000)^{5}. We thus use the method developed in Ogata (2005) to solve the Rintegration and integrate the remaining dimensions using the CUBATURE library^{6}. To model Γ^{(1)}, we apply the same transformations to Eq. (18) of S+05; for Γ^{(2)} and Γ^{(3)} we perform a cyclic permutation of the input variables as outlined in Eq. (12).
While a triangle of galaxy positions for which we evaluate the threepoint correlation function can be described by its three sidelengths x_{1}, x_{2} and x_{3}, it is certainly not a good idea to use these variables for a binning scheme; for example, for x_{1} > x_{2} + x_{3} a triangle can not be defined, which means that the 3pcf would not be defined for many bins. A better way to bin the triangles was introduced by Jarvis et al. (2004). Assuming x_{1} > x_{2} > x_{3} they defined a triangle via the values r ∈ [0, ∞], u ∈ [0, 1] and v ∈ [ − 1, 1] by
Here, v is positive for triangles where x_{1}, x_{2} and x_{3} are oriented clockwise and negative for a counterclockwise orientation. This binning choice allows us to bin the triangle size r in logarithmic steps without having bins where the 3pcf is not defined. In all cases, we bin the shear 3pcf logarithmic in r and linear in u and v.
We note that Eq. (12) implies that the four shear 3pcf for x_{1} > x_{2} > x_{3} (as in the binning scheme of Jarvis et al. 2004) already contain the entire information content of the thirdorder shear signal, as does knowledge of Γ^{(0)} and Γ^{(1)} for all combinations of x_{1}, x_{2} and x_{3}. In a similar manner, Eq. (12) implies that Γ^{(i)}(r, u, v) = Γ^{(i)*}(r, u, −v) holds, where the ‘^{*}’ denotes complex conjugation. To ensure compatibility with results from the measured 3pcf, we transform the modelled functions from the centroid (as used in S+05) to the orthocenter (as used in Jarvis et al. 2004, compare Sect. 4.3).
For a potential cosmological parameter analysis, the threepoint correlation functions face a few hurdles: Assuming we bin the threepoint correlation functions in 10 bins for each r, u and v, then our data vector consists of 8000 entries. This makes estimating a covariance matrix using simulations practically impossible and leads to a modelling time that is unfeasible even for a nontomographic analysis.
4.3. Measuring the shear threepoint correlation functions
We use the public treecode TREECORR (Jarvis et al. 2004) to measure the threepoint correlation functions Γ_{i}. This algorithm estimates the quantity
where the ε_{i} = ε_{t, i} + iε_{×, i} are the observed ellipticities of galaxies and w_{i} the associated weights. The other natural components are estimated in the same manner. This estimator has the advantage that it is not impacted by the survey geometry: As long as at least one galaxy triplet falls into each bin, it remains unbiased (Simon et al. 2008)^{7}. The disadvantage of the estimator is that its computational complexity scales with , which is not feasible to execute even for moderate values of N_{gal} ≳ 10^{6}. That is why TREECORR constructs a hierarchical ball tree out of the galaxy sample and calculates the correlation functions from this tree. This results in a remarkable speedup and allows us to calculate the shear 3pcf for an ensemble of about 10^{7} source galaxies distributed over a 10 × 10 deg^{2} field in about 1500 CPUh. A disadvantage of the treecode is that its execution time scales massively with the number of bins: If b is the logarithmic bin size, then the runtime scales roughly with b^{−4}.
The TREECORR algorithm also has a BINSLOP parameter, which allows balls of the KDtree to overlap the edges of a bin. This parameter heavily affects computation time, and while the expectation value is relatively stable between different values of BINSLOP, the covariance is subject to change (Secco et al. 2022).
4.4. Validation
We test our developed integration routine described in Sect. 4.2 with a lensing potential for which we can derive analytic expressions for both the convergence bispectrum and the shear 3pcf. As discussed in more detail in Appendix B we found that the analytic 3pcf and the one obtained by integration from the bispectrum agree to the subpercent level.
To validate our model for the shear 3pcf, we measure the shear signal at a redshift of z = 1 in the MS. We choose to measure the signal in 10^{3} bins (10 bins in each r, u and v), with logarithmic rbins from 0′.1 to 120′. To speed up computation time, we randomly select every tenth pixel of the 4096^{2} pixel grid. As we do not include shape noise, we expect the loss of signal to be small. The results can be seen in Fig. 4. We conclude that we can model the shear 3pcf reliably down to subarcminute scales for almost all triangle configurations. Only for almost degenerate, flattened triangle configurations (v > 0.9) do we see that the model and simulations differ significantly. This might signify that B IHALOFIT breaks down at the corresponding triangle configurations in Fourierspace. Alternatively, this might point towards a breakdown of the treecode’s accuracy at these very degenerate triangles. As we can see in Sect. 5.4.2, these points play a negligible role in the conversion to aperture mass statistics. Overall we see that the agreement between shear 3pcf is better than the one at the bispectrum level (compare Appendix A.2). We observe the same effects for the other natural components of the 3pcf (compare Fig. C.1).
Fig. 4. First natural component Γ^{(0)} of the shear threepoint correlation functions modelled from B IHALOFIT compared to measurements from the MS that were extracted using TREECORR. Each panel corresponds to one fixed triangle configuration; the u and v values are listed at the top, and the corresponding shape of the triangle is shown on the bottom right. The xaxis corresponds to the length of the red side of the triangle. We show both the real part (red) and the imaginary part (blue) for the model (dashed line) and the simulations (crosses). Under each panel we show the ratio between measurement and model. The error bars denote the error on the mean of the 64 lines of sight of the MS. 
5. Aperture mass statistics
5.1. Definition of aperture mass statistics
An alternative way to analyse cosmic shear is via aperture mass maps (Schneider 1996; Bartelmann & Schneider 2001). Their advantage is that they can separate the signal into socalled E and Bmodes (Schneider et al. 2002), where Bmodes can, to leading order, not be created by the weak gravitational lensing effect. Thus, the absence of Bmodes provides a crucial nulltest in a cosmic shear analysis (see, e.g., Hildebrandt et al. 2017; Asgari et al. 2019). Compared to convergence maps (Kaiser & Squires 1993; Seitz & Schneider 2001; Gatti et al. 2022), aperture mass maps are constructed in a way that they do not suffer from the wellknown masssheet degeneracy (Falco et al. 1985; Schneider & Seitz 1995).
The aperture mass ℳ_{ap} at position ϑ and filter radius θ are defined as
here, U_{θ}(ϑ) is a compensated filter (i.e. ∫dϑ ϑ U(ϑ) = 0). Given a shear field γ, the aperture mass ℳ_{ap} and its respective Bmode counterpart ℳ_{⊥} can be calculated as
where γ_{t} and γ_{×} are projected along the vector ϑ′−ϑ (see Eq. (10)) and Q_{θ} is related to U_{θ} via
For simplicity of notation, we define U_{θ}(ϑ) = θ^{−2}u(ϑ/θ) and denote by the Fourier transform of u. In this work, we opt for the filter function introduced by Crittenden et al. (2002),
While the construction of aperture mass maps has its uses (see for example HarnoisDéraps et al. 2021; Heydenreich et al. 2021), we do not care about the structure of an aperture mass map but rather about its statistical properties. We define for arbitrary combinations of E and Bmode aperture mass statistics
By construction, ⟨ℳ_{ap}⟩(θ) vanishes. In a paritysymmetric field, all odd powers of Bmode components vanish (Schneider 2003), meaning that the relevant Bmode counterparts to and are and , respectively.
5.2. Modelling aperture mass statistics
Given a model for the convergence power spectrum, the secondorder aperture statistics can be calculated as
As a model for the nonlinear power spectrum, we use the revised HALOFIT model of Takahashi et al. (2012). Equivalently, the thirdorder aperture statistics can be derived from a bispectrum model via (compare Jarvis et al. 2004; Schneider et al. 2005)^{8}
We use the public CUBATURE library to solve this integration. In our implementation, the integration kernel is executed on a graphics processing unit (GPU), yielding a significant speedup over parallelisation on central processing units (CPUs).
5.3. Measuring aperture mass statistics
There are three main methods to estimate the thirdorder aperture mass statistics , first, via the convergence field κ, second, via the shear field or, in practice, from the observed galaxy ellipticities, and third, via the thirdorder correlation functions Γ^{(i)}.
5.3.1. Measuring aperture mass statistics from convergence maps
The most straightforward way is to measure aperture mass maps directly on a convergence field using Eq. (25). In a real survey, this is difficult, as the convergence is not directly observable. In principle, one could compute the aperture mass statistics of a reconstructed convergence field, but this is not a good way to estimate aperture statistics, as the convergence reconstruction yields a convergence map that is necessarily smoothed and potentially also contains other systematic effects caused by masks or boundaries of the survey (Pires et al. 2020). While not really applicable to real data, this method yields a quick and unbiased way to estimate aperture statistics in lightcones from simulations, as for them, convergence maps are readily available. However, one faces the issue of boundary effects when the integral in Eq. (25) extends past the simulation boundary. To avoid this issue, we cut off a slice of width 4θ from the computed aperture mass maps^{9}.
To extract the data vectors from the (fullsky) convergence maps of the T17 simulations, we smooth the maps with HEALPY function SMOOTHING, with a given beam window function created by the function BEAM2BL, which in turn is determined by the corresponding U_{θ} filter. For each filter radius θ, this yields a fullsky aperture mass map , without the need to cut off boundaries.
5.3.2. Measuring aperture mass statistics from galaxy catalogues
Another way to estimate aperture statistics is from an ensemble of observed galaxy ellipticities, using
where ε_{t/×} are the observed galaxy ellipticities converted into their tangential and croos components according to Eq. (10); ϑ_{i} are their respective positions. Here, n_{gal} can be the global number density of galaxies (Bartelmann & Schneider 2001) or the number density of galaxies within the aperture radius (Martinet et al. 2018). For this work, we define n_{gal} as the number of galaxies weighted by the Qfilter function:
We tested all three definitions of n_{gal} using the SLICS and found that, for randomly distributed galaxies, setting n_{gal} as the number density within the aperture radius or the one weighted by the Qfilter function induces subpercent differences on the thirdorder aperture masses . However, setting n_{gal} as the global galaxy density can induce differences of about 5% in .
In the following, we adopt n_{gal} to be the number of galaxies weighted by the Qfilter function. Rewriting Eq. (10) as
where the vector ϑ − ϑ′ denotes the projection direction of the tangential shear in complex notation, we can rewrite Eqs. (26) and (32) as a convolution. To calculate aperture mass fields, we distribute galaxies on a grid using a cloudincell method. From this, we compute the aperture masses using a Fast Fourier Transform (FFT), allowing us to compute an aperture mass map in 𝒪(N_{pix} log N_{pix}) operations. To compute second and thirdorder aperture statistics, we apply the same cutoff of 4θ_{ap} to the aperture mass maps. From these aperture mass maps, we obtain estimates for the second and thirdorder aperture statistics by multiplication of the respective aperture mass maps on each pixel, and then taking the average of all pixel values.
5.3.3. Measuring aperture statistics from threepoint correlation functions
While the abovementioned method to estimate aperture statistics is extremely fast, it can not be applied to realistic survey data. Assuming a relatively large aperture radius of θ_{ap} = 30′, we would have to cut off a 2°strip around every edge or mask in the survey footprint, meaning that we would disregard most of the data. While active research is being conducted to circumvent these problems (Porth et al. 2020; Porth & Smith 2021), the arguably best method to estimate thirdorder aperture statistics from real data is to derive them from the measured 3pcf, as has been introduced in Jarvis et al. (2004), generalised in Schneider et al. (2005) and applied to survey data in Fu et al. (2014) and Secco et al. (2022). The shear 3pcf can be estimated straightforwardly from a survey with arbitrarily complex geometry, meaning that the converted aperture statistics are not biased by boundary effects. One caveat is that this conversion requires the knowledge of the 3pcf for all triangle configurations, particularly for infinitesimally small or extremely large ones, both of which can not be measured. The incomplete knowledge of the correlation functions can lead to a mixing of E and Bmodes for the aperture statistics (Kilbinger et al. 2006). However, for thirdorder aperture statistics, this effect appears to be not as severe as for their secondorder counterpart (at least for the diagonal part of the aperture statistics, this has been demonstrated in Shi et al. 2014, we are testing this assumption for nondiagonal aperture mass statistics in Sect. 5.4.1). Despite its advantages, this method comes at the price of computation time: Calculating the shear 3pcf for a realistic number of source galaxies takes orders of magnitude longer^{10} than the direct estimation of aperture statistics, so calculating the shear 3pcf of an ensemble of simulations (as would be necessary to estimate a covariance matrix) comes at a prohibitively high computational cost.
The method to transform shear threepoint correlation functions into aperture statistics is already implemented in TREECORR. To compute the thirdorder aperture statistics, the quantities
need to be computed. Here, A_{1, 2} and F_{1, 2} are the prefactors and filter functions, which are specified in S+05 (see Eqs. (62) and (71)). The aperture statistics and are linear combinations of these quantities:
We do not consider the quantities and , as they vanish for any paritysymmetric field (Schneider 2003). In TREECORR, this method is implemented in the following way: First the Γ^{(1)} and Γ^{(2)} are transformed into Γ^{(3)} via Eq. (12). For each bin (r, u, v), the transformation matrix is then calculated. The value of the respective integral is computed as the sum of the values of Γ^{(3)} multiplied by the determinant of the transformation matrix and the value of the filter functions F_{1, 2} at the bin centre. While numerically very cheap, this is probably not the most efficient way to compute that integral. In case the filter functions vary significantly over a bin (compare Figs. 2 and 3 of Schneider et al. 2005), it might be more appropriate to calculate the average of the filter function in a bin, for example. We tried to improve the integration results by interpolating the measured shear 3pcf and performing the same integral, achieving rather moderate improvements. We leave an optimisation of the conversion from shear 3pcf to aperture mass statistics for future work.
5.4. Validation
5.4.1. Binning choice of threepoint correlation functions
As a first step, we want to validate the conversion Γ^{(i)}⇢ performed by the TREECORR algorithm. In particular, we investigate the number of bins necessary to achieve an unbiased estimate for thirdorder aperture statistics and quantify the leakage of E and Bmodes. While the latter has already been investigated by Shi et al. (2014), we extend upon these results by using a realistic convergence bispectrum model and by taking into account nondiagonal aperture mass statistics.
Our modelling pipeline provides us with the ability to test this conversion. As we start from the same convergence bispectrum B_{κ}, the aperture mass statistics achieved by the conversion and by direct modelling have to be consistent. Furthermore, the modelled Γ^{(i)} are pure Emode functions, so any Bmodes that we observe have to be created by the transformation . These tests would be unfeasible to perform with simulations due to the prohibitively high computational cost of extracting the shear 3pcf from an extensive simulation set for different bin sizes.
The results of our tests can be seen in Fig. 5. We see that as long as the three filter radii θ_{1}, θ_{2} and θ_{3} are similar, the conversion appears to work reasonably well. Only when we bin Γ^{(i)} in 7^{3} bins do we get significant deviations, meaning that this is certainly not a sufficient number of bins. The conversion becomes less accurate when two filter radii are much smaller than the third one, as shown in the topright corner of Fig. 5. There also the results for 10^{3} bins show significant deviations, whereas the results for 15^{3} bins seem overall consistent with the ones from 20^{3} bins. We also observe a nonnegligible amount of Bmode leakage for these cases, even for 15^{3} and 20^{3} bins. We are planning to exclude all combinations of filter radii from cosmological parameter analyses where we observe a Bmode leakage of more than 10% for the 3pcf in 15^{3} bins.
Fig. 5. Thirdorder aperture mass statistics when modelled directly from the bispectrum (Eq. (31), blue line) compared to the ones constructed from modelled shear 3pcf (Eq. (35), coloured crosses). The dotted lines denote the respective Bmodes. We plot (θ_{1}, θ_{2}, θ_{3}), where θ_{1} is constant in each row, θ_{2} is constant in each column and θ_{3} varies along the xaxis in each panel. The yellow, green, red and purple dots denote the we get from the shear 3pcf when we bin r, u and v in 7 × 7 × 7, 10 × 10 × 10, 15 × 15 × 15 and 20 × 20 × 20 bins, respectively. 
Inspecting Fig. 6 we further note that the function Γ^{(1)}(r, u, v) is relatively smooth and wellbehaved with respect to r, but strongly varies as a function of u and v, especially when u ≈ 0 and v ≈ ±1. This implies that r can be binned rather coarsely, as long as u and v are finely binned. This is in contrast to binning choices in other studies, for example Secco et al. (2022), who preferred a fine binning in r (55 bins) and a coarser binning in u and v (10 bins).
Fig. 6. Threepoint correlation function Γ^{(1)} as a function of its three parameters r (left), u (middle) and v (right). We note that the jump in the imaginary part of Γ^{(1)} around v = 0 is due to the fact that Γ^{(i)}(r, u, v) = Γ^{(i)*}(r, u, −v) holds. 
5.4.2. Comparison to Nbody simulations
We compare the modelled aperture mass statistics with the ones we measure in the MS in Fig. 7. As expected from our discussions in Sect. 5.4.1, we observe that the conversion from shear 3pcf to thirdorder aperture masses fails when two aperture radii are small, and the third one is large (as observed in the topright panel). In these cases, we also register significant Bmodes. We also note that in most cases, the uncertainties on the direct measurements are significantly larger than those from the shear 3pcf. Cutting off a stripe around the boundary in order to avoid edge effect in the estimate of (as discussed in Sect. 5.3), leads to a loss of information compared to measuring the 3pcfs, for which all triplets of available points^{11} in the field are used.
Fig. 7. Thirdorder aperture masses from our modelling pipeline using B IHALOFIT (Eq. (31), blue) compared to the ones measured directly from the MS. The direct measurements on the convergence maps (Eq. (25)) are shown in green, and the measurements converted from shear 3pcf (Eq. (35)) are shown in orange. Both statistics were computed from 32 lines of sight of the MS, the error bars denote the error on the mean. We note that for the largest angle of 32′, a direct measurement could not be obtained due to the limited size of the individual light cones. 
6. Cosmological parameter estimation
To perform a cosmological parameter inference, we use our described pipeline to create a model vector. For the covariance, we rely on Nbody simulations, where we use the method of Percival et al. (2022) to debias the estimated covariance matrix . Given a data vector d and a covariance matrix measured from n_{r} simulated survey realisations, the posterior distribution of a model vector m(Θ) depending on n_{Θ} parameters, is
where
The powerlaw index m is
with n_{d} being the number of data points and
If m = n_{r} the formalism of Sellentin & Heavens (2016) is recovered.
Normally, one needs to evaluate this likelihood function in highdimensional parameter space to perform a cosmological parameter analysis with thirdorder statistics. The creation of the model vector for thirdorder aperture statistics with 35 combinations of aperture radii takes about one minute on an NVIDIA A40 GPU, and needs to be evaluated at about 10^{4} points even for sampling methods like POLYCHORD (Handley et al. 2015). This means that a complete likelihood analysis with thirdorder aperture statistics is possible but takes a long time, whereas a complete likelihood analysis with the 3pcf, where the modelling of the 3pcf in 10^{3} bins takes two to three hours, is not feasible. To solve this issue, we use a neural network emulator called COSMOPOWER (Spurio Mancini et al. 2022), which was first developed to emulate power spectra but can easily be adapted for arbitrary vectors. Since a neural network emulator needs as many as possible evaluation points, we calculate our model at 7500 points^{12} in a fourdimensional Latin hypercube describing a flat wCDM cosmological model, varying the parameters Ω_{m}, S_{8}, w_{0} and h with flat priors of Ω_{m} ∈ [0.1, 0.55], S_{8} ∈ [0.6, 0.9],w_{0} ∈ [ − 2, −0.5] and h ∈ [0.6, 0.82]. We leave all other parameters fixed at the values corresponding to the SLICS, which we introduced in Sect. 2. We then train the emulator with 6500 points and use the remaining 1000 as a validation test, as shown in Fig. C.2. This neural networkbased emulator performs extraordinarily well, modelling all statistics with a subpercent accuracy. As our bispectrum model is only accurate to about 10% (Takahashi et al. 2020), the emulator uncertainty plays a negligible role in the modelling process. This emulator yields extremely fast predictions for the thirdorder statistics, so that we can perform a full likelihood analysis using the MetropolisHastings sampling method.
Even though we trained our emulator on the cosmological parameters Ω_{m}, S_{8}, w_{0} and h, we notice that in the absence of tomography we can not constrain w_{0} or h. Therefore we fix these parameters in all subsequent MCMC analyses to their fiducial values w_{0} = −1 and h = 0.7. For Ω_{m} and S_{8} we use the same priors that we took for the creation of the training sample.
6.1. Shear threepoint correlation functions vs. thirdorder aperture masses
One aspect of this work is investigating the information loss using aperture mass statistics instead of the 3pcf itself. Although the aperture mass statistics are well suited for a cosmological parameter inference due to their E/Bmode decomposition and fast modelling times, they should not be used if the loss of information is too severe.
Unfortunately, we cannot quantify the full information content of the shear 3pcf, as the data vector contains about 10^{4} entries, and therefore a reliable covariance matrix is not accessible. To circumvent this problem, we model the shear 3pcf in 10^{3} bins in r, u and v, where we bin r logarithmically from 0.′1 to 100′, at 500 of the 7500 training nodes. We then perform a principal component analysis (PCA) to decide on the 20 most relevant principal components of the 3pcf data vector^{13}.
Using this PCA, we determine the covariance matrix for the shear 3pcf, which we measure from 200 10 × 10 deg^{2} tiles of the SLICS in the same configuration as the training data. For the model vector, we again use COSMOPOWER, which is trained on the PCA components of the 500 models used in determining the principal components.
In Fig. 8 we compare the constraining power of the PCA analysis of 3pcf to the analysis, where the covariance for PCA analysis of 3pcf is measured from 200 SLICS realisations; for we use all 927 available SLICS realisations and take into consideration all combinations of the filter radii 0′.5, 1′,2′,4′,8′,16′, and 32′. The different number of realisations are considered by Eq. (39). In both cases, the data vector was created by the COSMOPOWER Emulator. It is clearly seen that the constraining power from the PCA analysis of 3pcf is only slightly better than the one from , and this slight difference may very well be explained by the use of different scales between the 3pcf and aperture mass statistics, and the fact that the aperture mass maps are smaller due to the cutoff at the boundaries. Overall, the advantages of the justify their use, even considering their potentially slightly lower constraining power.
Fig. 8. Comparison of the posterior distributions of different thirdorder statistics using the SLICS simulations to determine the covariance. In red we show the posteriors of the PCA from the 3pcf and in blue the one resulting from the statistic. 
6.2. Combination of second and thirdorder aperture masses
To assess the constraining power for thirdorder aperture statistics, especially when combined with secondorder shear statistics, we perform a mock analysis for a nontomographic KiDS1000like setup. To estimate the covariance matrix, we made use of all 108 realisations of the T17 simulations with a resolution NSIDE = 4096 (corresponding to a pixel size of 0′.74^{2}). From each realisation, we extract 18 HEALPix squares of size ≈860 deg^{2} that do not share common borders. This results in 1944 independent realisations from which the covariance matrix is estimated, such that with a data vector size of ∼50, the statistical noise of the covariance matrix can be neglected. Since the area of the square patches is slightly larger than the one from KiDS1000 with ≈777.4 deg^{2}, the covariance needs to be rescaled by a factor of 1.11.
Furthermore, in order to have a data vector as unbiased and noisefree as possible, we estimate it with one fullsky realisation with a resolution NSIDE = 8192 (corresponding to a pixel size of 0.18 arcmin). Lastly, we use the filter scales of (4, 8, 16, 32) arcmin, as the model and simulations are inconsistent for smaller scales.
The resulting posterior distribution is shown in the left panel of Fig. 9, where we first notice that the combination of second and thirdorder statistics significantly increases the constraining power, especially in the Ω_{m}σ_{8} panel to the different degeneracy directions of the individual statistics (compare Table 1). Indeed, a joint analysis increases the constraints on S_{8} by 42% with respect to secondorder statistics; the constraints on Ω_{m} and σ_{8} increase by at least 68% and 54%, respectively^{14}. The figure of merit (Albrecht et al. 2006) in the Ω_{m}σ_{8} plane increases by a factor of 5.88. Finally, we note that the true cosmology is well within 1 σ of the expected KiDS1000 uncertainty for all three statistics.
Fig. 9. Left panel: Comparison of the posterior distributions of second and thirdorder aperture statistics with a joint analysis within a KiDS1000like setup. The data vector and covariance matrix are estimated by the T17 simulations, and we use filter scales of (4, 8, 16, 32) arcmin. The Hubble parameter h and the dark energy equationofstate are fixed to the T17 values. Right panel: Comparison of the combined second and thirdorder aperture statistics with the combination of secondorder aperture statistics with the equalscale thirdorder aperture statistics. The corresponding data vector and model vector can be found in Fig. C.3. 
Additionally, we investigate the constraining power if only equalscale aperture masses (θ, θ, θ) are used. As displayed in the right panel of Fig. 9, the loss of constraining power by the limitation to equalscale aperture masses is small, although not zero. This is in stark contrast to Kilbinger & Schneider (2005), who found a strong difference in constraining power using a Fisher forecast. However, their analysis was conducted using a covariance matrix from a significantly smaller set of simulations. Our results are roughly in line with the findings of Fu et al. (2014), who found rather marginal differences in an MCMC. We note that our analysis is quite simplified: We only vary two cosmological parameters and do not include any systematic effects. The introduction of additional parameters or tomography might lead to some degeneracies in equalscale aperture mass statistics that could be broken by combinations of different filter radii. Furthermore, future analyses might be able to use a larger dynamic range of scales. In this case the nonequal filter radii contain information about more flattened and squeezed bispectrum configurations, which can yield additional constraining power.
7. Discussion
In this work, which is the first of a series on cosmological analysis with thirdorder shear statistics, we have shown that a cosmological parameter analysis with thirdorder aperture mass statistics is feasible and beneficial for StageIII surveys. Both the shear 3pcf and the thirdorder aperture statistics can be modelled from the matter bispectrum, and we found that our models based on the B IHALOFIT bispectrum model are accurate enough for StageIII surveys. In particular, the flatsky and Limber approximations are valid for our selected range of scales, so the accuracy of our model is mainly limited by the accuracy of the bispectrum model. We note that we have not yet tested the impact of astrophysical or observational systematics.
We developed a test for binning strategies of the threepoint correlation functions in order to obtain unbiased estimates for aperture mass statistics and found that, at our selected scales, a measurement of the threepoint correlation functions in 15^{3} bins yields good results. In particular, we found that the leakage between E and Bmodes is at the percent level, and the bias in the aperture mass statistics is well below the sample variance for a StageIII survey. This extends upon the findings of Shi et al. (2014), who found a percentlevel leakage for diagonal aperture mass statistics utilising a simplified bispectrum model. We emphasise that, in addition to the effect of a minimum scale in the shear 3pcf investigated by Shi et al. (2014), our approach also quantifies the leakage that stems from the binning choices in the shear 3pcf, resulting in an inaccurate evaluation of the integral in Eq. (35).
We have tested the information loss when converting the shear threepoint correlation functions to thirdorder aperture mass statistics by performing a principal component analysis of the threepoint correlation functions and comparing the constraining power of the principal components to the one of thirdorder aperture statistics. We found comparable information content, suggesting that thirdorder aperture statistics constitute a good data compression method for the shear threepoint correlation functions. In addition to an easier modelling, the aperture statistics have the added advantage that they cleanly separate E and Bmodes.
We demonstrate that the computational load of a cosmological parameter analysis with thirdorder aperture statistics is manageable, particularly when utilising an emulator to speed up the MCMC. We make a COSMOSISmodule of our modelling algorithm publicly available^{15}.
Finally, we compare the constraining power between secondorder, thirdorder, and joint aperture statistics analysis. While secondorder aperture statistics are not being used in modern cosmological parameter analyses, we assume that all secondorder shear statistics exhibit similar constraining power and parameter degeneracies (compare Asgari et al. 2021). We find that thirdorder aperture statistics alone have a lower constraining power than their secondorder counterpart, but they exhibit a different degeneracy direction in the Ω_{m}–σ_{8}plane so that a joint analysis almost doubles the constraining power on the structure growth parameter S_{8} and increases the figure of merit in the Ω_{m}–σ_{8}plane by a factor of 5.9. However, the information gain predicted by Fisher analyses, especially for the difference between diagonal and full thirdorder aperture statistics (compare Kilbinger & Schneider 2005), appears overly optimistic. This suggests that a Fisher forecast might not be an optimal tool to forecast parameter constraints and mock sampling methods give more realistic constraints.
While we have demonstrated here that cosmological analyses with thirdorder shear statistics are feasible and promising, there are steps left to do before applying our methods to a concrete cosmological survey. The first of these is the development of a model for the covariance of thirdorder statistics, which is essential for a tomographic analysis and will be addressed in the following paper of this series (Linke et al., in prep). Additionally, as mentioned above, our model does not yet incorporate systematic and astrophysical effects, like baryonic feedback or intrinsic alignments of source galaxies (Semboloni et al. 2013; Halder & Barreira 2022; Pyne et al. 2022). We will develop and test strategies for treating these effects in future works of this series.
These maps are freely available for download at http://cosmo.phys.hirosakiu.ac.jp/takahasi/allsky_raytracing/
We note that due to the different definitions of the convergence bispectrum (Eq. (5) vs. Eq. (4) in S+05), we get a factor of 3 difference.
Even if certain bins remain empty, the estimated correlation function can be rebinned with a tesselation scheme to yield unbiased values of Γ^{(i)} for all bins, as was shown by Linke et al. (2020) for the related galaxygalaxygalaxylensing correlation function.
Acknowledgments
We thank Benjamin Joachimi and Mike Jarvis for providing valuable insights to this project. We would like to thank Joachim HarnoisDéraps for making public the SLICS mock data, which can be found at http://slics.roe.ac.uk/. This work was funded by the TRA Matter (University of Bonn) as part of the Excellence Strategy of the federal and state governments. This work has been supported by the Deutsche Forschungsgemeinschaft through the project SCHN 342/151. SH acknowledges support from the German Research Foundation (DFG SCHN 342/13), the International MaxPlanck Research School (IMPRS) and the German Academic Scholarship Foundation. Author contributions. All authors contributed to the development and writing of this paper. SH wrote the pipeline to model the bispectrum and shear 3pcf and the methods to measure the thirdorder statistics. LL implemented the modelling algorithm for aperture mass statistics, the GPUintegration, and the cosmosis module, aside from making various improvements to the codes. PB was responsible for everything regarding the T17 simulations and the MCMC runs, including the COSMOPOWER emulator. PS gave countless valuable insights into thirdorder shear statistics.
References
 Abbott, T. M. C., Aguena, M., Alarcon, A., et al. 2022, Phys. Rev. D, 105, 023520 [CrossRef] [Google Scholar]
 Albrecht, A., Bernstein, G., Cahn, R., et al. 2006, ArXiv eprints [arXiv:astroph/0609591] [Google Scholar]
 Asgari, M., Heymans, C., Hildebrandt, H., et al. 2019, A&A, 624, A134 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Asgari, M., Tröster, T., Heymans, C., et al. 2020, A&A, 634, A127 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Asgari, M., Lin, C.A., Joachimi, B., et al. 2021, A&A, 645, A104 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bartelmann, M. 2010, CQG, 27, 233001 [NASA ADS] [CrossRef] [Google Scholar]
 Bartelmann, M., & Schneider, P. 2001, Phys. Rep., 340, 291 [Google Scholar]
 Bernardeau, F., van Waerbeke, L., & Mellier, Y. 1997, A&A, 322, 1 [NASA ADS] [Google Scholar]
 Burger, P., Friedrich, O., HarnoisDéraps, J., & Schneider, P. 2022, A&A, 661, A137 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Crittenden, R. G., Natarajan, P., Pen, U.L., & Theuns, T. 2002, ApJ, 568, 20 [NASA ADS] [CrossRef] [Google Scholar]
 Deshpande, A. C., & Kitching, T. D. 2020, Phys. Rev. D, 101, 103531 [NASA ADS] [CrossRef] [Google Scholar]
 Di Valentino, E., Anchordoqui, L. A., Akarsu, Ö., et al. 2021a, Astropart. Phys., 131, 102605 [NASA ADS] [CrossRef] [Google Scholar]
 Di Valentino, E., Anchordoqui, L. A., Akarsu, Ö., et al. 2021b, Astropart. Phys., 131, 102604 [NASA ADS] [CrossRef] [Google Scholar]
 Eisenstein, D. J., & Hu, W. 1999, ApJ, 511, 5 [Google Scholar]
 Falco, E. E., Gorenstein, M. V., & Shapiro, I. I. 1985, ApJ, 289, L1 [Google Scholar]
 Flaugher, B. 2005, Int. J. Mod. Phys. A, 20, 3121 [Google Scholar]
 Fu, L., Kilbinger, M., Erben, T., et al. 2014, MNRAS, 441, 2725 [Google Scholar]
 Gatti, M., Jain, B., Chang, C., et al. 2022, Phys. Rev. D, 106, 083509 [NASA ADS] [CrossRef] [Google Scholar]
 GilMarín, H., Wagner, C., Fragkoudi, F., Jimenez, R., & Verde, L. 2012, JCAP, 2012, 047 [CrossRef] [Google Scholar]
 Gruen, D., Friedrich, O., Krause, E., et al. 2018, Phys. Rev. D, 98, 023507 [Google Scholar]
 Halder, A., & Barreira, A. 2022, MNRAS, 515, 4639 [NASA ADS] [CrossRef] [Google Scholar]
 Halder, A., Friedrich, O., Seitz, S., & Varga, T. N. 2021, MNRAS, 506, 2780 [CrossRef] [Google Scholar]
 Hamilton, A. J. S. 2000, MNRAS, 312, 257 [NASA ADS] [CrossRef] [Google Scholar]
 Handley, W. J., Hobson, M. P., & Lasenby, A. N. 2015, MNRAS, 453, 4384 [Google Scholar]
 HarnoisDéraps, J., Amon, A., Choi, A., et al. 2018, MNRAS, 481, 1337 [Google Scholar]
 HarnoisDéraps, J., Martinet, N., Castro, T., et al. 2021, MNRAS, 506, 1623 [CrossRef] [Google Scholar]
 Heydenreich, S., Brück, B., & HarnoisDéraps, J. 2021, A&A, 648, A74 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Heymans, C., Tröster, T., Asgari, M., et al. 2021, A&A, 646, A140 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hikage, C., Oguri, M., Hamana, T., et al. 2019, PASJ, 71, 43 [Google Scholar]
 Hilbert, S. J. 2008, Ph.D. Thesis, LMU Munich, Germany [Google Scholar]
 Hilbert, S., Hartlap, J., White, S. D. M., & Schneider, P. 2009, A&A, 499, 31 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hildebrandt, H., Viola, M., Heymans, C., et al. 2017, MNRAS, 465, 1454 [Google Scholar]
 Hildebrandt, H., Köhlinger, F., van den Busch, J. L., et al. 2020, A&A, 633, A69 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hoekstra, H., & Jain, B. 2008, Ann. Rev. Nucl. Part. Sci., 58, 99 [NASA ADS] [CrossRef] [Google Scholar]
 Huterer, D., Takada, M., Bernstein, G., & Jain, B. 2006, MNRAS, 366, 101 [Google Scholar]
 Ivezic, Z., Axelrod, T., Brandt, W. N., et al. 2008, Serb. Astron. J., 176, 1 [Google Scholar]
 Jarvis, M., Bernstein, G., & Jain, B. 2004, MNRAS, 352, 338 [Google Scholar]
 Jeffrey, N., Lanusse, F., Lahav, O., & Starck, J.L. 2020, MNRAS, 492, 5023 [Google Scholar]
 Joachimi, B., Shi, X., & Schneider, P. 2009, A&A, 508, 1193 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Joudaki, S., Hildebrandt, H., Traykova, D., et al. 2020, A&A, 638, L1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kaiser, N., & Jaffe, A. 1997, ApJ, 484, 545 [NASA ADS] [CrossRef] [Google Scholar]
 Kaiser, N., & Squires, G. 1993, ApJ, 404, 441 [Google Scholar]
 Kaiser, N., Squires, G., & Broadhurst, T. 1995, ApJ, 449, 460 [Google Scholar]
 Kayo, I., Takada, M., & Jain, B. 2013, MNRAS, 429, 344 [Google Scholar]
 Kilbinger, M., & Schneider, P. 2005, A&A, 442, 69 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kilbinger, M., Schneider, P., & Eifler, T. 2006, A&A, 457, 15 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Laureijs, R., Amiaux, J., Arduini, S., et al. 2011, ArXiv eprints [arXiv:1110.3193] [Google Scholar]
 Limber, D. N. 1954, ApJ, 119, 655 [NASA ADS] [CrossRef] [Google Scholar]
 Linke, L., Simon, P., Schneider, P., & Hilbert, S. 2020, A&A, 634, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Martinet, N., Schneider, P., Hildebrandt, H., et al. 2018, MNRAS, 474, 712 [Google Scholar]
 Munshi, D., Valageas, P., van Waerbeke, L., & Heavens, A. 2008, Phys. Rep., 462, 67 [NASA ADS] [CrossRef] [Google Scholar]
 Ogata, H. 2005, Publ. Res. Inst. Math. Sci., 41, 949 [Google Scholar]
 Peebles, P. J. E. 1980, The Largescale Structure of the Universe (Princeton University Press) [Google Scholar]
 Percival, W. J., Friedrich, O., Sellentin, E., & Heavens, A. 2022, MNRAS, 510, 3207 [NASA ADS] [CrossRef] [Google Scholar]
 Pires, S., Vandenbussche, V., Kansal, V., et al. 2020, A&A, 638, A141 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration VI. 2020, A&A, 641, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Porth, L., & Smith, R. E. 2021, MNRAS, 508, 3474 [NASA ADS] [CrossRef] [Google Scholar]
 Porth, L., Smith, R. E., Simon, P., Marian, L., & Hilbert, S. 2020, MNRAS, 499, 2474 [CrossRef] [Google Scholar]
 Pyne, S., & Joachimi, B. 2021, MNRAS, 503, 2300 [NASA ADS] [CrossRef] [Google Scholar]
 Pyne, S., Tenneti, A., & Joachimi, B. 2022, MNRAS, 516, 1829 [Google Scholar]
 Sato, M., & Nishimichi, T. 2013, Phys. Rev. D, 87, 123538 [NASA ADS] [CrossRef] [Google Scholar]
 Schneider, P. 1996, MNRAS, 283, 837 [Google Scholar]
 Schneider, P. 2003, A&A, 408, 829 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schneider, P., & Lombardi, M. 2003, A&A, 397, 809 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schneider, P., & Seitz, C. 1995, A&A, 294, 411 [NASA ADS] [Google Scholar]
 Schneider, P., van Waerbeke, L., Jain, B., & Kruse, G. 1998, MNRAS, 296, 873 [NASA ADS] [CrossRef] [Google Scholar]
 Schneider, P., van Waerbeke, L., & Mellier, Y. 2002, A&A, 389, 729 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schneider, P., Kilbinger, M., & Lombardi, M. 2005, A&A, 431, 9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Scoccimarro, R., & Couchman, H. M. P. 2001, MNRAS, 325, 1312 [NASA ADS] [CrossRef] [Google Scholar]
 Secco, L. F., Jarvis, M., Jain, B., et al. 2022, Phys. Rev. D, 105, 103537 [NASA ADS] [CrossRef] [Google Scholar]
 Seitz, S., & Schneider, P. 2001, A&A, 374, 740 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Seitz, S., Schneider, P., & Bartelmann, M. 1998, A&A, 337, 325 [NASA ADS] [Google Scholar]
 Sellentin, E., & Heavens, A. F. 2016, MNRAS, 456, L132 [Google Scholar]
 Semboloni, E., Hoekstra, H., & Schaye, J. 2013, MNRAS, 434, 148 [CrossRef] [Google Scholar]
 SevillaNoarbe, I., Bechtol, K., Carrasco Kind, M., et al. 2021, ApJS, 254, 24 [NASA ADS] [CrossRef] [Google Scholar]
 Shi, X., Joachimi, B., & Schneider, P. 2014, A&A, 561, A68 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Simon, P., Watts, P., Schneider, P., et al. 2008, A&A, 479, 655 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Springel, V., Yoshida, N., & White, S. D. M. 2001, Nature, 6, 79 [Google Scholar]
 Springel, V., White, S. D. M., Jenkins, A., et al. 2005, Nature, 435, 629 [Google Scholar]
 Spurio Mancini, A., Piras, D., Alsing, J., Joachimi, B., & Hobson, M. P. 2022, MNRAS, 511, 1771 [NASA ADS] [CrossRef] [Google Scholar]
 Takada, M., & Jain, B. 2004, MNRAS, 348, 897 [Google Scholar]
 Takahashi, R., Sato, M., Nishimichi, T., Taruya, A., & Oguri, M. 2012, ApJ, 761, 152 [Google Scholar]
 Takahashi, R., Hamana, T., Shirasaki, M., et al. 2017, ApJ, 850, 24 [Google Scholar]
 Takahashi, R., Nishimichi, T., Namikawa, T., et al. 2020, ApJ, 895, 113 [NASA ADS] [CrossRef] [Google Scholar]
 Troxel, M. A., MacCrann, N., Zuntz, J., et al. 2018, Phys. Rev. D, 98, 043528 [Google Scholar]
 Watkinson, C. A., Majumdar, S., Pritchard, J. R., & Mondal, R. 2017, MNRAS, 472, 2436 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Testing the BIHALOFIT bispectrum model
A.1. Measuring the bispectrum
To measure the convergence bispectrum B_{κ} from simulations, we adapt the estimator developed by Watkinson et al. (2017). While their algorithm has been presented for threedimensional density fields, it can be adapted to twodimensional convergence fields.
Given a convergence field κ(θ) and its Fourier transform , for an ℓbin we define as
and as its inverse Fourier transform. We also define as the inverse Fourier transform of with defined as in Eq. (A.1):
The estimator for the convergence bispectrum is then defined as
where Ω is the solid angle of the respective field and N_{pix} is the number of pixels. The advantage of this estimator is its speed: With seven Fourier transforms, we can extract the complete averaged bispectrum of a field, where the three Fourier transforms required for the computation of only need to be performed once, even when computing the bispectra of multiple fields. Furthermore, can be stored for computing bispectra of different triangle configurations containing the same bin. Nevertheless, the estimator suffers from one significant drawback: The Fourier transform assumes periodicity of the field κ, which is normally not given for convergence maps (in contrast to the threedimensional Nbody simulation cubes, which usually exhibit periodic boundary conditions). Therefore the estimator can be biased for ℓscales approaching the scales of either the field or individual pixels. We, therefore, discard scales smaller than 5 pixels or larger than a third of the field size.
A.2. Validation
To validate our implementation of the BIHALOFIT algorithm and the Limber integration, we compare our model bispectrum with one extracted from the MS. For all 64 linesofsight, we take a convergence map at redshift z = 1 and use the estimator described in Sect. A.1 to extract the bispectrum for a set of triangle configurations. We compare these to our model predictions in Fig. A.1. We recover the bispectrum signal quite well, although we fall short of the 1020% accuracy reported in Takahashi et al. (2020). However, that may very well be due to the limited sample size provided in the MS or the smoothing inside the MS induced by the raytracing (Hilbert et al. 2009).
Fig. A.1. Top panels: convergence bispectrum at z = 1 extracted from the MS (orange) and modelled from BIHALOFIT (blue). Bottom panels: ratio of measurement over the model. The error bars denote the error on the mean of all 64 lines of sight. For each point, we compute the average in a bin of width Δℓ = 0.13ℓ, both in the measurements and the model. For the model, we use equation (8) of Joachimi et al. (2009) as a weight for different triangle configurations within one bin. 
Appendix B: Testing the integration routine for Γ_{i}
When we define the threepoint correlation function of the deflection potential as
we can analytically compute both the threepoint correlation functions and the bispectrum. For this, we define
and use the relations
Defining x = X − Z and y = Y − Z, the following equations hold:
In the last equation, marked by ( * ), the variables x,y are interpreted as complex numbers x = x_{1} + ix_{2}; for their scalar product, x · y = x_{1}y_{1} + x_{2}y_{2} holds. The equation for ⟨γγγ^{*}⟩ was computed via MATHEMATICA and is too long to denote here. We can now set
We then transform the Cartesian shears in Eq. (B.8) to the orthocenters using Eq. (14), and use the fact that Γ^{(0)} = ⟨γγγ⟩ and Γ^{(3)} = ⟨γγγ^{*}⟩ (see Eq. 11), to test our integration routine using an analytic model. As can be seen in Fig. B.1, the integration routine is accurate to the subpercent level.
Fig. B.1. Real (red) and imaginary (blue) parts of the first natural component of the shear threepoint correlation function Γ^{(0)}. Each plot panel corresponds to one fixed triangle shape shown in the bottomright corner; the xaxis represents the length of the red triangle shape. We compare an integration of the analytic bispectrum model (B.9, dark solid) and the analytic solution (B.8, light dashed). Under each panel we show the difference between the integration and the analytic solution, multiplied by a factor of 10^{5} for visibility. 
Appendix C: Additional figures
Fig. C.2. Accuracy of the emulated aperture statistics model data vector m(Θ) for the KiDS1000 setup. Given the bispectrum accuracy of 10% (Takahashi et al. 2020), the emulator uncertainty plays a negligible role in the modelling process. 
Fig. C.3. Comparison of the measured KiDS1000like data vector in the T17 simulations to the modelled vector of the second and thirdorder aperture statistics. The orange band is the expected KiDS1000 uncertainty. 
All Tables
All Figures
Fig. 1. Schematic diagram of the modelling and validation pipeline introduced in this paper. The numbers on the arrows correspond to the respective section where this part of the pipeline is discussed. 

In the text 
Fig. 2. Redshift distribution constructed from the T17 simulation given the fiducial n(z) of the KiDS1000 data. 

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

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

In the text 
Fig. 5. Thirdorder aperture mass statistics when modelled directly from the bispectrum (Eq. (31), blue line) compared to the ones constructed from modelled shear 3pcf (Eq. (35), coloured crosses). The dotted lines denote the respective Bmodes. We plot (θ_{1}, θ_{2}, θ_{3}), where θ_{1} is constant in each row, θ_{2} is constant in each column and θ_{3} varies along the xaxis in each panel. The yellow, green, red and purple dots denote the we get from the shear 3pcf when we bin r, u and v in 7 × 7 × 7, 10 × 10 × 10, 15 × 15 × 15 and 20 × 20 × 20 bins, respectively. 

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

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

In the text 
Fig. 8. Comparison of the posterior distributions of different thirdorder statistics using the SLICS simulations to determine the covariance. In red we show the posteriors of the PCA from the 3pcf and in blue the one resulting from the statistic. 

In the text 
Fig. 9. Left panel: Comparison of the posterior distributions of second and thirdorder aperture statistics with a joint analysis within a KiDS1000like setup. The data vector and covariance matrix are estimated by the T17 simulations, and we use filter scales of (4, 8, 16, 32) arcmin. The Hubble parameter h and the dark energy equationofstate are fixed to the T17 values. Right panel: Comparison of the combined second and thirdorder aperture statistics with the combination of secondorder aperture statistics with the equalscale thirdorder aperture statistics. The corresponding data vector and model vector can be found in Fig. C.3. 

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

In the text 
Fig. B.1. Real (red) and imaginary (blue) parts of the first natural component of the shear threepoint correlation function Γ^{(0)}. Each plot panel corresponds to one fixed triangle shape shown in the bottomright corner; the xaxis represents the length of the red triangle shape. We compare an integration of the analytic bispectrum model (B.9, dark solid) and the analytic solution (B.8, light dashed). Under each panel we show the difference between the integration and the analytic solution, multiplied by a factor of 10^{5} for visibility. 

In the text 
Fig. C.1. Same as Fig. 4, just for the second natural component Γ^{(1)}. 

In the text 
Fig. C.2. Accuracy of the emulated aperture statistics model data vector m(Θ) for the KiDS1000 setup. Given the bispectrum accuracy of 10% (Takahashi et al. 2020), the emulator uncertainty plays a negligible role in the modelling process. 

In the text 
Fig. C.3. Comparison of the measured KiDS1000like data vector in the T17 simulations to the modelled vector of the second and thirdorder aperture statistics. The orange band is the expected KiDS1000 uncertainty. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.