© ESO 2023

1. Introduction

It is a well-established fact that the Universe is undergoing a phase of accelerated expansion (e.g., Riess et al. 1998; Perlmutter et al. 1999). Understanding what is driving this acceleration in the framework of a spatially flat universe is one of the (if not the) greatest challenges of modern-day cosmology. The concordance Λ cold dark matter (ΛCDM) model performs excellently in fitting the available data, yet the cosmological constant Λ is far from satisfactory from a theoretical point of view. To make things harder for ΛCDM, recent tensions have emerged due to an inconsistency between the values of some parameters from independent data measuring the same quantities in radically different ways. The most debated case is the discrepancy between the Hubble constant H₀ as measured from local probes and as inferred from cosmological data sets (see, e.g., Di Valentino et al. 2021, for a review). Another, albeit less significant, example is the disagreement between the cosmic microwave background (CMB) and lensing estimates of the growth of structure parameter, $S_8={\sigma }_8\sqrt{{\mathrm{\Omega }}_{\mathrm{m}}/0.3}$ (e.g., Hildebrandt et al. 2017; Heymans et al. 2021; Amon et al. 2022). Although some unknown systematic effects could have been missed in the analysis, such tensions may also be the first signs that alternative models are needed, relying on either dark energy in a general relativity framework or based on modified gravity (see, e.g., Joyce et al. 2016, and references therein). Discriminating among the plethora of viable candidates is the aim of Stage IV surveys such as the Dark Energy Spectroscopic Instrument (DESI, DESI Collaboration 2016), the Prime Focus Spectrograph (PFS, Takada et al. 2014), the Vera C. Rubin Observatory Legacy Survey of Space and Time (LSST, Ivezić et al. 2019), Euclid (Laureijs et al. 2011), SPHEREx (Doré et al. 2014), and the Nancy Grace Roman Space Telescope (Spergel et al. 2015).

In this context, the Euclid mission will play a pivotal role in cosmology by measuring the dark energy equation of state w and the S₈ parameter with exquisite precision and accuracy (see, e.g., Laureijs et al. 2011; Euclid Collaboration 2020). This will be achieved by exploiting the clustering of galaxies as well as cosmic shear: the Weak gravitational Lensing (WL) signal of galaxies due to the deflection of light rays by the large-scale structure. In this article, we focus on optimizing the extraction of cosmological information from the cosmic shear probe.

The classical cosmic shear analysis involves measuring the correlations between ellipticities of pairs of galaxies as a function of their separation; an estimator called the two-point correlation functions of the shear (γ-2PCF). Although it benefits from a comprehensive theoretical description, this estimator is only sensitive to the multiscale variance of the lensing field. However, the gravitational collapse of the matter perturbations introduces non-Gaussian features in the shear field. As a consequence, the γ-2PCF and any estimator probing the field up to second order do not contain all cosmological information. This is seen, for example, in the degeneracies between parameters, such as the one between Ω_m and σ₈, which means the γ-2PCF can only efficiently constrain their combination $S_8={\sigma }_8\sqrt{{\mathrm{\Omega }}_{\mathrm{m}}/0.3}$. To recover the extra information contained in nonlinear scales, many non-Gaussian estimators – also referred as higher-order statistics (HOS) in contrast to two-point statistics – have been introduced in the literature. Such HOS include higher-order moments (e.g., Van Waerbeke et al. 2013; Gatti et al. 2022; Porth & Smith 2021), peak counts (e.g., Marian et al. 2009; Dietrich & Hartlap 2010; Kacprzak et al. 2016; Martinet et al. 2018, 2021b; Harnois-Déraps et al. 2021), one-point probability distributions (e.g., Barthelemy et al. 2020; Boyle et al. 2021; Liu & Madhavacheril 2019; Thiele et al. 2020), Minkowski functionals (e.g., Kratochvil et al. 2012; Petri et al. 2015; Vicinanza et al. 2019; Parroni et al. 2020), Betti numbers (e.g., Feldbrugge et al. 2019; Parroni et al. 2021), persistent homology (e.g., Heydenreich et al. 2021, 2022), scattering transform coefficients (e.g., Cheng et al. 2020; Cheng & Ménard 2021b), as well as map-level inference (Porqueres et al. 2022; Boruah et al. 2022). Despite their increased complexity, which often requires resorting to numerical simulations to model their cosmology dependence, all the references above have demonstrated that these new statistics have superior constraining power compared to the γ-2PCF. However, each of these new HOS is usually developed and studied by independent teams, which renders a fair comparison between them extremely difficult.

The Higher-Order Weak Lensing Statistics (HOWLS) project has been initiated to remedy this situation. One of its main aims is, indeed, to test HOS probes by relying on the same mock data, here mimicking those that Euclid will make available. In contrast to some early (e.g., Pires et al. 2009; Hilbert et al. 2012) and recent efforts in the literature (e.g., Zürcher et al. 2022), HOWLS was designed as a challenge to the community, thus attracting contributions from the largest team of HOS experts ever. Individual teams within the Euclid community have applied 24 different algorithms to the same mocks for a total of two second-order statistics (the shear and convergence two-point correlation functions γ-2PCF and κ-2PCF) and ten different HOS: convergence one-point probability distribution (κ-PDF), higher-order convergence moments (HOM), n-th order aperture mass moments $⟨M_{\mathrm{ap}}^n⟩$, aperture mass peak counts (peaks), convergence Minkowski functionals (MFs), convergence Betti numbers (BNs), aperture mass persistent homology Betti numbers (pers. BNs) and heatmap (pers. heat.), and convergence scattering transform coefficients (ST). Such a large number is unprecedented and offers the possibility of investigating which one (or which combination) is best suited to be coupled with the standard γ-2PCF probe to narrow down the constraints on cosmological parameters (CPs). Different HOS are in fact sensitive to different scales and features in the convergence (κ) maps and thus they couple to the γ-2PCF in their own way. Moreover, HOWLS can also check for correlations among the various HOS probes, revealing which ones are sufficiently uncorrelated such that their combination does indeed improve the total constraining power. It is also worth stressing that the present paper is only the first in a series. HOWLS will actually serve as a preparation for the application of WL HOS to the Euclid Survey, defining common tools and pipelines for the consortium.

The HOWLS data set is based on the DUSTGRAIN-pathfinder simulations (Giocoli et al. 2018a), designed to model the cosmological dependence of every statistic, and on the Scinet LIght-Cones Simulations (SLICS, Harnois-Déraps et al. 2018) for estimating covariances. We have built realistic Euclid mocks out of these simulations, in particular mimicking the expected galaxy density, intrinsic ellipticities, and redshift distribution. For every mock, we built a convergence map following the Kaiser & Squires (1993) implementation described in Pires et al. (2020). We measured the γ-2PCF in the ellipticity catalogs; the κ-2PCF, κ-PDF, MFs, BNs, ST, and HOM from the convergence maps; $⟨M_{\mathrm{ap}}^3⟩$, $⟨M_{\mathrm{ap}}^n⟩$, and pers. BNs and pers. heat. from the aperture mass calculated from the shear field; and peaks of aperture mass maps calculated from the reconstructed convergence fields. We then developed two independent analysis pipelines to compute the Fisher information and thus forecast the constraining power of HOS compared to two-point statistics.

This first paper in the HOWLS series is intended to introduce the data set (Sect. 2) and HOS (Sect. 3), and to conduct a Fisher analysis (Sect. 4) to compare them. Forecasts are presented and discussed in Sect. 5. We conclude in Sect. 6 by listing the refinements that we will include in the following HOWLS publications. The HOWLS data set and applied statistics are publicly released with this article¹.

2. HOWLS data set

To perform a Fisher analysis, one needs to compute data vector (DV) derivatives with respect to individual CPs. We therefore run the DUSTGRAIN-pathfinder simulations, varying one parameter at a time among Ω_m, σ₈, and w, for four different values around the fiducial ones. We additionally used the DUSTGRAIN-pathfinder and SLICS to build the covariance matrix necessary to forecast parameter constraints. The simulations are summarized in Table 1 and described in detail below in Sect. 2.1 for DUSTGRAIN-pathfinder and Sect. 2.2 for SLICS.

2.1. DUSTGRAIN-pathfinder simulations

The DUSTGRAIN-pathfinder suite consists of N-body simulations of volume (750 h⁻¹ Mpc)³ filled with N_p = 768³ particles, corresponding to a particle mass resolution m_p of approximately 8 × 10¹⁰h⁻¹ M_⊙ (Giocoli et al. 2018a). The standard reference CPs have been set to values consistent with the results of the Planck-2015 (Planck Collaboration I 2016) cosmological data analysis, namely a matter density Ω_m of 0.31345, a baryon density Ω_b of 0.0491, Hubble constant H₀ of 67.31 km s⁻¹ Mpc⁻¹, a scalar spectral index n_s of 0.9658, and mean amplitude σ₈ of the linear density fluctuations on the 8 h⁻¹ Mpc scale of 0.842. Since the DUSTGRAIN-pathfinder simulations are part of a cosmological data set that also accounts for modified gravity models, the simulations (including those assuming standard general relativity) were carried out with the MG-Gadget code (Puchwein et al. 2013).

For the analyses performed for this work, in addition to the reference ΛCDM, we used 12 other cosmological runs. In particular, we considered cosmological simulations where only one of the CPs Ω_m, σ₈, and w was varied, either by ±4% and ±16% for σ₈ and w, or by +28%/−36% for Ω_m to allow for existing data to be reused. When varying Ω_m the value of the physical baryon density Ω_bh² was kept fixed in the computation of the linear matter power spectrum adopted in the initial conditions, which was performed by means of the Boltzmann code CAMB (Lewis et al. 2000).

For each cosmological simulation, we built up mass density planes and then shooting-rays in light cones for 128 different line-of-sight realizations (256 in the case of the fiducial cosmology). The past light cones were built using the MapSim routine (Giocoli et al. 2015) following a pyramidal geometry. This method has been used and tested on a variety of cosmological simulations (Tessore et al. 2015; Castro et al. 2018; Giocoli et al. 2018b) and recently compared with other algorithms (Hilbert et al. 2020) finding only percent-level differences for both cosmic shear two-point and peak statistics. The approach we follow in MapSim is based on using several snapshots from a single realization of an N-body simulation to build a light-cone up to a redshift of 4. For the DUSTGRAIN-pathfinder suite there are 21 snapshots available in this redshift range. Given the box length of 750 h⁻¹Mpc, roughly 7 (5) boxes are needed to cover the comoving distance of about 5 (3.6) h⁻¹ Gpc to a source redshift z_s of 4 (2). To obtain better redshift sampling, the volume required to construct the light-cone is divided along the line-of-sight into multiple contiguous redshift slices obtained from the individual snapshots. If the redshift slice reaches beyond the boundary of a single box, two lens planes are constructed from a single snapshot. The total number of lens planes up to z_s = 4 (z_s = 2) is 27 (19). To avoid replicating the same structure along the line of sight, the 7 boxes needed to cover the light-cone are randomized. This randomization procedure allows us to extract multiple realizations from a single simulation. Randomization is achieved by using seeds that act on the simulation boxes based on: (i) changing the location of the observer, typically placed on the center of one of the faces of the box, (ii) redefining the center of the box (taking advantage of periodic boundary conditions), and (iii) changing the signs of the box axes. From each line of sight realization, we then project, using the Born approximation, to construct convergence and shear maps of 5 × 5 deg² at various source redshifts.

In Fig. 1 we exhibit the κ maps for our 13 simulations considering sources at redshift z_s = 2. As can be seen, we display the same line-of-sight realization for the different cosmological simulations: in each subpanel we can recognize by eye the same large-scale structure distribution. More detailed and quantitative statistical information can now be extracted from these maps, for example in Fig. 2 we show the corresponding convergence power spectra, averaged over 128 lines of sight. The left, central and right panels display the Ω_m, σ₈, and w variations, respectively. The black solid line and gray shaded region – the same in all panels – exhibit the average convergence power spectra and the corresponding dispersion for the reference ΛCDM run.

Fig. 1. Convergence maps for sources at zs = 2 of the same light-cone random realization constructed from the DUSTGRAIN-pathfinder simulations used in this work. The region displayed covers a region of approximately 2.5 × 2.5 deg2.

Fig. 2. Convergence power spectra for sources at zs = 2. From left to right, we exhibit the runs that account for Ωm, σ8 and w variations, respectively. The curves show the average over all 128 line of sight realizations, and the gray areas show the scatter of the convergence power spectrum around the mean value for the reference ΛCDM simulation (256 lines of sight) considering a field of view of 5 deg on a side.

2.2. SLICS simulations

The SLICS are a suite of 924 fully independent N-body simulations specifically designed for the estimation of covariance matrices describing WL observables. They were produced by cubep³m, a Poisson solver that computes the nonlinear evolution of 1536³ particles starting from initial conditions created at z = 120 under the Zeldovich approximation (Harnois-Déraps et al. 2013), in boxes of 505 h⁻¹ Mpc on a side. Every run shares the same CPs² but embodies a unique noise realization, thereby providing a large ensemble ideally suited for sample variance estimation. The matter power spectrum agrees to within 2 percent with the COSMIC EMULATOR (Heitmann et al. 2014) up to k = 3.0 h⁻¹ Mpc at z = 0.6 (Harnois-Déraps & van Waerbeke 2015).

As detailed in Harnois-Déraps et al. (2018), the particles were collapsed on-the-fly into mass sheets at 18 predetermined redshifts, from which WL light-cones were constructed up to a redshift of 3.0 following the standard multiple-plane technique. Specifically, convergence and shear maps of 100 deg² were constructed under the Born approximation at 18 source planes, and subsequently sampled to generate Euclid-like galaxy mocks with properties listed in Sect. 2.3, similarly to the methods presented in Sect. 2.1. Each plane has a thickness of 257.5 h⁻¹ Mpc, which, according to Zorrilla Matilla et al. (2020), results in sub-dominant biases on cosmic shear statistics for upcoming surveys.

