# The Ray Bundle method for calculating weak magnification by gravitational lenses

## 1 Introduction

Gravitational lensing is the study of the effects of matter on the propagation of light. The most obvious observational results are the production of multiple images (as first seen with the multiply imaged quasar 0957+561 by Walsh, Carswell & Weymann \[Walsh et al. 1979\]), the creation of giant luminous arcs (first identified in the galaxy clusters Abell 370 and Cl 2244 by Lynds & Petrosian \[Lynds & Petrosian 1986\] and, independently, in Abell 370 by Soucail et al. \[Soucail et al. 1987\]) and the large magnifications of source flux seen in microlensing events (for example, brightening of a single image of the multiply-imaged quasar 2237+0305 due to compact objects in the lensing galaxy, first detected by Irwin et al. \[Irwin et al. 1989\]).

The paths of light rays from a source to the observer are conveniently described with the gravitational lens equation (Equation (1) below). This equation is highly non-linear, so that, except for a small number of specific cases, there are no analytic solutions. In particular, there is no straightforward result which determines the image locations or magnifications for an ensemble of many lenses. This presents a serious problem when we wish to study the lensing properties of a complex lensing structure such as the Universe.

Fortunately, a number of approximate numerical methods have been developed which allow us to calculate magnifications and other properties of a collection of lenses. Foremost amongst these are the Ray Shooting methods, introduced by Paczyński \[Paczyński 1986\] and Kayser, Refsdal & Stabell \[Kayser et al. 1986\] and developed by Schneider & Weiss \[Schneider & Weiss 1987\], Kayser et al. \[Kayser et al. 1989\] and Lewis et al. \[Lewis et al. 1993\]. For a discussion of other methods, such as the use of a scalar deflection potential, and the optical scalar equations, see Schneider, Ehlers & Falco \[Schneider et al. 1992\].

In this paper, we present an alternative technique for calculating the magnification properties of an ensemble of lenses – the Ray Bundle method (RBM). The RBM is particularly well suited to studies of the weak lensing limit, where we are not concerned with the creation of multiple (comparably bright) images. Like the Ray Shooting method, the RBM uses backwards propagation of light rays from the observer to the source, which are deflected by the distribution of lenses, and are mapped to the source plane. Whereas the Ray Shooting method collects the deflected light rays within a rectangular grid of pixels, we consider an infinitesimal bundle of rays which form a circular image, and maintain this association to the source plane.

At first sight, the Ray Bundle method may seem to be less computationally efficient than the Ray Shooting method, as a grid based technique lends itself to (fast) hierarchical calculations \[Wambsganss 1990\]. However instead of requiring $`\stackrel{>}{}\mathrm{\hspace{0.17em}100}`$ light rays per source grid to reduce the statistical error, we need only use a bundle of $`N_{\mathrm{ray}}=8`$ rays to obtain magnifications which are correct to better than 5 per cent for an equivalent source size.

With a Ray Shooting method, we have both image and source pixels, but cannot easily determine the correspondence between them. By keeping track of the individual light bundles with the RBM, we are able to monitor the shape distortions of the beam caused by the shear and convergence of a lens ensemble. This is of particular interest when the RBM is applied to multiple lens plane geometries (as are conventionally used for studying cosmologically distributed lenses). Details of the beam shape allows for the opportunity to make comparisons between results calculated with the gravitational lens equation, and those using the optical scalar equations (which are more easily applied to smooth mass distributions, or approximate mass distributions such as the Swiss cheese model \[Einstein & Straus 1945\]).

In Section 2 we introduce various basic results of gravitational lensing, particularly with regards to magnification and the magnification probability distribution. The Ray Bundle method is introduced in Section 3 and compared in detail with analytic solutions for the Schwarzschild lens model. By obtaining the magnification probability distribution with both the Ray Bundle and Ray Shooting methods for a variety of lens geometries, we demonstrate the general applicability of the RBM. We do not discuss applications of the RBM here, but reserve details for Fluke, Webster & Mortlock (in preparation) where the Ray Bundle method is used to investigate weak lensing within realistic cosmological models (generated with N-body simulations).

## 2 Gravitational Lensing

We present here a number of important results from gravitational lensing which we will require: the gravitational lens equation, the magnification of source flux and the magnification probability distribution. Several excellent sources exist which describe gravitational lensing far more comprehensively than may be discussed here, for example Schneider et al. \[Schneider et al. 1992\] and Narayan & Bartelmann \[Narayan & Bartelmann 1996\].

### 2.1 The gravitational lens equation

The deflection of a light ray by a massive object is conveniently expressed with the gravitational lens equation (GLE), which may be derived from simple geometrical arguments (as demonstrated in Figure 1).

The GLE relates the impact parameter, $`𝝃`$, in the lens (deflector) plane, the source position, $`𝜼`$, in the source plane and the deflection angle, $`\widehat{𝜶}(𝝃)`$, of the light ray:

$$𝜼=\frac{D_{\text{os}}}{D_{\text{od}}}𝝃D_{\text{ds}}\widehat{𝜶}(𝝃).$$
(1)   
The distances ($`D_{\text{ij}}`$) in Equation (1) are angular diameter distances between the \[o\]bserver, \[d\]eflector and \[s\]ource planes. For a given source position, an image will occur for each value of $`𝝃`$ which is a solution of Equation (1). We define the lens axis as the line from the observer through a single lens (at the origin of the lens plane), and which is perpendicular to both the lens and source planes.

The GLE may be recast into a dimensionless form by introducing the scaling lengths $`\xi _0`$ and $`\eta _0=\xi _0D_{\text{os}}/D_{\text{od}}`$. Defining $`𝒙=𝝃/\xi _0`$, $`𝒚=𝜼/\eta _0`$ and the dimensionless deflection angle

$$𝜶(𝒙)=\frac{D_{\text{od}}D_{\text{ds}}}{D_{\text{os}}\xi _0}\widehat{𝜶}(\xi _0𝒙)$$
(2)   
it follows that

$$𝒚=𝒙𝜶(𝒙).$$
(3)   
### 2.2 Magnification

For a bundle of light rays passing through a transparent lens, the number of photons is conserved, ie. gravitational lensing does not change the specific intensity of the source. A change in flux, however, can occur as the cross-sectional area of a bundle of light rays will be affected by a gravitational lens. The change in the apparent luminosity is entirely due to the change in the solid angle that the image covers (at the expense of the rest of the sky) – with a lens present, more photon trajectories are brought to the observer’s eye than if there were no lens \[Dyer & Roeder 1981\].

If the flux at a frequency $`\nu `$ is $`S_\nu =I_\nu d\mathrm{\Omega }_{\mathrm{obs}}`$, where $`I_\nu `$ is the specific intensity and $`d\mathrm{\Omega }_{\mathrm{obs}}`$ is the solid angle subtended by the source at the observer’s location, then the magnification is

$$|\mu |=\frac{S_\nu }{S_\nu ^{}}=\frac{\mathrm{d}\mathrm{\Omega }_{\mathrm{obs}}}{\mathrm{d}\mathrm{\Omega }_{\mathrm{obs}}^{}}$$
(4)   
with primes denoting quantities when the lens is absent.

It is usual to measure the magnification with respect to either a ‘full beam’ or an ‘empty beam’. The full beam refers to the case where matter is smoothly distributed everywhere (inside and outside of a bundle of light rays), so that there is no shear and the magnification is due entirely to convergence. An empty beam corresponds to the case where there is a smooth distribution of matter (on average) external to the beam, and no matter within the beam. The empty beam is the maximally divergent beam, as there is now no magnification due to convergence, and the minimum magnification is $`\mu _{\mathrm{empty}}=1`$. The addition of material outside of the beam will result in a total magnification $`\mu 1`$ due to shear \[Schneider 1984\]. In the work that follows, we will calculate magnifications with respect to an empty beam.

### 2.3 The magnification probability distribution

The (differential) magnification probability, $`p(\mu ,z)\mathrm{d}\mu `$, is the probability that the total magnification of a source at redshift $`z`$ will lie in the range $`\mu `$ to $`\mu +\mathrm{d}\mu `$. The probability is subject to the constraints

$$_1^{\mathrm{}}p(\mu ,z)d\mu =1$$
(5)   
and

$$_1^{\mathrm{}}p(\mu ,z)\mu d\mu =\mu (z),$$
(6)   
where $`\mu (z)`$ is the mean magnification. For a universe where the matter is smoothly distributed, $`\mu (z)=1`$, while for the empty beam magnifications considered here, $`\mu (z)>1`$. Equation (5) provides the normalisation of the probability, while Equation (6) expresses the conservation of flux \[Weinberg 1976\].

These constraints are only strictly true when we consider the magnification probability over the whole sky – in the tests we conduct here, we use artificial lens distributions that do not cover the entire celestial sphere.

The probability distribution provides a straightforward way of making a comparison between different techniques for solving the GLE. Although we might expect the calculation of the magnification along a particular line-of-sight to vary slightly (depending on source geometry, etc), the statistical properties of a distribution of sources should be essentially independent of the method.

In practice, we will use a histogram of the probability as a function of $`\mu `$. This involves summing over a range of magnifications from $`\mu _1`$ to $`\mu _2=\mu _1+\mathrm{\Delta }\mu `$, so that each histogram bin represents

$$p(\mu _1,z)\mathrm{\Delta }\mu =_{\mu _1}^{\mu _2}p(\mu ,z)d\mu .$$
(7)   
In the following, we will refer to this as the magnification probability histogram (MPH).

### 2.4 Solving the gravitational lens equation

The GLE is highly non-linear, and in general there are multiple solutions for the image locations for a single source position. As a result, analytic solutions (where we can invert the GLE to solve for all $`𝝃`$ as a function of $`𝜼`$) exist for only special lens geometries and models, which may not always be realistic or practical. Instead, various numerical methods have been developed which enable solutions to the lens equation to be found accurately and efficiently.

