Noiseless Lensing Simulations

Noiseless Gravitational Lensing Simulations


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 Recursive-TCM, 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 phase-space volume elements. We discuss the advantages of this method over the widely used cloud-in-cells and adaptive-kernel smoothing density estimators. Applying the new method to a cluster-sized 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 high-precision theoretical predictions required to interpret and fully exploit strong gravitational lensing observations.

gravitational lensing: strong – gravitational lensing: weak – cosmology: theory – dark matter – large-scale structure of the Universe – methods: numerical

March 6, 2018

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 large-scale 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 wide-field surveys (e.g DES3, J-PAS4, Euclid5, LSST6), and by Cosmic Microwave Background (CMB) experiments (Planck7, SPT8, ACT9). From correlations in the distortions, one can infer the amplitude, shape, and redshift evolution of the matter power spectrum – quantities sensitive to the initial density perturbations, the law of gravity and the cosmic expansion. Therefore, gravitational lensing measurements are expected to contribute significantly to our understanding of the Dark Energy and the physics of the early Universe (e.g. Huterer, 2010; Marian et al., 2011; Oguri & Takada, 2011; Hilbert et al., 2012; Giannantonio et al., 2012).

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 small-scale 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 self-interaction cross-sections 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 Poisson-Vlasov 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 small-scale noise that is very often larger than the small-scale 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 Recursive-TCM10, a method to create gravitational lensing simulations free of discreteness noise, with no tunable parameters nor additional smoothing scales. Our procedure builds on a recently proposed method to solve for the collisionless dynamic of the DM fluid (Abel et al., 2012; Shandarin et al., 2012; Hahn et al., 2013; Kaehler et al., 2012; Angulo et al., 2013). The novel approach considers simulation particles as the vertices of Lagrangian phase-space volume elements, not mass carriers as in the usual interpretation of numerical simulations. The evolution and distortion of these volume elements is described by the Eulerian coordinates of simulation particles. Consequently, the DM density field is determined by the spatially overlapping phase-space elements, these can be exactly deposited onto a target grid using a recursive algorithm. The result is a continuous and smooth density field ideal for small-scale lensing simulations.11

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 cluster-size 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)


where denotes an angular position on the (plane) sky, and the convergence is defined as


Here, denotes the projected angular surface mass density of the lens mass concentration. The critical angular surface mass density is defined as


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


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:


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.12 To recover the correct mean convergence [which is lost in the FFT methods due to setting to zero], the potential from the FFT is corrected by a term . Finally, the lensing deflection, shear, and magnification can be obtained by computing derivatives of either in Fourier space, or in real space by using finite difference methods (Hilbert et al., 2009).

2.2 Recursive-TCM: 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 three-dimensional (3D) distribution of simulation particles onto a two-dimensional (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, large-scale -body simulations, with volumes comparable to that of future wide-field galaxy surveys, have typically low number densities of simulation particles, which adds a Poisson shot-noise that dominates the small-scale 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 higher-order or more subtle features of the DM halo substructures such as streams, caustics, etc. For weak lensing, maps are often smoothed with a fixed-size kernel, which decreases the particle noise but also erases actual small-scale structure (Hilbert et al., 2009; Takahashi et al., 2011).

Here, we discuss Recursive-TCM, a mass-depositing 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 phase-space 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 Recursive-TCM are:

  • 1) Creating the Initial Tessellation: First, we need to define a set of disjoint Lagrangian phase-space 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 Recursive-TCM, 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 phase-space volumes, our method provides a very reliable proxy for the density field (Abel et al., 2012). However, tetrahedron-based density calculations are biased if the distortion of an initial phase-space 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 galaxy-sized DM haloes, where the observational search for substructures is focused), and also are affected by baryonic processes absent in DM-only -body simulations (such as feedback, adiabatic compression, etc). Hence, any DM-only-based 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 Recursive-TCM, in bins of Recursive-TCM densities. This aims to correct the overestimation of central densities.

    • We apply an additional correction factor to account for spatially coherent, large-scale biases related to tidal stripping of substructures. This extra factor is set to the ratio of the density maps using the traditional and the Recursive-TCM estimator (after the above correction is applied), both Gaussian smoothed to keep only large-scale 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 Recursive-TCM 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 Warm-DM (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 Tree-PM method, as implemented in the L-Gadget3 code (Angulo et al., 2012; Springel et al., 2005). Gravitational forces are smoothed using a Plummer-equivalent 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 self-bound 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 cluster-sized 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).14 The force resolution is times smaller than the halo’s virial radius, and thus the halo structure is resolved adequately for our purposes.

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.