The SLICS simulations are publicly available³ and were used in a number of cosmic shear data analyses including CFHTLenS (e.g., Joudaki et al. 2017), Kilo Degree Survey (Hildebrandt et al. 2017) and Dark Energy Survey (Harnois-Déraps et al. 2021) data, as well as clustering data analyses including 2dFLenS (Blake et al. 2016), GAMA (van Uitert et al. 2018) and BOSS (Xia et al. 2020). Notably, the covariance matrix estimates of two-point functions have been shown to match well the analytical calculations in Hildebrandt et al. (2017) and Harnois-Déraps et al. (2019), leading to comparable constraints on CPs.

Some dissimilarities between the SLICS and the DUSTGRAIN-pathfinder simulations are worth noting here, as they might have a small but nonnegligible impact on the results presented in this paper. The two suites are based on distinct N-body codes, which results in residual differences in the nonlinear clustering. In addition, we note differences in particle count and mass resolution, in the light-cone opening angles, in pixel sizes⁴ and in the lens randomization procedure. Finally, the shear maps construction pipelines differ in that they are computed directly from the convergence field of view for DUSTGRAIN-pathfinder, but over the full periodic box for SLICS, which eliminates residual edge effects and B-mode leakage. It is worth mentioning that for the DUSTGRAIN-pathfinder runs, we assume void boundary condition, with κ = 0 outside the field of view, when computing the shear field. The only scales that are affected are angular modes l ∼ 10⁵ and change is below one percent. This last caveat has been also highlighted in Hilbert et al. (2020), for a lower resolution map, to have little impact on two-point statistics as also suggested from the good agreement between numerical and theoretical forecasts for γ-2PCF and κ-PDF in the present analysis (see Sect. 4.5).

2.3. Mock galaxy catalogs

We generate mock galaxy catalogs by sampling the shear and convergence planes defined earlier from the DUSTGRAIN-pathfinder and SLICS simulations. This is to achieve a high degree of realism, allowing us to reproduce Euclid survey properties such as the redshift distribution, galaxy density, and shape noise. The procedure closely follows the implementation of Martinet et al. (2021b), but additionally considers magnitudes in the I_E band of the Euclid VIS instrument for sampling the redshift distribution in the case of DUSTGRAIN-pathfinder. In contrast to the mentioned article, we do not present results including tomography in the present work. This is delayed to a future HOWLS paper but our mocks already support any tomographic slicing by construction.

The redshift distributions n(z) of the DUSTGRAIN-pathfinder and SLICS simulations are shown in Fig. 3. They are built from the COSMOS2015 photometric redshift catalogs (Laigle et al. 2016) after a cut in magnitudes I_E ≤ 24.5 in the case of DUSTGRAIN-pathfinder and i′≤24.5 for SLICS. This difference is explained by the fact the SLICS mocks were built earlier for Martinet et al. (2021b) than DUSTGRAIN-pathfinder and did not yet include the information about the Euclid VIS magnitudes (computed as a combination of the r, i′, and z′ magnitudes in the present work). Therefore, the DUSTGRAIN-pathfinder n(z) is closer to the expected Euclid n(z). The difference is, however, less than 5%, except in the redshift range z ≤ 0.2, which contains little information with regard to lensing. Additionally, the model of the dependence on cosmology is built from DUSTGRAIN-pathfinder while SLICS enters in the computation of the covariance matrix. This distinction lowers the risk of any impact of the n(z) difference between the two sets of data on the cosmological inference. After the magnitude cut, the COSMOS2015 n(z) is smoothed by fitting the parametrization from Fu et al. (2008) as it was shown in Martinet et al. (2021b) to capture the high-redshift tail better than the standard Smail et al. (1994) fit. The redshift distributions of the DUSTGRAIN-pathfinder and SLICS mocks are fully characterized by

$$\begin{array}{c}\hfill n(z)=A\frac{z^a+z^{\mathit{ab}}}{z^b+c},\end{array}$$(1)

Fig. 3. Top: Redshift distributions of the mocks built from the DUSTGRAIN-pathfinder (red dots) and SLICS (blue dots) simulations, normalized to 30 galaxies per arcmin2. These correspond to a Fu et al. (2008) fit to the Laigle et al. (2016) COSMOS2015 catalog after removing galaxies with magnitudes IE ≥ 24.5 and i′≥24.5 for DUSTGRAIN-pathfinder and SLICS, respectively (red and blue histograms). Bottom: Fractional difference between the DUSTGRAIN-pathfinder and SLICS redshift distributions. The difference is always below 5% except for redshifts lower than 0.2.

and the parameters listed in Table 2.

Shape noise is included by assigning an intrinsic ellipticity to each galaxy. Specifically, we draw each ellipticity component (ϵ_i, i = {1, 2}) from a Gaussian random distribution centered on 0 and with a dispersion σ_ϵi = 0.26. This reference value (e.g., Euclid Collaboration 2019b) has been measured for a sample of galaxies observed with the Hubble Space Telescope and with similar photometric properties to the expected VIS sample (magnitudes I₈₁₄ ∼ 24.5, Schrabback et al. 2018).

The impact of shape noise on the cosmological model is minimized by using the same random realization of galaxy intrinsic ellipticities and positions across all cosmologies for a given mock. This means that we have 128 independent realizations of shape noise for the DUSTGRAIN-pathfinder mocks, but these are identical for the 12 cosmologies probed. Conversely, the positions and ellipticities are fully random for any realization used in the covariance matrix computation, either with the DUSTGRAIN-pathfinder or SLICS simulations. This ensures the shape noise contribution to the error budget is faithfully captured. This process of fixing shape noise in the model and leaving it free in the covariance has become standard practice for simulation-based inference with higher-order mass map estimators (e.g., Kacprzak et al. 2016; Martinet et al. 2018; Harnois-Déraps et al. 2021).

These galaxy catalogs are then fed to the Euclid convergence map reconstruction pipeline described in Pires et al. (2020) and in Sect. 2.4 to produce the convergence maps on which most HOS will be measured. These catalogs are also used to compute direct statistics from the shear, specifically γ-2PCFs (Sect. 3.1) and aperture masses (M_ap, see Sects. 3.4 and 3.8).

2.4. Mass mapping

The statistical properties of the WL field can be assessed by a statistical analysis either of the shear field or of the convergence field. Many HOS are traditionally computed from the κ field. This requires solving a mass inversion problem that consists of reconstructing the convergence κ from the measured shear field γ. Using complex notation, the shear field is written as γ = γ₁ + iγ₂, and the convergence field as κ = κ_E + iκ_B, with κ_E and κ_B respectively corresponding to the E- and B-mode components of the field, by analogy with the electromagnetic field.

We can derive the relation between the shear field γ and the convergence field κ in the Fourier domain (Kaiser & Squires 1993) with

$$\begin{array}{c}\hfill \widehat{\kappa }={\widehat{P}}^{\ast }\widehat{\gamma },\end{array}$$(2)

where the hat symbol denotes Fourier transforms, ${\widehat{P}}^{\ast }$ is the complex conjugate, and $\widehat{P}={\widehat{P}}_1+\mathrm{i}{\widehat{P}}_2$ with

$$\begin{array}{cc}& \widehat{P_1}(\mathit{\ell })=\frac{{\ell }_1^2-{\ell }_2^2}{{\ell }^2},\hfill \\ \hfill & \widehat{P_2}(\mathit{\ell })=\frac{2{\ell }_1{\ell }_2}{{\ell }^2},\hfill \end{array}$$(3)

with ${\mathcal{l}}^2\equiv {\mathcal{l}}_1^2+{\mathcal{l}}_2^2$ and ℓ_i the Fourier counterparts of the angular coordinates θ_i. The convergence can only be determined up to an additive constant because there is a degeneracy when ℓ₁ = ℓ₂ = 0 (see e.g., Bartelmann 1995). In practice, we impose that the mean convergence vanishes across the field by setting the reconstructed ℓ = 0 mode to zero.

Assuming the mass inversion is conducted without noise regularization, the same information is contained in the shear field as in the convergence maps. However, it is well known that the Kaiser-Squires inversion creates undesirable artifacts at the borders of the reconstructed convergence maps (see e.g., Seitz & Schneider 1996, 2001; Pires et al. 2020). This is due to the fact that the discrete Fourier transform implicitly assumes periodicity of the image along both dimensions. In a future work, mass mapping systematic effects will be further mitigated and their impact quantified by propagating the errors into CP forecasts using HOS. In addition to border effects and masks, we will test the impact of reduced shear. In this article, however, we assume that the mean ellipticity is an unbiased estimator of the mean shear instead of reduced shear such that we can replace γ by ϵ in Eq. (2), an assumption only correct in the weak regime.

The shear is sampled only at the positions of the galaxies. Therefore, the first step of the mass inversion method is to bin the observed ellipticities of galaxies on a regular pixel grid to create the shear maps. In practice, we bin the galaxies in pixels of size of 0$\stackrel{\prime }{.}$59, resulting in shear maps of 512 × 512 pixels (for DUSTGRAIN-pathfinder) and 1024 × 1024 pixels (for SLICS) that are then converted into convergence maps using Eq. (2). The statistical analysis is performed only on the E-modes convergence maps because WL only produces E-modes. However, the application of the HOS to the B-modes map can be used to test for residual systematic effects.

3. Statistics

Many higher-order probes have been proposed, tested, and measured on present-day Stage III lensing maps with promising preliminary results (e.g., Martinet et al. 2018; Gatti et al. 2022; Harnois-Déraps et al. 2021; Heydenreich et al. 2022; Zürcher et al. 2022; Burger et al. 2023). This consideration motivated us to focus our analysis on HOS of scalar fields derived from the shear: convergence and aperture mass. In the following paragraphs, we review the probes we consider, referring the interested reader to the quoted papers for further details. Far from being fully exhaustive, this short review aims at presenting the tools we have used and giving the reader an overview of the many roads that open up when going beyond second-order statistics. A short description of theoretical predictions is given for the 2PCF and the κ-PDF as they are used to validate the simulated derivatives entering the Fisher analysis, see Sect. 4.5. The impatient reader can directly look at Table 3, where we list the HOS we have used together with their abbreviations, the number of independent teams that applied them, and the subsections in which they are described. The fact that each statistic was computed by a different team led to a variety of choices in terms of the filtering of the shear or convergence field. For this publication we do not try to homogenize these choices as we consider it part of each method and list these differences in Table 3.

3.1. Two-point correlation functions

Although useful on their own, HOS are at their best when used in combination with standard second-order statistics, coupling the typically larger signal-to-noise ratio (S/N) of lower-order statistics with the degeneracy-breaking power of HOS. In the following, we therefore quantify both the constraints from each HOS probe alone and the improvement of the constraints from joint second- and higher-order statistics with respect to the second-order-only case.

The most basic cosmic shear observable is the real-space shear two-point correlation functions (γ-2PCF), since it can be estimated by simply multiplying the ellipticities of galaxy pairs and averaging. The shear can conveniently be decomposed into a tangential, γ_t, and cross – component, γ_×, such that

$$\begin{array}{c}\hfill {\gamma }_{\mathrm{t}}({\mathit{\vartheta }}^{\prime },\mathit{\vartheta })=-\mathfrak{R}[\gamma ({\mathit{\vartheta }}^{\prime }){\mathrm{e}}^{-2\mathrm{i}\varphi }]{\gamma }_{\times }({\mathit{\vartheta }}^{\prime },\mathit{\vartheta })=-\mathfrak{I}[\gamma ({\mathit{\vartheta }}^{\prime }){\mathrm{e}}^{-2\mathrm{i}\varphi }],\end{array}$$(4)

with ℜ and ℑ the real and imaginary parts, ϕ the polar angle of the direction vector between the galaxy position ϑ′ and a reference point ϑ, and the minus sign a convention to have the tangential shear positive around a mass overdensity. The two shear components can then be combined to get two 2PCFs (see, e.g., Kilbinger 2015, and references therein)

