GLAMER Part II: Multiple Plane Gravitational Lensing
Abstract
We present an extension to multiple planes of the gravitational lensing code GLAMER. The method entails projecting the mass in the observed lightcone onto a discrete number of lens planes and inverse rayshooting 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, powerlaw, etc.), haloes extracted from numerical simulations and clusters constructed from semianalytic models (MOKA). Likewise, there are several different options for modeling the source(s) which can be distributed throughout the lightcone. The distribution of matter in the lightcone can be either taken from a preexisting Nbody 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.
keywords:
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 energy^{1}^{1}1see DES: http://www.darkenergysurvey.org, PanStarrs :http://panstarrs.ifa.hawaii.edu 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.^{2}^{2}2HST Frontier Fields: http://www.stsci.edu/hst/campaigns/frontierfields/ 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 fourimages 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 lineofsight 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 lineofsight 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 quasistrong 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 Nbody 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 raybundles 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.
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 lightcone in angle and redshift. We also discuss some of the ways a lightcone 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 lineofsight objects to strong lenses, the importance of strongly lensed individual galaxies to weak lensing surveys, the statistical properties of strong lenses, the importance of multiplane 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 lightcones implemented in the code. In Sec. 5.2 we discuss the use a simulated halo catalog to populate the lightcone and in Sec. 5.3, the generation of a lightcone 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 galaxygalaxy 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 RobertsonWalker metric. The line element of this metric in the Newtonian gauge is
(1) 
where is the cosmological scale factor, is the Newtonian potential, is the speed of light, is the conformal time time, and
(2)  
(3) 
Here is the comoving distance and , the comoving angular diameter distance, is defined for different curvatures as follows:
(4) 
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
(5) 
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
(6) 
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,
(7) 
where is now the surface potential.
The surface potential is related to the surface density by
(8) 
where is the surface density within the slab which can be solved to get
(9) 
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
(10) 
This matrix is often decomposed into the convergence, and the components of shear . In the notation of appendix A, this decomposition is
(11) 
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
(12) 
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 multiplane lensing.
In this case, the source and image positions are related by
(13) 
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
(14) 
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 Bmodes. Also every point on the lens plane has a onetoone correspondence to a point on the image plane. As we shall see these things are not guaranteed when there are multiple lens planes.
3 multiplane formalism
In the multiplane 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
(15) 
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:
(16)  
(17) 
where is the redshift of the plane in consideration.
It is important to note here that since we calculate the angular diameter distances between planes assuming the RobertsonWalker metric, the mass added to those lenses planes will cause our lightcone to effectively contain more mass than the average for the FriedmanLemaitre 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)^{3}^{3}3This could be avoided by using a DryerRoeder 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
(18) 
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
(19) 
with the initial conditions
(20) 
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
(21) 
This is what the normal magnification matrix would be if the th plane where the source plane. Differentiating equation (19) yields
(22)  
(23) 
with the initial conditions
(24) 
The forcing term is
(25) 
where the notation of Appendix A is used. The coefficients are related to the Newtonian surface potential on the planes by
(26) 
The term (23) can be written out as
(27) 
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
Multiplane 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 Nbody simulations and semianalytic 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 – NavarroFrenkWhite (NFW) (Navarro et al., 1997), nonsingular isothermal sphere (NSIE), powerlaw, Hernquist (1990), pointmass, 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.
5 Generating Light Cones
In this section, we describe the different types of lightcones that are implemented within or can be input into the GLAMER code. The mass within the lightcone can be generated in four different ways: () particles can be taken directly from a Nbody simulation of a lightcone and sorted onto planes, () halos can be identified in the Nbody 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 Nbody simulation or semianalytic simulation, can be input into GLAMER. A lightcone 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 LightCones 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 stateoftheart numerical simulations the number of particles within a complete lightcone 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 Nbody 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 lightcone that is represented in a different way.
5.2 LightCones from Halo Catalog
GLAMER is able to read in a catalog of halos and their properties. These halos can be taken from an Nbody 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 Nbody 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 Nbody 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 Nbody 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 LightCones 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 Nbody simulation. No attempt has been made yet to cluster these halos or to create subhalos within halos for the full lightcone. 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);

