GLAMER Part II: Multiple Plane Gravitational Lensing

GLAMER Part II: Multiple Plane Gravitational Lensing

Margarita Petkova,R. Benton Metcalf, Carlo Giocoli

Dipartimento di Fisica e Astronomia, Università di Bologna, viale Berti Pichat 6/2, 40127, Bologna, Italy
Department of Physics, Ludwig-Maximilians-Universität, Scheinerstr. 1, D-81679 München, Germany
Excellence Cluster Universe, Boltzmannstr. 2, D-85748 Garching, Germany
INAF - Osservatorio Astronomico di Bologna, via Ranzani 1, 40127, Bologna, Italy
INFN - Sezione di Bologna, viale Berti Pichat 6/2, 40127, Bologna, Italy
September 30, 2019

We present an extension to multiple planes of the gravitational lensing code GLAMER. The method entails projecting the mass in the observed light-cone onto a discrete number of lens planes and inverse ray-shooting from the image to the source plane. The mass on each plane can be represented as halos, simulation particles, a projected mass map extracted form a numerical simulation or any combination of these. The image finding is done in a source oriented fashion, where only regions of interest are iteratively refined on an initially coarse image plane grid. The calculations are performed in parallel on shared memory machines. The code is able to handle different types of analytic halos (NFW, NSIE, power-law, etc.), haloes extracted from numerical simulations and clusters constructed from semi-analytic models (MOKA). Likewise, there are several different options for modeling the source(s) which can be distributed throughout the light-cone. The distribution of matter in the light-cone can be either taken from a pre-existing N-body numerical simulations, from halo catalogs, or are generated from an analytic mass function. We present several tests of the code and demonstrate some of its applications such as generating mock images of galaxy and galaxy cluster lenses.


1 Introduction

Gravitational lensing has become a tool of greater and greater importance to cosmology and the study of galaxy structure. Large scale weak lensing surveys are being used to measure cosmological parameters and study dark energy111see DES:, Pan-Starrs : and larger ones are planned in space and with purpose built telescopes (Laureijs et al., 2011; LSST Science Collaborations et al., 2009). At the same time strong lensing is being used to study cosmology (Suyu et al., 2012; Treu et al., 2013) and the mass distribution within galaxies down to unprecedentedly small scales (Metcalf & Amara, 2012; Xu et al., 2013; Vegetti et al., 2012). Long term monitoring of strong lenses is being carried out with networks of telescopes spanning the globe to obtain time delays (Sluse et al., 2012; Treu et al., 2013). The Hubble Space Telescope (HST) has spent many orbits observing lenses (Postman et al., 2012; Coe et al., 2012; Coe et al., 2013; Zitrin et al., 2013) and many more are scheduled.222HST Frontier Fields: Getting the most out of these huge observational efforts will require very good simulations of the lensing phenomena for interpretation of the data and to control systematic effects. These simulations should include all the physical and instrumental effects that we know will affect the data. Currently, weak lensing simulations have limitations in resolution and the contribution of baryonic physics is usually neglected. In addition, they often use unperturbed light rays instead of fully tracing their paths through the simulation. Individual galaxies and the distortions they experience are not resolved in these simulations. Strong lensing and cluster lens simulations often use highly idealized lens models and usually neglect the environment of the lens and contributions to the lens that lie at different redshifts along the line of sight. The GLAMER code is an attempt to improve this situation in a way that can be applied to many different lensing situations on different scales.

Typically, simulations of strong gravitational lensing treat the lens as a single plane without foreground and background structure. Systematic studies, on clusters extracted from a cosmological numerical simulation, have been conducted by Meneghetti et al. (2010a, 2011) With the aim of better understanding of the properties of galaxy clusters that can potentially act as strong lenses. In these works however the simulated clusters are analyzed using the single lens approximation. The same method has also been used by Giocoli et al. (2012a); Boldrin et al. (2012), where the individual cluster components are analytically modeled. However, it is probable that structures along the line of sight have a significant effect on the observed properties of the lens: both in the weak and in the strong lensing regime. For example, the inconsistencies in the flux ratios between four-images quasars (Metcalf & Madau, 2001; Dalal & Kochanek, 2002; Metcalf & Zhao, 2002; Amara et al., 2006; Macciò & Miranda, 2006; Chen, 2009; Xu et al., 2009; Chen et al., 2011) can be at least partially explained by line-of-sight objects (Metcalf, 2005; Puchwein & Hilbert, 2009; Wambsganss et al., 2005; Xu et al., 2010). Furthermore, in cluster lensing, the mass and concentration of halos derived from lensing properties can be affected by the line-of-sight contribution. Puchwein & Hilbert (2009) showed that the presence of uncorrelated large scale structures can boost the Einstein radius by as much as 50%.

Weak lensing simulations are usually done by creating a shear map on a fixed resolution grid via FFT of the potential on one uniform grid (Jain et al., 2000; Vale & White, 2003; Pace et al., 2007; Angulo et al., 2013; Hamana & Mellier, 2001; Sato et al., 2009; Takahashi et al., 2011) or on separate grids for long and short range (Hilbert et al., 2007, 2009). The resolution of the grid is much larger than the image of a typical galaxy. The shear at the positions of individual galaxies is found by interpolating off this grid and is assumed to be uniform over the image of the source. The lensing mass in the light -cone is projected onto multiple lensing planes and rays are traced from the observer to the source. Also the simulation used to represent the mass distribution typically does not have high enough resolution to represent the inner regions of halos or any of the baryonic physics happening there. Because of these limitations no strong or quasi-strong lensing of individual galaxies will occur in these simulations. Although the contribution from smaller scales to the weak lensing might be small it is a potential source of systematic error and it is important to characterize it in the future high precision weak lensing surveys.

