Raytracing through the Millennium Simulation: Born corrections and lenslens coupling in cosmic shear and galaxygalaxy lensing
S. Hilbert^{1,2}  J. Hartlap^{2}  S. D. M. White^{1}  P. Schneider^{2}
1  Max Planck Institute for Astrophysics, KarlSchwarzschildStr. 1, 85741 Garching, Germany
2  ArgelanderInstitut für Astronomie, Universität Bonn, Auf dem Hügel 71, 53121 Bonn, Germany
Received 29 September 2008 / Accepted 13 March 2009
Abstract
Context. Weaklensing surveys need accurate theoretical predictions for interpretation of their results and cosmologicalparameter estimation.
Aims. We study the accuracy of various approximations to cosmic shear and weak galaxygalaxy lensing and investigate effects of Born corrections and lenslens coupling.
Methods. We use raytracing through the Millennium Simulation, a large Nbody simulation of cosmic structure formation, to calculate various cosmicshear and galaxygalaxylensing statistics. We compare the results from raytracing to semianalytic predictions.
Results. (i) We confirm that the firstorder approximation (i.e. neglecting lensing effects beyond first order in density fluctuations) provides an excellent fit to cosmicshear power spectra as long as the actual matter power spectrum is used as input. Common fitting formulae, however, strongly underestimate the cosmicshear power spectra (by
on scales
). Halo models provide a better fit to cosmic shearpower spectra, but there are still noticeable deviations (
). (ii) Cosmicshear Bmodes, which are induced by Born corrections and lenslens coupling, are at least three orders of magnitude smaller than cosmicshear Emodes. Semianalytic extensions to the firstorder approximation predict the right order of magnitude for the Bmode. Compared to the raytracing results, however, the semianalytic predictions may differ by a factor two on small scales and also show a different scale dependence. (iii) The firstorder approximation may under or overestimate the galaxygalaxylensing shear signal by several percent due to the neglect of magnification bias, which may lead to a correlation between the shear and the observed number density of lenses.
Conclusions. (i) Current semianalytic models need to be improved in order to match the degree of statistical accuracy expected for future weaklensing surveys. (ii) Shear Bmodes induced by corrections to the firstorder approximation are not important for future cosmicshear surveys. (iii) Magnification bias can be important for galaxygalaxylensing surveys.
Key words: gravitational lensing  dark matter  largescale structure of Universe  cosmology: theory  methods: numerical
1 Introduction
During the past few years, weak gravitational lensing has developed rapidly from mere detection to an important cosmological tool (Munshi et al. 2008). Measurements of cosmic shear help us to constrain the properties of the cosmic matter distribution (e.g. Semboloni et al. 2006; Benjamin et al. 2007; Hoekstra et al. 2006; Simon et al. 2007; Fu et al. 2008; Massey et al. 2007b), the growth of structure (e.g. Massey et al. 2007c; Bacon et al. 2005), and the nature of the dark energy (e.g. Schimd et al. 2007; Amendola et al. 2008; Taylor et al. 2007). Weak galaxygalaxy lensing can be used to study the properties of galactic darkmatter halos and the relation between luminous and dark matter (e.g. Gavazzi et al. 2007; Mandelbaum et al. 2006; Simon et al. 2007).
The accuracy that can be reached in weaklensing surveys is determined by several factors. On the observational side, high accuracy requires large field sizes and deep observations with a high number density of galaxies with measurable shapes. Moreover, it is crucial to obtain an accurate and unbiased measurement of galaxy ellipticities. Finally, for the interpretation of the resulting data and the inference of cosmological parameters, an accurate theoretical model is needed. A thorough understanding of systematic effects in weak lensing will become particularly important with the advent of very large weaklensing surveys such as CFHTLS^{}, KIDS^{}, PanSTARRS^{}, and LSST^{}, or the planned Dark Energy Survey^{}, DUNE (Réfrégier et al. 2006), and SNAP^{}. For these surveys, the statistical uncertainties will be very small, so the accuracy will be limited by the remaining systematics in the data reduction and theoretical modeling.
While significant improvement on imageellipticity measurements are expected in the near future (Massey et al. 2007a), one still needs to investigate, how uncertain current theoretical predictions are, and how much improvement can be expected for these. Presently, the most accurate way to obtain predictions for weaklensing surveys is to perform raytracing through large highresolution Nbody simulations of cosmic structure formation (see, e.g., Van Waerbeke et al. 2001; Vale & White 2003; Jain et al. 2000; White 2005; White & Hu 2000; Wambsganss et al. 1998; Hamana & Mellier 2001). The drawback of this approach is that large Nbody simulations are computationally demanding, so using them to explore the whole parameter space of cosmological models is currently unrealistic. On the other hand, raytracing simulations enable one to check the approximations and assumptions made in computationally less demanding (semi)analytic models, and adjust and extend these models where necessary.
Numerous raytracing methods have been developed to study the many aspects of gravitational lensing. Treebased raytracing methods (Aubert et al. 2007) that adapt to the varying spatial resolution of Nbody simulations have been used to study the impact of substructure on strong lensing by dark matter halos (Peirani et al. 2008). Cluster strong lensing simulations, which require good mass modelling of galaxy clusters, usually ignore the matter distribution outside clusters and use the thinlens approximation in the raytracing (e.g. Bartelmann & Weiss 1994; Meneghetti et al. 2007; Rozo et al. 2008).
Many simulations of weak lensing by clusters and largescale structure (e.g. Pace et al. 2007; Vale & White 2003; Jain et al. 2000; Wambsganss et al. 1998) employ algorithms that are based on the multiplelensplane approximation (Blandford & Narayan 1986) to trace light rays through cosmological Nbody simulations. Others (e.g. Carbone et al. 2008; Couchman et al. 1999) perform raytracing though the threedimensional gravitational potential. In a simpler approach (e.g. Hilbert et al. 2007a; White & Vale 2004; Heymans et al. 2006), the matter in the Nbody simulation is projected along unperturbed light paths onto a single lens plane, which is then used to calculate lensing observables. Recent simulations of CMB lensing use generalisations of the single or multipleplane approximation that take the curvature of the sky into account (e.g. Das & Bode 2008; Fosalba et al. 2008; Teyssier et al. 2008).
In this work, we employ multiplelensplane raytracing through the Millennium Simulation (Springel et al. 2005) to study weak lensing. One of the largest Nbody simulations available today, the Millennium Simulation provides not only a much larger volume, but also a higher spatial and mass resolution than simulations used for earlier weaklensing studies. In order to take full advantage of the large simulation volume and high resolution, the raytracing algorithm used here differs in several aspects from algorithms used in previous works (e.g. Jain et al. 2000). Here, we give a detailed description of our raytracing algorithm.
Semianalytic weaklensing predictions are usually based on the firstorder approximation, in which light deflections are only considered to first order in the peculiar gravitational potential and hence, to first order in the matter fluctuations. The raytracing approach allows us to look at effects neglected in the firstorder approximation such Born corrections and lenslens coupling. Here, we investigate the cosmicshear Bmodes induced by these effects and compare the raytracing results to semianalytic estimates (Shapiro & Cooray 2006; Cooray & Hu 2002; Hirata & Seljak 2003), whose accuracy has not been confirmed by numerical simulations yet. Moreover, we investigate how well fitting formulae (Peacock & Dodds 1996; Smith et al. 2003; Eisenstein & Hu 1999) and halo models (Cooray & Sheth 2002; Seljak 2000) reproduce cosmicshear power spectra. Finally, we investigate the accuracy of the firstorder approximation for weak galaxygalaxy lensing.
The paper is organised as follows. In the next section, we introduce the theoretical background and notation used in our lensing analysis. In Sect. 3, we discuss our raytracing algorithm. The results from our raytracing analysis are presented in Sect. 4. We conclude our paper with a summary in Sect. 5.
2 Theory
2.1 Gravitational light deflection
In this section, we introduce the formulae relating the ``apparent'' positions of distant light sources to their ``true'' positions. In order to label spacetime points in a model universe with a weakly perturbed FriedmannLemaîtreRobertsonWalker (FLRW) metric, we choose a coordinate system
based on physical time t, two angular coordinates
,
and the lineofsight comoving distance w relative to the observer. The spacetime metric of the model universe is then given by (see, e.g., Bartelmann & Schneider 2001):
Here, c denotes the speed of light, a=a(t) denotes the scale factor, and denotes the peculiar (Newtonian) gravitational potential. The comoving angular diameter distance is defined as:
where K denotes the curvature of space. The particular choice for the angular coordinates is convenient for the application of the ``flatsky'' approximation, where the metric near is approximated using .
Consider the path, parametrised by comoving distance w, of a photon eventually reaching the observer from angular direction
.
The angular position
of the photon at comoving distance w is then given by (see, e.g., Jain & Seljak 1997, for a sketch of a derivation):
with , and t(w') denoting the cosmic time of events at lineofsight comoving distance w' from the observer. By differentiation this equation w.r.t. , we obtain the distortion matrix , i.e. the Jacobian of the lens mapping :
Due to the matrix products in Eq. (4), the distortion matrix is generally not symmetric. However, it can be decomposed into a rotation matrix (related to a usually unobservable rotation in the source plane) and a symmetric matrix (Schneider et al. 1992):
The decomposition defines the rotation angle , the convergence , and the two components and of the shear, which may be combined into the complex shear .
The shear field
can be decomposed into a rotationfree part
and a divergencefree part
.
For infinite fields, the decomposition into these E/Bmodes is most easily written down in Fourier space:
Here, hats denote Fourier transforms, denotes the Fourier wave vector, and . Care must be taken when decomposing the shear in fields of finite size, where the field boundaries can cause artifacts (Seitz & Schneider 1996). These artifacts can be avoided by using aperture masses to quantify the shear E and Bmode contributions (Schneider et al. 2002; Crittenden et al. 2002).
Equations (3) and (4) are implicit relations for the light path and the Jacobian. The solution of Eq. (3) to first order in the potential is obtained by integrating along undisturbed light paths:
The distortion to first order reads:
The firstorder approximation to the distortion contains the Born approximation, which ignores deviations of the actual light path from the undisturbed path on the r.h.s. of Eq. (4). Moreover, lenslens coupling is neglected, i.e. the appearance of the distortion on the r.h.s. of Eq. (4). The neglected lenslens coupling and corrections to the Born approximation account for the effect that light from a distant source ``sees'' a distorted image of the lowerredshift matter distribution due to higherredshift matter inhomogeneities along the lineofsight. Thus, the firstorder approximation works well in regions where larger matter inhomogeneities are absent or confined to a small redshift range, but fails in regions where noticeable distortions arise from matter inhomogeneities at multiple redshifts.
Born corrections and lenslens coupling effects may create shear Bmodes. The perturbative calculation of the shear Bmodes by iteratively solving Eq. (4) is possible (Hirata & Seljak 2003; Cooray & Hu 2002), but tedious, and the accuracy of this approach is not known. However, multiple deflections and lenslens coupling effects are fully included in the multiplelensplane approximation as described below. We will thus use this approximation to investigate these effects and assess the quality of perturbative calculations of these effects.
2.2 The multiplelensplane approximation
In the multiplelensplane approximation (see, e.g., Blandford & Narayan 1986; Schneider et al. 1992; Seitz et al. 1994; Jain et al. 2000), a series of lens planes perpendicular to the central lineofsight is introduced into the observer's backward light cone. The continuous deflection that a light ray experiences while propagating through the matter inhomogeneities in the light cone is then approximated by finite deflections at the lens planes. The deflections are calculated from a projected matter distribution on the lens planes. This corresponds to solving the integral Eqs. (3) and (4) by discretisation (and using the impulse approximation).
The deflection
of a light ray intersecting the
lens plane (here, we count from the observer to the source) at angular position
can be expressed as the gradient of a lensing potential
:
The differential deflection is then given by higher derivatives of the lensing potential. The second derivatives can be combined into the shear matrix
(10) 
The lensing potential is a solution of the Poisson equation:
(11) 
The dimensionsless surface mass density is given by a projection of the matter distribution in a slice around lens plane:
(12) 
Here, H_{0} denotes the Hubble constant, the mean matter density in terms of the critical density, f_{K}^{(k)}=f_{K}(w^{(k)}) and a^{(k)}=a(w^{(k)}), with w^{(k)} denoting the lineofsight comoving distance of the plane. Furthermore, denotes the threedimensional density contrast at comoving position relative to the mean matter density. The slice boundaries and have to satisfy and . They are usually chosen to correspond to the mean redshifts (e.g. Jain et al. 2000) or comoving distances (e.g. Wambsganss et al. 2004) of successive planes^{}. These conditions ensure that every region of the light cone contributes exactly to one lens plane, which is the closest plane in redshift or comoving distance.
Figure 1: Schematic view of the observer's backward light cone in the multiplelensplane approximation. A light ray (red line) experiences a deflection only when passing through a lens plane (solid blue lines). The deflection angle of a ray passing through the lens plane at distance f_{K}^{(k1)} from the observer is obtained from the matter distribution between and projected onto the plane. Using the deflection angle of the light ray at the previous lens plane and the ray's angular positions and on the two previous planes, the angular position on the current plane can be computed. 

Open with DEXTER 
Given the deflection angles on the lens planes, one can trace back a light ray reaching the observer from angular position
on the first lens plane to the other planes:
Here, .
Equation (13) is not practical for tracing rays through many lens planes. An alternative expression is obtained as follows see, e.g., Hartlap 2005; or Seitz et al. 1994 for a different derivation: the angular position of a light ray on the lens plane k is related to its positions and on the two previous lens planes by (see Fig. 1):
Hence,
For a light ray reaching the observer from angular position on the first lens plane, one can compute its angular position on the other lens planes by iterating Eq. (15) with initial values .
Differentiating Eq. (15) with respect to
,
we obtain a recurrence relation for the distortion matrix:
With the knowledge of the involved distances and shear matrices, this equation allows us to iteratively compute the distortion matrix of a light ray from the observer to any lens plane. This equation requires in practice much fewer arithmetic operations and memory than the commonly used relations (e.g. by Jain et al. 2000) based on Eq. (13).
For comparison and testing, we will also use the multiplelensplane algorithm to calculate the distortion in the firstorder approximation by:
3 The raytracing algorithm
The methods we use for raytracing through Nbody simulations to study lensing are generally similar to those used by, e.g., Jain et al. (2000) or Vale & White (2003). First, the matter distribution on the past light cone of a fiducial observer is constructed from the simulation data. Then, the past light cone is partitioned into a series of redshift slices. The content of each slice is projected onto a lens plane. Finally, the multiplelensplane approximation is used to trace back light rays from the observer through the series of lens planes to the sources.
The purpose of our raytracing algorithm is to simulate strong and weak lensing in a way that takes full advantage of the unprecedented statistical power offered by the large volume and high spatial and mass resolution of the Millennium Simulation^{}. Therefore, our raytracing method differs in many details from previous works. Most notably, we use a multiplemesh method and adaptive smoothing to calculate light deflections and distortions from the projected matter distribution on the lens planes. This allows us to simulate lensing on the full range of scales covered by the Millennium Simulation, ranging from strong lensing on scales to cosmic shear on scales . A brief outline of our algorithms for the construction of the past light cones and the lens planes has been given in an earlier work (Hilbert et al. 2007b). Here, we extend the discussion and provide a more detailed description.
3.1 The Millennium Simulation
The Millennium Simulation (Springel et al. 2005) is a large Nbody simulation of cosmic structure formation in a flat CDM universe. The following cosmological parameters were assumed for the simulation: a matter density of in units of the critical density, a cosmological constant with , a Hubble constant h=0.73 in units of , a spectral index n=1 and a normalisation parameter for the primordial linear density power spectrum. These chosen parameters are consistant with the 2dF (Colless et al. 2003) and WMAP 1styear data analysis (Spergel et al. 2003). The simulation followed the evolution of the matter distribution in a cubic region of comoving side length from redshift z=127 to the present using a TreePM version of GADGET2 (Springel 2005) with 2160^{3} particles of mass and a force softening length of comoving.
Snapshots of the simulation were stored on disk at 64 output times. These snapshots contain, among other data, the positions, SPH smoothing lengths and friendoffriend (FoF) group data of the particles. The storage order for the particle data is based on a spatial octtree decomposition of the simulation cube (see Springel et al. 2005; Springel 2005, for details), which facilitates access to the particle data for small subvolumes of the simulation.
Complex physical processes of baryonic matter such as the formation and evolution of stars in galaxies has not been incorporated directly into the Millennium Simulation. However, several galaxyformation models have been used to predict the properties of galaxies in the simulation (Springel et al. 2005; Bower et al. 2006; De Lucia & Blaizot 2007; Croton et al. 2006). The raytracing methods presented in this paper will allow us (in future work) not only to study cosmic shear in great detail, but also to make predictions for galaxygalaxy lensing (and related higherorder statistics) for the various galaxyformation models.
3.2 The construction of the matter in the backward light cones
Even with a comoving size of , the simulation box is too small to trace back light rays within one box. We therefore exploit the periodic boundary conditions of the simulation by arranging replicas of the simulation box in a simple cubic lattice with a lattice constant equal to the box size L to fill space. We refrain from randomly shifting or rotating the content of the lattice cells, because the simulation box is far too large to be projected onto a single lens plane. In addition, this allows us to keep the matter distribution continuous across the cell boundaries.
In this periodic matter distribution, light rays would encounter the same structures many times at different epochs before reaching relevant source redshifts if one chose the line of sight (LOS) parallel to the box edges. Hence, the LOS must be chosen at a skewed angle relative to the box axes. On the other hand, the application of Fourier methods for the calculation of the light deflection at the lens planes requires a matter density that is periodic perpendicular to the LOS. Choosing a LOS parallel to with suitable coprime , one can obtain a large enough repetition length of along the LOS (see Appendix A). At the same time, the matter distribution is periodic perpendicular to the LOS with an area of periodicity given by . Our choice for yields a LOS periodicity of (corresponds to z=3.87) and a rectangular unit cell of for the lens planes. Moreover, any directions with shorter periodicity are at least away from , and a light cone with a field of view does not intersect with itself up to redshift z=3.87 when folded back into the simulation cube^{}. The resulting orientation of the LOS and the lens planes w.r.t. the simulation box are illustrated in Fig. 2.
We partition the observer's backward light cone into redshift slices such that each slice contains the part of the light cone that is closer in redshift to one of the snapshots than any other snapshot (with exceptions near the boundaries discussed below). The boundary between two redshift slices with snapshot redshifts z^{(k)} and z^{(k+1)} is thus a plane at comoving distance . In addition, . The particle data of the snapshot closest in redshift is then used to approximate the matter distribution in each of these slices. Fast boxintersection tests (Gottschalk et al. 1996) and the spatialocttree storage order of the simulation are utilised to minimise reading of the particle data (which reduces run time by factors 510).
In the construction of the matter distribution in the light cone, special care is taken for the particles near the boundary of two slices. In the simulation, particle concentrations representing dark matter halos of galaxies or clusters were identified by a friendoffriend (FoF) group finding algorithm. Some of these halos are located on the slice boundaries with particles on either side^{}. In order to avoid that such a halo is only partially included into a slice (and hence would be only partially projected onto a lens plane), a halo is either included as a whole if its central particle is inside the slice as defined by boundary planes, or completely excluded otherwise.
If the matter structure in the simulation were static, this procedure would suffice to prevent parts of the same halo from being projected onto adjacent lens planes, which would create artificial close pairs of halos on the sky. Halos, however, may have moved across a slice boundary between two snapshots. We therefore amend the above inclusion criterion for halos near the boundary: if a halo is included in (excluded from) the slice of the later snapshot, its progenitors in the earlier snapshot are excluded from (included in) the ``earlier'' slice even if their centres lie on the ``early'' (``late'') side of the slice boundary. These inclusion criteria for halos are illustrated in Fig. 3.
Figure 2: Schematic view of the orientation of the lineofsight (red line) and the lens planes (blue area) relative to the simulation box (indicated by black lines). 

Open with DEXTER 
Figure 3: Schematic view of the adaptive slice boundaries to avoid the truncation or double inclusion of halos that are located near a slice boundary. Halos near the boundary of slice k and k+1 are eitherincluded as a whole in slice k or completely excluded depending on the positions of their centres (a). Halos that are included (excluded) in slice k, are excluded (included) from slice k+1 even if they have crossed the slice boundary between redshift k and k+1 (b). 

Open with DEXTER 
3.3 The lens planes
The matter content of each redshift slice of the backward light cone is projected along the LOS direction onto a lens plane. Each lens plane is placed at the comoving distance of the corresponding snapshot's redshift. The lens planes serve also as source planes for the raytracing. The resulting number of lens planes as a function of the source redshift is shown in Fig. 4.
Figure 4: The number of lens planes used for the raytracing as afunction of the source redshift . 

Open with DEXTER 
The light deflection angles and distortions resulting from the projected matter density on the lens planes are computed by particlemesh (PM) methods (Hockney & Eastwood 1981). Mesh methods have the advantage that, once the deflection and distortion are computed on a mesh (e.g. by Fast Fourier methods), the computation of the deflections and distortions for many light rays intersecting the plane is very fast (compared to, e.g., directsummation or tree methods). One disadvantage is that the used mesh spacing limits the spatial resolution of the projected matter distribution. However, any Nbody simulation providing the matter distribution for the raytracing has a limited resolution as well. In dense regions, the spatial resolution of the Millennium Simulation is effectively determined by the force softening, which is comoving. Thus, a mesh spacing of comoving is required to avoid resolution degradation for the projected matter density. However, a single mesh covering the full periodic area of the lens plane (i.e. comoving) with such a small mesh spacing would be too demanding, in particular regarding the memory required both for its computation and storage. We therefore use a hierarchy of meshes instead.
The lensing potential is split into a longrange part and a shortrange part . The split is defined in Fourier space by:
(18)  
(19) 
The splitting angle , with comoving splitting length and comoving angular diameter distance of the lens plane f_{K}(w), quantifies the spatial scale of the split. Different meshes are then used to calculate and .
First, the particles in each slice are projected onto a coarse mesh of 16 384 16 384 points covering the whole periodic area of the lens plane using cloudsincells (CIC) assignment (Hockney & Eastwood 1981). The longrange potential is then calculated on this mesh by means of fast Fourier transform (FFT) techniques (Cooley & Tukey 1965; Frigo & Johnson 2005). The splitting length is chosen slightly larger than the coarse mesh spacing ( and comoving, respectively), so the coarse mesh samples with sufficient accuracy. For each lens plane, the longrange potential is calculated once, and the result is stored on disk for later use during the raytracing.
The shortrange potential is calculated ``on the fly'', i.e. during the actual raytracing. The area where the light rays intersect the plane is determined and, if larger than comoving, subdivided into several patches up to that size. Each patch is covered by a fine mesh with a mesh spacing of comoving and up to 16 384 16 384 mesh points. The fine meshes are chosen slightly larger than the patches in order to take into account all matter within the effective range of , for which we assume ( ). The limited range of ensures that the matter distribution outside the mesh affects only mesh points close to its boundary (i.e. within the effective range), but not the interior mesh points used for the subsequent analysis. Periodic boundary conditions can therefore be used for the FFT on the patches without ``zero padding''.
In order to reduce the shot noise from the individual particles, either a fixed or an adaptive smoothing scheme is used for the matter distribution on the fine meshes. In case of the fixed smoothing, the particles in the slice are projected onto the fine mesh using CIC. The resulting matter density on the fine mesh is then smoothed in Fourier space with a Gaussian lowpass filter
(20) 
whose filter scale is determined by the lens plane's comoving distance w and a fixed comoving filter scale . This is done during the calculation of the shortrange potential with FFT methods.
In case of the adaptive smoothing, the mass associated with each simulation particle contributes
(21) 
to the surface mass density on the fine mesh. Here, denotes comoving position on the lens plane, is the projected comoving particle position, and denotes the comoving distance to the 64th nearest neighbour particle in three dimensions (i.e. before projection). The adaptive smoothing is essentially equivalent to the assumption that, in threedimensional space, each simulation particle represents a spherical cloud with a Gaussian density profile and an rms radius that is half the distance to its 64th nearest neighbour. From the resulting surface mass density on the fine mesh, the shortrange potential is then calculated by FFT methods.
Figure 5: The deflection angle as a function of the comoving distance r to a point mass computed by our mesh method. a) Shortrange component (dashed line) and longrange component (dotted line) of the deflection angle (full line). b) Fractional difference between the value calculated by our mesh method and the theoretical value . For the plots, we placed 10 point masses of 8.6 randomly on the lens plane at z=0.5 and calculated the resulting deflection at 1000 random positions around each of them. 

