text:
# Simulating weak lensing by large scale structure

## 1. Introduction

Weak gravitational lensing is becoming an indispensable tool in modern cosmology. Building on early work by Gunn (hhm10 (1967)), Blandford et al. (hhm6 (1991)) and Miralda-Escude (hhm23 (1991)) showed that the shear and magnification are on the order of a few percent in popular cosmologies, which is now within reach of observations. Indeed detections of shear correlations by several groups (Bacon et al. hhm1 (2000); Van Waerbeke et al. 2000a ; Rhodes et al. hhm27 (2001); Hoekstra et al. hhm15 (2002); Brown et al. hhm7 (2002); Jarvis et al. hhm20 (2003)) indicate the dramatic progress being made in observational power and precision. Other projects, such as the CFHT legacy survey111http://www.cfht.hawaii.edu/Science/CFHLS, the Deep Lens Survey222http://dls.bell-labs.com/, and NOAO deep survey333http://www.noao.edu/noao/noaodeep/ and proposed projects such as SNAP444http://snap.lbl.gov, DMT/LSST555http://www.dmtelescope.org/dark\_home.html, and Pan$``$STARRS666http://www.ifa.hawaii.edu/pan-starrs/ will continue this trend as they map the shear on large fractions of the sky. In order for the full scientific return from lensing to be realized, both fundamental theory and data analysis techniques must keep pace with these observational advances.

In principle the effect exploited by these surveys is simple: gravitational lensing of light rays by foreground large-scale structure magnifies and shears the images of background galaxies, and this distortion can be mapped to provide information on the matter distribution and cosmological model. Like the cosmic microwave background (CMB), the theory of weak lensing is simple and clean (see Bartelmann & Schneider hhm4 (2001) or Mellier hhm22 (1999) for a review). Unlike the the CMB, however, lensing on sub-degree scales probes large-scale structure in the non-linear regime and a full description has not been achieved through purely analytic means. Fortunately, on scales relevant to weak lensing, the hierarchical growth of structure can be accurately simulated using N-body methods. The lensing effect of this structure can then be computed directly by ray tracing through the simulation. In this paper we describe an implementation of such a method, building upon the work of (Blandford et al. hhm6 (1991); Wambsganss, Cen & Ostriker hhm36 (1998); Couchman, Barber, & Thomas hhm8 (1998); Fluke, Webster, & Mortlock hhm9 (1998); Hamana, Martel & Futamase hhm11 (2000); Jain, Seljak, & White hhm19 (2000); White & Hu hhm38 (2000); Barber, Thomas, Couchman, & Fluke hhm2 (2000); Hamana & Mellier hhm12 (2001); Hennawi et al. hhm14 (2001); Padmanabhan, Seljak, & Pen hhm25 (2002)). Our goal is to quantify the numerical precision of such algorithms, and to investigate the computational cost of achieving a given theoretical precision.

The outline of the paper is as follows. In §2 we introduce the weak lensing formalism and in §3 we describe our implementation of a multi-plane ray-tracing algorithm. In §4 we present some basic results for comparison to previous work while in §5 we evaluate the numerical convergence of the algorithm. Some interesting physical results are given in §6 before we conclude in §7.

## 2. Weak lensing formalism

We make use of a multiple lens plane algorithm in order to simulate the distortion and magnification effect of foreground matter on background light rays. In the following section, we provide a summary of equations directly relevant to this paper (see Jain, Seljak, & White (hhm19 (2000)) or Schneider, Ehlers, & Falco (hhm28 (1992)) for a more thorough treatment) and then describe our ray tracing method. Note that we adopt units where $`c=1`$, and we will work in comoving coordinates.

In a universe governed by the Robertson-Walker metric, the change in direction of a light ray propagating through space is:

$`d\stackrel{}{\alpha }=2_{}\varphi d\chi `$ (1)
where $`d\stackrel{}{\alpha }`$ is the bend angle, $`_{}`$ denotes the spatial gradient perpendicular to the path of the light ray, $`\varphi `$ is the 3 dimensional peculiar gravitational potential, and $`\chi `$ is the radial comoving coordinate.

The change in position on a plane perpendicular to the line of sight at a position $`\chi `$ due to a deflection $`d\stackrel{}{\alpha }`$ at $`\chi ^{}`$ is:

$`d\stackrel{}{x}(\chi )=r(\chi \chi ^{})d\stackrel{}{\alpha }(\chi ^{})`$ (2)
where $`r(\chi )`$ is the comoving angular diameter distance. Integrating to get $`\stackrel{}{x}(\chi )`$ and then dividing by $`r(\chi )`$ yields:

$$\stackrel{}{\theta }(\chi )=\frac{2}{r(\chi )}_0^\chi 𝑑\chi ^{}r(\chi \chi ^{})_{}\varphi +\stackrel{}{\theta }(0)$$
(3)
where $`\stackrel{}{\theta }(\chi )`$ is the angular position of the light ray at comoving coordinate $`\chi `$. The $`2\times 2`$ distortion matrix, $`A_{ij}\theta _i(\chi )/\theta _j(0)`$, is given by:

$`A_{ij}=2{\displaystyle _0^\chi }𝑑\chi ^{}g(\chi ^{},\chi )_i_j\varphi +\delta _{ij}`$ (4)
where we make use of the definition:

