Issue 
A&A
Volume 646, February 2021



Article Number  A62  
Number of page(s)  16  
Section  Cosmology (including clusters of galaxies)  
DOI  https://doi.org/10.1051/00046361/202039679  
Published online  05 February 2021 
Probing dark energy with tomographic weaklensing aperture mass statistics
^{1}
AixMarseille Univ, CNRS, CNES, LAM, Marseille, France
email: nicolas.martinet@lam.fr
^{2}
Institute for Astronomy, University of Edinburgh, Royal Observatory, Blackford Hill, Edinburgh EH9 3HJ, UK
^{3}
Astrophysics Research Institute, Liverpool John Moores University, 146 Brownlow Hill, Liverpool L3 5RF, UK
^{4}
ArgelanderInstitut für Astronomie, Universität Bonn, Auf dem Hügel 71, 53121 Bonn, Germany
Received:
14
October
2020
Accepted:
20
November
2020
We forecast and optimize the cosmological power of various weaklensing aperture mass (M_{ap}) map statistics for future cosmic shear surveys, including peaks, voids, and the full distribution of pixels (1D M_{ap}). These alternative methods probe the nonGaussian regime of the matter distribution, adding complementary cosmological information to the classical twopoint estimators. Based on the SLICS and cosmoSLICS Nbody simulations, we build Euclidlike mocks to explore the S_{8} − Ω_{m} − w_{0} parameter space. We develop a new tomographic formalism that exploits the crossinformation between redshift slices (crossM_{ap}) in addition to the information from individual slices (autoM_{ap}) probed in the standard approach. Our autoM_{ap} forecast precision is in good agreement with the recent literature on weaklensing peak statistics and is improved by ∼50% when including crossM_{ap}. It is further boosted by the use of 1D M_{ap} that outperforms all other estimators, including the shear twopoint correlation function (γ2PCF). When considering all tomographic terms, our uncertainty range on the structure growth parameter S_{8} is enhanced by ∼45% (almost twice better) when combining 1D M_{ap} and the γ2PCF compared to the γ2PCF alone. We additionally measure the first combined forecasts on the dark energy equation of state w_{0}, finding a factor of three reduction in the statistical error compared to the γ2PCF alone. This demonstrates that the complementary cosmological information explored by nonGaussian M_{ap} map statistics not only offers the potential to improve the constraints on the recent σ_{8}–Ω_{m} tension, but also constitutes an avenue to understanding the accelerated expansion of our Universe.
Key words: gravitational lensing: weak / cosmology: observations / dark energy / largescale structure of Universe / surveys
© N. Martinet et al. 2021
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1. Introduction
The coherent distortion between galaxy shapes due to the gravitational lensing by largescale structure has emerged as one of the most powerful cosmological probes today. Socalled cosmic shear analyses have started to unlock their potential thanks to large observational surveys of hundreds of square degrees at optical wavelength: for example, CFHTLenS (Heymans et al. 2012), KiDS (de Jong et al. 2013), DES (Flaugher 2005), and 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: for example, Euclid (Laureijs et al. 2011), Vera C. Rubin Observatory (formerly LSST, Ivezić et al. 2019), and Nancy Grace Roman Space Telescope (formerly WFIRST, Spergel et al. 2015).
So far, these surveys have mostly relied on twopoint estimators for their cosmological analysis, for example the shear twopoint correlation functions (γ2PCF, e.g., Kilbinger et al. 2013; Troxel et al. 2018; Hikage et al. 2019; Asgari et al. 2021). These estimators are inherited from cosmic microwave background (CMB) analyses, which probe the matter distribution of the early Universe. However, cosmic shear probes the recent Universe, where the matter distribution is more complex due to the nonlinear accretion of structures that creates nonGaussian features in the matter field on small scales (e.g., Codis et al. 2015; Barthelemy et al. 2020). Twopoint statistics fail to capture this nonGaussian information and thus yield an incomplete description of the matter distribution at low redshift. To close this gap, the community has recently started to explore nonGaussian cosmic shear estimators: for example weaklensing peaks (e.g., Kruse & Schneider 1999, 2000; Dietrich & Hartlap 2010; Kratochvil et al. 2010; Fan et al. 2010; Yang et al. 2011; Maturi et al. 2011; Hamana et al. 2012; Hilbert et al. 2012; Marian et al. 2012, 2013; Shan et al. 2014, 2018; Lin & Kilbinger 2015; Martinet et al. 2015, 2018; Liu et al. 2015a,b; Kacprzak et al. 2016; Petri et al. 2016; Zorrilla Matilla et al. 2016; Giocoli et al. 2018; Peel et al. 2018; Davies et al. 2019; Fong et al. 2019; Li et al. 2019; Weiss et al. 2019; Yuan et al. 2019; Coulton et al. 2020; Ajani et al. 2020; Zürcher et al. 2021), Minkowski functionals (e.g., Kratochvil et al. 2012; Petri et al. 2015; Vicinanza et al. 2019; Parroni et al. 2020; Zürcher et al. 2021), higherorder moments (e.g., Van Waerbeke et al. 2013; Petri et al. 2015; Peel et al. 2018; Vicinanza et al. 2018; Chang et al. 2018; Gatti et al. 2020), threepoint statistics (e.g., Schneider & Lombardi 2003; Takada & Jain 2003, 2004; Semboloni et al. 2011; Fu et al. 2014), density split statistics (e.g., Friedrich et al. 2018; Gruen et al. 2018; Burger et al. 2020), persistent homology (e.g., Heydenreich et al. 2021), scattering transform (e.g., Cheng et al. 2020), and machine learning (e.g., Merten et al. 2019; Ribli et al. 2019; Peel et al. 2019; Shirasaki et al. 2019; Fluri et al. 2019). Although these new statistics have shown a great potential, they have not yet been fully generalized to data analyses, mostly because they need a large number of computationally expensive Nbody simulations to calibrate their dependence on cosmology, while this dependence is accurately predicted by theoretical models in the case of twopoint estimators.
In particular, weaklensing 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 simulationbased (e.g., Zürcher et al. 2021) 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 VI 2020; Joudaki et al. 2020; Heymans et al. 2021). In this article, we generalize the peak formalism to aperturemass statistics, which provides additional constraining power on S_{8}.
The aperture mass map (M_{ap} map in the following) is essentially 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} mapsampling methods exist beside peaks: for example voids (e.g., Gruen et al. 2016; Coulton et al. 2020) and the full distribution of pixels, also known as the lensing probability distribution function (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 are likely better for 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 determined.
In the context where StageIV weak lensing experiments are about to see first light, one of the key questions 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. 2021), less so with a realistic redshift distribution (e.g., Martinet et al. 2015; Petri et al. 2016; Zürcher et al. 2021). 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 autocorrelation 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. 2021). In this article, we develop a new mass map tomographic approach to include the crossinformation 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} nonGaussian 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 twopoint 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 Nbody simulations.
The covariance matrix is built from the SLICS simulations (HarnoisDéraps et al. 2015) while the cosmology dependence is modeled with the cosmoSLICS (HarnoisDé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 cosmoSLICS Nbody 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.
2. 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).
2.1. 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 spin2 quantities. Assuming that the mean intrinsic ellipticity within the aperture averages to 0, Eq. (1) is an unbiased estimator of Eq. (3) in the weaklensing regime () where the reduced shear 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 through (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 to 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, that is only due to galaxy intrinsic ellipticities:
The σ(M_{ap}(θ_{0})) accounts for local variations of the galaxy density that arise from masks or different magnitude depths within a survey.
From Eqs. (1) & (6), one can derive the local signaltonoise ratio (S/N), 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 smallscale 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. The 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 S/N 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 six different effective scales (0.5′,1.0′,1.5′,2.0′,2.5′,3′) 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. This optimal scale is closely followed by the 1′ and 2′ scales in term of forecast precision, in good agreement with other stage IV forecasts (for example Li et al. 2019 found an optimal smoothing scale of 2′ for their LSST mocks). The fact that we find the same optimal θ_{ap} in the different tomographic cases suggests that the choice of aperture size is dominated by the scale of the matched structures rather than the noise smoothing thanks to the high galaxy density. 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.
2.2. 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 min 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 directspace 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 Fourierspace and realspace 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 S/N values of the M_{ap} map, we find that the Fourier approach introduces a loss of power for the high S/N, with 5% fewer pixels with S/N ≥ 3.0 compared to the realspace 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 two 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 Fourierspace map computations. In the rest of the paper we use the Fourierspace approach and do not apply any observational mask.
3. Building mock catalogs
3.1. SLICS and cosmoSLICS
The analysis presented in this work relies on two suites of numerical weak lensing simulations, the Scinet LIght Cones Simulations (SLICS; HarnoisDéraps et al. 2015, 2018) and the cosmoSLICS (HarnoisDéraps et al. 2019), which are respectively used for estimating the covariance matrix and modeling the signal of our different statistics. Both are based on a series of 100 deg^{2} lightcones extracted from Nbody simulations with the multiple plane technique described in the above references. The underlying gravityonly calculations were carried out by CUBEP^{3}M (HarnoisDéraps et al. 2013) and evolve 1536^{3} 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 lightcone was constructed from each Nbody run, ensuring a complete independence between the 928 lightcones.
The cosmoSLICS were run with a similar but complementary Nbody 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 4D space follows a Latin hypercube, which maximizes the interpolation accuracy. Additionally, at each of these nodes, a pair of simulations was produced, in which the sampling variance was highly suppressed from a careful selection of the initial conditions (see HarnoisDéraps et al. 2019) such that any measurement averaged over the pair nearly follows the ensemble mean. In total, ten pseudoindependent lightcones were constructed at every cosmoSLICS cosmology (five per pair member) by randomly rotating the planes and shifting the origin^{2}. A 26^{th} cosmoSLICS simulation was run with the same setup 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 cosmoSLICS Nbody simulations was quantified with twopoint correlation functions and power spectra in HarnoisDéraps et al. (2015, 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 HarnoisDé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.
3.2. Redshift distribution
From the shear and convergence planes described in the last section, we then constructed 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 built a representative redshift distribution from the COSMOS 2015 catalog (Laigle et al. 2016) and used a galaxy density of 30 arcmin^{−2} (Laureijs et al. 2011). We first selected 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 VIS band, a large filter covering the r, i, and z bands. We then built 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 highz tail, and smooths isolated peaks that are 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 highredshift 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, and c = 0.7259.
Fig. 1.
Redshift distribution normalized to 30 gal.arcmin^{−2}. The red histogram corresponds to the COSMOS2015 photometric redshift distribution after a cut at i′≤24.5, the green curve to the fit of a Smail et al. (1994) function to the COSMOS2015 histogram, and the blue curve to a Fu et al. (2008) fit. In this paper we use the Fu et al. (2008) fit that captures the highredshift tail in contrast to the other function. 
3.3. 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. However, when using simulations in combination with observations, 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 that corresponds 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 simulationonly 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 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 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 cosmoSLICS 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.
4. 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.
4.1. 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 that have larger S/N values than their n neighboring pixels. In particular, n = 8 corresponds to maxima, and n = 0 to minima. Maxima, or peaks, trace the matter overdensities and in the highS/N regime are a good representation of massive halos. Peaks therefore share some information with galaxy cluster counts and both probes aim at exploring the nonGaussian tail of the matter distribution. Peaks, however, present the advantages of being insensitive to the classical selection function and massobservable relation issues inherent to cluster count studies (e.g., Sartoris et al. 2016; Schrabback et al. 2018; Dietrich et al. 2019; McClintock et al. 2019). Minima probe the voids of the matter distribution. Alternatively, one can use the full distribution of pixel values (hereafter 1D M_{ap} or lensing PDF).
The number count distribution of these estimators as a function of S/N 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 cosmoSLICS fiducial cosmology (2 different initial conditions for the Nbody run, five different raytracing through the lightcones, and five 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 S/N pixels while the distribution of voids is centered on lower S/N. The dotted lines in Fig. 2 show the DVs measured from noiseonly maps computed for the five 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 S/N values well above that of the noiseonly DVs. These broader wings also imply that the DVs are below the noise near their maximum to preserve normalization.
Fig. 2.
S/N distribution (or data vectors) of three different M_{ap} samplings: voids (blue), peaks (green), and full 1D distribution (red). Dotted lines correspond to the S/N distributions measured in noiseonly maps. Numbers are quoted for a 1024 × 1024 pixels map of 100 deg^{2} with a filter size θ_{ap} = 10′. 
In Fig. 3, we exhibit the cosmology dependence by comparing 1D M_{ap} from Fig. 2 (in black) with that from the cosmoSLICS 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 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.
Fig. 3.
Variation of the excess number of M_{ap} pixels over noise with cosmology. The black curve corresponds to the excess 1D M_{ap} in the observation mock, with error bars from the diagonal of the covariance matrix. Colorcoded are variations in Ω_{m} (top left), S_{8} (top right), h (bottom left), and w_{0} (bottom right). A smooth color gradient indicates a correlation between the data vector and the probed cosmological parameter. No tomography has been performed here. 
In Fig. 3 and all the following results we consider only the S/N ranges that can be accurately modeled with our emulator (see Sect. 6.3): −2.5 < S/N < 5.5 for peaks, −5 < S/N < 3 for voids, and −4 < S/N < 5 for 1D M_{ap}, and use a bin size ΔS/N = 1. Figure 2 suggests that bins of larger S/N could be used in the case of peaks and 1D M_{ap} when no tomography is applied. We find however that these bins only mildly improve the forecast precision (by less than 5% and 2% for peaks and 1D M_{ap} respectively) as the associated pixel counts are particularly low while they introduce biases due to the large uncertainty of the emulator model in this S/N range. The chosen bin size corresponds to the largest width 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, for example through 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.
4.2. Shear twopoint correlation functions
To assess the complementary information from nonGaussian 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 0th/4th order Bessel function to build a theoretical model of its cosmology dependence (see e.g., Kilbinger 2015, for a review), we use the cosmoSLICS 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 cosmoSLICS mocks does not impact the γ2PCF, which are insensitive to source positions.
We measure ξ_{±} with the tree code Athena (Kilbinger et al. 2014). 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 Nbody simulations allows us to include scales below 0.5′; these probe the deeply nonlinear 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 Nbody 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 surveys (e.g., in KiDS: Hildebrandt et al. 2017) with the addition of the small nonlinear scales for ξ_{+} as they are expected to be used in Euclid (also directly calibrated from Nbody simulations, as in Euclid Collaboration 2019).
5. 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 ten slices, the target for twopoint 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 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. 2021). 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 autoM_{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} S/N map reconstructed for each individual redshift slice and for a few combinations of redshifts in a five tomographic bin setup. A striking point is the repetition of highS/N 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 higherredshift slices contain more information than lower ones. In particular the lowest tomographic slice (z_{1} < 0.4676) is almost consistent with a noiseonly 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 S/N is larger in the nontomographic case, which probes five 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.
Fig. 4.
Aperture mass maps of 1 deg^{2} obtained from the full redshift distribution (0 < z_{all} < 3), for five 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 crosscombinations 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 crossM_{ap} terms allow us to retain the information about correlations of structures in between slices. 
For γ2PCF, the tomographic DV is defined as the concatenation of the autocorrelations (ξ_{±} 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 crossterms 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 ten tomographic bins, the γ2PCF DV would contain 880 bins requiring several thousand simulations for an accurate estimate of the covariance matrix (and its inverse). While we study the previously used tomography of autoM_{ap} up to ten slices, we do not include tomography above five slices for the γ2PCF and for the crossM_{ap}. This corresponds to DVs of, respectively, 240 and 279 elements for a total of 519 elements in the combined γ2PCF and crossM_{ap} analysis. For larger number of slices, the estimation of the covariance starts to become noisy, especially when combining γ2PCF and M_{ap}.
6. 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 labeled 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 cosmoSLICS simulations (DVs labeled x_{m}) in Sect. 6.3. The mock observation (DVs labeled x) consists in one of the cosmoSLICS model that was not used to calibrate the cosmological dependence.
6.1. 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,
The 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 S/N and up to ∼20% at S/N of 0. The position of sources has a weaker impact on large S/N values that correspond to massive halos for which the shear is larger.
Using a finite suite of simulations results in uncertainties in the covariance matrix which propagates into an error on the constraints on the cosmological parameters (Hartlap et al. 2007; Dodelson & Schneider 2013; Taylor & Joachimi 2014). This loss of information can be quantified by employing a Studentt distribution to describe the distribution of values in each S/N bins (Sellentin & Heavens 2017). This distribution is the result of the central limit theorem for a finite ensemble of independent realizations, and converges to the Gaussian distribution for an infinite number of samples. The systematically lost information in the variance of each cosmological parameter marginalized over all other parameters can be computed from the number of realizations used in the computation of the covariance matrix (N_{s}), the size of the data vector (N_{d}) and the number of cosmological parameters inferred (N_{p}), via Eq. (42) of Sellentin & Heavens (2017). 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 5slice tomography including all cross terms, we find that the accuracy on our errors on each individual parameter is better than 1%.
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 S/N 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 crosscorrelations: ξ_{±}(z_{1}), ξ_{±}(z_{1}, z_{2}),…, ξ_{±}(z_{2}), ξ_{±}(z_{2}, z_{3}),…,ξ_{±}(z_{5}). The upper right quadrant shows 1D M_{ap} first for autoM_{ap}, and then by increasing the combination size of crossM_{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 nontomographic 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 crossM_{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.
Fig. 5.
Correlation matrix computed from 1D M_{ap} in 928 simulations with independent noise realizations (sample variance, shape noise, and position noise). 
Fig. 6.
Correlation matrix computed from the combined γ2PCF and 1D M_{ap} DVs with five tomographic slices and including crossterms. ξ_{±} 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. 
6.2. Model
The model of the cosmology dependence of each data vector is built from the cosmoSLICS simulations. For each of the 25 nodes in the parameter space, we compute M_{ap} maps and measure the distribution of S/N 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 raytracing. As detailed in HarnoisDéraps et al. (2019), the cosmoSLICS Nbody 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, ten pseudoindependent lightcones are generated for each of these simulation pairs at every cosmology, using the same random rotations and shifting. Averaging over these further suppresses the sampling variance associated with the observer’s position. To minimize shape noise, galaxy intrinsic ellipticities and positions are kept identical for every cosmology and lightcone, and we generate five 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: (two different initial conditions for the Nbody simulations) × (five different raytracings) × (five 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.
6.3. 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 weaklensing 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 S/N 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 that 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 S/N to retain an accuracy of better than 5% in the interpolation. We also display the relative difference between the 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.
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). 
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 cosmoSLICS 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.
6.4. Likelihood
The model built from the cosmoSLICS 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 cosmoSLICS simulations, and p(x) is a normalization factor.
We use a Studentt likelihood, which is a variation to the Gaussian likelihood that accounts for the uncertainty in the inverse of the covariance matrix due to the finite number of simulations used in its estimation (Sellentin & Heavens 2016) that would otherwise artificially improve the precision of the forecasts. When neglecting the dependence of the covariance matrix on cosmology, this likelihood reads (Sellentin & Heavens 2016):
The N_{s} = 928 is the number of SLICS simulations used to compute the covariance matrix. We verified that every bin of our DVs are compatible with a multivariate Gaussian distribution by comparing their distribution across the SLICS to a Gaussian model drawn from the mean and dispersion of the same bin values in the SLICS. This shows that we have a sufficiently large ensemble of independent realizations that we could use a Gaussian likelihood, but we prefer to apply the Studentt which is more accurate for finite samples (it converges to the Gaussian likelihood for N_{s} → ∞). χ^{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 that is most likely invalid in the case of nonGaussian statistics. Our approach is also closer to an observationbased analysis and allows us to estimate not only the precision on the cosmological parameters but also potential biases.
6.5. Forecasts
The cosmological forecasts are obtained by finding the cosmology π_{best} that maximizes the posterior likelihood. Although we follow a full 4D 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 Nbody simulations.
Forecasts are computed for a 100 deg^{2} survey at Euclid depth. Given the few percent accuracy of the cosmoSLICS 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.
7. 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.
7.1. 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}Euclidlike survey for 1D M_{ap} without tomography (blue contours and curves) and with a tomographic analysis with 5 redshift slices including both auto and crossterms (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 5slice tomography by 32% on S_{8}, 49% on w_{0}, and 46% on Ω_{m} for 1D M_{ap}. Using only the autoM_{ap} terms would reduce this gain to 13%, 34%, and 20% on S_{8}, w_{0}, and Ω_{m}, respectively. These results show that crossM_{ap} terms contain significant additional information to the autoterms, notably improving w_{0} forecasts by an extra ∼50%. The reported values assume the same S/N range for the nontomographic and tomographic cases. We note that with a larger area one might be able to also accurately emulate these statistics for larger S/N bins in the nontomographic case, such that the gain brought by the tomographic approach could be slightly lower.
Fig. 8.
Forecast of cosmological parameters from 1D M_{ap} in a 100 deg^{2} survey at Euclid depth. Marginalized 2D (1 and 2σ contours) and 1D (full likelihood) constraints are displayed in blue without tomography and in green for a tomographic analysis with 5 redshift slices and including auto and crossM_{ap} terms between redshift slices. Dashed lines correspond to the 1σ constraints in the 1D marginalized likelihood. The red crosses and lines indicate the input values while the blue and green crosses indicate the best estimate in the 2D constraints. Gray crosses correspond to parameters of the 25 cosmoSLICS simulations that are used to estimate the cosmology dependence of 1D M_{ap}. 
Cosmological predictions for 100 deg^{2}Euclidlike mock without tomography and with a fiveslice tomography setup, including crossM_{ap} terms, for voids, peaks, 1D M_{ap}, γ2PCF, and for the combination of M_{ap} map estimators with the γ2PCF.
The likelihood contours appear noisy because of the interpolation of the model DVs that loses accuracy when further away from one of the 25 cosmoSLICS simulations, represented as gray 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 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 setup including only the M_{ap} maps of individual redshift slices (autoM_{ap}), and the solid lines to our refined setup with combinations of redshift slices (auto and crossM_{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).
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 setup (autoM_{ap} only) and to the improved tomography (auto and crossM_{ap}). Black dots represent the γ2PCF (both auto and crosscorrelations) 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 crossterms. 
When considering only autoM_{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 Euclidlike weaklensing 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 crossM_{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 autoM_{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 crossterms bring new information by probing correlations between structures across the different redshift slices. This extra information also removes the doublepeak in the w_{0} likelihood, that explains the strange behavior of the autoM_{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 five redshift slices. However, identifying the number of zslices at which the information content saturates will require a larger number of Nbody 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 zbins 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.
7.2. 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 setup. 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 five 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 Nbody simulations, as studied in Appendix B.
Comparing weaklensing peak constraints to current surveys, we find a dramatic improvement due to the larger galaxy density and to the fact that higherredshift 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. It should then 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.
8. Combining M_{ap} statistics and shear twopoint correlation functions
8.1. Comparison of 1D M_{ap} with γ2PCF
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 nonGaussian 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.
Fig. 10.
Forecast of cosmological parameters from 1D M_{ap}, γ2PCF and their combination in a 100 deg^{2}Euclidlike survey without tomography. Marginalized 2D (1 and 2σ contours) and 1D (full likelihood) constraints are displayed in black for the γ2PCF, green for 1D M_{ap}, and purple for their combination. The red crosses and lines indicate the truth. Gray crosses correspond to the 25 simulation points from which the model is built. 
When applying tomography with 5 redshift slices and including both auto and crossterms, 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 setup, 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 maximize 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 S/N 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 setup (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. (2021) that the γ2PCF benefits more from tomography than mass map estimators in their analysis because they include crosscorrelations between redshift slices in the γ2PCF but not in their mass reconstruction. In contrast, our addition of the crossM_{ap} terms allows us to better extract the tomographic information similarly to the γ2PCF crosscorrelations.
We note here that our constraints on the γ2PCF are representative of the Euclid survey setup, 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 five tomography setup 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 ten tomographic slices (Laureijs et al. 2011). We also verified that ξ_{±} measured in our simulationbased 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 highresolution Nbody simulations to account for the nonlinear 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σ.
8.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 nonGaussian complementary information to the γ2PCF is contained in all M_{ap} DVs.
We also show in Figs. 10 and 11 the 2D and 1D likelihoods for 1D M_{ap} (in green), γ2PCF (in black), and their combination (in purple), with one and five 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}.
Fig. 11.
Same as Fig. 10 for a tomographic analysis with five redshift slices, including auto and crossM_{ap} terms. 
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 nonGaussian and traditional twopoint estimators.
8.3. Comparison with literature on weaklensing peak statistics
We now compare our forecasts with recent results in the literature of weaklensing 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 crossM_{ap}, we find an improvement of 36% on S_{8} and 51% on Ω_{m} with peaks and γ2PCF compared to γ2PCF alone. In order to compare with the literature, we also compute the combined forecasts for autoM_{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. (2021) on DES mocks including tomography.
The similar forecast improvements on S_{8} and Ω_{m} between our analysis with autoM_{ap} terms and the recent literature comforts the robustness of our results. The fact that our constraints are also systematically better when including the crossM_{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 setup (and 37% when including only autoM_{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} tension 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.
9. 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 Euclidlike galaxy mocks of 100 deg^{2} based on the SLICS Nbody simulations to build the covariance matrix, and on 25 different cosmologies, the cosmoSLICS simulations, to infer the model dependence of different estimators on cosmological parameters S_{8}, w_{0}, and Ω_{m}. We computed aperturemass maps from the shear instead of traditional massmap reconstructions that require the calculation of the convergence first, and verified that realistic observational masks are not biasing the forecasts. We developed a new tomographic approach that includes crossinformation between redshift slices in the range 0 < z < 3. We emulated the cosmology dependence of three nonGaussian 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 nonGaussian 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 crossM_{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 (autoM_{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 nonGaussian massmap 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 parts of the matter distribution that are probed: mostly Gaussian information from the growth of structure on large scales and at higher redshift from twopoint estimators, and the nonGaussian part from small scales and the 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 crossM_{ap} terms of our analysis.

For the first time, we also estimated joint tomographic forecasts on w_{0} finding a decrease in 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 nonGaussian part of the matter distribution that is probed by M_{ap} map statistics. The high sensitivity of the halomass 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 highorder moments could constitute possible ways of explaining this large gain of information compared to twopoint statistics.
These results highlight the feasibility and potential of M_{ap} statistics (including weaklensing 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, such as 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 twopoint estimators inherited from CMB analyses. While the observational systematics of γ2PCF are very well understood, combining both estimators in observational data requires us to improve our modeling of the biases inherent to mass map estimators. This will be done by further developing Nbody 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.
Acknowledgments
We thank Tiago Castro, Carlo Giocoli, the anonymous referee, and the many KiDS collaborators for useful discussions. NM acknowledges support from a fellowship of the Centre National d’Etudes Spatiales (CNES). JHD acknowledges support from an STFC Ernest Rutherford Fellowship (project reference ST/S004858/1). Computations for the Nbody simulations were enabled by Compute Ontario (www.computeontario.ca), Westgrid (www.westgrid.ca) and Compute Canada (www.computecanada.ca). The SLICS numerical simulations can be found at http://slics.roe.ac.uk/, while the cosmoSLICS can be made available upon request. This work has been carried out thanks to the support of the OCEVU Labex (ANR11LABX0060) and of the Excellence Initiative of AixMarseille University – A*MIDEX, part of the French “Investissements d’Avenir” program.
References
 Aihara, H., Arimoto, N., Armstrong, R., et al. 2018, PASJ, 70, S4 [NASA ADS] [Google Scholar]
 Ajani, V., Peel, A., Pettorino, V., et al. 2020, Phys. Rev. D, 102, 103531 [CrossRef] [Google Scholar]
 Asgari, M., Lin, C.A., Joachimi, B., et al. 2021, A&A, 645, A104 [CrossRef] [EDP Sciences] [Google Scholar]
 Barthelemy, A., Codis, S., Uhlemann, C., Bernardeau, F., & Gavazzi, R. 2020, MNRAS, 492, 3420 [Google Scholar]
 Burger, P., Schneider, P., Demchenko, V., et al. 2020, A&A, 642, A161 [EDP Sciences] [Google Scholar]
 Carron, J. 2013, A&A, 551, A88 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Chang, C., Pujol, A., Mawdsley, B., et al. 2018, MNRAS, 475, 3165 [CrossRef] [Google Scholar]
 Cheng, S., Ting, Y.S., Ménard, B., & Bruna, J. 2020, MNRAS, 499, 5902 [CrossRef] [Google Scholar]
 Codis, S., Pichon, C., & Pogosyan, D. 2015, MNRAS, 452, 3369 [Google Scholar]
 Coulton, W. R., Liu, J., McCarthy, I. G., & Osato, K. 2020, MNRAS, 495, 2531 [CrossRef] [Google Scholar]
 Coupon, J., Czakon, N., Bosch, J., et al. 2018, PASJ, 70, S7 [Google Scholar]
 Davies, C. T., Cautun, M., & Li, B. 2019, MNRAS, 488, 5833 [Google Scholar]
 de Jong, J. T. A., Verdoes Kleijn, G. A., Kuijken, K. H., & Valentijn, E. A. 2013, Exp. Astron., 35, 25 [Google Scholar]
 Dietrich, J. P., & Hartlap, J. 2010, MNRAS, 402, 1049 [NASA ADS] [CrossRef] [Google Scholar]
 Dietrich, J. P., Bocquet, S., Schrabback, T., et al. 2019, MNRAS, 483, 2871 [NASA ADS] [CrossRef] [Google Scholar]
 Dodelson, S., & Schneider, M. D. 2013, Phys. Rev. D, 88, 063537 [NASA ADS] [CrossRef] [Google Scholar]
 Euclid Collaboration (Knabenhans, M., et al.) 2019, MNRAS, 484, 5509 [NASA ADS] [CrossRef] [Google Scholar]
 Euclid Collaboration (Martinet, N., et al.) 2019, A&A, 627, A59 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Euclid Collaboration (Desprez, G., et al.) 2020, A&A, 644, A31 [EDP Sciences] [Google Scholar]
 Fan, Z., Shan, H., & Liu, J. 2010, ApJ, 719, 1408 [NASA ADS] [CrossRef] [Google Scholar]
 Flaugher, B. 2005, Int. J. Mod. Phys. A, 20, 3121 [Google Scholar]
 Fluri, J., Kacprzak, T., Lucchi, A., et al. 2019, Phys. Rev. D, 100, 063514 [NASA ADS] [CrossRef] [Google Scholar]
 Fong, M., Choi, M., Catlett, V., et al. 2019, MNRAS, 488, 3340 [Google Scholar]
 Friedrich, O., Gruen, D., DeRose, J., et al. 2018, Phys. Rev. D, 98, 023508 [CrossRef] [Google Scholar]
 Fu, L., Semboloni, E., Hoekstra, H., et al. 2008, A&A, 479, 9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fu, L., Kilbinger, M., Erben, T., et al. 2014, MNRAS, 441, 2725 [NASA ADS] [CrossRef] [Google Scholar]
 Gatti, M., Chang, C., Friedrich, O., et al. 2020, MNRAS, 498, 4060 [CrossRef] [Google Scholar]
 Giocoli, C., Moscardini, L., Baldi, M., Meneghetti, M., & Metcalf, R. B. 2018, MNRAS, 478, 5436 [NASA ADS] [CrossRef] [Google Scholar]
 Gruen, D., Friedrich, O., Amara, A., et al. 2016, MNRAS, 455, 3367 [NASA ADS] [CrossRef] [Google Scholar]
 Gruen, D., Friedrich, O., Krause, E., et al. 2018, Phys. Rev. D, 98, 023507 [NASA ADS] [CrossRef] [Google Scholar]
 Hamana, T., Oguri, M., Shirasaki, M., & Sato, M. 2012, MNRAS, 425, 2287 [NASA ADS] [CrossRef] [Google Scholar]
 HarnoisDéraps, J., Amon, A., Choi, A., et al. 2018, MNRAS, 481, 1337 [NASA ADS] [CrossRef] [Google Scholar]
 HarnoisDéraps, J., Pen, U.L., Iliev, I. T., et al. 2013, MNRAS, 436, 540 [NASA ADS] [CrossRef] [Google Scholar]
 HarnoisDéraps, J., van Waerbeke, L., Viola, M., & Heymans, C. 2015, MNRAS, 450, 1212 [NASA ADS] [CrossRef] [Google Scholar]
 HarnoisDéraps, J., Giblin, B., & Joachimi, B. 2019, A&A, 631, A160 [CrossRef] [EDP Sciences] [Google Scholar]
 Hartlap, J., Simon, P., & Schneider, P. 2007, A&A, 464, 399 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Heitmann, K., Lawrence, E., Kwan, J., Habib, S., & Higdon, D. 2014, ApJ, 780, 111 [Google Scholar]
 Hetterscheidt, M., Erben, T., Schneider, P., et al. 2005, A&A, 442, 43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Heydenreich, S., Brück, B., & HarnoisDéraps, J. 2021, A&A, in press, https://doi.org/10.1051/00046361/202039048 [Google Scholar]
 Heymans, C., Van Waerbeke, L., Miller, L., et al. 2012, MNRAS, 427, 146 [NASA ADS] [CrossRef] [Google Scholar]
 Heymans, C., Tröster, T., Asgari, M., et al. 2021, A&A, in press, https://doi.org/10.1051/00046361/202039063 [Google Scholar]
 Hikage, C., Oguri, M., Hamana, T., et al. 2019, PASJ, 71, 43 [NASA ADS] [CrossRef] [Google Scholar]
 Hilbert, S., Marian, L., Smith, R. E., & Desjacques, V. 2012, MNRAS, 426, 2870 [NASA ADS] [CrossRef] [Google Scholar]
 Hildebrandt, H., Viola, M., Heymans, C., et al. 2017, MNRAS, 465, 1454 [Google Scholar]
 Ivezić, Ž., Kahn, S. M., Tyson, J. A., et al. 2019, ApJ, 873, 111 [NASA ADS] [CrossRef] [Google Scholar]
 Joudaki, S., Hildebrandt, H., Traykova, D., et al. 2020, A&A, 638, L1 [CrossRef] [EDP Sciences] [Google Scholar]
 Kacprzak, T., Kirk, D., Friedrich, O., et al. 2016, MNRAS, 463, 3653 [NASA ADS] [CrossRef] [Google Scholar]
 Kaiser, N., & Squires, G. 1993, ApJ, 404, 441 [NASA ADS] [CrossRef] [Google Scholar]
 Kilbinger, M. 2015, Rep. Prog. Phys., 78, 086901 [Google Scholar]
 Kilbinger, M., Benabed, K., Guy, J., et al. 2009, A&A, 497, 677 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kilbinger, M., Fu, L., Heymans, C., et al. 2013, MNRAS, 430, 2200 [NASA ADS] [CrossRef] [Google Scholar]
 Kilbinger, M., Bonnett, C., & Coupon, J. 2014, Astrophysics Source Code Library [record ascl:1402.026] [Google Scholar]
 Kratochvil, J. M., Haiman, Z., & May, M. 2010, Phys. Rev. D, 81, 043519 [NASA ADS] [CrossRef] [Google Scholar]
 Kratochvil, J. M., Lim, E. A., Wang, S., et al. 2012, Phys. Rev. D, 85, 103513 [CrossRef] [Google Scholar]
 Kruse, G., & Schneider, P. 1999, MNRAS, 302, 821 [NASA ADS] [CrossRef] [Google Scholar]
 Kruse, G., & Schneider, P. 2000, MNRAS, 318, 321 [NASA ADS] [CrossRef] [Google Scholar]
 Laigle, C., McCracken, H. J., Ilbert, O., et al. 2016, ApJS, 224, 24 [NASA ADS] [CrossRef] [Google Scholar]
 Laureijs, R., Amiaux, J., Arduini, S., et al. 2011, ArXiv eprints [arXiv:1110.3193] [Google Scholar]
 Li, Z., Liu, J., Matilla, J. M. Z., & Coulton, W. R. 2019, Phys. Rev. D, 99, 063527 [NASA ADS] [CrossRef] [Google Scholar]
 Lin, C.A., & Kilbinger, M. 2015, A&A, 576, A24 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Liu, J., & Madhavacheril, M. S. 2019, Phys. Rev. D, 99, 083508 [NASA ADS] [CrossRef] [Google Scholar]
 Liu, J., Petri, A., Haiman, Z., et al. 2015a, Phys. Rev. D, 91, 063507 [NASA ADS] [CrossRef] [Google Scholar]
 Liu, X., Pan, C., Li, R., et al. 2015b, MNRAS, 450, 2888 [NASA ADS] [CrossRef] [Google Scholar]
 Marian, L., Smith, R. E., Hilbert, S., & Schneider, P. 2012, MNRAS, 423, 1711 [NASA ADS] [CrossRef] [Google Scholar]
 Marian, L., Smith, R. E., Hilbert, S., & Schneider, P. 2013, MNRAS, 432, 1338 [NASA ADS] [CrossRef] [Google Scholar]
 Martinet, N., Bartlett, J. G., Kiessling, A., & Sartoris, B. 2015, A&A, 581, A101 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Martinet, N., Schneider, P., Hildebrandt, H., et al. 2018, MNRAS, 474, 712 [NASA ADS] [CrossRef] [Google Scholar]
 Maturi, M., Fedeli, C., & Moscardini, L. 2011, MNRAS, 416, 2527 [NASA ADS] [CrossRef] [Google Scholar]
 McClintock, T., Varga, T. N., Gruen, D., et al. 2019, MNRAS, 482, 1352 [NASA ADS] [CrossRef] [Google Scholar]
 Merten, J., Giocoli, C., Baldi, M., et al. 2019, MNRAS, 487, 104 [CrossRef] [Google Scholar]
 Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493 [NASA ADS] [CrossRef] [Google Scholar]
 Patton, K., Blazek, J., Honscheid, K., et al. 2017, MNRAS, 472, 439 [Google Scholar]
 Parroni, C., Cardone, V. F., Maoli, R., & Scaramella, R. 2020, A&A, 633, A71 [CrossRef] [EDP Sciences] [Google Scholar]
 Peel, A., Pettorino, V., Giocoli, C., Starck, J.L., & Baldi, M. 2018, A&A, 619, A38 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Peel, A., Lalande, F., Starck, J.L., et al. 2019, Phys. Rev. D, 100, 023508 [CrossRef] [Google Scholar]
 Petri, A., Liu, J., Haiman, Z., et al. 2015, Phys. Rev. D, 91, 103511 [NASA ADS] [CrossRef] [Google Scholar]
 Petri, A., May, M., & Haiman, Z. 2016, Phys. Rev. D, 94, 063534 [CrossRef] [Google Scholar]
 Pires, S., Vandenbussche, V., Kansal, V., et al. 2020, A&A, 638, A141 [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration VI. 2020, A&A, 641, A6 [Google Scholar]
 Ribli, D., Pataki, B. Á., & Csabai, I. 2019, Nat. Astron., 3, 93 [Google Scholar]
 Sartoris, B., Biviano, A., Fedeli, C., et al. 2016, MNRAS, 459, 1764 [Google Scholar]
 Schirmer, M., Erben, T., Hetterscheidt, M., & Schneider, P. 2007, A&A, 462, 875 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schneider, P. 1996, MNRAS, 283, 837 [Google Scholar]
 Schneider, P., & Bartelmann, M. 1997, MNRAS, 286, 696 [Google Scholar]
 Schneider, P., & Lombardi, M. 2003, A&A, 397, 809 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schneider, P., van Waerbeke, L., Kilbinger, M., & Mellier, Y. 2002, A&A, 396, 1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schrabback, T., Applegate, D., Dietrich, J. P., et al. 2018, MNRAS, 474, 2635 [Google Scholar]
 Seitz, C., & Schneider, P. 1997, A&A, 318, 687 [NASA ADS] [Google Scholar]
 Sellentin, E., & Heavens, A. F. 2016, MNRAS, 456, L132 [Google Scholar]
 Sellentin, E., & Heavens, A. F. 2017, MNRAS, 464, 4658 [Google Scholar]
 Semboloni, E., Schrabback, T., van Waerbeke, L., et al. 2011, MNRAS, 410, 143 [Google Scholar]
 Shan, H. Y., Kneib, J.P., Comparat, J., et al. 2014, MNRAS, 442, 2534 [Google Scholar]
 Shan, H., Liu, X., Hildebrandt, H., et al. 2018, MNRAS, 474, 1116 [Google Scholar]
 Shirasaki, M., Yoshida, N., Ikeda, S., Oogi, T., & Nishimichi, T. 2019, ArXiv eprints [arXiv:1911.12890] [Google Scholar]
 Smail, I., Ellis, R. S., & Fitchett, M. J. 1994, MNRAS, 270, 245 [Google Scholar]
 Smith, R. E., Peacock, J. A., Jenkins, A., et al. 2003, MNRAS, 341, 1311 [Google Scholar]
 Spergel, D., Gehrels, N., Baltay, C., et al. 2015, ArXiv eprints [arXiv:1503.03757] [Google Scholar]
 Takada, M., & Jain, B. 2003, MNRAS, 344, 857 [Google Scholar]
 Takada, M., & Jain, B. 2004, MNRAS, 348, 897 [Google Scholar]
 Takahashi, R., Sato, M., Nishimichi, T., Taruya, A., & Oguri, M. 2012, ApJ, 761, 152 [Google Scholar]
 Taylor, A., & Joachimi, B. 2014, MNRAS, 442, 2728 [Google Scholar]
 Troxel, M. A., MacCrann, N., Zuntz, J., et al. 2018, Phys. Rev. D, 98, 043528 [Google Scholar]
 Van Waerbeke, L., Benjamin, J., Erben, T., et al. 2013, MNRAS, 433, 3373 [Google Scholar]
 Vicinanza, M., Cardone, V. F., Maoli, R., Scaramella, R., & Er, X. 2018, Phys. Rev. D, 97, 023519 [Google Scholar]
 Vicinanza, M., Cardone, V. F., Maoli, R., et al. 2019, Phys. Rev. D, 99, 043534 [Google Scholar]
 Weiss, A. J., Schneider, A., Sgier, R., et al. 2019, J. Cosmology Astropart. Phys., 2019, 011 [Google Scholar]
 Yang, X., Kratochvil, J. M., Wang, S., et al. 2011, Phys. Rev. D, 84, 043529 [Google Scholar]
 Yuan, S., Pan, C., Liu, X., Wang, Q., & Fan, Z. 2019, ApJ, 884, 164 [Google Scholar]
 Zorrilla Matilla, J. M., Haiman, Z., Hsu, D., Gupta, A., & Petri, A. 2016, Phys. Rev. D, 94, 083506 [Google Scholar]
 Zürcher, D., Fluri, J., Sgier, R., Kacprzak, T., & Refregier, A. 2021, J. Cosmol. Astropart. Phys., 01, 028 [Google Scholar]
Appendix A: Effect of masks
In this section we explore the impact of observational masks on the forecasts. This is important to verify that masks do not introduce biases when applying M_{ap} statistics to observations.
Similarly to Laureijs et al. (2011), we build a mask that captures prime features of the expected masks in Euclid, as a collection of circular star masks with random positions. We therefore neglect any contribution from noncircular artifacts, for example bad pixels, cosmic rays, and charge transfer inefficiency. Those generally occupy a smaller fraction of the image such that only accounting for stars is a good first order approximation. We draw random radii for the stellar masks from a uniform distribution between 0.29′ and 8.79′, corresponding to 0.5 and 15 pixels. The largest radius is chosen to match that of the bright Gaia star masks used in the HSC survey (Coupon et al. 2018). This value can be seen as an upper limit as the Euclid point spread function will be at least four times smaller than that of the Subaru groundbased observation. Using a uniform distribution for star radii does also not reflect the star luminosity function of the Milky Way, overestimating the number of bright stars. With this setup, we build the mask out of 1000 random stars, covering about 20% of the total 100 deg^{2} area of our simulation mocks, and resembling the Euclid mask used in Pires et al. (2020). We use the same mask for all mocks. This is correct for the model determination but slightly underestimates the impact of the mask in the covariance matrix. We recall that our goal is not to apply an exact masking of the data but rather to include the main mask features in order to measure the impact of masking on the M_{ap} map computation. Our mask is presented in Fig. A.1 for the same 1 deg^{2} as in Fig. 4 and for the real and Fourier space reconstructions. The two reconstructions are visually almost impossible to differentiate one from another, and no remaining structures can be seen in their difference.
Fig. A.1.
Aperture mass maps of 1 deg^{2} reconstructed in real (left) and Fourier space (middle) after applying the star mask, and their difference (right). 
We find that our mask has a negligible impact on the forecasts, both in real and Fourier space, and both without and with tomography. The only difference with the unmasked case is a decrease in the precision consistent with the loss of 20% of the probed area.
These results are good news for M_{ap} statistics as they tend to show that we can safely use the Fourier approach to speed up computation without introducing any bias to the cosmological constraints. We however recall that we used several approximations in the mask building and that we also cannot forecast biases for the full coverage of future Stage IV surveys, such that a more complex analysis will be required to verify the possible impact of masks on these surveys to greater precision.
Appendix B: Effect of the cosmological parameters sampling
In this section we test the robustness of our forecasts to the number of points explored in the simulated parameter space. This is done by recomputing the forecast when removing the n closest points to the observational mocks in the cosmological parameter space. We stress that this represents a worstcase scenario. We perform this test on 1D M_{ap} without including tomography. The former case is representative of most bias behaviors in our study as we found the bias to be similar for the different M_{ap} estimators.
We find that removing up to the five closest points to the truth in the studied parameter space weakly affects the forecast precision, while the bias increases for all parameters compared to our fiducial analysis: by 60% for S_{8}, 15% for w_{0}, and 63% for Ω_{m}. Removing more points significantly increases the bias, that becomes larger than the precision and would dominate the forecasts.
These results suggest that the 25points modeled by the cosmoSLICS are sufficient and well adapted to the purpose of this article. However, we cannot find a simple trend of the bias with respect to the number of points included in the parameter space. This prevents us from predicting the necessary number of points to be explored to reach a 1% accuracy on w_{0} from M_{ap} statistics, which will be mandatory to exploit the full potential of future cosmic shear surveys. A more quantitative study with a highly sampled parameter space is required here. In the case of γ2PCF, subpercent level accuracy could be obtained with 250 cosmoSLICS nodes, as shown in Appendix A of HarnoisDéraps et al. (2019).
All Tables
Cosmological predictions for 100 deg^{2}Euclidlike mock without tomography and with a fiveslice tomography setup, including crossM_{ap} terms, for voids, peaks, 1D M_{ap}, γ2PCF, and for the combination of M_{ap} map estimators with the γ2PCF.
All Figures
Fig. 1.
Redshift distribution normalized to 30 gal.arcmin^{−2}. The red histogram corresponds to the COSMOS2015 photometric redshift distribution after a cut at i′≤24.5, the green curve to the fit of a Smail et al. (1994) function to the COSMOS2015 histogram, and the blue curve to a Fu et al. (2008) fit. In this paper we use the Fu et al. (2008) fit that captures the highredshift tail in contrast to the other function. 

In the text 
Fig. 2.
S/N distribution (or data vectors) of three different M_{ap} samplings: voids (blue), peaks (green), and full 1D distribution (red). Dotted lines correspond to the S/N distributions measured in noiseonly maps. Numbers are quoted for a 1024 × 1024 pixels map of 100 deg^{2} with a filter size θ_{ap} = 10′. 

In the text 
Fig. 3.
Variation of the excess number of M_{ap} pixels over noise with cosmology. The black curve corresponds to the excess 1D M_{ap} in the observation mock, with error bars from the diagonal of the covariance matrix. Colorcoded are variations in Ω_{m} (top left), S_{8} (top right), h (bottom left), and w_{0} (bottom right). A smooth color gradient indicates a correlation between the data vector and the probed cosmological parameter. No tomography has been performed here. 

In the text 
Fig. 4.
Aperture mass maps of 1 deg^{2} obtained from the full redshift distribution (0 < z_{all} < 3), for five 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 crosscombinations 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 crossM_{ap} terms allow us to retain the information about correlations of structures in between slices. 

In the text 
Fig. 5.
Correlation matrix computed from 1D M_{ap} in 928 simulations with independent noise realizations (sample variance, shape noise, and position noise). 

In the text 
Fig. 6.
Correlation matrix computed from the combined γ2PCF and 1D M_{ap} DVs with five tomographic slices and including crossterms. ξ_{±} 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. 

In the text 
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). 

In the text 
Fig. 8.
Forecast of cosmological parameters from 1D M_{ap} in a 100 deg^{2} survey at Euclid depth. Marginalized 2D (1 and 2σ contours) and 1D (full likelihood) constraints are displayed in blue without tomography and in green for a tomographic analysis with 5 redshift slices and including auto and crossM_{ap} terms between redshift slices. Dashed lines correspond to the 1σ constraints in the 1D marginalized likelihood. The red crosses and lines indicate the input values while the blue and green crosses indicate the best estimate in the 2D constraints. Gray crosses correspond to parameters of the 25 cosmoSLICS simulations that are used to estimate the cosmology dependence of 1D M_{ap}. 

In the text 
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 setup (autoM_{ap} only) and to the improved tomography (auto and crossM_{ap}). Black dots represent the γ2PCF (both auto and crosscorrelations) 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 crossterms. 

In the text 
Fig. 10.
Forecast of cosmological parameters from 1D M_{ap}, γ2PCF and their combination in a 100 deg^{2}Euclidlike survey without tomography. Marginalized 2D (1 and 2σ contours) and 1D (full likelihood) constraints are displayed in black for the γ2PCF, green for 1D M_{ap}, and purple for their combination. The red crosses and lines indicate the truth. Gray crosses correspond to the 25 simulation points from which the model is built. 

In the text 
Fig. 11.
Same as Fig. 10 for a tomographic analysis with five redshift slices, including auto and crossM_{ap} terms. 

In the text 
Fig. A.1.
Aperture mass maps of 1 deg^{2} reconstructed in real (left) and Fourier space (middle) after applying the star mask, and their difference (right). 

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