To ensure better resolution and less computational demand, other codes have used the tree force algorithm (Barnes & Hut, 1986) to calculate the deflection and distortion of rays. Wambsganss et al. (1998) used this method particularly for simulations of quasar microlensing. And Aubert et al. (2007) use it for lensing by N-body particles. Both of these use the single plane approximation.

Another set of codes (Fluke et al., 1999; Killedar et al., 2012) use the ray bundle method, where a collection of rays is traced from the image to the observer, integrating the and computing the shear and convergence of the bundle along the way. A third class of codes (Couchman et al., 1999; Carbone et al., 2008) calculate the distortion to ray-bundles along unperturbed paths that pass through a three dimensional simulated potential of a cosmological simulation. Still others use the unperturbed paths, also know as the the Born approximation, but put the mass from the simulations on multiple planes (White & Vale, 2004; Hilbert et al., 2007) The Born approximation can cause systematic errors and is clearly not sufficient in the case of strong lensing.

Figure 1: A schematic view of the multi-plane formalism, as described in eqn. 19.

In Metcalf & Petkova (2013) (here after MP) we describe the GLAMER code’s structure, its adaptive mesh refinement capabilities and how it calculates the deflection of light rays by a single lens plane. These capabilities alleviate many of the limitations codes have with resolution and make it possible to construct fully resolved images of millions of galaxies in a single simulation. In this paper, we describe the code’s ability to simulate lensing in three dimensions; that is with the lensing mass and also the sources distributed throughout the light-cone in angle and redshift. We also discuss some of the ways a light-cone can be generated and input into the code.

There are many interesting questions that can be addressed with a lensing code like the one we have developed. Among them are: the importance of line-of-sight objects to strong lenses, the importance of strongly lensed individual galaxies to weak lensing surveys, the statistical properties of strong lenses, the importance of multi-plane effects on precision weak lensing, and many more.

We start this paper by presenting the basics of lensing in Sec. 2. We elaborate on the multiple plane formalism in Sec. 3. We continue in section. 4 by describing the implementation of this formalism in the GLAMER code. In Sec. 5 we describe the different types of light-cones implemented in the code. In Sec. 5.2 we discuss the use a simulated halo catalog to populate the light-cone and in Sec. 5.3, the generation of a light-cone from a mass function. Then in Sec. 6 we perform several tests of the code, determining the optimal number of lensing planes to be used. We also show some future applications of the code, such as Einstein radius mapping of galaxy clusters and producing galaxy-galaxy strong lens catalogs. We conclude our discussions in Sec. 8. For more extended information we have also included several appendices.

2 lensing basics

We will assume a weak field approximation to General Relativity. In this case, light propagates through the universe on the geodesics of the perturbed Robertson-Walker metric. The line element of this metric in the Newtonian gauge is


where is the cosmological scale factor, is the Newtonian potential, is the speed of light, is the conformal time time, and


Here is the comoving distance and , the comoving angular diameter distance, is defined for different curvatures as follows:


where is the curvature scale ( in the standard cosmological model before radiation domination).

The geodesic equation for a null geodesic on small enough scale to ignore the curvature is


neglecting higher order terms in and being an affine distance along the path. From the prospective of the light ray this can be written more simple as


where and are the coordinates perpendicular and parallel to the direction of propagation respectively.

If we imagine a thin slab that the ray enters perpendicularly we can integrate equation (6) with a unperturbed path within the slab. We can then calculate the deflection angle of the ray on leaving the slab,


where is now the surface potential.

The surface potential is related to the surface density by


where is the surface density within the slab which can be solved to get


The position on the sky will be . The position of this point where there no lensing will be , the source position. The magnification matrix, or the Jacobian matrix of the map between these coordinates, is


This matrix is often decomposed into the convergence, and the components of shear . In the notation of appendix A, this decomposition is


We have included a third component of “shear”, , which is not usually included, but, as we will see, is required when there are more than one lens plane. The magnification of an infinitesimal point on the sky is then


2.1 single lensing plane

It is instructive to consider the case of a single lens plane where it is assumed that the light rays follow unperturbed geodesics outside of the lens plane. This is the case usually considered in the literature so expressing it here in our notation will serve to clarify our approach and what new phenomena arise from multi-plane lensing.

In this case, the source and image positions are related by


where is the angular diameter distance to the source and , between lens plane and the source (). In this case, the convergence can be written in terms of the surface density on the plane


so in this case can be considered a dimensionless surface density. This is not the case for multiple lens planes. The deflection field is a potential field which mean that , the field has no curl and the shear field has no B-modes. Also every point on the lens plane has a one-to-one correspondence to a point on the image plane. As we shall see these things are not guaranteed when there are multiple lens planes.

Figure 2: Halo mass function, recreated in a light -cone of , for different redshift bins. The left panel shows the Press-Schechter mass function, the middle one – the Sheth-Tormen one, and the right one – a power law with a slope of . The distribution of masses (histograms), for which we put together 64 different light-cone realizations, is in very good agreement with the theoretical mass functions (lines) for all redshifts.

3 multi-plane formalism

In the multi-plane formalism, the mass in the universe is projected onto multiple lensing planes, perpendicular to the line of sight. As light travels from the lensed object to the observer, it experiences deflection from each plane and passes unperturbed through the space between the planes. The multiple plane lens equation has the following form


where is the number of planes, is the angular diameter distance to the source and – between plane and the source. See figure 1.

For the purpose of our work, from now on we will use only comoving positions and angular distances:


where is the redshift of the plane in consideration.