Figure 1: Projected DM density, as computed by Recursive-TCM, for the most massive halo in our simulations at (). Each image corresponds to a square region of side length and projection depth around the cluster center. The white circle indicates the virial radius . The left panel shows the halo in a WDM scenario, the right panel assumes a CDM cosmology.

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 small-scale 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 low-mass 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 Recursive-TCM 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 Recursive-TCM algorithm using a maximum of ten recursion levels, i.e. every top-level 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 non-periodic 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 Recursive-TCM 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 GPU-rendering routines can, in principle, achieve significantly better performances (Kaehler et al., 2012).

For comparison, we computed densities using two other techniques, in addition to Recursive-TCM. 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 most-common projection method in cosmology (Hockney & Eastwood, 1981). The second method, referred to as Sph, employs a spherically-symmetric 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 32-th nearest neighbour. This approach is the core of the Smoothed-Particle-Hydrodynamics (SPH, Monaghan, 1992) numerical formalism.

3.4 Bias correction

As discussed in Section 2.2, Recursive-TCM densities can be biased in regions where heavy winding of the primordial phase-space 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 two-step bias correction.

Figure 2: Relation between the densities estimated using Recursive-TCM and Sph, for pixels covering the inner region centred in our WDM cluster. The red line display the mean value in logarithmic bins of .

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 Recursive-TCM, 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.

Figure 3: Fractional differences between the convergence maps built using the Sph and Recursive-TCM methods. The region displayed correspond to the inner centered in our WDM halo.

The second correction step accounts for spatially coherent biases. In Fig. 3, we show the ratio between convergence maps in Sph and in Recursive-TCM for our WDM halo and after the first correction has been applied. Ideally, this figure would be a pure white-noise field, however, we can clearly see that there is a large-scale component in this field. This arises partially because the 2D projected density field is not a one-to-one function of the full 6D phase-space 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 tetrahedron-based densities to represent the tidal stripping of infalling DM halos. In order to correct for this, we further divide the Recursive-TCM 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

Figure 4: The convergence in the central region about our WDM (top row) and CDM clusters (bottom row), using CiC, Sph, or Recursive-TCM density estimates. The colour scale is identical for all six sub-images. Note the different noise levels present in the different projection methods. Recursive-TCM displays the least noise.

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 small-scale 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 discretisation-related noise in traditional density estimators.

In contrast, the new Recursive-TCM 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 Recursive-TCM 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.

Figure 5: The power spectra of the 2D over-density , produced by the three projection algorithms when applied to a WDM halo. Wavelengths are shown in units of the Nyquist frequency of the density mesh . The method proposed here, Recursive-TCM, yields the lowest amplitude on small-scale fluctuations. Also shown are expectations for a smooth NFW profile without and with central core. The dotted horizontal line indicates the white-noise level.

Fig. 5 shows the power spectrum of the surface over-density, , 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 white-noise 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 spherically-averaged 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 Recursive-TCM 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) small-scale 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 Recursive-TCM 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 Recursive-TCM 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 Recursive-TCM 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 down-sampled 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 down-sampled 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 down-sampled 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).