Open with DEXTER 
Figure 6: The magnification for sources at z=2 as a function of the angular separation to a point mass of at redshift z=1 computed by our mesh method. a) Numerical values (symbols) compared to the analytical result (solid line). b) Fractional difference between the measured magnification and its theoretical value . The magnification diverges at the Einstein radius (dashed vertical line). For the plots, we placed 10 point masses of randomly on the lens plane at z=1 and calculated the resulting magnification at 1000 random positions around them. 

Open with DEXTER 
The long and shortrange contributions to the deflection angles and shear matrices are calculated on the coarse and fine mesh by finite differencing of the potentials^{}. The values between mesh points are obtained by bilinear interpolation. The resulting deflection angles and shear matrices at the ray positions are then used to advance the rays and their associated distortion matrices from one plane to the next.
The raytracing algorithm reproduces the deflection angles and distortions caused by a single point mass very accurately, as is shown in Figs. 5 and 6. For the deflection angle, the relative deviation from the known analytical result is at most one percent, with a much smaller rms error. Apart from scales below comoving, where discreteness of the fine mesh becomes important, the largest relative errors occur on the scale of the split between short and longrange potential. If desired, an even smaller error in this region could be obtained by increasing the splitting scale.
In this work, we do not consider the effects of the stellar mass in galaxies. Note, however, that the raytracing algorithm can be extended to include the gravitational effects of the stars in galaxies as described in Hilbert et al. (2008): the positions and stellar masses of the galaxies are inferred from semianalytic galaxyformation models implemented within the evolving darkmatter distribution of the Millennium Simulation (Springel et al. 2005; De Lucia & Blaizot 2007; Croton et al. 2006). The light deflection induced by the stellar matter is then calculated using analytic expressions for the projected stellar mass profiles on the lens planes. The error made by adding the stellar matter onto the lens planes, although the darkmatter particles of the simulation represent the total matter, can then be compensated using extended analytic profiles with negative masses (as was done, e.g., by Wambsganss et al. 2008).
3.4 Lensing maps and image positions
To produce lensing maps, we set up rays starting from a fiducial observer on a regular grid in the image plane with an angular field size and resolution suitable for the particular application. The resulting shear and convergence maps may be used directly to obtain, e.g., the shear correlation functions and power spectra.
We also wish to perform simulations of galaxygalaxy lensing. Not only are the image positions of distant source galaxies affected by lensing. Also the apparent positions of galaxies and halos that act as lenses for background galaxies (and are to be correlated with the shear field) are affected by lensing due to foreground matter. We therefore have to compute the image positions given the galaxies' source positions (i.e. the projected galaxy positions on the lens planes) and the lens mapping sampled on the grid of light rays in the raytracing algorithm. To reach this, we make use of a triangle technique described, e.g., in Schneider et al. (1992). We partition the region of the image plane that is covered by the grid of rays into triangles formed by rays of adjacent grid points (see Fig. 7). On each lens plane, we identify for each such triangle all galaxies with source position inside the backtraced triangle. The image positions of these galaxies are then computed by linear interpolation between the rays. This method takes into account strong lensing, as a galaxy might lie in more than one triangle on the lens plane, resulting in multiple images of that galaxy.
Figure 7: Interpolation scheme used for determining image positions of galaxies. The regular grid of rays in the image plane (left filled circles) is used to partition the image plane into triangles (right blue lines). The image position (left open circle) of a source inside a triangle (right blue lines) formed by the backtraced rays on the source plane (right filled circles) is then determined by linear interpolation. 