Some of the most widely used methods are those based on the concept of Ray Shooting. The basic principle of the Ray Shooting method (RSM) is to use Equation (1) to propagate light rays backwards from the observer through a sequence of one or more lens planes to the source plane, where the rays are collected on a rectangular pixel grid. Typically $`10^6`$ pixels are required, with an average of $`\overline{N}100`$ rays per source plane pixel. The magnification in each source pixel $`(i,j)`$ is then proportional to the density of rays collected therein:

$$\mu _{\mathrm{pixel}}(i,j)=N_{\mathrm{collected}}(i,j)/\overline{N}.$$
(8)   
A shooting region in the lens plane is chosen such that only a few rays from outside of this grid are mapped onto one of the sources, otherwise there would be missing flux. By using a uniform grid of image rays, the relative error is approximately $`\overline{N}^{3/4}`$ \[Kayser et al. 1986\], which is better than the Poisson error if a random distribution of image rays was used (however a regular grid may introduce systematic errors). As a consequence of this error, it is possible to get magnifications $`\mu _{\mathrm{pixel}}(i,j)<1`$, violating the condition on the total magnification (Section 2.2), which is entirely a numerical effect.

Early versions of this method were introduced by Paczyński \[Paczyński 1986\] and Kayser et al. \[Kayser et al. 1986\] and developed by Schneider & Weiss \[Schneider & Weiss 1987\] and Kayser et al. \[Kayser et al. 1989\]. Hierarchical tree methods were applied to microlensing scenarios by Wambsganss \[Wambsganss 1990\] and Wambsganss, Paczyński & Katz \[Wambsganss et al. 1990\]. The hierarchical methods approximate the effects of lenses which are far from a light ray, and allow the inclusion of many thousands of lenses at a low computational cost (O($`N\mathrm{log}_2N`$) for a tree-code with N lenses, versus O($`N^2`$) when the contribution of every lens is explicitly calculated). An improvement on conventional Ray Shooting methods for obtaining statistical properties of microlensing light curves can be made with the efficient one-dimensional contour following algorithm of Lewis et al. \[Lewis et al. 1993\]. We now introduce the Ray Bundle method as an alternative to the RSM.

## 3 The Ray Bundle Method

The Ray Bundle method (RBM) is similar to the RSM, in that the lens equation is used to propagate light rays backwards from the observer to the source plane. However we now consider a bundle consisting of a central ray (the null geodesic) surrounded by $`N_{\mathrm{ray}}`$ light rays, which create an image shape (usually circular). As the ray bundle passes through the lens distribution, its shape will be distorted due to shear (stretching along an axis) and convergence (focusing due to matter within the beam). For an ‘infinitesimal’ ray bundle, the magnification is determined from Equation (4) by calculating the area of the bundle in the image plane ($`d\mathrm{\Omega }_{\mathrm{obs}}^{}`$) and the source plane ($`d\mathrm{\Omega }_{\mathrm{obs}}`$).

Since we are using backwards ray-tracing through a single image position, we do not know where other images may occur – hence any measurement of the magnification using the RBM will underestimate the total magnification. However, this is only a significant problem when we consider images located near the critical curves, when the contribution to the total magnification due to any other images becomes important (see below). The RBM was developed for applications in the weak lensing limit, and should be used with caution for strong lensing cases.

### 3.1 Comparison with analytic solutions

We can investigate the validity of the RBM by comparison with the various analytic solutions which exist for the Schwarzschild lens (see Appendix A for a summary).

Consider first a circular source of radius $`R_\mathrm{s}`$ with centre at $`𝒚_\mathrm{c}=(y_{1,\mathrm{c}},y_{2,\mathrm{c}})`$. The circumference of the source is then described by the set of vectors $`𝒚=(y_1,y_2)`$ with

$`y_1`$ $`=`$ $`y_{1,\mathrm{c}}+R_\mathrm{s}\mathrm{cos}(\varphi )`$   
$`y_2`$ $`=`$ $`y_{2,\mathrm{c}}+R_\mathrm{s}\mathrm{sin}(\varphi )`$ (9)   
where $`0\varphi <2\pi `$.

For each $`𝒚`$, we can solve for the two solutions, $`𝒙_\pm `$, with Equation (19). In this case we are using the GLE to map from the source plane to the image plane. A source far from the lens axis produces one highly demagnified image ($`\mu _{\mathrm{faint}}`$) located near the lens axis (at $`𝒙_{\mathrm{faint}}`$). The second image will have a magnification $`\mu _{\mathrm{bright}}1`$, and an angular position near the source (at $`𝒙_{\mathrm{bright}}`$). As the source is moved towards the lens axis, the images are stretched in the tangential direction<sup>1</sup><sup>1</sup>1The Schwarzschild lens has a single (degenerate) caustic point at $`y=0`$, which corresponds to a tangential critical curve at $`x=1`$ (the Einstein ring). and become comparable in brightness ($`|\mu _{\mathrm{faint}}||\mu _{\mathrm{bright}}|1`$ as $`𝒚_c0`$). When $`𝒚_c=0`$, the two images merge into a highly magnified ring (the Einstein ring) with total magnification given by Equation (24).