$`g(\chi ^{},\chi ){\displaystyle \frac{r(\chi \chi ^{})r(\chi ^{})}{r(\chi )}}`$ (5)
This matrix is normally decomposed as:

$`𝐀=\left(\begin{array}{cc}1\kappa \gamma _1& \gamma _2\omega \\ \gamma _2+\omega & 1\kappa +\gamma _1\end{array}\right)`$ (6)
where $`\kappa `$ is the convergence, $`\gamma `$ is the shear, and $`\omega `$ is the rotation.

The gravitational potential, $`\varphi ^{}`$, is related to the mass density, $`\rho `$, through Poisson’s equation:

$$^2\varphi ^{}=4\pi Ga^2\rho $$
(7)
where the gradient is with respect to comoving coordinates (that’s why there’s an $`a^2`$ on the right hand side) and $`\varphi ^{}`$ is the sum of a smooth background potential $`\overline{\varphi }`$ and the local peculiar potential $`\varphi `$. We can subtract the background terms and use $`^2\overline{\varphi }=4\pi G\overline{\rho }`$, where $`\overline{\rho }`$ is the mean density in the region, and substitute $`a^3\overline{\rho }=\overline{\rho }_0`$ to find

$$^2\varphi =\frac{4\pi G\overline{\rho }_0}{a}\left(\frac{\rho }{\overline{\rho }}1\right).$$
(8)
If we then define the relative mass overdensity, $`\delta \rho /\overline{\rho }1`$, and substitute $`\overline{\rho }_0=3H_0^2\mathrm{\Omega }_m/8\pi G`$ we obtain

$$^2\varphi =\frac{3H_0^2\mathrm{\Omega }_m\delta }{2a}.$$
(9)
In the discrete approximation, the interval between the source and the observer is divided into distinct regions separated by a comoving distance $`\mathrm{\Delta }\chi `$. The relative mass overdensity in a given region is then projected onto a lens plane within the region and perpendicular to the line of sight. An effective two dimensional gravitational potential, $`\psi =2\varphi 𝑑\chi `$, is defined on the this plane, enabling a two dimensional Poisson equation relating $`\psi `$ and the matter content in the region to be written as

$$^2\psi =3H_0^2\mathrm{\Omega }_m\mathrm{\Sigma }$$
(10)
where $`\mathrm{\Sigma }=\delta 𝑑\chi `$ is the projected two dimensional relative mass overdensity. Eqs (3) and (4) then become

$`\stackrel{}{\theta }_n={\displaystyle \underset{m=1}{\overset{n1}{}}}{\displaystyle \frac{r(\chi _n\chi _m)}{r(\chi _n)}}_{}\psi _m+\stackrel{}{\theta }_1`$
$`𝐀_n=𝐈{\displaystyle \underset{m=1}{\overset{n1}{}}}g(\chi _m,\chi _n)𝐔_m𝐀_m`$ (11)
where $`𝐔_m`$ is the shear tensor in the $`m^{th}`$ region defined by

$$U_{ij}\frac{^2\psi _m}{x_ix_j}$$
(12)
It will be useful to decompose the distortion matrix at a given plane into two components (Seitz, Schneider, and Ehlers hhm29 (1994)) $`𝐁_n`$ and $`𝐂_n`$ such that for a zero curvature universe (assumed hereafter):

$`𝐀_{n+1}𝐈𝐁_n+𝐂_n`$
$`𝐁_n=𝐁_{n1}+\chi _n𝐔_n𝐀_n`$
$`𝐂_n={\displaystyle \frac{\chi _n}{\chi _{n+1}}}(𝐂_{n1}+\chi _n𝐔_n𝐀_n)`$ (13)
## 3. ray-tracing algorithm

The goal of our weak lensing algorithm is to simulate the deflection and distortion that light rays would experience as they propagate through an intervening matter distribution that is statistically similar to that of the real universe. In this section we describe two different implementations of the weak lensing algorithm: first, a standard implementation based on multiple projected-mass lens planes, and then a second, three-dimensional version. The modeling of the mass distribution is discussed in §5.

### 3.1. Standard (2D) algorithm

The standard implementation consists of creating a two-dimensional projected mass plane from an N-body simulation, computing the position, propagation direction, and distortion matrix on the plane for each ray, and then repeating the process until the rays have been traced from the observer to the source. Here is a more detailed description:

1. We begin by using an N-body simulation (which we describe later) to create the structure responsible for the lensing. Since lensing is only sensitive to density fluctuations perpendicular to the line of sight, the simulation must resolve structures on length scales $`\chi \theta `$, where $`\chi `$ is a typical survey depth, and $`\theta `$ is the desired angular resolution. Light-rays from redshifts of interest must therefore travel thousands of Mpc through structures resolved below Mpc scales.

It is currently computationally very expensive to simulate both scales simultaneously, but, fortunately, this is not necessary. Since there is little lensing power on very large scales, it makes no sense to simulate huge boxes whose pieces are almost independent. We employ a technique (common in the literature) in which we simulate a volume which is large enough to be a fair sample of the universe, yet small enough to provide sufficient small scale resolution at reasonable cost. In this simulation, mass particles in a box evolve in time under the influence of gravity, and the positions of the particles are recorded at time intervals corresponding to regular intervals in comoving distance. These output boxes are then used to create the mass distribution to be traced through. Note that since each box is actually the same matter distribution at different stages of evolution, care must be taken to avoid tracing over the same structures more than once.

