Probing dark energy with tomographic weak-lensing aperture mass statistics

We forecast and optimize the cosmological power of various weak-lensing aperture mass ($M_{\rm ap}$) map statistics for future cosmic shear surveys, including peaks, voids, and the full distribution of pixels (1D $M_{\rm ap}$). These alternative methods probe the non-Gaussian regime of the matter distribution, adding complementary cosmological information to the classical two-point estimators. Based on the SLICS and cosmo-SLICS $N$-body simulations, we build Euclid-like mocks to explore the $S_8 - \Omega_{\rm m} - w_0$ parameter space. We develop a new tomographic formalism which exploits the cross-information between redshift slices (cross-$M_{\rm ap}$) in addition to the information from individual slices (auto-$M_{\rm ap}$) probed in the standard approach. Our auto-$M_{\rm ap}$ forecast precision is in good agreement with the recent literature on weak-lensing peak statistics, and is improved by $\sim 50$% when including cross-$M_{\rm ap}$. It is further boosted by the use of 1D $M_{\rm ap}$ that outperforms all other estimators, including the shear two-point correlation function ($\gamma$-2PCF). When considering all tomographic terms, our uncertainty range on the structure growth parameter $S_8$ is enhanced by $\sim 45$% (almost twice better) when combining 1D $M_{\rm ap}$ and the $\gamma$-2PCF compared to the $\gamma$-2PCF alone. We additionally measure the first combined forecasts on the dark energy equation of state $w_0$, finding a factor of three reduction of the statistical error compared to the $\gamma$-2PCF alone. This demonstrates that the complementary cosmological information explored by non-Gaussian $M_{\rm ap}$ map statistics not only offers the potential to improve the constraints on the recent $\sigma_8$ - $\Omega_{\rm m}$ tension, but also constitutes an avenue to understand the accelerated expansion of our Universe.


Introduction
The coherent distortion between galaxy shapes due to the gravitational lensing by large-scale structure has emerged as one of the most powerful cosmological probes today. So-called cosmic shear analyses have started to unlock their potential thanks to large observational surveys of hundreds of square degrees at optical wavelength: e.g. CFHTLenS (Heymans et al. 2012), KiDS (de Jong et al. 2013), DES (Flaugher 2005), HSC (Aihara et al. 2018). In the meantime, the community is preparing for the next generation of cosmic shear surveys that will probe thousands of square degrees and extend to large redshifts: e.g. Euclid (Laureijs et al. 2011), Rubin Observatory (formerly LSST, Ivezić et al. 2019), Nancy Grace Roman Space Telescope (formerly WFIRST, Spergel et al. 2015).
In particular, weak-lensing peak statistics have proven to be one of the most powerful such new probes. Recent observation-(e.g. Liu et al. 2015a;Kacprzak et al. 2016;Martinet et al. 2018) and simulation-based (e.g. Zürcher et al. 2020) analyses show that peak statistics improve the constraints on S 8 = σ 8 (Ω m /0.3) 0.5 by 20% to 40% compared to the γ-2PCF alone. This particular combination of the matter density Ω m and clustering amplitude σ 8 runs across the main degeneracy direction obtained from cosmic shear constraints. It has notably been used to assess the potential tension between KiDS γ-2PCF cosmic shear and Planck CMB analyses (e.g. Hildebrandt et al. 2017;Planck Collaboration et al. 2020;Joudaki et al. 2020;Heymans et al. 2020). In this article we generalize the peak formalism to aperture-mass statistics, which provides additional constraining power on S 8 .
The aperture mass map (M ap map in the following) essentially is a map of convergence, convolved with a 'Mexican hat'like filter, which can be directly obtained from the estimated shear field (see Sect. 2). Other interesting M ap map-sampling methods exist beside peaks: e.g. voids (e.g. Gruen et al. 2016;Coulton et al. 2020) and the full distribution of pixels (aka the lensing PDF, e.g. Patton et al. 2017;Shirasaki et al. 2019;Liu & Madhavacheril 2019). Because they probe more massive structures, peaks are less sensitive to noise, and likely better at studying cosmological parameters sensitive to matter (Ω m , σ 8 ) and the sum of neutrino mass (Coulton et al. 2020). Voids, however, could offer a better sensitivity to dark energy. The optimal choice of lensing estimator is currently unknown and needs to be found.
In the context where Stage-IV weak lensing experiments are about to see first light, one of the key question that needs to be addressed is how much improvement can we expect from higherorder statistics when it comes to constraining the dark energy equation of state w 0 . In this paper, we present the first forecasts of this parameter from joint M ap and γ-2PCF statistics, and provide recommendations for applying tomography in a way that optimizes the total constraining power.
Indeed, while routinely used in γ-2PCF analyses, tomography has rarely been applied in mass map cosmological analyses (e.g. Yang et al. 2011;Martinet et al. 2015;Petri et al. 2016;Li et al. 2019;Coulton et al. 2020;Ajani et al. 2020;Zürcher et al. 2020), less so with a realistic redshift distribution (e.g. Martinet et al. 2015;Petri et al. 2016;Zürcher et al. 2020). One caveat of the current mass map tomography is that it only probes the information from individual redshift slices. This is similar to using only the auto-correlation terms in the case of the γ-2PCF and likely explains why the γ-2PCF benefits more from tomography than mass map estimators (Zürcher et al. 2020). In this article, we develop a new mass map tomographic approach to include the cross-information between redshift slices, and apply it to simulated lensing catalogs matched in depth to the future Euclid survey, with a broad redshift distribution between 0 < z < 3.
The aim of this article is to serve as a basis for the application of M ap non-Gaussian statistics to future cosmic shear surveys. In particular we review and optimize the statistical precision gain on S 8 and w 0 brought by different samplings of the M ap map, either using the full distribution of pixels, the peaks, or the voids. We also optimize the tomographic approach in the context of a joint analysis with standard two-point estimators. We finally test the impact of performing the M ap map computation in real or Fourier space, in presence of observational masks, and study the effect of parameter sampling explored with the N-body simulations.
The covariance matrix is built from the SLICS simulations (Harnois-Déraps et al. 2015) while the cosmology dependence is modeled with the cosmo-SLICS (Harnois-Déraps et al. 2019), which sample the parameters S 8 , Ω m , w 0 and the reduced Hubble parameter h at 25 points. This allows us to model the signals everywhere inside a broad parameter space and therefore to surpass the Fisher analyses often performed in the literature.
We first describe the mass map calculation that we use (M ap ) in Sect. 2, and investigate its sensitivity to observational masks in Appendix A. We then present in Sect. 3 the SLICS and cosmo-SLICS N-body simulation mocks. We detail the different data vectors (DVs) in Sect. 4 and our new tomographic approach in Sect. 5. We present the likelihood pipeline used to make our cosmological forecasts in Sect. 6, examining the cosmological parameter sampling in Appendix B. We forecast the constraints on S 8 and w 0 from M ap statistics in Sect. 7, focusing on the tomography improvement (Sect. 7.1) and on the sampling of M ap (Sect. 7.2), and finally combine it with the γ-2PCF (Sect. 8). We summarize our results in Sect. 9.

Aperture mass computation
The M ap statistic, first defined by Schneider (1996), presents several advantages over classical mass reconstructions that particularly fit the purpose of cosmological analyses from mass maps. We first describe the general expression for M ap in the shear and convergence space and discuss the pros and cons of this method (Sect. 2.1). We then develop a Fourier space approach that increases the speed of the reconstruction (Sect. 2.2).