Now consider a circular image of radius $`R_\mathrm{i}`$ centred on the location of the bright image, $`𝒙_\mathrm{c}=𝒙_{\mathrm{bright}}(𝒚_\mathrm{c})`$, with circumferential points $`𝒙=(x_1,x_2)`$

$`x_1`$ $`=`$ $`x_{1,\mathrm{c}}+R_\mathrm{i}\mathrm{cos}(\varphi )`$   
$`x_2`$ $`=`$ $`x_{2,\mathrm{c}}+R_\mathrm{i}\mathrm{sin}(\varphi ).`$ (10)   
This time, the GLE maps in the opposite direction – from the image plane to the source plane. The source shape we obtain is stretched along the radial direction, and differentially compressed in the tangential direction.

Figure 2 demonstrates the differences between the shape and locations of a circular source (solid line), for which the circumferential points may be determined for both images, and a circular image and its corresponding source (short dashed line). In this example, we have considered a source which is near the Einstein radius where strong lensing effects dominate, and there may be a significant contribution to the flux from the second image. In all cases where we will apply the RBM, the image is chosen to be well away from the Einstein radius (critical curve), and so the flux lost from the second image is not important, as we now show.

### 3.2 The magnification deficit

We now look at how accurately the RBM approximates the total magnification, even though it includes the contribution of only one image. If we set $`R_\mathrm{s}=R_\mathrm{i}`$, as we expect only small changes to the shape and hence radius of the image in the weak lensing limit ($`|𝒙_\mathrm{c}|1`$), then this is a two parameter problem ($`R_\mathrm{i},N_{\mathrm{ray}}`$).

Defining the RBM magnification in terms of the ray bundle image and source areas ($`A_{\mathrm{i},\mathrm{RBM}},A_{\mathrm{s},\mathrm{RBM}}`$) as

$$\mu _{\mathrm{RBM}}=\frac{A_{\mathrm{i},\mathrm{RBM}}}{A_{\mathrm{s},\mathrm{RBM}}}$$
(11)   
and the true (total) magnification as

$$\mu _{\mathrm{true}}=\frac{A_{\mathrm{faint}}+A_{\mathrm{bright}}}{A_\mathrm{s}},$$
(12)   
where $`A_{\mathrm{faint}}`$, $`A_{\mathrm{bright}}`$ are the areas of the two images, then the relative error in $`\mu _{\mathrm{RBM}}`$ is

$$\frac{\mathrm{\Delta }\mu _{\mathrm{RBM}}}{\mu _{\mathrm{true}}}=\frac{|\mu _{\mathrm{true}}\mu _{\mathrm{RBM}}|}{\mu _{\mathrm{true}}}.$$
(13)   
Due to the circular symmetry of the Schwarzschild lens model, we need only determine the radius, $`x_{\mathrm{cut}}`$, within which the RBM produces a relative error $`\frac{\mathrm{\Delta }\mu _{\mathrm{RBM}}}{\mu _{\mathrm{RBM}}}>p`$ per cent.

By using $`N_{\mathrm{ray}}`$ rays in the image and source bundles, we are approximating the shape of a circular image/source by a polygon with $`N_{\mathrm{ray}}`$ sides. Clearly, when $`N_{\mathrm{ray}}1`$, we will have a reasonable approximation to the true shape of the image/source. However, to improve the speed of the ray bundle method (at the cost of a small error), we ideally want to select a small value of $`N_{\mathrm{ray}}`$ ($`\stackrel{<}{}\mathrm{\hspace{0.17em}20}`$). The areas are calculated as a sum of triangular components within the image/source polygon, where each triangle has a common vertex at $`𝒚_\mathrm{c}`$ or $`𝒙_\mathrm{c}`$ (ie. the null geodesic).

We need to first check that the calculated $`\mu _{\mathrm{true}}`$ is not significantly in error using a particular value of $`N_{\mathrm{ray}}`$. This was achieved by numerically integrating Equation (23), and comparing with $`\mu _{\mathrm{true}}`$ for $`4N_{\mathrm{ray}}256`$. The relative error in $`\mu _{\mathrm{true}}`$ was found to be independent of the choice of $`N_{\mathrm{ray}}`$, and was well below the percentage cut-off level selected for $`\mu _{\mathrm{RBM}}`$ at the same radii.

For source positions which are near the lens axis, it is difficult to keep track of which solutions $`x_\pm `$ belong to which of the images, particularly when merging of images is taking place. In these cases, $`\mu _{\mathrm{true}}`$ is a misnomer, as it can be significantly higher than the total magnification from Equation (23), as may be seen in Figure 3. $`\mu _{\mathrm{RBM}}`$ is subject to a similar error, as the areas here are also being calculated on the basis of $`N_{\mathrm{ray}}`$ rays. We can proceed under the assumption that $`\mu _{\mathrm{true}}`$ is sufficiently accurate for comparison with $`\mu _{\mathrm{RBM}}`$ when the two images are separable.