$$\begin{array}{c}\hfill \{\begin{array}{c}{\displaystyle {\xi }_{+}(\theta )=⟨\gamma {\gamma }^{\ast }⟩(\theta )=⟨{\gamma }_{\mathrm{t}}{\gamma }_{\mathrm{t}}⟩(\theta )+⟨{\gamma }_{\times }{\gamma }_{\times }⟩(\theta )}\hfill \\ \hfill \\ {\displaystyle {\xi }_{-}(\theta )=\mathfrak{R}[⟨\gamma \gamma ⟩(\theta )e^{-4\mathrm{i}\varphi }]=⟨{\gamma }_{\mathrm{t}}{\gamma }_{\mathrm{t}}⟩(\theta )-⟨{\gamma }_{\times }{\gamma }_{\times }⟩(\theta )}\hfill \end{array},\end{array}$$(5)

where the dependence is only on the angular separation θ on the sky because under the Cosmological Principle cosmic fields are statistically invariant under translation and rotation.

The main virtue of these γ-2PCF is that they can be straightforwardly estimated from the measured ellipticities ε_i as

$$\begin{array}{c}\hfill {\xi }_{\pm }(\theta )=\frac{{\sum }_{i,j}{w}_i{w}_j({\epsilon }_{\mathrm{t},i}{\epsilon }_{\mathrm{t},j}\pm {\epsilon }_{\times ,i}{\epsilon }_{\times ,j})}{{\sum }_{i,j}{w}_i{w}_j},\end{array}$$(6)

where the sum runs over the galaxy pairs with positions on the sky (θ_i, θ_j) having angular separation |θ_i − θ_j| in a bin centered on θ. The weight w_i of the ellipticity ε_i can be related to the measurement error, and set to zero if the galaxy is in a masked region. In the present article the lensing weights are set to 1, as the impact of masks and shear measurement methods are beyond the scope of this analyis.

The γ-2PCF can be easily computed for a given cosmological model by first going to Fourier space and then converting back into real space. The final result is

$$\begin{array}{c}\hfill \{\begin{array}{c}{\displaystyle {\xi }_{+}(\theta )=\frac{1}{2\pi }{\displaystyle \int }[P_{\kappa }^E(\ell )+P_{\kappa }^B(\ell )]J_0(\ell \theta )\ell \mathrm{d}\ell }\hfill \\ {\displaystyle {\xi }_{-}(\theta )=\frac{1}{2\pi }{\displaystyle \int }[P_{\kappa }^E(\ell )-P_{\kappa }^B(\ell )]J_4(\ell \theta )\ell \mathrm{d}\ell }\hfill \end{array},\end{array}$$(7)

where $P_{\kappa }^E(\ell )$ and $P_{\kappa }^B(\ell )$ are the power spectra of the convergence E and B modes, while J_n(x) is the n-th order spherical Bessel function of the first kind. We note that, in the absence of systematics and neglecting higher-order effects, WL does not produce B modes so that $P_{\kappa }^B(\ell )=0$, and ξ_± reduce to different Hankel transforms of the same quantity, which can be derived from the matter power spectrum P_δ(ℓ).

As done for the shear, we can similarly define the κ-2PCF for the convergence. In particular, this will inherit the same properties of translation and rotation invariance as in the shear case. The main difference is that, with the convergence being a scalar quantity, there is only one single correlation function. In Fourier space, this can be related to the convergence power spectrum as

$$\begin{array}{c}\hfill ⟨\widehat{\kappa }(\mathit{\ell }){\widehat{\kappa }}^{\ast }(\mathit{\ell }\prime )⟩={(2\pi )}^2{\delta }_{\mathrm{D}}(\mathit{\ell }-{\mathit{\ell }}^{\prime })P_{\kappa }(\ell ),\end{array}$$(8)

where δ_D is the Dirac-δ function and the convergence power spectrum depends only on the modulus of ℓ due to the statistical homogeneity and isotropy. It is then possible to show that $P_{\kappa }(\ell )=P_{\kappa }^E(\ell )$, and the κ-2PCF is given by

$$\begin{array}{c}\hfill {\xi }_{\kappa }(\theta )=\frac{1}{2\pi }{\displaystyle \int }{P}_{\kappa }(\ell )J_0(\ell \theta )\ell \mathrm{d}\ell ,\end{array}$$(9)

which is the same as ξ₊(θ). In the following, we keep this different label in order to distinguish between the γ-2PCF and κ-2PCF.

In this article, the γ-2PCF has been computed with the ATHENA software (Kilbinger et al. 2014) using 10 logarithmic bins between 0$\stackrel{\prime }{.}$1 and 300′. The κ-2PCF has been measured using the public code TreeCorr (Jarvis et al. 2004). The minimum separation considered is 0$\stackrel{\prime }{.}$59 (corresponding to the pixel scale) and the maximum separation is approximately 424′ (corresponding to the DUSTGRAIN-pathfinder map diagonal), with 25 bins. As discussed in Sect. 4 several bins are later removed to pass our quality criteria, with the final scale range described in Table 3. We also show these DVs as well as their standard Fisher derivatives with respect to Ω_m, σ₈, and w in Figs. 4 and 5 for γ-2PCF and κ-2PCF respectively. Finally, we note some difference in the chosen binning by the different teams computing the 2PCFs which leads to some artificial difference between ξ₊ and ξ_κ, for example.

Fig. 4. γ-2PCF ξ+ (top) and ξ− (bottom) DVs and derivatives. The gray lines, whose scale is given by the second axis, correspond to the average DVs computed from the 924 SLICS realizations. The Fisher derivatives (defined in Eq. (47)) are computed from the DUSTGRAIN-pathfinder simulations with large variations of Ωm (orange), σ8 (green), and w (blue). The solid, dashed, and dotted lines respectively correspond to the average over 128, 64, and 32 realizations. The shaded areas represent the uncertainty computed from the 128 DUSTGRAIN-pathfinder realizations, and the gray error bars those of the 924 SLICS realizations. The inclusion of the dashed and dotted lines within the shaded areas highlights the low numerical noise.

Fig. 5. Fiducial DV for κ-2PCF (gray) along with its derivatives with respect to Ωm (orange), σ8 (green), and w (blue).

3.2. Convergence PDF

The one-point probability distribution function (PDF) of the WL convergence κ encodes vital information about the non-Gaussian late-time density field. In the presence of shape noise, which dominates on small scales, this information can be extracted most conveniently from the PDF of the convergence after smoothing on an angular scale big enough that the variance set by the gravitational clustering is larger than the shape noise contribution. As a one-point statistic, the PDF is straightforward to measure from simulated and real data, as evidenced by recent analyses of density-split statistics and weak lensing moments in the Dark Energy Survey (Friedrich et al. 2018; Gruen et al. 2018; Gatti et al. 2020, 2022), which carry PDF information in a compressed way and can deal with survey mask effects.

When focusing on mildly nonlinear scales, typically corresponding to smoothing scales of about 10′ at low redshifts, this one-point PDF also admits accurate theoretical predictions (see for example Barthelemy et al. 2020; Boyle et al. 2021). The main idea of this theoretical model relies on the fact that, within the Limber approximation, since the WL convergence probes the matter density in a cone whose opening angle is set by the smoothing scale, the κ-PDF can be predicted using cylindrical collapse applied to the 2D slices of the density in the circular cross-sections. Hence, for a set of sources regrouped in redshift bins and with a certain source distribution across the survey n(z), we recall that the theoretical cumulant generating function (CGF) ϕ_κ, θ at a given angular smoothing scale θ can be expressed as

$$\begin{array}{c}\hfill {\varphi }_{\kappa ,\theta }(\lambda )={\displaystyle \int }\frac{\mathrm{d}zc}{H(z)}{\varphi }_{\mathrm{cyl}}[{\omega }_{n(z)}(z)\lambda ,z],\end{array}$$(10)

with ϕ_cyl the CGF of the density in each individual 2D slice of radius χ(z)θ, and ω_n(z) the generalized lensing kernel given a wide distribution of sources following the normalized distribution n(z)

$$\begin{array}{c}\hfill {\omega }_{n(z)}(z)=\frac{3{\mathrm{\Omega }}_{\mathrm{m}}{H}_0^2}{2c^2}{\displaystyle \int }\mathrm{d}{z}_{\mathrm{s}}n(z_{\mathrm{s}})\frac{[\chi (z_{\mathrm{s}})-\chi (z)]\chi (z)}{\chi (z_{\mathrm{s}})}\mathcal{H}(z_{\mathrm{s}}-z)(1+z),\end{array}$$(11)

where the Heaviside ℋ ensures that the integrand vanishes for z ≥ z_s. ϕ_cyl is obtained from the cylindrical collapse dynamics as a proxy for the whole nonlinear evolution in cylindrically symmetric configurations (top-hat smoothing) and enforces a specific hierarchy of cumulants. The PDF is then recovered through the inverse Laplace transform of the exponential of the κ CGF as

$$\begin{array}{c}\hfill {\mathcal{P}}_{\theta }(\kappa )={\displaystyle {\int }_{-i\infty }^{+i\infty }}\frac{\mathrm{d}\lambda }{2\pi \mathrm{i}}exp[-\lambda \kappa +{\varphi }_{\kappa ,\theta }(\lambda )].\end{array}$$(12)

Note that the cylindrical collapse model is accurate enough to predict the so-called reduced cumulants of the (2D) density field $(S_n={⟨{\delta }^n⟩}_{\mathrm{c}}/{\sigma }_{\delta }^{2(n-1)})$, where ⟨δⁿ⟩_c are the cumulants of the matter density. This means that the model also requires an external input for the prediction of the nonlinear variance of all the 2D slices of the density along the line of sight. Fortunately, the emulation of the matter power spectrum has received lots of attention these past years and can be estimated rather accurately, for example by using the revised Halofit model (Takahashi et al. 2012) or the Euclid Emulator (Euclid Collaboration 2019a). Galaxy shape noise can be included in the theoretical prediction through a convolution of the noiseless PDF as described in Boyle et al. (2021) or directly at the level of the CGF by simple addition of – for example – an associated shape noise variance ${\sigma }_{\text{SN}}^2{\lambda }^2/2$.

We extracted the κ-PDF from the simulated maps for two top-hat filters with smoothing scales θ_s ∈ {4$\stackrel{\prime }{.}$69, 9$\stackrel{\prime }{.}$37}, corresponding to 8/16 pixels in the DUSTGRAIN-pathfinder simulations. After smoothing the maps with the appropriate top-hat filter, we excluded all pixels whose smoothing circle of radius θ_s would intersect the patch boundary.⁵ The PDF was obtained from a histogram using 201 linearly spaced bins in the range κ ∈ [ − 0.1, 0.1]. Cuts were made in the considered range of κ to exclude the extreme tails and keep the DV close to Gaussian as discussed in Sect. 4.2. The precise κ values are given in Fig. 6. The cuts exclude 0.05% and 0.25% of the cumulative probability from the low-κ and high-κ tails, respectively, for the smaller smoothing scale. We subsequently compare the theoretical predictions to the measurements in the simulation in Sect. 4.5.

Fig. 6. Fiducial DV for κ-PDF (gray) along with its Fisher derivatives with respect to Ωm (orange), σ8 (green), and w (blue). The κ-PDF was computed on κ maps smoothed with a top-hat filter with a scale of 4$\stackrel{\prime }{.}$69.

One alternative to the κ-PDF is to turn to the (slightly) more involved – but also more easily accessible through observations – modeling of the one-point PDF of the aperture mass from Eq. (13) discussed in Sect. 3.4. The analogous theoretical model was developed in Barthelemy et al. (2021). The next extension would be to include some tomographic information in the analysis. At the level of the PDF one could rely on the ideas developed in Barthelemy et al. (2022) for the joint PDF between – built to be – independent lensing kernels that probe structures along the line of sight, but the generalization to any other lensing observables should be straightforward.

3.3. Higher-order convergence moments

The shear and convergence 2PCF and their harmonic space counterparts, the power spectra, are related to the variance of the lensing fields. It is then natural to ask whether additional information is encoded in moments higher than two, given that they are indicators of the non-Gaussianity of the field. We consider here the second, third, and fourth moments of the convergence field.

Moments of the smoothed WL convergence field κ_θ can be calculated from weighted averages of the lensing PDF as $⟨{\kappa }^n⟩(\theta )=\int \mathrm{d}\kappa {\kappa }^n{\mathcal{P}}_{\theta }(\kappa )$, assuming a lensing convergence field of zero mean. The variance ⟨κ²⟩=σ² fully characterizes a Gaussian one-point distribution, while the skewness ⟨κ³⟩/σ³ and kurtosis ⟨κ⁴⟩/σ⁴ − 3 are the most common examples encoding non-Gaussian information. While an analysis of the variance, skewness, and kurtosis promises a significant compression of the κ-PDF DV, higher-order cumulants are known to be very sensitive to the tails of the distribution, rendering their signal and likelihood hard to predict.

An approximated yet accurate formula for the moments can be obtained under the hierarchical ansatz (Bernardeau & Schaeffer 1992; Szapudi & Szalay 1993), using this ansatz to express the bispectrum and trispectrum in terms of the power spectrum and plugging them into the general expression for the convergence moments (see Munshi & Jain 2001; Vicinanza et al. 2019, for details).

It is worth noting that moments are a function of smoothing radius, so one typically considers as a cosmological probe their scaling with θ. In contrast, here we use them for a fixed θ since, in a preliminary analysis, we have found moments obtained from different smoothing radii to be highly correlated. Although some information is still present in the correlated moments, we prefer here to first focus on the choice of the filter radius in our two parameter analysis for which combining moments at different radii is not essential. We, therefore, measure them like the κ-PDF using a top-hat filter with two different values of the aperture radius, namely θ ∈ {4$\stackrel{\prime }{.}$69, 9$\stackrel{\prime }{.}$37}, and take {⟨κ²⟩(θ),⟨κ³⟩(θ),⟨κ⁴⟩(θ)} for a fixed θ as our DV (see Fig. 7). To this end, we first define

$$\begin{array}{c}\hfill {\kappa }_{\theta }(\mathit{\vartheta })={\displaystyle \int }{\mathrm{d}}^2\vartheta \prime {\mathcal{H}}_{\theta }({\mathit{\vartheta }}^{\prime }-\mathit{\vartheta })\kappa ({\mathit{\vartheta }}^{\prime })\end{array}$$

Fig. 7. Fiducial DV (gray) along with its derivatives with respect to Ωm (orange), σ8 (green), and w (blue) for HOM computed on κ maps smoothed with a top-hat filter with a scale of 4$\stackrel{\prime }{.}$69.

as the smoothed convergence field. We then take it to the second, third, fourth power and compute its average value over the pixels to compute the data vector for every single map. The mean over the different realizations of the simulations makes up our final DV.

3.4. Higher-order aperture mass moments

Higher-order moments as introduced before were extracted for the smoothed convergence maps. However, it is also possible to compute them for the aperture mass M_ap map, at sky location ϑ and scale θ given by (Schneider et al. 1998):

$$\begin{array}{c}\hfill M_{\mathrm{ap}}(\mathit{\vartheta },\theta )={\displaystyle \int }{\mathrm{d}}^2{\vartheta }^{\prime }{U}_{\theta }(|\mathit{\vartheta }-{\mathit{\vartheta }}^{\prime }|)\kappa ({\mathit{\vartheta }}^{\prime }),\end{array}$$(13)

where U_θ is a compensated filter function with

$$\begin{array}{c}\hfill {\displaystyle \int }{\mathrm{d}}^2\vartheta U_{\theta }(|\mathit{\vartheta }|)=0.\end{array}$$(14)

The higher-order aperture maps statistics $⟨M_{\mathrm{ap}}^n⟩$ correlate the aperture masses on n different scales θ₁, ⋯, θ_n, and are defined as

$$\begin{array}{c}\hfill ⟨M_{\mathrm{ap}}^n⟩({\theta }_1,\cdots ,{\theta }_n)=⟨M_{\mathrm{ap}}(\mathit{\vartheta },{\theta }_1)\cdots M_{\mathrm{ap}}(\mathit{\vartheta },{\theta }_n)⟩.\end{array}$$(15)

The $⟨M_{\mathrm{ap}}^3⟩$ is closely related to the third-order moment of the convergence maps since both of them are sensitive to the matter bispectrum (Schneider et al. 2005). However, in contrast to other HOS considered in this work, the $⟨M_{\mathrm{ap}}^3⟩$ can be directly inferred from the shear and are therefore not affected by systematics induced by the convergence reconstruction. Additionally $⟨M_{\mathrm{ap}}^3⟩$ can be inferred from the shear three-point correlation functions Γ₀ = ⟨γ γ γ⟩ and Γ₁ = ⟨γ^* γ γ⟩, which can be easily estimated even for irregular survey geometries and in the presence of masks (Schneider & Lombardi 2003; Heydenreich et al. 2023). Consequently, $⟨M_{\mathrm{ap}}^3⟩$ is more straightforward to measure in a realistic survey setting than HOS of κ.

We infer the $⟨M_{\mathrm{ap}}^3⟩$ from the tangential shear γ_t(ϑ′;ϑ) at position ϑ′ with respect to ϑ. This inference is possible because each filter U_θ is associated with a function Q_θ for which

$$\begin{array}{c}\hfill M_{\mathrm{ap}}(\mathit{\vartheta },\theta )={\displaystyle \int }{\mathrm{d}}^2{\vartheta }^{\prime }{Q}_{\theta }(|\mathit{\vartheta }-{\mathit{\vartheta }}^{\prime }|){\gamma }_{\mathrm{t}}({\mathit{\vartheta }}^{\prime };\mathit{\vartheta }).\end{array}$$(16)

The function Q_θ is given by

$$\begin{array}{c}\hfill Q_{\theta }(\vartheta )=\frac{2}{{\vartheta }^2}{\displaystyle {\int }_0^{\vartheta }}\mathrm{d}{\vartheta }^{\prime }{\vartheta }^{\prime }{U}_{\theta }({\vartheta }^{\prime })-U_{\theta }(\vartheta ).\end{array}$$(17)

To mimic the application of the $⟨M_{\mathrm{ap}}^3⟩$ statistics to a survey, we use Eq. (16) and work directly with the simulated shear catalogs. We use the filter function U_θ from Crittenden et al. (2002),

$$\begin{array}{c}\hfill U_{\theta }(\vartheta )=\frac{1}{2\pi {\theta }^2}(1-\frac{{\vartheta }^2}{2{\theta }^2})exp(-\frac{{\vartheta }^2}{2{\theta }^2}),\end{array}$$(18)

and the scale radii θ ∈ {1$\stackrel{\prime }{.}$17, 2$\stackrel{\prime }{.}$34, 4$\stackrel{\prime }{.}$69, 9$\stackrel{\prime }{.}$37}.

Our measurement of $⟨M_{\mathrm{ap}}^3⟩$ proceeds in two steps following the procedure in Sect. 5.3.1 of Heydenreich et al. (2023). First, we measure M_ap(ϑ, θ). For this, we employ the convolution theorem to solve the convolution of Q_θ and γ_t in Fourier space. Since the tangential shear γ_t can be written as

$$\begin{array}{c}\hfill {\gamma }_{\mathrm{t}}({\mathit{\vartheta }}^{\prime },\mathit{\vartheta })=-\mathfrak{R}\left[\gamma ({\mathit{\vartheta }}^{\prime })\frac{{(\mathit{\vartheta }-{\mathit{\vartheta }}^{\prime })}^{\ast }}{\mathit{\vartheta }-{\mathit{\vartheta }}^{\prime }}\right],\end{array}$$(19)

where the vectors are interpreted as complex numbers ϑ = ϑ₁ + iϑ₂, Eq. (16) transforms into

$$\begin{array}{c}\hfill M_{\mathrm{ap}}(\mathit{\vartheta },\theta )=-{\displaystyle \int }{\mathrm{d}}^2{\vartheta }^{\prime }{Q}_{\theta }(|\mathit{\vartheta }-{\mathit{\vartheta }}^{\prime }|)\frac{{(\mathit{\vartheta }-{\mathit{\vartheta }}^{\prime })}^{\ast }}{\mathit{\vartheta }-{\mathit{\vartheta }}^{\prime }}\gamma ({\mathit{\vartheta }}^{\prime }).\end{array}$$(20)

Therefore, we can calculate both $Q_{\theta }(|\mathit{\vartheta }|)\frac{{\mathit{\vartheta }}^{\ast }}{\mathit{\vartheta }}$ and γ(ϑ) on a grid, and solve the convolution in Eq. (20) using the Fast Fourier Transform (FFT). To avoid border effects, we cut off a strip of width of 4θ from each border. This large cut-off is needed because the exponential aperture filter is not exactly zero for ϑ > θ, so at distance θ, we still experience border effects from each side of the M_ap(ϑ, θ) map. Our cut-off means we neglect 0.07% of the total filter power, which we deem acceptable. Second, we measure $⟨M_{\mathrm{ap}}^3⟩$(θ₁, θ₂, θ₃). For this, we multiply M_ap(ϑ, θ₁), M_ap(ϑ, θ₂), and M_ap(ϑ, θ₃) for each line-of-sight and each position ϑ. Then, we average over ϑ, which gives $⟨M_{\mathrm{ap}}^3⟩$(θ₁, θ₂, θ₃) for each line-of-sight. This DV is shown in Fig. 8.

Fig. 8. Fiducial DV (gray) along with its derivatives with respect to Ωm (orange), σ8 (green), and w (blue) for $⟨M_{\mathrm{ap}}^3⟩$ of aperture mass maps computed from shear fields smoothed with a Crittenden et al. (2002) filter with scales (θ1, θ2, θ3, θ4) = (1$\stackrel{\prime }{.}$17, 2$\stackrel{\prime }{.}$34, 4$\stackrel{\prime }{.}$69, 9$\stackrel{\prime }{.}$37).

We estimate the n-th order moments $⟨M_{\mathrm{ap}}^n⟩$ (Fig. 9) with a different approach. For this we cover the survey footprint with apertures and estimate $⟨M_{\mathrm{ap}}^n⟩$ within each aperture (Schneider et al. 1998) as

$$\begin{array}{cc}\hfill \widehat{M_{\mathrm{ap}}^n}(\mathit{\vartheta };{\theta }_1,\cdots ,{\theta }_n)& =\left(\pi {\theta }_1^2\right)\cdots \left(\pi {\theta }_n^2\right)\hfill \\ \hfill & \times \frac{{\sum }_{i_1\ne i_2\ne \cdots \ne i_n}{w}_{i_1}{Q}_{{\theta }_1;i_1}{\epsilon }_{\mathrm{t},i_1}\cdots w_{i_n}{Q}_{{\theta }_n;i_n}{\epsilon }_{\mathrm{t},i_n}}{{\sum }_{i_1\ne i_2\ne \cdots \ne i_n}{w}_{i_1}\cdots w_{i_n}},\hfill \end{array}$$(21)

Fig. 9. Fiducial DV (gray) along with its derivatives with respect to Ωm (orange), σ8 (green), and w (blue) for $⟨M_{\mathrm{ap}}^n⟩$ of aperture mass maps computed from shear fields smoothed with a polynomial filter with scales θi, j, k ∈ {1$\stackrel{\prime }{.}$17, 2$\stackrel{\prime }{.}$34, 4$\stackrel{\prime }{.}$69, 9$\stackrel{\prime }{.}$37}. We first display second, fourth, and fifth identical scale moments, followed by the third order moments ordered identically to Fig. 8.

where each w_i is the weight of the ellipticity ε_i and we abbreviated Q_θj(|ϑ_ij − ϑ|) ≡ Q_θj; ij for an aperture centered at ϑ. Note that the index i_k in the kth sum in (21) only considers galaxies that lie within the support of Q_θk. By averaging over all apertures in the footprint with appropriate weights w_ap, one then obtains an unbiased estimate for $⟨M_{\mathrm{ap}}^n⟩$. As shown in Porth et al. (2020) and Porth & Smith (2021), one can decompose the nested sums appearing in (21), so that the full estimation process scales linearly with the number of galaxies. In this work, we estimate the connected parts of the aperture-mass statistics for n ∈ {2, 3, 4, 5}, where for n = 3 we take into account all of the different combinations of aperture radii and in the other cases only consider the components for which θ₁ = … = θ_n. In contrast to the FFT-based estimation procedure described above, we employ a polynomial filter function introduced in Schneider et al. (1998) for the direct estimator method,

$$\begin{array}{c}\hfill Q_{\theta }(\vartheta )=\frac{6}{\pi {\theta }^2}{\left(\frac{\vartheta }{\theta }\right)}^2[1-{\left(\frac{\vartheta }{\theta }\right)}^2]\mathcal{H}(\theta -\vartheta ),\end{array}$$(22)

where ℋ(x) denotes the Heaviside function. For our choice of w_ap, we follow Porth & Smith (2021) who propose a form that approximates an inverse variance weighting scheme. As the mocks used in this work have w_i ≡ 1, it can be shown that in this case the weight of an aperture centered at ϑ is equal to the number of multiplet counts within its configuration of aperture radii,

$$\begin{array}{c}\hfill w_{\mathrm{ap}}(\mathit{\vartheta };{\theta }_1,\cdots ,{\theta }_n)={\displaystyle \sum _{i_1\ne i_2\ne \cdots \ne i_n}}1,\end{array}$$(23)

where the indices are bound to the same constraints as in (21). Analogous to the FFT-based estimation procedure we need to avoid border effects; with our choice (22) for the filter function this results in cutting off a strip of width of max({θ₁,⋯,θ_n}) from each border. In the resulting DV (Fig. 9) the first twelve elements correspond to the equal scale statistics of order {2, 4, 5} while the final 20 elements correspond to the third order statistics and are ordered as the DV in Fig. 8.

3.5. Aperture mass peak counts

A different way to probe the convergence field at all statistical orders in a single step is to consider peak counts. As the name itself implies, one is now searching for peaks (i.e., local over-densities) on the smoothed convergence maps.

Some studies only focus on peaks with very large S/N because they are good tracers of massive galaxy clusters (Kruse & Schneider 1999; Marian & Bernstein 2006; Gavazzi & Soucail 2006; Hamana et al. 2015; Miyazaki et al. 2017). The advantage of this approach is that the dependence on cosmology of the WL cluster abundance can be accurately predicted by theoretical models (see e.g., Kruse & Schneider 2000; Bartelmann etal. 2001; Hamana et al. 2004; Marian et al. 2009). Another approach consists of also considering low-amplitude peaks. Since there is no analytical prediction for the full range of S/N peaks, it is necessary to run a large number of N-body simulations to calibrate the dependence of the peak count statistics on cosmology. However, a significant fraction of these peaks arise from large-scale structure projections, and as such, carry additional cosmological information (see e.g., Yang et al. 2013; Lin et al. 2016). Over the past two decades, many cosmological studies have been performed based on the second approach (see e.g., Pires et al. 2009; Dietrich & Hartlap 2010; Martinet et al. 2018, 2021a,b; Peel et al. 2018; Li et al. 2019; Ajani et al. 2020; Harnois-Déraps et al. 2021) and have shown the strength of peak counts in discriminating between cosmological models. In a preliminary Fisher analysis, we tested the two approaches and, as expected, the best results were obtained when including the full range of S/N peaks, which is therefore what we present in the following.

The computation of the peak count statistics requires filtering the convergence maps because of galaxy shape noise. In practice, the peak count can be evaluated at a given scale by convolving the convergence maps with a filter function of a specific scale (i.e., aperture radius). We choose to filter the maps using compensated filters defined in Eq. (14). Compensated filters (e.g., aperture mass filters or wavelet filters) have been preferred over low-pass filters (e.g., Gaussian filters) due to their shapes, which reduce the overlap between different scales and then the correlations between them (see e.g., Lin et al. 2016; Ajani et al. 2020). In Leonard et al. (2012), it is demonstrated that the aperture mass is formally identical to a wavelet transform at a specific scale.

As such, we use the starlet transform (Starck et al. 2007; Leonard et al. 2012) to simultaneously compute five aperture mass maps corresponding to scales of {1$\stackrel{\prime }{.}$17, 2$\stackrel{\prime }{.}$34, 4$\stackrel{\prime }{.}$69, 9$\stackrel{\prime }{.}$37, 18$\stackrel{\prime }{.}$74}. The starlet transform decomposes the convergence as follows:

$$\begin{array}{c}\hfill \kappa (\mathit{\vartheta })=C_J(\mathit{\vartheta })+{\displaystyle \sum _{i=1}^J}{W}_{\theta }^i(\mathit{\vartheta }),\end{array}$$(24)

where C_J is a smooth version of the convergence κ and $W_{\theta }^i$ are the wavelet maps corresponding to a scale of θ = 2ⁱ pixels.

The starlet transform is equivalent to applying the following aperture mass filter to the convergence map (Leonard et al. 2012):

$$\begin{array}{c}\hfill U_{\theta }(\vartheta )=\frac{31}{3}{\left|\frac{\vartheta }{\theta }\right|}^3-\frac{64}{9}({|\frac{1}{2}-\frac{\vartheta }{\theta }|}^3+{(\frac{1}{2}+\frac{\vartheta }{\theta })}^3)\\ \hfill +2({|1-\frac{\vartheta }{\theta }|}^3+{|1+\frac{\vartheta }{\theta }|}^3)-\frac{1}{18}({|2-\frac{\vartheta }{\theta }|}^3+{|2+\frac{\vartheta }{\theta }|}^3).\end{array}$$(25)

In the aperture mass maps, the peaks are identified as individual pixels higher than their eight neighbors. The edges are discarded from the computation. Once the peaks are detected on each aperture mass map, they are classified depending on their amplitudes with respect to the shape noise. Several implementations of the peaks have been tested. The main differences between the implementations are in the range of the amplitudes of the peaks that are considered, the number of equidistant bins and the size of the edges to be discarded.

The results for the peak counts presented in this study are obtained by sorting the peaks into 14 equidistant bins between −1 and 6 times the dispersion value (see Table 4 and Fig. 10) to ensure it is Gaussian distributed for the Fisher analysis and by discarding a stripe of 1 pixel on each side of the M_ap map before counting peaks.

Fig. 10. Fiducial DV (gray) along with its derivatives with respect to Ωm (orange), σ8 (green), and w (blue) for peaks of aperture mass maps, computed from κ-fields smoothed with a starlet filter with a scale of 2$\stackrel{\prime }{.}$34.

3.6. Convergence Minkowski functionals

The HOS probes previously described can be considered extensions of the standard second-order probes since they all aim at probing the convergence field at higher orders to sort out the information contained in its non-Gaussianity. A different way to access this additional constraining power is represented by topological indicators.

Topology is the branch of mathematics that addresses the shapes, boundaries and connectivity of structural features in a field. More precisely, it is concerned with the properties of a geometric object that are preserved under continuous deformations, such as stretching, twisting, crumpling, and bending. Basic examples of topological properties are: the dimension distinguishing between a line and a surface; the compactness differentiating between a straight line and a circle; and connectedness, which separates a circle from two nonintersecting circles.

The first topological statistics we use as a cosmological tool are the Minkowski functionals (MFs). According to Hadwiger’s theorem (Hadwiger 1957), under a few simple requirements, any morphological descriptor of a d-dimensional scalar field is a linear combination of d + 1 MFs V_n with n ranging from 0 to d. The geometrical interpretation of each functional depends on the considered dimension d. Therefore, the morphology of the 2-dimensional convergence scalar field κ(ϑ) = κ(ϑ₁, ϑ₂) with variance σ² will be described by the first three MFs⁶ as

$$\begin{array}{cc}& V_0(\nu )=\frac{1}{A}{\displaystyle {\int }_{Q_{\nu }}\mathrm{d}a,}\hfill \\ \hfill & V_1(\nu )=\frac{1}{4A}{\displaystyle {\int }_{\partial Q_{\nu }}\mathrm{d}l,}\hfill \\ \hfill & V_2(\nu )=\frac{1}{2\pi A}{\displaystyle {\int }_{\partial Q_{\nu }}\mathcal{K}\mathrm{d}l,}\hfill \end{array}$$(26)

where A represents the map area, ∂Q_ν is the boundary of the excursion set Q_ν = {(ϑ₁, ϑ₂)∣κ(ϑ₁, ϑ₂)/σ ≥ ν} at a given threshold ν, da and dl are respectively the surface and the line element along ∂Q_ν, and 𝒦 is the local geodesic curvature of ∂Q_ν. Qualitatively, one can say that V₀(ν) and V₁(ν) quantify, respectively, the area and the perimeter of the excursion set Q_ν, while V₂(ν) gives the Euler characteristic, which is the topological quantity measuring the connectivity of the field. For an orientable surface, the latter is defined as χ = 2(𝒜 − ℬ), that is the difference between the number 𝒜 of disconnected regions above the threshold ν and the number ℬ of those below ν, that is the number of holes. We note that in Morse theory the Euler characteristics are also given by the alternating sum of critical point counts and hence related to the peak counts described in the previous section.

Equation (26), applied to the convergence field, can be rewritten in an operative form as

$$\begin{array}{cc}& V_0(\nu )=\frac{1}{A}{\displaystyle {\int }_A{\mathrm{d}}^2\mathit{\vartheta }\mathcal{H}(\kappa -\nu \sigma ),}\hfill \\ \hfill & V_1(\nu )=\frac{1}{4A}{\displaystyle {\int }_A{\mathrm{d}}^2\mathit{\vartheta }{\delta }_{\mathrm{D}}(\kappa -\nu \sigma )g(\kappa ),}\hfill \\ \hfill & V_2(\nu )=\frac{1}{2\pi A}{\displaystyle {\int }_A{\mathrm{d}}^2\mathit{\vartheta }{\delta }_{\mathrm{D}}(\kappa -\nu \sigma )h(\kappa ),}\hfill \end{array}$$(27)

where

$$\begin{array}{cc}& g(\kappa )=\sqrt{{\kappa }_{,1}^2+{\kappa }_{,2}^2}\hfill \\ \hfill & h(\kappa )=\left(\frac{2{\kappa }_{,1}{\kappa }_{,2}{\kappa }_{,12}-{\kappa }_{,1}^2{\kappa }_{,22}-{\kappa }_{,2}^2{\kappa }_{,11}}{{\kappa }_{,1}^2+{\kappa }_{,2}^2}\right),\hfill \end{array}$$(28)

with κ_,i and κ_,ij the first- and second-order partial derivatives of the convergence field with respect to ϑ_i and ϑ_iϑ_j and for i, j = 1, 2.

Equation (27) suggests that the zeroth MF V₀ can be evaluated by the integration of the Heaviside step function over the whole excursion set Q_ν. The other MFs V₁, V₂ have been transformed from line integrals into surface integrals by inserting a delta function δ_D and the appropriate Jacobian. The integrands of the first and second MFs, including the geodesic curvature 𝒦, are then given by second-order invariants that depend solely on the threshold ν, the field value κ and functions of its first- and second-order derivatives κ_,i, κ_,ij defined in Eq. (28), (see Sect. 2.3 of Schmalzing & Górski 1998, for a detailed calculation)⁷.

Equation (27) is the basis for developing an algorithm able to numerically measure MFs for a given threshold on a given convergence map. For our analysis, we make use of the code proposed and implemented in Vicinanza et al. (2019), Parroni et al. (2020). We apply it to the convergence map after performing a Gaussian smoothing with four different apertures, namely (1$\stackrel{\prime }{.}$17, 2$\stackrel{\prime }{.}$34, 4$\stackrel{\prime }{.}$69, 9$\stackrel{\prime }{.}$37). Figure 11 shows the corresponding DV for the 2$\stackrel{\prime }{.}$34 smoothing. The code has been validated by measuring MFs on simulated Gaussian maps since, in this case, it is possible to analytically predict the expected values of V_n(ν). This is no longer the case for non-Gaussian fields. A perturbative approach can be used when deviations from Gaussianity are very small (see e.g., Pogosyan et al. 2016), but it gives approximate results that require the evaluation of high-order polyspectra to be more accurate. Moreover, this formulation does not account for the presence of noise and systematics from the convergence reconstruction from biased shear data. We refer the reader to Parroni et al. (2020) and references therein for the derivation and testing of an approximate formulation dealing with these issues. However, in that same paper, it has been shown that a large number of nuisance parameters must be introduced, so it is desirable to rely on direct measurements from simulated convergence maps as we do here, which can be used to develop future emulators.

Fig. 11. Fiducial DV (gray) along with its derivatives with respect to Ωm (orange), σ8 (green), and w (blue) for MFs V0 (top), V1 (center), and V2 (bottom) computed on κ maps smoothed with a Gaussian filter with a scale of 2$\stackrel{\prime }{.}$34.

3.7. Convergence Betti numbers

In order to individually study the different topological features in the maps (holes, islands), one can rely on the Betti Numbers (BNs) β_i(ν), (whose alternating sum gives back the Euler characteristics). BNs consist of topological invariants that formalize the pure topological information in a given field through the relationship between singularity structure and connectivity. They represent homology measures since they describe the spatial distribution of the critical points of the field and their connectivity. The homology of a 3-dimensional manifold M of a 2-dimensional scalar field such as the convergence κ(ϑ) = κ(ϑ₁, ϑ₂) is characterized by the first three BNs⁸, β_i:ℝ↦ℤ ∣ i = 0, 1, 2. The manifold M of κ is a subset of ℝ³ defined by the relation

$$\begin{array}{c}\hfill M=\{(\mathit{\vartheta },\kappa (\mathit{\vartheta }))\mid \mathit{\vartheta }\in {\mathbb{R}}^2\}\subset {\mathbb{R}}^3.\end{array}$$(29)

In particular, similarly to the excursion set Q_ν introduced in Sect. 3.6, it is useful to construct a superlevel set filtration at a threshold ν formally defined as

$$\begin{array}{c}\hfill M(\nu )=\{(\mathit{\vartheta },\kappa (\mathit{\vartheta }))\in M\mid \kappa (\mathit{\vartheta })/\sigma \ge \nu \},\end{array}$$(30)

with σ² the variance of the convergence field. For all the thresholds ν₁ ≥ ν₂ made inline M(ν₁)⊆M(ν₂) Indeed, lowering the threshold ν from ∞ to −∞, new points are included into the manifold M(ν). It is trivial to show that M(∞) = ∅ and M(−∞) = M.

Following the Morse–Smale theoretical approach introduced in Feldbrugge et al. (2019), it is possible to infer integral relations for BNs valid for an arbitrary 2-dimensional random field (not necessarily Gaussian):

$$\begin{array}{cc}\hfill {\beta }_0(\nu )& ={\displaystyle {\int }_{\nu }^{\infty }\{{\mathcal{N}}_{\max }(\nu )-[1-g(\nu )]{\mathcal{N}}_{\mathrm{saddle}}(\nu )\}\mathrm{d}\nu ,}\hfill \\ \hfill {\beta }_1(\nu )& ={\displaystyle {\int }_{\nu }^{\infty }[g(\nu ){\mathcal{N}}_{\mathrm{saddle}}(\nu )-{\mathcal{N}}_{\min }(\nu )]\mathrm{d}\nu ,}\hfill \\ \hfill {\beta }_2(\nu )& =0,\hfill \end{array}$$(31)

with 𝒩_min,𝒩_saddle,𝒩_max:ℝ↦ℝ respectively the minima, saddle points and maxima density at threshold ν. The function g(ν):ℝ↦[0,1] and the opposite 1 − g(ν) respectively represent the probability of two different topological transitions caused by the introduction of a saddle point of the field at threshold ν. These latter affect the value of β₁ or β₀ as described in the incremental algorithm introduced in Feldbrugge et al. (2019). As is evident in Eq. (31), BNs depend solely on the diagnostic parameter ν, similarly to MFs. Furthermore, in the 2 - dimensional case, only the first two BNs β₀(ν),β₁(ν) are relevant since β₂(ν) will increase only when the lowest minimum exceeds the threshold ν. Unfortunately, there is no analytical expression for g(ν) in either the Gaussian or non-Gaussian case. In fact, only an approximated form has been found by Feldbrugge et al. (2019) for Gaussian 2 - dimensional fields. On the other hand, by relying on the incremental algorithm presented in Feldbrugge et al. (2019), codes can be developed in order to numerically count BNs. For our analysis, we measure BNs through a code developed in Parroni et al. (2021) and tested against Gaussian maps. As done for the MFs, the convergence field is smoothed with the same Gaussian filter and with identical apertures {1$\stackrel{\prime }{.}$17, 2$\stackrel{\prime }{.}$34, 4$\stackrel{\prime }{.}$69, 9$\stackrel{\prime }{.}$37}. The 2$\stackrel{\prime }{.}$34 scale DV used for the Fisher forecasts is displayed in Fig. 12.

Fig. 12. Fiducial DV (gray) along with its derivatives with respect to Ωm (orange), σ8 (green), and w (blue) for BNs β0 (top) and β1 (bottom) computed on κ maps smoothed with a Gaussian filter with a scale of 2$\stackrel{\prime }{.}$34.

3.8. Persistent homology of aperture mass

Persistent homology is a tool from topological data analysis to quantify the topological structure of data in the presence of noise. We briefly introduce the method here. Further details can be found in Heydenreich et al. (2021, 2022).

As in Sect. 3.7, we construct excursion sets on a manifold M; here, M is represented by the S/N map of aperture masses M_ap/σ(M_ap) on a square patch on the sky, seen as a subset of ℝ². The filter function of the aperture mass map was optimized to detect dark matter halos and reads (Schirmer et al. 2007)

$$\begin{array}{cc}\hfill Q_{\theta }(\vartheta )=& {[1+exp(6-150\frac{\vartheta }{\theta })+exp(-47+50\frac{\vartheta }{\theta })]}^{-1}\hfill \\ \hfill & \times {\left(\frac{\vartheta }{x_{\mathrm{c}}\theta }\right)}^{-1}tanh\left(\frac{\vartheta }{x_{\mathrm{c}}\theta }\right).\hfill \end{array}$$(32)

Here we use the filter scale of θ = 2$\stackrel{\prime }{.}$34 and a concentration index of x_c = 0.15 (Hetterscheidt et al. 2005). By varying the threshold, topological features such as connected components and holes emerge. The number of these features at varying threshold levels depends on the topological structure of the aperture mass map. The topological features also have a physical interpretation: Connected components correspond to local minima of the aperture mass map, and thus to underdensities in the integrated matter distribution. On the other hand, holes correspond to local maxima of the aperture mass map, or overdensities in the integrated matter distribution.

In contrast to BNs, which just count the number of features at each filtration level, persistent homology tracks the birth b and death d of each topological feature. The “persistence” of a feature, defined as its “lifetime” d − b, is a measure of how far a topological feature protrudes from its surroundings. In other words, features with low persistence are more likely to be attributed to noise (Sousbie 2011).

To visualize a persistent homology statistic, the features (b, d) are shown in a scatter plot, which is called a “persistence diagram”. Here, we distinguish between two summary statistics that we generate from the persistence diagram. The first are the “persistent Betti numbers” (pers. BNs) β_i(ν_b, ν_d). In contrast to regular Betti numbers β_i(ν), which just count the number of features alive at threshold ν, persistent Betti numbers β_i(ν_b, ν_d) count all features that were born before ν_b and die after ν_d. This allows us to take into account the persistence of topological features and disregard the features that are only short-lived (for more detail, see Heydenreich et al. 2021). As a more robust summary statistic, we generate “persistent heatmaps” (pers. heat.) by replacing each point in the persistence diagram by a Gaussian of width σ = 0.2.⁹ For both summary statistics, we then construct a data vector by employing the χ²-maximiser method discussed in Heydenreich et al. (2022). This method evaluates how much a pers. BNs or pers. heat. varies between different cosmologies with respect to the expected standard deviation across realizations, and then picks a sample of evaluation points that maximize the cosmological information content. For both pers. BNs and pers. heat., we generate a DV containing 30 entries each and show them in Fig. 13. The jagged features in these DVs are due to the ordering of the x-axis which is dictated by the cosmological information content of the respective point in the heatmap, and not by its value.

Fig. 13. Fiducial DV (gray) along with its derivatives with respect to Ωm (orange), σ8 (green), and w (blue) for pers. BNs (top) and pers. heat. (bottom) computed on aperture mass maps from γ, smoothed with a 2$\stackrel{\prime }{.}$34 Schirmer et al. (2007) filter. The x-axis corresponds to compression coefficients with no physical meaning.

3.9. Scattering transform coefficients

The scattering transform (Mallat 2012) was invented as a novel summary statistic which borrows ideas from convolutional neural networks but does not need training. It was first applied in computer vision, and shortly after introduced to cosmology in Cheng et al. (2020), it has gained an increasing attention for various applications (e.g., Cheng & Ménard 2021b; Valogiannis & Dvorkin 2022a,b; Chung 2022; Delouis et al. 2022; Greig et al. 2023). Below, we present a brief description of the ST. A nontechnical guide to this new statistic and more intuitive understandings and details can be found in Cheng & Ménard (2021a) and references therein.

The ST statistics is similar to a lower order and compressed version of the N-point polyspectra. It keeps the ability to describe different morphological configurations without raising the power of variables by multiplication, thus numerically much more stable. In addition, it also uses logarithmic binning of frequencies, which significantly reduces the size without losing much information for fields with multiscale structures in comparison to polyspectra.

In this study, we use up to second order ST, which resembles up to some binned tri-spectrum and is defined in the following way:

$$\begin{array}{c}\hfill \{\begin{array}{c}{\displaystyle S_0=⟨\kappa ⟩=0}\hfill \\ \hfill \\ {\displaystyle S_1^{j_1,l_1}=⟨|\kappa \star {\psi }^{j_1,l_1}|⟩}\hfill \\ \hfill \\ {\displaystyle S_2^{j_1,l_1,j_2,l_2}=⟨||\kappa \star {\psi }^{j_1,l_1}|\star {\psi }^{j_2,l_2}|⟩}\hfill \end{array},\end{array}$$(33)

where we apply the ST to the lensing convergence field κ. In the above equations, κ ⋆ ψ is a convolution with wavelets ψ labeled by a size index j and orientation index l. Wavelets are localized band-pass filters, so convolution with a wavelet selects Fourier modes in a wide frequency bin. The wavelets chosen here are qualitatively similar to the aperture mass filters, but they are oriented instead of isotropic and thus also complex-valued, and we use wavelets with a wide range of sizes. | ⋅ | denotes the modulus of the complex-valued fields, which is the key operation in the design of the ST. Finally, S_0, 1, 2 are the zeroth, first, and second order scattering coefficients, which are translation invariant.

Following Cheng et al. (2020), we use Morlet wavelets and sample L = 4 different orientations. Morlet wavelets is a sinusoidal wave multiplied by a Gaussian profile and modified to be strictly band-pass:

$$\begin{array}{c}\hfill {\psi }_0(x,y)=\frac{1}{\sigma \sqrt{2}}{\mathrm{e}}^{-(x^2+y^2)/2{\sigma }^2}({\mathrm{e}}^{i{k}_{0}x}-\beta ),\end{array}$$(34)

where $\beta ={\text{e}}^{-k_0^2\sigma /2}$, σ = 0.8 pixel, $k_0=\frac{3\pi }{4}$ pixel⁻¹. The whole library of wavelet ψ is obtained by dilating this prototype wavelet by factors of 2 for j times and rotating it by l × 45 degrees.

Now we explain the meaning of the ST in some detail. Those coefficients are obtained by first nonlinearly transforming the input field in a nested way and then taking the spatial average. In Fourier space, wavelets are wide window functions with different sizes (labeled by j) and orientations (labeled by l), so wavelet convolution will select Fourier modes around a frequency. The nonlinear operation, modulus, is not common in traditional statistics but proves useful in computer vision. Note that the modulus is applied in real space to each pixel (not to Fourier coefficients), similar to the pointwise nonlinearity used in convolutional neural networks. One intuitive way to understand the role of modulus is by noticing that for any field $|f(x,\mathit{y})|=\sqrt{f(x,\mathit{y})f^{\ast }(x,\mathit{y})}$ is the square root of a quadratic form of f. Cheng & Ménard (2021a) explained that the field |κ ⋆ ψ| is similar to a map of locally measured square root of the power spectrum averaged over the pass-band of filter ψ. So S₁ resembles a binned power spectrum, and S₂ resembles the spatial variance of power spectrum (4-point information), but all in a low-order form.

The total number of scattering coefficients depends on the number of distinct wavelets used. We find that the orientation dependence and the coefficients with j₁ ≤ j₂ contain little cosmological information for weak lensing, so the scattering coefficients can be reduced to a smaller set by averaging over all orientations:

$$\begin{array}{c}\hfill s_1^{j_1}={⟨S_1^{j_1,l_1}⟩}_{l_1}\\ \hfill s_2^{j_2}={⟨S_2^{j_1,l_1,j_2,l_2}⟩}_{l_1,l_2},\end{array}$$(35)

with j₁ ≤ j₂. If in total J dyadic scales (spaced by power of 2 in size) are used, the first and second-order scattering coefficients have sizes of J and J(J + 1)/2, respectively. In our study, we explore J = 8 dyadic scales, resulting in a compact set of 44 scattering coefficients used as the summary statistic. In Fig. 14, we show the mean and dispersion of the scattering transform coefficients in fiducial cosmology and its cosmological sensitivity. The jagged features are again due to the ordering of the x-axis. ST is a function of two scales, and would be a smooth function on a two-dimension space; however to display it in one dimension, we need to slice it along one direction, creating these wiggles. Besides compactness, the ST statistics also have the desirable properties of being informative and numerically stable, which is confirmed in the cosmological forecast analysis and reliability tests of the following sections.

Fig. 14. Fiducial data vector (gray) along with its derivatives with respect to Ωm (red), σ8 (green), and w (blue) for ST computed on κ maps smoothed with a Gaussian filter with a scale of 2$\stackrel{\prime }{.}$34. The first 8 points correspond to the first-order scattering coefficients, followed by 8 series of second-order scattering coefficients, one per scale number j1 and spanning all j2 values for each.

4. Fisher analysis

The aim of this paper is to estimate the constraints on (a subset of) CPs for each of the HOS probes described in Sect. 3. To this end, one could perform a Monte Carlo Markov chain (MCMC) analysis fitting the data as measured on, for example, the SLICS simulated maps. For this approach to be feasible, a mandatory ingredient would be a way to estimate the expected HOS value for any given set of CPs – in other words, a theoretical model and/or an emulator-like method. For most of the probes of interest here, these are not available, and so we need to look for an alternative approach that can be implemented for all the statistics. This is the Fisher matrix method, which is routinely used to forecast the accuracy that a probe can achieve on constraining CPs given the specifics of the survey at hand. In the following, we first give a brief introduction to the method, highlighting the key quantities, and then discuss in some more detail how we estimate them and the necessary caveats.

4.1. Fisher formalism

According to Bayes’ theorem, the posterior distribution P(p | D_obs) of the model parameters p given the observed data D_obs is given by

$$\begin{array}{c}\hfill P(\mathbf{p}|{\mathbf{D}}_{\mathrm{obs}})=\frac{\mathcal{L}({\mathbf{D}}_{\mathrm{obs}}|\mathbf{p})\mathrm{Pr}(\mathbf{p})}{\mathcal{E}({\mathbf{D}}_{\mathrm{obs}})},\end{array}$$(36)

where ℒ(D_obs | p) is the likelihood of the data D_obs for the parameters p, Pr(p) is the prior on the model parameters, and ℰ(D_obs) is the evidence. The Fisher matrix F is defined as the expectation value of the second derivative of the likelihood function, that is its elements are (Bunn 1995; Vogeley & Szalay 1996; Tegmark et al. 1997)

$$\begin{array}{c}\hfill F_{\alpha \beta }=-\langle \frac{{\partial }^{2}ln\mathcal{L}({\mathbf{D}}_{\mathrm{obs}}|\mathbf{p})}{\partial p_{\alpha }\partial p_{\beta }}\rangle ,\end{array}$$(37)

where (α, β) label the parameters of interest, and the derivatives are evaluated at a fiducial point p_fid assumed to be the maximum of the likelihood function. The Fisher matrix is then the Hessian of the logarithmic likelihood, thus estimating how fast it decreases from the maximum.

Under the usual assumption of a Gaussian likelihood function, the Fisher matrix elements reduce to

$$\begin{array}{cc}\hfill F_{\alpha \beta }& =\frac{1}{2}\mathrm{tr}\left[\frac{\partial {\mathsf{C}}_{\mathrm{d}}}{\partial p_{\alpha }}{\mathsf{C}}_{\mathrm{d}}^{-1}\frac{\partial {\mathsf{C}}_{\mathrm{d}}}{\partial p_{\beta }}{\mathsf{C}}_{\mathrm{d}}^{-1}\right]\hfill \\ \hfill & +\frac{\partial {\mathbf{D}}_{\mathrm{fid}}(\mathbf{p})}{\partial p_{\alpha }}{\mathsf{C}}_{\mathrm{d}}^{-1}\frac{\partial {\mathbf{D}}_{\mathrm{fid}}(\mathbf{p})}{\partial p_{\beta }},\hfill \end{array}$$(38)

where

$$\begin{array}{c}\hfill {\mathsf{C}}_{\mathrm{d}}=\langle ({\mathbf{D}}_{\mathrm{obs}}-⟨\mathbf{D}⟩){({\mathbf{D}}_{\mathrm{obs}}-⟨\mathbf{D}⟩)}^T\rangle \end{array}$$(39)

is the data covariance matrix, and we have assumed that the mean of the data ⟨D⟩ coincides with the theoretical DV D_fid computed for the fiducial parameters. As in most lensing studies, we assume that the covariance matrix is model independent so that the Fisher matrix reduces to the second term only in Eq. (38). Strictly speaking, the Fisher matrix is the expectation of the likelihood Hessian, while in Eq. (38), we have implicitly assumed that this is the same as computing the derivatives at the fiducial point. Because of the Cramer–Rao inequality, this means that the inverse of the Fisher matrix gives a lower limit to the full expected parameter covariance matrix of the CPs. That is to say, the marginalized errors (i.e., having included all of the degeneracies with respect to other parameters) estimated from the diagonal elements C_αα of C as

$$\begin{array}{c}\hfill \sigma (p_{\alpha })=C_{\alpha \alpha }^{1/2}={\left[{({\mathsf{F}}^{-1})}_{\alpha \alpha }\right]}^{1/2},\end{array}$$(40)

are a lower limit to the actual errors on the CPs.

The Fisher matrix is also useful for estimating the degeneracy directions in parameter space. The correlation coefficient among the parameters (p_α, p_β) is given by

$$\begin{array}{c}\hfill {\rho }_{\alpha \beta }=C_{\alpha \beta }/\sqrt{C_{\alpha \alpha }{C}_{\beta \beta }},\end{array}$$(41)

while the angle ϕ_αβ defining the degeneracy direction in the 2D plane (p_α, p_β) may be estimated as

$$\begin{array}{c}\hfill {\varphi }_{\alpha \beta }=\frac{1}{2}arctan\left(\frac{C_{\alpha \beta }}{C_{\alpha \alpha }-C_{\beta \beta }}\right).\end{array}$$(42)

Both Eqs. (41) and (42) show that two parameters are completely independent if C_αβ = 0, that is the corresponding off-diagonal term of the Fisher matrix vanishes.

The above description of the Fisher matrix formalism makes it clear that three assumptions are implicitly used when computing the Fisher matrix: i. the likelihood may be approximated as Gaussian; ii. a reliable estimate of the covariance matrix is available; iii. one knows how to compute the derivatives of the DV with respect to the CPs. We subsequently discuss in the next few paragraphs the caveats associated with each of these points and how we deal with them.

4.2. Gaussian likelihood and data selection

Equation (38) has been obtained under the assumption that the likelihood is Gaussian. However, such an assumption must be tested for each particular probe at hand since it is manifestly violated for some of them under certain conditions. As an example, one can think of the number of peaks at very large κ, which is expected to be zero in the absence of noise. Similarly, the value of the zeroth order MF V₀(κ) is identically 1 for very small κ. As a consequence, every measurement of these quantities will give the same value so that the distribution of repeated measurements would definitely be non-Gaussian. Should we include them in the DV, the Gaussian assumption of the likelihood would break down, introducing a systematic error that is hard to quantify. We therefore need to perform data selection to be confident that the likelihood can be well approximated by a Gaussian.

To this end, we note that if the likelihood is indeed Gaussian, the quantity

$$\begin{array}{c}\hfill y_i=({\mathbf{D}}_i-⟨\mathbf{D}⟩){\mathsf{C}}_{\mathrm{d}}^{-1}{({\mathbf{D}}_i-⟨\mathbf{D}⟩)}^{\mathrm{T}},\end{array}$$(43)

with D_i the DV measured on the i-th SLICS map, and ⟨D⟩ the mean over the 924 realizations, must be distributed as a χ² distribution with N_d degrees of freedom (N_d being the length of the DV). For each HOS probe, we therefore compare the y_i distribution to the χ² one, and quantify the agreement by computing the weighted average SMAPE (Symmetrized Mean Absolute Percentage Deviation) defined as (Rizzato & Sellentin 2023)

$$\begin{array}{c}\hfill \mathcal{S}=\frac{\sum {\omega }_i{\mathcal{S}}_i}{\sum {\omega }_i},\end{array}$$(44)

with

$$\begin{array}{c}\hfill {\mathcal{S}}_i=\frac{|P_{\mathrm{obs}}(y_i)-P_{{\chi }^2}(y_i)|}{\left|P_{\mathrm{obs}}(y_i)\right|+\left|P_{{\chi }^2}(y_i)\right|},\end{array}$$(45)

and

$$\begin{array}{c}\hfill {\omega }_i=P_{{\chi }^2}(y_i),\end{array}$$(46)

where P_obs(y_i) and P_χ²(y_i) are the measured and χ² distribution in y_i, and the sum is over the points used to sample the y distribution.¹⁰ Note that we have defined the weights ω_i in such a way that we ask for a better match in the central regions of the distributions, reducing the sensitivity to outliers. There is no general rule to define a limiting value for 𝒮, but it is clear that the closer it is to zero, the more the two distributions match. For each probe and configuration, we therefore implement a four-step procedure to decide whether the distribution is compatible with the one expected for a Gaussian probe according to the following scheme: i. We first remove all the elements in the DV that are consistent with zero within 2σ, with σ estimated from the corresponding diagonal element of the covariance matrix. Note that this step is equivalent to saying that we are removing points with S/N < 2, so that we are actually deleting those elements carrying no valuable information; ii. We recompute the weighted average SMAPE defined above using only the surviving elements and comparing them with a χ² distribution with N_obs degrees of freedom, and with N_obs ≤ N_d the number of surviving elements of the full DV; iii. We generate 924 mock realizations of the DV sampling from a Gaussian distribution with a mean and covariance matrix equal to the observed one, and compute the corresponding 𝒮_mock value; iv. We repeat the above step 500 times, and estimate the 95% confidence limit (CL) range of the distribution of 𝒮_mock, finally considering the observed DV Gaussian if its 𝒮_obs is within this range.¹¹ This was empirically found to roughly correspond to a limit of 𝒮_obs ≤ 0.15.

Applying this procedure to the full list of probes is, however, nontrivial given the different peculiarities. Roughly speaking, we can divide them into two classes. First, we have the κ-PDF as a function of κ, the number of peaks as a function of S/N ν, the BNs and the MFs versus the threshold ν. All are measured by first choosing a range (x_min, x_max), and splitting it into N_x bins (with x = κ, ν). Second, we have other probes that are measured for a given fixed configuration so that they are set once for all by the measurement pipeline. In the first case, we can vary (x_min, x_max, N_x) until the Gaussianity criterion has been successfully passed. We investigate, in particular, configurations with (κ_min, κ_max) in the range ( − 0.1, 0.1) and (ν_min, ν_max) in the range ( − 5, 10), varying the number of bins from 10 to 70 in steps of 5. We then retain the configurations passing the SMAPE test according to the above procedure. This is illustrated in Fig. 15, where we can appreciate the correspondence between the measured and theoretical χ² distribution for all probes but moments in configurations where 𝒮 ≤ 0.15.

Fig. 15. Verification of the Gaussian hypothesis for all statistics computed in this analysis. For a Gaussian-distributed DV, the histogram of χ2 values from the SLICS simulations (red) should match the theoretical prediction (black) modulo the sampling noise across the finite set of 924 realizations. The corresponding SMAPE value (Eq. (45)) is given in the top left quadrant of each panel. 𝒮 ≤ 0.15 indicates compatibility with the Gaussian hypothesis, which is verified for all statistics but the three bottom panels with gray shaded background. The y-axes correspond to the frequency of each χ2 bin value across the SLICS and are different for each panel.

For the other family of HOS probes, we could only compute the value of the average SMAPE, and checked whether it passed our Gaussianity criterion. We computed the Fisher matrix forecasts for all the probes, but left a warning if the SMAPE value is outside the expected 95% CL range, that is if the y_i distribution is too different from a χ² one. This is the case for all moments-based statistics: HOM, $⟨M_{\mathrm{ap}}^3⟩$, $⟨M_{\mathrm{ap}}^n⟩$ (gray shaded panels in Fig. 15), which have significant non-Gaussian distributions across the different realizations, rendering them unsuitable for a Fisher analysis. This issue, which seems partially related to the small patch sizes considered here, could be mitigated by considering "bulk" moments obtained from integrating the κ-PDF over a limited range using the tail cuts employed for the PDF. We nonetheless prefer to report the forecasts for these probes too since the fact that their likelihood is non Gaussian does not prevent one from using them in the future. Indeed, should the corresponding data be available, one can simply run an MCMC analysis based on a non-Gaussian likelihood. The classification of Gaussian and non-Gaussian DVs is given in Table 4, together with the selected range of scales, convergence, or S/N that allows each probe to pass the Gaussianity test.

As a final remark, we note that the implemented criterion looks at the overall properties of the DV. However, one could also apply a more stringent criterion looking at the distribution of each single element of the DV. A way to quantitatively implement this approach is based on the D’Agostino–Pearson test. It turns out that this approach is, however, overly restrictive, cutting a large part of the DV for most probes. This is because this method looks for a perfect match with a Gaussian likelihood even in the tail of the distributions, which are actually not of interest for the computation of the Fisher matrix. Indeed, the Fisher matrix requires that the likelihood can be approximated as Gaussian in the neighborhood of its peak, so one is not really interested in what happens in the tails. This is why we finally prefer to rely on the y_i distribution rather than the D’Agostino–Pearson test. This choice is also supported by the results of Lin et al. (2020), who found that the skewness and kurtosis of the γ-2PCF in LSST-like mocks have little impact on forecasts when computed with a Gaussian or an Edgeworth non-Gaussian likelihood.

4.3. Covariance matrices

The covariance matrix for each HOS probe can be estimated using the SLICS ellipticity mocks and convergence maps. To this end, we simply measure each probe on all realizations, take the mean as our fiducial DV, and use the definition (39) to compute the data covariance matrix C_d. We also use this simple method when investigating the joint use of second and HOS, or of more HOS probes at the same time. We then just concatenate the individual DVs into a single DV, which we then measure on the different mock realizations for the estimate of the covariance matrix. Unless otherwise stated, however, we focus first on single HOS probes and later combine them with the γ-2PCF. Moreover, to reduce the number of possible combinations, we choose a reference smoothing angle of 4 pixels for all probes, except for the convergence PDF and moments, which we evaluate from maps smoothed with an 8 pixel radius.¹² We have checked that, for each probe, there is indeed a strong correlation between values estimated from the same map with different smoothing radii, so that using one smoothing angle per statistics captures most of the cosmological information.

Using a finite number of realizations alters the accuracy of the inverse covariance matrix (e.g., Hartlap et al. 2007). This arises from the fact that the DV cannot be considered as Gaussian distributed (this is only the case for an infinite number of realizations according to the central limit theorem) but rather follows a student-t distribution. Sellentin & Heavens (2017) propagated the impact of using a student-t likelihood in the Fisher formalism instead of the Gaussian likelihood, deriving a simple correction factor to this effect to be applied to the inverse Fisher matrix, which depends on the number of realizations used in the computation of the covariance matrix N_f, the number of degrees of freedom, that is the DV size N_d, and the number of CPs to be determined N_p: c = (N_f − 1)/(N_f − N_d + N_p − 1), (see Eq. (28) of Sellentin & Heavens 2017).

Although we use the N_f = 924 SLICS maps to compute the covariance in our analysis, we also use the 256 DUSTGRAIN-pathfinder realizations as a cross-check. However, the fact that the DUSTGRAIN-pathfinder realizations are pseudo-independent (i.e., derived from a single simulation) renders the comparison with the truly independent SLICS difficult. This nonetheless allows us to quantify the impact of the choice of a simulation set on the covariance computation and to a lesser extent on the possible cosmology dependence of the covariance as the SLICS and DUSTGRAIN-pathfinder have slightly different CPs. This comparison is done at the level of the Fisher forecasts, keeping the derivatives fixed and using three different covariance matrices: computed from the 256 pseudo-independent DUSTGRAIN-pathfinder simulations, and 256 and 924 SLICS simulations. We find perfect agreement between the SLICS for the different number of realizations when factoring in the Sellentin & Heavens (2017) correction factor, highlighting the robustness of our covariance estimate in terms of the number of realizations used. However, we see some variations between SLICS and DUSTGRAIN-pathfinder. Overall, the order of magnitude of the forecasts remains the same and it is difficult to attribute the difference we measure to a particular effect. This could stem from the simulations themselves, in particular the differences in the computation of the lensing quantities discussed in Sect. 2, the differing fiducial cosmologies, the different area of the mocks (25 deg² for DUSTGRAIN-pathfinder and 100 deg² for SLICS), or to the level of independence of the realizations in the DUSTGRAIN-pathfinder simulations. To specifically assess whether this difference could be due to the different field of view of the SLICS and DUSTGRAIN-pathfinder mocks, we also tried the following experiment. We computed the covariance from the central 5 × 5 deg² region of each of the 924 SLICS mocks, scaled the covariance matrix by four, and compared it with the 10 × 10 deg² result. We found excellent agreement showing that the difference with the DUSTGRAIN-pathfinder 5 × 5 deg² covariance is not due to the change in area. The good agreement between the SLICS covariance for both areas also supports our decision to scale up the survey area to the nominal 15 000 deg² expected for Euclid. In the case of the κ-PDF, it has also been shown that the numerical DUSTGRAIN-pathfinder covariances agree with those obtained from shifted lognormal maps created with FLASK (Xavier et al. 2016) for a fixed source redshift, which can be created in large quantities and tuned to replicate the desired skewness thus exquisitely matching the more expensive simulated covariances (Boyle et al. 2021).

Ultimately we chose to use the SLICS covariance because of the larger field-of-view and number of realizations, and the independence of each realization of this simulation set. We also rescale the covariance to 15 000 deg². We note that our simulated covariances from small patches with vanishing mean density neglect the super-sample covariance effect. As its impact depends on the specific summary statistic, including it in the analysis is beyond the scope of this work. In the future, those effects can be estimated by treating the full data covariance as a sum of the simulated covariance and a super-sample term built from the response of the summary to a background density, as done for angular power spectra (Lacasa & Grain 2019) and the κ-PDF (Uhlemann et al. 2023). Alternatively, an analytic covariance model including super-sample effects can be derived from the estimator, as done for the κ-PDF (Uhlemann et al. 2023) and $⟨M_{\mathrm{ap}}^3⟩$ (Linke et al. 2023).

4.4. Derivatives and numerical noise

Computing the derivatives of the observable is a key step in the estimate of the Fisher matrix. This would be a trivial task should a theoretical procedure be available to compute the HOS corresponding to a given set of CPs. Unfortunately, this is not the case for most of the probes of interest here. For some of them, a theoretical formulation has been developed, but this typically does not account for noise in the data, and assumes that the κ map is perfectly reconstructed from the γ data.

Here, we must therefore compute derivatives numerically using the DUSTGRAIN-pathfinder simulations for shifted values of (Ω_m, σ₈, w). For each given probe, we first measure the DV as the mean over the N_c = 128 mocks available for each model. We then use four different methods to estimate numerical derivatives. First, we consider the simple 3 - point stencil derivatives defined as

$$\begin{array}{c}\hfill \frac{\partial \mathbf{D}}{\partial p_{\alpha }}=\frac{\mathbf{D}(p_{\alpha }+\epsilon p_{\alpha }^{\mathrm{fid}})-\mathbf{D}(p_{\alpha }-\epsilon p_{\alpha }^{\mathrm{fid}})}{2\epsilon |p_{\alpha }^{\mathrm{fid}}|},\end{array}$$(47)

with D(p_α) the DV measured on the maps generated from simulations with all CPs but p_α set to the fiducial values. Two choices are possible for the shift ε, the small and large increments, respectively corresponding to 4 and 16% except for the large increment of Ω_m for which ε takes the values +28%/−36%. Alternatively, we can use all four values as entries in a linear fit, taking the best fit slope as an estimate of the first derivative. We perform both an unweighted fit and a weighted fit with weights given by

$$\begin{array}{c}\hfill w_i^{\mathrm{fit}}={[1-{\left(\frac{|{\epsilon }_i-{\epsilon }_{\mathrm{ref}}|}{1.1{\epsilon }_{\mathrm{ref}}}\right)}^3]}^3,\end{array}$$(48)

with ε_ref = 0.12, and ε_i = ( − 0.16, −0.04, 0.04, 0.16), again with adapted values for the large increment of Ω_m. In practice, we find that the weighted and unweighted fits are both equal to the large increment derivatives so we only discuss the large and small increment cases in the following.

As an example, we show in Figs. 4–14 the numerical derivatives for the large CP variation of all probes with respect to Ω_m, σ₈, and w in red, green, and blue, respectively. Shaded areas correspond to the noise computed as the dispersion over the N_c = 128 DUSTGRAIN-pathfinder realizations at each cosmology. We also show in gray the mean DV over the N_f = 924 SLICS realizations with error bars corresponding to the diagonal elements of the SLICS covariance matrix to better visualize the trend of the derivatives on the DV. We note from these figures that the derivatives with respect to w are close to 0, suggesting that we actually measure very little information on this parameter. This could be alleviated in the future for instance by increasing the number of realizations in the simulation, by adding simulations with a larger step in w in order to strengthen its effect, or by including a tomographic decomposition. On the other hand, there is a strong signal for the other two parameters and we therefore focus the forecasts on Ω_m and σ₈ in the remainder of the paper.

As a further check on the reliability of the numerical derivatives, one can investigate whether they are stable against the number of DUSTGRAIN-pathfinder maps used to compute them. In Figs. 4–14, we therefore add the derivatives computed for N_c = 32 and 64 realizations as dotted and dashed lines, respectively. For the large variation derivatives, these lines are within the error bars computed from the 128 realizations, highlighting the robustness against numerical noise. This is, however, not always the case at the level of the Fisher forecasts. In Fig. 16, the covariance matrix is held fixed to the 924 SLICS realizations and the dotted, dashed, and solid ellipses correspond to the forecasts computed with Fisher derivatives of 32, 64, and 128 DUSTGRAIN-pathfinder realizations, respectively. We note that several probes suffer from a slow convergence, for example κ-PDF, MFs, BNs, and pers. heat.. The fact that the forecasts are not yet fully converged is due to the complex interplay between the errors on the derivatives and that of the covariance matrix, which can amplify some noise fluctuations in the derivatives. We demonstrate the trustworthiness of the not fully converged κ-PDF derivatives through a comparison with theoretical predictions that lead to marginally wider and slightly tilted Fisher contours as shown in Sect. 4.5.2. Although one could try to develop some numerical methods to correct for this bias, we did not do so because of the three following reasons. First, the convergence of the numerical derivatives is different for each probe with only half of them being affected with the current number of realizations of the DUSTGRAIN-pathfinder simulations, suggesting a nontrivial correction scheme that would introduce individual differences in the analysis of each probe. Second, the effect remains small as shown by the smaller difference between the 64 and 128 realizations than between the 32 and 64. Finally, the comparison with theoretical predictions in the case of the κ-PDF strongly suggests that convergence has been reached in spite of the small variation still being seen between the 64 and 128 realizations in Fig. 16. Although we are confident in the forecasts when using the large CP increment to compute the derivatives, when considering the small increment derivatives, these are not even converged at the level of the Fisher derivatives. This is due to the CP step being too small, so that the change in the observable due to the shift from the fiducial CPs gets partially lost into the numerical noise. Apart from a few additional tests when comparing numerical and theoretical derivatives in Sect. 4.5, we therefore only consider the ±16% derivatives in the rest of this paper.

Fig. 16. Fisher forecasts for all probes computed for a Euclid-like survey without tomography. The different ellipses refer to different numbers of DUSTGRAIN-pathfinder realizations Nc used to compute the DV derivatives: dotted, dashed, and solid lines, respectively for 32, 64, and 128 realizations. The covariance matrix is held fixed to SLICS with Nf = 924 realizations. For γ-2PCF and κ-PDF we separately assess the convergence of derivatives by a comparison with theoretical predictions in Figs. 17 and 18.

4.5. Theory versus numerical forecasts

For assessing the robustness of our method in creating reliable Fisher forecasts from the simulated maps, we compare them to the theoretical results for a set of statistics for which such a description is available, namely the γ-2PCF, κ-2PCF, κ-PDF, MFs, HOM, and $⟨M_{\mathrm{ap}}^3⟩$. We here only report on two of them: the γ-2PCF in Sect. 4.5.1 and κ-PDF in Sect. 4.5.2, but the general conclusions are the following. Numerical noise in the Fisher derivatives can generate some artificial degeneracy breaking in the forecasts from the simulation-based model compared to the theoretical ones. Except for moments where we could not reconcile the theory and simulations due to the small field-of-view of our mocks, we find a fair agreement for other probes when applying some restrictions to our analysis:

• we only consider Fisher forecasts for the two parameters Ω_m and σ₈ as the derivatives with respect to w are too noisy in our nontomographic setup;

• we use only the large variations of the CPs (±16%) as the small variation (±4%) is too small to capture the signal compared to the numerical noise;

• we further restrict the range of scales used in the 2PCFs because of the resolution and field-of-view of the DUSTGRAIN-pathfinder mocks.

4.5.1. Shear two-point correlation functions

As a first example, we choose the γ-2PCF introduced in Sect. 3.1. For computing the theoretical prediction, we use a modified version of the public software NICAEA (Kilbinger et al. 2009) and choose the transfer function from Eisenstein & Hu (1998) and the revised Halofit prescription of Takahashi et al. (2012) for the modeling of the nonlinear power spectrum. We also account for finite field effects in the ray-tracing of the DUSTGRAIN-pathfinder mocks by restricting the support of the convergence power spectrum in Eq. (7) to the accessible field-of-view; in particular we reweight P_κ by the fraction of inverse ℓ-modes that can be confined to a 5 × 5 deg² region. This modification introduces a strong ringing effect in the γ-2PCF which becomes significant at around 5′/150′ for ξ₋/ξ₊. In the upper panels of Fig. 17, we compare the theoretical and numerical derivatives and we see that there is a nonnegligible discrepancy on small scales for both ξ₊ and ξ₋ and an additional discrepancy on large scales for ξ₊.

Fig. 17. Upper row: Comparison of the derivatives of the γ-2PCF as predicted from theory (solid lines) and as measured from the DUSTGRAIN-pathfinder simulations (dots). The errorbars are scaled to a single line-of-sight in the DUSTGRAIN-pathfinder ensemble. The shaded regions display the scales that were discarded from the Fisher analysis presented in this work. In the lower panel we plot the relative deviation between the measured quantities and the theoretical predictions when accounting for (solid lines) or neglecting (dotted lines) finite field effects. The gray dashed lines display the 10 percent errorband. Lower row: Comparison of the Fisher forecast using the theoretical model (orange dashed) or the simulation measurements (blue dashed). The thin solid lines correspond to simulated analyses that are used to test the stability of the ellipses given the numerical noise in the derivatives. For the panel on the left we use all available scales on the γ-2PCF while for the figure on the right only the scales with realistic numerical derivatives (i.e., removing the gray shaded area of the upper panel) are included in the forecast. In the figure on the right we show for reference again the constraints of the theoretical model with all scales included.

In the lower-left panel of Fig. 17, we show how these biases propagate in the forecasts with fixed w. In particular, when including a broad range of scales, we find that the simulation-based constraints seem much stronger as the different level of bias in the derivatives artificially breaks parameter degeneracies. When using only the scales for which the numerical derivatives are in reasonable agreement with their theoretical counterparts, this effect disappears at the cost of overall reduced constraining power. From the theoretical point of view this reduction remains small as shown by comparing the orange and black ellipses in Fig. 17. Besides the bias, noise in the simulation-based derivatives can also lead to artificial degeneracy breaking. To assess the magnitude of the latter effect, we first estimate the covariance matrix of the simulation-based derivatives and, from there, generate 100 noise realizations. When adding these realizations to the theoretical derivatives, we can assess the stability of our forecast in the ideal case, that is when no bias is present in the simulations. On the other hand, when adding the noise on top of the numerical derivatives, we get an upper bound on the instability of the actual forecast, even without a theoretical description of the observables underlying the analysis. For the γ-2PCF, we find that the ellipses are stable in all cases. For our subsequent analysis we limit ourselves to those scales in which the mean bias on the derivatives in (Ω_m, σ₈) is less than 10%. This range corresponds to the scales displayed in Table 4. We note that it is more restrictive than the range used in Martinet et al. (2021b,a) based on the same SLICS mocks for the covariance because of the better resolution and larger field-of-view of the cosmo-SLICS used for the model compared to the DUSTGRAIN-pathfinder simulations here.

4.5.2. Convergence PDF

We can use the theoretical model for the κ-PDF discussed in Sect. 3.2 to obtain the PDF at different cosmologies and cross-validate it with the PDFs measured from the simulations, which are also used to obtain the PDF covariance matrix entering the Fisher forecast.

We compare the mean of the measured κ-PDF for the fiducial cosmology with the theoretical predictions (including shape noise) in the top part of Fig. 18 (gray points and solid line) for a smoothing scale of θ = 4$\stackrel{\prime }{.}$69. We see that for the chosen cuts in the tails of the PDF and the measured error bars (standard deviation across all realizations) the theory does perform pretty well, which is reinforced by the fact that the derivatives with respect to the CPs are also well captured. Measured derivatives for both smaller increment sizes (shown as stars) and larger increment sizes (shown as points) in the CP variation are compared with our theoretical model for the κ-PDF, where the nonlinear variances along the line of sight are set externally. For σ₈ (green), we see very good agreement between the theoretical and simulated derivatives. We also show the results for Ω_m (red). In this case, a more substantial degree of discrepancy is seen between the theoretical and measured derivatives. The underprediction of the amplitude of the derivative stems from an underprediction of the difference in variances between the two cosmologies as predicted by Halofit. This is consistent with the simulations predicting a stronger response to Ω_m than the theory in Fig. 17. We found a similar discrepancy in the Ω_m derivatives when cross-checking with a larger smoothing scale of 9$\stackrel{\prime }{.}$37. We observe a similar underprediction for the w derivative (blue), which is also consistent with results for the γ-2PCF in Fig. 17. We also checked that replacing the input Halofit nonlinear variance with that from the Euclid emulator in each case resulted in minor changes of the predictions.

Fig. 18. Top: Overall shape of the lensing κ-PDF, smoothed with a top-hat filter of radius 4$\stackrel{\prime }{.}$69 and including shape noise. The fiducial cosmology is shown in gray as predicted from large deviation theory (line) and measured in the DUSTGRAIN-pathfinder simulations (data points with error bars indicating the standard deviation across the 256 realizations). Derivatives of the κ-PDF with respect to the CPs are shown in color, with theoretical predictions using nonlinear variances predicted by Halofit shown as solid lines. The markers represent the derivatives obtained from the simulations using finite differences based on larger and smaller increments (points and stars, respectively). Bottom: Fisher forecast constraints on Ωm and σ8 for a Euclid-like survey from the κ-PDF shown above. Numerically measured (orange) and theoretically predicted (pink) derivatives find good agreement, with differences in the final 1σ parameter constraints of approximately 20%.

Having discussed the theoretical modeling and measurement of the WL κ-PDF, we perform a Fisher forecast for the two parameters Ω_m and σ₈ using the central region of the WL κ-PDF at a single scale assuming a Euclid-like survey area, shape and source redshift distribution as discussed before. The result shown in the bottom part of Fig. 18 shows reasonably good convergence between the theoretical and numerical results, which reinforces the validity of both sets of forecasts. When adding w to the forecast, the simulated DVs would lead to unrealistically tight constraints compared to the predicted DVs. To avoid an artificial degeneracy breaking due to numerical noise, we limit the following analysis to two CPs.

5. Results

We now discuss the results of the Fisher analysis. We start by considering each probe independently (Sect. 5.1), then quantify the gain on CPs from our statistics when combined with the γ-2PCF alone (Sect. 5.2), and finally investigate the correlations between probes and their combinations (Sect. 5.3). For all results in this section of the paper, we consider a Euclid-like nontomographic analysis of 15 000 deg² and vary only Ω_m and σ₈, as the Fisher derivatives are too noisy for constraining all three parameters (i.e., including w) in the nontomographic setup. We also present forecasts for the growth of structure parameter Σ₈ = σ₈(Ω_m/0.3)^α (equivalent to S₈ when α = 0.5) in Appendix A to ease comparison with other analyses; however these results lead to similar conclusions as for σ₈. The mass maps are smoothed with a 2$\stackrel{\prime }{.}$34 Gaussian, wavelet or M_ap filter except for moments and the κ-PDF where a 4$\stackrel{\prime }{.}$69 top-hat filter is used, thus probing larger scales. The derivatives are computed with the large variation of CPs in the DUSTGRAIN-pathfinder simulations and the covariance from the 924 SLICS simulations. We also note that these results were cross-checked by running two independent Fisher analysis pipelines that obtained the same output. Finally, we make a distinction between results for Gaussian-distributed DVs, for which the Fisher forecasts can be robustly computed, and results for non-Gaussian DVs, where the Fisher formalism cannot be trusted. The latter are included for consistency of the HOWLS project and because they can be used in future analyses by using a non-Gaussian likelihood.

5.1. Individual forecasts

Individual Fisher forecasts are displayed in columns 2 and 3 of Table 5 and correspond to the 1σ marginalized parameter errors given as percentages of the CP fiducial values: δσ₈/σ₈ and δΩ_m/Ω_m. These forecasts can also be visually appreciated as the orange Fisher ellipses in the Ω_m − σ₈ planes of Fig. 19.

Fig. 19. Individual Fisher forecasts in the σ8 − Ωm plane for a nontomographic Euclid-like survey for the 11 statistics (orange) and γ-2PCF (blue), as well as their combination (green). The corresponding marginalized precision on CPs can be found in Table 5. The black dashed line in the first quadrant indicates constant $S_8={\sigma }_8\sqrt{{\mathrm{\Omega }}_{\mathrm{m}}/0.3}$. The bottom row with the gray shaded backgrounds displays probes that are not Gaussian distributed and cannot be robustly interpreted with Fisher forecasts. The abbreviated name of each summary statistic is displayed in the top-right part of each panel: γ-2PCF and κ-2PCF for the shear and convergence two-point correlation functions, κ-PDF for the convergence one-point probability distribution, peaks for aperture mass peak counts, MFs for convergence Minkowski functionals, BNs for convergence Betti numbers, pers. BNs and pers. heat. for aperture mass persistent homology Betti numbers and heatmap, ST for convergence scattering transform coefficients, HOM for higher-order convergence moments, and $⟨M_{\mathrm{ap}}^3⟩$ and $⟨M_{\mathrm{ap}}^n⟩$ for third and n-th order aperture mass moments.

If we first focus on second-order statistics, we see that the forecasts are on the order of one percent on both parameters. This is the same order of magnitude as the forecasts from Euclid Collaboration (2020). A quantitative comparison is, however, not possible due to large differences in both analyses: we vary only 2 CPs in a nontomographic setup while they use a minimum of 5 parameters and a tomographic setup. Additionally, their analysis is based on theoretical predictions, while we are affected by numerical noise because we need to rely on simulations for a fair comparison among HOS. If we compare to analyses more similar to ours, for instance, the cosmo-SLICS based likelihood inference of Martinet et al. (2021b) for 4 CPs in their nontomographic setup, we also find a good agreement: rescaling the forecasts of their Table 1 to 15 000 deg² gives a 1.9% precision on Ω_m to be compared with the 1.2% of the present analysis. Our Fisher forecasts are tighter because we neglect the degeneracies between our two parameters and the additional two (w and h) probed in the aforementioned analysis. Additionally, we know that degeneracies between pairs of CPs, which are treated as simple ellipses in the Fisher formalism, are usually more complex. Looking at the orientation of the blue Fisher ellipse in the top-left quadrant of Fig. 19, we also see the usual lensing degeneracy shown by the dashed black line of constant S₈. As expected, the κ-2PCF forecasts are of the same order as the γ-2PCF, but weaker due to the pixelization of the κ map, which destroys signal below the pixel size of 0$\stackrel{\prime }{.}$59.

The κ-PDF forecasts are shown as the orange ellipse in the top-right quadrant of Fig. 19. The forecasts in this case are on the order of half a percent and are therefore tighter than those of the γ-2PCF thanks to the non-Gaussian small-scale information that this statistic is able to probe. The forecasts in this case also show different degeneracy directions. This is also the case for other HOS, with some diversity between DVs ranging from forecasts on the order of a third of a percent to about one percent, and various degeneracy orientations as shown by the different Fisher ellipses of Fig. 19. In the left part of Fig. 20 we also make an attempt to rank the statistics in terms of constraining power. We see a plateau in this figure where most of the statistics reach similar forecast values, about twice as good as that from the γ-2PCF. It is therefore hard to identify a single superior statistic, with notably pers. BNs, MFs, ST, κ-PDF, peaks, and pers. heat. all performing very well.

Fig. 20. Fisher forecasts for σ8 and Ωm for a nontomographic Euclid-like survey for all statistics. The black crosses (left) show the individual forecasts and the gray circles (right) the forecasts combined with the γ-2PCF. Forecasts are ranked within 3 categories: two-point statistics, higher-order Gaussian statistics, and higher-order non-Gaussian statistics. The latter are displayed on a gray background as they should be confirmed with a non-Gaussian likelihood framework. The corresponding marginalized precision on CPs can be found in Table 5.

Finally, we also display the forecasts for the statistics which are not Gaussian distributed at the end of Table 5 and with a gray shaded background in Figs. 19 and 20: HOM, $⟨M_{\mathrm{ap}}^3⟩$, and $⟨M_{\mathrm{ap}}^n⟩$. These are only displayed for the consistency of the HOWLS challenge but must be interpreted with caution as the Gaussian likelihood assumption is required for the Fisher formalism. There is possible strong constraining power but a refined analysis that goes beyond the simple Fisher method is needed to assess the robustness of these results. We however note that in the case of $⟨M_{\mathrm{ap}}^3⟩$ our results are consistent with that of the full MCMC implementation of Heydenreich et al. (2022).