Figure 3: Position of the four images (different symbols) of a lensed quasar at redshift by a galaxy at redshift in simulations with different numbers of lensing planes (6, 8, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 32, 44, 56, 68, 80, 92, 104). The color of the symbols darkens with increasing number of planes. The positions clearly converge as the number of planes increases.
Figure 4: Image positions and magnification as a function of the comoving distance between the lensing planes (different number of planes) for the lensing simulation of a quasar at redshift by a galaxy at redshift . Separation of approximately 300 comoving Mpc  corresponds to 22 lensing planes in this case.

It is important to note here that since we calculate the angular diameter distances between planes assuming the Robertson-Walker metric, the mass added to those lenses planes will cause our light-cone to effectively contain more mass than the average for the Friedman-Lemaitre universe used to calculate the distances if it is not compensated for – the average angular size through the planes will the smaller on average. This is not a small effect if a significant fraction of the mass in the universe is represented on the planes. When needed we correct for this by add a uniform negative mass density to each of the planes that compensates for the clustered mass so that the average surface density is zero on each plane (Schneider & Weiss, 1988)333This could be avoided by using a Dryer-Roeder angular diameter distances.. This will be further discussed in section 5.

The deflection angles at each plane add to give the angle with respect to the radial direction at the th plane. From equation (15) it follows then that the position on the th plane is


Applying this to the th plane in terms of the th plane and subtracting the equations gives a recursion relation for the positions on the planes


with the initial conditions


The first of these is the requirement that all the rays converge at the observer.

The magnification matrix defined in (10) can be expressed on each plane as


This is what the normal magnification matrix would be if the th plane where the source plane. Differentiating equation (19) yields


with the initial conditions


The forcing term is


where the notation of Appendix A is used. The coefficients are related to the Newtonian surface potential on the planes by


The term (23) can be written out as


It can be seen in equation (3) that there is a component to the magnification matrix which does not exist in the case of one lens. In appendix A, it is shown that this implies a rotation in the image (Pen & Mao, 2006). The case of two lens planes is worked out explicitly for pedagogical purposes in appendix B.

4 code implementation

Multi-plane lensing is implemented as an extension to the preexisting GLAMER code. Other aspects of the GLAMER code are described in more detail in Metcalf & Petkova (2013) (MP). The code is written in C++ in an object oriented manor so that characteristics of the lenses and sources can be chosen by the user in a very flexible fashion. There are a number of options allowed which will be described briefly here.

Rays are shot from the observer to the source plane using equations (19) and (3). In each plane the deflection, convergence and shear are calculated using the methods described in MP. The mass distribution on each plane can be represented in several different ways. A surface density map can be input in FITS format. This option is useful for representing the output of N-body simulations and semi-analytic methods for constructing galaxies and galaxy clusters such as MOKA (Giocoli et al., 2012). The lens planes can also contain analytic “halos”. The halos can have a variety of different mass profiles – Navarro-Frenk-White (NFW) (Navarro et al., 1997), non-singular isothermal sphere (NSIE), power-law, Hernquist (1990), point-mass, etc. These “halos” can represent the baryonic component of galaxies, the dark matter halos or even stars. The contributions of the halos to the lensing quantities are calculated using a modified tree algorithm when there are enough of them to warrant it (see MP). Finally, the mass on a plane can be represented by simulation particles with an adaptive smoothing in which case the lensing quantities are calculated by tree algorithm. The code allows for any combination of these representations. For example, the dark matter halos could be represented by NFW profiles, the baryonic galaxy by an NSIE, the mass outside of halos by a mass map and the stars by point masses all at the same time.

Once the lens has been constructed rays are shot back to the source plane in a uniform grid on the image plane. This initial gridding can be used to make shear or magnification map with uniform resolution if that is desired. As described in MP there are a number of options for refining the grid of rays. The code can be made to find critical curves and caustics and increase the resolution of them to the desired level. It can also be made to refine the grid to resolve a particular source. The ray shooting is parallelized with POSIX threads, which increases the performance significantly especially when the number of planes is large.

There are also a variety of options for representing the sources. The simplest source is circular with uniform surface brightness. The source can also have an analytic surface brightness distribution or a pixelized surface brightness distribution that can be input in FITS format. There can be any number of sources and they can be of mixed type. Each source can also have a different redshift. Because of the adaptive gridding of the rays each source can be resolved by shooting a relatively small number of rays to the particular redshift of the source.

Figure 5: Relative error (31) of the lensing properties – inverse of the magnification (top), magnitude of the deflection angle, magnitude of the shear and convergence (bottom) – for simulations with the same light-cone, but differing comoving distance between the lensing planes (different number of planes). Separation of approximately 300 comoving Mpc  corresponds to 18 lensing planes. The field of view is , the mass function goes down to halos and up to redshift . The true solution in this case is assumed to be the one with one plane for each halo at precisely its redshift.
Figure 6: Distribution of the deflection angle, convergence, and shear of a set of 32 random and uncorrelated fields of view. The area is and the NFW halos are from a ST mass function down to , distributed on 20 planes. The source is at redshift . The results agree very well with simulations by Pace et al. (2011)

5 Generating Light Cones

In this section, we describe the different types of light-cones that are implemented within or can be input into the GLAMER   code. The mass within the light-cone can be generated in four different ways: () particles can be taken directly from a N-body simulation of a light-cone and sorted onto planes, () halos can be identified in the N-body simulation and inputted into GLAMER as a halo catalog, () random halos can be generated within GLAMER from a redshift dependent mass function and other analytic descriptions of their internal parameters, () two dimensional maps of surface density, which are generated from a N-body simulation or semi-analytic simulation, can be input into GLAMER. A light-cone can have lensing planes of any mixture of these types – for example a galaxy cluster might be represented by a mass map and the foreground and background objects generated randomly from a mass function. We will discuss each of these cases below.

5.1 Light-Cones from Simulation Particles