Open with DEXTER 
Figure 8 illustrates how well the image positions obtained by the triangle interpolation method are mapped back onto the source positions by the raytracing. Shown is the difference between the true source positions and the positions obtained by tracing back light rays starting from the interpolated image positions to the source plane. The difference is always much smaller than the resolution of the matter distribution on the lens planes. The slight anisotropy seen as a larger spread along one diagonal is caused by the particular way the image plane was partitioned into triangles. The diagonal coincides with the diagonal chosen for splitting the square mesh pixels into triangles. If one used the other diagonal for splitting the mesh pixels into triangles, a larger spread along that diagonal would be seen instead.
Figure 8: Comparison of the true source positions and the source positions obtained by tracing back light rays through the Millennium Simulation starting from the image positions computed by the interpolation method. Shown is the difference between true and tracedback comoving source position for 1000 galaxy centres at z=1. At this redshift, comoving correspond to an angle of . The right and upper frames sides are labelled by the corresponding angular difference between true and tracedback source position in units of the mesh spacing used for the rays. 

Open with DEXTER 
4 Results
We compute various weaklensing twopoint statistics from a set of raytracing simulations and compare the results to semianalytic predictions. If not stated otherwise, adaptive snoothing of the matter distribution on the lens planes is applied for the raytracing.
4.1 Power spectra
We start our discussion with the convergence power spectrum
.
In the firstorder approximation (8), the convergence power spectrum is given by (see, e.g., Schneider 2006):
Here,
with the probability distribution of visible sources in comoving distance. Furthermore, denotes the threedimensional power spectrum of the matter contrast at cosmic time t and comoving wave vector k. For an accurate comparison with the results obtained by our raytracing algorithm, we use Eq. (22) together with the threedimensional power spectra measured from the Millennium Simulation (see also Vale & White 2003). In the following, we will call the resulting power spectra firstorder prediction for brevity.
Figure 9: Convergence power spectra for sources at redshift z=1 (lower curves) and z=2 (upper curves). The simulation results (from random fields of 3 ) are shown as diamonds with errorbars (indicating standard deviation calculated from the fieldtofield variance), the corresponding firstorder predictions as solid lines. The predictions using the Peacock & Dodds (1996) prescription together with the transfer function from Eisenstein & Hu (1999) are given as dotted lines, those obtained from the Smith et al. (2003) fitting formula as dashed lines. The predictions of a halo model using the concentration parameters of Neto et al. (2007) are shown as dashdotted lines. 