2. A selected number of light rays to be traced are then initialized at the observer for a square field of view. For each ray, we track the position $`\stackrel{}{x}`$, the propagation direction $`\stackrel{}{\alpha }_n`$ relative to the line of sight, and the matrices $`𝐁`$ and $`𝐂`$ from equation (13). The recursion relations (Eq. 11) for $`𝐀_n`$ and $`\stackrel{}{\theta }_n`$ require as inputs the previous values for $`_{}\psi `$ and $`𝐔`$, all of which need to be stored. However, we make use of a less memory intensive algorithm for the distortion matrix based on Eq. (13) in which we track the values the decomposed distortion matrix $`𝐁`$ and $`𝐂`$ for each ray. Each iteration requires as a new input only the current value of $`𝐔`$ at each ray location, which may then be discarded. Similarly, by tracking $`\stackrel{}{\alpha }_n`$, we require only one value of $`_{}\psi `$ at a time.

3. Next, we select the source redshift and the number of lens planes we wish to use. Note that the sources are all initialized at the same redshift, and the planes are equally spaced in comoving distance. This sets the comoving distance between planes, $`\mathrm{\Delta }\chi `$. We create a lens plane from an N-body simulation box using the following technique. We select the N-body particle dump closest in comoving distance to that of the desired lens plane, and then randomly choose the $`x`$, $`y`$, or $`z`$ axis as the propagation axis, $`\chi `$. A two-dimensional $`\mathrm{N}\times \mathrm{N}`$ lens plane grid is then created by projecting mass particles from the N-body simulation in a direction parallel to the $`\chi `$ axis onto a plane at the desired comoving distance, and then recording the mass onto the grid locations using a cloud in cell (CIC) assignment scheme. All of the particles in the box within $`\frac{1}{2}\mathrm{\Delta }\chi `$ of the desired plane are used to make the grid, which is then normalized by the mean density. The projected field is thus a square of side length equal to that of the box.

The origin of the lens plane is randomly selected. This origin shift is facilitated by the fact that the N-body simulation uses periodic boundary conditions, and we can use this periodicity to re-map the mass density about any origin we choose. This process of rotation and origin shift makes it unlikely that the same structures will be traced repeatedly. Note that both the origin and orientation of the axes are maintained until the box is traced all the way through, at which point the box is re-randomized, so that consecutive lens planes are made from a matter distribution that is continuous everywhere except at the box boundaries.

A variation of this technique, to our knowledge presented here for the first time, and which we make use of for comparison, avoids introducing these discontinuities entirely. The propagation direction is selected at an angle relative to one of the box axes. The rays are initialized about a central line of sight along this direction, and the lens planes are created perpendicular to the central line of sight and with their origins on it. The angle is chosen so that the rays, traced from the observer to the source, never propagate through the same volume twice. However, since the boxes are never rotated, nor is there any random origin shift, the lens planes are made from a matter distribution that is continuous everywhere.

4. The two component gradient of the potential $`\psi `$ and the four components of the matrix $`𝐔`$ are then computed on the grid. This is done by first taking the Fourier transformation of the mass density plane, then taking the first or second derivatives, respectively, and computing the inverse transformation. We do this by making use of the relation of real space derivatives in Fourier space, for example

$`{\displaystyle \frac{\psi (x)}{x}}=i{\displaystyle k_x\stackrel{~}{\psi }(k)e^{i\stackrel{}{k}\stackrel{}{x}}d^2k}`$ (14)
5. The position $`\stackrel{}{x}_n`$ for each ray on the plane is computed using the stored values $`\stackrel{}{x}_{n1}`$ and $`\stackrel{}{\alpha }_{n1}`$. Values for $`𝐔`$ and $`\psi `$ are then assigned at each ray position, again using a CIC assignment scheme, and the new deflection angle $`\stackrel{}{\alpha }_n=\stackrel{}{\alpha }_{n1}+\mathrm{\Delta }\stackrel{}{\alpha }_n`$ is then computed where

$`\mathrm{\Delta }\stackrel{}{\alpha }=\psi `$ (15)
follows from Eq. (1). The distortion matrix $`𝐀_n`$ is temporarily created from stored values of $`𝐁_{n1}`$ and $`𝐂_{n1}`$ and used along with $`𝐔_n`$ to create the new values $`𝐁_n`$ and $`𝐂_n`$ for each ray on the plane.

6. The process is repeated until each ray is traced all the way back to its source, thus ensuring that all the rays will meet at the observer, and the final angular position and distortion matrix are recorded for each ray.

Clearly, one limitation of any algorithm based on lens planes is that the light rays themselves do not all propagate in precisely the same direction and cannot all be perpendicular to a given lens plane. The weak lensing formalism of the previous section calls for the calculation of derivatives perpendicular to the light rays propagation direction, so the multiple lens plane algorithm is necessarily only approximate in this respect.

### 3.2. Extended (3D) algorithm

We have created a second version of the ray-tracing algorithm that is fully three-dimensional, and therefore doesn’t suffer from this drawback. It is quite similar to the standard implementation in most respects, and makes use of the same N-body particle dumps. Here is a brief description.