GLAMER allows the user to input a list of particle positions and smoothing scales for a light -cone. The smoothing is best done adaptively according to particle density in three dimensions – as in Smooth Particle Hydrodynamics (SPH). The particles are sorted onto the desired number of lens planes and then a tree structure is constructed on each plane to calculate the deflection and other lensing quantities (see MP for details). For state-of-the-art numerical simulations the number of particles within a complete light-cone can be very large which makes memory management difficult. For this reason, we favor using methods () or () or a combination of them to represent the output of a N-body simulation, but the ability to use the particles directly is maintained. A simulation of a single object with relatively few particles can be implanted within a light-cone that is represented in a different way.

5.2 Light-Cones from Halo Catalog

GLAMER   is able to read in a catalog of halos and their properties. These halos can be taken from an N-body simulation or any other source. The halos are projected onto the lensing planes keeping their angular positions fixed and a tree structure is constructed on each plane. The force and deflections are calculated by tree algorithm (see MP).

To represent an N-body halo as an analytic halo, parameters of the halo model must be matched to quantities calculated from the particle distribution – virial mass, half mass radius, maximum circular velocity, moments of inertia, etc. This procedure is somewhat dependent on the halo catalog used and what information it provides. One can also attempt to represent the baryons that might not be in the simulation with an additional “halo” within the dark matter halo. This requires further modeling of the relation between halo mass and stellar mass and the distribution of baryonic mass within a N-body halo of a given type. These are interesting subjects that are being actively investigated. Lensing has great potential for investigating these relationships. We will not discuss them in detail here since our purpose is to introduce a tool for investigating them and not to advocate any particular model.

We have implemented one full pipeline for N-body to sky images using the Millennium Simulation catalog of halos and with sources from the Millennium Run Observatory (Overzier et al., 2013). This project will be discussed in a separate paper.

5.3 Light-Cones from Generated Halos

The halos catalog can also be created internally and then sorted onto lens planes as in the case a catalog from a N-body simulation. No attempt has been made yet to cluster these halos or to create subhalos within halos for the full light-cone. The halos masses are drawn randomly from a mass function. The mass functions implemented in our code so far are:

  • Press & Schechter (1974) (PS);

  • Sheth & Tormen (1999) (ST);

  • Power-Law (PL), given by the following equation:


    where represents the slope of the power-law, , and the non-linear mass at redshift given such that , with the overdensity required for spherical collapse at , and is the variance in the linear fluctuation field when smoothed with a top-hat filter of scale where is the comoving density of the background. The Sheth & Tormen (1999) mass function in CD on cluster mass scales is well approximated by .

The three different mass function are plotted in Fig. 2 . For a test, we have populated light-cones of and measured the mass function in different redshift bins. The histograms of the halo masses, for which we put together different realizations of the light-cone, are also plotted in Fig. 2. The agreement is very good.

In order to model the halo lensing properties we need to make some assumption about the dark matter distribution inside the virial radius of each system - the virial radius defines the scale inside which the overdensity exceeds the critical value predicted for collapse (Gunn, 1977; Bond et al., 1991; Lacey & Cole, 1993). Numerical simulations predict that the dark matter follows an NFW profile (Navarro et al., 1997, 2004) or something very close to it. Using this model we must also fix the concentration (the ratio between the scale and the virial radius: ) for each halo. In GLAMER   we have implemented four model for the mass-concentration relation (we refer the reader to the cited papers for more details): Zhao et al. (2009), Muñoz-Cuartas et al. (2011), Giocoli et al. (2012b) and a power-law relation:


where the redshift evolution factor is motivated by the change of the normalization as seen in cosmological numerical simulations.

This method of populating the light-cone with halos has several limitations with respect to using a halos catalog derived from an simulation. Since the positioning of the halos is random, we discard any clustering or particular environmental effects. Halos do not have subhalos within them. The advantage is that there is no resolution limit and no limit to how many light-cones can be created.

As in the halos from N-body case, the baryonic component can be modeled by adding an additional halo to the center of each NFW halo. So far we have implemented one simple way of doing this. The galaxy is modeled as a NSIE with a random orientation. The mass of the NSIE is calculated using the stellar mass/halo mass relation of Moster et al. (2010):


with , , , and . This is derived from matching the abundance of observed galaxies to the number of predicted halos in mass bins. The NFW halos mass is reduced accordingly. Much more sophisticated ways of doing this are possible and will be implemented in the future.

The centers of some halos will lie outside the light-cone, but they will intersect the cone’s boundary. This can result in boundary effects because the density out the boundary of the field is artificially low. A buffering option is implemented to avoid this. All the halos are created in the region outside of the cone but within a fixed proper distance () from the boundary. The result is a funnel shape instead of a cone.

5.4 Mass Maps

Maps of surface density in FITS format can be read into GLAMER and placed on a lens plane at the desired redshift. The deflection and shear are found by Fourier transforming the surface density, solving equation (9) for the potential, evaluating equations (26) and then Fourier transforming back to real space. Maps of shear and deflection are stored so that this need be done only once. The lensing quantities are then evaluated by interpolating from these maps during the ray shooting procedure.

When using a light-cone taken from a large-scale simulation it is very advantageous to project the particles in shells around the observer onto planes and then construct two dimensional density maps by CIC (clouds in cell) or some other method. Because the shells can be quite thick, as we will see in section 6, this can save an enormous amount of memory and computational time over using the particles directly without significant lose in accuracy. We have also used this capability to implant simulated galaxy clusters into light-cones that are otherwise populated with halos. Fig. 7 shows and example of this.

5.5 Sources

