Noiseless Gravitational Lensing Simulations
Abstract
The microphysical properties of the DM particle can, in principle, be constrained by the properties and abundance of substructures in DM halos, as measured through strong gravitational lensing. Unfortunately, there is a lack of accurate theoretical predictions for the lensing signal of substructures, mainly because of the discreteness noise inherent to body simulations. Here we present RecursiveTCM, a method that is able to provide lensing predictions with an arbitrarily low discreteness noise, without any free parameters or smoothing scale. This solution is based on a novel way of interpreting the results of body simulations, where particles simply trace the evolution and distortion of Lagrangian phasespace volume elements. We discuss the advantages of this method over the widely used cloudincells and adaptivekernel smoothing density estimators. Applying the new method to a clustersized DM halo simulated in warm and cold DM scenarios, we show how the expected differences in their substructure population translate into differences in the convergence and magnification maps. We anticipate that our method will provide the highprecision theoretical predictions required to interpret and fully exploit strong gravitational lensing observations.
keywords:
gravitational lensing: strong – gravitational lensing: weak – cosmology: theory – dark matter – largescale structure of the Universe – methods: numericalAugust 9, 2019
1 Introduction
Gravitational Lensing has become a powerful and robust technique to explore the “dark side” of our Universe (see Bartelmann, 2010, for a recent review). In the near future, it is expected to probe the accelerated cosmic expansion, and to constrain the properties of the dark matter (DM) particle.
In the weak regime, lensing by the largescale structure of the Universe causes
small distortions in the apparent shape of galaxies and in the temperature
fluctuations maps at the last scattering surface. This effect can be detected
statistically by future widefield surveys (e.g
DES
In the strong regime, efficient lensing configurations can produce multiple images of the same background galaxy or quasar. Each of these images is further distorted by intervening smallscale structures, thus the differences in their shape and/or flux can be used to constrain the substructure content of galaxy and cluster haloes (Mao & Schneider, 1998; Metcalf & Madau, 2001; Dalal & Kochanek, 2002; Kochanek & Dalal, 2004; Natarajan & Springel, 2004; Natarajan et al., 2007; Vegetti et al., 2010). This method is, in fact, the only way of detecting substructures in distant galaxies (Vegetti et al., 2012). The amount and compactness of halo substructures depends strongly on the nature of the DM particle: colder candidates produce more and denser substructures (e.g. Moore et al., 1999; Klypin et al., 1999; Diemand et al., 2007; Springel et al., 2008); particles with larger selfinteraction crosssections produce shallower and more spherical density profiles (Meneghetti et al., 2001; Peter et al., 2013). Therefore, strong gravitational lensing can probe the microphysical properties of the DM particle and thus provide a direct test of the Cold Dark Matter (CDM) paradigm.
In order to fully exploit gravitational lensing measurements in both strong and weak regimes, it is essential to have accurate predictions for the nonlinear state of the mass distribution in the Universe. In particular, for the abundance, spatial distribution and internal properties of DM halos and their substructure. Among the different theoretical approaches available, cosmological body simulations appear as the only robust and accurate method that meets these requirements.
Moreover, cosmological body simulations (e.g. Peebles, 1971; Efstathiou & Eastwood, 1981; Efstathiou et al., 1985; Springel et al., 2005; Angulo et al., 2012) are also invaluable cosmological tools: i) They are the most reliable and precise method to follow the highly nonlinear evolution of primordial density fluctuations (e.g. Kuhlen et al., 2012, for a recent review). ii) They provide virtual universes with which we can test, predict, and interpret astronomical observations (e.g. Overzier et al., 2013). iii) They allow us to experiment with the laws of physics and the background cosmological model (e.g. Fontanot et al., 2012, 2013). Thus, numerical simulations not only can provide the theoretical predictions required by gravitational lensing, but also can be particularly useful for testing analysis algorithms and for exploring the connection between lensing observations and the underlying cosmological model (Bartelmann et al., 1998; Jain et al., 2000; Vale & White, 2003; Meneghetti et al., 2007; Hilbert et al., 2009).
Unfortunately, numerical simulations have a serious limitation that is inherent to the formulation of the body problem: in order to efficiently solve the PoissonVlasov equation, the initial cosmic density field must be represented by a set of discrete bodies. This discretization is necessary to efficiently follow the dynamics and evolution of the DM fluid, but it introduces a smallscale noise that is very often larger than the smallscale lensing signal itself. The noise decreases on large scales and/or with better mass resolution. However, it is still comparable to the strong lensing signal from most of the substructure population, even with the highest resolution simulations to date (Xu et al., 2009; Rau et al., 2013). In other words, the substructure lensing properties that could allow us to constrain the DM particle mass remain buried beneath the discreteness noise. Hence, current theoretical predictions are not sufficiently accurate for upcoming lensing measurements.
In this paper, we propose RecursiveTCM
We devote this paper to the presentation and testing of the algorithm. We start in Section 2 by describing how we compute the gravitational lensing signal of a set of simulation particles. We then apply our method to a clustersize halo simulated in CDM and WDM cosmologies. These simulations are described in Section 3. In Section 4, we compare our method with standard approaches, and show how the noise in the surface density and magnification maps is greatly reduced. This allows us to explore the impact of substructure on the strong lensing magnification fields for our CDM and WDM haloes. We present our conclusions and a discussion of possible future work in Section 5.
2 Lensing Simulations
We start by describing how the gravitational lensing signal of a set of simulation particles is computed, including details of our method to estimate the respective surface density maps.
2.1 Gravitational lensing
Within the plane lens approximation, the lensing distortions produced by a concentrated mass distribution can be derived from a lensing potential, , (e.g. Bartelmann & Schneider, 2001)
(1) 
where denotes an angular position on the (plane) sky, and the convergence is defined as
(2) 
Here, denotes the projected angular surface mass density of the lens mass concentration. The critical angular surface mass density is defined as
(3) 
with the speed of light , gravitational constant , scale factor at lens redshift, and comoving angular diameter distances , , and from the observer to the lens, from the observer to the source, and between the source and the lens, respectively.
The deflection angle , the complex shear , and the magnification , are given by
(4)  
(5)  
(6) 
where the subscripts refer to partial derivatives with respect to one of the angular coordinates.
There are several ways of computing the lensing signal from numerical simulations (e.g Wambsganss et al., 1998; Couchman et al., 1999; Jain et al., 2000; Aubert et al., 2007; Hilbert et al., 2009). Here we choose one of the simplest, which consists in computing the surface density on a regular lattice and then solving for the lensing potential in Fourier space:
(7)  
(8) 
where the superscript indicates a Fourier transform. These
expressions can be readily evaluated by using Fast Fourier Transforms (FFT).
However, this requires additional corrections, because FFT algorithms implicitly
assume periodic boundary conditions, while the appropriate conditions should be
vacuum boundary conditions. To suppress shear artefacts induced by periodic
images of the mass distribution, generous zero padding is employed.
2.2 RecursiveTCM: A new density estimator
The problem is now reduced to obtain the surface mass density, , on a uniform grid from which the respective lensing potential can be computed. Essentially, this step consists in mapping a threedimensional (3D) distribution of simulation particles onto a twodimensional (2D) grid. Although it is in principle a simple task, in practice it is rather difficult to accurately carry out the mapping. Several authors have explored different projection methods and have concluded that all of them gave rise to a noise field of amplitude comparable to the strong lensing signal produced by real DM substructures (Bradač et al., 2004; Li et al., 2006; Xu et al., 2009; Rau et al., 2013). This is true even for the highest mass resolution simulations of DM haloes available to date. Similarly, largescale body simulations, with volumes comparable to that of future widefield galaxy surveys, have typically low number densities of simulation particles, which adds a Poisson shotnoise that dominates the smallscale predictions for weak gravitational lensing (e.g. Jain et al., 2000; Vale & White, 2003; Sato et al., 2009; Hilbert et al., 2009).
There are several proposed ways of dealing with this problem. For strong lensing, one of the most common is to model the smooth component of a DM halo with an analytic expression (e.g. a Single Isothermal Sphere), and then add on top the substructure population (e.g. Xu et al., 2009). Although, it is possible to incorporate the correct density profile and the triaxiality of the DM halo, this method washes out all other higherorder or more subtle features of the DM halo substructures such as streams, caustics, etc. For weak lensing, maps are often smoothed with a fixedsize kernel, which decreases the particle noise but also erases actual smallscale structure (Hilbert et al., 2009; Takahashi et al., 2011).
Here, we discuss RecursiveTCM, a massdepositing scheme that captures a simulated density field in all its complexity without the noise introduced by the finite number of particles, and without any smoothing parameter. The method extends the techniques proposed by Abel et al. (2012); Hahn et al. (2013); Kaehler et al. (2012); Angulo et al. (2013), and thus we refer to them for an extensive discussion of the method. The key idea is to consider simulation particles as vertices of Lagrangian phasespace tetrahedra. At any redshift, the particles indicate the current positions of these vertices. To create surface density maps, the matter represented by these tetrahedra is distributed on a target mesh using a recursive splitting scheme.
One way of interpreting our method is that it assigns to each particle a smoothing kernel whose size and shape are given by their Lagrangian (not Eulerian as in most smoothing methods) neighbours. In particular, this kernel is anisotropic and not even uniquely defined in an Eulerian space. We also note that our method is conceptually different to those that project a Delaunay or Voronoi tessellation built from the Eulerian particle distribution (Schaap & van de Weygaert, 2000; Bradač et al., 2004).
The four main steps of RecursiveTCM are:

1) Creating the Initial Tessellation: First, we need to define a set of disjoint Lagrangian phasespace elements that fully fill the volume of a body simulation. In three dimensions, the most natural choice is a Delaunay tessellation of the unperturbed particle distribution. The result is a set of tetrahedra (six times more abundant than the number of particles) whose corners are given by the simulation particles.
^{13} The connectivity of each tetrahedra is fixed and stored (it can also be trivially recovered from the particles’ ID number in case the ID is related to the position of a particle in an unperturbed lattice). 
2) Reconstructing the Evolved Tessellation: After the simulation has been evolved and the particles moved to different locations, the initial set of tetrahedra (which therefore also moved) is reconstructed using the stored connectivity. The internal density of each tetrahedron is assumed to be uniform, and the density field at any given location is simply given by all those tetrahedra that intersect the target location. We note that it is also possible to compute, at any point of space, other quantities besides the density, such as the number of streams, the velocity dispersion tensor, vorticity, etc. (Hahn et al in prep).

3) Projecting the Density Field: The next step is to compute the projected density field on a grid, i.e, to map the tetrahedra onto a 2D regular mesh. The simplest way, called TCM by Hahn et al. (2013), is to represent each tetrahedron by a single point mass located at the center of mass. Another option is to represent each tetrahedra by 4 particles, preserving the monopole and quadrupole of the parent polyhedra (Hahn et al., 2013). Here we propose a more exact deposit, referred to as RecursiveTCM, which consists in recursively biparting each tetrahedron along its longest edge. The process continues until all the child tetrahedra are completely contained inside one grid cell, or a maximum number of levels in the recursion is reached. Then, each child tetrahedra is subsequently represented using a single particle of mass (where is the recursion level and is the mass of the top tetrahedron) that is deposited using a Nearest Grid Point (NGP) assignment scheme.

