Optimizing simulation parameters for weak lensing analyses involving non-Gaussian observables
José Manuel Zorrilla Matilla Stefan Waterval Zoltán Haiman




José Manuel Zorrilla Matilla
jzorrilla@astro.columbia.edu

0000-0002-6267-716X]José Manuel Zorrilla Matilla
Department of Astronomy, Columbia University, New York, USA


Department of Physics, ETH, Zürich, Switzerland
Department of Physics, New York University Abu Dhabi, Abu Dhabi, United Arab Emirates


Department of Astronomy, Columbia University, New York, USA



We performed a series of numerical experiments to quantify the sensitivity of the predictions for weak lensing statistics obtained in raytracing DM-only simulations, to two hyper-parameters that influence the accuracy as well as the computational cost of the predictions: the thickness of the lens planes used to build past light-cones and the mass resolution of the underlying DM simulation. The statistics considered are the power spectrum and a series of non-Gaussian observables, including the one-point probability density function, lensing peaks, and Minkowski functionals. 
Counter-intuitively, we find that using thin lens planes (< 60 h^-1 Mpc on a 240 h^-1 Mpc simulation box) suppresses the power spectrum over a broad range of scales beyond what would be acceptable for an LSST-type survey. A mass resolution of 7.2× 10^11 h^-1 M_⊙ per DM particle (or 256^3 particles in a (240 h^-1 Mpc)^3 box) is sufficient to extract information using the power spectrum and non-Gaussian statistics from weak lensing data at angular scales down to 1 arcmin with LSST-like levels of shape noise.











§ INTRODUCTION
 

Weak gravitational lensing (WL) enables the mapping of the distribution of dark matter (DM) in the universe on large scales and as a result is a powerful probe to infer cosmological parameters such as σ_8 and Ω_m (see comprehensive reviews by, e.g. ,   and ). In practice, the lensing signal can be extracted from statistical measurements of the shapes of background galaxies, distorted by deflections in the path of light rays as they traverse the vicinity of matter over- and under-densities.

Upcoming surveys such as DESI <cit.>, LSST <cit.>, Euclid <cit.>, WFIRST-AFTA <cit.> and SKA <cit.>, will provide WL data of unprecedented quality and quantity and will require correspondingly accurate and precise models to extract information from these datasets. WL observables delve into the non-linear regime, which can be forward-modeled by ray-tracing photons through high-resolution simulated dark matter (DM) density fields <cit.>. Simulating the large volumes encompassed by future surveys is computationally expensive, especially when a high-dimensional parameter space needs to be explored. Different ideas to reduce the cost of forward modeling WL observables have been put forward and tested. These include using approximate codes to simulate the evolution of the matter density field – for example ICE-COLA <cit.>, L-PICOLA <cit.> or FastPM <cit.>–, analytic or semi-analytic models – for example Camelus <cit.> and machine learning approaches – for example generative adversarial networks (GANs; e.g. ). While analytic models can predict two-point statistics with sufficient accuracy <cit.>, higher-order statistics which capture non-Gaussian information from the non-linear regime require a numerical approach <cit.>.

In this paper we study, within the framework of raytracing N-body simulations, the sensitivity of the power spectrum (PS) and several non-Gaussian statistics commonly used in WL studies: the one-point probability density function (PDF), lensing peak counts, and the full set of Minkowski functionals (MFs) to the two hyper-parameters with the highest impact on the computational cost of the simulations: (i) the thickness of the lens planes used to construct the past light cones in raytracing and (ii) the mass resolution of the N-body simulations used to model the underlying matter density field. Previous studies have already tackled some of these or related aspects. For instance, <cit.> verified the effect on the measured PS of the mass resolution of the underlying DM simulation, the grid size and interpolation method used when raytracing and also studied the contribution of super-sample modes to the PS variance. <cit.> looked at the effect of the resolution of the 2D lens planes on the measured convergence power spectrum, and evaluated the non-Gaussian contribution to its covariance matrix. <cit.> investigated the effect of mass resolution and comoving distance between lens planes on the convergence power spectrum, skewness and kurtosis. This work revisits the sensitivity of the measured convergence power spectrum to the mass resolution and distance between lens planes, and extends the analysis of numerical convergence analyses to non-Gaussian statistics within the framework of a tomographic WL analysis of a LSST-like survey.