Sources can be distributed throughout the light-cone. They can be represented by a number of analytic forms including a circular uniform surface brightness profile, a Sérsic (1963) profile and an inclined exponential disk profile. A pixelized image of the source can also be used. A catalog of sources can be read in from an external source or individual sources can be created randomly. GLAMER can adaptively find and refine around the image even when they are highly distorted and there are multiple images of the source. The rays are shot only back to the redshift of this source.

A numerical effect can occur if the source is placed in the center of a halo and then the halo is projected to the nearest lens plane that is at lower redshift than the source. This will produce alignments with significant differences in radial distance and spurious lensing. The distance between planes should be small enough to make the lensing by the nearest plane small, but the perfect alignment can have unwanted effects. To avoid this, the deflection is turned off for the lens plane that contains the source’s halo.

Figure 7: Convergence map of a cluster without (left) and with (right) structure along the line of sight. The cluster lens is located at redshift , it posses a mass of and a concentration of . The green contours in the two figures show the location of radial and tangential critical lines. The cluster on the left possess an Einstein radius of arcsec while in the right panel, where the structure along the line of sight are included, the Einstein radius is arcsec.
Figure 8: The number of planes needed in a simulation versus source redshift. The values are chosen for separation distance between lensing planes of 300 comoving Mpc.

6 Tests

We perform several tests to check the accuracy and convergence of our multi-plane gravitational lensing implementation in the code GLAMER. In our first test we determine the convergence of our code by comparing the image positions and magnification of a four-imaged quasar, as a function of the number of lensing planes. In Fig. 3 we plot the position of the four images of the quasar at redshift , lensed by a galaxy at redshift . The field of view of with a buffer of 1 Mpc, filled with NFW halos from a ST mass function down to . We compute the lensing properties and change the number of planes in two plane interval as the main galaxy lens stays fixed on a plane at redshift .

The changes in the image positions are no greater than 0.4 arcsec  in one x-direction and 1.2 arcsec  in the y-direction. This can also be observed in Fig. 4, where plot the image positions and inverse magnifications as a function of the distance between the planes. The positions show very small ( 0.2 arcsec) scatter for 300 comoving Mpc  or less, equivalent to 22 lensing planes. The same holds for the inverse magnification.

Our next test aims to examine how other lensing properties, such as shear and convergence, converge with increasing the number of lensing planes. For the study we generate a light-cone with a field of view of that extends to redshift . We populate it with NFW halos from a ST mass function down to . As a “true” result we consider the case where an infinite number of planes is used, i.e. each lensing halo occupies a single lensing plane. We will denote this results as the “true” solution.

We compare the relative error of the lensing properties of simulations with different number of planes to the true solution. To do this we compare the image maps pixel by pixel and we formulate the the relative error in the following way:




and is the value of the quantity at the respective pixel .

In Fig. 5 we plot the relative error for the convergence, shear, deflection angle, and inverse magnification as a function of the comoving separation between the lensing planes. The relative error for the deflection angle and the inverse magnification is very small – for less than 300 comoving Mpc  between the planes it is around 0.1%. The error for the convergence and shear is ten times larger and drops down to less than 5% for less than 300 comoving Mpc, which supports the conclusions from our previous test. This separation corresponds to 18 lensing planes for a source at redshift .

In the next test, we investigate the agreement of our distributions of the lensing quantities with already published results. We compute histograms of the lensing properties of a field of view with the same parameters as the one from the previous test (field of view of , extending to redshift , populated with NFW halos from a ST mass function down to ), but keeping the number of planes fixed to 20. We change randomly the distribution of the NFW halos in the light-cone and obtain a set of 32 different and uncorrelated fields of view. In Fig. 6 we show the average of the convergence, shear, and deflection angle for the 32 realizations of the image plane and we overplot the standard deviation. The results agree very well with the ones from Pace et al. (2011), despite the much smaller field of view.

Finally, we show how structure along the line of sight may change the Einstein radius of a galaxy cluster. In our simulation we define the Einstein radius as the median distance of the tangential central critical points from the cluster center Meneghetti et al. (2010a). In the left panel of Fig. 7 we show convergence map of a halo created with MOKA (Giocoli et al., 2012a) at redshift . The cluster possess a virial mass of and a concentration of . In the right panel we include in the map the effect of uncorrelated structures along the line of sight, assuming an empty light-cone. The green lines in the maps represent the radial and tangential critical curves, where the magnification of the background sources, located at , is infinite. The green connected points in the center of the maps indicate the Einstein radius of the cluster, measured as the median distance of the tangential central critical points to the halo center. From the analysis we notice that while the Einstein radius in the left figure is arcsec, in the right increases to arcsec, where LSS are included. We will discuss motivate and better quantify this effect in a future paper Giocoli et al. in preparation.

We have also performed further small test to ensure the accuracy of our multi-plane lensing scheme. For example, we have tested that the average lensing mass on each plane is zero and that an empty light-cone with an arbitrary number of lensing planes always produces zero deflection. We also checked that an analytic lens placed in an empty light-cone with an arbitrary number of planes agrees fully with the single lens plane calculation.

From our tests we conclude that a separation of 300 comoving Mpc  or less is sufficient for accurate representation of the lensing field. In Fig. 8 we show the correspondence between number of planes and source redshift, that gives approximately 300 comoving Mpc  plane separation. For redshift , 18 planes are sufficient. For much higher redshifts, , the number increases to 32.

7 Some applications

In this section, we present some example applications to demonstrate GLAMER’s applicability to interface with other lensing codes and to produce simulated observations. We do this without fully describing the parts of the simulations that are not directly related to the functionality of GLAMER.