1. The source redshift is chosen and light rays in a square field of view are initialized at the observer. As before, we will track the position and propagation direction of each ray, and also the matrices $`𝐁`$ and $`𝐂`$. Note that these are still $`2\times 2`$ matrices, since, like the matrix $`𝐀`$, they still describe the distortion of a two-dimensional image.

2. The matter distribution from the particle dump nearest in redshift to the current position is mapped onto a three-dimensional regular grid using a CIC assignment scheme, and the grid is then normalized by the mean density. As before, it’s important not to re-trace the same structures, so the box orientations and origins are randomly reassigned, when necessary, to avoid this.

3. The relationship in Fourier space between $`\varphi `$ and $`\delta `$ is used to calculate the three first derivatives and six independent second derivatives of $`\varphi `$ at each point on the grid. These derivatives are then evaluated using the CIC method at each ray position. The components perpendicular to each ray direction are calculated at each ray position and used to update the values of $`𝐁`$, $`𝐂`$ and $`\stackrel{}{\alpha }_n`$. The rays are then projected forward by a pre-selected comoving distance, and the derivatives are then evaluated at the new positions. This process is repeated until it becomes time to use the next particle dump.

## 4. basic results

Before we discuss the numerical convergence of our technique, we present some basic science results. These are broadly in agreement with similar studies presented by (Jain, Seljak, & White hhm19 (2000); Hamana, Martel, & Futamase hhm11 (2000); White & Hu hhm38 (2000)) as we elucidate later. Fig. 1 is an example of the convergence ($`\kappa `$) and shear ($`\gamma `$) maps that result from the ray-tracing procedure. The maps are $`2^{}`$ on a side and contain $`2048^2`$ lines of sight. The average magnitude of the convergence and shear is about $`1.8\%`$, close to the level predicted by Blandford et al. ( hhm6 (1991)) and Miralda-Escude ( hhm23 (1991)).

In the following section we describe these fields. We begin by presenting the 1-point distribution, and demonstrate that there are not likely to be any unlensed images at high redshift for this cosmology. We then present the power spectrum, which we find to be consistent with the Born and Limber approximations. We then discuss the skewness ($`S_3`$) and the kurtosis ($`S_4`$), which are given as a function of angular smoothing scale. We conclude with a description of the convergence peaks in our maps.

The lowest order statistic for the maps is the distribution of the shear amplitude or convergence. We show two typical examples in Fig. 2 in terms of the magnification, which in the weak lensing limit is simply $`\mu =1+2\kappa `$. This is compared to the ‘universal’ form of Wang, Holz, & Munshi (hhm37 (2002)) in §6. The average magnitude of the convergence over 10 realizations is 0.018, consistent with the expectation of a weak lensing effect on the order of a few percent. The mean and variance of the convergence are -0.002 and $`7.5\times 10^4`$, respectively. The minimum $`\kappa `$ for any of our realizations was $`0.049`$, which is greater than the ‘empty beam’ value (an empty beam is a path from the source to the observer which is sufficiently void of matter that gravitational interactions may be ignored) of $`0.064`$ for our cosmology, consistent with the fact that all the light rays in our simulations are lensed to some extent.

The shear correlations, from which cosmological inferences are usually drawn, are shown in Fig. 3 in terms of angular power spectra. Figure 3 shows the power spectrum of ($`\kappa `$) vs. the spectra of the E-mode of the shear ($`\gamma _E`$) and of the rotation ($`\omega `$) for one of our simulated maps. The shear E-mode is defined to be the curl free component of the shear tensor, and is easily calculated in Fourier space using

$$\stackrel{~}{\gamma }_E(\stackrel{}{k})=\frac{(k_x^2k_y^2)\stackrel{~}{\gamma }_1(\stackrel{}{k})+2k_xk_y\stackrel{~}{\gamma }_2(\stackrel{}{k})}{k_x^2+k_y^2}$$
(16)
where the tilde denotes the Fourier transform. The ‘excess’ power in $`\kappa `$ and $`\gamma _E`$ at $`\mathrm{}200`$ is just a fluctuation in this particular map.

One of the most important predictions of the weak lensing approximation is that the shear matrix $`𝐀`$ can be well approximated using the 2nd derivatives of a scalar potential, as would be the case if only a single lens plane were used to compute it. Then $`𝐀`$ would be symmetric, $`\omega `$ would be zero, and $`\gamma _E`$ would equal $`\kappa `$. Thus, in the weak lensing regime, the spectra of $`\kappa `$ and $`\gamma _E`$ are predicted to be nearly equal to each other, and to be much larger than that of $`\omega `$. We agree with earlier work (e.g. JSW) in verifying this basic assumption of weak lensing (Fig. 3).

Figure 4 shows the angular power spectrum of $`\kappa `$ (at various levels angular resolution) compared to the semi-analytic result which makes both the Born and Limber approximations. This result was first derived by Kaiser (Kaiser (1992)) and extended by Jain & Seljak (hhm18 (1997)). We present it here in the form of White & Hu (hhm38 (2000))