Open with DEXTER 
In Fig. 9, we compare the firstorder prediction to the convergence power spectra obtained from the raytracing. As has already been observed by Jain et al. (2000) and Vale & White (2003), the power spectra from the raytracing are in very good agreement with the firstorder prediction. For the considered source redshifts, the difference for . The larger deviations at wave numbers 10^{4} are due to smoothing effects discussed below.
In the firstorder approximation, the power spectra of the convergence and the shear are identical. As has already been found, e.g., by Jain et al. (2000), the convergence and shear power spectra from raytracing agree very well, too. On scales , the difference between both is well below one percent in our raytracing results.
If the firstorder prediction for the convergence powerspectra is assumed to be correct to very high accuracy, the smoothing tests can be considered as a test of the accuracy of our raytracing algorithm. Then the results shown in Fig. 9 suggest that the raytracing is able to reproduce weaklensing effects within accuracy on scales .
The comparison of the raytracing power spectra with some of the popular fitting formulae is less encouraging: Both the prescriptions by Peacock & Dodds (1996) (with the transfer function by Eisenstein & Hu 1999) and Smith et al. (2003) strongly underpredict the power on intermediate and small scales. These fitting formulae are based on older simulations, whose matter power spectra are noticeably different from the power spectra of more recent, higherresolution simulations. The deviations from the simulated convergence power spectra exceed 30% for , so these fitting formulae seem to be of limited use for the interpretation of data from future weaklensing surveys.
A prediction based on the popular halo model (Cooray & Sheth 2002; Seljak 2000) and the halo concentrationmass relation of Neto et al. (2007) provides a better fit to the convergence power spectrum. There are, however, still deviations ( ), in particular for higher source redshifts and intermediate scales (i.e. 10^{3}). This coincides with the transition region of the one and twohalo terms, which is difficult to model accurately due to halo exclusion effects (see e.g. Tinker et al. 2005, and references therein), which are not included in our prediction.
Figure 10: Convergence power spectra for sources at redshift z=1. Compared are the results from raytracing (symbols) using various smoothing schemes (none/Gaussian with fixed scale /adaptive) and the corresponding firstorder prediction (lines) obtained projecting and smoothing the measured 3D power spectra of the actual mass distribution in the simulation. 