As previously discussed, Fig. 7 shows the convergence map for a simulated galaxy cluster. The cluster lens located at redshift was created using MOKA (Giocoli et al., 2012), and has as input parameters a mass of and a concentration of , and sources are located at redshift . The axial ratios of the cluster ellipsoid are randomly drown from the Jing & Suto (2002) distributions, they are and where , and represent the smallest, intermediate and largest axis of the halo ellipsoid, respectively; the random orientation of the systems, with respect to the line of sight, gives in the plane of the sky an axial ratio of the ellipse that is . As a result the cluster on the left is characterized by an Einstein radius arcsec. In GLAMER we have created an interface that allows MOKA outputs to be reprocessed including background and foreground objects by constructing a random light-cone as described in section 5.3. The MOKA map was in this way inserted into the light-cone. As in all these examples, the ray tracing was done without the Born approximation. Similar map of the shear and deflection are created. The critical lines and caustics can be mapped and the images of lensed background objects can be found. The created light-cone realization modify the size of the Einstein radius to arcsec.

Fig. 9 shows a map of the inverse of the magnification for a random 8.4 by 8.4 arcminute field. The light-cone was constructed from a halo catalog extracted from the Millennium simulation. Each halo is represented by a NFW profile for the dark matter and a NSIE profile for the inner regions. The grid in which the rays have been shot was refined to give higher angular resolution where the magnification is high.

Fig. 10 shows twelve examples of simulated galaxy-galaxy strong lenses. The ray shooting grid is refined to adapt to the surface brightness of the images. Each source is ray traced back to its own redshift and has a NFW+NSIE halo. Here the caustics in a random light-cone were found and an extra source was put close to the caustic. Many thousands of such lenses can be created easily. Such images are being used to test lens finding robots, lens modeling algorithms and predict the statistics of detectable lenses.

Figure 9: Map of the inverse of the magnification, , for a 8.4 by 8.4 arcmin random field with the sources plane at z=2.5. The image is on a regular grid, but the rays have been refined to resolve the critical curves so the high magnification ( and the ) regions had higher resolution before this image was interpolated to a regular grid. In some of the regions where the resolution of the adaptive grid is lower then the final image resolution.
Figure 10: Twelve simulated galaxy-galaxy lenses. The mass distributions of the galaxies and their dark matter halos are taken to be NSIE within NFW profiles. The visible galaxies are modeled as Sérsic (1963) profiles at the centers of the halos. Noise has been added and a Gaussian point spread function (psf) applied. Two of these are the same lens with different sources behind it. The output images are in standard fits files.

8 discussion and conclusions

We have presented an extension of the ray-shooting lensing code GLAMER to multiple planes and its ability to perform lensing simulations in a light-cone, rather than using a single lens. This allows us to perform more accurate statistical strong and weak lensing simulations and study properties of the strong lensing galaxies and galaxy clusters, as well as the structure along the line of sight outside of the primary lens.

We have presented in detail the methodology and numerical implementation that we have use. The mass in the universe, in the form of halos, is projected onto a discrete number of lensing planes and the lensing properties, such as deflection and convergence, are computed on each of these planes and summed up according to the lens equations. The summation of the contributions from the different halos is performed via a spacial tree algorithm in order to improve accuracy and speed. Details of the performance and limitations of the tree can be found in MP.

The halos, populating a light-cone, can be either taken from a dark matter simulation or randomly generated from a mass function at different redshifts. For the second one we have adopted three different mass functions – Press-Schechter, Sheth- Tormen and a simple power law. In this context more general and useful relations can be implemented and inputed as desired.

We have presented various tests and applications of our code. We have shown that the optimal number of planes per light-cone is achieved for a separation of circa 300 comoving Mpc. We have presented that the image position of a lensed quasar by a galaxy converges for this value. Moreover, the relative uncertainty between the lensing properties in a simulation with that given plane separation and an infinitesimal one is less than 3%. We have also shown that the statistical distribution for the lensing properties in our simulations agree very well with previously published results. Finally, we have demonstrated some of the applications that the code is already being applied to.