$$\mathrm{\Delta }_\kappa ^2(l)=\frac{9\pi }{4l}\mathrm{\Omega }_m^2H_0^4\chi ^{}𝑑\chi ^{}\left[\frac{g(\chi ^{},\chi )}{a(\chi ^{})}\right]^2\mathrm{\Delta }_\mathrm{m}^2(k,a)$$
(17)
where $`g(\chi ^{},\chi )`$ is defined in Eq. (5), $`\mathrm{\Delta }_\mathrm{m}^2(k)=k^3P(k)/(2\pi ^2)`$ is the contribution to the mass variance per logarithmic interval in wavenumber and $`\mathrm{\Delta }_\kappa ^2(\mathrm{})=\mathrm{}^2C_{\mathrm{}}/(2\pi )`$ is the contribution to $`\kappa _{\mathrm{rms}}^2`$ per logarithmic interval in angular wavenumber $`\mathrm{}`$.

The agreement between the numerical and semi-analytic predictions is quite good, but not perfect. One reason for this could be that in evaluating Eq. (17) we used the method of Peacock & Dodds (hhm26 (1996)) to compute the non-linear power spectrum $`P(k)`$. As can be seen in Figure 5, where we compare this prediction with the measured 3-dimensional matter power spectrum of our N-body simulation, the Peacock & Dodds hhm26 (1996)) fitting formula is not exact. If we replace the $`\mathrm{\Delta }_\mathrm{m}^2(k)`$ in Eq. (17) with the values actually measured in our simulation, the discrepancy disappears. We conclude therefore that the Born and Limber approximations are consistent with the power spectrum results from our simulations (we shall return to this issue in §6).

The maps in Fig. 1 are noticeably non-Gaussian, with identifiable objects accounting for large positive values of the convergence. The 2-point statistics, such as the power spectra, therefore cannot contain all of the information in the map. Further cosmological information is contained in the higher order moments, however the number of degrees of freedom rapidly becomes large making it difficult to present all of the information. For this reason we shall present only the ‘collapsed’ N-point functions, and since the non-Gaussianity is most clearly manifest in real (rather than Fourier) space we shall present only real-space results. In Figure 6 we show the skewness and kurtosis of the convergence on an angular scale $`\theta `$. These are obtained by smoothing the $`\kappa `$ maps using the IDL command SMOOTH, which performs a real space boxcar average using square apertures with an integer number of pixels on a side. We then calculate the skewness

$`S_3(\theta )={\displaystyle \frac{\kappa _\theta ^3}{\kappa _\theta ^2^2}}`$ (18)
and the kurtosis

$`S_4(\theta )={\displaystyle \frac{\kappa _\theta ^43\kappa _\theta ^2^2}{\kappa _\theta ^2^3}}`$ (19)
by averaging over the pixels in a given field and then averaging the results for ten realizations. Individual realizations vary substantially, suggesting error bars of $`5\%`$ for $`S_3`$ and $`10\%`$ for $`S_4`$, in line with the results of White & Hu (hhm38 (2000)). Figure 7 shows the dependence on scale of the lower order moments for two cosmological models which differ only in the clustering strength, $`\sigma _8`$.

An accurate measurement of the higher order moments, $`S_3`$ and $`S_4`$, from lensing maps may provide important constraints on the cosmological parameters $`\mathrm{\Omega }_\mathrm{m}`$ and $`\sigma _8`$ (Bernardeau, van Waerbeke, & Mellier hhm5 (1997); Jain & Seljak hhm18 (1997)). In perturbation theory, the variance of the convergence exhibits a strong dependence on both quantities, while the skewness and kurtosis are independent of the latter. Although the independence of $`S_3`$ and $`S_4`$ from $`\sigma _8`$ is not expected to hold on scales of interest where the lensing signal can be measured, a joint measurement of these quantities should still be useful for breaking the degeneracy between $`\mathrm{\Omega }_\mathrm{m}`$ and $`\sigma _8`$ if the actual dependence can be understood. This has motivated efforts to extend skewness and kurtosis predictions to smaller angular scales.

We compare some of these calculations with a variety of numerical results to assess the degree of convergence between groups and methods. Table 1 presents $`S_3`$ and $`S_4`$ reported by several groups, both from simulations and by analytic means. All are for cosmological models with $`\mathrm{\Omega }_\mathrm{m}=0.3=1\mathrm{\Omega }_\mathrm{\Lambda }`$ and $`h=0.7`$. The moments from simulations are from JSW, HM (Hamana & Mellier hhm12 (2001)), WH (White & Hu hhm38 (2000)), H2 (Hamana hhm13 (2003)), and this paper (VW), while the last two rows are analytic predictions. The first of these, v.W (van Waerbeke 2000b ), is for a semi-analytic fit to N-body simulations which tries to interpolate between the Hyper-Extended Perturbation Theory (HEPT) prediction (Scoccimarro & Frieman hhm31 (1999); Hui hhm17 (1999)) and the perturbation theory result. The second, TJ (Takada & Jain hhm33 (2002)), uses the Halo model approach.

We note that the skewness in the highly non-linear regime is predicted by HEPT to vary as $`\sigma _8^{0.4}`$. Lowest order perturbation theory predicts that $`S_3`$ and $`S_4`$ would be independent of $`\sigma _8`$. If we adjust all of the results to $`\sigma _8=0.9`$ the disagreements are lessened (for example using HEPT $`S_3(1^{})`$ would scale to 144 for VWa and to 148 for VWb). Keeping in mind that the simulation results have uncertainties of 10% or so, the agreement between groups is quite good, though more work needs to be done if we are to obtain precision cosmology from the higher order moments.