The manuscript is organized as follows. In  <ref> we describe our simulation pipeline ( <ref>), how we quantify the impact of the hyper-parameters ( <ref>), and the statistics used to assess that impact ( <ref>). In  <ref> we report and discuss the results of the analysis. Finally, we summarize our conclusions in  <ref>.







§ METHODS




 §.§ Simulating convergence maps
 

Our analysis is based on lensing statistics measured on convergence maps obtained from ray-tracing dark matter only N-body simulations. Since weak lensing probes large-scale structures in the non-linear regime, direct simulations offer a way to characterize the WL signal.

The N-body simulations are run using the publicly available Tree-PM code Gadget2 <cit.>,  which evolves Gaussian initial conditions generated with NGenIC <cit.>. The initial conditions are defined by power spectra computed with CAMB <cit.>. The positions of the particles at different redshifts are used to build past light-cones. The trajectory of a bundle of rays is followed along past light-cones according to the multi-plane algorithm to generate synthetic convergence maps. We refer the reader to <cit.> for a detailed step-by-step description of our pipeline and its implementation in LensTools,  and to <cit.> for a review of the theoretical basis of the algorithms used. 

All the simulations have the same comoving volume of (240h^-1 Mpc)^3, phases for the initial conditions, and underlying cosmology. The cosmological parameters are consistent with <cit.>: (H_0, Ω_m, Ω_Λ, Ω_b, w_0, σ_8, n_s) = (67.7 km s^-1Mpc^-1, 0.309, 0.691, 0.0486, -1.0, 0.816, 0.967). As an illustration, we compare the convergence power spectrum measured in our highest-resolution simulation (1024^3 DM-only particles, or a particle mass of 1.1×10^10  M_⊙), with flat-sky approximation predictions (see Eq. (<ref>)) computed using the halofit model <cit.>, as implemented in the
Core Cosmology Library (CCL; ). In data smoothed at 1  arcmin resolution with the presence of the shape noise levels considered here, the agreement between theory and simulations stays within 20% for ℓ < 12,000. In noiseless data, a loss of power in the simulations is noticeable, especially for sources at low redshift. The worst case corresponds to the bin with sources at z=0.5, which shows a 50% loss in power at ℓ≈7,500. In all cases we measure excess power at large angular scales relative to the theoretical expectation (up to ≈20% for the noiseless case and ≈13% in the presence of shape noise). The reason of this excess is the sample variance in the matter power spectrum, as low multipoles probe scales that are comparable with the size of our simulated volume. Since the purpose of this study is the calibration of hyper-parameters, we are interested in relative changes and the sample variance at large scales should not impact our conclusions. Furthermore, we are mostly interested in relative changes on small scales, where the use of non-Gaussian observables is relevant. To perform inference from data, this sample variance would need to be tackled, for instance, by simulating a larger volume (or many independent N-body realizations of our simulation volume).



For each configuration, a single volume is simulated and reused through random shifts and rotations to generate as many as 𝒪(10^4) pseudo-independent past light-cones <cit.>. Since each convergence map covers only 3.5×3.5 deg^2 on the sky at this redshift, the flat-sky approximation holds.  To account for their intrinsic ellipticity, Gaussian random shape noise is added independently at each of the 1024×1024 pixels in each convergence map. The standard deviation of this noise,


    σ_pix = √(σ_ε^2/2n_gA_pix),


depends on the variance of the galaxies' intrinsic ellipticity (assumed to be σ_ϵ=0.35), the solid area covered by a pixel (A_pix=0.04 arcmin^2), and the effective galaxy number density. 

We consider an LSST-like survey with a footprint of 20,000  deg^2 and  galaxy density distribution n(z)∝ z^2 exp[-2z] <cit.>. Because a small number of redshift bins suffices to extract most of the tomographic information encoded in weak lensing data sets <cit.>, we use five tomographic bins. In each bin, we assume all lensed galaxies are at a fixed redshift of [0.5, 1.0, 1.5, 2.0, 2.5]. The galaxy densities we consider in each redshift bin are [8.83, 13.25, 11.15, 7.36, 4.26]  arcmin^-2.