By randomly selecting $`5000`$ source (RBM image) locations, we calculate $`\mu _{\mathrm{true}}`$ and $`\mu _{\mathrm{RBM}}`$. Table 1 shows the dimensionless radius, $`x_{\mathrm{cut}}`$, within which $`\mathrm{\Delta }\mu _{\mathrm{RBM}}/\mu _{\mathrm{true}}>p`$ per cent for a range of source radii ($`10^5<R_\mathrm{i}<0.01`$). For a given value of the relative error, these results are essentially independent of both the source radii and $`N_{\mathrm{ray}}`$ in the range $`4N_{\mathrm{ray}}256`$.

The choice of $`N_{\mathrm{ray}}=8`$ appears to be the best compromise between accuracy and speed for our later investigations of the RBM. The smaller $`N_{\mathrm{ray}}`$, the fewer total deflection calculations which need to be made, yet we retain high accuracy for $`\mu _{\mathrm{RBM}}`$ in the weak lensing limit. With 8 rays, we have a symmetric image with a quadrupole component, so it is possible to determine the distribution of image ellipticities.

In Figure 3, we plot $`\mu _{\mathrm{RBM}}`$ (thick solid line), $`\mu _{\mathrm{true}}`$ (thin solid line) and $`\mu _{\mathrm{integral}}`$, the solution of Equation (23) (dashed line), as functions of the source impact parameter ($`y=|𝒚|`$). The source/image bundle radius in this case was $`R_\mathrm{s}=R_\mathrm{i}=0.4`$. We have used a larger radius than would be practical for the RBM in order to show in more detail what the high magnification behaviour is like near $`y=0`$. As the source is comparable in size to the Einstein radius, the calculated magnifications near $`y0.4`$ are noisy. In this case, a circular source should produce a pair of highly distorted arcs, but this is not well represented by the use of triangular components within the ray bundle (ie. parity changes can occur in the individual triangles which make up the bundle). The vertical dotted line shows where the relative error in $`\mu _{\mathrm{RBM}}`$ is accurate to $`10`$ per cent, at $`y_{\mathrm{cut}}=1.17`$ or $`x_{\mathrm{cut}}=1.74`$ (note that these values are higher than in Table 1 due to the larger source radius).

This test has shown that the Ray Bundle method can give a highly accurate value of the magnification (relative error better than $`5`$ per cent) in the weak lensing limit. However, we have modified the Ray Bundle method slightly by requiring that the image position agrees with the brighter image of the corresponding circular source. In a true application of the RBM, we randomly select image positions without the a priori knowledge of source locations.

Figure 4 shows the RBM magnification for bundle positions near, and within, the critical curve of a single Schwarzschild lens (thick solid line). For comparison, the total magnification of a point source with one image at the same location as the RBM image is shown as the dashed line. Bundle positions within $`x=0.5`$ produce magnifications $`\mu _{\mathrm{RBM}}1`$. For these bundles we actually selected the fainter of the images, as opposed to the case considered previously, where we purposefully selected the bundle corresponding to the brighter image. Solutions of the lens equation, Equation (19), which produce an image near the lens ($`x_{}\stackrel{<}{}\mathrm{\hspace{0.17em}1}`$) correspond to total magnifications near $`\mu =1`$, as the source is far from the caustic point ($`y1/x_{}>1`$). This is demonstrated with the thin solid line in Figure 4, where we plot ray bundles with $`\mu 1`$ within the Einstein radius as $`\mu ^{}=\mu +1`$.

It is clear that a small strip of image locations near the critical curve is responsible for the high magnification region of the MPH, for which the Ray Bundle method is not well suited.

Images well within the Einstein radius correspond to sources far from the lens axis. Such images are the faint images described in Section 3.1 (at $`𝒙_{\mathrm{faint}}`$), and so there will be a second image at $`𝒙_{\mathrm{bright}}`$ which contributes the majority of the flux. Although it is not possible to solve for $`𝒙_{\mathrm{bright}}`$ given $`𝒙_{\mathrm{faint}}`$, if the entire image plane is well sampled with ray bundles, then we can expect that another bundle will pass through $`𝒙_{\mathrm{bright}}`$ and the source magnification will be calculated on the basis of this second bundle only.

Two restrictions are now imposed on the Ray Bundle method for its later application to ensembles of lenses, which serve to complete the definition of the method. Firstly, image positions within the Einstein radius, or equivalently, any magnifications which are calculated to be $`\mu <1`$ are discarded. Secondly, after selecting a relative error for $`\mu _{\mathrm{RBM}}`$, we do not include images which fall within $`x_{\mathrm{cut}}`$ Einstein radii of any given lens.