Finally it is interesting to look at the extrema of the lensing maps. We located the convergence peaks for our maps by selecting all pixels with $`\kappa `$ greater a selected threshold convergence $`\kappa _{\mathrm{min}}`$, and then directly comparing the convergence of nearby pixels within a radius of 0.3 arcminutes. The largest value was then selected as the peak’s center. Table 2 shows the average number of peaks in a $`1^{}\times 1^{}`$ field for selected values of $`\kappa _{\mathrm{min}}`$. The average profile of these peaks is roughly consistent with that of a projected NFW profile (Bartelmann hhm3 (1996); Navarro, Frenk, & White hhm24 (1997)).

## 5. Numerical issues

The numerical convergence of the algorithm described in the previous section is limited by the size and resolution of the N-body simulation and by the finite grids used to compute the distortion and deflection of the light rays. In the following section, we make use of the power spectrum as a metric in order to test these limits under various conditions, and to evaluate the basic ray-tracing implementation against the less approximate versions. We then comment on the numerical issues involving the skewness and kurtosis.

### 5.1. N-body Simulations

We use N-body simulations to calculate the evolution of the dark matter which governs the formation of large-scale structure. We use 3 simulations each of a $`\mathrm{\Lambda }`$CDM model. The first simulation evolves $`512^3`$ particles in a $`300h^1`$Mpc box, using the TreePM-SPH code (see the appendix of White hhm39 (2002)) but without any SPH particles. The cosmological model is chosen to provide a reasonable fit to a wide range of observations with $`\mathrm{\Omega }_\mathrm{m}=0.3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$, $`H_0=100h\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$ with $`h=0.7`$, $`\mathrm{\Omega }_Bh^2=0.02`$, $`n=1`$ and $`\sigma _8=1`$ (corresponding to $`\delta _H=5.3\times 10^5`$). The simulation was started at $`z=60`$ and evolved to the present with the full phase space distribution dumped every $`100h^1`$Mpc from $`z2`$ to $`z=0`$. The gravitational force softening was of a spline form, with a “Plummer-equivalent” softening length of $`20h^1`$kpc comoving, and the particle mass is $`1.7\times 10^{10}h^1M_{}`$. The 3-d mass power spectrum of this simulation is is given for selected redshifts in Figure 5.

We use this simulation as our ‘fiducial’ model unless otherwise specified. To test the dependence on the amplitude of the power spectrum we also used a second simulation which is identical to the first except that the amplitude of clustering is reduced: $`\sigma _8=0.8`$. We also make use of a smaller, lower resolution simulation for certain of the numerical tests described below. This final simulation is similar to the larger ones except that it uses fewer ($`128^3`$) particles and a smaller ($`128h^1`$Mpc) box and it dumps the output more frequently (every $`25h^1`$Mpc). The softening length in this simulation is $`36h^1`$Mpc.

### 5.2. Large scale resolution

Our ability to resolve large scale features is limited by the field of view of our ray tracing, and therefore by the size of the box of the N-body simulation and the distance to the furthest lens plane. At $`z=1`$, this limit is 7.5 degrees for our larger simulation and 3.2 degrees for our smaller one. However, most of our ray tracing runs were made at a field of view of 2 degrees. Only a few modes are available to contribute to the power spectrum on large scales for any given realization, leading to substantial variation for different realizations. We therefore generate multiple realizations of ray tracing for an N-body simulation and average them to obtain estimates on large scales. For each realization, we project mass randomly along one of three different axes, and then select the origin of the projected plane at random. Depending on the distance from the observer, only a fraction of this projected plane is used in any given ray tracing realization. For example, a 2 degree run on our larger simulation uses less than 2% of a lens plane at the peak of the lensing kernel. Multiple nearly independent realizations are therefore possible using this method, and the results can be averaged to obtain lensing statistics.

### 5.3. Small scale resolution

The discrete nature of our ray tracing procedure introduces numerical smoothing when we interpolate mass density values from the N-body simulation onto the Fourier grid. The N-body simulation is characterized by both a finite spatial resolution and by shot noise due to the finite particle number. To investigate the effect of this limitation on the angular resolution of the power spectrum, we model the finite resolution effect with Gaussian filters and thereby modify Eq. (17) so that $`\mathrm{\Delta }_\mathrm{m}^2(k)`$ is replaced with

$$\mathrm{\Delta }_e^2(k)\left(\mathrm{\Delta }_\mathrm{m}^2(k)e^{\sigma _n^2k^2}+\frac{k^3}{2\pi ^2\overline{n}}\right)e^{\sigma _g^2k^2}$$
(20)
where $`\mathrm{\Delta }_e^2(k)`$ is an effective 3-d mass power spectrum, $`k^3/2\pi ^2\overline{n}`$ is the analytic expression for the shot noise, $`\overline{n}`$ is the mean particle number density, and $`\sigma _n`$ and $`\sigma _g`$ are characteristic resolution limits of the N-body simulation and the Fourier grid, respectively. We did not find it necessary to introduce a term for shot-noise on each individual lens plane.