4) Removing Density Biases: Over the range of scales in which the mass resolution of a given simulation is adequate to describe the evolution and distortion of Lagrangian phasespace volumes, our method provides a very reliable proxy for the density field (Abel et al., 2012). However, tetrahedronbased density calculations are biased if the distortion of an initial phasespace volume can not be represented by linear transformations. This happens in two situations. One is at the center of DM haloes, which have high densities and short dynamical times. As discussed in Hahn et al. (2013), this has the net effect of densities being overestimated at the halo center, and underestimated at slightly outer regions. The second situation regards the tidal stripping of substructures, where some vertices of a given tetrahedron are stripped while others might still be attached to the substructure. This has the net effect of underestimating the mass associated to substructures.
Fortunately, these biases in the density are small and can be identified and corrected for. Moreover, the centers of haloes are typically dominated by a stellar component (specially in galaxysized DM haloes, where the observational search for substructures is focused), and also are affected by baryonic processes absent in DMonly body simulations (such as feedback, adiabatic compression, etc). Hence, any DMonlybased predictions need to be altered to account for these and produce realistic lensing efficiencies (e.g. Xu et al., 2009), so an additional correction that remove biases of our density estimator can be easily included.
Here we propose and use a simple way to remove the bias:

We first compute an unbiased 2D density map using a traditional (noisier) estimator applied to the same object.

We then apply a correction factor, defined as the average ratio between densities computed using a traditional estimator and the RecursiveTCM, in bins of RecursiveTCM densities. This aims to correct the overestimation of central densities.

We apply an additional correction factor to account for spatially coherent, largescale biases related to tidal stripping of substructures. This extra factor is set to the ratio of the density maps using the traditional and the RecursiveTCM estimator (after the above correction is applied), both Gaussian smoothed to keep only largescale modes.
As we will show in the next section (cf.Section 3.4), this simple correction procedure eliminates most biases in the surface density maps, preserving the reduced noise properties. We note, however, that more sophisticated correction methods are possible and should decrease the biases even further. Some possible extensions are applying corrections to 3D densities instead of projected ones, and/or applying separate correction factors for different substructures.