Open with DEXTER 
As mentioned above, the deviations of the measured power spectra and the firstorder predictions at large are due to smoothing effects. In Fig. 10, we present the convergence power spectra from raytracing runs of the same set of fields (with a cumulative area of and sources at z=1), but with different smoothing schemes. In addition to adaptive smoothing, which is intractable analytically, we also employ smoothing with a Gaussian kernel of fixed comoving size on the lens planes. The raytracing simulations with Gaussian smoothing on the lens planes show  apart from sampling variance  perfect agreement with the firstorder prediction if the smoothing is into taken into account there. Only the spectrum for the smallest smoothing length shows some aliasing effects on very small scales. The spectrum of the adaptivesmoothing runs happens to match the spectrum for a Gaussian smoothing length of comoving quite well, but one should be cautious when considering this as an ``effective'' smoothing length in a different context.
4.2 Aperturemass statistics
A suitable cosmicshear measure that allows one to decompose the shear signal in a finitesized field into E and Bmodes is the aperture mass dispersion (Schneider et al. 1998,2002). The E and Bmode aperture mass at position
on the sky and scale
are defined by:
In this work we use the polynomial filter function Q proposed by Schneider et al. (1998):
(25) 
The tangential and cross components of the shear are defined by
(26a)  
(26b) 
where is the polar angle for the direction defined by .
Figure 11: Aperture mass dispersion a) and b) as a function of filter scale for sources at z=1 and z=2. Compared are the firstorder prediction (solid lines, Emode only) and the results from raytracing (symbols with error bars indicating standard deviation, obtained from 7 fields of 5 ). 

Open with DEXTER 
Figure 12: Bmode aperture mass dispersion for sources at z=2 decomposed into the contributions by lenslens coupling and Born corrections. Raytracing results (symbols with error bars indicating standard deviation, calculated from 7 fields of 1 ): full raytracing (squares), only Born corrections (diamonds), only lenslens coupling (triangles). Predictions by Cooray & Hu (2002): full signal (dashed line), only Born corrections (dotted line), only lenslens coupling (dashdotted line). Prediction by Hirata & Seljak (2003): full signal (solid line). 

Open with DEXTER 
An estimate for the aperture mass dispersion
as a function of the filter scale
can be computed from a given shear field by a spatial average. Figure 11a shows the Emode aperture mass dispersion measured from our set of simulations. The dispersion measured from the raytracing is in very good agreement with the firstorder prediction (Schneider et al. 1998):
where is given by Eq. (22), and J_{4} denotes a Bessel function of the first kind. The deviations of the measured aperture mass dispersion from the firstorder prediction seen on scales can be attributed to smoothing.
In the firstorder approximation, the Bmode aperture mass dispersion vanishes. The measured Bmode dispersion from the full raytracing is shown in Fig. 11b. The Bmode signal is at least 3 orders of magnitude smaller than the Emode. On larger scales their ratio even drops below 10^{5}.
Theoretical predictions of the amplitude of the lensinginduced Bmodes have been made by Cooray & Hu (2002) and Hirata & Seljak (2003), who calculated corrections to the E and Bmode shear power spectra by expanding Eq. (4) to second order in the gravitational potential. As Fig. 12 illustrates, the predictions based on their methods (and the measured threedimensional power spectra of the Millennium Simulation) are of the correct order of magnitude and reproduce some qualitative features of the raytracing simulations, but the match is far from being perfect. While the Bmode predictions are lower by a factor of on small scales, the signal measured from the raytracing declines much more quickly on larger scales. However, the discrepancies are not large enough to challenge the finding of Shapiro & Cooray (2006) that the lensinginduced Bmode is unimportant even for an allsky survey.
In order to determine their individual contributions to the total Bmode signal, we switch off ray deflections (i.e. we employ the Born approximation by setting in Eq. (15)) and/or lenslens coupling (i.e. we set in the third term on the r.h.s. of Eq. (16)). Again the predictions and the measured signal differ by factors on small scales, and the measured signal decreases much stronger with increasing scale.
We closely examined various steps involved in the calculations to exclude numerical artifacts as the reason for the discrepancy. The smoothing tests show that the smoothing we applied to the matter distribution on the lens planes can only account for deviations on scales . Examining the variance between the different raytraced fields, we can exclude ``cosmic variance'' as a major source of the discrepancy.
The raytracing results did not change when different ways of estimating in real and Fourier space, as well as different methods of numerical integration for the theoretical curves were used. Furthermore, only a tiny Bmode (at least 6 orders of magnitudes smaller than the Emode) remained, when both ray deflections and lenslens coupling were switched off in the simulation (which is essentially equivalent to the firstorder approximation). The origin of this tiny signal is found to be the interpolation of the Jacobian matrix between the grid points to obtain their values at the light ray positions. Sampling a Bmodefree, continuous shear field on a grid and subsequent interpolation yields again a continuous shear field. This, however, agrees exactly with the original field only at the grid points. Therefore, it may in general contain a small Bmode contribution, depending on the grid resolution and the interpolation scheme used.
4.3 Galaxygalaxy lensing
We test the effect of Born corrections and lenslens coupling on galaxygalaxy lensing (GGL) by producing a catalogue of unbiased mock galaxies. We achieve this by first drawing a number of simulation particles on each lens plane at random and using their positions as lens galaxy positions in the algorithm described in Sect. 3.4. We then obtain a catalogue of source galaxies by randomly sampling positions in the image plane assuming a uniform image distribution over the fieldofview.
The GGL signal we are interested in is given by the mean tangential shear
at the image positions of the source galaxies as a function of angular separation
to the positions of the lens galaxies. In the simple case of unbiased galaxies considered here, the expected GGL signal can be computed in the firstorder approximation by:
where J_{2} is a Bessel function of the first kind, is the probability distribution of the lens galaxies' distances, the lensing weight q(w) is given by Eq. (23), and denotes again the 3D matter power spectrum. For simplicity, we will consider a volumelimited sample of lens galaxies with constant comoving density in the following.
Due to statistical parity invariance, the cross component is expected to vanish when averaged over many sourcelens pairs. The observed mean cross component can therefore be used as a test for systematic effects and ``cosmic variance''. As shown in Fig. 13, is consistent with zero in our raytracing.
While the cross component provides a test for systematic effects, the tangential shear contains the desired information about the matter and galaxy distribution. As can be seen in Fig. 13, the mean tangential shear is significantly smaller ( at an angular separation of ) in the raytracing than expected from the firstorder prediction (28).
The reason for this discrepancy is magnification bias: lenses, i.e. dense matter structures such as galaxies or clusters with their dark matter halos, magnify the regions behind them. The magnification reduces the apparent number density of higherredshift lens galaxies around lowerredshift lenses in a volume limited survey (as has been simulated here). Underdense regions, on the other hand, demagnify the regions behind them, thereby increasing the apparent number density of lens galaxies behind them. The de/magnification leads to an anticorrelation between the positions of highredshift lens galaxies and the tangential shear induced by lowredshift structures. The anticorrelation reduces the signal compared to the firstorder approximation^{}. We can suppress the magnification bias in the raytracing by switching off the deflections and using Eq. (17) to calculate the distortions. In this case our simulations are fully consistent with the firstorder prediction, as is shown in Fig. 13.
Figure 13: Galaxygalaxylensing signal for sources at redshift z=1 and unbiased lens galaxies with a constant comoving mean densitybetween z=0 and z=1. a) Shown are the measured tangential component of the shear from full raytracing (diamonds) and raytracing using the firstorder approximation (17) (squares), and the firstorder prediction (28) (solid line). b) Measured cross component from full raytracing (diamonds) and firstorder raytracing (squares). Error bars denote the standard deviation calculated from a set of 24 simulated fields of 3 . 