Having shown that in the ‘weak’ lensing limit ($`|𝒙_\mathrm{c}|\stackrel{>}{}\mathrm{\hspace{0.17em}2}`$), we can be sure of calculating the magnification to within $`5`$ per cent (or better) using $`N_{\mathrm{ray}}4`$ we can now proceed to a statistical comparison between the Ray Shooting and Ray Bundle methods.

### 3.3 Comparison of magnification probability histograms

Next we compare the MPH obtained with the RBM and the RSM. A number of subtle differences exist between the two methods, even when applied to the same lens model. The source size investigated is limited by the number of pixels in the source plane for the RSM. For a grid of $`N_{\mathrm{pix}}\times N_{\mathrm{pix}}`$ covering a square region $`2y_{\mathrm{max}}\times 2y_{\mathrm{max}}`$, the source ‘radius’ is $`R_\mathrm{s}=y_{\mathrm{max}}/N_{\mathrm{pix}}`$.

We choose the RSM sources to be squares with a side-length equal to the diameter of the circular RBM image bundles. Since the magnification decreases as the source area is increased, we expect each RSM source to have a systematically smaller magnification than the corresponding RBM image. It is possible to use a much finer resolution grid for the ray shooting, and then integrate over a larger source size, but we have elected not to do this. This decision was based on a comparison of the computation time: for our implementation of the RSM, a grid of $`1000\times 1000`$ pixels, with an average of $`\overline{N}=250`$ rays per pixels took approximately eight hours of computation time. An equivalent number of images (where $`N_{\mathrm{ray}}=8`$) is completed in 1 minute with the RBM<sup>2</sup><sup>2</sup>2More computationally efficient versions of the RSM are available, which can produce the same level of resolution in a time comparable to that of the RBM (J. Wambsganss, private communication). 

Due to the distortion of rays near the boundary, as shown in Figure 5, we actually shoot rays with the RSM through a larger angular region in the image plane. This prevents us from including source pixels which are not well sampled by rays.

The accuracy of the RSM magnifications depends on the average number of rays collected in each pixel, $`\overline{N}`$. Since we are not implementing a (fast) hierarchical method, the time required to obtain the magnification distribution is proportional to the total number of rays, $`N_{\mathrm{RSM}}=\overline{N}N_{\mathrm{pix}}^2`$. The RSM does not fully sample the highest magnification region $`(\mu \mu _{\mathrm{max}})`$ due to the regular placement of sources on a grid.

With the RBM image size fixed by $`R_\mathrm{i}=R_\mathrm{s}`$, we are left to choose the number of rays which make up the ray bundle, $`N_{\mathrm{ray}}`$, and the number of images. Ideally we want the images to completely cover the image plane (which again has a slightly larger angular size than the source plane due to the deflection of light rays near the boundaries) which requires

$$N_{\mathrm{image}}f\frac{\pi R_\mathrm{i}^2}{(2y_{\mathrm{max}})^2},$$
(14)   
where $`f`$ is (approximately) the fraction of the source plane covered by sources. For the RBM then, the total number of rays required is $`N_{\mathrm{RBM}}=N_{\mathrm{image}}\times N_{\mathrm{ray}}`$.

The comparative computational speed and resulting magnification accuracy may be obtained by requiring $`N_{\mathrm{RBM}}=N_{\mathrm{RSM}}`$. Setting $`f=1`$ and $`R_\mathrm{i}=R_\mathrm{s}`$, we have

$$\overline{N}=\frac{4}{\pi }N_{\mathrm{ray}}.$$
(15)   
We have already seen that $`N_{\mathrm{ray}}=8`$ is a suitable choice, which gives $`\overline{N}10.2`$. This corresponds to a relative error $`\overline{N}^{3/4}17`$ per cent. Although we can relax the constraint on the fraction of the source plane covered somewhat, and still have a well sampled MPH with the Ray Bundle method, we do not have this flexibility with the grid based Ray Shooting method. In addition, as we decrease the source size, the number of pixels required for the RSM increases, and a higher density of rays is necessary.

Figure 6 shows the Magnification Probability histograms obtained for a single Schwarzschild lens using the Ray Shooting (thin line) and Ray Bundle (thick line) methods. The parameters for each method are listed in Table 2. The vertical axis of this (and later histograms) is the normalised number of bundles/sources in each magnification bin, $`N(\mu )`$, which is equivalent to the definition of $`p(\mu _1)\mathrm{\Delta }\mu `$ in Equation (7).

A cut-off was imposed on image locations for the RBM at $`x_{\mathrm{cut}}=1.01`$ Einstein radii. The two distributions are qualitatively very similar, when the various caveats described above are considered. Imposing a larger value of $`x_{\mathrm{cut}}`$ serves to reduce the maximum magnification with the RBM. The poor sampling in the highest $`\mu `$ bins for the RSM is clearly demonstrated. The dashed line shows the expected $`\mu ^2`$ power law slope of the MPH at large $`\mu `$ (see Appendix A), and both distributions have this approximate form (although for the RSM the statistical significance of the histogram bins with $`\mu >10`$ is low).