Real space computation
The M ap at a position θ 0 consists in a convolution between the tangential ellipticity t (θ, θ 0 ) with respect to that point and a circularly symmetrical aperture filter Q(θ) (with θ = |θ|), where the sum is carried out over observed galaxies and normalized by the galaxy density n gal within the aperture. The tangential ellipticity is defined from the two components of the ellipticity and the polar angle φ(θ, θ 0 ) of the separation vector θ − θ 0 with respect to the center of the aperture θ 0 , A hat denotes spin-2 quantities. Assuming that the mean intrinsic ellipticity within the aperture averages to 0, Eq. (1) is an unbiased estimator of Eq. (3) in the weak-lensing regime (|γ|, |κ| 1) where the reduced shearĝ =γ/ (1 − κ) is equal to the shearγ: The tangential shear γ t (θ, θ 0 ) is computed by replacingˆ byγ in Eq.
(2), and is related to the intrinsic galaxy ellipticityˆ S througĥ = (ˆ S +ĝ)/(1 +ĝ * ˆ S ) (Seitz & Schneider 1997). Throughout this paper, we use the reduced shear, as this corresponds to the quantity observed in data. We also tested our whole pipeline on the shearγ instead and verified that it does not impact our cosmological forecasts, validating the use of M ap with reduced shear.
As demonstrated in Schneider & Bartelmann (1997), the aperture mass can equivalently be written as a convolution of the convergence κ(θ) and a compensated aperture filter U(θ): As U(θ) is a compensated filter, dθ θ U(θ) = 0, which means that the residual constant from integrating the shear into the convergence vanishes. This is not the case for classical mass reconstructions (e.g. Kaiser & Squires 1993) which usually resort to averaging the mass map to zero to lower the impact of that constant, referred to as mass sheet degeneracy. As the integration constant will vary from field to field and between the observational data and the modeled ones, it introduces random variations in the mass levels that could bias the cosmological constraints from mass map statistics. Dealing with the mass sheet degeneracy issue makes M ap a well suited method for cosmology analyses (e.g. peak statistics).
Another important advantage of M ap is that one can define a local noise as the dispersion of M ap assuming no shear, i.e. only due to galaxy intrinsic ellipticities: σ(M ap (θ 0 )) accounts for local variations of the galaxy density, arising from masks or different magnitude depths within a survey. From Eqs.
(1) & (6), one can derive the local signal-to-noise ratio (hereafter SNR), which will be used to measure the height of M ap pixels: In contrast to Martinet et al. (2018), we do not include galaxy weights and shear multiplicative biases as this article is primarily intended to assess the potential of M ap map cosmology for future cosmic shear surveys in terms of gain in precision. The impact of the shear measurement systematics, while being an important topic, is beyond the scope of this paper.
The formalism above does not assume any particular function for Q(θ) besides the circular symmetry. In this work, we use the Schirmer et al. (2007) filter which is a simpler form of an NFW (Navarro et al. 1997) halo mass profile, and is hence optimized to detect dark matter halos associated with small-scale structures in M ap maps, where θ ap is the size of the aperture and x c indicates a change of slope at the angular scale x c θ ap and is analogous to the r s parameter in the NFW profile. x c θ ap is referred to as the effective size of the aperture, and corresponds to the smoothing scale of the mass map in other traditional mass reconstruction (e.g. Kaiser & Squires 1993). The first term in Eq. (8) introduces a smooth transition to 0 in the core and edge of the profile. These cuts avoid the singularity in the aperture center, where γ t is not defined, and the contamination from the strong shear regime when centered on a massive halo. Following Hetterscheidt et al. (2005), we use x c = 0.15 as an optimal value for halo detection. Although a particular effective aperture size increases the SNR of halos with the same physical scale in the M ap map, it also acts as a smoothing for shape noise and depends on the galaxy density. We tested several effective scales and found that θ ap x c = 1.5 (corresponding to an aperture size θ ap = 10 ) maximizes the cosmological forecasts in the case of our mocks in all tomographic configurations, and therefore only present results for this aperture size. We do not explore the combination of M ap statistics measured from multiple smoothing scales as Martinet et al. (2018) showed that this combination of scales only marginally improves the cosmological constraints in the case of M ap map peaks.

Fourier space computation
Aperture masses are conveniently computed in Fourier space, where the convolution becomes a product: where a tilde denotes Fourier transformed quantities. Carrying out the convolution with FFT significantly reduces the computational time. In our case, for a 1024 × 1024 pixel map of 100 deg 2 with a galaxy density of 30 arcmin −2 and a filter scale of 10 arcmin, the Fourier space approach is ∼ 25 times faster than in real space, dropping from ∼ 50 to ∼ 2 minutes on a regular desktop computer. Although they are equivalent in theory, the real and Fourier space approaches differ in practice: the Fourier approach suffers from a loss of resolution and periodic boundary effects. The latter effect is suppressed by removing a stripe of width θ ap at the field boundary after the convolution (as is also done in the direct-space computation). While we use exact galaxy positions in real space, we have to bin them on a map before performing the Fourier transform, reducing the resolution to that of the final map.
Despite these limitations, our Fourier-space and real-space approaches agree very well in the absence of masks, producing visually identical maps with a null mean residual per pixel. Exploring the impact on the distribution of SNR values of the M ap map, we find that the Fourier approach introduces a loss of power for the high SNR, with 5% fewer pixels with SNR ≥ 3.0 compared to the real-space computation. This power leakage due to the discretization of galaxy positions in Fourier space depends on the sampling frequency: it is of 14% if we decrease the pixel resolution by a factor 2 and drops to only 2% if we increase the resolution by the same factor (i.e. reconstructing a 512 × 512 and 2048 × 2048 pixels map respectively). We however keep our fiducial resolution (1024 × 1024 pixels) to lower the computational time. We also verified that the Fourier reconstruction does not significantly bias the cosmological forecasts. We find almost no difference (below 0.3% on any parameter precision and bias) between the forecasts in Fourier and real space, validating the use of the Fourier approach.
In the M ap formalism, observational masks affect all pixels closer to the mask than the aperture size both in the direct and Fourier approaches. It is then mandatory to apply the same mask to the simulations used to build the cosmology dependence of M ap . In Appendix A, we investigate the possible impact of a representative Euclid mask on the cosmological forecasts. We find a decrease in the forecast precision consistent with the loss of area due to masking and no bias for both the direct-and Fourier-space map computations. In the rest of the paper we use the Fourierspace approach and do not apply any observational mask.