Open with DEXTER 
The effect of the magnification bias on the GGL depends on the redshift distribution of the sources and the lenses. Moreover, the shape of the lens luminosity function may be important if the lens population is selected using a magnitude limit. For example, the firstorder approximation may underestimate for a lens population with a very steep luminosity function near the survey magnitude limit. We reserve a more detailed investigation of this effect with realistic source and lens distributions for future work.
5 Summary
In this work, we have described a new variant of the multiplelensplane algorithm, which is particularly suited for raytracing through very large cosmological Nbody simulations. The algorithm differs in some important details from previous works. This allows us to take full advantage of the unprecedented statistical power offered by the large volume and high spatial and mass resolution of the Millennium Simulation. The features discussed include: a tilted lineofsight (to avoid periodic repetition of structures along the lineofsight), adaptive slice boundaries (to avoid the slicing and duplication of bound structures), adaptive smoothing of the projected matter distribution on the lens planes (to reduce shot noise from the particles), a mutliplemesh method for calculating the light deflections and distortions at the lens planes (which takes into account the smallscale and largescale structure simultaneously), and a method to include galaxies (as lenses and sources) from semianalytic galaxyformation models in the raytracing process.
We have used the raytracing code and the Millennium Simulation to investigate the impact of lenslens coupling and multiple ray deflections on various cosmic shear twopoint statistics. We have computed convergence power spectra from a set of raytracing realisations. For testing and comparison, we have also computed a firstorder prediction of the convergence power spectrum using the measured threedimensional power spectra of the mass distribution in the Millennium Simulation. We find that this firstorder prediction agrees very well with the raytracing results except for very small scales (the difference is only for ), where smoothing on the lens planes becomes important.
Comparing the convergence power spectrum from the raytracing to the predictions based on the fitting formulae for the matter power spectrum by Peacock & Dodds (1996) and Smith et al. (2003), we find significant discrepancies ( for ), casting the usefulness of these fitting formulae for cosmological parameter estimation for future surveys into doubt. A prediction based on the popular halo model and the halo concentrationmass relation of Neto et al. (2007) fits better, but there are still noticeable deviations, in particular for higher source redshifts ( for sources at ). This indicates a need for more accurate descriptions of matter power spectra.
Furthermore, we have computed the E and Bmode aperture mass dispersion using our raytracing algorithm. We find the Bmode to be finite, but at least three orders of magnitude smaller than the Emode. The amplitude of the Bmode is slightly larger and shows a different scale dependence than the predictions of Cooray & Hu (2002) and Hirata & Seljak (2003). We have performed various tests to exclude numerical artifacts as the origin of the deviations. Despite these discrepancies, we can confirm the finding of Shapiro & Cooray (2006) that the lensinginduced Bmode can be safely neglected even in an allsky survey.
Corrections to the firstorder approximation can have a considerable impact on galaxygalaxy lensing. In the simple case of a volumelimited sample of unbiased lens galaxies and all sources at redshifts z=1, the firstorder approximation overestimates the mean tangential shear around lenses by at an angular separation of due to its failure to incorporate the magnification bias. The impact of the magnification bias on the galaxygalaxy lensing signal depends on the survey selection criteria and the luminosity and redshift distribution of the sources and the lenses. A detailed investigation of this effect should be carried out in future work.
Note added in proof. We became aware recently that Ziour & Hui (2008) discuss the impact of magnification bias on galaxygalaxy lensing power spectra.
Acknowledgements
We thank Volker Springel, Jeremy Blaizot, Ole Möller, Martin Kilbinger, Emilio PastorMira, and Uros Seljak for helpful discussions, Jasmin Pielorz for providing the halo model power spectra, and the referee for very useful comments. This work was supported by the DFG within the Priority Programme 1177 under the projects SCHN 342/6 and WH 6/3.
Appendix A: Lattice planes
The periodicity of our matter distribution along and perpendicular to the lineofsight can be studied within the theory of crystal lattices (see, e.g., Ashcroft & Mermin 1976). Here, we give a practical explanation rather than a rigorous proof.
Consider an array of unit cubes forming a simple cubic lattice with lattice constant unity. Choose two linearly independent lattice vectors and with and and . These two vectors span a plane which is perpendicular to the lattice vector with (q_{1},q_{2},q_{3}), .
Since the planespanning vectors are lattice vectors, the plane is itself periodic and therefore represents a plane lattice. With and as basis vectors of the plane lattice, the plane is periodic along the direction of and with periodicity length and , respectively. The parallelogram constructed from and represents a unit cell of the plane lattice with a cell area of . One can show that there is no smaller unit cell if the integer coefficients n_{1}, n_{2}, and n_{3} are coprime. Furthermore, there is no shorter nonzero lattice vector perpendicular to the plane than in this case, and hence, the shortest periodicity along the normal direction is .
For the computational cube of the Millennium Simulation with side length , the lengths and areas above have to be multiplied by L and L^{2}, respectively. Our choices and yield a LOS vector with and a rectangular area of for the lens planes.
References
 Amendola, L., Kunz, M., & Sapone, D. 2008, J. Cosmol. AstroPart. Phys., 4, 13 [NASA ADS] [CrossRef]
 Ashcroft, N. W., & Mermin, N. D. 1976, Solid State Physics (Philadelphia: Saunders College) (In the text)
 Aubert, D., Amara, A., & Metcalf, R. B. 2007, MNRAS, 376, 113 [NASA ADS] [CrossRef] (In the text)
 Bacon, D. J., Taylor, A. N., Brown, M. L., et al. 2005, MNRAS, 363, 723 [NASA ADS] [CrossRef]
 Bartelmann, M., & Weiss, A. 1994, A&A, 287, 1 [NASA ADS]
 Bartelmann, M., & Schneider, P. 2001, Phys. Rep., 340, 291 [NASA ADS] [CrossRef] (In the text)
 Benjamin, J., Heymans, C., Semboloni, E., et al. 2007, MNRAS, 381, 702 [NASA ADS] [CrossRef]
 Blandford, R., & Narayan, R. 1986, ApJ, 310, 568 [NASA ADS] [CrossRef] (In the text)
 Bower, R. G., Benson, A. J., Malbon, R., et al. 2006, MNRAS, 370, 645 [NASA ADS]
 Carbone, C., Springel, V., Baccigalupi, C., Bartelmann, M., & Matarrese, S. 2008, MNRAS, 388, 1618 [NASA ADS] [CrossRef]
 Colless, M., Dalton, G., Maddox, S., et al. 2003, VizieR Online Data Catalog, 7226, 0 (In the text)
 Cooley, J. W., & Tukey, J. W. 1965, Math. Comput., 19, 297 [CrossRef]
 Cooray, A., & Hu, W. 2002, ApJ, 574, 19 [NASA ADS] [CrossRef]
 Cooray, A., & Sheth, R. 2002, Phys. Rep., 372, 1 [NASA ADS] [CrossRef]
 Couchman, H. M. P., Barber, A. J., & Thomas, P. A. 1999, MNRAS, 308, 180 [NASA ADS] [CrossRef]
 Crittenden, R. G., Natarajan, P., Pen, U.L., & Theuns, T. 2002, ApJ, 568, 20 [NASA ADS] [CrossRef]
 Croton, D. J., Springel, V., White, S. D. M., et al. 2006, MNRAS, 365, 11 [NASA ADS] [CrossRef]
 Das, S., & Bode, P. 2008, ApJ, 682, 1 [NASA ADS] [CrossRef]
 De Lucia, G., & Blaizot, J. 2007, MNRAS, 375, 2 [NASA ADS] [CrossRef]
 Eisenstein, D. J., & Hu, W. 1999, ApJ, 511, 5 [NASA ADS] [CrossRef]
 Faure, C., Kneib, J. P., Hilbert, S., et al. 2009, ApJ, 695, 1233 [NASA ADS] [CrossRef]
 Fosalba, P., Gaztañaga, E., Castander, F. J., & Manera, M. 2008, MNRAS, 391, 435 [NASA ADS] [CrossRef]
 Frigo, M., & Johnson, S. G. 2005, Proc. IEEE, 93, 216 [CrossRef]
 Fu, L., Semboloni, E., Hoekstra, H., et al. 2008, A&A, 479, 9 [NASA ADS] [CrossRef] [EDP Sciences]
 Gavazzi, R., Treu, T., Rhodes, J. D., et al. 2007, ApJ, 667, 176 [NASA ADS] [CrossRef]
 Gottschalk, S., Lin, M. C., & Manocha, D. 1996, in SIGGRAPH '96: proceedings of the 23rd annual conference on Computer graphics and interactive techniques (New York, NY, USA: ACM Press), 171 (In the text)
 Hamana, T., & Mellier, Y. 2001, MNRAS, 327, 169 [NASA ADS] [CrossRef]
 Hartlap, J. 2005, Studying GalaxyGalaxyLensing Using RayTracing Simulations, Diplomarbeit (In the text)
 Heymans, C., White, M., Heavens, A., Vale, C., & van Waerbeke, L. 2006, MNRAS, 371, 750 [NASA ADS] [CrossRef]
 Hilbert, S., Metcalf, R. B., & White, S. D. M. 2007a, MNRAS, 382, 1494 [NASA ADS] [CrossRef]
 Hilbert, S., White, S. D. M., Hartlap, J., & Schneider, P. 2007b, MNRAS, 382, 121 [NASA ADS] [CrossRef] (In the text)
 Hilbert, S., White, S. D. M., Hartlap, J., & Schneider, P. 2008, MNRAS, 386, 1845 [NASA ADS] [CrossRef] (In the text)
 Hirata, C. M., & Seljak, U. 2003, Phys. Rev. D, 68, 083002 [NASA ADS] [CrossRef]
 Hockney, R. W., & Eastwood, J. W. 1981, Computer Simulation Using Particles (New York: McGrawHill) (In the text)
 Hoekstra, H., Mellier, Y., van Waerbeke, L., et al. 2006, ApJ, 647, 116 [NASA ADS] [CrossRef]
 Jain, B., & Seljak, U. 1997, ApJ, 484, 560 [NASA ADS] [CrossRef] (In the text)
 Jain, B., Seljak, U., & White, S. 2000, ApJ, 530, 547 [NASA ADS] [CrossRef]
 Mandelbaum, R., Seljak, U., Kauffmann, G., Hirata, C. M., & Brinkmann, J. 2006, MNRAS, 368, 715 [NASA ADS] [CrossRef]
 Massey, R., Heymans, C., Bergé, J., et al. 2007a, MNRAS, 376, 13 [NASA ADS] [CrossRef] (In the text)
 Massey, R., Rhodes, J., Ellis, R., et al. 2007b, Nature, 445, 286 [NASA ADS] [CrossRef]
 Massey, R., Rhodes, J., Leauthaud, A., et al. 2007c, ApJS, 172, 239 [NASA ADS] [CrossRef]
 Meneghetti, M., Argazzi, R., Pace, F., et al. 2007, A&A, 461, 25 [NASA ADS] [CrossRef] [EDP Sciences]
 Munshi, D., Valageas, P., van Waerbeke, L., & Heavens, A. 2008, Phys. Rep., 462, 67 [NASA ADS] [CrossRef] (In the text)
 Neto, A. F., Gao, L., Bett, P., et al. 2007, MNRAS, 381, 1450 [NASA ADS] [CrossRef] (In the text)
 Pace, F., Maturi, M., Meneghetti, M., et al. 2007, A&A, 471, 731 [NASA ADS] [CrossRef] [EDP Sciences]
 Peacock, J. A., & Dodds, S. J. 1996, MNRAS, 280, L19 [NASA ADS]
 Peirani, S., Alard, C., Pichon, C., Gavazzi, R., & Aubert, D. 2008, MNRAS, 390, 945 [NASA ADS] [CrossRef] (In the text)
 Réfrégier, A., Boulade, O., Mellier, Y., et al. 2006, in Space Telescopes and Instrumentation I: Optical, Infrared, and Millimeter, ed. J. C. Mather, H. A. MacEwen, & M. W. M. de Graauw, Proc. SPIE, 6265, 62651Y (In the text)
 Rozo, E., Nagai, D., Keeton, C., & Kravtsov, A. 2008, ApJ, 687, 22 [NASA ADS] [CrossRef]
 Schimd, C., Tereno, I., Uzan, J.P., et al. 2007, A&A, 463, 405 [NASA ADS] [CrossRef] [EDP Sciences]
 Schneider, P. 2006, in SaasFee Advanced Course 33: Gravitational Lensing: Strong, Weak and Micro, ed. G. Meylan, P. Jetzer, P. North, et al. (Berlin: Springer), 269 (In the text)
 Schneider, P., Ehlers, J., & Falco, E. E. 1992, Gravitational Lenses (Berlin, Heidelberg, New York: SpringerVerlag) (In the text)
 Schneider, P., van Waerbeke, L., Jain, B., & Kruse, G. 1998, MNRAS, 296, 873 [NASA ADS] [CrossRef]
 Schneider, P., van Waerbeke, L., & Mellier, Y. 2002, A&A, 389, 729 [NASA ADS] [CrossRef] [EDP Sciences]
 Seitz, S., & Schneider, P. 1996, A&A, 305, 383 [NASA ADS] (In the text)
 Seitz, S., Schneider, P., & Ehlers, J. 1994, Classical and Quantum Gravity, 11, 2345 [NASA ADS] [CrossRef]
 Seljak, U. 2000, MNRAS, 318, 203 [NASA ADS] [CrossRef]
 Semboloni, E., Mellier, Y., van Waerbeke, L., et al. 2006, A&A, 452, 51 [NASA ADS] [CrossRef] [EDP Sciences]
 Shapiro, C., & Cooray, A. 2006, J. Cosmol. AstroPart. Phys., 3, 7 [NASA ADS] [CrossRef]
 Simon, P., Hetterscheidt, M., Schirmer, M., et al. 2007, A&A, 461, 861 [NASA ADS] [CrossRef] [EDP Sciences]
 Smith, R. E., Peacock, J. A., Jenkins, A., et al. 2003, MNRAS, 341, 1311 [NASA ADS] [CrossRef]
 Spergel, D. N., Verde, L., Peiris, H. V., et al. 2003, ApJS, 148, 175 [NASA ADS] [CrossRef] (In the text)
 Springel, V. 2005, MNRAS, 364, 1105 [NASA ADS] [CrossRef] (In the text)
 Springel, V., White, S. D. M., Jenkins, A., et al. 2005, Nature, 435, 629 [NASA ADS] [CrossRef] (In the text)
 Taylor, A. N., Kitching, T. D., Bacon, D. J., & Heavens, A. F. 2007, MNRAS, 374, 1377 [NASA ADS] [CrossRef]
 Teyssier, R., Pires, S., Prunet, S., et al. 2008, ArXiv eprints
 Tinker, J. L., Weinberg, D. H., Zheng, Z., & Zehavi, I. 2005, ApJ, 631, 41 [NASA ADS] [CrossRef] (In the text)
 Vale, C., & White, M. 2003, ApJ, 592, 699 [NASA ADS] [CrossRef]
 Van Waerbeke, L., Hamana, T., Scoccimarro, R., Colombi, S., & Bernardeau, F. 2001, MNRAS, 322, 918 [NASA ADS] [CrossRef]
 Wambsganss, J., Cen, R., & Ostriker, J. P. 1998, ApJ, 494, 29 [NASA ADS] [CrossRef]
 Wambsganss, J., Bode, P., & Ostriker, J. P. 2004, ApJ, 606, L93 [NASA ADS] [CrossRef] (In the text)
 Wambsganss, J., Ostriker, J. P., & Bode, P. 2008, ApJ, 676, 753 [NASA ADS] [CrossRef] (In the text)
 White, M. 2005, Astropart. Phys., 23, 349 [NASA ADS] [CrossRef]
 White, M., & Hu, W. 2000, ApJ, 537, 1 [NASA ADS] [CrossRef]
 White, M., & Vale, C. 2004, Astropart. Phys., 22, 19 [NASA ADS] [CrossRef]
 Ziour, R. & Hui, L. 2008, Phys. Rev. D, 78, 123517 [NASA ADS] [CrossRef] (In the text)