3 RecursiveTCM in action
For illustrative purposes, we now apply our new method to numerical simulations of cold and warm DM cosmologies. We start by presenting these body simulations, together with one particular DM halo on which we focus our analysis. Then, we provide details of our density estimator when applied to these simulated objects.
3.1 Parent body runs
We employ two of the cosmological body calculations presented in Angulo et al. (2013), that simulate two different cosmological scenarios: i) a standard CDM, and ii) a WarmDM (WDM) model with a DM particle mass. In the latter, fluctuations below are suppressed, which translates into a lack of collapsed structures below , and consequently into a strong suppression of the subhalo population of massive haloes. Although this WDM model is ruled out by observations (Viel et al., 2013), we will consider it for illustrative purposes. The cosmological parameters for the simulations are consistent with the WMAP7 data release (Komatsu et al., 2011): , , , , and spectral index .
Each of these two simulations corresponds to a cubic region of side length, containing simulation particles of mass . The initial conditions were created using the MUSIC code (Hahn & Abel, 2011) at . The particles were subsequently evolved using a TreePM method, as implemented in the LGadget3 code (Angulo et al., 2012; Springel et al., 2005). Gravitational forces are smoothed using a Plummerequivalent softening length set to . Additionally, we have located DM haloes using a FoF algorithm (Davis et al., 1985) (using a standard value for the linking length ), and identified selfbound substructures (or subhalos) within these haloes using the SUBFIND algorithm (Springel et al., 2001).
The numerical simulations were started using identical phases and evolved with the same numerical parameters, which allows a direct comparison of structure formation in general, and of the gravitational lensing signatures explored in our paper.
3.2 Target clustersized DM halo
For our strong gravitational lensing analysis, we will focus on the most massive
cluster present in our simulations at . This object has a mass of , and it is resolved with more than million
particles. For comparison, this corresponds roughly to the lowest resolution
runs of the clusters in the Phoenix project (Gao et al., 2012), and it is a factor
of coarser than the main cluster employed by Rau et al. (2013).
The spherically averaged density profile of the halo is well fitted by a cored NFW profile (Navarro et al., 1996, 1997), , with concentration , and core radius, , both in CDM and WDM. We note that the core is a numerical artifact, and arises due to the lack of force resolution on scales smaller than the simulation softening length.
Fig. 1 shows an image of the selected halo in our two cosmological scenarios. The DM halo displays the same overall structure in WDM and in CDM, and the differences caused by the DM particle mass are evident only in their smallscale properties. In CDM, the halo contains a wealth of substructure: a large number of small clumps that are the remnants of previously accreted DM haloes. These, in contrast, are almost absent when WDM is adopted, but caustics, streams and filaments are instead much better defined. Inside of the CDM halo, we find substructures with more than particles, which corresponds to a minimum subhalo mass of . Contrasting this, we found only substructures inside the WDM halo – which are mostly a result of numerical fragmentation of filaments (Wang et al., 2007; Angulo et al., 2013). The substructure population contributes and of the mass inside , respectively for our WDM and CDM halo.
Considering only substructure with masses above , the subhalo mass function in CDM follows a power law . However, the slope decreases to when we consider all the subhalos detected. These values are shallower that the average slope found in other simulations (0.9, e.g. Angulo et al., 2008; Gao et al., 2012). The discrepancy is most likely caused by our low force resolution compared to our mass resolution (many lowmass haloes are tidally disrupted too efficiently due to our low force resolution, which makes our subhalo mass function being incomplete up to higher subhalo masses than in simulations with smaller softening lengths), which could explain the change in slope in the subhalo mass function. Although these discrepancies are not important for our work, we caution the reader that the amount of substructure in our work is underestimated compared to other simulations of similar mass resolution.
3.3 RecursiveTCM lensing simulations
We are now in position of applying our method to the WDM and CDM halo described before. We artificially place the haloes at , where the most massive galaxy clusters are expected to be observed, and assume a background source population located at . We consider the 3D particle distribution inside a region of dimensions (equal to ) and deep centred on our halo, and project it onto a pixels mesh. This yields a spatial resolution of , sufficient to resolve the smallest structures present in our simulation given our gravitational resolution limit ().
We apply our full RecursiveTCM algorithm using a maximum of ten recursion levels, i.e. every toplevel tetrahedron is split into smaller tetrahedra, at most. Using these maps, we create convergence fields, compute the lensing potential, and derive the associated , , and , as described in Section 2.1. We use a FFT mesh (this is approximately a factor of thousand larger than the density mesh to allow for nonperiodic boundary conditions), and compute the spatial derivatives in Fourier space. We note that by considering only a small inner region, we are neglecting the weak lensing effects of structures further away from the cluster center. However, we have explicitly checked that this is a good approximation for the quantities explored here.
The computational cost of our RecursiveTCM algorithm is higher than usual projection methods, but it is still negligible compared to the time employed in carrying out an body simulation. For our particular data structure, data access and target grid, and maximum recursion level, it took minutes using processors, i.e. CPU hours. It is important to remember that this corresponds to an implementation in software for CPUs, and that less recursion levels significantly reduce the execution time. In addition, alternative algorithms based on GPUrendering routines can, in principle, achieve significantly better performances (Kaehler et al., 2012).
For comparison, we computed densities using two other techniques, in addition to RecursiveTCM. The first one, referred to as CiC, represents each particle by a cube of uniform density and size equal to the cell size of the target grid. This is the mostcommon projection method in cosmology (Hockney & Eastwood, 1981). The second method, referred to as Sph, employs a sphericallysymmetric polynomial kernel (Springel et al., 2005) to project each particle onto our 2D grid. The characteristic size of the kernel is given by the local density about each particle, explicitly, by the 3D distance to the 32th nearest neighbour. This approach is the core of the SmoothedParticleHydrodynamics (SPH, Monaghan, 1992) numerical formalism.
3.4 Bias correction
As discussed in Section 2.2, RecursiveTCM densities can be biased in regions where heavy winding of the primordial phasespace sheet occurs. Fortunately, we can use a noisy estimator to correct for such biases, as we will show. In practice, we follow the procedure described in §2.2 and apply a twostep bias correction.
The first step corrects for the overestimation of central densities by an average factor, shown in Fig. 2. This factor was computed as the geometric mean of the ratio between densities estimated using Sph and RecursiveTCM, in logarithmic bins of . Because the densities produced by our tetrahedral approach are biased high in central regions of DM haloes (see Section 2.2), the average ratio progressively decreases at higher densities. The largest correction factor is at the very center of our halo.
The second correction step accounts for spatially coherent biases. In Fig. 3, we show the ratio between convergence maps in Sph and in RecursiveTCM for our WDM halo and after the first correction has been applied. Ideally, this figure would be a pure whitenoise field, however, we can clearly see that there is a largescale component in this field. This arises partially because the 2D projected density field is not a onetoone function of the full 6D phasespace structure (which determines the amount of winding and overestimation). Another aspect contributing to this map is related to mass accretion history and shortcomings of the tetrahedronbased densities to represent the tidal stripping of infalling DM halos. In order to correct for this, we further divide the RecursiveTCM map by the a version of the map shown in Fig. 3, but smoothed with a Gaussian kernel of size . We note that this scale is set to be larger than those dominated by discreteness noise in Sph. As we will see later, the simple procedure described here, is successful in producing accurate lensing maps.
4 Results
We now present and discuss the results of our lensing simulations. We first focus on the input surface density maps, and then on the lensing magnification. In particular, we discuss the performance of our algorithm and the role of discreteness noise compared to the signal of DM substructures.
4.1 The surface densities and lensing convergence
Fig. 4 shows maps of the surface density in the inner regions of our WDM (top row) and CDM haloes (bottom row). Each column shows the result of one of our three projection methods, as indicated by the legend. Note that the colour scale is identical in all six panels.
In both cosmological scenarios, we can appreciate how the smallscale noise is decreased from left to right. The CiC method shows the largest fluctuations, though we note that the visual impression of the noise depends on the target mesh resolution, as this sets the size of the CiC smoothing kernel. The Sph method reduces the noise significantly in this case, though still a considerable amount of spurious fluctuations remains. These two images illustrate the discretisationrelated noise in traditional density estimators.
In contrast, the new RecursiveTCM method, does not present any appreciable noise thanks to the extra sophistication in the mass deposit, nor presents appreciable biases thanks to the correction method described in Section 2.2. We note that despite being much smoother, all those peaks seen clearly in the CiC and Sph maps also appear in the RecursiveTCM maps. It is important to note that our method is the only one that could in principle distinguish random fluctuations in surface density maps from those produced by halo substructures: compare the differences between the CDM and WDM halo in rightmost column, with the differences seen in the leftmost column. In both SPH and CiC it is almost impossible to visually distinguish CDM from WDM.
Fig. 5 shows the power spectrum of the surface overdensity, , as given by the three projection methods applied to the WDM halo. We note that the measurements are robust only until , where denotes the Nyquist frequency of the surface density mesh. On smaller scales, aliasing and the mass assignment window introduce visible artefacts, damping the measured power spectra. We also display the expectations of a whitenoise field with the same number of point particles as those projected in our surface density maps. For comparison, we also show the expectations for a (cored) NFW halo with the same mass and concentration as the sphericallyaveraged density profile of our DM haloes, but without noise.
Comparing all the measured power spectra, we observe a situation consistent with the visual impression provided before. On large and intermediate scales, all methods provide essentially identical power spectra decreasing as , as expected for an NFW profile. This incidentally supports the validity of our simple approach to correct the biases in RecursiveTCM densities.
On scales smaller than or roughly () – a scale much larger than the typical size of substructures – all methods begin to differ. The CiC spectrum follows the value expected for a 2D Poisson field: . Interestingly, the Sph method levels to this expectation at roughly the same scale as the CiC spectrum, but then decays much more quickly. This is because the Sph method heavily smooths the field on scales smaller than the kernel size. This smoothing also erases true (specially anisotropic) smallscale density fluctuations present in our DM haloes. For instance, it can be trivially seen that structures resolved with less than particles will be smoothed out.
The performance of RecursiveTCM appears to be considerably better. The noise level is a few orders of magnitude below that of the other projection methods. The noise measurements are consistent with our expectations of reducing the noise in a way proportional to the maximum level of recursion, , in the adaptive mass deposit: , where is the number density of bodies used in the map creation. The prefactor of six is understood in terms of the six times more particles (one per tetrahedron) employed to describe the mass field. This scaling predicts the Poisson noise in our measurements to be a factor of smaller than the CiC method. This value is very close to the actual differences (measured at ): .
Finally, we can see that RecursiveTCM creates a projected mass power spectrum that is very close to that of an ideal NFW halo, differing only on the slope at high wavenumbers. On these scales, the core introduced by the softening length in our simulations becomes relevant, and the RecursiveTCM power spectrum follows that of a cored NFW profile.
4.2 Mass resolution dependence
We now explore how the noise of our convergence maps vary with the mass resolution of the underlying body simulation. For this, we have downsampled the field by factors of , and in each coordinate, or equivalently, reducing the total number of particles in our body simulations by factors of , , and . The most downsampled case is equivalent to a particle simulation, where our WDM halo would be resolved with only particles. For each case, we have created convergence maps with maximum recursion levels set at for the original maps, and to , and , respectively for the downsampled versions. The increased recursion level compensates the sparser particle data, so that in all four cases the maps are created with roughly the same number of tetrahedra (i.e. keeping constant).
In Fig. 6 we show the four resulting convergence maps. In all subpanels we see that our method produces extremely smooth surface density maps. Naturally, as we decrease the effective resolution, smallscale features slowly disappear, for instance, the three substructures located on the bottomright side of the image. With low mass resolution, orbits and accretion become discrete and subtle radial features appear. However, even in the lowerright case, which contains almost orders of magnitude less particles than in our original simulation, the smallscale noise is considerably smaller than with the CiC method (compare to the leftmost panel of the top row in Fig. 4). This shows that the limitation of our maps is in the actual amount of structure that the parent body simulation tracks correctly, and not in the discreteness noise associated with the finite number of particles.
A quantitative comparison can be obtained from Fig. 7, which shows the 2D power spectra of our four test cases, with and without applying our density correction (cf. Section 3.4). The overall shape of the power spectra is very similar among all resolutions, as expected from the previous images, but small differences arise due to the different amount of structure resolved in the different cases. Nevertheless, the discreteness noise appears at the same level and is set by the maximum amount of deposited tetrahedra. This again shows that the limitation of our gravitational lensing simulations reside on the ability of the parent body simulation to represent DM structures properly, and not in an artificial noise introduced by our massprojection method. In the next subsection we will show that this has positive repercussions on simulated lensing signals.
4.3 The lensing magnification
The magnification field, , which gives the ratio of the area of the lensed image to the original area of the source, depends on secondorder derivatives of the lensing potential (whereas depends only on firstorder derivatives). Thus, is very sensitive to small perturbations to the lensing potential (such as those caused by subhalos). This is also the reason why is very susceptible to noise introduced by discretisation in body simulations.
In Fig. 8 we display maps of the inverse of the magnification field, , created from each of the three projection methods we consider. We highlight two places i) where the magnification is formally infinite, , which are known as critical lines; and ii) where , which are useful to explore the impact of noise and substructures in the topology of the magnification field. Note that our simulated cluster is not a particularly efficient lens, partly due to its particular dynamical state, and also because of our modest force resolution and the lack of a modelling of a central stellar component.
While the magnification maps from CiC, Sph, and RecursiveTCM agree on large scales, they differ substantially in the amount of associated smallscale noise. In the CiC case, the noise fluctuations make it almost impossible to distinguish CDM from WDM based solely on the topology of iso lines. The same is true to some degree in for Sph. The method RecursiveTCM is superior: The contours are not disturbed by discreteness noise. This allows us to explore the magnification field with great detail. When applied to the WDM case, we see contour lines that are extremely smooth. For the CDM halo, the contour lines are also very smooth in most parts of the map, but display many small protuberances with high significance. As we will see in the next section, these are caused by the rich substructure population of this halo and thus are a distinctive signature of the DM candidate properties.
Figure 9 shows frequency of magnification values predicted by the different methods for the WDM (results for the CDM halo are very similar). The noise in the CiC magnification maps leads to substantially broader magnification distributions compared to those for Sph and RecursiveTCM. In contrast, the distributions predicted by Sph and RecursiveTCM are very similar. However, the probability densities predicted by Sph display more features, i.e. local deviations from a smooth density function, which we attribute to residual noise in the Sph maps.
4.4 The impact of substructures
In order to further explore the capabilities of our method, we now focus on the impact of halo substructure on the magnification and convergence maps. In Fig. 10 we show isomagnification lines, as computed in RecursiveTCM maps, together with the substructures identified with more than particles in our simulated halo. Note that the small amount of noise visible in the outer contour is a result of the maximum number of recursion levels () employed in our method. This noise can also be seen in the power spectra of the projected density field shown in Fig. 5, and, as we discussed earlier, it can be reduced further by simply increasing at the expense of more CPU time.
We can clearly see the differences between CDM and WDM. In the WDM case, there are only substructures in the field, and consequently iso lines are mostly smooth, showing only a few notorious protuberances. In contrast, in the CDM case, there are substructures, and the iso lines show many protuberances, but are almost smooth otherwise.
As Fig. 10 illustrates, all protuberances in the isomagnification contours are associated with nearby substructures. However, the relation is not simple, and different substructures produce perturbations of different importance. Moreover, in some cases fluctuations are not caused by a single substructure, but by a group of them. This case is seen, for example, in the lower right section of the CDM . On the contrary, some substructures near contours do not strongly perturb the magnification field. These objects have typically surface mass densities below average. For instance, the substructure located at in the WDM case, has a projected density a factor of ten smaller than the substructure found at .
In Fig. 11 we compare the signal of identified substructures among the convergence maps. The xaxis shows the true subhalo mass . The yaxis shows the ratio of the excess convergence with respect to a simple expectation based on the substructure properties. The measured value, , is defined as , where is the mean convergence within the halfmass radius , and the background convergence is the given by the mean convergence in an annulus with . The expected value, , is defined as .
The deviations of from unity seen in Fig. 11 can be explained, e.g., by particle noise (though not for RecursiveTCM), projections effects, triaxiality, or inaccuracies in the background estimation. Note that the scatter in is much lower in our method than in the other two. This is thanks to a less noisy estimation of both the signal itself and the background. However, the average value of is roughly a factor of two smaller in RecursiveTCM than in the other two methods. The discrepancy is originated by two factors. First, it is due to an overestimation of the bias correction factor: since this factor is essentially set by the background halo, it does not capture the exact density biases for substructures, which have different central densities and dynamical times. The second aspect is an intrinsic underestimation of the mass associated to subhalos in RecursiveTCM. This has an origin in tetrahedra being stretched along a subhalo’s orbit by tidal stripping. We estimated this effect to cause an underestimation of about 30% in the mass inside subhalos for resolved substructures in our CDM halo. It remains to be explored whether more sophisticated correction procedures, or modifications to the RecursiveTCM algorithm, will alleviate these discrepancies.
In order to quantify the performance of RecursiveTCM concerning substructure lensing signals, we have implemented a hierarchical peak finder algorithm. This proceeds as follows: We start by smoothing the convergence field with a Gaussian kernel of size and then identify local peaks in the smoothed field. Then, we progressively reduce the kernel size and search for new peaks, discarding those that are inside a larger peak. We repeat this procedure for different scales, uniformly spaced in , down to . Finally, we compute the signal associated to each peak as , where is the average convergence within and is the local background value defined as the average of the map in an annulus of .
The results are shown in Fig. 12, which displays the cumulative number of peaks detected by our algorithm when applied to CiC, Sph and RecursiveTCM maps, as a function of their local convergence excess . In addition, we plot as a black line the substructures detected in 3D by SUBFIND. The associated is computed in the same way as that of our peaks, but using as the peak scale. Our algorithm finds (CiC), (Sph) and (RecursiveTCM) peaks in the CDM map, and (CiC), (Sph) and (RecursiveTCM), peaks in the WDM map.
It is clear that the CiC and Sph maps contain a large amount of spurious peaks produced by the discreteness noise. At all there are between one and two orders of magnitude more detected peaks than real substructures. Moreover, the number of detected peaks is almost identical between CDM and WDM, even though the actual amount of substructure is very different. This further exemplifies that in current lensing simulations the impact of noise is comparable or larger than that of real DM substructures.
In contrast, our method recovers roughly the correct amount of peaks, which is an order of magnitude larger in CDM than in WDM. Furthermore, and of the substructures can be matched to a peak in CDM and WDM, respectively. Among those substructures not detected as peaks, we mostly find objects with a negative or very small value, which are also not detected in the CiC or Sph. This suggests that these might indeed be the falsepositives in SUBFIND or special cases of projection effects. However, it might also be a consequence of the simplicity of our peak detection algorithm. In order to quantify the implications for observational constraints, this issue requires further study considering realistic input signals and analysis procedures, and perhaps more sophisticated procedures to correct for RecursiveTCM biases.
There are also peaks in the convergence maps that are not related to any identified 3D substructure. In many cases, these are due our 3D substructure finder algorithm: a density peak not found by SUBFIND, one that fell below our massresolution, or one that is not a selfbound object. For instance, the large peak located at . An object of a different nature is shown in Fig. 13. It does not correspond to a selfbound spherical overdensity, but to a DM stream of a tidally disrupted substructure. Such features are expected in hierarchical structure formation scenarios, and perhaps they could be eventually detected trough their lensing signal. (Note that this feature is barely distinguished over the noise in Sph). For the moment, this detection serves as a further example of the potential accuracy and precision of the lensing maps created by the method presented and discussed here.
5 Conclusions
The next generation of gravitational lensing observations might help us to decipher the mysteries of the Dark Universe: Dark Energy and Dark Matter. However, highprecision theoretical predictions are essential to ensure the correct interpretation of future datasets.
We presented RecursiveTCM, a method aimed at predicting the lensing signal from cosmological simulations with the required accuracy. This algorithm originates from a novel way of interpreting the results of body simulations and overcomes one of the most serious limitations of current lensing simulations: the noise introduced by the discrete nature of the particle representation of the density field.
We applied RecursiveTCM to cluster sizedhaloes simulated in WDM and CDM universes. We showed that the method produces convergence maps with noise several orders of magnitude below that of traditional methods (here, a factor of smaller than the formal shotnoise limit), and that it recovers the underlying power spectrum of fluctuations well below the particle shot noise of the simulations.
With traditional projection methods, the discreteness noise in lensing maps are comparable to the signal from DM substructure. This is not true for RecursiveTCM, where the features associated with real overdensities are preserved, but the discreteness noise is absent. All this comes without free parameters or additional smoothing scales. We also showed that there are density biases associated to RecursiveTCM, which, however, can be mostly eliminated by simple correction procedures. Therefore, this method is well suited for creating highprecision predictions for the relation between the underlying cosmological model and the expected signatures of smallcale structure in strong gravitational lensing observations.
With RecursiveTCM, we were able to clearly show the differences in the lensing properties between CDM and WDM. The differences come mainly from their substructure population, and thus lensing might be able to constrain the DM particle mass. However, we found that the relation between substructures and the associated lensing effects is not trivial: some substructures do not affect the convergence noticeably; many lensing perturbations are caused by more than a single structure; and a few perturbations are not associated with any selfbound substructure, but, e.g. with tidal debris. This suggests that in order to interpret correctly the lensing measurements of substructures, a rigorous study of the detectability of substructures needs to be carried out.
In this paper, we have shown the feasibility of our method, providing examples of its potential when applied to a rather modest simulation. In the future, we expect our method to enable many detailed theoretical studies, exploiting stateoftheart simulations of much higher force and mass resolution, also simulating more realistic WDM scenarios, and even taking advantage of hydrodynamical simulations. The presented method will also be very useful for creating largescale weak lensing shear and magnification maps with high fidelity and low particle noise. Moreover, the method will be crucial for testing and characterising the performance of algorithms of extracting substructure information from observed lensed galaxies, in particular methods for constraining the DM particle mass from image perturbations or fluxratio anomalies in multipleimage systems. All this together will allow us to understand better the impact of the underlying cosmological model in lensing observations, and therefore help to unleash the full potential of gravitational lensing.
Acknowledgements
We warmly thank Oliver Hahn and Ralf Kaehler for many enlightening discussions. We also thank Simon White for valuable suggestions that improved our paper. We acknowledge useful comments on the manuscript from Carlo Giocoli, Yashar Hezaveh, Phil Marshal, Ben Metcalf, Stefan Rau and Simona Vegetti. T.A. gratefully acknowledges support by the National Science foundation through award number AST0808398 and the LDRD program at the SLAC National Accelerator Laboratory as well as the Terman Fellowship at Stanford University. We gratefully acknowledge the support of Stuart Marshall and the SLAC computational team, as well as the computational resources at SLAC. Parts of the computations were carried out at the Rechenzentrum Garching (RZG) of the Max Planck Society.
Footnotes
 pagerange: Noiseless Gravitational Lensing Simulations–Noiseless Gravitational Lensing Simulations
 pubyear: 2013
 http://www.darkenergysurvey.org/
 http://jpas.org/
 http://sci.esa.int/euclid/
 http://www.lsst.org/
 http://www.esa.int/Our_Activities/Space_Science/Planck/
 http://pole.uchicago.edu/
 http://www.princeton.edu/act/
 The term RecursiveTCM abbreviates for ”Recursive deposit of Tethrahedra approximated by their Center of Mass”.
 The reduction of discreteness noise also helps to suppress the artificial fragmentation of filaments seen in warm Dark Matter (WDM) simulations (Hahn et al., 2013; Angulo et al., 2013).
 For simplicity, we refrain from also applying a forcerange cutoff.
 Constructing the tessellation can be a computationally expensive task for stateoftheart simulations (e.g. Pandey et al., 2013). However, this is trivial if the particle distribution is arranged in a regular lattice (as opposite to a glasslike distribution): each set of 8 grid points defines a cube that is subdivided into six disjoint tetrahedra.
 Our force resolution is also much lower.
