Open Access
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/0004-6361/202244820
Published online 27 March 2023

© The Authors 2023

Licence Creative CommonsOpen 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, six-parameter model describes a wide range of observations, from the cosmic microwave background (CMB) to the observed large-scale structure (LSS) of galaxies, with remarkable accuracy. As the reported uncertainties on cosmological parameters are approaching the per-cent level, a few tensions arise between CMB observations of the early Universe and observations of the local Universe that quantify the LSS (for example, in the Hubble parameter H0, 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, two-point statistics have been established as the main analysis tool for cosmic shear (Schneider et al. 1998; Troxel et al. 2018; Hildebrandt et al. 2017, 2020; Hikage et al. 2019; Asgari et al. 2020). These statistics capture the entire information content of a Gaussian random field. Since the initial density field of the Universe is believed to be Gaussian, two-point statistics capture a large amount of cosmological information. However, in late times non-linear structure formation has introduced non-Gaussian features at the smaller scales of the matter distribution, whose information content cannot be captured by two-point statistics. To use this information, a variety of higher-order statistics has been introduced in recent years, including peak count statistics (Martinet et al. 2018; Harnois-Déraps et al. 2021), persistent homology (Heydenreich et al. 2021), density split statistics (Gruen et al. 2018; Burger et al. 2022), and many others. In this work, we consider third-order shear statistics, which measure the skewness of the LSS at various scales.

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

This article aims to prepare a cosmological parameter analysis with cosmic shear three-point statistics by developing a pipeline to measure and model both the three-point correlation functions Γ(i) and third-order 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 Stage-III survey, which constitutes a significant advantage over their Fourier-space counterpart, the convergence bispectrum. The aperture mass statistics have several additional advantages: They offer good data compression, are not subject to the mass-sheet degeneracy, and, most importantly, they decompose the signal into E- and B-modes, where to leading order only E-modes can be created by gravitational lensing.

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

thumbnail Fig. 1.

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

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

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

2.1. T17 simulations

The T17 are used in this work to perform a realistic analysis of a survey that mimics the KiDS-1000 data and are constructed from a series of nested cubic boxes with side lengths of L, 2L, 3L... placed around a fixed vertex representing the observer’s position, with L = 450 Mpc/h. Each box is replicated eight times and placed around the observer using periodic boundary conditions. With the N-body code GADGET2 (Springel et al. 2001) the gravitational evolution of 20483 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 GRAYTRIX2 to trace the light-ray trajectories from the observer to the last scattering surface3. The cosmological parameters of the simulation are Ωm = 1 − ΩΛ = 0.279, Ωb = 0.046, h = 0.7, σ8 = 0.82, and ns = 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 KiDS-1000 n(z) – see Fig. 2.

thumbnail Fig. 2.

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

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

(1)

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

2.2. Millennium simulations

The MS were run with 21603 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 ns = 1. Subsequently, shear- and convergence-maps of 64 independent lines of sight with an area of 4 × 4 deg2 each were created at 36 different redshifts between z = 0 and z = 3 (Hilbert 2008; Hilbert et al. 2009). Each map is calculated on a grid of 4096 × 4096 pixel. For each line of sight, we use the full shear- and convergence maps at redshift z = 1. As we use the MS solely to validate our model, we do not add any noise to the maps.

2.3. Scinet-LIghtCone Simulations

The SLICS were run with 15363 particles in a 505 h−1 Mpc box, filling up 10 × 10 deg2 light-cones 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 ns = 0.969. From these simulations, we use convergence maps as well as galaxy catalogues with a redshift distribution of

(2)

with z0 = 0.637, β = 1.5 and the overall proportionality constant given by normalising the distribution to 30 gal/arcmin2. We use mock galaxy catalogues provided for 924 (pseudo-)independent lines of sight. The SLICS are useful to estimate the constraining power of third-order statistics since the three-point correlation function can be calculated relatively quickly on the 100 deg2 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 line-of-sight integration, weighted by the lensing efficiency

(3)

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

3.1. Definition of power- and bispectrum

The matter power spectrum Pδ(k, z) and bispectrum Bδ(k1, k2, k3, z) can be defined as

(4)

(5)

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

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

(6)

(7)

Here,

(8)

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

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

Instead of the three -values 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

(9)

3.2. Modelling the bispectrum

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

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

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

4. Shear three-point correlation functions

4.1. Definition of the shear three-point correlation functions

Shear three-point correlation functions (3pcf) are the natural extension to the widely used two-point correlation functions. Let γc = γ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 cross-components with respect to a point fixed with respect to the triangle, for example one of its centres,

(10)

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

(11)

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

thumbnail Fig. 3.

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

Parametrising the shear 3pcf by the triangle side-lengths, x1, x2 and x3, where the indices 1, 2 and 3 are ordered in a counter-clockwise direction, the natural components exhibit a nice behaviour concerning cyclic permutation of arguments. While the first natural component Γ(0) is invariant under cyclic permutations, the other three components transform into each other,

(12)

Similar behaviour can be observed for parity transformations:

(13)

4.2. Modelling the shear three-point correlation functions

We model the natural components of the three-point correlation functions of cosmic shear using the methods described in Schneider et al. (2005, hereafter S+05): When we project the shear of all three galaxies to the orthocenter of the triangle, the projection direction at the triangle vertex Xi is orthogonal to the orientation φi of the triangle side xi. Thus, the shear transforms as

(14)

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

(15)

where β is the polar angle of , to write

(16)

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

(17)

with

(18)

The quantities A1, 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

(19)

with

(20)

Defining

(21)

and E1 and E2 via cyclic permutations of indices, we can write

(22)

The R-integration filters the bispectrum with a 6-th order Bessel function, making the integration routine difficult to solve numerically, as the functional form of the bispectrum prevents the application of fast Hankel transform algorithms like FFTLog (Hamilton 2000)5. We thus use the method developed in Ogata (2005) to solve the R-integration and integrate the remaining dimensions using the CUBATURE library6. 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 three-point correlation function can be described by its three side-lengths x1, x2 and x3, it is certainly not a good idea to use these variables for a binning scheme; for example, for x1 > x2 + x3 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 x1 > x2 > x3 they defined a triangle via the values r ∈ [0, ∞], u ∈ [0, 1] and v ∈ [ − 1, 1] by

(23)

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

We note that Eq. (12) implies that the four shear 3pcf for x1 > x2 > x3 (as in the binning scheme of Jarvis et al. 2004) already contain the entire information content of the third-order shear signal, as does knowledge of Γ(0) and Γ(1) for all combinations of x1, x2 and x3. 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 three-point correlation functions face a few hurdles: Assuming we bin the three-point correlation functions in 10 bins for each r, u and v, then our data vector consists of 8000 entries. This makes estimating a covariance matrix using simulations practically impossible and leads to a modelling time that is unfeasible even for a non-tomographic analysis.

4.3. Measuring the shear three-point correlation functions

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

(24)

where the εi = εt, i + iε×, i are the observed ellipticities of galaxies and wi 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 Ngal ≳ 106. That is why TREECORR constructs a hierarchical ball tree out of the galaxy sample and calculates the correlation functions from this tree. This results in a remarkable speed-up and allows us to calculate the shear 3pcf for an ensemble of about 107 source galaxies distributed over a 10 × 10 deg2 field in about 1500 CPUh. A disadvantage of the tree-code is that its execution time scales massively with the number of bins: If b is the logarithmic bin size, then the run-time scales roughly with b−4.

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

4.4. Validation

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

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

thumbnail Fig. 4.

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

5. Aperture mass statistics

5.1. Definition of aperture mass statistics

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

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

(25)

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

(26)

where γt and γ× are projected along the vector ϑ′−ϑ (see Eq. (10)) and Qθ is related to Uθ via

(27)

For simplicity of notation, we define Uθ(ϑ) = θ−2u(ϑ/θ) and denote by the Fourier transform of u. In this work, we opt for the filter function introduced by Crittenden et al. (2002),

(28)

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

(29)

By construction, ⟨ℳap⟩(θ) vanishes. In a parity-symmetric field, all odd powers of B-mode components vanish (Schneider 2003), meaning that the relevant B-mode counterparts to and are and , respectively.

5.2. Modelling aperture mass statistics

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

(30)

As a model for the non-linear power spectrum, we use the revised HALOFIT model of Takahashi et al. (2012). Equivalently, the third-order aperture statistics can be derived from a bispectrum model via (compare Jarvis et al. 2004; Schneider et al. 2005)8

(31)

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

5.3. Measuring aperture mass statistics

There are three main methods to estimate the third-order aperture mass statistics , first, via the convergence field κ, second, via the shear field or, in practice, from the observed galaxy ellipticities, and third, via the third-order correlation functions Γ(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 maps9.

To extract the data vectors from the (full-sky) convergence maps of the T17 simulations, we smooth the maps with HEALPY function SMOOTHING, with a given beam window function created by the function BEAM2BL, which in turn is determined by the corresponding Uθ filter. For each filter radius θ, this yields a full-sky aperture mass map , 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

(32)

where εt/× are the observed galaxy ellipticities converted into their tangential and croos components according to Eq. (10); ϑi are their respective positions. Here, ngal 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 ngal as the number of galaxies weighted by the Q-filter function:

(33)

We tested all three definitions of ngal using the SLICS and found that, for randomly distributed galaxies, setting ngal as the number density within the aperture radius or the one weighted by the Q-filter function induces sub-percent differences on the third-order aperture masses . However, setting ngal as the global galaxy density can induce differences of about 5% in .

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

(34)

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

5.3.3. Measuring aperture statistics from three-point correlation functions

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

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

(35)

need to be computed. Here, A1, 2 and F1, 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:

(36)

We do not consider the quantities and , as they vanish for any parity-symmetric 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 F1, 2 at the bin centre. While numerically very cheap, this is probably not the most efficient way to compute that integral. In case the filter functions vary significantly over a bin (compare Figs. 2 and 3 of Schneider et al. 2005), it might be more appropriate to calculate the average of the filter function in a bin, for example. We tried to improve the integration results by interpolating the measured shear 3pcf and performing the same integral, achieving rather moderate improvements. We leave an optimisation of the conversion from shear 3pcf to aperture mass statistics for future work.

5.4. Validation

5.4.1. Binning choice of three-point correlation functions

As a first step, we want to validate the conversion Γ(i) performed by the TREECORR algorithm. In particular, we investigate the number of bins necessary to achieve an unbiased estimate for third-order aperture statistics and quantify the leakage of E- and B-modes. While the latter has already been investigated by Shi et al. (2014), we extend upon these results by using a realistic convergence bispectrum model and by taking into account non-diagonal aperture mass statistics.

Our modelling pipeline provides us with the ability to test this conversion. As we start from the same convergence bispectrum Bκ, the aperture mass statistics achieved by the conversion and by direct modelling have to be consistent. Furthermore, the modelled Γ(i) are pure E-mode functions, so any B-modes 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 73 bins do we get significant deviations, meaning that this is certainly not a sufficient number of bins. The conversion becomes less accurate when two filter radii are much smaller than the third one, as shown in the top-right corner of Fig. 5. There also the results for 103 bins show significant deviations, whereas the results for 153 bins seem overall consistent with the ones from 203 bins. We also observe a non-negligible amount of B-mode leakage for these cases, even for 153 and 203 bins. We are planning to exclude all combinations of filter radii from cosmological parameter analyses where we observe a B-mode leakage of more than 10% for the 3pcf in 153 bins.

thumbnail Fig. 5.

Third-order 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 B-modes. We plot (θ1, θ2, θ3), where θ1 is constant in each row, θ2 is constant in each column and θ3 varies along the x-axis in each panel. The yellow, green, red and purple dots denote the 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 well-behaved with respect to r, but strongly varies as a function of u and v, especially when u ≈ 0 and v ≈ ±1. This implies that r can be binned rather coarsely, as long as u and v are finely binned. This is in contrast to binning choices in other studies, for example Secco et al. (2022), who preferred a fine binning in r (55 bins) and a coarser binning in u and v (10 bins).

thumbnail Fig. 6.

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

5.4.2. Comparison to N-body simulations

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

thumbnail Fig. 7.

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

6. Cosmological parameter estimation

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

(37)

where

(38)

The power-law index m is

(39)

with nd being the number of data points and

(40)

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

Normally, one needs to evaluate this likelihood function in high-dimensional parameter space to perform a cosmological parameter analysis with third-order statistics. The creation of the model vector for third-order aperture statistics with 35 combinations of aperture radii takes about one minute on an NVIDIA A40 GPU, and needs to be evaluated at about 104 points even for sampling methods like POLYCHORD (Handley et al. 2015). This means that a complete likelihood analysis with third-order aperture statistics is possible but takes a long time, whereas a complete likelihood analysis with the 3pcf, where the modelling of the 3pcf in 103 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 points12 in a four-dimensional Latin hypercube describing a flat wCDM cosmological model, varying the parameters Ωm, S8, w0 and h with flat priors of Ωm ∈ [0.1, 0.55], S8 ∈ [0.6, 0.9],w0 ∈ [ − 2, −0.5] and h ∈ [0.6, 0.82]. We leave all other parameters fixed at the values corresponding to the SLICS, which we introduced in Sect. 2. We then train the emulator with 6500 points and use the remaining 1000 as a validation test, as shown in Fig. C.2. This neural network-based emulator performs extraordinarily well, modelling all statistics with a sub-per-cent accuracy. As our bispectrum model is only accurate to about 10% (Takahashi et al. 2020), the emulator uncertainty plays a negligible role in the modelling process. This emulator yields extremely fast predictions for the third-order statistics, so that we can perform a full likelihood analysis using the Metropolis-Hastings sampling method.

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

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

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

Unfortunately, we cannot quantify the full information content of the shear 3pcf, as the data vector contains about 104 entries, and therefore a reliable covariance matrix is not accessible. To circumvent this problem, we model the shear 3pcf in 103 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 vector13.

Using this PCA, we determine the covariance matrix for the shear 3pcf, which we measure from 200 10 × 10 deg2 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.

thumbnail Fig. 8.

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

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

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

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

The resulting posterior distribution is shown in the left panel of Fig. 9, where we first notice that the combination of second- and third-order statistics significantly increases the constraining power, especially in the Ωm-σ8 panel to the different degeneracy directions of the individual statistics (compare Table 1). Indeed, a joint analysis increases the constraints on S8 by 42% with respect to second-order statistics; the constraints on Ωm and σ8 increase by at least 68% and 54%, respectively14. 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 KiDS-1000 uncertainty for all three statistics.

thumbnail Fig. 9.

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

Table 1.

Marginalised one-dimensional parameter constraints from the left part of Fig. 9.

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

7. Discussion

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

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

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

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

Finally, we compare the constraining power between second-order, third-order, and joint aperture statistics analysis. While second-order aperture statistics are not being used in modern cosmological parameter analyses, we assume that all second-order shear statistics exhibit similar constraining power and parameter degeneracies (compare Asgari et al. 2021). We find that third-order aperture statistics alone have a lower constraining power than their second-order counterpart, but they exhibit a different degeneracy direction in the Ωmσ8-plane so that a joint analysis almost doubles the constraining power on the structure growth parameter S8 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 third-order aperture statistics (compare Kilbinger & Schneider 2005), appears overly optimistic. This suggests that a Fisher forecast might not be an optimal tool to forecast parameter constraints and mock sampling methods give more realistic constraints.

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


3

These maps are freely available for download at http://cosmo.phys.hirosaki-u.ac.jp/takahasi/allsky_raytracing/

4

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.

5

FFTLog requires the integral to be of the form g(k) = ∫drf(r)Jν(kr), whereas in our case we have to compute a three-dimensional integral where the bispectrum depends non-linearly on all three integration variables.

7

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 galaxy-galaxy-galaxy-lensing correlation function.

8

Again, due to the different definitions of the convergence bispectrum, we get a factor of 3 difference with respect to S+05. We also use the symmetry of the bispectrum to only integrate from 0 to π in φ, introducing a prefactor of 2.

9

Both the Q- and the u-filter function have 99.9% of their power within this range, meaning that boundary effects beyond this cut-off are negligible.

10

For a 10 × 10 deg2 field of a Stage-IV survey, direct estimation of aperture mass statistics takes a few min. vs. ∼1500 CPUh for the estimation of the 3pcf.

11

Again, we only use every tenth pixel to calculate the 3pcf, but do not expect the a strong loss of signal to noise from this.

12

Due to the long modelling time, we arbitrarily pick 500 of the 7500 sampled points for the 3pcf and verify by eye that they uniformly sample the prior range.

13

We found that using more principal components does not improve the constraints and just leads to a noisier covariance matrix.

14

As the constraints for second-order statistics are at least partly dominated by the prior, the true increase is likely to be much greater.

Acknowledgments

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

References

  1. Abbott, T. M. C., Aguena, M., Alarcon, A., et al. 2022, Phys. Rev. D, 105, 023520 [CrossRef] [Google Scholar]
  2. Albrecht, A., Bernstein, G., Cahn, R., et al. 2006, ArXiv e-prints [arXiv:astro-ph/0609591] [Google Scholar]
  3. Asgari, M., Heymans, C., Hildebrandt, H., et al. 2019, A&A, 624, A134 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  4. Asgari, M., Tröster, T., Heymans, C., et al. 2020, A&A, 634, A127 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  5. Asgari, M., Lin, C.-A., Joachimi, B., et al. 2021, A&A, 645, A104 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  6. Bartelmann, M. 2010, CQG, 27, 233001 [NASA ADS] [CrossRef] [Google Scholar]
  7. Bartelmann, M., & Schneider, P. 2001, Phys. Rep., 340, 291 [Google Scholar]
  8. Bernardeau, F., van Waerbeke, L., & Mellier, Y. 1997, A&A, 322, 1 [NASA ADS] [Google Scholar]
  9. Burger, P., Friedrich, O., Harnois-Déraps, J., & Schneider, P. 2022, A&A, 661, A137 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  10. Crittenden, R. G., Natarajan, P., Pen, U.-L., & Theuns, T. 2002, ApJ, 568, 20 [NASA ADS] [CrossRef] [Google Scholar]
  11. Deshpande, A. C., & Kitching, T. D. 2020, Phys. Rev. D, 101, 103531 [NASA ADS] [CrossRef] [Google Scholar]
  12. Di Valentino, E., Anchordoqui, L. A., Akarsu, Ö., et al. 2021a, Astropart. Phys., 131, 102605 [NASA ADS] [CrossRef] [Google Scholar]
  13. Di Valentino, E., Anchordoqui, L. A., Akarsu, Ö., et al. 2021b, Astropart. Phys., 131, 102604 [NASA ADS] [CrossRef] [Google Scholar]
  14. Eisenstein, D. J., & Hu, W. 1999, ApJ, 511, 5 [Google Scholar]
  15. Falco, E. E., Gorenstein, M. V., & Shapiro, I. I. 1985, ApJ, 289, L1 [Google Scholar]
  16. Flaugher, B. 2005, Int. J. Mod. Phys. A, 20, 3121 [Google Scholar]
  17. Fu, L., Kilbinger, M., Erben, T., et al. 2014, MNRAS, 441, 2725 [Google Scholar]
  18. Gatti, M., Jain, B., Chang, C., et al. 2022, Phys. Rev. D, 106, 083509 [NASA ADS] [CrossRef] [Google Scholar]
  19. Gil-Marín, H., Wagner, C., Fragkoudi, F., Jimenez, R., & Verde, L. 2012, JCAP, 2012, 047 [CrossRef] [Google Scholar]
  20. Gruen, D., Friedrich, O., Krause, E., et al. 2018, Phys. Rev. D, 98, 023507 [Google Scholar]
  21. Halder, A., & Barreira, A. 2022, MNRAS, 515, 4639 [NASA ADS] [CrossRef] [Google Scholar]
  22. Halder, A., Friedrich, O., Seitz, S., & Varga, T. N. 2021, MNRAS, 506, 2780 [CrossRef] [Google Scholar]
  23. Hamilton, A. J. S. 2000, MNRAS, 312, 257 [NASA ADS] [CrossRef] [Google Scholar]
  24. Handley, W. J., Hobson, M. P., & Lasenby, A. N. 2015, MNRAS, 453, 4384 [Google Scholar]
  25. Harnois-Déraps, J., Amon, A., Choi, A., et al. 2018, MNRAS, 481, 1337 [Google Scholar]
  26. Harnois-Déraps, J., Martinet, N., Castro, T., et al. 2021, MNRAS, 506, 1623 [CrossRef] [Google Scholar]
  27. Heydenreich, S., Brück, B., & Harnois-Déraps, J. 2021, A&A, 648, A74 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  28. Heymans, C., Tröster, T., Asgari, M., et al. 2021, A&A, 646, A140 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  29. Hikage, C., Oguri, M., Hamana, T., et al. 2019, PASJ, 71, 43 [Google Scholar]
  30. Hilbert, S. J. 2008, Ph.D. Thesis, LMU Munich, Germany [Google Scholar]
  31. Hilbert, S., Hartlap, J., White, S. D. M., & Schneider, P. 2009, A&A, 499, 31 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  32. Hildebrandt, H., Viola, M., Heymans, C., et al. 2017, MNRAS, 465, 1454 [Google Scholar]
  33. Hildebrandt, H., Köhlinger, F., van den Busch, J. L., et al. 2020, A&A, 633, A69 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  34. Hoekstra, H., & Jain, B. 2008, Ann. Rev. Nucl. Part. Sci., 58, 99 [NASA ADS] [CrossRef] [Google Scholar]
  35. Huterer, D., Takada, M., Bernstein, G., & Jain, B. 2006, MNRAS, 366, 101 [Google Scholar]
  36. Ivezic, Z., Axelrod, T., Brandt, W. N., et al. 2008, Serb. Astron. J., 176, 1 [Google Scholar]
  37. Jarvis, M., Bernstein, G., & Jain, B. 2004, MNRAS, 352, 338 [Google Scholar]
  38. Jeffrey, N., Lanusse, F., Lahav, O., & Starck, J.-L. 2020, MNRAS, 492, 5023 [Google Scholar]
  39. Joachimi, B., Shi, X., & Schneider, P. 2009, A&A, 508, 1193 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  40. Joudaki, S., Hildebrandt, H., Traykova, D., et al. 2020, A&A, 638, L1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  41. Kaiser, N., & Jaffe, A. 1997, ApJ, 484, 545 [NASA ADS] [CrossRef] [Google Scholar]
  42. Kaiser, N., & Squires, G. 1993, ApJ, 404, 441 [Google Scholar]
  43. Kaiser, N., Squires, G., & Broadhurst, T. 1995, ApJ, 449, 460 [Google Scholar]
  44. Kayo, I., Takada, M., & Jain, B. 2013, MNRAS, 429, 344 [Google Scholar]
  45. Kilbinger, M., & Schneider, P. 2005, A&A, 442, 69 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  46. Kilbinger, M., Schneider, P., & Eifler, T. 2006, A&A, 457, 15 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  47. Laureijs, R., Amiaux, J., Arduini, S., et al. 2011, ArXiv e-prints [arXiv:1110.3193] [Google Scholar]
  48. Limber, D. N. 1954, ApJ, 119, 655 [NASA ADS] [CrossRef] [Google Scholar]
  49. Linke, L., Simon, P., Schneider, P., & Hilbert, S. 2020, A&A, 634, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  50. Martinet, N., Schneider, P., Hildebrandt, H., et al. 2018, MNRAS, 474, 712 [Google Scholar]
  51. Munshi, D., Valageas, P., van Waerbeke, L., & Heavens, A. 2008, Phys. Rep., 462, 67 [NASA ADS] [CrossRef] [Google Scholar]
  52. Ogata, H. 2005, Publ. Res. Inst. Math. Sci., 41, 949 [Google Scholar]
  53. Peebles, P. J. E. 1980, The Large-scale Structure of the Universe (Princeton University Press) [Google Scholar]
  54. Percival, W. J., Friedrich, O., Sellentin, E., & Heavens, A. 2022, MNRAS, 510, 3207 [NASA ADS] [CrossRef] [Google Scholar]
  55. Pires, S., Vandenbussche, V., Kansal, V., et al. 2020, A&A, 638, A141 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  56. Planck Collaboration VI. 2020, A&A, 641, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  57. Porth, L., & Smith, R. E. 2021, MNRAS, 508, 3474 [NASA ADS] [CrossRef] [Google Scholar]
  58. Porth, L., Smith, R. E., Simon, P., Marian, L., & Hilbert, S. 2020, MNRAS, 499, 2474 [CrossRef] [Google Scholar]
  59. Pyne, S., & Joachimi, B. 2021, MNRAS, 503, 2300 [NASA ADS] [CrossRef] [Google Scholar]
  60. Pyne, S., Tenneti, A., & Joachimi, B. 2022, MNRAS, 516, 1829 [Google Scholar]
  61. Sato, M., & Nishimichi, T. 2013, Phys. Rev. D, 87, 123538 [NASA ADS] [CrossRef] [Google Scholar]
  62. Schneider, P. 1996, MNRAS, 283, 837 [Google Scholar]
  63. Schneider, P. 2003, A&A, 408, 829 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  64. Schneider, P., & Lombardi, M. 2003, A&A, 397, 809 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  65. Schneider, P., & Seitz, C. 1995, A&A, 294, 411 [NASA ADS] [Google Scholar]
  66. Schneider, P., van Waerbeke, L., Jain, B., & Kruse, G. 1998, MNRAS, 296, 873 [NASA ADS] [CrossRef] [Google Scholar]
  67. Schneider, P., van Waerbeke, L., & Mellier, Y. 2002, A&A, 389, 729 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  68. Schneider, P., Kilbinger, M., & Lombardi, M. 2005, A&A, 431, 9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  69. Scoccimarro, R., & Couchman, H. M. P. 2001, MNRAS, 325, 1312 [NASA ADS] [CrossRef] [Google Scholar]
  70. Secco, L. F., Jarvis, M., Jain, B., et al. 2022, Phys. Rev. D, 105, 103537 [NASA ADS] [CrossRef] [Google Scholar]
  71. Seitz, S., & Schneider, P. 2001, A&A, 374, 740 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  72. Seitz, S., Schneider, P., & Bartelmann, M. 1998, A&A, 337, 325 [NASA ADS] [Google Scholar]
  73. Sellentin, E., & Heavens, A. F. 2016, MNRAS, 456, L132 [Google Scholar]
  74. Semboloni, E., Hoekstra, H., & Schaye, J. 2013, MNRAS, 434, 148 [CrossRef] [Google Scholar]
  75. Sevilla-Noarbe, I., Bechtol, K., Carrasco Kind, M., et al. 2021, ApJS, 254, 24 [NASA ADS] [CrossRef] [Google Scholar]
  76. Shi, X., Joachimi, B., & Schneider, P. 2014, A&A, 561, A68 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  77. Simon, P., Watts, P., Schneider, P., et al. 2008, A&A, 479, 655 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  78. Springel, V., Yoshida, N., & White, S. D. M. 2001, Nature, 6, 79 [Google Scholar]
  79. Springel, V., White, S. D. M., Jenkins, A., et al. 2005, Nature, 435, 629 [Google Scholar]
  80. Spurio Mancini, A., Piras, D., Alsing, J., Joachimi, B., & Hobson, M. P. 2022, MNRAS, 511, 1771 [NASA ADS] [CrossRef] [Google Scholar]
  81. Takada, M., & Jain, B. 2004, MNRAS, 348, 897 [Google Scholar]
  82. Takahashi, R., Sato, M., Nishimichi, T., Taruya, A., & Oguri, M. 2012, ApJ, 761, 152 [Google Scholar]
  83. Takahashi, R., Hamana, T., Shirasaki, M., et al. 2017, ApJ, 850, 24 [Google Scholar]
  84. Takahashi, R., Nishimichi, T., Namikawa, T., et al. 2020, ApJ, 895, 113 [NASA ADS] [CrossRef] [Google Scholar]
  85. Troxel, M. A., MacCrann, N., Zuntz, J., et al. 2018, Phys. Rev. D, 98, 043528 [Google Scholar]
  86. 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 three-dimensional density fields, it can be adapted to two-dimensional convergence fields.

Given a convergence field κ(θ) and its Fourier transform , for an -bin we define as

(A.1)

and as its inverse Fourier transform. We also define as the inverse Fourier transform of with defined as in Eq. (A.1):

(A.2)

The estimator for the convergence bispectrum is then defined as

(A.3)

where Ω is the solid angle of the respective field and Npix 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 three-dimensional N-body simulation cubes, which usually exhibit periodic boundary conditions). Therefore the estimator can be biased for -scales approaching the scales of either the field or individual pixels. We, therefore, discard scales smaller than 5 pixels or larger than a third of the field size.

A.2. Validation

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

thumbnail Fig. A.1.

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

Appendix B: Testing the integration routine for Γi

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

(B.1)

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

(B.2)

and use the relations

(B.3)

(B.4)

(B.5)

Defining x = XZ and y = YZ, the following equations hold:

(B.6)

(B.7)

(B.8)

In the last equation, marked by ( * ), the variables x,y are interpreted as complex numbers x = x1 + ix2; for their scalar product, x · y = x1y1 + x2y2 holds. The equation for ⟨γγγ*⟩ was computed via MATHEMATICA and is too long to denote here. We can now set

(B.9)

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 sub-percent level.

thumbnail Fig. B.1.

Real (red) and imaginary (blue) parts of the first natural component of the shear three-point correlation function Γ(0). Each plot panel corresponds to one fixed triangle shape shown in the bottom-right corner; the x-axis 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 105 for visibility.

Appendix C: Additional figures

thumbnail Fig. C.1.

Same as Fig. 4, just for the second natural component Γ(1).

thumbnail Fig. C.2.

Accuracy of the emulated aperture statistics model data vector m(Θ) for the KiDS-1000 setup. Given the bispectrum accuracy of 10% (Takahashi et al. 2020), the emulator uncertainty plays a negligible role in the modelling process.

thumbnail Fig. C.3.

Comparison of the measured KiDS-1000-like data vector in the T17 simulations to the modelled vector of the second- and third-order aperture statistics. The orange band is the expected KiDS-1000 uncertainty.

All Tables

Table 1.

Marginalised one-dimensional parameter constraints from the left part of Fig. 9.

All Figures

thumbnail 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
thumbnail Fig. 2.

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

In the text
thumbnail Fig. 3.

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

In the text
thumbnail Fig. 4.

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

In the text
thumbnail Fig. 5.

Third-order 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 B-modes. We plot (θ1, θ2, θ3), where θ1 is constant in each row, θ2 is constant in each column and θ3 varies along the x-axis in each panel. The yellow, green, red and purple dots denote the 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
thumbnail Fig. 6.

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

In the text
thumbnail Fig. 7.

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

In the text
thumbnail Fig. 8.

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

In the text
thumbnail Fig. 9.

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

In the text
thumbnail 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
thumbnail Fig. B.1.

Real (red) and imaginary (blue) parts of the first natural component of the shear three-point correlation function Γ(0). Each plot panel corresponds to one fixed triangle shape shown in the bottom-right corner; the x-axis 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 105 for visibility.

In the text
thumbnail Fig. C.1.

Same as Fig. 4, just for the second natural component Γ(1).

In the text
thumbnail Fig. C.2.

Accuracy of the emulated aperture statistics model data vector m(Θ) for the KiDS-1000 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
thumbnail Fig. C.3.

Comparison of the measured KiDS-1000-like data vector in the T17 simulations to the modelled vector of the second- and third-order aperture statistics. The orange band is the expected KiDS-1000 uncertainty.

In the text

Current usage metrics show cumulative count of Article Views (full-text 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 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.