We generate spectra using Eq. (20), normalizing the height to equal that of the simulation at $`\mathrm{}=300`$. As shown in Fig. 8, the spectra generated in this manner using $`\sigma _g=0.54L_{\mathrm{box}}/N_{\mathrm{grid}}`$ (corresponding to a full width half max resolution of 1.3 grid cells) and $`\sigma _n`$ set equal to zero are a good fit with those from our high resolution simulation. This suggests that for this simulation the power spectrum resolution is limited by the grid and not by the N-body resolution. It is also consistent with our expectation that the Fourier grid’s resolution is roughly equal to its spacing.

The limits imposed by the N-body simulation can be seen in the spectra of our low resolution simulation (Fig. 9), which are well approximated by those generated using Eq. (20), with $`\sigma _g=0.54L_{\mathrm{box}}/N_{\mathrm{grid}}`$ as before, and $`\sigma _n=50h^1\mathrm{kpc}`$, or 5% of the mean inter-particle spacing $`\overline{n}^{\frac{1}{3}}`$. The three curves are for different values of N. The curve for $`N=1024`$ is similar to those of Figure 8, however the roll over is now caused by both the grid resolution and the N-body resolution. The curves at $`N=2048`$ and $`N=4096`$ show substantially more power at high angular resolution. This is because the noise term, which is larger than the signal from this N-body simulation for $`\mathrm{}10^4`$ and above, is no longer suppressed by the grid resolution.

The value for $`\sigma _n`$ is independent of grid spacing and is instead related to the resolution of the simulation. However, it is not equal to the “Plummer equivalent” softening length, which is $`36h^1`$kpc. To check this, we computed the 3-d mass power spectrum $`\mathrm{\Delta }_\mathrm{m}^2(k)`$ using the rebinning technique of Smith et al. (hhm32 (1999)) for three simulations with identical initial conditions but varying particle number. We find that $`\sigma _n=0.05\overline{n}^{\frac{1}{3}}`$ does roughly reproduce the roll off in power between the models, though it clearly doesn’t get the whole form right because the models have more power at high-$`k`$ than this filter suggests.

We make use of our low resolution simulation to investigate the convergence of the algorithm with respect to the number of lens planes used and to determine the effect of decreasing the redshift/time interval between particle dumps of the N-body simulation. Figure 10 shows the power spectrum, averaged over 10 realizations, for different values of inter-plane spacing ($`\mathrm{\Delta }\chi `$) and redshift/time interval ($`\mathrm{\Delta }t`$). A change in $`\mathrm{\Delta }\chi `$ from $`125h^1`$Mpc to $`25h^1`$Mpc produced a change in the power spectrum roughly at the $`10\%`$ level for this simulation (for our high resolution simulation this value was $`5\%`$), however, additional decreases in $`\mathrm{\Delta }\chi `$ did not significantly improve the resolution. For a similar change in $`\mathrm{\Delta }t`$, the spectra are within $`1\%`$ of each other. This is not surprising, since structures evolve on a scale roughly equal to the Hubble time, which is much longer than the time it takes for light to travel $`100h^1`$Mpc.

We also wish to understand the numerical resolution of the weak lensing algorithm as it relates to the higher order moments. Figure 6 shows the effect on $`S_3`$ and $`S_4`$ of changing the Fourier grid resolution and the inter-plane spacing $`\mathrm{\Delta }\chi `$. Increasing the resolution of the grid increases the small angle power in both the skewness and kurtosis, while decreasing $`\mathrm{\Delta }\chi `$ increases their power over the full range of resolution. However, this effect seems to be due entirely to the decrease in the convergence variance $`\kappa ^2`$, while $`\kappa ^3`$ and $`\kappa ^4`$ are essentially unaffected.

### 5.4. Discontinuities

The process of rotation and origin shifting for each box necessarily leads to discontinuities in the 3-dimensional mass distributions, which may, in principle, introduce artifacts during ray-tracing runs. To investigate this possibility we traced through our large simulation using the “angle” method described above. We randomly select the origin and orientation of the first box and then trace the rays through all the stacked boxes without any further rotations or origin shifts. The rays are then initialized in a field of view about a central line of sight that is angled with respect to the boxes to avoid repeating any structures, and the lens planes are made so that they are perpendicular this line of sight. We averaged the power spectra of ten runs created using this method and find that it is the same as for the standard “rotate and shift” method to better than $`1\%`$.

### 5.5. 3-dimensional ray tracing

One approximation inherent in the lens plane method is that since the light rays are all traveling in different directions they are only roughly perpendicular to the lens planes. The weak lensing formalism is based on Eq. (1), $`d\stackrel{}{\alpha }=2_{}\varphi d\chi `$, which requires that we use the gradient perpendicular to the path of the light-rays.

We used the 3-dimensional ray-tracing algorithm described previously to test the validity of this approximation, and we find that it holds extremely well for two degree fields. We used three-dimensional Fourier grids with $`512\times 512\times 32`$ points, with the latter number applying to the direction of the line of sight. We find that convergence maps made with the two different methods agree to better than $`0.1\%`$, and measures such as the power spectra, skewness, and kurtosis are virtually indistinguishable.

## 6. Physical Issues

### 6.1. The Born approximation near peaks in the convergence

An approximation that is often useful in weak lensing analysis is the Born approximation, in which the lensing events for a given light ray are completely decoupled from one another. The distortion and deflection contributions for each event are computed at points along an undeflected path and as if for an undistorted ray. These are then summed to produce the total deflection and distortion. The Born approximation can be calculated by tracing rays along straight paths without deflections and approximating $`𝐀_n`$ from Eq. (11) as