As expected, the RSM produces magnifications which are $`\mu <1`$ (when $`\mu _{\mathrm{empty}}1`$ is expected) due to numerical effects. The sample mean and variance of the two distributions are $`\mu =0.98`$ and $`\sigma _\mu ^2=0.03`$ for Ray Shooting, and $`\mu =1.02`$ and $`\sigma _\mu ^2=0.08`$ for the Ray Bundle method<sup>3</sup><sup>3</sup>3We reiterate that in the definition of the RBM, we do not include images with $`\mu <1`$, see end of Section 3.2). . The Ray Shooting method provides higher accuracy at high magnifications, but at magnifications $`\mu 1`$, the Ray Bundle method is more accurate even though flux from additional images is neglected.

One aspect of the MPH we have not yet discussed has to do with the weighting we apply to each ray bundle. For the RSM, every source is a pixel with the same area. For the RBM, the initial ray bundles (images) have the same area, but the resulting sources must have different areas (by the definition of a magnification). It is sufficient to weight each ray bundle in the MPH by the area of the resulting source. Figure 7 shows the RBM (thin line) with no weighting compared with the correct area weighting (thick line).

### 3.4 More complex lens distributions

For an ensemble of $`N`$ lenses in the lens plane, each with mass $`M_j`$, the total deflection angle generalises to

$$\widehat{𝜶}(𝝃)=\underset{j=1}{\overset{N}{}}\frac{4GM_j}{c^2}\frac{(𝝃𝝃_j)}{|𝝃𝝃_j|^2},$$
(16)   
where the $`(𝝃𝝃_j)`$ are the impact parameters to each lens.

Consider the case where lenses are restricted to lie within a rectangular region<sup>4</sup><sup>4</sup>4This is the case in common implementations of Ray Shooting as it allows for the easy implementation of fast, hierarchical methods \[Wambsganss 1990\].. Light rays passing through one of the corners of the shooting region will necessarily be deflected inwards by the mass distribution. This is appropriate for an isolated configuration of lenses, such as in studies of the microlensing effect of many stars which make up a galaxy (where the contribution to the deflection by an external mass distribution may be modelled by adding a shear term to the lens equation). However, for an investigation of the lensing due to large scale structure, where the mass distribution is assumed to be continuous and homogeneous in all directions about a ray, we may introduce an artificial shear on rays near the shooting boundary.

For the RBM, we choose to calculate each of the deflection angles explicitly, with an increase in the computational time over an equivalent hierarchical method. By making the direct calculation of the deflection, we are free to choose the geometry of the region within which we include lenses. The most natural choice for a distribution which is homogeneous and isotropic beyond some length scale $`R_\mathrm{H}`$ is to include lenses within a circular region around each ray out to the radius $`R_{\mathrm{lens}}=R_\mathrm{H}`$. For the isolated lens geometries we examine here, we need only set $`R_{\mathrm{lens}}`$ to encompass the RSM image plane. A discussion of the appropriate choice for $`R_\mathrm{H}`$ in cosmological lensing scenarios is reserved for Fluke et al. (in preparation).

As a final test here, we now consider a random distribution of $`N_{\mathrm{lens}}=5`$ lenses (with the same lens positions for both the RSM and RBM). The resulting MPH is shown in Figure 8 and various parameters in Table 2. For both methods we use a total of $`10^6`$ bundles/sources, but discarding bundles within $`x_{\mathrm{cut}}=1.01`$ Einstein radii for the RBM reduces this to $`9.3\times 10^5`$ bundles for the histogram. The sample mean and variance are $`\mu =1.09`$, $`\sigma _\mu ^2=2.09`$ for the RBM and $`\mu =0.98`$, $`\sigma _\mu ^2=0.23`$ for the RSM.

Using an ensemble of lenses introduces a new length scale to the problem, so that sub-structure at intermediate magnifications ($`3\mu 30`$) is seen in the MPH using the RSM (but not with the RBM). The ‘bump’ in the MPH occurs for a planar distribution of lenses, and has been studied both numerically \[Rauch et al. 1992\] and analytically \[Kofman et al. 1997\], and is believed to be due to the caustic patterns of pairs of point lenses. If this is a caustic-induced feature, ie. a region of high magnification, then the RBM will not provide the ‘correct’ magnification. We feel the feature may be in part due to the low resolution with which the complex caustic structure was mapped with the RSM on a regular grid ($`1000\times 1000`$ pixels), however a discussion of this effect is beyond the scope of this paper.

## 4 Conclusions

The Ray Bundle method provides a computationally fast, accurate and flexible alternative to the Ray Shooting method for studies of the weak gravitational lensing limit. A wide variety of lens models are easily incorporated – we have considered here only the case of the Schwarzschild lens, but changing to a different model involves modifying only $`\widehat{𝜶}(𝝃)`$ in the GLE. One alternative to ray based methods requires solving for the scalar lensing potential, but it may not always be possible to find an analytic solution for all lens models. The RBM also allows us to avoid artificial shear introduced by grid based methods for light rays which pass near the shooting boundary, when a small portion of an otherwise homogeneous and isotropic distribution is used.