Figure 6: Comparison of the convergence maps created by the same particle distribution of the WDM halo down-sampled by factors of (top-right), (bottom-left) or (bottom-right). Therefore, these would correspond to haloes simulated with , , and , respectively, as indicated by the legend. The original map is shown on the top-left panel. We use an identical logarithmic colour scale in all subpanels.

In Fig. 6 we show the four resulting convergence maps. In all sub-panels we see that our method produces extremely smooth surface density maps. Naturally, as we decrease the effective resolution, small-scale features slowly disappear, for instance, the three substructures located on the bottom-right side of the image. With low mass resolution, orbits and accretion become discrete and subtle radial features appear. However, even in the lower-right case, which contains almost orders of magnitude less particles than in our original simulation, the small-scale 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.

Figure 7: Comparison between the power spectra of projected overdensity maps created from the WDM halo with the Recursive-TCM method applied on progressively sparser data. The solid red line show the result for our original dataset and a maximum level of recursion set to , whereas the dot-dashed magenta line shows the result for a case using less particles but allowing seven further levels of refinements.
Figure 8: Map of the inverse of the magnification field, , at the central region of our WDM (top) and CDM (bottom). The region displayed matches that shown in Fig. 4. White and black lines shows contours where and , respectively. Note we use the same colour scale in all panels, ranging from (white) to (light yellow).

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 mass-projection 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 second-order derivatives of the lensing potential (whereas depends only on first-order 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.

Figure 9: Probability density of the magnification modulus computed with CiC (red dotted), Sph (blue dashed), and Recursive-TCM (black solid lines) for the WDM halo.

Figure 10: The relation between substructure and perturbations in lensing magnification for our simulated CDM (left) and WDM (right) halo. Black lines denote iso-magnification contours at and inwards. Red circles indicate the positions where substructures were identified, and their radii is equal to the half-mass radius of the respective subhalo. Note the reduced number of substructures in the WDM case, which result from the initial suppression of small-scale fluctuations.

While the magnification maps from CiC, Sph, and Recursive-TCM agree on large scales, they differ substantially in the amount of associated small-scale 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 Recursive-TCM 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 Recursive-TCM. In contrast, the distributions predicted by Sph and Recursive-TCM 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 iso-magnification lines, as computed in Recursive-TCM 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 iso-magnification 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 .

Figure 11: Ratio of the measured excess convergence of substructures and a simple expectation , as a function of the subhalo mass, , in maps constructed using the CiC, Sph and Recursive-TCM. Solid lines show the median values in seven logarithmic bins in .

In Fig. 11 we compare the signal of identified substructures among the convergence maps. The x-axis shows the true subhalo mass . The y-axis 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 half-mass 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 Recursive-TCM), 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 Recursive-TCM 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 Recursive-TCM. 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 Recursive-TCM algorithm, will alleviate these discrepancies.

Figure 12: The cumulative number of local peaks in convergence maps detected in our WDM (top) and CDM (bottom) cluster-sized halo. In each case we show the results for a map created using three density projection methods: CiC (red lines), Sph (green line) and Recursive-TCM blue line. In addition, we show the abundance of substructures detected in 3D by the SUBFIND algorithm.

In order to quantify the performance of Recursive-TCM 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 Recursive-TCM 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 (Recursive-TCM) peaks in the CDM map, and (CiC), (Sph) and (Recursive-TCM), 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 false-positives 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 Recursive-TCM 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 mass-resolution, or one that is not a self-bound 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 self-bound 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.

Figure 13: Zoom into a region, centred at relative to the main halo, showing the convergence field associated to a tidal feature. The left panel shows a map computed using a Recursive-TCM projection method, whereas the right panels shows one computed using the Sph method.

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, high-precision theoretical predictions are essential to ensure the correct interpretation of future datasets.

We presented Recursive-TCM, 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 Recursive-TCM to cluster sized-haloes 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 shot-noise 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 Recursive-TCM, 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 Recursive-TCM, which, however, can be mostly eliminated by simple correction procedures. Therefore, this method is well suited for creating high-precision predictions for the relation between the underlying cosmological model and the expected signatures of small-cale structure in strong gravitational lensing observations.