This research is part of the project GLENCO, funded under the Seventh Framework Programme, Ideas, Grant Agreement n. 259349.


  • Amara et al. (2006) Amara A., Metcalf R. B., Cox T. J., Ostriker J. P., 2006, MNRAS, 367, 1367
  • Angulo et al. (2013) Angulo R. E., Chen R., Hilbert S., Abel T., 2013, ArXiv e-prints
  • Aubert et al. (2007) Aubert D., Amara A., Metcalf R. B., 2007, MNRAS, 376, 113
  • Barnes & Hut (1986) Barnes J., Hut P., 1986, Nature, 324, 446
  • Bartelmann & Schneider (2001) Bartelmann M., Schneider P., 2001, Physics Reports, 340, 291
  • Boldrin et al. (2012) Boldrin M., Giocoli C., Meneghetti M., Moscardini L., 2012, MNRAS, 427, 3134
  • Bond et al. (1991) Bond J. R., Cole S., Efstathiou G., Kaiser N., 1991, ApJ, 379, 440
  • Carbone et al. (2008) Carbone C., Springel V., Baccigalupi C., Bartelmann M., Matarrese S., 2008, MNRAS, 388, 1618
  • Chen (2009) Chen J., 2009, AAP, 498, 49
  • Chen et al. (2011) Chen J., Koushiappas S. M., Zentner A. R., 2011, ApJ, 741, 117
  • Coe et al. (2012) Coe D., Umetsu K., Zitrin A., Donahue M., Medezinski E., Postman M., Carrasco M., Anguita T., et al. 2012, ApJ, 757, 22
  • Coe et al. (2013) Coe D., Zitrin A., Carrasco M., Shu X., Zheng W., Postman M., Bradley L., Koekemoer et al. 2013, ApJ, 762, 32
  • Couchman et al. (1999) Couchman H. M. P., Barber A. J., Thomas P. A., 1999, MNRAS, 308, 180
  • Dalal & Kochanek (2002) Dalal N., Kochanek C. S., 2002, ApJ, 572, 25
  • Fluke et al. (1999) Fluke C. J., Webster R. L., Mortlock D. J., 1999, MNRAS, 306, 567
  • Giocoli et al. (2012) Giocoli C., Meneghetti M., Bartelmann M., Moscardini L., Boldrin M., 2012, MNRAS, 421, 3343
  • Giocoli et al. (2012a) Giocoli C., Meneghetti M., Bartelmann M., Moscardini L., Boldrin M., 2012a, MNRAS, 421, 3343
  • Giocoli et al. (2012b) Giocoli C., Tormen G., Sheth R. K., 2012b, MNRAS, 422, 185
  • Gunn (1977) Gunn J. E., 1977, ApJ, 218, 592
  • Hamana & Mellier (2001) Hamana T., Mellier Y., 2001, MNRAS, 327, 169
  • Hernquist (1990) Hernquist L., 1990, ApJ, 356, 359
  • Hilbert et al. (2009) Hilbert S., Hartlap J., White S. D. M., Schneider P., 2009, AAP, 499, 31
  • Hilbert et al. (2007) Hilbert S., Metcalf R. B., White S. D. M., 2007, MNRAS, 382, 1494
  • Hilbert et al. (2007) Hilbert S., White S. D. M., Hartlap J., Schneider P., 2007, MNRAS, 382, 121
  • Jain et al. (2000) Jain B., Seljak U., White S., 2000, ApJ, 530, 547
  • Jing & Suto (2002) Jing Y. P., Suto Y., 2002, ApJ, 574, 538
  • Kaiser & Squires (1993) Kaiser N., Squires G., 1993, ApJ, 404, 441
  • Killedar et al. (2012) Killedar M., Lasky P. D., Lewis G. F., Fluke C. J., 2012, MNRAS, 420, 155
  • Lacey & Cole (1993) Lacey C., Cole S., 1993, MNRAS, 262, 627
  • Laureijs et al. (2011) Laureijs R., Amiaux J., Arduini S., Auguères J. ., Brinchmann J., Cole R., Cropper M., Dabin C., Duvet L., et al. 2011, ArXiv e-prints
  • LSST Science Collaborations et al. (2009) LSST Science Collaborations Abell P. A., Allison J., Anderson S. F., Andrew J. R., Angel J. R. P., Armus L., Arnett D., Asztalos S. J., Axelrod T. S., et al. 2009, ArXiv e-prints
  • Macciò & Miranda (2006) Macciò A. V., Miranda M., 2006, MNRAS, 368, 599
  • Meneghetti et al. (2010a) Meneghetti M., Fedeli C., Pace F., Gottlöber S., Yepes G., 2010a, AAP, 519, A90+
  • Meneghetti et al. (2011) Meneghetti M., Fedeli C., Zitrin A., Bartelmann M., Broadhurst T., Gottlöber S., Moscardini L., Yepes G., 2011, AAP, 530, A17+
  • Metcalf & Petkova (2013) Metcalf R., Petkova M., 2013, arXiv:1312.1128
  • Metcalf (2005) Metcalf R. B., 2005, ApJ, 629, 673
  • Metcalf & Amara (2012) Metcalf R. B., Amara A., 2012, MNRAS, 419, 3414
  • Metcalf & Madau (2001) Metcalf R. B., Madau P., 2001, MNRAS, 563, 9
  • Metcalf & Zhao (2002) Metcalf R. B., Zhao H., 2002, ApJL, 567, L5
  • Miralda-Escude (1991) Miralda-Escude J., 1991, ApJ, 380, 1
  • Moster et al. (2010) Moster B. P., Somerville R. S., Maulbetsch C., van den Bosch F. C., Macciò A. V., Naab T., Oser L., 2010, ApJ, 710, 903
  • Muñoz-Cuartas et al. (2011) Muñoz-Cuartas J. C., Macciò A. V., Gottlöber S., Dutton A. A., 2011, MNRAS, 411, 584
  • Navarro et al. (1997) Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ, 490, 493
  • Navarro et al. (2004) Navarro J. F., Hayashi E., Power C., Jenkins A. R., Frenk C. S., White S. D. M., Springel V., Stadel J., Quinn T. R., 2004, MNRAS, 349, 1039
  • Overzier et al. (2013) 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
  • Pace et al. (2007) Pace F., Maturi M., Meneghetti M., Bartelmann M., Moscardini L., Dolag K., 2007, AAP, 471, 731
  • Pace et al. (2011) Pace F., Moscardini L., Bartelmann M., Branchini E., Dolag K., Grossi M., Matarrese S., 2011, Monthly Notices of the Royal Astronomical Society, 411, 595
  • Pen & Mao (2006) Pen U.-L., Mao S., 2006, MNRAS, 367, 1543
  • Postman et al. (2012) Postman M., Coe D., Benítez N., Bradley L., Broadhurst T., Donahue M., Ford H., Graur O., Graves G., Jouvel S., et al. 2012, ApJ Sup., 199, 25
  • Press & Schechter (1974) Press W. H., Schechter P., 1974, ApJ, 187, 425
  • Puchwein & Hilbert (2009) Puchwein E., Hilbert S., 2009, MNRAS, 398, 1298
  • Sato et al. (2009) Sato M., Hamana T., Takahashi R., Takada M., Yoshida N., Matsubara T., Sugiyama N., 2009, ApJ, 701, 945
  • Schneider & Weiss (1988) Schneider P., Weiss A., 1988, ApJ, 327, 526
  • Seitz & Schneider (1997) Seitz C., Schneider P., 1997, AAP, 318, 687
  • Sérsic (1963) Sérsic J. L., 1963, Boletin de la Asociacion Argentina de Astronomia La Plata Argentina, 6, 41
  • Sheth & Tormen (1999) Sheth R. K., Tormen G., 1999, MNRAS, 308, 119
  • Sluse et al. (2012) Sluse D., Chantry V., Magain P., Courbin F., Meylan G., 2012, AAP, 538, A99
  • Suyu et al. (2012) Suyu S. H., Auger M. W., Hilbert S., Marshall P. J., Tewes M., Treu T., Fassnacht C. D., Koopmans L. V. E., Sluse D., Blandford R. D., Courbin F., Meylan G., 2012, ArXiv e-prints
  • Takahashi et al. (2011) Takahashi R., Oguri M., Sato M., Hamana T., 2011, ApJ, 742, 15
  • Treu et al. (2013) Treu T., Marshall P. J., Cyr-Racine F.-Y., Fassnacht C. D., Keeton C. R., Linder E. V., Moustakas L. A., Bradac M., Buckley-Geer E., et al. 2013, ArXiv e-prints
  • Treu et al. (2013) Treu T., Marshall P. J., Cyr-Racine F.-Y., Fassnacht C. D., Keeton C. R., Linder E. V., Moustakas L. A., et al. 2013, ArXiv e-prints
  • Vale & White (2003) Vale C., White M., 2003, ApJ, 592, 699
  • Vegetti et al. (2012) Vegetti S., Lagattuta D. J., McKean J. P., Auger M. W., Fassnacht C. D., Koopmans L. V. E., 2012, Nature, 481, 341
  • Wambsganss et al. (2005) Wambsganss J., Bode P., Ostriker J. P., 2005, ApJL, 635, L1
  • Wambsganss et al. (1998) Wambsganss J., Cen R., Ostriker J. P., 1998, ApJ, 494, 29
  • White & Vale (2004) White M., Vale C., 2004, Astroparticle Physics, 22, 19
  • Xu et al. (2010) Xu D. D., Mao S., Cooper A. P., Wang J., Gao L., Frenk C. S., Springel V., 2010, MNRAS, 408, 1721
  • Xu et al. (2009) Xu D. D., Mao S., Wang J., Springel V., Gao L., White S. D. M., Frenk C. S., Jenkins A., Li G., Navarro J. F., 2009, MNRAS, 398, 1235
  • Xu et al. (2013) Xu D. D., Sluse D., Gao L., Wang J., Frenk C., Mao S., Schneider P., 2013, ArXiv e-prints
  • Zhao et al. (2009) Zhao D. H., Jing Y. P., Mo H. J., Bnörner G., 2009, ApJ, 707, 354
  • Zitrin et al. (2013) Zitrin A., Meneghetti M., Umetsu K., Broadhurst T., Bartelmann M., Bouwens R., Bradley L., Carrasco M. e. a., 2013, ApJL, 762, L30