A final smoothing is applied to both noiseless and noisy maps with a Gaussian kernel of 1 arcmin standard deviation. While only results from maps with shape noise are relevant for the analysis of future survey data, we show also results for smoothed, noiseless convergence maps, since the effect of different simulation choices are often more discernible in those.





 §.§ Assessing the impact of hyper-parameters
 

The accuracy of the forward model as a function of different values of the hyper-parameters is assessed by comparing the statistics of observables measured over 10,048 convergence maps simulated for each configuration. As explained in  <ref>, for each set of hyper-parameters, all 10,048 maps are generated from a single, recycled, N-body simulation.

We chose a "fiducial" configuration as a reference. The difference between an observable's mean for all cases and the fiducial model is compared with a standard error. As standard error, we consider 3 standard deviations measured on the fiducial model's maps, scaled to a survey sky area of 2×10^4 deg^2 (commensurate with LSST's). This scaled standard deviation represents a lower bound on the uncertainty expected in future surveys, as it includes only the statistical error from intrinsic ellipticities. For a review of all sources of error in WL surveys, see for instance <cit.>.

Ultimately, the most relevant metric is how much the lens plane thickness affects the inference of cosmological parameters. A definite answer to this question requires the calculation of the credible contours for the parameters of interest. Within the scope of our single-cosmology numerical experiments, we instead look at the reduced χ^2, which combines differences in the mean observables with their covariance matrix.

In general, a given observable is a vector 𝐬 with components corresponding to bins of spherical harmonic index ℓ for the power spectrum, or κ thresholds for the other statistics, such as the one-point PDF, lensing peaks or Minkowski functionals. We compute the χ^2 statistic for each configuration, considering the mean of the fiducial model as ground truth:


    χ^2 = (𝐬-𝐬_fid)^T𝖢^-1(𝐬-𝐬_fid),


where 𝖢^-1 is an unbiased estimator of the precision matrix. We compute (and report) results based both on the precision matrix estimated at the fiducial model, and the precision matrix estimated at each specific configuration. We used the prescription from <cit.> to de-bias the estimator of the precision matrix. For each observable, we report the χ^2 per degree of freedom. Furthermore, to account for the effect of a change in configuration on the covariance matrix, we report two values of the reduced χ^2, one based on the covariance matrix for the reference configuration and another based on the configuration-specific covariance. We consider two configurations as statistically equivalent when their χ^2 per degree of freedom, in the presence of noise, is less than or equal to unity.



 §.§ Hyper-parameter configurations


We analyzed different configurations for two hyper-parameters: the thickness of the lensing planes used in the ray-tracing algorithm, and the mass resolution (number of particles) of the N-body simulations. For each configuration, we ray-traced 10,048 convergence maps, each of these maps with an area of 12.25  deg^2.

The positions of the DM particles were saved at redshifts that allowed the construction of light-cones with lens planes at comoving distances of 20 h^-1, 40 h^-1, 60 h^-1 80 h^-1, and 120 h^-1 Mpc. Thinner lens planes can potentially capture more accurately the evolution of the matter density field with redshift. However the number of planes needed to cover a redshift range increases as the plane thickness decreases, and so do the computational and storage requirements for the simulations. In particular, the two tasks that account for the largest increase in computation time are the calculation of the gravitational potential at the planes (solving a 2D Poisson equation in Fourier space) and the computation of the Jacobian matrix that determines the light ray's deflections at each point on the planes <cit.>. The fiducial case corresponds to a lens plane thickness of 80 h^-1 Mpc, a value that has been typically used in prior work <cit.>, provides 9 independent lens planes per simulation snapshot (increasing the number of pseudo-independent κ maps that can be generated from a single N-body simulation), and is not large enough to show discreteness effects with lensed galaxies at z=2 <cit.>.

The minimum angular scale at which cosmological information can be extracted is limited by the depth of the survey, which determines the number density of galaxies whose shape can be measured, and the accuracy of the forward models used to predict the signal. Baryonic physics <cit.> and intrinsic alignments <cit.> restrict the accuracy of current models at small scales. Matching the mass resolution of the underlying N-body simulations to the scales at which the analysis of the data is reliable saves computational resources.

In our fiducial N-body simulations we used 1024^3 particles, and we ran additional simulations with 128^3, 256^3 and 512^3 particles, which yield mass resolutions per DM particle of 5.7× 10^12 h^-1 M_⊙, 7.2× 10^11 h^-1 M_⊙, 9.0× 10^10 h^-1 M_⊙ and 1.1× 10^10 h^-1 M_⊙. We note that due to memory limitations, we increased the number of CPUs allocated to our higher-resolution simulation, resulting  in a worse scaling than the one that can be achieved in an optimized setup <cit.>.

Tables <ref> and <ref> summarize the computational cost of the main steps involved in our simulation pipeline, and the disk space required for storing the different data products (in practice, not all need to be saved). Performance benchmarks are based on runs using Intel Knights Landing nodes from TACC's Stampede2 supercomputer at the NSF XSEDE facility.










 §.§ Observables


We analyze the effect of changes in the hyper-parameters used to generate the synthetic data on the convergence power spectrum and a series of non-Gaussian observables: the one-point probability density function (PDF), lensing peak counts, and the full set of Minkowski functionals (MFs).

The power spectrum is the Fourier transform of the two-point correlation function, encodes all the information available in a Gaussian random field, and can be accurately predicted from theory on large angular scales. Therefore, it is commonly used in WL analyses <cit.>. We measured it in 20 logarithmic bins between ℓ_min=200 and ℓ_max=12,000, extending slightly above the smoothing scale of 1  arcmin which corresponds to ℓ≈ 10,800. The relatively featureless power spectra can be properly characterized with 20 bins, without yielding a very large data vector when combining the 5 auto-spectra with the 15 cross-spectra for all 5 redshift bins. Since the synthetic maps cover a small field of view, we can use the flat-sky approximation to model the expected power spectrum from measuring the ellipticities of lensed sources at a fixed radial comoving distance, χ_s:


    P_κ(ℓ) = (3H_0^2/2c^2Ω_m)^2 ∫_0^χ_sdχ/a^2(χ)(1-χ/χ_s)^2 P_δ(k=ℓ/χ; χ).


Here P_δ is the (non-linear) matter power spectrum, which depends on the wavenumber k and evolves with redshift, and H_0, Ω_m and c are the Hubble constant, matter density parameter and speed of light, respectively. Shape noise contributes white noise with total power of P_ϵ=σ_ϵ^2/n_g. Smoothing with a Gaussian kernel of width σ_s multiplies the convergence power spectrum by the beam function b_ℓ=exp[-ℓ(ℓ+1)σ_s^2].

The effect that the finite thickness of the lens planes has on the measured convergence power spectrum can be modelled as convolving the matter density power spectrum with a window function. Following <cit.>, it suffices to substitute  P_δ in Eq. (<ref>) by


    P'_δ(k; χ) = Δχ/π∫_0^k^ max_∥ dk_∥P_δ(k) sinc^2(k_∥Δχ/2),


where k=√(k_∥^2+k_⊥^2), k_⊥=ℓ/χ, sinc(x)=sin x/x, and Δχ is the thickness of the lens planes used to build past light-cones in our simulations.

Changes in the mass resolution of the N-body simulations will also affect the measured κ power spectrum. On one hand, the matter power will be suppressed on scales below the simulation resolution. On the other hand, an additional shot-noise component will result from the discreteness of the simulated matter density field. Both effects can be accounted for as in <cit.>, by smoothing the matter power spectrum with a Gaussian kernel of with σ_N and adding a Poisson noise contribution corresponding to the finite number density n_p of simulated particles,


    P”_δ(k;χ) = P_δ(k;χ)exp(-σ_N^2k^2)+1/n_p.


Most scales in our simulated maps probe non-linear structures, and as a result, there is non-Gaussian information not captured by the convergence power spectrum. To assess the impact of the different simulation parameters on the non-Gaussian information content of convergence maps, we analyze the effect on the one-point probability density function, the distribution of lensing peaks, and the Minkowski functionals.

The convergence PDF has been shown to help break the power spectrum degeneracy in the Ω_m-σ_8 plane, improving cosmological constraints from the power spectrum alone (e.g. ). It has also been used in order to assess the ability of weak lensing measurements to infer the neutrino masses <cit.>. We use the binned PDF measured on our (smoothed) simulated maps within the range [-3σ_κ, 5σ_κ], where σ_κ is the standard deviation of the maps corresponding to the highest mass resolution (1024^3 particles), and the fiducial, 80 h^-1  Mpc, lens planes. For each redshift bin in the noiseless case, we use σ_κ = [0.007, 0.014, 0.019, 0.023, 0.027 ], and for the maps that include shape noise σ_κ = [0.034, 0.031, 0.035, 0.043, 0.055]. We compare the measurements over simulations with the analytic expectation for Gaussian random fields (GRFs) with variance


    σ_0^2=∫_0^∞ℓ dℓ/2πP_κ(ℓ).


Lensing peaks are local maxima of the convergence field; their distribution as a function of their height, κ, was proposed to constrain cosmological parameters <cit.>, and has been used extensively and successfully ever since <cit.>. We measured peak histograms in 50 linear bins covering the range κ∈[-0.5, 5.0] in units of the r.m.s. shape noise in each redshift bin. For the maps smoothed on 1  arcmin scales, the r.m.s. of the shape noise for the redshift bins [0.5, 1.0, 1.5, 2.0, 2.5] is σ_ n=[0.033, 0.027, 0.030, 0.036, 0.048 ]. As we did for the PDF, we compare measurements with the predictions for GRFs with the same convergence power spectrum as the simulated κ maps. The analytic expressions for the number of peaks in GRFs was derived in <cit.>. We reproduce here the relevant equations for convenience. The number of lensing peaks in a GRF for a given height ν=κ/σ_0 is given by:


    n_ peaks(ν) = 1/(2π)^3/2θ^2exp(-ν^2/2)G(γ, γν),


where θ=√(2)σ_1/σ_2 and γ=σ_1^2/(σ_0σ_2) depend on the moments of the convergence power spectrum


    σ_p^2∫_0^∞ℓ d ℓ/2πℓ^2pP_κ(ℓ)


and


    G(γ, x) = (x^2-γ^2)[1-1/2 erfc(x/√(2(1-γ^2)))]  + 
    
         x(1-γ^2)exp(-x^2/2(1-γ^2))/√(2π(1-γ^2)) + 
    exp(-x^2/3-2γ^2)/√(3-2γ^2)[1-1/2 erfc(x/√(2(1-γ^2)(3-2γ^2)))].


Finally, the three Minkowski functionals that can be defined in a two-dimensional random field, {V_0, V_1, V_2}, measure the area, length, and genus of the subset of points in the random field whose value exceeds a given threshold <cit.>. While V_0 conveys the same information as the PDF, we include it in our measurements for completeness. They have been extensively used to measure or extract non-Gaussian information from cosmological datasets, including weak lensing <cit.>.

We measured the three MFs in 28 linear bins covering the range κ∈[-2, 5] in units of the r.m.s. shape noise for each redshift bin, and compared it with analytical predictions for GRFs. The MFs of two-dimensional GRFs are determined by the moments of their power spectrum, and their expression as a function of ν can be found in <cit.>:


    V_0(ν)    =   1/2 erfc(ν/√(2)) 
    
        V_1(ν)    =   1/8√(2)σ_1/σ_0exp(-ν^2/2)
    
        V_2(ν)    =   1/(2π)^3/2(σ_1/√(2)σ_0)^2νexp(-ν^2/2).


While a perturbative expansion has been found for MFs for non-Gaussian random fields <cit.>, they have been shown to converge too slowly at the relevant ∼arcmin smoothing scales <cit.>, so we do not investigate them here.









§ RESULTS AND DISCUSSION
 

We discuss the effect on lensing statistics of the comoving distance between lens planes (plane thickness) in   <ref> and mass resolution in  <ref>.



 §.§ Lens plane thickness
 


We show the mean percentage difference in the power spectrum for noiseless and noisy convergence maps for different lens plane thicknesses (or comoving distance between lens planes) relative to the fiducial case of 80 h^-1 Mpc in Figure <ref>. For each plane thickness, the power spectrum is measured and averaged over 10,048 convergence maps. For clarity, we do not display the results for the 10 cross-spectra between redshift bins, but the conclusions are not altered. The gray band in each panel represents a standard error corresponding to 3 times the standard deviation measured over the fiducial maps, scaled to a 2×10^4 deg^2 LSST-like survey. This does not incorporate the effect of off-diagonal terms in the covariance matrix, which are fully utilized (including auto and cross spectra) in the χ^2 statistic described in  <ref> to assess the significance of the differences between each configuration and the fiducial case.



The effect of the finite width of lens planes, measured on the simulations, matches well the expectation from the model described in  <cit.>, in which the matter power spectrum is convolved with the window function corresponding to a spherical lens shell. Since the differences between the fiducial and the 60 h^-1 Mpc case are in general sub-percent level, they fall in the linear range of the y-scale and appear to be larger than those for other cases, but they are not. As the lens planes cut structures, there is an appreciable loss of power for thin planes throughout the multipole range considered in this study. The presence of shape noise renders those differences statistically insignificant for the range of thicknesses 60 h^-1 - 120 h^-1 Mpc (see Table <ref>).



The effect on the PDF is shown in Figure <ref>, where the PDF change measured from maps built using different lens plane thicknesses is plotted (thick lines) together with predictions for Gaussian random fields (thin lines). The GRFs have the same power spectrum as the simulated κ maps. The overall impact on the GRFs is straightforward: the loss in power translates into a smaller variance for the random field (see Eq. (<ref>)). As a result we would expect, for planes thinner than the fiducial, a suppression in the tails of the PDF compensated by an enhancement near the peak of the PDF (and the opposite for a the case with thicker lens planes). Qualitatively, the measurements in the presence of shape noise follow that trend, which is to be expected since the shape noise model considered is Gaussian. In the limit of noiseless data, while there is still some qualitative agreement, the quantitative differences with the GRF model are evident.



From a practical perspective, lens planes whose thickness lays in the range of 60-120 h^-1 Mpc yield convergence maps with statistically indistinguishable PDFs, even in the absence of shape noise (see Table <ref>).


The sensitivity of the distribution of lensing peak counts to the width of the lens planes used during raytracing is shown in Figure <ref>. For noiseless data, lens planes thinner than 60 h^-1 Mpc overproduce low-significance peaks while underproducing higher-significance ones. Table <ref> shows how the results from using lens planes within the range 60-120 h^-1 Mpc are statistically indistinguishable. The addition of shape noise reduces the differences further, to the point where using thinner planes, at 40 h^-1 Mpc, could be safe. As seen from the power spectrum, there is not a strong case to move from the fiducial thickness of 80 h^-1 Mpc.



As for the PDF, the changes in the lensing peak counts in κ maps are qualitatively similar to those in GRFs with the same power spectrum, but there are quantitative differences. The effect is also similar to the one seen for the PDF, and can be understood intuitively: thicker lens planes provide more power on larger scales that modulates the heights of local maxima, increasing the variance of their distribution (and conversely for thinner planes).

Past studies have found that the use of the Born approximation, in which the convergence is directly estimated weighting the projected matter density, can induce a significant bias to predictions of non-Gaussian observables such as the skewness and kurtosis <cit.>. Since in the limit of very thick planes, the multi-plane ray-tracing algorithm is similar to the Born approximation (except for the lack of redshift evolution in the matter density field within each plane), any attempt to predict peak histograms beyond the range of thickness explored in this study, would need to be validated.




The effect of lens plane thickness on the MFs is different for each functional (see Figure <ref>). As expected, shape noise "gaussianizes" the convergence field and as a result measured differences follow qualitatively the predictions for GRFs. This is not the case in the limit of noiseless data. The κ field is less sensitive to changes in the plane thickness than what would be expected for GRFs, especially V_1. For V_2, the region of κ≈0, where the functional changes sign, yields noisy measurements.



In terms of χ^2 (Table <ref>), the MFs are more demanding than other non-Gaussian observables. While in the presence of shape noise, lens planes in the range of 60-120 h^-1 Mpc yield measurements that are statistically equivalent, that is not the case anymore in the limit of noiseless data for V_1 and V_2.

As explained in  <ref>, the reduced χ^2 allows us to assess the impact on inference, of the differences in the measured statistics from synthetic maps generated with different hyper-parameters. Changes in the covariance matrix are taken into consideration when computing χ^2. We use either the covariance matrix from maps in the reference configuration, or an alternative covariance from maps recomputed after changing the thickness of the lens planes. In both cases, we use the full covariance to compute χ^2, including, for instance, the 10 cross power spectra not shown in Figure <ref>. As an illustration, we display the full covariance matrix for the power spectrum measured on noisy maps in Figure <ref>.



The effect of lens planes in the covariance is small (for example, less than 20% for the dominant, diagonal elements for the power spectrum). As a result, the conclusions are the same regardless of which covariance is used. Also,  <cit.> have analyzed the impact of  non-Gaussian contributions to the lensing power spectrum covariance on parameter inference. They find that mean changes in the covariance matrix of ≈ 20% translate to ≤ 5% changes in the parameter uncertainties.


Table <ref> collects the results for the observables under study. With the shape noise levels considered, lens planes whose thickness is in the range of 60-120 h^-1 Mpc yield statistically indistinguishable measurements for all the statistics we have studied, and are therefore safe to use. For a scheme, like ours, that recycles N-body snapshots to build different past light-cones, there is a trade-off between a loss in power induced by the lens plane's window function, and the number of lens planes that can be built from a single N-body simulation. The fiducial choice of 80 h^-1 Mpc seems reasonable, since it allows to generate as many as a few× 10^4 pseudo-independent convergence maps from a single simulation <cit.>. In the limit of noiseless data this is not the case anymore for the power spectrum and MFs. For the power spectrum, the effect of the lens plane thickness can be incorporated analytically.

These results do not take into consideration the uncertainty in the measured covariance matrices: if they are overestimated, our results could change. Nevertheless, the required errors must be substantial: the true covariance would need to be overestimated by more than 150% (its true elements be smaller than 0.4× the elements of the matrix measured on our simulations) for the Minkowski functionals (the most sensitive statistic to this hyper-parameter) measured on maps generated with 80 h^-1 Mpc lensing planes to be statistically distinguishable from those measured on maps generated with 120 h^-1 Mpc lensing planes.










 §.§ Mass resolution
 


Figure <ref> displays the mean percentage difference in the power spectrum for noiseless and noisy maps for different mass resolutions, compared to the highest-resolution case. The number of particles in the four configurations we tested are 128^3, 256^3, 512^3 and 1024^3. As done  in Figure <ref> with the sensitivity to the lens plane thickness, a standard error corresponding to 3 standard deviations for the fiducial case is shaded for reference in Figure <ref>.



The main difference between the power spectra is an increase in power on small scales, with a relative loss of power on intermediate scales, as the mass resolution decreases. Qualitatively, this behavior matches what is expected from the matter power spectrum: a loss of power at intermediate scales as the mass resolution decreases (captured by the first term in Eq. (<ref>)) followed by a sharp increase at small scales due to shot noise. Similar effects have been described in past studies <cit.>. We cannot fit all the data with the simple model in Eq. (<ref>), though. For example, a smoothing scale of σ_N=5% of the average inter-particle separation gives results that are qualitatively in line with the relative power change when going from 256^3 to 1024^3 particles, but not for the other resolutions tested. An effect not taken into consideration in that simple model is that of the softening length on the relaxation time of the simulations. All of our simulations kept that length constant, and the softening for the lower resolution runs may have been insufficient. Another effect that can induce a loss in power is the reduction in the linear growth factor due to the discreteness of the initial conditions, as described in  <cit.>. The presence of shape noise, which dominates at small scales, mitigates both effects, rendering the measured statistic more forgiving to the resolution of the underlying N-body simulation.




The effect of mass resolution on the one-point probability density function is not as straightforward as for the lens plane thickness. The GRF model does not work particularly well for noiseless data, as can be seen in Figure <ref>. Except for the lowest-resolution simulation, shape noise still keeps any differences within the standard error of the simulations themselves. In the absence of noise, 512^3 and 1024^3 particles seem to yield equivalent results.








Lensing peaks are a more robust statistic to the mass resolution (see Figure <ref>), and the configurations with 256^3 and 512^3 particles are statistically indistinguishable from the 1024^3 fiducial for both noisy and noiseless data. The model with the lowest mass resolution yields histograms whose differences from the fiducial case clearly exceed the statistical uncertainty. The difference is more important for low-significance peaks, which can be induced by the Poisson shot noise present in the underlying N-body simulation. Shape noise reduces the differences at the low significance tail, where peaks are noise-dominated. The GRF model for peak counts seems to work better than for the PDF.








The effect of mass resolution on the measured MFs is diplayed in Figure <ref>. As with the PDF, a GRF model does not fit as well the impact of changes in mass resolution as it did for changes in lens plane thickness. This may be partly due to the relatively larger change in power, as well as the ℓ-dependence of that change. In general, measurements performed on simulations are less sensitive than what would be expected from GRFs, and their changes remain within errors for most cases.








Table <ref> collects the statistical significance of changes in the mass resolution of the N-body simulations used to build synthetic κ maps on the observables under study. Most of the considerations discussed in  <ref> about differences in the covariance matrices for the observables apply here as well. For example, the differences for the power spectrum covariance are modest with differences in the dominant, diagonal terms, of less than 3% between the covariances for the 256^3 and 1024^3 particle simulations. For analyses in the presence of shape noise levels commensurate with the ones considered in this study, 256^3 particles suffice. If the multipole range is limited to ≈ 3500, even lower resolutions can be used for the power spectrum (or alternatively the same 256^3 resolution can be used for negligible levels of shape noise). The non-Gaussian statistics that are the most sensitive to resolution effects are the MFs.

As with the results for the thickness of the lensing planes, the figures in Table <ref> do not incorporate the possible effect of errors in the covariance matrices used in the calculation of χ^2. In this case, the Minkowski functionals are also the most sensitive statistic. For their measurements on the 256^3 particles' simulation to be statistically different than the measurements on the 1024^3 particles' maps, the true covariance would need to be overestimated by more than 117%.






§ CONCLUSIONS
 
We performed a series of numerical experiments to test the influence of the lens plane thickness and the mass resolution of ray-traced N-body simulations on commonly used WL statistics: the convergence power spectrum, the one-point probability density function, lensing peak counts, and Minkowski functionals. While our simulations cannot be used directly to analyze survey data, they can serve to guide design choices in studies of non-Gaussian statistics on small scales. They set some minimal requirements to avoid significant biases in predictions obtained from the application of the multiple lens plane algorithm. 

We have found that in the multi-plane ray-tracing algorithm, lens planes in the range 60-120 h^-1 Mpc are safe to use in analysis of WL data with galaxy densities commensurate with those in LSST-like surveys. Thinner planes will result in a loss of power across a wide range of multipoles, that can in principle be accounted for analytically for the power spectrum but not for non-Gaussian statistics. On the thick-plane end, there are no biases induced in the observables, and the computational time for raytracing past light cones can be reduced (e.g. by ≈ 33% when 120 h^-1 Mpc planes are used instead of 80 h^-1 Mpc). However, such a choice implies a larger number of N-body simulations to generate the same number of pseudo-independent realizations of the κ field, increasing the overall computation budget.

In order to analyze WL data sets at angular resolutions of 1 arcmin with LSST levels of shape noise, simulations with mass resolutions of 7.2× 10^11 M_⊙ per DM particle are sufficient (corresponding to 256^3 particles in a 240 h^-1 Mpc simulation box), even if non-Gaussian statistics are included in the analysis. Moving, for instance from 512^3 to 256^3 particles could bring computational time savings of ≈ 77%.




§ ACKNOWLEDGMENTS

We thank the anonymous referee for providing insightful feedback that has resulted in an improved manuscript. We also thank Alexandre Refregier for useful comments on the original manuscript and Takashi Hamana for pointing us to an analytic model for the loss of power due to finite-thickness lens planes. The simulations and data analysis for this work were performed on NSF's XSEDE facility. This work was also supported by NASA ATP grant 80NSSC18K1093. 