The point source limit may be approached as values of $`R_\mathrm{i}`$ may be selected without the restriction of introducing a finer source grid, and a corresponding increase in the density of rays required (provided we can relax the constraint that the source plane must be completely covered by sources with the RBM).

The RBM should only be applied with caution to strong lensing scenarios, where the RSM is far superior. As only one image is followed to the source, there will be an error in the total magnification when the image is near a critical curve.

The RBM is, however, particularly well suited to problems were we want to investigate in detail individual lines-of-sight for various lens geometries and models. An important advantage of the RBM is that we can associate a particular image position and shape with the corresponding source position and shape. This provides us with the opportunity of following the development of the shape of a ray bundle through a sequence of lens planes, as used in models of cosmological lensing (for example, Wambsganss, Cen & Ostriker \[Wambsganss et al. 1998\]).

An important application of weak lensing is to determine the effect of small changes in the magnification of standard candle sources (such as Type Ia Supernovae) on the derived values of cosmological parameters \[Wambsganss et al. 1997\]. The high accuracy of the RBM in the weak lensing limit makes it a value tool for such studies.

## Acknowledgments

The authors would like to thank Peter Thomas and Andrew Barber (University of Sussex), and Hugh Couchman (University of Western Ontario) for helpful discussions. The authors are grateful to the referee, Joachim Wambsganss, for his insightful comments. CJF and DJM are funded by Australian Postgraduate Awards. CJF is grateful for financial assistance from the University of Melbourne’s Melbourne Abroad scholarships scheme, and the Astronomical Society of Australia’s travel grant scheme.

## Appendix A The Schwarzschild Lens

The simplest lens model is the point-mass or Schwarzschild lens (for example Schneider et al. \[Schneider et al. 1992\], Narayan & Bartelman \[Narayan & Bartelmann 1996\]) , for which the deflection angle due to a mass, $`M`$, is

$$\widehat{𝜶}(𝝃)=\frac{4GM}{c^2|𝝃|^2}𝝃.$$
(17)   
The dimensionless lens equation for a point source with this model is

$$y=x1/x$$
(18)   
so that there are two images (one located on either side of the lens, and co-linear with the lens) at

$$x_\pm =\frac{1}{2}\left(y\pm \sqrt{y^2+4}\right)$$
(19)   
with corresponding magnifications

$$\mu _\pm =\frac{1}{2}\left(\frac{y^2+2}{y\sqrt{y^2+4}}\pm 1\right).$$
(20)   
The total magnification is the sum of the absolute values of the individual image magnifications: $`\mu _\mathrm{p}=|\mu _+|+|\mu _{}|`$. When the source is far from the lens axis ($`y1`$), one of the images (say, $`x_{}`$) will be significantly demagnified ($`\mu _{}1`$). The total magnification is then $`\mu _\mathrm{p}\mu _+`$.

It is possible to derive an analytic form for the magnification probability, $`p(\mu ,z)`$, for the Schwarzschild lens for large values of $`\mu _\mathrm{p}`$ \[Schneider et al. 1992\], which is reasonably generic for most lens models \[Peacock 1982\]:

$$p(\mu )\mu _\mathrm{p}^3.$$
(21)   
On integrating to form the MPH we have

$$p(\mu )\mathrm{\Delta }\mu =_{\mu _1}^{\mu _2}p(\mu )d\mu =\frac{1}{\mu _1^2}\left(1\frac{\mu _1^2}{\mu _2^2}\right).$$
(22)   
If the histogram bins are equally spaced logarithmically, $`\mu _1/\mu _2`$ is constant, and so we have an additional constraint that the MPH must have a power law slope $`2`$ for $`\mu \stackrel{>}{}\mathrm{\hspace{0.17em}10}`$.

For an extended source, the total magnification is the integral of $`\mu _\mathrm{p}`$ over the source, weighted by the intensity profile, $`(𝒚)`$. For a circular source with dimensionless radius $`R_\mathrm{s}=R_\mathrm{s}/\eta _0`$, and a uniform intensity profile we have (eg. Pei \[Pei 1993\], or Schneider et al. \[Schneider et al. 1992\] for an alternative formulation)

$$\mu _\mathrm{e}(y)=\frac{1}{\pi R_\mathrm{s}^2}_{|yR_\mathrm{s}|}^{y+R_\mathrm{s}}dt\frac{\sqrt{t^24}\left(R_\mathrm{s}^2y^2+t^2\right)}{\sqrt{R_\mathrm{s}^2(yt)^2}\sqrt{(y+t)^2R_\mathrm{s}^2}}.$$
(23)   
Equation (23) approaches the point source solution as $`y\mathrm{}`$ or $`R_\mathrm{s}0`$, and the maximum magnification is

$$\mu _{\mathrm{e},\mathrm{max}}=\frac{\sqrt{R_\mathrm{s}^2+4}}{R_\mathrm{s}}$$
(24)   
at $`y=0`$.