With Recursive-TCM, 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 self-bound 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 state-of-the-art 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 large-scale 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 flux-ratio anomalies in multiple-image 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.


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 AST-0808398 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.


  1. pagerange: Noiseless Gravitational Lensing SimulationsNoiseless Gravitational Lensing Simulations
  2. pubyear: 2013
  10. The term Recursive-TCM abbreviates for ”Recursive deposit of Tethrahedra approximated by their Center of Mass”.
  11. 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).
  12. For simplicity, we refrain from also applying a force-range cut-off.
  13. Constructing the tessellation can be a computationally expensive task for state-of-the-art 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 glass-like distribution): each set of 8 grid points defines a cube that is subdivided into six disjoint tetrahedra.
  14. Our force resolution is also much lower.


  1. Abel T., Hahn O., Kaehler R., 2012, MNRAS, 427, 61
  2. Angulo R. E., Baugh C. M., Lacey C. G., 2008, MNRAS, 387, 921
  3. Angulo R. E., Hahn O., Abel T., 2013, arXiv:1304.2406
  4. Angulo R. E., Springel V., White S. D. M., Jenkins A., Baugh C. M., Frenk C. S., 2012, MNRAS, 426, 2046
  5. Aubert D., Amara A., Metcalf R. B., 2007, MNRAS, 376, 113
  6. Bartelmann M., 2010, Classical and Quantum Gravity, 27, 233001
  7. Bartelmann M., Huss A., Colberg J. M., Jenkins A., Pearce F. R., 1998, A&A, 330, 1
  8. Bartelmann M., Schneider P., 2001, Phys. Rep., 340, 291
  9. Bradač M., Schneider P., Lombardi M., Steinmetz M., Koopmans L. V. E., Navarro J. F., 2004, A&A, 423, 797
  10. Couchman H. M. P., Barber A. J., Thomas P. A., 1999, MNRAS, 308, 180
  11. Dalal N., Kochanek C. S., 2002, ApJ, 572, 25
  12. Davis M., Efstathiou G., Frenk C. S., White S. D. M., 1985, ApJ, 292, 371
  13. Diemand J., Kuhlen M., Madau P., 2007, ApJ, 667, 859
  14. Efstathiou G., Davis M., White S. D. M., Frenk C. S., 1985, ApJS, 57, 241
  15. Efstathiou G., Eastwood J. W., 1981, MNRAS, 194, 503
  16. Fontanot F., Puchwein E., Springel V., Bianchi D., 2013, arXiv:1307.5065
  17. Fontanot F., Springel V., Angulo R. E., Henriques B., 2012, MNRAS, 426, 2335
  18. Gao L., Navarro J. F., Frenk C. S., Jenkins A., Springel V., White S. D. M., 2012, MNRAS, 425, 2169
  19. Giannantonio T., Porciani C., Carron J., Amara A., Pillepich A., 2012, MNRAS, 422, 2854
  20. Hahn O., Abel T., 2011, MNRAS, 415, 2101
  21. Hahn O., Abel T., Kaehler R., 2013, MNRAS
  22. Hilbert S., Hartlap J., White S. D. M., Schneider P., 2009, A&A, 499, 31
  23. Hilbert S., Marian L., Smith R. E., Desjacques V., 2012, MNRAS, 426, 2870
  24. Hockney R. W., Eastwood J. W., 1981, Computer Simulation Using Particles. Computer Simulation Using Particles, New York: McGraw-Hill, 1981
  25. Huterer D., 2010, General Relativity and Gravitation, 42, 2177
  26. Jain B., Seljak U., White S., 2000, ApJ, 530, 547
  27. Kaehler R., Hahn O., Abel T., 2012, arXiv:1208.3206
  28. Klypin A., Kravtsov A. V., Valenzuela O., Prada F., 1999, ApJ, 522, 82
  29. Kochanek C. S., Dalal N., 2004, ApJ, 610, 69
  30. Komatsu E. et al., 2011, ApJS, 192, 18
  31. Kuhlen M., Vogelsberger M., Angulo R., 2012, Physics of the Dark Universe, 1, 50
  32. Li G.-L., Mao S., Jing Y. P., Kang X., Bartelmann M., 2006, ApJ, 652, 43
  33. Mao S., Schneider P., 1998, MNRAS, 295, 587
  34. Marian L., Hilbert S., Smith R. E., Schneider P., Desjacques V., 2011, ApJ, 728, L13
  35. Meneghetti M., Argazzi R., Pace F., Moscardini L., Dolag K., Bartelmann M., Li G., Oguri M., 2007, A&A, 461, 25
  36. Meneghetti M., Yoshida N., Bartelmann M., Moscardini L., Springel V., Tormen G., White S. D. M., 2001, MNRAS, 325, 435
  37. Metcalf R. B., Madau P., 2001, ApJ, 563, 9
  38. Monaghan J. J., 1992, ARA&A, 30, 543
  39. Moore B., Ghigna S., Governato F., Lake G., Quinn T., Stadel J., Tozzi P., 1999, ApJ, 524, L19
  40. Natarajan P., De Lucia G., Springel V., 2007, MNRAS, 376, 180
  41. Natarajan P., Springel V., 2004, ApJ, 617, L13
  42. Navarro J. F., Frenk C. S., White S. D. M., 1996, ApJ, 462, 563
  43. Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ, 490, 493
  44. Oguri M., Takada M., 2011, Phys. Rev. D, 83, 023008
  45. 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
  46. Pandey B., White S. D. M., Springel V., Angulo R., 2013, arXiv:1301.3789
  47. Peebles P. J. E., 1971, A&A, 11, 377
  48. Peter A. H. G., Rocha M., Bullock J. S., Kaplinghat M., 2013, MNRAS, 430, 105
  49. Rau S., Vegetti S., White S. D. M., 2013, MNRAS, 430, 2232
  50. Sato M., Hamana T., Takahashi R., Takada M., Yoshida N., Matsubara T., Sugiyama N., 2009, ApJ, 701, 945
  51. Schaap W. E., van de Weygaert R., 2000, A&A, 363, L29
  52. Shandarin S., Habib S., Heitmann K., 2012, Phys. Rev. D, 85, 083005
  53. Springel V. et al., 2008, MNRAS, 391, 1685
  54. Springel V. et al., 2005, Nature, 435, 629
  55. Springel V., Yoshida N., White S. D. M., 2001, New Astronomy, 6, 79
  56. Takahashi R., Oguri M., Sato M., Hamana T., 2011, ApJ, 742, 15
  57. Vale C., White M., 2003, ApJ, 592, 699
  58. Vegetti S., Koopmans L. V. E., Bolton A., Treu T., Gavazzi R., 2010, MNRAS, 408, 1969
  59. Vegetti S., Lagattuta D. J., McKean J. P., Auger M. W., Fassnacht C. D., Koopmans L. V. E., 2012, Nature, 481, 341
  60. Viel M., Becker G. D., Bolton J. S., Haehnelt M. G., 2013, Phys. Rev. D, 88, 043502
  61. Wambsganss J., Cen R., Ostriker J. P., 1998, ApJ, 494, 29
  62. Wang H. Y., Mo H. J., Jing Y. P., 2007, MNRAS, 375, 633
  63. Xu D. D. et al., 2009, MNRAS, 398, 1235
This is a comment super asjknd jkasnjk adsnkj
The feedback cannot be empty
Comments 0
The feedback cannot be empty
Add comment

You’re adding your first comment!
How to quickly get a good reply:
  • Offer a constructive comment on the author work.
  • Add helpful links to code implementation or project page.