Footnotes
 ... CFHTLS^{}
 http://www.cfht.hawaii.edu/Science/CFHLS
 ... KIDS^{}
 http://http://www.astrowise.org/projects/KIDS
 ... PanSTARRS^{}
 http://panstarrs.ifa.hawaii.edu
 ... LSST^{}
 http://www.lsst.org
 ... Survey^{}
 http://www.darkenergysurvey.org
 ... SNAP^{}
 http://snap.lbl.gov
 ... planes^{}
 The exact choice for the projection boundaries becomes unimportant for sufficiently small spacings between the lens planes.
 ... Simulation^{}
 This work concentrates on weak lensing, but the algorithm is also used for stronglensing studies (Faure et al. 2009; Hilbert et al. 2007b,2008).
 ... cube^{}
 We often use a larger field of view  in particular, if only lower source redshifts are considered. Even for high source redshifts, where the resulting light cone may cover the same simulation region more than once, a large field of view can be used with due care (Hilbert et al. 2007a).
 ... side^{}
 Approx. 0.5% (5%) of halos with virial masses ( ) are affected by this procedure. Though not essential for cosmic shear simulations (test show a relative difference for the shear power spectra), the proper treatment of halos near slice boundaries is important for groupgalaxy lensing and strong lensing simulations.
 ... potentials^{}
 The raytracing algorithm has to compute 5 derivatives of the lensing potential (2 deflection angle components and 3 shear matrix components) starting from the matter distribution. On large meshes, lowerorder finite differencing operations (FDs) are much faster than FFTs. Using FFT derivatives would require 6 FFTs, whereas our approach requires only 2 FFTs and 5 FDs.
 ... approximation^{}
 Note that in the firstorder approximation, magnification effects are neglected. Thus, the positions of galaxies at any given redshift are uncorrelated with the shear induced by galaxies at different redshifts.