Appendix A ellipticity in tensor representation

We have found the following formalism useful in dealing with propagating magnification matrices and dealing with galaxy ellipticities.

The magnification matrix, ellipticity or any 2-by-2 matrix can be represented as


where , and are scalar coefficients which may depend on position on the sky. The matrices are


These matrices obey the following algebraic relations


Some of the standard quantities and transformations are easily expressed in this basis:

trace (38)
determinant (39)
transpose (40)

The rotation operator in this space is


Using the algebraic relations, the behavior of these matrices under rotations can be written


Thus and are rotationally invariant.

The ellipticity of a galaxy can be represented by a matrix with only and components


Using equation (43) it can be seen that an ellipticity with only an component will become an ellipticity with a component when rotated. In this way an angle can be associated with any ellipticity – the rotation angle required to rotate the ellipticity into having only an component. This angle is defined only up to an additive integer multiple of .

We can multiply ellipticities together using the algebraic relations (37)


where . This product is rotationally invariant and describe the relative orientation of the ellipticities. Half its trace can be used as a scalar product.

Appendix B Two plane lens

To get a better intuitive feel for what is happening we can apply the equations in section 3 to the case where there are two lens planes. In this case, the position on the source plane () is


and the magnification is


The nonlinear term written out is


It is seen in these formulae that the magnification matrices of each of the lens planes acting alone add (50), but there is an additional nonlinear contribution (51) which gives rise to a rotation.

Appendix C transformation of the moments

For completeness we include here the transformation of an image’s moments by the gravitational lensing distortion in the tensor formalism. We feel that this is a more intuitive way to express these transformations than the more common ways in the literature (see Miralda-Escude, 1991; Kaiser & Squires, 1993; Bartelmann & Schneider, 2001).

The moments tensor of an image with surface brightness can be define as


where is the centroid of the image. The surface brightness is conserved but the coordinates are transformed


where it is assumed that the magnification matrix is constant over the source. Primes will indicate quantities for the original, pre-lensed source. Written out in our tensor notation this is


This can be readily worked out in terms of convergence and shear using equations (39), (40) and (41) along with the decomposition of  (11) (note the negative signs),


where is the reduced shear,


Equation (61) is the full transformation without averaging over the unknown orientation of the source or making the weak lensing approximation (while it is assumed that the shear and convergence are constant over the whole image). Averaging over the orientation of the source gives and, by the same argument, and . By expanding to first order in shear and convergence (the weak lensing approximation) and then averaging over the source’s orientation one obtains,


where has been used. For convenience and to comply with the common definition of ellipticity, we define and . 444It is common to define the components of ellipticity as and .

It is convenient to renormalize moments matrix and subtract the trace


where the average in the denominator is now over both the magnitude and the orientation on the source ellipticity. When is average in this way we get


Thus is an estimator of the weak shear. We find this expression simpler and more intuitive than the more common transformations using “polars” (Kaiser & Squires, 1993) or complex ellipticities (Seitz & Schneider, 1997).

Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description