SLICS and cosmo-SLICS
The analysis presented in this work relies on two suites of numerical weak lensing simulations, the Scinet LIght Cones Simulations (Harnois-Déraps et al. 2015  particles inside a box of comoving side set to 505 h −1 Mpc. Multiple mass sheets were written to disk at regular comoving intervals, and were subsequently used to generate convergence (κ) and shear (γ 1/2 ) maps for a series of source planes (between 15 and 28, depending on the cosmology) from which the lensing quantities can be interpolated given a position and a redshift.
The SLICS were specifically designed for the estimation of covariance matrices: they consist of 928 fully independent Nbody ΛCDM runs in which the cosmological parameters are fixed to Ω m = 0.2905, σ 8 = 0.826, h = 0.6898, n s = 0.969 and Ω b = 0.0473, while the random seeds in the initial conditions are varied, allowing us to estimate covariance matrices with the 'brute force' ensemble approach. A single light-cone was constructed from each N-body run, ensuring a complete independence between the 928 light-cones.
The cosmo-SLICS were run with a similar but complementary N-body configuration: the random seeds are fixed while the cosmological parameters Ω m , S 8 , h, and w 0 are sampled at 25 points 1 . The parameters n s and Ω b were fixed to the value used in the SLICS. The sampling of this four-dimensional space follows a Latin hypercube, which maximizes the interpolation accuracy. Additionally, at each of these nodes, a pair of simulations were produced in which the sampling variance was highly suppressed from a careful selection of the initial conditions (see Harnois-Déraps et al. 2019), such that any measurement averaged over the pair nearly follows the ensemble mean. In total, ten pseudo-independent light-cones were constructed at every cosmo-SLICS cosmology (five per pair member) by randomly rotating the planes and shifting the origin 2 . A 26 th cosmo-SLICS simulation is run with the same set-up but with cosmological parameters fixed to the SLICS ΛCDM model, and is used to mock the observational data from which we try to estimate the cosmology.
The accuracy of the SLICS and cosmo-SLICS N-body simulations has been quantified with two-point correlation functions and power spectra in Harnois-Déraps et al. (2015) and Harnois-Déraps et al. (2019), and is basically within 2% of the Cosmic Emulator (Heitmann et al. 2014) up to k = 2.0 h Mpc −1 , with a gradual degradation of the accuracy at smaller scales. The covariance matrix measured from the SLICS is also in excellent agreement with the analytical predictions, at least for γ-2PCF and power spectra, and was shown in Harnois-Déraps et al. (2019) to contain more than 80% of the 'Super Sample Covariance'. This is a strong indicator of the high accuracy we achieve on the covariance matrix estimated with the SLICS.

Redshift distribution
From the shear and convergence planes described in the last section, we then construct a mock galaxy catalog for each lightcone. Because the shear and convergence depend both on the foreground mass and the source redshifts, it is mandatory to reproduce the redshift distribution of the target survey when building the mocks. In order to mimic a survey at the expected depth of the Euclid Wide Survey (Laureijs et al. 2011), we build a representative redshift distribution from the COSMOS 2015 catalog (Laigle et al. 2016) and use a galaxy density of 30 arcmin −2 (Laureijs et al. 2011). We first select COSMOS galaxies with an i -band magnitude lower than 24.5 to account for the expected limiting magnitude of the Euclid lensing sample in the V IS band, a large filter covering the r, i, and z bands. We then build the distribution of photometric redshifts in the range 0 < z < 3, with a large bin size of 0.1 in redshift in order to smooth the individual structures in the COSMOS field. We fit this redshift distribution with the Fu et al. (2008) function, We see in Fig. 1 that this function is able to reproduce the full redshift distribution including the high-z tail and smooths isolated peaks due to cosmic variance. This would not be the case for the classical Smail et al. (1994) function, n(z) ∝ (z/z 0 ) α exp [−(z/z 0 ) β ], which fails to capture the high-redshift tail above z ∼ 2. The parameters of the fit, that are used to generate the redshift distributions in our simulations are: A = 1.7865 arcmin −2 , a = 0.4710, b = 5.1843, c = 0.7259.

Galaxy positions and ellipticities
Galaxy 2D projected positions are randomly drawn from a uniform distribution. The correlation between the density of galaxies and matter is therefore not captured in our galaxy mocks.
Since the average value of the lensing estimator is unaffected by the source galaxy positioning, this choice does not bias our cosmological predictions. When using simulations in combination with observations though, the positions of galaxies in the simulations are chosen to mimic that of the observation, so as to reduce the impact of masks and shape noise. As a side effect, this positioning introduces wrong correlations between galaxy positions and matter overdensities in the simulations. This is often corrected for through the computation of a boost factor corresponding to the ratio of the radial number density profiles of galaxies around observed and simulated mass peaks (e.g. Kacprzak et al. 2016;Shan et al. 2018). Although we can neglect this effect in simulation-only analyses, a careful measurement of this bias will be necessary to unlock the full potential of M ap statistics in observational analyses. When assigning an intrinsic ellipticityˆ S to the mock galaxies, each component is drawn from a Gaussian distribution with zero mean and a dispersion σ = 0.26. This latter value corresponds to Hubble Space Telescope measurements for galaxies of magnitudes i ∼ 24.5 (Schrabback et al. 2018) which will constitute most of the Euclid lensing sample and was used previously to investigate the impact of faint blends on shear estimates (e.g. Euclid Collaboration: Martinet et al. 2019).
Galaxy positions and intrinsic ellipticities constitute the primary source of noise for the extraction of the cosmological signal. To prevent any noise effect on the model of the cosmological dependence of the mass map, we fix the positions, redshifts and intrinsic ellipticities of every galaxy in the cosmo-SLICS mocks for all the different cosmologies, in which only the cosmic shear is changed. These sources of noise are however included in the covariance matrix: each mass map is computed from an individual realization of the SLICS mocks with different galaxy positions and ellipticities. See Sect. 6 for more details.

Data vectors
In this section we define several estimators built from the M ap maps (Sect. 4.1) as well as the classical γ-2PCF (Sect. 4.2) which serves as our baseline comparison method.

Sampling the M ap distribution
Once M ap maps have been generated, there are several choices to sample the 1D distribution of M ap . In this article we investigate the use of distributions of M ap pixels with larger SNR value than their n neighboring pixels. In particular, n = 8 corresponds to maxima, and n = 0 to minima. Maximia, or peaks, trace the matter overdensities and in the high-SNR regime are a good representation of massive halos. Peaks therefore share some infor- The number count distribution of these estimators as a function of SNR constitute our data vectors -DVs for short. They are shown in Fig. 2 for a 100 deg 2 of mock data without tomographic decomposition. 1D M ap is displayed in red, peaks in green, and voids in blue. The curves correspond to the mean over 50 mock observations from the cosmo-SLICS fiducial cosmology (2 different initial conditions for the N-body run, 5 different ray-tracing through the light-cones, and 5 different shape noise realizations) and the error bars are estimated from the diagonal elements of the covariance matrix computed from the SLICS. As expected from their definition, peaks correspond to higher SNR pixels while the distribution of voids is centered on lower SNR. The dotted lines in Fig. 2 show the DVs measured from noise-only maps computed for the 5 shape noise realizations we use. We see that M ap maps are dominated by noise, with only a small fraction of the DVs carrying the cosmic shear information. However, the DVs extend to SNR values well above that of the noise-only DVs. These broader wings also imply that the DVs are below the noise near their maximum to preserve normalization.
In Fig. 3, we exhibit the cosmology dependence by comparing 1D M ap from Fig. 2 (in black) with that from the cosmo-SLICS simulations probing a range of cosmological parameters (blue to red curves). As for Fig. 2, we show the mean over 50 M ap maps of 100 deg 2 each, from which we subtracted the noise contribution to better visualize the cosmology dependence. The dependence on Ω m , S 8 , h, and w 0 are shown in different panels and colored from blue to red for increasing parameter values. We find a smooth gradient of 1D M ap with S 8 , and a somewhat A&A proofs: manuscript no. spe smaller variation with Ω m and w 0 , highlighting a strong dependence on these parameters. There however does not seem to be any particular dependence on h.
In Fig. 3 and all the following results we consider only the SNR ranges that can be accurately modeled with our emulator (see Sect. 6.3): −2.5 < SNR < 5.5 for peaks, −5 < SNR < 3 for voids, and −4 < SNR < 5 for 1D M ap , and use a bin size ∆SNR = 1. This corresponds to the largest bin size for which the cosmological forecasts are not degraded in our analysis. The high correlation between close bins in Fig. 3 however tends to show that the information is partly redundant. DV compression, e.g. though principal component analysis, would probably allow one to reduce the size of the M ap DVs to a few interesting numbers, although we do not explore this possibility in the present article.

Shear two-point correlation functions
To assess the complementary information from non-Gaussian estimators to traditional ones, we measure the γ-2PCF in the same mocks on which we compute the aperture mass maps. We use the classical ξ + /ξ − definition of γ-2PCFs (Schneider et al. 2002) which sums the ellipticity correlations over N pairs pairs of galaxies in bins of separation θ = |θ a − θ b |, where t and × correspond to the tangential and cross ellipticity components with respect to the separation vector θ a −θ b between the two galaxy positions θ a and θ b . As for M ap we do not include observational weights in the former equation.
Although ξ + /ξ − can be linked to the integrals of the matter power spectrum weighted by a 0 th /4 th order Bessel function to build a theoretical model of its cosmology dependence (see e.g. Kilbinger 2015, for a review), we use the cosmo-SLICS mocks to model the dependence of the γ-2PCF on cosmology. This ensures a fair comparison between the different probes and also allows for a full combination of DVs when measuring joint constraints. We note that fixing or varying galaxy positions in the cosmo-SLICS mocks does not impact the γ-2PCF, which are insensitive to source positions.
We measure ξ ± with the tree code Athena . We use an opening angle of 1.7 degrees to merge trees and speed up the computation. This approximation mainly affects the computation at large scales which are not included in our analysis given the size of our simulation mocks (10 × 10 deg 2 ), and we verified that ξ ± do not change for a smaller opening angle of 1.13 degrees (as advised in the Athena readme). We compute the correlations for 10 separation bins logarithmically spaced between 0.1 and 300 . We note that the use of N-body simulations allows us to include scales below 0.5 ; these probe the deeply non-linear regime of structure formation, for which the dependence on cosmology is difficult to model correctly with theoretical approaches or fit functions. We nevertheless need to cut large scales due to the size of our simulation mocks, and therefore discard the two last bins of ξ + . We also remove the two first bins of ξ − where the information content is low and which correspond to scales where the matter power spectrum is not fully resolved in the N-body simulations. We are left with 8 bins centered on (0.15 , 0.33 , 0.74 , 1.65 , 3.67 , 8.17 , 18.20 , 40.54 ) for ξ + , and (0.74 , 1. 65 , 3.67 , 8.17 , 18.20 , 40.54 , 90.27 , 201.03 ) for ξ − . These scales are similar to that used in current cosmic shear survey (e.g. in KiDS: Hildebrandt et al. 2017) with the addition of the small non-linear scales for ξ + as they are expected to be used in Euclid (also directly calibrated from N-body simulations, as in Euclid Collaboration: Knabenhans et al. 2019).

Tomography
Tomography is key to improving the constraints on the dark energy equation of state w 0 . It is traditionally implemented by slicing the redshift distribution of sources and computing the estimators for the different slices. In principle, a large number of slices allows one to better sample the information of the third dimension. This however decreases the density of galaxies per slice such that measurements become noisier and it is not yet clear what is the ideal trade off between these two effects.
We therefore test various tomographic configurations from 1 slice (i.e. no tomography) up to 10 slices, the target for 2-point statistics in Euclid (Laureijs et al. 2011). The boundaries of each tomographic bin are chosen to have the same galaxy density in each slice. In this way, shape noise affects all slices similarly, which allows us to fairly compare the signal from different redshift bins. All redshifts are assumed to be perfectly known in this article as we do not intend to assess the impact of photometric redshift uncertainties (e.g. Euclid Collaboration: Desprez et al. 2020) or biases in the mean of the redshift distribution (e.g. Joudaki et al. 2020).
This tomographic approach to mass maps is however suboptimal as it does not exploit the correlations between different redshift bins. This loss of information likely explains why the γ-2PCF is better improved by tomography than mass map estimators (Zürcher et al. 2020). Indeed, the former includes correlation functions in between slices (cross terms) in addition to within slices (auto terms). Following the example of γ-2PCF we develop a new approach in which we explicitly compute cross terms for M ap statistics. This is performed by computing M ap maps for all combinations of redshift slices (M ap (z i ∪ z j ∪ ...)) in addition to the auto-M ap terms (M ap (z i )). For example, with 3 tomographic slices, we use the cross terms M ap (z 1 ∪ z 2 ), M ap (z 1 ∪ z 3 ), M ap (z 2 ∪ z 3 ), M ap (z 1 ∪ z 2 ∪ z 3 ), in addition to the usual auto terms M ap (z 1 ), M ap (z 2 ), M ap (z 3 ). The DVs measured in these extra maps cannot be reconstructed from those measured in the M ap maps of the individual slices and contain information about the relative positions of structures in the different redshift slices. This is illustrated in Fig. 4 where we show a 1 deg 2 M ap SNR map reconstructed for each individual redshift slice and for a few combinations of redshifts in a 5 tomographic bin set up. A striking point is the repetition of high-SNR patterns across slices, which likely originate from the same massive structure. Indeed, a foreground massive halo will distort galaxy shapes in several background source planes and higher redshift slices contain more information than lower ones. In particular the lowest tomographic slice (z 1 < 0.4676) is almost consistent with a noise-only map. As our DVs are made of pixels counts, the information about the position of structures is lost. The addition of combined redshift bins can compensate (at least partially) for this loss as we include the information about the relative positions of structures probed by the combined redshift planes. We also note that the SNR is larger in the non-tomographic case, which probes 5 times as many galaxies as a single slice. This z all map exhibits several structures that likely originate from chance superposition of distinct smaller halos (Yang et al. 2011) that are often seen in the tomographic slices, and appear as a single large feature in the z all map. The disentangling of these structures with the tomographic approach improves the recovery of the cosmological information embedded in the 3D matter distribution.
For γ-2PCF, the tomographic DV is defined as the concatenation of the auto-correlations (ξ ± in one slice) and crosscorrelations (ξ ± between pairs of slices). For M ap statistics, we concatenate the DVs measured in every combination of redshift slices, including combinations of more than two slices. The cross-terms significantly increase the size of the DV and with the 928 realizations of the SLICS simulations used to compute the covariance matrix, the latter becomes too noisy for large number of tomographic slices. For 10 tomographic bins, the γ-2PCF DV would contain 880 bins requiring several thousands simulations for an accurate estimate of the covariance matrix (and its inverse). While we study the previously used tomography of auto-M ap up to 10 slices, we do not include tomography above 5 slices for the γ-2PCF and for the cross-M ap . This corresponds to DVs of respectively 240 and 279 elements for a total of 519 elements in the combined γ-2PCF and cross-M ap analysis. For larger number of slices, the estimation of the covariance starts to become noisy, especially when combining γ-2PCF and M ap .

Computing cosmological forecasts
Forecasts are obtained by finding the cosmological parameters that maximize the likelihood (see Sect. 6.4) of the model given an observation. Noise is accounted for through the covariance matrix computed from the SLICS simulations (DVs labelled x s ) in Sect. 6.1. The model that characterizes the dependence of a given estimator on cosmology is described in Sect. 6.2 and is emulated for any cosmology from the cosmo-SLICS simulations (DVs labelled x m ) in Sect. 6.3. The mock observation (DVs labelled x) consists in one of the cosmo-SLICS model that was not used to calibrate the cosmological dependence. Fig. 4. Aperture mass maps of 1 deg 2 obtained from the full redshift distribution (0 < z all < 3), for 5 continuous redshift slices (0 < z 1 < 0.4676, 0.4676 < z 2 < 0.7194, 0.7194 < z 3 < 0.9625, 0.9625 < z 4 < 1.3319, 1.3319 < z 5 < 3) and for a few cross-combinations of redshift slices (z 2 ∪ z 3 , z 2 ∪ z 4 , z 2 ∪ z 3 ∪ z 4 ). Massive foreground halos are traced in multiple redshift planes while the noise fluctuates randomly. The cross-M ap terms allow us to retain the information about correlations of structures in between slices.

Noise
The noise is accounted for through the covariance matrix, that estimates the correlation between the data vectors x s,i for N s different mock realizations i, (x s,i (π 0 ) −x s (π 0 )) (x s,i (π 0 ) −x s (π 0 )) T . (12) x s denotes the average value of the data vector over its different realizations. We neglect the dependence of the covariance matrix on cosmology and estimate it from the ΛCDM SLICS simulations at cosmology π 0 . This has been shown to be the correct approach in a Gaussian likelihood framework, since varying the cosmology in the covariance can introduce unphysical information (Carron 2013).
Each of the N s = 928 mock realizations is fully independent in terms of all three noise components: sample variance, shape noise, and position noise. We find that the shape noise, introduced by a particular realization of galaxy intrinsic ellipticities, is of the same order as the sample variance. The fact that galaxy positions are random gives rise to 'position noise', which adds to the shape noise, but only increases the covariance by a few percents at high SNR and up to ∼ 20% at SNR of 0. The position of sources has a weaker impact on large SNR values that correspond to massive halos for which the shear is larger.
The fact that we estimate our DV covariance matrix from a finite suite of simulations results in residual noise, which propagates into an error on the constraints on the cosmological parameter (Taylor & Joachimi 2014) and into an additional error on the parameters (Dodelson & Schneider 2013). The accuracy ν of the parameter covariance depends on the number of realizations used in its computation (N s ), on the size of the data vector (N d ) and on the number of cosmological parameters inferred (N p ), via N s N d + N p ν −1 + 2ν −2 (from Eq. 13 of Taylor & Joachimi 2014). With N s = 928, N p = 4 and N d varying between 8 bins in the simplest case to 519 bins for the combination of 1D Map and γ-2PCF with a 5-slice tomography including all cross terms, we find that the accuracy on our errors on π = (S 8 , Ω m , w 0 , h) is at worst 7.5%. The additional error prescribed by Dodelson & Schneider (2013) is obtained by rescaling the covariance by a factor (1 + N d /N s ), a 0.4% correction that we neglect in this work.
We show the covariance matrix normalized by its diagonal elements (the correlation matrix) for 1D M ap in Fig. 5. It displays high correlations between close SNR bins, and presents three different regimes as already seen from the DV excess in Fig. 3. Figure 6 is similar to Fig. 5 but for a tomographic analysis with 5 redshift slices and combining with γ-2PCF. The lower left quadrant of the matrix displays ξ ± ordered from low to high redshifts including both auto-and cross-correlations: ξ ± (z 1 ), ξ ± (z 1 , z 2 ), ..., ξ ± (z 2 ), ξ ± (z 2 , z 3 ), ..., ξ ± (z 5 ). The upper right quadrant shows 1D M ap first for auto-M ap , and then by increasing the combination size of cross-M ap : M ap (z 1 ), ..., M ap (z 5 ), M ap (z 1 ∪ z 2 ), ..., M ap (z 4 ∪ z 5 ), ..., M ap (z 1 ∪ z 2 ∪ z 3 ), ..., M ap (z 1 ∪ z 2 ∪ z 3 ∪ z 4 ), ..., M ap (z all ). First looking at M ap , we note that the non-tomographic correlation pattern is seen in every tomographic auto-and crossinformation term. We additionally see correlations in between redshift slices and for combined slices that fade for larger distance between slices. As already noted from the M ap maps (see Fig. 4) this is probably due to massive halos being observed from multiple background source planes. The high correlations between cross-M ap terms also suggest that one could find a more compact representation of the tomographic DV. Now focusing on the combination of the two estimators, we find somewhat lower correlations than for individual DVs. This shows that γ-2PCF and M ap maps indeed probe partially independent information such that their combination should result in a gain of precision in the forecasts.

Model
The model of the cosmology dependence of each data vector is built from the cosmo-SLICS simulations. For each of the 25 nodes in the parameter space, we compute M ap maps and measure the distribution of SNR of each estimator and then interpolate the DV at any cosmology.
In contrast to the computation of the covariance matrix where all sources of noise must be accounted for, we need the model to be as insensitive to noise as possible to ensure the model is accurate in the parameter space close to that of the observational data.
Sample variance is reduced at the level of the simulations and of the ray-tracing. As detailed in Harnois-Déraps et al. (2019), the cosmo-SLICS N-body runs were produced in pairs for which the sample variance is already highly suppressed, allowing us to rapidly approach the ensemble mean by averaging our estimated signal over these pairs. Additionally, 10 pseudoindependent light-cones are generated for each of these simulation pairs at every cosmology, using the same random rotations and shifting. Averaging over these further suppresses the sam- Fig. 6. Correlation matrix computed from the combined γ-2PCF and 1D M ap DVs with 5 tomographic slices and including cross-terms. ξ ± is ordered by increasing redshift from ξ ± (z 1 ), ξ ± (z 1 , z 2 ), ..., to ξ ± (z 5 ). 1D M ap is ordered by increasing size of the combination of slices: M ap (z i ), ..., M ap (z i ∪ z j ), ..., M ap (z i ∪ z j ∪ z k ), ..., M ap (z i ∪ z j ∪ z k ∪ z l ), ..., M ap (z all ). The lower correlations between the two probes indicate a potential gain of cosmological information from their combination.
pling variance associated with the observer's position. To minimize shape noise, galaxy intrinsic ellipticities and positions are kept identical for every cosmology and light-cone, and we generate 5 different shape noise realizations to lower the noise due to a particular choice of ellipticities. In practice we only change the orientation of galaxies and keep the same ellipticity amplitude, as those should be determined from the observations. This results in a total of 50 mocks of 100 deg 2 per cosmology: (2 different initial conditions for the N-body simulations) × (5 different ray-tracings) × (5 different shape noise realizations). These numbers are chosen such that adding extra sample variance or shape noise realizations does not affect the cosmological forecasts any more.

Emulator
We then emulate the DVs measured from the 25 nodes at any point in the parameter space through radial basis functions. We use the scipy.interpolate.rbf python module that was shown to perform well in interpolating weak-lensing peak distributions (e.g. Liu et al. 2015a;Martinet et al. 2018). We improve from these studies by using a cubic function instead of multiquadric, as we find it to give more accurate results in the case of our mocks. Each bin of SNR of the DV is interpolated independently.
A classical method to verify the accuracy of the interpolation is to remove one node from the cosmological parameter space and to compare the emulation with the measured DV at this node. The proposed verification would however overestimate the errors by measuring them at more than the largest possible distance to a point in the latin hypercube parameter space. We choose instead to evaluate the results of the emulation at the only data point which was not part of the latin hypercube: our observation mock. In Fig. 7, we show the relative error due to the interpolation of 1D M ap at this point (solid green curve). The DV is cut at low and high SNR to retain an accuracy of better than 5% in the interpolation. We also display the relative difference between the Fig. 7. Testing the emulation of 1D M ap . The solid green curve corresponds to the relative excess of the interpolated DV with respect to the measured DV. This error on the interpolation remains below 5% (dashed black lines) and small in comparison to the relative excess interpolated for a variation of 2% (doted green curve) and 5% (dashed green curve) of S 8 , and to the 25 DVs used to build the emulator (blue to red curves). data vector for the 25 values of S 8 and the observation (blue to red curves) to compare the cosmology dependence of the model with the error on the interpolation. As can be seen, the error on the interpolation is low compared to the variation due to cosmology. The interpolated DV around the input S 8 value (doted and dashed green curves) fits well in the cosmology dependence of the model. Similar figures can be produced for the other 3 cosmological parameters only with a weaker cosmology dependence as already seen from Fig. 3. The stability of the emulator is tested in more details in Appendix B where we also vary the number of points used in the initial parameter space.
Using the emulator we then generate DVs for a grid of cosmological parameters spanning the full parameter space with 40 points for each parameter. This might introduce a bias on the measured parameter values of up to the step of the interpolation: 0.011 for Ω m , 0.007 for S 8 , 0.037 for w 0 , and 0.005 for h. We find this limitation to be subdominant in the case of the forecasts for 100 deg 2 , and close to the 2% uncertainties from the cosmo-SLICS simulations. We also apply a Gaussian smoothing to the interpolation results to improve the rendering of the likelihood contours generated for this sparse ensemble of points. In further studies it would be interesting to develop a refining mesh approach or a Markov Chain Monte Carlo sampler to improve the sampling of the parameter space.

Likelihood
The model built from the cosmo-SLICS simulations allows us to estimate the likelihood: the probability p(x|π) to measure a DV x for a given set of cosmological parameters π = (S 8 , w 0 , Ω m , h). This likelihood is related to the quantity of interest, the posterior likelihood or the probability of a cosmology given the observed DV p(π|x) through Bayes' theorem, where p(π) is a uniform prior within the parameter range probed by the cosmo-SLICS simulations, and p(x) is a normalization factor. We use a Student-t likelihood, which is a variation to the Gaussian likelihood which includes a marginalization over the true covariance matrix (Sellentin & Heavens 2016). This marginalization accounts for the uncertainty in the inverse of the covariance matrix (the precision matrix) due to the finite number of simulations used in its estimation that would otherwise artificially improve the precision of the forecasts (Hartlap et al. 2007). When neglecting the dependence of the covariance matrix on cosmology, this likelihood reads (Sellentin & Heavens 2016): N s = 928 is the number of SLICS simulations used to compute the covariance matrix. We note that for large N s , Eq. (14) converges to the Gaussian likelihood. For the likelihood to be correct, we need to ensure that the DV behaves as a multivariate Gaussian distribution. We therefore discard SNR bins with less than 10 M ap pixel counts, although this criterion is almost always already satisfied by the limitation we impose on the DV interpolation. χ 2 (x, π) is a measure of the deviation of a measured DV x from a DV x m modeled at cosmology π, and accounting for the measurement errors through the precision matrix Σ −1 at a fixed cosmology π 0 . We note that this Bayesian approach is superior to the Fisher formalism sometimes used in the literature (e.g. Martinet et al. 2015), as it does not assume a linear dependence of the DV on the cosmological parameters considered here, an assumption which most likely is invalid in the case of non-Gaussian statistics. Our approach is also closer to an observation-based analysis and allows us to estimate not only the precision on the cosmological parameters but also potential biases.

Forecasts
The cosmological forecasts are obtained by finding the cosmology π best that maximizes the posterior likelihood. Although we follow a full 4-dimensional analysis we display 2D and 1D forecasts by marginalizing over the other parameters. The 1σ and 2σ contours are calculated as the contours enclosing respectively 68% and 95% of the marginalized likelihood.
Throughout this article, reported forecast values correspond to the 1σ limit in the 1D marginalized likelihood. Likewise, the biases are computed as the difference between the best estimate in the 1D likelihood and the true parameter value in input of the observation mock. Percentage values are in percentage of the true parameter value. We study both the expected precision and bias on S 8 , w 0 , and Ω m , but do not present results for h. Indeed, our precision forecasts on h are unreliable because the likelihood extends up to the prior limit of that parameter. This advocates for using a wider prior on h when defining the parameter space to run the N-body simulations.
Forecasts are computed for a 100 deg 2 survey at Euclid depth. Given the few percent accuracy of the cosmo-SLICS simulations, it is not possible to study potential biases on the recovered parameter values for the full 15 000 deg 2 Euclid survey area, since the statistical precision would exceed that of the simulations, and therefore we also refrain from computing accurate forecasts in the latter case. Although our analysis is well suited to explore the potential of M ap map statistics as a cosmological probe, simulations with improved accuracy will be needed for future observations.

M ap statistics forecasts
We concentrate first on results obtained for the M ap statistics and discuss those from the γ-2PCF in the next section.

1D M ap tomographic forecasts
We show in Fig. 8 and in the upper part of Table 1 the forecasts for a 100 deg 2 Euclid-like survey for 1D M ap without tomography (blue contours and curves) and with a tomographic analysis with 5 redshift slices including both auto-and cross-terms (green contours and curves). We find a good sensitivity of 1D M ap on S 8 and Ω m . The sensitivity to w 0 is lower but improves more drastically than other parameters when including tomography, as seen especially in the 1D likelihoods. All constraints become tighter in the tomographic case, highlighting the gain of information from probing along the redshift direction. The 1σ credible intervals are reduced by including 5-slice tomography by 32% on S 8 , 49% on w 0 , and 46% on Ω m for 1D M ap . Using only the auto-M ap terms would reduce this gain to 13%, 34%, and 20% on S 8 , w 0 , and Ω m respectively. These results show that cross-M ap terms contain significant additional information to the auto-terms, notably improving w 0 forecasts by an extra ∼ 50%.
The likelihood contours appear noisy because of the interpolation of the model DVs that loses accuracy when further away from one of the 25 cosmo-SLICS simulations, represented as grey crosses in Fig. 8. We also note a small bias when comparing the best estimates to the true values but it remains within the 1σ error bars, except for Ω m when including tomography. Although the precision of the forecasts is greatly enhanced by the tomographic approach, the biases are less affected by it. Only the bias on w 0 is significantly reduced with tomography because of the disappearance of the double peak in the 1D likelihood (see Fig. 8). The fact that the bias remains constant suggests that it is driven by the simulation parameter space Table 1. Cosmological predictions for 100 deg 2 Euclid-like mock without tomography and with a 5-slice tomography set-up, including cross-M ap terms, for voids, peaks, 1D M ap , γ-2PCF, and for the combination of M ap map estimators with the γ-2PCF. The 'δ' values refer to the 1σ precision forecasts, while '∆' measure the bias between the best fit value and the input. Numbers in parenthesis show the results in percentage of the input value. rather than the DV measurements. We explore this further in Appendix B, where we vary the number of points in the parameter space used to perform the model interpolation. We find that removing up to 5 points out of the 25 we use in the fiducial analysis only slightly decreases the forecast precision while the bias increases dramatically. These results confirm the high impact of the number of training nodes on this bias, and demonstrate that a densely sampled parameter space will be necessary when applying M ap statistics to observations.
We now investigate how the precision gain varies with the considered number of tomographic slices, and try to identify the optimal number of slices for M ap map cosmological analyses. We show this variation of the forecast precision in Fig. 9 for S 8 , w 0 , and Ω m for all DVs. The dotted lines correspond to the previous tomography set-up including only the M ap maps of individual redshift slices (auto-M ap ), and the solid lines to our refined set-up with combinations of redshift slices (auto-and cross-M ap ). For the former we explore the dependence on up to 10 redshift slices, and for the latter up to 5 because of the larger DV size and the corresponding restriction of the calculation of the covariance from SLICS (see Sect. 6.1).
When considering only auto-M ap , the gain brought by tomography rapidly converges to a fixed value for all three parameters probed. The largest gain is observed when going from 1 (no tomography) to 2 slices, and does not evolve any more above ∼ 5 slices. This is in good agreement with results from Yuan et al. (2019) who found for a Euclid-like weak-lensing peak analysis that the constraints on Ω m and σ 8 are only marginally improved when using 8 tomographic slices compared to 4. This is no longer the case when we additionally include cross-M ap terms. We find that the forecasts are greatly enhanced compared to the previous implementation of mass map tomography, and that the gain keeps increasing up to the largest number of tomographic bins that we could test. In the auto-M ap approach, most of the tomographic gain is likely brought by separating the unlensed foreground galaxies from the background galaxies that carry the cosmic shear information. This is supported by the observation that the gain is mainly concentrated on the low number of slices in this case. In contrast, the cross-terms bring new information by probing correlations between structures across the different redshift slices. This extra information also removes the double-peak in the w 0 likelihood, that explains the strange behavior of the auto-M ap constraints between 2 and 4 tomographic slices: a small shift in the position of the likelihood maximum with respect to this double peak can significantly reduce the error bars by excluding one of the peaks from the ±1σ uncertainty interval.
Our new tomographic approach shows great potential and will allow us to exploit the 3D cosmic shear information with improved constraints at least up to 5 redshift slices. However, identifying the number of z-slices at which the information content saturates will require a larger number of N-body simulations in order to improve the accuracy of the model dependence of M ap on cosmology and of the covariance matrix. Finally, we note that the number of z-bins that can realistically be employed might also be set by our ability to recover unbiased mean source redshifts in these slices, an effect that is not studied in this article.

Comparison of 1D M ap with peaks and voids
In Fig. 9 and Table 1, we also present the forecasts for peaks and voids measured in the M ap maps.
1D M ap performs better than both peaks and voids regardless of the considered cosmological parameter and tomographic set-up. Peaks show a higher sensitivity to S 8 and Ω m than voids and both estimators perform as well when focusing on w 0 . This can be explained by the fact that peaks probe matter overdensities, while dark energy also significantly affects the properties of voids. 1D M ap has a superior constraining power because it includes both peaks and voids. The likelihood is however noisier for 1D M ap as many M ap pixels contain no signal, but this does not degrade the performance of this estimator. The forecast precision for 1D M ap is 7%, 31%, and 23% better than for peaks on S 8 , w 0 , and Ω m respectively without tomography, and 4%, 26%, and 21% when including tomography with 5 redshift slices. We also note that the relative gain of information from tomography is larger for voids than for peaks or 1D M ap .
With regards to biases in the recovered cosmology, we see that they are quite similar for the different estimators, another hint that they are primarily due to the sampling of the parameter space by the N-body simulations, as studied in Appendix B. Fig. 9. Dependence of the 1σ error on the number of tomographic bins for S 8 (top), w 0 (middle), and Ω m (bottom). Constraints are represented in blue for voids, green for peaks, and red for the full 1D M ap distribution. Dotted and solid lines respectively correspond to the previous tomography set-up (auto-M ap only) and to the improved tomography (auto-and cross-M ap ). Black dots represent the γ-2PCF (both auto-and cross-correlations) in configurations with 1 and 5 tomographic slices, and color dots the combination between M ap map DV of the same color and γ-2PCF including all auto-and cross-terms.
Comparing weak-lensing peak constraints to current surveys, we find a dramatic improvement due to the larger galaxy density and to the fact that higher-redshift galaxies carry more lensing information. With 100 deg 2 at Euclid depth without including tomography, we forecast a precision of 4.1% on S 8 which is almost twice as good as the currently best observational constraints from peak statistics (7.1% using the first 450 deg 2 of KiDS; Martinet et al. 2018). This highlights the potential of forthcoming Stage IV cosmic shear surveys.
Finally, we also experimented with other DVs: the distribution of M ap pixels higher than their n neighbors, for n varying between 0 and 8. n = 4 is particularly interesting, as it includes (but is not restricted to) the saddle points that can be linked to the presence of filaments (Codis et al. 2015). This estimator behaves very similarly to 1D M ap but has a lower constraining power. It however exceeds the forecast precision of peaks and voids. The other n values present a smooth degradation of the constraints from n = 4 to the peaks (n = 8) and n = 4 to the voids (n = 0).
Overall, 1D M ap performs better than any other estimator and should be used to obtain the best possible cosmological constraints from M ap map statistics, unless it is found in the future that it is more severely affected by systematic uncertainties than the other probes. We first compare the statistical power of M ap DVs and γ-2PCF before combining them. This comparison is based on Fig. 9 and the values of the top part of Table 1. When no tomography is included, the γ-2PCF is outperformed by all M ap DVs for all parameters, except by voids on Ω m . The best estimator, 1D M ap , presents constraints better than the γ-2PCF by 24%, 48%, and 47% on S 8 , w 0 , and Ω m respectively, highlighting the greater power of non-Gaussian estimators. We also note that the γ-2PCF presents less noisy likelihoods (see Fig. 10), as this estimator is less impacted by shape noise than M ap maps.

Combining
When applying tomography with 5 redshift slices and including both auto-and cross-terms, 1D M ap still outperforms the γ-2PCF for w 0 and Ω m , but presents a similar constraining power than the γ-2PCF on S 8 . In this set-up, the constraints are still 7%, 47% and 44% better with 1D M ap than γ-2PCF on S 8 , w 0 and Ω m . The smaller difference on S 8 with the tomographic γ-2PCF could partially be explained by the favorable definition of this parameter. S 8 corresponds to the particular value of α = 0.5 in Σ 8 = σ 8 (Ω m /0.3) α which is chosen to maximise the precision of the γ-2PCF on that parameter. Since M ap DVs present different degeneracies in the σ 8 -Ω m plane than the γ-2PCF (in particular large SNR peaks; e.g. Shan et al. 2018), they are probably more sensitive to a Σ 8 defined with a different α parameter.
For w 0 and Ω m , we find the 1D M ap forecasts to be superior to the γ-2PCF by roughly the same amount without and with tomography. This is in contrast with the previous M ap tomography set-up (doted lines in Fig. 9) which is closer to and sometimes even outperformed by the γ-2PCF for larger number of redshift slices. These results confirm the hypothesis of Zürcher et al. (2020) that the γ-2PCF benefits more from tomography than mass map estimators in their analysis because they include cross-correlations between redshift slices in the γ-2PCF but not in their mass reconstruction. In contrast, our addition of the cross-M ap terms allows us to better extract the tomographic information similarly to the γ-2PCF cross-correlations.
We note here that our constraints on the γ-2PCF are representative of the Euclid survey set up, and can therefore be safely compared to other estimators in this analysis. Assuming a precision dependence proportional to the square root of the number of galaxies considered, the 18.6% constraints on w 0 in the 5 tomography set up extrapolates to a 1.5% forecasts on this parameter for the 15 000 deg 2 final area, in good agreement with the expected 1% precision for γ-2PCF when including 10 tomographic slices (Laureijs et al. 2011). We also verified that ξ ± measured in our simulation-based mocks is consistent with theoretical predictions from Takahashi et al. (2012) using the nicaea software (Kilbinger et al. 2009). This model is based on halofit (Smith et al. 2003) with a calibration on high-resolution N-body simulations to account for the non-linear power spectrum contribution on small scales. The Takahashi et al. (2012) model is within the 1σ uncertainty of our DV (computed from the covariance matrix) on every scale except on the smallest bin where it is within 2σ.

Combination of 1D M ap and γ-2PCF
The combination of M ap estimators with the γ-2PCF is very powerful. It is shown as colored dots in Fig. 9. One striking point is the smaller difference between the three M ap DVs forecasts when combined with the γ-2PCF. This suggests that a large part of the non-Gaussian complementary information to the γ-2PCF is contained in all M ap DVs.
We also show in Figs. 10 & 11 the 2D and 1D likelihoods for 1D M ap (in green), γ-2PCF (in black), and their combination (in purple), with 1 and 5 tomographic slices. We see a large gain on all parameters when combining estimators. Without tomography, the combined forecasts are largely dominated by the 1D M ap constraints. With tomography, both 1D M ap and γ-2PCF contribute to the better forecasts, especially for w 0 .
The quantitative forecasts are displayed in the lower part of Table 1 and can be compared with the upper part values for the individual estimators. We focus here on the combination of the γ-2PCF with 1D M ap , which is the best M ap estimator. Without tomography the combined forecasts for S 8 , w 0 , and Ω m are respectively 34%, 56%, and 55% better than that of the γ-2PCF alone. The gain is lower when compared to 1D M ap alone: 13% on S 8 , 16% on w 0 , and 15% on Ω m . When applying a 5 bin tomography with all cross terms, the gain is of 46%, 68%, and 57% on S 8 , w 0 , and Ω m respectively, compared to the γ-2PCF alone, and 42%, 39%, and 24% compared to 1D M ap alone. These improved forecasts confirm the good complementarity between non-Gaussian and traditional 2-point estimators.

Comparison with literature on weak-lensing peak statistics
We now compare our forecasts with recent results in the literature of weak-lensing peak statistics. Without tomography we find an improvement of 34% on S 8 when combining peaks and γ-2PCF compared to the γ-2PCF alone. This is somewhat higher than the value of 20% found in the case of the KiDS data in Martinet et al. (2018), and could be explained by the higher galaxy density of our Stage IV mocks which decreases the impact of shape noise in the M ap maps.
When including tomography with auto-and cross-M ap , we find an improvement of 36% on S 8 and 51% on Ω m with peaks Fig. 11. Same as Fig. 10 for a tomographic analysis with 5 redshift slices, including auto-and cross-M ap terms. and γ-2PCF compared to γ-2PCF alone. In order to compare with the literature, we also compute the combined forecasts for auto-M ap only. In that case we find an improvement of 21% on S 8 and 38% on Ω m .
The latter constraints are in very good agreement with recent results from Li et al. (2019) who found an improvement of 32% on Ω m when combining convergence peaks and power spectrum in their tomographic analysis with LSST mocks. These authors also computed forecasts on the amplitude of the primordial matter power spectrum A s (similar to σ 8 ) and the sum of the neutrinos that we do not explore in our analysis.
Our forecasts on S 8 are also close to the 25% improvement from the combination of peaks and power spectrum compared to power spectrum alone performed by Zürcher et al. (2020) on DES mocks including tomography.
The similar forecast improvements on S 8 and Ω m between our analysis with auto-M ap terms and the recent literature comforts the robustness of our results. The fact that our constraints are also systematically better when including the cross-M ap terms highlights again the superiority of this new tomographic approach for cosmological analyses with mass maps.
For the first time, we forecast constraints on w 0 through a combination of tomographic M ap methods and γ-2PCF. We find that the combination of peaks and γ-2PCF improves the constraints on this parameter by 52% compared to the γ-2PCF alone in our tomography set-up (and 37% when including only auto-M ap ). The gain from M ap -map peaks is lower than for 1D M ap which also probes the information from other structures. In addition to improve our chances of solving the recent σ 8 -Ω m ten-sion by increasing the precision on S 8 , we now know that weaklensing peak statistics (and more so 1D M ap ) also probe significant complementary information to the γ-2PCF on w 0 , and are therefore a powerful tool to help understand the nature of dark energy from cosmic shear surveys.

Conclusion
In this article, we optimized the method to apply M ap statistics to future cosmic shear surveys (e.g. Euclid, LSST, WFIRST). We created Euclid-like galaxy mocks of 100 deg 2 based on the SLICS N-body simulations to build the covariance matrix, and on 25 different cosmologies, the cosmo-SLICS simulations, to infer the model dependence of different estimators on cosmological parameters S 8 , w 0 , and Ω m . We computed aperture-mass maps from the shear instead of traditional mass-map reconstructions that require to calculate the convergence first, and verified that realistic observational masks are not biasing the forecasts. We developed a new tomographic approach which includes cross-information between redshift slices in the range 0 < z < 3. We emulated the cosmology dependence of three non-Gaussian statistics and of the γ-2PCF. We forecast their precision and bias on cosmology, maximizing the likelihood between the model, accounting for sample variance and shape noise in the covariance matrix. We then combined these non-Gaussian estimators with the γ-2PCF, including the correlations between the two probes, to assess the statistical gain. The main results of our analysis are the following: • We built different DVs from the M ap maps: the full distribution of pixels (1D M ap ), the local maxima (peaks), and minima (voids). We found that 1D M ap outperforms other estimators, followed by peaks, and finally voids, both with and without tomography. These estimators are also superior to the γ-2PCF: in particular, the uncertainties on w 0 from the 1D M ap are smaller than for the γ-2PCF by a factor of ∼ 2, irrespective of the tomographic configuration.
• The explicit computation of the cross-M ap terms from the combinations of redshift slices increases the tomographic forecast precision by an extra ∼ 50% compared to the previous tomographic approach which only relies on individual redshift slices (auto-M ap terms in our framework). Only these latter terms have been used in the literature up to now, explaining why authors find, in contrast to us, that non-Gaussian mass-map estimators profit less from tomography than the γ-2PCF.
• The combination of M ap maps with γ-2PCF significantly improves the forecasts with respect both to γ-2PCF alone and M ap map alone. This important gain is explained by the different part of the matter distribution that are probed: mostly Gaussian information from the growth of structure on large scale and at higher redshift from 2-point estimators, and the non-Gaussian part from small scale and lower redshift evolution of structures for mass maps. Uncertainties in S 8 are reduced by ∼ 35% and ∼ 45% (almost twice better), respectively without and with tomography, compared to the γ-2PCF alone. These numbers are larger than in the recent literature due to our new tomographic approach but are in good agreement when we exclude the cross-M ap terms of our analysis.
• For the first time, we also estimated joint tomographic forecasts on w 0 finding a decrease of the uncertainty by about a factor of three when combining 1D M ap and γ-2PCF compared to the γ-2PCF alone. Further studies are required to better understand how dark energy affects the non-Gaussian part of the matter distribution which is probed by M ap map statistics. The high sensitivity of the halo-mass function to the growth of structures at small scales, or the fact that the M ap statistics studied in this paper are probing a collection of high-order moments could constitute possible ways of explaining this large gain of information compared to 2-point statistics.
These results highlight the feasibility and potential of M ap statistics (including weak-lensing peaks and voids) as complementary probes to the γ-2PCF in cosmic shear surveys. The large gain on both S 8 and w 0 offers a greater chance to solve some of the most puzzling current questions in cosmology, e.g. the recent σ 8 -Ω m tension and the nature of the accelerated expansion of our Universe. Relying solely on γ-2PCF to answer these questions would deprive us of a large statistical power already present in our observations, and which is both larger and complementary to that probed by standard 2-point estimators inherited from CMB analyses. While the observational systematics of γ-2PCF are very well understood, combining both estimators in observational data however requires to improve our modeling of the biases inherent to mass map estimators. This will be done by further developing N-body simulations specifically tailored to study the impact of the model emulation, preferentially in a Markov Chain Monte Carlo framework, and by carefully investigating the impact of classical cosmic shear systematics on M ap map estimators: shear multiplicative bias, uncertainty on the mean redshift, baryon feedback, and intrinsic alignments.