All Figures
Figure 1: Schematic view of the observer's backward light cone in the multiplelensplane approximation. A light ray (red line) experiences a deflection only when passing through a lens plane (solid blue lines). The deflection angle of a ray passing through the lens plane at distance f_{K}^{(k1)} from the observer is obtained from the matter distribution between and projected onto the plane. Using the deflection angle of the light ray at the previous lens plane and the ray's angular positions and on the two previous planes, the angular position on the current plane can be computed. 

Open with DEXTER  
In the text 
Figure 2: Schematic view of the orientation of the lineofsight (red line) and the lens planes (blue area) relative to the simulation box (indicated by black lines). 

Open with DEXTER  
In the text 
Figure 3: Schematic view of the adaptive slice boundaries to avoid the truncation or double inclusion of halos that are located near a slice boundary. Halos near the boundary of slice k and k+1 are eitherincluded as a whole in slice k or completely excluded depending on the positions of their centres (a). Halos that are included (excluded) in slice k, are excluded (included) from slice k+1 even if they have crossed the slice boundary between redshift k and k+1 (b). 

Open with DEXTER  
In the text 
Figure 4: The number of lens planes used for the raytracing as afunction of the source redshift . 

Open with DEXTER  
In the text 
Figure 5: The deflection angle as a function of the comoving distance r to a point mass computed by our mesh method. a) Shortrange component (dashed line) and longrange component (dotted line) of the deflection angle (full line). b) Fractional difference between the value calculated by our mesh method and the theoretical value . For the plots, we placed 10 point masses of 8.6 randomly on the lens plane at z=0.5 and calculated the resulting deflection at 1000 random positions around each of them. 

Open with DEXTER  
In the text 
Figure 6: The magnification for sources at z=2 as a function of the angular separation to a point mass of at redshift z=1 computed by our mesh method. a) Numerical values (symbols) compared to the analytical result (solid line). b) Fractional difference between the measured magnification and its theoretical value . The magnification diverges at the Einstein radius (dashed vertical line). For the plots, we placed 10 point masses of randomly on the lens plane at z=1 and calculated the resulting magnification at 1000 random positions around them. 

Open with DEXTER  
In the text 
Figure 7: Interpolation scheme used for determining image positions of galaxies. The regular grid of rays in the image plane (left filled circles) is used to partition the image plane into triangles (right blue lines). The image position (left open circle) of a source inside a triangle (right blue lines) formed by the backtraced rays on the source plane (right filled circles) is then determined by linear interpolation. 

Open with DEXTER  
In the text 
Figure 8: Comparison of the true source positions and the source positions obtained by tracing back light rays through the Millennium Simulation starting from the image positions computed by the interpolation method. Shown is the difference between true and tracedback comoving source position for 1000 galaxy centres at z=1. At this redshift, comoving correspond to an angle of . The right and upper frames sides are labelled by the corresponding angular difference between true and tracedback source position in units of the mesh spacing used for the rays. 

Open with DEXTER  
In the text 
Figure 9: Convergence power spectra for sources at redshift z=1 (lower curves) and z=2 (upper curves). The simulation results (from random fields of 3 ) are shown as diamonds with errorbars (indicating standard deviation calculated from the fieldtofield variance), the corresponding firstorder predictions as solid lines. The predictions using the Peacock & Dodds (1996) prescription together with the transfer function from Eisenstein & Hu (1999) are given as dotted lines, those obtained from the Smith et al. (2003) fitting formula as dashed lines. The predictions of a halo model using the concentration parameters of Neto et al. (2007) are shown as dashdotted lines. 

Open with DEXTER  
In the text 
Figure 10: Convergence power spectra for sources at redshift z=1. Compared are the results from raytracing (symbols) using various smoothing schemes (none/Gaussian with fixed scale /adaptive) and the corresponding firstorder prediction (lines) obtained projecting and smoothing the measured 3D power spectra of the actual mass distribution in the simulation. 

Open with DEXTER  
In the text 
Figure 11: Aperture mass dispersion a) and b) as a function of filter scale for sources at z=1 and z=2. Compared are the firstorder prediction (solid lines, Emode only) and the results from raytracing (symbols with error bars indicating standard deviation, obtained from 7 fields of 5 ). 

Open with DEXTER  
In the text 
Figure 12: Bmode aperture mass dispersion for sources at z=2 decomposed into the contributions by lenslens coupling and Born corrections. Raytracing results (symbols with error bars indicating standard deviation, calculated from 7 fields of 1 ): full raytracing (squares), only Born corrections (diamonds), only lenslens coupling (triangles). Predictions by Cooray & Hu (2002): full signal (dashed line), only Born corrections (dotted line), only lenslens coupling (dashdotted line). Prediction by Hirata & Seljak (2003): full signal (solid line). 

Open with DEXTER  
In the text 
Figure 13: Galaxygalaxylensing signal for sources at redshift z=1 and unbiased lens galaxies with a constant comoving mean densitybetween z=0 and z=1. a) Shown are the measured tangential component of the shear from full raytracing (diamonds) and raytracing using the firstorder approximation (17) (squares), and the firstorder prediction (28) (solid line). b) Measured cross component from full raytracing (diamonds) and firstorder raytracing (squares). Error bars denote the standard deviation calculated from a set of 24 simulated fields of 3 . 

Open with DEXTER  
In the text 
Copyright ESO 2009