PowerLaw (PL), given by the following equation:
(28) where represents the slope of the powerlaw, , and the nonlinear 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 tophat 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 lightcones 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 lightcone, 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 massconcentration relation (we refer the reader to the cited papers for more details): Zhao et al. (2009), MuñozCuartas et al. (2011), Giocoli et al. (2012b) and a powerlaw relation:
(29) 
where the redshift evolution factor is motivated by the change of the normalization as seen in cosmological numerical simulations.
This method of populating the lightcone 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 lightcones can be created.
As in the halos from Nbody 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):
(30) 
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 lightcone, 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 lightcone taken from a largescale 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 lightcones that are otherwise populated with halos. Fig. 7 shows and example of this.
5.5 Sources
Sources can be distributed throughout the lightcone. 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.
6 Tests
We perform several tests to check the accuracy and convergence of our multiplane 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 fourimaged 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 xdirection and 1.2 arcsec in the ydirection. 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 lightcone 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:
(31) 
where
(32) 
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 lightcone 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 lightcone. 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 multiplane lensing scheme. For example, we have tested that the average lensing mass on each plane is zero and that an empty lightcone with an arbitrary number of lensing planes always produces zero deflection. We also checked that an analytic lens placed in an empty lightcone 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 lightcone as described in section 5.3. The MOKA map was in this way inserted into the lightcone. 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 lightcone 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 lightcone 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 galaxygalaxy 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 lightcone 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.
8 discussion and conclusions
We have presented an extension of the rayshooting lensing code GLAMER to multiple planes and its ability to perform lensing simulations in a lightcone, 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 lightcone, 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 – PressSchechter, 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 lightcone 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.
Acknowledgments
This research is part of the project GLENCO, funded under the Seventh Framework Programme, Ideas, Grant Agreement n. 259349.
References
 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 eprints
 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 eprints
 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 eprints
 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
 MiraldaEscude (1991) MiraldaEscude 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ñozCuartas et al. (2011) MuñozCuartas 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 eprints
 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., CyrRacine F.Y., Fassnacht C. D., Keeton C. R., Linder E. V., Moustakas L. A., Bradac M., BuckleyGeer E., et al. 2013, ArXiv eprints
 Treu et al. (2013) Treu T., Marshall P. J., CyrRacine F.Y., Fassnacht C. D., Keeton C. R., Linder E. V., Moustakas L. A., et al. 2013, ArXiv eprints
 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 eprints
 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 2by2 matrix can be represented as
(33) 
where , and are scalar coefficients which may depend on position on the sky. The matrices are
(34)  
(35)  
(36) 
These matrices obey the following algebraic relations
(37)  
Some of the standard quantities and transformations are easily expressed in this basis:
trace  (38)  
determinant  (39)  
transpose  (40)  
(41) 
The rotation operator in this space is
(42) 
Using the algebraic relations, the behavior of these matrices under rotations can be written
(43) 
(44) 
(45) 
Thus and are rotationally invariant.
The ellipticity of a galaxy can be represented by a matrix with only and components
(46) 
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)
(47)  
(48) 
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
(49) 
and the magnification is
(50)  
(51) 
The nonlinear term written out is
(52) 
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 MiraldaEscude, 1991; Kaiser & Squires, 1993; Bartelmann & Schneider, 2001).
The moments tensor of an image with surface brightness can be define as
(53)  
(54) 
where is the centroid of the image. The surface brightness is conserved but the coordinates are transformed
(55)  
(56)  
(57)  
(58) 
where it is assumed that the magnification matrix is constant over the source. Primes will indicate quantities for the original, prelensed source. Written out in our tensor notation this is
(59)  
(60) 
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),
(61)  
where is the reduced shear,
(62) 
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,
(63)  
(64)  
(65) 
where has been used. For convenience and to comply with the common definition of ellipticity, we define and . ^{4}^{4}4It is common to define the components of ellipticity as and .
It is convenient to renormalize moments matrix and subtract the trace
(66) 
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
(67) 
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).