References
 Abel T., Hahn O., Kaehler R., 2012, MNRAS, 427, 61
 Angulo R. E., Baugh C. M., Lacey C. G., 2008, MNRAS, 387, 921
 Angulo R. E., Hahn O., Abel T., 2013, arXiv:1304.2406
 Angulo R. E., Springel V., White S. D. M., Jenkins A., Baugh C. M., Frenk C. S., 2012, MNRAS, 426, 2046
 Aubert D., Amara A., Metcalf R. B., 2007, MNRAS, 376, 113
 Bartelmann M., 2010, Classical and Quantum Gravity, 27, 233001
 Bartelmann M., Huss A., Colberg J. M., Jenkins A., Pearce F. R., 1998, A&A, 330, 1
 Bartelmann M., Schneider P., 2001, Phys. Rep., 340, 291
 Bradač M., Schneider P., Lombardi M., Steinmetz M., Koopmans L. V. E., Navarro J. F., 2004, A&A, 423, 797
 Couchman H. M. P., Barber A. J., Thomas P. A., 1999, MNRAS, 308, 180
 Dalal N., Kochanek C. S., 2002, ApJ, 572, 25
 Davis M., Efstathiou G., Frenk C. S., White S. D. M., 1985, ApJ, 292, 371
 Diemand J., Kuhlen M., Madau P., 2007, ApJ, 667, 859
 Efstathiou G., Davis M., White S. D. M., Frenk C. S., 1985, ApJS, 57, 241
 Efstathiou G., Eastwood J. W., 1981, MNRAS, 194, 503
 Fontanot F., Puchwein E., Springel V., Bianchi D., 2013, arXiv:1307.5065
 Fontanot F., Springel V., Angulo R. E., Henriques B., 2012, MNRAS, 426, 2335
 Gao L., Navarro J. F., Frenk C. S., Jenkins A., Springel V., White S. D. M., 2012, MNRAS, 425, 2169
 Giannantonio T., Porciani C., Carron J., Amara A., Pillepich A., 2012, MNRAS, 422, 2854
 Hahn O., Abel T., 2011, MNRAS, 415, 2101
 Hahn O., Abel T., Kaehler R., 2013, MNRAS
 Hilbert S., Hartlap J., White S. D. M., Schneider P., 2009, A&A, 499, 31
 Hilbert S., Marian L., Smith R. E., Desjacques V., 2012, MNRAS, 426, 2870
 Hockney R. W., Eastwood J. W., 1981, Computer Simulation Using Particles. Computer Simulation Using Particles, New York: McGrawHill, 1981
 Huterer D., 2010, General Relativity and Gravitation, 42, 2177
 Jain B., Seljak U., White S., 2000, ApJ, 530, 547
 Kaehler R., Hahn O., Abel T., 2012, arXiv:1208.3206
 Klypin A., Kravtsov A. V., Valenzuela O., Prada F., 1999, ApJ, 522, 82
 Kochanek C. S., Dalal N., 2004, ApJ, 610, 69
 Komatsu E. et al., 2011, ApJS, 192, 18
 Kuhlen M., Vogelsberger M., Angulo R., 2012, Physics of the Dark Universe, 1, 50
 Li G.L., Mao S., Jing Y. P., Kang X., Bartelmann M., 2006, ApJ, 652, 43
 Mao S., Schneider P., 1998, MNRAS, 295, 587
 Marian L., Hilbert S., Smith R. E., Schneider P., Desjacques V., 2011, ApJ, 728, L13
 Meneghetti M., Argazzi R., Pace F., Moscardini L., Dolag K., Bartelmann M., Li G., Oguri M., 2007, A&A, 461, 25
 Meneghetti M., Yoshida N., Bartelmann M., Moscardini L., Springel V., Tormen G., White S. D. M., 2001, MNRAS, 325, 435
 Metcalf R. B., Madau P., 2001, ApJ, 563, 9
 Monaghan J. J., 1992, ARA&A, 30, 543
 Moore B., Ghigna S., Governato F., Lake G., Quinn T., Stadel J., Tozzi P., 1999, ApJ, 524, L19
 Natarajan P., De Lucia G., Springel V., 2007, MNRAS, 376, 180
 Natarajan P., Springel V., 2004, ApJ, 617, L13
 Navarro J. F., Frenk C. S., White S. D. M., 1996, ApJ, 462, 563
 Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ, 490, 493
 Oguri M., Takada M., 2011, Phys. Rev. D, 83, 023008
 Overzier R., Lemson G., Angulo R. E., Bertin E., Blaizot J., Henriques B. M. B., Marleau G.D., White S. D. M., 2013, MNRAS, 428, 778
 Pandey B., White S. D. M., Springel V., Angulo R., 2013, arXiv:1301.3789
 Peebles P. J. E., 1971, A&A, 11, 377
 Peter A. H. G., Rocha M., Bullock J. S., Kaplinghat M., 2013, MNRAS, 430, 105
 Rau S., Vegetti S., White S. D. M., 2013, MNRAS, 430, 2232
 Sato M., Hamana T., Takahashi R., Takada M., Yoshida N., Matsubara T., Sugiyama N., 2009, ApJ, 701, 945
 Schaap W. E., van de Weygaert R., 2000, A&A, 363, L29
 Shandarin S., Habib S., Heitmann K., 2012, Phys. Rev. D, 85, 083005
 Springel V. et al., 2008, MNRAS, 391, 1685
 Springel V. et al., 2005, Nature, 435, 629
 Springel V., Yoshida N., White S. D. M., 2001, New Astronomy, 6, 79
 Takahashi R., Oguri M., Sato M., Hamana T., 2011, ApJ, 742, 15
 Vale C., White M., 2003, ApJ, 592, 699
 Vegetti S., Koopmans L. V. E., Bolton A., Treu T., Gavazzi R., 2010, MNRAS, 408, 1969
 Vegetti S., Lagattuta D. J., McKean J. P., Auger M. W., Fassnacht C. D., Koopmans L. V. E., 2012, Nature, 481, 341
 Viel M., Becker G. D., Bolton J. S., Haehnelt M. G., 2013, Phys. Rev. D, 88, 043502
 Wambsganss J., Cen R., Ostriker J. P., 1998, ApJ, 494, 29
 Wang H. Y., Mo H. J., Jing Y. P., 2007, MNRAS, 375, 633
 Xu D. D. et al., 2009, MNRAS, 398, 1235