$`𝐀_n=𝐈{\displaystyle \underset{m=1}{\overset{n1}{}}}g(\chi _m,\chi _n)𝐔_m`$ (21)
The Born approximation has been shown before to be valid for whole convergence maps using measures such as the power spectrum. However, a breakdown of this approximation, if one were to occur, would almost certainly take place where lensing effects are strongest, which could in principle be an issue of relevance to lensing by clusters of galaxies in the real universe. In Fig. 11 we compare 300 peaks with convergence values of at least $`0.35`$ computed using our standard ray-tracing vs. the same peaks using the Born approximation.

We find that the Born approximation is nearly always accurate to within 10% even for these large peaks, and the single biggest departure was 20%. In addition, the mean of the ratios is very nearly equal to one, so deviations appear to be unbiased. We expect that the Born approximation becomes increasingly accurate as $`\kappa `$ is decreased, but matching the ‘Born peaks’ with those of the full ray tracing (including deflections) becomes increasingly difficult so we have not been able to quantify this.

### 6.2. Comparison to a universal magnification PDF

In Figure 2, we compare the magnification probability distribution function of a typical run, $`P(\mu )`$, with the analytical prediction using the universal probability distribution function $`P(\eta )`$ (Wang, Holz, & Munshi hhm37 (2002)), where

$`P(\mu )=P(\eta )/2|\kappa _{\mathrm{min}}|.`$ (22)
The reduced convergence is defined as $`\eta \frac{1}{2}(\mu 1)/|\kappa _{\mathrm{min}}|+1`$, and $`\kappa _{\mathrm{min}}`$ is the “empty beam” convergence.

The top panel shows a typical PDF from our simulations vs. the prediction from the model. The shape of the model PDF agrees with the simulation value to within $`5\%`$ near the PDF peak and $`20\%`$ in the high $`\mu `$ tail. However, we found that the model PDF was often shifted (in either direction) on the horizontal axis with respect to that of the simulation. The magnification $`\mu `$ at the PDF peaks was typically within $`0.5\%`$ but disagreed by as much as $`1.5\%`$ for extreme cases, such as the one shown here. This is a substantial disagreement given the few percent amplitude of the lensing signal. Even so, the shape of the PDF is still well fit, as can be easily seen if the model PDF is shifted to lower $`\mu `$ so that the curves overlap.

### 6.3. Multiple lensing events

During the ray-tracing runs, we tracked the number of lensing events above selected threshold magnitudes $`\kappa _m`$ for each ray. A lensing event is defined here as occurring when a ray experiences a shear on a single lens plane that causes a change in the convergence at the observer in excess of $`\kappa _m`$, so that $`g(\chi _m,\chi _n)𝐔_m>\kappa _m`$ in Eq. (11). In Table 3, we present the results from ten simulations with a total of more than $`4\times 10^7`$ lines of sight. It is quite uncommon for a given light-ray to experience multiple lensing events of a large magnitude, which suggests that a given large peak in the convergence maps is nearly always associated with a single localized massive object in the N-body simulation.

## 7. Conclusions

Weak gravitational lensing has become a powerful tool for observing the large-scale matter distribution in the universe, with rapid advances in observational capabilities hinting at even greater things to come. As observers achieve ever greater precision, it’s important that the theory keep pace.

In this paper, we have described a multi lens-plane algorithm for generating maps of weak lensing distortion from structure generated by N-body simulations paying particular attention to the numerical convergence of the algorithm(s). While our findings are substantially in line with previous works, we have noticed that biases can creep in if not enough lensing planes or time dumps are used. As concerns the range of validity of the generated maps, we find that the small-scale spatial resolution is well described by Eq. (20). With modern, parallel, codes the N-body simulation is not a limiting factor in lensing studies, suggesting that grids of models could be run relatively inexpensively.

The basic multi lens-plane algorithm introduces two numerical artifacts to the calculation of distortion maps (other than finite resolution of the grids used in computations). One is the introduction of structure discontinuities at the box boundaries in the 3-dimensional matter distribution used to make the lens planes, and the other is that light-rays are not truly perpendicular to the lens planes. We have examined the effect of both of these artifacts using techniques specifically designed to avoid them, by tracing through the boxes at an angle in the first case, and using a 3-dimensional simulation in the second, and we have shown that they do not introduce significant changes in the maps.

We examined the applicability of the Born approximation to weak lensing using measures of the statistical properties of whole maps, such as the power spectrum and the skewness, and confirm that they are in generally good agreement with semi-analytic predictions. We also confirm that the Born approximation is valid even near peaks in the convergence, where the most non-linear lensing events in our simulations are likely to occur.

Weak gravitational lensing is a tool that is providing insight into the nature of our universe at a rate that is expected to increase dramatically in coming years, and simulations will play a useful role. We have demonstrated that simulations of weak lensing based on a simple implementation of the multiple lens-plane algorithm generate maps of lensing distortion that are a good approximation of maps created using more computationally intensive methods. We have also shown that the numerical resolution of these simulations is well approximated by a simple formula. This suggests that a relatively modest investment in computation can provide the highly accurate, simulated maps which are crucial to developing algorithms and understanding data in this field.

The simulations used here were performed on the IBM-SP2 at the National Energy Research Scientific Computing Center. This research was supported by the NSF and NASA.