text:
# Numerical Methods in Gravitational Lensing

## 1 Introduction

Only for very special lens models can numerical methods be avoided in gravitational lensing studies. There are three essential reasons for that. One is the non-linearity of gravitational lensing, i.e. the fact that image and source positions are related to one another in a non-linear fashion. This gives rise to the well-known phenomena of mutiple imaging, strong image distortions, and so forth. The second reason is that lenses exist which are themselves best described by numerical models. Galaxy clusters are one example, lensing by large-scale structures is another. Although it is true that many aspects of gravitational lensing by large-scale structures can be derived analytically, detailed simulations require numerical techniques. The third reason is that the interpretation of gravitational lensing effects or events often require the application of sophisticated algorithms to ever growing amounts of data. One example is the reconstruction of the projected mass density distribution of a galaxy cluster from the observed image distortions due to gravitational shear.

Needless to say, there are many more aspects of numerical methods related to gravitational lensing than I can cover in this review. An outstanding example are the highly elaborate methods that have been developed over recent years for determining image shapes of faint background galaxies on CCD frames, and for extracting the gravitational shear signal from them. This is a whole branch of data analysis on its own. Here, I can only deal with numerical methods for relating mass distributions to their gravitational lensing effects.

Consequently, the outline of this lecture is as follows: First, I shall discuss methods for studying individual lenses, i.e. their imaging properties, their critical curves and caustics. In particular, the use of adaptive grids and techniques for searching and characterising images will be discussed. Second, I shall describe how extended lenses can be treated numerically using the multiple-lens plane theory. This will lead to the basic equations for tracing light rays through (simulated) cosmological volumes. A large fraction of the discussion will be devoted to issues of resolution and noise, and to spurious effects in simulated lensing statistics. Finally, third, I shall describe inversion techniques, i.e. methods for reconstructing the projected mass distribution of lenses whose distortion has been measured. The classic Kaiser-Squires method will be described, and also maximum-likelihood techniques and maximum-entropy methods.

General lensing theory and the theory of weak lensing are covered by Koenraad Kuijken’s and Peter Schneider’s lectures in this volume. Basic references on lensing include the textbook by Schneider et al. (1992) and the lecture by Narayan & Bartelmann (1999), reviews of weak lensing are Mellier (1999) and Bartelmann & Schneider (2001).

## 2 Individual Lenses

### 2.1 Assumptions

A brief reminder of the basic assumptions underlying the theory of individual lenses may be in order. There are three main assumptions. First, the Newtonian gravitational potential of the lens be small, $`|\mathrm{\Phi }|c^2`$. Second, velocities in the gravitational lens system, both of constituents within the lenses and of the lenses with respect to the rest frame of the microwave background, be small $`vc`$. Third, the extent of the lenses along the line-of-sight be small compared to the other distances in the system, which are usually cosmological and thus comparable to the Hubble radius, $`c/H_0=3h^1\mathrm{Gpc}`$, with $`H_0`$ being the Hubble constant and $`h=H_0/100\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$.

It is worth noting how well these assumptions are satisfied in ordinary lensing situations. Consider a galaxy cluster with mass $`M=10^{15}h^1M_{}`$. Assuming spherical symmetry, the Newtonian potential at a distance $`R=1h^1\mathrm{Mpc}`$ from its centre is

$$|\mathrm{\Phi }|\frac{GM}{R}(2\times 10^3\mathrm{km}\mathrm{s}^1)^2,$$
(1)
evidently much smaller than the speed of light squared. A typical length scale for the radius of a galaxy cluster is $`11.5h^1\mathrm{Mpc}`$, which is several hundred times smaller than typical distances in a cluster-lensing system. Finally, peculiar velocities of galaxy clusters with respect to the Hubble flow are of order several hundred $`\mathrm{km}\mathrm{s}^1`$, and typical velocities of galaxies within galaxy clusters reach of order $`10^3\mathrm{km}\mathrm{s}^1`$, but both velocities are way below the speed of light. The above assumptions hold even better for lensing by galaxies, of course.

We can thus safely assume the above conditions to be satisfied. It is then possible to project the lensing mass distribution onto a plane perpendicular to the line-of-sight, the lens plane, and describe it by its surface mass density $`\mathrm{\Sigma }`$. Sources are assumed to be located on a corresponding plane, the source plane. A typical lens system is sketched in Fig. 1.

The three distances $`D_{\mathrm{d},\mathrm{s},\mathrm{ds}}`$ shown in Fig. 1 and explained in its caption are generally not additive because of space-time curvature, thus $`D_\mathrm{s}D_\mathrm{d}+D_{\mathrm{ds}}`$ in contrast to flat space-time.

### 2.2 Coordinates and Notation

Let us now introduce physical coordinates $`\stackrel{}{\xi }`$ and $`\stackrel{}{\eta }`$ on the lens and source planes, respectively. Alternatively, it is often convenient to introduce angular coordinates $`\stackrel{}{\theta }`$ and $`\stackrel{}{\beta }`$, which are obviously related to $`\stackrel{}{\xi }`$ and $`\stackrel{}{\eta }`$ through

$$\stackrel{}{\xi }=D_\mathrm{d}\stackrel{}{\theta },\stackrel{}{\eta }=D_\mathrm{s}\stackrel{}{\beta }.$$
(2)
Dimensional coordinates are of course not suitable for numerical calculations, which can only handle numbers. We thus have to introduce a length scale $`\xi _0`$, or alternatively an angular scale $`\theta _0`$, in the lens plane. This length scale is so far *arbitrary*. It implies a length or angular scale

$$\eta _0=\frac{D_\mathrm{s}}{D_\mathrm{d}}\xi _0,\text{or}\beta _0=\frac{\eta _0}{D_\mathrm{s}}=\frac{\xi _0}{D_\mathrm{d}}=\theta _0$$
(3)
in the source plane. Dimension-less coordinates are then defined by

$$\stackrel{}{x}=\frac{\stackrel{}{\xi }}{\xi _0}=\frac{\stackrel{}{\theta }}{\theta _0},\text{or}\stackrel{}{y}=\frac{\stackrel{}{\eta }}{\eta _0}=\frac{\stackrel{}{\beta }}{\theta _0}$$
(4)
in the lens and source planes, respectively. The numerical code will have to deal with the dimension-less vectors $`\stackrel{}{x}`$ and $`\stackrel{}{y}`$. It helps numerical accuracy greatly if these numbers are of order unity. Thus, the first challenge in setting up a lensing simulation is to choose an appropriate length- or angular scale $`\xi _0`$ or $`\theta _0`$, which should both be adapted to the physical problem at hand, and to the requirement that numerical codes work most accurately if the numbers they are dealing with are neither too large nor too small, compared to machine accuracy. Choosing unappropriate length scales can, for instance, render image searches unsuccessful.

### 2.3 The Lensing Potential

It will be convenient for the following discussion to introduce the lensing potential $`\psi `$ as the basic physical quantity for lensing studies. It is the scaled, projected Newtonian gravitational potential of the lens,

$$\psi (\stackrel{}{x})=\frac{2}{c^2}\frac{D_\mathrm{d}D_{\mathrm{ds}}}{\xi _0^2D_\mathrm{s}}\mathrm{\Phi }(\xi _0\stackrel{}{x},l)dl.$$
(5)
The so-called *reduced* (i.e. appropriately scaled) deflection angle is the gradient of the potential,

$$\stackrel{}{\alpha }(\stackrel{}{x})=_\stackrel{}{x}\psi (\stackrel{}{x}),$$
(6)
and the lensing convergence (i.e. the appropriately scaled surface-mass density) is

$$\kappa (\stackrel{}{x})=\frac{1}{2}_\stackrel{}{x}^2\psi (\stackrel{}{x})=\frac{1}{2}_\stackrel{}{x}\stackrel{}{\alpha }(\stackrel{}{x}).$$
(7)
Finally, the gravitational tidal field is described by the two-component shear,

$$\gamma _1(\stackrel{}{x})=\frac{1}{2}\left(\psi _{,11}\psi _{,22}\right)=\frac{1}{2}\left(\alpha _{1,1}\alpha _{2,2}\right),\gamma _2(\stackrel{}{x})=\psi _{,12}=\alpha _{1,2},$$
(8)
where the convention was used that $`f_{i,j}`$ is the derivative of the $`i`$-th component of $`\stackrel{}{f}`$ with respect to the coordinate $`x_j`$. It is important to note that the fact that all lensing quantities can be derived from the scalar lensing potential establishes relations between all of them. This will be exploited several times later.

Note that the lensing quantities must be rescaled in case the coordinate scale $`\xi _0`$ is changed. Suppose $`\xi _0^{}`$ is introduced instead of $`\xi _0`$. Since the physical surface-mass density of the lens must remain the same at any given physical location, the reduced deflection angle must transform as

$$\stackrel{}{\alpha }(\stackrel{}{x}^{})=\frac{\xi _0^{}}{\xi _0}\stackrel{}{\alpha }(\stackrel{}{x}),$$
(9)
and convergence and shear transform as

$$[\kappa (\stackrel{}{x}^{}),\gamma _i(\stackrel{}{x}^{})]=\left(\frac{\xi _0^{}}{\xi _0}\right)^2[\kappa (\stackrel{}{x}),\gamma _i(\stackrel{}{x})].$$
(10)
### 2.4 Imaging

Suppose now we were given some description of the lensing potential $`\psi (\stackrel{}{x})`$, or of the deflection angle $`\stackrel{}{\alpha }(\stackrel{}{x})`$. This description could be an analytical formula, or it could be in form of an array, i.e. a set of numbers given at grid points $`(x_i,x_j)`$. We wish to know how the given lens images its background.

We introduce a coordinate grid $`\stackrel{}{x}_{ij}`$ on the lens plane subject to the condition that it be sufficiently well resolved. This means that the smallest features in the lens must be covered by at least a few grid points. Since we are given the deflection angle as a function of position, we can compute a deflection-angle grid, $`\stackrel{}{\alpha }_{ij}=\stackrel{}{\alpha }(\stackrel{}{x}_{ij})`$. The mapped grid on the source plane is then simply $`\stackrel{}{y}_{ij}=\stackrel{}{x}_{ij}\stackrel{}{\alpha }_{ij}`$. This mapped grid will appear as a distorted image of the regular grid in the lens plane, as the example in the left panel of Fig. 2 shows.

The mapping process must now be reversed in order to obtain an image created by the lens. For doing so, the *source* plane is first covered with a regular grid, $`\stackrel{}{y}_{ij}^{}`$. Next, we loop over all grid points $`\stackrel{}{x}_{ij}`$ in the *lens* plane and find its mapped source point $`\stackrel{}{y}_{ij}`$ in the source plane, and search for the nearest neightbours $`\stackrel{}{y}_{kl}^{}`$ surrounding $`\stackrel{}{y}_{ij}`$ in the source plane. This is illustrated in the right panel of Fig. 2. The surface brightness of the source, known at the positions $`\stackrel{}{y}_{kl}^{}`$, can then be interpolated to $`\stackrel{}{y}_{ij}`$ and the result assigned to the image point $`\stackrel{}{x}_{ij}`$. That way, the surface brightness at all points in the lens plane can be determined, and thus the lensed image be constructed. Fig. 3 shows an example.

The left panel of the figure shows a simulated CMB temperature fluctuation field of $`10^{}\times 10^{}`$ angular size. The temperature increases from white to red. In essence, the temperature fluctuation corresponds to a fairly smooth gradient across the field. The right panel shows the gravitational lensing signature imprinted on the CMB at such angular scales by a galaxy cluster. The temperature visible at an angular position $`\stackrel{}{\theta }`$ on the sky, $`T^{}(\stackrel{}{\theta })`$, is related to the intrinsic temperature $`T`$ through $`T^{}(\stackrel{}{\theta })=T[\stackrel{}{\theta }\stackrel{}{\alpha }(\stackrel{}{\theta })]`$. Thus, the light deflection by the cluster causes the visible temperature distribution to be rearranged, yielding a highly specific pattern (Seljak & Zaldarriaga 2000).

### 2.5 Critical Curves and Caustics

As mentioned in the introduction, the deflection-angle field contains full information on the lensing mass distribution. All other quantities like convergence and shear, but also image magnifications, follow from the deflection angle via differentiation. It is thus a common task to compute numerical derivatives.

Suppose a function $`f(\stackrel{}{x})`$ is tabulated on a grid, so that we are given the values $`f_{ij}`$ at the grid points $`\stackrel{}{x}_{ij}`$. The derivative of $`f(\stackrel{}{x})`$ at a particular point $`\stackrel{}{x}_{00}`$ in the first coordinate direction is approximated by

$$\frac{f(\stackrel{}{x})}{x_1}|_{\stackrel{}{x}_{00}}=\frac{1}{2h}\left(f_{10}f_{10}\right)+𝒪(h^2)$$
(11)
where $`h`$ is the separation of the grid points in the chosen direction; cf. the left panel of Fig. 4. This centred difference has the advantage compared to the more straightforward one-sided differences $`f_{10}f_{00}`$ or $`f_{00}f_{10}`$ of being second-order in the grid separation $`h`$. There are higher-order differencing schemes using function values at more than two grid points, but the second-order scheme is usually sufficient. No lensing quantity should vary strongly between two adjacent grid points because otherwise the resolution of the grid would be grossly insufficient.

We will typically need derivatives of the deflection angle $`\stackrel{}{\alpha }`$. Since $`\stackrel{}{\alpha }`$ is itself the gradient of a scalar potential, its derivatives must satisfy $`\alpha _{1,2}=\psi _{,12}=\psi _{,21}=\alpha _{2,1}`$. It is thus usually preferable to check that this relation is satisfied within numerical accuracy, and to use $`(\alpha _{1,2}+\alpha _{2,1})/2`$ instead of either $`\alpha _{1,2}`$ or $`\alpha _{2,1}`$ alone.

Critical curves in the lens plane are defined by the condition that the Jacobian determinant of the lens mapping vanish there, $`det𝒜(\stackrel{}{x})=0`$. The elements of the Jacobian matrix are $`𝒜_{ij}=\delta _{ij}\alpha _{i,j}`$, thus the Jacobian determinant is

$$Ddet𝒜=(1\alpha _{1,1})(1\alpha _{2,2})\alpha _{1,2}^2.$$
(12)
It can be computed once the (numerical) derivatives of the both deflection-angle components have been determined.

One method of identifying grid points in the lens plane next to the critical curve proceeds as follows. Let $`S=\mathrm{sign}(D)`$, and consider one particular grid point $`\stackrel{}{x}_{00}`$ in the lens plane. The point is next to the critical curve if, and only if, the sign of the Jacobian determinant changes between it and one or more of its nearest neighbours. Hence, if the condition

$$S_{00}(S_{10}+S_{10}+S_{01}+S_{01})<4$$
(13)
is satisfied, the grid point $`\stackrel{}{x}_{00}`$ is next to a critical curve (cf. the right panel of Fig. 4). Of course, $`\stackrel{}{x}_{00}`$ is not itself *on* the critical curve, but to the positional accuracy determined by the grid resolution, the position of the critical curve can be constrained that way. Points on the source plane next to the caustic curve are then easily found via the lens equation, $`\stackrel{}{y}_{\mathrm{C}ij}=\stackrel{}{x}_{\mathrm{C}ij}\stackrel{}{\alpha }(\stackrel{}{x}_{\mathrm{C}ij})`$, where the $`\stackrel{}{x}_{\mathrm{C}ij}`$ are the grid points in the lens plane next to critical curves.

As an example, consider a lens model for a spiral galaxy, consisting of a spherical halo and a flat disk seen almost edge-on (Bartelmann & Loeb 1998). The deflection-angle field of such a lens can be given analytically (cf. Keeton & Kochanek 1998). Convergence and total shear $`(\gamma _1+\gamma _2)^{1/2}`$ as determined by numerical differentiation are shown together with the modulus of the deflection angle in Fig. 5.

The critical curves and caustics of that lens model as determined with the method described above are shown in Fig. 6.

### 2.6 Adaptive Source Grids

One of the most prominent goals of gravitational lensing studies with individual strong lenses is to determine the imaging statistics of a given lens model, for example the abundance of highly magnified events, the occurrence of multiple imaging with the images satisfying certain conditions, and the like. This is done in principle by distributing many sources across the source plane, imaging them as described before, and determining the image properties. However, such events are rare. If one were to cover the entire source plane with a regular grid of sources, this grid would have to have a very high resolution for rare events to be reliably found. In turn, most of the sources probed would produce images failing the criteria imposed, so by far the largest fraction of the CPU time used would be wasted.

This situation calls for adaptive grids. We know in advance that any strongly lensed image will occur near a critical curve, or any strongly lensed source near a caustic. It is those sources that we need to treat in detail, while those far from caustic curves are usually only required to normalise the statistics properly.

One approach for defining an adaptive grid, and there may be others more suitable for a particular lensing situation, proceeds as follows. Again, we assume that we know the deflection angle of the lens, either because it was provided numerically or because it is described by a known analytic formula. Then, we saw in the preceding subsection how grid points can easily be identified which are close to a critical curve in the lens plane, or a caustic curve in the source plane.

In order to save computational time, the source plane is first covered with a coarse grid. This grid should obviously be fine enough for the caustics to be properly resolved; for instance, it must not be so coarse that the typically two types of caustic curve, the radial and the tangential one, are closer than a few times the grid separation.

Next, those points on that coarse grid are identified and saved which are next to a caustic curve. This can, for instance, be done by masking, i.e. by attaching a logical variable to all grid points and setting it to either *true* or *false* depending on whether it is or is not next to a critical curve. One can then cover the source plane with a grid whose resolution is doubled in both dimensions, and keep only those points which are identical with, or surrounded by, points of the coarse grid which were masked in the preceding step. This procedure can be repeated as often as desired, i.e. until the finest grid level reaches the ultimately required resolution. Note that it is not the grids and their masks which need to be saved, but only the coordinates of those grid points which are either part of the coarse initial grid, or whose logical mask values are *true*. That way, lists of source positions can be constructed which are to be probed later for the images they give rise to.

Naturally, this can only be a basic recipe which needs to be adapted to the situation at hand. For instance, the condition that grid points need to be next to a caustic can be replaced by the condition that the absolute value of the Jacobian determinant be less than a given threshold which can be lowered at each step of grid refinement. Such a criterion would naturally increase the grid resolution near such grid positions where sources are certain to be highly magnified.

Of course, if statistics is the ultimate goal, one has to take into account that sources near caustic curves were positioned such as to have an unfair advantage over sources far from caustics. Since we have chosen to double the grid resolution at each refinement step, each source on a refined grid represents only a quarter of the area on the source plane represented by a source on the next coarser grid. Assigning a statistical weight of unity to the sources on the finest grid, the weight must quadruple for each coarser level. If the grid was refined $`N`$ times, the weight of sources on the coarsest grid is thus $`w_i=2^{2N}`$. Each source is assigned a statistical weight $`w_i`$ in that way, and counts $`w_i`$ times in the final statistical evaluation.

The left panel in Fig. 7 shows the source locations chosen for evaluating image statistics of the spiral lens model illustrated in Fig. 6.

### 2.7 Finding Images

The principle of finding the images of a given source is simple: Given a source at position $`\stackrel{}{y}_\mathrm{s}`$, find those grid points $`\stackrel{}{x}_{ij}`$ on the lens plane which are mapped sufficiently close to $`\stackrel{}{y}_\mathrm{s}`$, i.e. whose mapped points $`\stackrel{}{y}_{ij}`$ are within a specified distance from $`\stackrel{}{y}_\mathrm{s}`$.

The problem with this approach is that a square-shaped or rectangular grid cell from the image plane is mapped onto a distorted figure in the source plane. In most cases, this figure will be a parallelogram, but in rare cases, opposing corners of the original rectangle may even be interchanged on the source plane. How can it then be decided whether a given point in the source plane is inside or outside the mapped grid cell, or in other words, whether the image of the given source falls within that particular grid cell on the lens plane?

The solution is to split each grid cell in the lens plane into two triangles, because a mapped triangle always remains a triangle, which always has a well-defined interiour (cf. Schneider et al. 1992).

Consider Fig. 8. The three grid points marked on the lens plane in the left panel of the figure are mapped to the distorted triangle shown on the right panel, which contains the source position. Call $`\stackrel{}{d}_{1,2,3}`$ the three vectors from the mapped triangle’s corners towards the source position. It can be shown that the source is inside the mapped triangle if the three vector products

$$\stackrel{}{d}_1\times \stackrel{}{d}_2,\stackrel{}{d}_1\times \stackrel{}{d}_3,\stackrel{}{d}_2\times \stackrel{}{d}_3$$
(14)
are all positive, with the vector product in two dimensions being defined as

$$\stackrel{}{a}\times \stackrel{}{b}a_1b_2a_2b_1.$$
(15)
One straightforward way to verify this condition is to convince one’s self that the source point is inside the triangle if all vectors $`\stackrel{}{d}_i`$ point within the angles spanned by the adjacent sides of the triangle, and that this condition translates to Eq. (14) above.

This algorithm for finding images works well as long as the separation between images is larger than the size of the grid cells in the lens plane. Very close images can be contained within the same grid cell, in which case the algorithm would find only one. Of course, this potential problem can be remedied by increasing the grid resolution, but then a very large number of grid cells would have to be checked in vain for containing an image.

Again, a viable solution uses adaptive grids. One can start with a coarse grid on the lens plane. Searching for images on that coarse grid will almost certainly not yield all images of a multiply imaged source, but those missed will be closer than the grid separation to those found. Then, those grid cells containing images can individually be covered with a highly resolved grid, and the image search repeated on those sub-grids. Hence, the first step represents a coarse scan of the lens plane for grid cells containing at least one image, and the second step scans only those regions on the lens plane in detail where images are sure to be found. If needed, further sub-grids can be similarly nested.

Of course, even though this procedure is highly adaptive and efficient, it always has a remaining resolution limit, and images closer than that will not be resolved. It is then important to adapt the resolution of the finest sub-grid to the situation at hand, for instance such that remaining unresolved images would neither be resolved by observations. The right panel in Fig. 7 illustrates the result of an adaptive image search for all sources at the positions shown in the figure’s left panel. Colours denote image numbers: Black means one image, blue three, and red five, while green shows source positions for which an even number of images has been found, in contradiction to the necessarily odd image number produced by non-singular lenses. Such events are rare, but they do occur because of the finite resolution limit of the algorithm applied.

Figure 9 gives an example for possible results of that adaptive technique for finding images. Colour-coded is the total magnification of point sources in the source plane behind the almost edge-on spiral lensing galaxy introduced above. The increasing spatial resolution towards the caustic curves is evident. The panel inserted into the figure shows caustics (blue) and critical curves (red) of the lens, the source position as a blue dot just inside the right-hand “naked” cusp, and the three images as red hexagons whose size logarithmically encodes the image magnification.

### 2.8 Asymmetric Lenses

So far, we have used a model for a spiral galaxy as an example for a complex lens whose properties need to be determined numerically. Despite its complexity, the model is still highly symmetric; and what is more, its deflection angle is given as an analytic formula. Sources were so far assumed to be point-like.

Let us now increase the level of complexity and use a numerically simulated galaxy cluster to gravitationally lens extended sources. Again, we assume the deflection angle to be given and postpone the question as to how it can be determined from an $`N`$-body simulation.

All techniques described above for computing convergence and shear from the deflection angle, for finding critical curves and caustics, for placing sources on an adaptive grid, and for finding images within grid cells split into triangles remain valid unchanged. Figure 10 shows an example.

The modulus of the cluster’s deflection angle is shown as the colour plot in the left panel. The right panel shows a section of the source plane with the dots marking source positions, and their colour illustrating the image number. Black, blue and red means one, three, or five images, respectively. The caustic structures can clearly be identified as the boundaries between black and blue and between blue and red, respectively.

### 2.9 Imaging Extended Sources

Extended sources can be described in a variety of ways. What follows is a simple description for elliptical sources, but alternative source models can easily be constructed along similar lines. We assume that source positions $`\stackrel{}{y}_\mathrm{s}`$ have already been found, preferentially on an adaptive grid as described before. Also, we need to be sure that the grid resolution in the source plane is sufficiently high as to resolve the smallest sources to be considered.

Elliptical sources are described by three more parameters, viz. their size, their ellipticity, and their position angle $`\varphi `$. Let us describe the ellipticity by $`e=b/a`$, with $`a`$ and $`b`$ being the semi-major and semi-minor axes of the ellipse, respectively. Finally, we introduce an effective radius $`r`$ by demanding that a circle of radius $`r`$ have the same area as the ellipse, hence $`r=\sqrt{ab}`$. By rotating by an angle $`\varphi `$ an ellipse centred on the coordinate origin whose axes are aligned with the coordinate axes, it can straightforwardly be shown that a grid point $`\stackrel{}{y}_{ij}`$ is enclosed by the ellipse if the condition

$`\mathrm{cos}^2\varphi \left({\displaystyle \frac{\delta y_1^2}{e}}+e\delta y_2^2\right)`$ $`+`$ $`\mathrm{sin}^2\varphi \left({\displaystyle \frac{\delta y_2^2}{e}}+e\delta y_1^2\right)`$ (16)
$`+`$ $`2\delta y_1\delta y_2\mathrm{sin}\varphi \mathrm{cos}\varphi \left({\displaystyle \frac{1}{e}}e\right)r^2`$
is satisfied, where $`\delta \stackrel{}{y}\stackrel{}{y}_{ij}\stackrel{}{y}_\mathrm{s}`$. If the grid point $`\stackrel{}{x}_{ij}`$, whose image in the source plane is $`\stackrel{}{y}_{ij}`$, satisfies Eq. (16), the image point $`\stackrel{}{x}_{ij}`$ is part of the source, and the image can be constructed by assigning the source’s surface brightness at $`\stackrel{}{y}_{ij}`$ to the image point $`\stackrel{}{x}_{ij}`$. By mapping the entire lens plane onto the source plane and checking Eq. (16) for each individual imaged grid point $`\stackrel{}{y}_{ij}`$, all image points belonging to the given source can be identified.

It is often desired for statistical purposes to automatically characterise a large number of images. An example is the determination of cross sections for the formation of large gravitational arcs by a numerically simulated galaxy cluster, for which a large number of sources need to be imaged and the image properties automatically quantified to search for the rare “giant” arcs. Most of the methods described here have been introduced and used extensively e.g. by Bartelmann & Weiss (1994), Bartelmann et al. (1995, 1998), Meneghetti et al. (2000, 2001); see also the contribution by Massimo Meneghetti to this volume.

A source may have multiple images, thus the point sets in the lens plane found by imaging extended sources need not be connected. The first step is therefore to group the image points into images. This can be done with a variant of the classical friends-of-friends algorithm: Pick one arbitrary point out of any given set of image points and search for another image point which is at most $`\sqrt{2}h`$ grid units away from the first point; $`h`$ is the grid size in the lens plane. If there is such a point, it is called a “friend” and grouped into the same image as the first point. Now take the “friend” and repeat until no further “friends” can be found and the image is complete. If more image points are left on the lens plane, pick one of those and repeat the process until all image points have been grouped. If the image is large enough, and the grid resolution on the lens plane is high enough for the image to consist of many points, the image magnification is simply the ratio between the numbers of pixels covered by an image and the number of pixels covered by the source.

Once all image points belonging to a single image have been identified, it is often useful to determine the boundary points of that image, e.g. by identifying those points inside an image which have a neighbour outside the image. By suitably ordering the boundary points, a boundary line can be found whose length can be measured and used in further steps of the automatic image classification. Next, the curvature of the image can be found by first identifying the image point which of the source centre, then search for the boundary point most distant from the so-defined image centre, and finally searching for the boundary point most distant from the first boundary point. These three points uniquely define a circle whose radius can be used as an approximation for the arc radius. And so on, you get the drift: Once image points are grouped into individual images and boundary curves have been determined, images can be classified by adapting elementary geometrical figures to them.

### 2.10 Deflection Angles of Asymmetric Lenses

So far we have assumed to be given the deflection angle either as an analytic expression or as two two-dimensional arrays of numbers giving its two components as a function of position in the lens plane. We now need to describe methods for obtaining the deflection angle of a numerically simulated lens.

The first issue to be discussed is the spatial resolution. Since the simulated lens is composed of discrete particles which represent a smooth mass distribution in reality, the deflection angle must not be computed by simply summing up the deflection angles of the individual particles: The result would be a collection of microlenses rather than a single macrolens, having many spurious and undesired imaging features.

Rather, the collection of particles has to be projected onto a lens plane, on which it needs to be smoothed in some way. We will return later to the issue of how particles should be sorted into grid cells. An important point to be addressed before is how large the grid cells should be chosen. They should be small enough for important features of the lens to remain identifiable; they should be large enough for the surface density to lose the “graininess” due to its being composed of individual particles, and they should be large enough so that Poisson errors are smaller than a certain threshold. If the number of particles per grid cell is $`nh^2`$, its Poisson fluctuation is $`\sqrt{nh^2}`$, thus the discreteness of the particles gives rise to fluctuations in the surface-mass density. Demanding that the relative fluctuations of the density should be smaller than $`ϵ1`$, the cell size $`h`$ has to be chosen such as to satisfy $`(nh^2)^{1/2}ϵ`$. It is impossible to give a general rule applicable to the majority of lensing situations, but it is clear that resolution, smoothing and particle noise have to be carefully balanced by choosing the grid cell size appropriately.

Assigning particle masses to grid points in order to obtain a smooth density distribution is an art of its own (cf. Hockney & Eastwood 1988). In principle, the particle mass could simply be attributed to the single grid point next to its position. This “nearest grid point” (NGP) method is appropriate for particles near the centre of a cell, but particles near cell boundaries should be attributed to the cell and its neighbour(s) in order to avoid boundary effects like density discontinuities. Numerous schemes for interpolating particles across cells have been proposed. They are generally of the form

$$Q(\stackrel{}{x})=\underset{i}{}W\left(\stackrel{}{x}\stackrel{}{x}_i\right)Q(\stackrel{}{x}_i),$$
(17)
where $`Q`$ is the quantity to be interpolated onto a point $`\stackrel{}{x}`$, e.g. the particle mass, the sum extends over all particles sufficiently close to the point of interest $`\stackrel{}{x}`$, and $`W(\stackrel{}{x}\stackrel{}{x}_i)`$ is a smoothing or interpolation kernel depending on the separation vector between the particle position $`\stackrel{}{x}_i`$ and $`\stackrel{}{x}`$. The kernel is decomposed into three factors directions,

$$W(\delta \stackrel{}{x})=w(\delta x_1)w(\delta x_2)w(\delta x_3),$$
(18)
one for each dimension, the $`i`$-th of which depends only on the $`i`$-component of the separation vector. Interpolation methods can now be classified according to the kernel factors $`w(\delta x)`$ and their width.

The “cloud-in-cell” (CIC) scheme uses the kernel factors

$$w_{\mathrm{CIC}}(\delta x)=\{\begin{array}{cc}1|\delta x|/h\hfill & \text{for}|\delta x|<h\hfill \\ 0\hfill & \text{otherwise}\hfill \end{array},$$
(19)
which implies that the particle is distributed over the four nearest grid points. A more elaborate scheme is the “triangular shaped cloud” (TSC) method, which uses the kernel factors

$$w_{\mathrm{TSC}}(x)=\{\begin{array}{cc}3/4\delta x^2/h^2\hfill & \text{for}|\delta x|h/2\hfill \\ (3/2|\delta x|/h)^2/2\hfill & \text{for}h/2|\delta x|<3h/2\hfill \\ 0\hfill & \text{otherwise}\hfill \end{array}.$$
(20)
The CIC and TSC interpolation schemes are illustrated for two dimensions in Fig. 11. For all schemes, the kernel has to be normalised such that all particle mass fractions add up to unity.

Suppose now we have obtained the surface mass density on a grid $`\kappa _{ij}=\kappa (\stackrel{}{x}_{ij})`$, then the deflection angle can most straightforwardly be determined by direct summation as

$$\stackrel{}{\alpha }_{ij}=\frac{1}{\pi }\underset{kl}{}\kappa _{kl}\frac{\stackrel{}{x}_{ij}\stackrel{}{x}_{kl}}{\left|\stackrel{}{x}_{ij}\stackrel{}{x}_{kl}\right|^2}.$$
(21)
Depending on the number of grid cells, the direct summation can be prohibitively slow. In many circumstances of astrophysical interest, fast-Fourier techniques can then be applied. In order to see how this works, note that the deflection angle can be written as a convolution of the convergence $`\kappa (\stackrel{}{x})`$ with a kernel

$$\stackrel{}{K}(\stackrel{}{x})=\frac{1}{\pi }\frac{\stackrel{}{x}}{\left|\stackrel{}{x}\right|^2}.$$
(22)
This allows the Fourier convolution theorem to be applied, which holds that the Fourier transform of a convolution is the product of the Fourier transforms of the functions to be convolved, hence

$$\widehat{\stackrel{}{\alpha }}(\stackrel{}{k})=\widehat{\kappa }(\stackrel{}{k})\widehat{\stackrel{}{K}}(\stackrel{}{k}).$$
(23)
The Fourier transform of the kernel $`\stackrel{}{K}`$ can be determined and tabulated once. Using fast-Fourier techniques to determine the Fourier transform of the convergence $`\widehat{\kappa }(\stackrel{}{k})`$ requires the convergence to be periodic on the lens plane. In many cases, this can be safely assumed or arranged. Often, lens planes are constructed from large-scale $`N`$-body simulations which have periodic boundary conditions by design, or the lens is an isolated object like a galaxy cluster, which can be surrounded by a sufficiently large field for the convergence to drop near zero everywhere around the edges of the field. Fast-Fourier methods speed up the computation of the deflection angle considerably.

If necessary, derivatives of the deflection angle field can also be determined in Fourier space. Once the convergence has been Fourier transformed, one can employ the two-dimensional Poisson equation to compute the Fourier transform of the lensing potential,

$$\widehat{\psi }=\frac{2}{k^2}\widehat{\kappa },$$
(24)
from which the Fourier transforms of the deflection angle and the shear components can easily be determined,

$$\widehat{\stackrel{}{\alpha }}=\mathrm{i}\stackrel{}{k}\widehat{\psi },\widehat{\gamma }_1=\frac{1}{2}\left(k_1^2k_2^2\right)\widehat{\psi },\widehat{\gamma }_2=k_1k_2\widehat{\psi }.$$
(25)
Relations like those and the exploitation of fast-Fourier methods are particularly relevant for simulating gravitational lensing by large-scale structures.

## 3 Lensing by Large-Scale Structures

### 3.1 Resolution Issues

Obviously, the thin-lens approximation that we have been using so far breaks down if one wishes to study gravitational lensing by large-scale structures. The solution then is to cover the complete cosmic volume whose lensing effects one wants to simulate with simulation boxes stacked along the line-of-sight, to project suitable slices on individual lens planes, and to use multiple lens-plane theory for describing light propagation.

The multiplicity of lens planes, and the general weakness of lensing by large-scale structures, make questions of angular and mass resolution particularly relevant for cosmic lensing. For instance, lens planes close to the observer are typically poorly resolved because even small grid cells span a large solid angle near the observer, and making grid cells smaller is not generally an acceptable solution because then the number of particles per grid cell becomes small, and the shot noise possibly unacceptably large. However, lens planes near the observer are less efficient than lens planes approximately half-way to the source because the lensing efficiency function is zero at the observer and source redshifts and peaks in between. Yet, structures grow over time, thus the lensing efficiency function is skewed towards lower redshifts because structures are geometrically less efficient lenses, but their density contrast keeps growing. By a related argument, sources at very high redshifts do not require the entire cosmological volume between them and the observer to be filled with lens planes because lens planes at very high redshift are geometrically inefficient and have a low density contrast. The left panel of Fig. 12 shows two examples for the lensing efficiency function times the linear growth factor, which is the relevant quantity combining structure growth with geometrical efficiency.

Similarly, the effective angular resolution of the simulation is dominated by the angular resolution of those lens planes near the peak in the combined efficiency function, i.e. the product of geometrical efficiency and linear growth factor.

The shot noise caused by the discretisation of mass into particles is particularly important for studies of weak lensing by large-scale structures. Even in absence of density inhomogeneities, shot noise leads to density fluctuations. They need to be sufficiently smaller than the signal, i.e. the convergence fluctuations which cause weak lensing.

In essence, this requirement also imposes a resolution limit. Suppose we wish to quantify the weak-lensing signal within a solid angle $`\delta \mathrm{\Omega }`$. The volume spanned by $`\delta \mathrm{\Omega }`$ within redshifts $`z`$ and $`z+\mathrm{d}z`$ is

$$\mathrm{d}V(z)=\delta \mathrm{\Omega }D^2(z)\left|\frac{\mathrm{d}D_{\mathrm{prop}}}{\mathrm{d}z}\right|\mathrm{d}z,$$
(26)
where $`D(z)`$ and $`D_{\mathrm{prop}}(z)`$ are the angular diameter and proper distances to redshift $`z`$. In absence of density inhomogeneities, this volume element contains $`\mathrm{d}N(z)`$ particles, with

$$\mathrm{d}N(z)=\mathrm{d}V(z)\frac{\overline{\rho }(z)}{m_\mathrm{p}},$$
(27)
where $`\overline{\rho }(z)`$ is the mean matter density at redshift $`z`$, and $`m_\mathrm{p}`$ is the mass of an $`N`$-body particle in the simulation. The contribution to the lensing convergence by these particles has to be weighted by the effective lensing distance, $`D_{\mathrm{eff}}(z,z_\mathrm{s})`$, and by numerical factors. Poisson fluctuations in the particle number thus cause convergence fluctuations whose variance is

$$\delta ^2\kappa _0^{z_\mathrm{s}}dzD_{\mathrm{eff}}^2(z,z_\mathrm{s})dN(z).$$
(28)
These fluctuations need to be compared with, and smaller than, the convergence fluctuations due to large-scale structure, which are typically of order $`\kappa ^2^{1/2}5\%`$ for sources near redshift unity and angular scales of order $`1^{}`$. According to Eqs. (26) through (28), the *rms* shot noise $`\delta ^2\kappa ^{1/2}`$ scales like $`\delta \mathrm{\Omega }^{1/2}`$, thus the requirement that the signal-to-noise ratio

$$\frac{\mathrm{S}}{\mathrm{N}}=\left(\frac{\kappa ^2}{\delta ^2\kappa }\right)^{1/2}$$
(29)
exceed a specified threshold translates into a lower limit to the solid angle $`\delta \mathrm{\Omega }`$ which can reasonably be resolved by the simulation. The smallness of the *rms* cosmic convergence $`\kappa _{\mathrm{rms}}=\kappa ^2^{1/2}`$ implies that many particles need to be enclosed by the “cone” spanned by $`\delta \mathrm{\Omega }`$ for the simulation to be reliable. The right panel of Fig. 12 shows an example. The *rms* cosmic convergence in per cent and the *noise-to-signal* ratio are plotted as functions of angular scale. The noise level was adapted to an $`N`$-body simulation with particle mass $`m_\mathrm{p}=6.8\times 10^{10}h^1M_{}`$. The curves show that the noise-to-signal ratio drops below unity for sources at redshift $`z_\mathrm{s}=1`$ only if the angular resolution is lowered to $`5^{}`$, while an angular resolution of $`0.8^{}`$ can be achieved for $`z_\mathrm{s}=1000`$ (i.e. for weak lensing of the CMB; Pfrommer 2002).

### 3.2 Multiple Lens-Plane Theory

Weak lensing by large-scale structures requires the cosmic volume to be split into multiple lens planes rather than a single one (for general reference on multiple lens-plane theory, see Schneider et al. 1992). The lens plane closest to the observer is the image plane which represents the observer’s sky. A light ray piercing the image plane at a physical coordinate $`\stackrel{}{\xi }_1`$ is mutiply deflected on $`N`$ lens planes and finally reaches the source plane at the physical coordinate

$$\stackrel{}{\eta }(\stackrel{}{\xi }_1)=\frac{D_\mathrm{s}}{D_1}\stackrel{}{\xi }_1+\underset{i=1}{\overset{N}{}}D_{i\mathrm{s}}\stackrel{}{\widehat{\alpha }}(\stackrel{}{\xi }_i),$$
(30)
where the $`D_i`$ and $`D_{i\mathrm{s}}`$ are the angular diameter distances from the observer to the $`i`$-the lens plane, and from the $`i`$-th lens plane to the source, respectively. The light ray passes the $`i`$-th plane at $`\stackrel{}{\xi }_i`$, where it is deflected by $`\stackrel{}{\widehat{\alpha }}(\stackrel{}{\xi }_i)`$. Similarly, the $`\stackrel{}{\xi }_i`$ are determined by

$$\stackrel{}{\xi }_j(\stackrel{}{\xi }_1)=\frac{D_j}{D_1}\stackrel{}{\xi }_1+\underset{i=1}{\overset{j1}{}}D_{ij}\stackrel{}{\widehat{\alpha }}(\stackrel{}{\xi }_i),$$
(31)
where $`D_{ij}`$ is the angular diameter distance from the $`i`$-th to the $`j`$-th lens plane.

Introducing angular coordinates $`\stackrel{}{\theta }_i=\stackrel{}{\xi }_i/D_i`$ yields

$$\stackrel{}{\theta }_j(\stackrel{}{\theta }_1)=\stackrel{}{\theta }_1+\underset{i=1}{\overset{j1}{}}\frac{D_{ij}D_\mathrm{s}}{D_jD_{\mathrm{is}}}\stackrel{}{\alpha }(\stackrel{}{\theta }_i),$$
(32)
where we have introduced the *reduced* deflection angle $`\stackrel{}{\alpha }=(D_{i\mathrm{s}}/D_\mathrm{s})\widehat{\stackrel{}{\alpha }}`$. We now define the matrices

$$𝒜_i=\frac{\stackrel{}{\theta }_i}{\stackrel{}{\theta }_1},𝒰_i=\frac{\stackrel{}{\alpha }_i}{\stackrel{}{\theta }_i}.$$
(33)
Clearly, $`𝒜_i`$ is the Jacobian matrix of the lens mapping between the $`i`$-th lens plane and the image plane, thus $`𝒜_N`$ is the Jacobian matrix of the mapping between the source and image planes. The goal is thus to determine $`𝒜_N`$ in order to obtain convergence, shear, and magnification for a light ray starting out into direction $`\stackrel{}{\theta }_1`$. The ray-tracing equation (32) implies the recursion relation

$$𝒜_j=\underset{i=1}{\overset{j1}{}}\frac{D_{ij}D_\mathrm{s}}{D_jD_{\mathrm{is}}}𝒰_i𝒜_i,$$
(34)
starting with $`𝒜_1=`$, the identity matrix. In summary, the deflection-angle fields $`\stackrel{}{\alpha }_i`$ on the $`N`$ lens planes can be used to construct the matrices $`𝒰_i`$ according to Eq. (33), then Eq. (34) can be used to determine the lensing experienced by a light ray starting out into any direction on the image plane. The left panel in Fig. 13 shows the total convergence experienced by sources at $`z_\mathrm{s}=5`$ on a lens plane with a side length of $`4.25^{}`$, obtained from an $`N`$-body simulation (Pfrommer 2002).

The right panel in Fig. 13 shows numerically determined power spectra for the effective convergence as functions of wave number $`l`$, which is the Fourier conjugate variable to the angular scale. The lines in this figure show the theoretically expected power spectra. The agreement between the numerical and theoretical results is very good over a limited range of wave numbers. Once the wave numbers increase beyond the limit set by the angular resolution, the simulated convergence fields lack power and the numerical results fall below the theoretical ones. This happens at lower $`l`$ for smaller source redshifts, because a fixed angular scale, and thus wave number $`l`$, corresponds to smaller physical scales at lower distances. On the low-$`l`$ end, i.e. for large structures, the errors on the numerically determined power spectra increase because the number of independent modes in the simulated convergence field decreases as the modes increase. This example should suffice to demonstrate that numerical simulations of gravitational lensing by large-scale structures should be carefully designed to match their final purpose.

## 4 Inversion Techniques

Let us conclude with a brief discussion of inversion techniques. They are typically less demanding numerically, but the methods which have been developed for this purpose are interesting in their own right.

### 4.1 Shear Deconvolution

We have seen before in Eqs. (24) and (25) that convergence and shear are related because they are both linear combinations of second derivatives of the scalar lensing potential $`\psi `$. In Fourier space, the relations are algebraic and can easily be combined to eliminate the Fourier transform $`\widehat{\psi }`$ of the potential. Transforming back into configuration space, the convergence turns out to be a convolution of the shear with a well-known kernel,

$$\kappa (\stackrel{}{\theta })=\frac{1}{\pi }\mathrm{d}^2\theta ^{}\left[𝒟_1(\stackrel{}{\theta }\stackrel{}{\theta }^{})\gamma _1(\stackrel{}{\theta }^{})+𝒟_2(\stackrel{}{\theta }\stackrel{}{\theta }^{})\gamma _2(\stackrel{}{\theta }^{})\right],$$
(35)
with

$$𝒟_1(\stackrel{}{\theta })=\frac{\theta _1^2\theta _2^2}{|\stackrel{}{\theta }|^4},𝒟_2(\stackrel{}{\theta })=\frac{2\theta _1\theta _2}{|\stackrel{}{\theta }|^4}.$$
(36)
This is the classic Kaiser & Squires (1993) shear inversion equation. Its limitations have been discussed in detail and removed to satisfaction by modifying e.g. the kernel components $`𝒟_{ij}`$; they are not of interest for the discussion here (cf. Peter Schneider’s lecture in this volume).

A suitable practical approximation of (35) using measured galaxy ellipticities $`ϵ_i`$ ($`i=1,2`$) is

$$\kappa (\theta )\frac{1}{n\pi }\underset{i=1}{\overset{N}{}}\left[𝒟_1ϵ_{1,i}+𝒟_2ϵ_{2,i}\right],$$
(37)
where $`n`$ is the number density of lensed galaxies on the sky. In practice, however, it turns out that an approximation like (37) would have infinite noise because of the random sampling of the shear components $`\gamma _i`$ by $`N`$ galaxy ellipticities $`ϵ_i`$. This can be remedied by introducing a smoothed kernel $`𝒟^{}`$ instead of $`𝒟`$, e.g.

$$𝒟^{}=\left[1\left(1+\frac{\theta ^2}{\theta _\mathrm{s}^2}\right)\mathrm{exp}\left(\frac{\theta ^2}{\theta _\mathrm{s}^2}\right)\right]𝒟,$$
(38)
where $`\theta _\mathrm{s}`$ is the angular smoothing scale (Seitz & Schneider 1995). The noise convariance matrix between the convergence values at two different grid points $`\stackrel{}{\theta }_i`$ and $`\stackrel{}{\theta }_j`$ is then

$$\kappa (\stackrel{}{\theta }_i)\kappa (\stackrel{}{\theta }_j)=\frac{\sigma _ϵ^2}{4\pi \theta _\mathrm{s}^2n}\mathrm{exp}\left[\frac{(\stackrel{}{\theta }_i\stackrel{}{\theta }_j)^2}{2\theta _\mathrm{s}^2}\right],$$
(39)
where $`\sigma _ϵ`$ is the scatter of the intrinsic galaxy ellipticities (van Waerbeke 2000). This expression demonstrates that smoothing introduces correlations on the convergence map on the angular scale $`\theta _\mathrm{s}`$, but the variance of $`\kappa `$ can become very high if $`\theta _\mathrm{s}`$ is chosen too small. A careful balance between the local variance and non-local correlations is necessary in order to arrive at a convergence map with the required properties.

### 4.2 Maximum-Likelihood Lens Inversion

An entirely different approach to lens inversion uses the maximum-likelihood technique (Bartelmann et al. 1996). Each lensed background galaxy $`i`$ provides a measurement of two ellipticity components $`(ϵ_{1,i},ϵ_{2,i})`$ and its angular size. Comparing the size of a galaxy behind a galaxy cluster to the average size of unlensed galaxies of the same surface brightness, an estimate $`r_i`$ of the inverse magnification of the lensed galaxy can be derived. Thus $`N`$ galaxies provide a $`3N`$-dimensional data vector

$$\stackrel{}{d}=(ϵ_{1,1},ϵ_{2,1},r_1,\mathrm{},ϵ_{1,N},ϵ_{2,N},r_N).$$
(40)
The goal of the lens inversion is then to find a two-dimensional array $`\psi _{jk}`$ of lensing potential values such that the ellipticities and inverse magnifications caused by that potential at the positions $`\stackrel{}{\theta }_i`$ of the real galaxies optimally reproduce the measured ellipticities and inverse magnifications. In other words, the potential values $`\psi _{jk}`$ have to be determined such as to minimise the mean-square difference between the data vector $`\stackrel{}{d}`$ and the model data vector $`\stackrel{}{d}[\psi _{jk}(\stackrel{}{x}_i)]`$,

$$\chi ^2(\psi _{jk})=\underset{i=1}{\overset{3N}{}}\left\{\frac{[d_id_i(\psi _{jk})]^2}{\sigma _i^2}\right\},$$
(41)
where the errors $`\sigma _i`$ can be estimated from the data themselves. The minimisation of $`\chi ^2`$ with respect to the potential values $`\psi _{jk}`$ can be done with any minimisation algorithm like, e.g. the downhill simplex method. For large fields, the number of potential values can become very large. In that case, conjugate-gradient methods are preferred, which make use of the fact that the derivatives of $`\chi ^2`$ with respect to the $`\psi _{jk}`$ are known analytically. Such methods can speed up the minimisation sufficiently to render it feasible even for large potential arrays (cf. Press et al. 1992).

### 4.3 Maximum-Entropy Methods

The minimisation of $`\chi ^2`$ is a special case of the maximum-likelihood technique for assumed Gaussian deviations of the measured data around the model values. Improvements of the maximum-likelihood technique can be derived starting from Bayes’ theorem,

$$P(\psi |\stackrel{}{d})=\frac{P(\stackrel{}{d}|\psi )}{P(\stackrel{}{d})}P(\psi ),$$
(42)
which states that the probability $`P(\psi |\stackrel{}{d})`$ of finding the potential $`\psi `$ given the data $`\stackrel{}{d}`$ is proportional to the probability $`P(\stackrel{}{d}|\psi )`$ of obtaining the data given the potential, times the probability $`P(\psi )`$ for finding the potential. The denominator $`P(\stackrel{}{d})`$ is called the *evidence* and simply normalises Eq. (42). $`P(\psi )`$ is called the prior, quantifying any *a priori* information one has or assumes on the potential $`\psi `$, $`P(\stackrel{}{d}|\psi )`$ is called the likelihood, and $`P(\psi |\stackrel{}{d})`$ is the *posterior* probability. The goal is now to maximise the latter, which is equivalent to maximising the product $`P(\stackrel{}{d}|\psi )P(\psi )`$ of likelihood and prior. If we have or can assume Gaussian noise and a diagonal noise correlation matrix, the likelihood reduces to $`P(\stackrel{}{d}|\psi )=\mathrm{exp}(\chi ^2/2)`$.

It can now be shown that in absence of any further information, the best, i.e. least prejudiced, prior is the maximum-entropy prior,

$$P(\psi )\mathrm{exp}\left[\alpha S(\psi ,\stackrel{}{m})\right],$$
(43)
with the *cross entropy*

$$S(\psi ,\stackrel{}{m})=\underset{i=1}{\overset{3N}{}}\psi _im_i\psi _i\mathrm{ln}\frac{\psi _i}{m_i},$$
(44)
where $`\stackrel{}{m}`$ is a model vector for the potential which can encode expectations on the potential, or simply be chosen to be uniform for all $`i`$. The potential array is then determined by maximising $`\mathrm{exp}(\chi ^2/2+\alpha S)`$, or equivalently by minimising

$$F\frac{1}{2}\chi ^2\alpha S$$
(45)
instead of the simple $`\chi ^2`$ in Eq. (41). The parameter $`\alpha `$ can be included into the minimisation. Bayesian theory implies that a good approximation to the optimal choice for $`\alpha `$ is determined such that $`F3N/2`$ at the potential minimum $`\overline{\psi }`$. The error covariance matrix for the potential $`\psi `$ is given by the inverse curvature matrix of $`F`$,

$$(\psi \overline{\psi })(\psi \overline{\psi })^\mathrm{T}\left(\frac{^2F}{\psi _i\psi _j}\right)^1.$$
(46)
Maximum-entropy methods have been suggested and used for regularising shear-inversion techniques such that their spatial resolution is adapted to the strength of the lensing signal (Bridle et al. 1998; Seitz et al. 1998).

## 5 Concluding Remarks

Many numerical methods have been used for gravitational lensing studies which I was not able to cover during the limited time of the lecture. Among them are the hierarchical tree-code methods introduced into microlensing by Wambsganss et al. (1990) and the methods for constraining cluster mass distributions from multiple arc systems (e.g. Kneib et al. 1993; see also Jean-Paul Kneib’s presentation in this volume). Despite this unavoidable incompleteness, I hope to have given a flavour of how numerical methods can be used for lensing, and what the main problem areas are.

###### Acknowledgements.

I wish to thank the organisers of the lensing school in Aussois which provided a nice opportunity to discuss a broad variety of lensing topics, and a scenic frame.


format: text

text:
# Mixture of multiple copies of maximally entangled states is quasi-pure

## Abstract

Employing the general BXOR operation and local state discrimination, the mixed state of the form

$$\rho _d^{(k)}=\frac{1}{d^2}\underset{m,n=0}{\overset{d1}{}}(|\varphi _{mn}\varphi _{mn}|)^k$$

is proved to be quasi-pure, where $`\{|\varphi _{mn}\}`$ is the canonical set of mutually orthogonal maximally entangled states in $`d\times d`$. Therefore irreversibility does not occur in the process of distillation for this family of states. Also, the distillable entanglement is calculated explicitly.

Entanglement distillation plays a crucial role in quantum information theory Bennett1 . As well-known, entanglement is responsible for many quantum tasks and pure entangled states are required in most cases. Unfortunately, due to the interaction with the environment, pure entanglement is fragile and easy to be blurred by the noise. So distillation of entanglement is of importance. Though many distillation protocols and upper bounds are known, distillable entanglement are known in few nontrivial cases Rains . Actually, entanglement distillation is closely related to local state discrimination. Two counter-intuitive facts are found about local state discrimination. One is that there exist product orthogonal states that could not be discriminated exactly by LOCC operations Bennett . The other is that any two orthogonal multipartite states could be discriminated with certainty by only LOCC operations Walgate1 . Based on the results of Walgate1 , the distillable entanglement Bennett1  is calculated for a class of mixed states comprising of Bell basis Ghosh ; Yang . For generic bipartite mixed states, entanglement distillation is an irreversible process in the sense that more pure-state entanglement is required to create it than can be distilled from it in the asymptotic limit Vidal . However, the process of distillation is reversible in the so-called quasi-pure mixed state Vollbrecht . In this paper, the mixture of multiple copies of maximally entangled states with equal probability Hamieh  is proved to be quasi-pure. So the distillation process is reversible that is not so obvious at first glance. Furthermore, the distillable entanglement is calculated explicitly.

We consider the mixed state of the form,

$$\rho _d^{(k)}=\frac{1}{d^2}\underset{m,n=0}{\overset{d1}{}}(|\varphi _{mn}\varphi _{mn}|)^k,$$
(1)
where $`\{|\varphi _{mn}\}`$ is the canonical set of mutually orthogonal maximally entangled states(MES) in $`d\times d`$ defined as,

$$|\varphi _{mn}=\frac{1}{\sqrt{d}}\underset{j=0}{\overset{d1}{}}e^{2\pi ijn/d}|j|jm,n,m=0,1,\mathrm{},d1.$$
(2)
where $``$ means addition modulo $`d`$.

The situation occurs where Alice and Bob share $`k`$ pairs of qudits, each pair in the same maximally entangled state. In order to share entanglement, they should store the classical information in classical memory that identifies in which state the $`k`$ pairs of qudits are. Now imagine that by some misfortune, the data in the classical memories are deleted, how many ebits can be distilled? As well-known, it is difficult to calculate the distillable entanglement of the mixed states. The main result of the paper is that the mixed state of $`Eq.(1)`$ is quasi-pure Vollbrecht . For quasi-pure states, irreversibility does not occur in the entanglement distillation. And the distillable entanglement is easy to calculated. Before the proof of the main result, let us explain what is the meaning of ”quasi-pure”.

A mixed state is quasi-pure if it can be reversibly transformed to the following form under the local operation and classical communication (LOCC),

$$\rho =p_i\rho _i|\varphi _i\varphi _i|,$$
(3)
where $`\rho _i`$ are separable and orthogonal states that can be discriminated exactly by LOCC operations. The physical meaning of the quasi-pure state is that every pure entangled state is identified by a separable ”tag state” that can be distinguished by LOCC operations. As the transformation is reversible, all the entanglement in the original mixed state is contained in the pure states and can be distilled. As Bennett et. al Bennett  showed that even separable orthogonal states might contain nonlocal property that can not be identified by LOCC operations, it is necessary to introduce the ”tag state” distinguishable by LOCC operations. Similarly, we can define separable-quasi-pure and PPT-quasi-pure with respect to the operations allowed. For instance, consider the constructed state $`\rho =1/9|\psi _i\psi _i||\varphi _i\varphi _i|`$, where $`|\psi _i`$ is the separable pure orthogonal basis in $`3\times 3`$ in Bennett  and $`|\varphi _i`$ is orthogonal MES in $`3\times 3`$. Explicitly, this state is separable-quasi-pure. And the distillation is a reversible process if separable operations are allowed. Is it also LOCC-quasi-pure? Does irreversibility occur in this mixed state? If it does, we can conclude that the irreversibility originates from the nonlocal property of the separable basis. Though we conjecture it is the case, the strict proof escape from us. However, we will show that the state of $`Eq.(1)`$ is quasi-pure though it looks unlike.

Main Result: The mixed state of maximally entangled states in $`Eq.(1)`$ is quasi-pure.

First we introduce the operations that are employed. The control-not operation $`C`$ is defined as $`C|i|j=|i|ji`$, and the bilateral control-not operation (BXOR) $`B`$ is $`B|i_{A1}|r_{B1}|j_{A2}|s_{B2}=|i_{A1}|r_{B1}|ji_{A2}|sr_{B2}`$, where $`\{|i\}`$ is the computational basis set. Denote $`B(m,n)`$ as the BXOR operation performed on the mth pair (source) and the nth pair (target), and $``$ act on $`k`$ copies of $`\varphi _{mn}`$ as

$`\varphi _{mn}^k=B(1,k)B(2,k)\mathrm{}B(k1,k)\varphi _{mn}^k=\varphi _{m,0}^{(k1)}\varphi _{km,n}.`$ (4)
For brevity, we denote the projector $`|\varphi _{mn}\varphi _{mn}|`$ as $`P_{mn}`$ and omit the normalization for the state unless confusion occurs. Three important properties about $`|\varphi _{mn}`$ are utilized. First, $`_mP_{mn}`$ and $`_nP_{mn}`$ are separable. For example, the subspace spanned by $`\{\varphi _{0n}\}`$ can also be spanned by $`\{|ii,i=0,\mathrm{},d1\}`$. For the subspace spanned by $`\{\varphi _{m0}\}`$, it is not so obvious, but can be proved to be separable in the dual basis $`\{|e_k^A\}`$, $`\{|e_k^B\}`$, where $`|k^A=1/\sqrt{d}_{j=0}^{d1}e^{2\pi ijk/d}|e_j^A`$ and $`|e_k^B=1/\sqrt{d}_{j=0}^{d1}e^{2\pi ijk/d}|j^B`$. Second, $`_mP_{mn}`$ are orthogonal to each other and can be distinguished exactly by LOCC. So do $`_nP_{mn}`$. Third, for fixed $`n`$ and varying $`m`$, $`|\varphi _{mn}`$ can be transformed to $`|\varphi _{0m}`$ by local basis transformation.

Proof:

$$\rho _d^{(k)}\stackrel{}{}\underset{m}{}P_{m0}^{(k1)}\underset{n}{}P_{km,n}.$$
(5)
Denote the maximal common factor between $`k`$ and $`d`$ as $`g(k,d)`$. When $`g(k,d)=1`$, We know that $`\{km,misvarying\}`$ is a permutation of the numbers $`\{0,1,\mathrm{},d1\}`$. The $`d`$ mixed states $`_nP_{km,n}`$ are separable and orthogonal to each other. So $`\rho _d^{(k)}`$ is quasi-pure. When $`g(k,d)>1`$, the problem is a little complicated. Suppose $`k=g\stackrel{~}{k},d=g\stackrel{~}{d}`$, where $`\stackrel{~}{k},\stackrel{~}{d}`$ are integers and satisfy $`g(\stackrel{~}{k},\stackrel{~}{d})=1`$. For $`m=0,1,\mathrm{},d1`$, it can be written as $`m=s\stackrel{~}{d}+t,s=0,1,\mathrm{},g1,t=0,1,\mathrm{},\stackrel{~}{d}1`$. Then

$$_d(km)=_{(g\stackrel{~}{d})}g\stackrel{~}{k}(s\stackrel{~}{d}+t)=_{(g\stackrel{~}{d})}g\stackrel{~}{k}t=g_{\stackrel{~}{d}}\stackrel{~}{k}t.$$
(6)
Because $`g(\stackrel{~}{k},\stackrel{~}{d})=1`$, $`\{_{\stackrel{~}{d}}\stackrel{~}{k}t,t=0,1,\mathrm{},\stackrel{~}{d}1\}`$ is a permutation of the numbers $`\{0,1,\mathrm{},\stackrel{~}{d}1\}`$. For brevity, set $`T(t)=g_{\stackrel{~}{d}}\stackrel{~}{k}t`$ that is different for $`t=0,1,\mathrm{},\stackrel{~}{d}1`$. So

$$\underset{m}{}P_{m0}^{(k1)}\underset{n}{}P_{km,n}=\underset{t}{}\{\underset{s}{}P_{s\stackrel{~}{d}+t,0}^{k1}\}\{\underset{n}{}P_{T(t),n}\}.$$
(7)
For every $`t`$, $`_nP_{T(t),n}`$ is separable and orthogonal and can be discriminated by local operations. Now, we show that $`_sP_{s\stackrel{~}{d},0}^{(k1)}`$ is quasi-pure since for the other cases $`t0`$, $`_sP_{s\stackrel{~}{d}+t,0}^{(k1)}`$ can be transformed to $`_sP_{s\stackrel{~}{d},0}^{(k1)}`$ by the permutation of the basis states.

$$P_{s\stackrel{~}{d},0}(d)=P_{s,0}(g)P_{00}(\stackrel{~}{d}),$$
(8)
where $`P_{s,0}(g)`$ represents the projector of the MES in $`g\times g`$. So

$$\underset{s}{}P_{s\stackrel{~}{d},0}^{(k1)}=\underset{s}{}P_{s,0}^{(k1)}(g)P_{00}^{(k1)}(\stackrel{~}{d}).$$
(9)
Under the dual basis, $`_sP_{s,0}^{(k1)}(g)`$ can be written as $`_sP_{0,s}^{(k1)}(g)`$. Further employing the operation $``$

$$\underset{s}{}P_{0,s}^{(k1)}(g)\stackrel{}{}P_{0,s}^{(k2)}(g)\underset{s}{}P_{0,s}(g)$$
(10)
So $`_sP_{s\stackrel{~}{d},0}^{(k1)}`$ is quasi-pure. As all the operations are local reversible operations, the original state is quasi-pure. Though measurements are performed during the process, they are non-demolation measurements and do not destroy the state.

In the following, we give some examples and calculate the distillable entanglement $`E`$ that is equal to the entanglement cost as the mixed state is quasi-pure.

Ex1 $`d=2`$

$`\rho _2^{(2)}`$ $`\stackrel{}{}(P_{00}+P_{10})(P_{00}+P_{01}),`$ (11a)
$`\rho _2^{(3)}`$ $`\stackrel{}{}P_{00}^2(P_{00}+P_{01})+P_{10}^2(P_{10}+P_{11}),`$ (11b)
$`\rho _2^{(4)}`$ $`\stackrel{}{}(P_{00}^3+P_{10}^3)(P_{00}+P_{01})`$ (11c)
Because $`P_{00}+P_{10}`$, $`P_{00}+P_{01}`$ are separable, $`\rho _2^{(2)}`$ is separable. The same conclusion is obtained through direct calculation Smolin . $`E(\rho _2^{(2)})=0`$. Similarly it is easy to show $`\rho _2^{(3)}`$ is quasi-pure and $`E(\rho _2^{(3)})=2`$. It is required to explain for $`\rho _2^{(4)}`$. $`P_{00}^3+P_{10}^3`$ can be transformed to $`P_{00}^3+P_{01}^3`$ by local unitary operation. Performing $``$ again, $`P_{00}^3+P_{01}^3\stackrel{}{}P_{00}^2(P_{00}+P_{01})`$ that is quasi-pure and $`E(\rho _2^{(4)})=2`$. For larger $`k`$, the cases above occur similarly.

Ex2 $`d=3`$

$`\rho _3^{(2)}`$ $`\stackrel{}{}P_{00}{\displaystyle \underset{n}{}}P_{0n}+P_{10}{\displaystyle \underset{n}{}}P_{2n}+P_{20}{\displaystyle \underset{n}{}}P_{1n},`$ (12a)
$`\rho _3^{(3)}`$ $`\stackrel{}{}{\displaystyle \underset{m}{}}P_{m0}^2{\displaystyle \underset{n}{}}P_{0n},`$ (12b)
$`\rho _3^{(4)}`$ $`\stackrel{}{}P_{00}^3{\displaystyle \underset{n}{}}P_{0n}+P_{10}^3{\displaystyle \underset{n}{}}P_{1n}+P_{20}^3{\displaystyle \underset{n}{}}P_{2n}`$ (12c)
It is required to show $`_mP_{m0}^2`$ is quasi-pure. $`\varphi _{m0}`$ can be transformed to $`\varphi _{0m}`$ in the dual basis. Performing $``$ again, $`_mP_{0m}^2\stackrel{}{}P_{00}_mP_{0m}`$ that is quasi-pure and $`E(\rho _3^{(2)})=\mathrm{log}3`$. For larger $`k`$, the cases above occur similarly.

Ex3 $`d=4`$

$`\rho _4^{(2)}`$ $`\stackrel{}{}(P_{00}+P_{20}){\displaystyle \underset{n}{}}P_{0n}+(P_{10}+P_{30}){\displaystyle \underset{n}{}}P_{2n},`$ (13a)
$`\rho _4^{(3)}`$ $`\stackrel{}{}P_{00}^2{\displaystyle \underset{n}{}}P_{0n}+P_{10}^2{\displaystyle \underset{n}{}}P_{3n}+P_{20}^2{\displaystyle \underset{n}{}}P_{2n}+P_{30}^2{\displaystyle \underset{n}{}}P_{1n},`$ (13b)
$`\rho _4^{(4)}`$ $`\stackrel{}{}{\displaystyle \underset{m}{}}P_{m0}^3{\displaystyle \underset{n}{}}P_{0n},`$ (13c)
$`\rho _4^{(5)}`$ $`\stackrel{}{}P_{00}^4{\displaystyle \underset{n}{}}P_{0n}+P_{10}^4{\displaystyle \underset{n}{}}P_{1n}+P_{20}^3{\displaystyle \underset{n}{}}P_{2n}+P_{30}^4{\displaystyle \underset{n}{}}P_{3n}`$ (13d)
$`\rho _4^{(3)}`$ and $`\rho _4^{(5)}`$ are exlicitly quasi-pure and $`E(\rho _4^{(3)})=4`$, $`E(\rho _4^{(5)})=8`$. It is required to explain for $`\rho _4^{(2)}`$ and $`\rho _4^{(4)}`$. As $`\varphi _{00}=(|00+|11)(|00+|11)`$ and $`\varphi _{20}=(|01+|10)(|00+|11)`$, so $`P_{00}+P_{20}=(P(|00+|11)+P(|01+|10))P(|00+|11)`$ that is quasi-pure. $`P_{10}+P_{30}`$ can be transformed to $`P_{00}+P_{20}`$ by local permutation the basis. So $`\rho _4^{(2)}`$ is quasi-pure and $`E(\rho _4^{(2)})=1`$. In the dual basis, $`_mP_{m0}^3`$ is represented as $`_nP_{0n}^3`$. Further applying the operation $``$, $`_nP_{0n}^3\stackrel{}{}P_{00}^2_nP_{0n}`$ that is quasi-pure. For larger $`k`$, the cases above occur similarly.

Ex4 d=5

$`\rho _5^{(2)}`$ $`\stackrel{}{}P_{00}{\displaystyle \underset{n}{}}P_{0n}+P_{10}{\displaystyle \underset{n}{}}P_{2n}+P_{20}{\displaystyle \underset{n}{}}P_{4n}+P_{30}{\displaystyle \underset{n}{}}P_{1n}+P_{40}{\displaystyle \underset{n}{}}P_{3n},`$ (14a)
$`\rho _5^{(3)}`$ $`\stackrel{}{}P_{00}^2{\displaystyle \underset{n}{}}P_{0n}+P_{10}^2{\displaystyle \underset{n}{}}P_{3n}+P_{20}^2{\displaystyle \underset{n}{}}P_{1n}+P_{30}^2{\displaystyle \underset{n}{}}P_{4n}+P_{40}^2{\displaystyle \underset{n}{}}P_{2n},`$ (14b)
$`\rho _5^{(4)}`$ $`\stackrel{}{}P_{00}^3{\displaystyle \underset{n}{}}P_{0n}+P_{10}^3{\displaystyle \underset{n}{}}P_{4n}+P_{20}^3{\displaystyle \underset{n}{}}P_{3n}+P_{30}^3{\displaystyle \underset{n}{}}P_{2n}+P_{40}^3{\displaystyle \underset{n}{}}P_{1n},`$ (14c)
$`\rho _5^{(5)}`$ $`\stackrel{}{}{\displaystyle \underset{m}{}}P_{m0}^4{\displaystyle \underset{n}{}}P_{0n},`$ (14d)
$`\rho _5^{(6)}`$ $`\stackrel{}{}{\displaystyle \underset{m}{}}P_{m0}^5\{{\displaystyle \underset{n}{}}P_{mn}\}.`$ (14e)
It is obvious that $`\rho _5^{(k)},k=2,3,4,6`$ are quasi-pure and $`E(\rho _5^{(k)})=(k1)\mathrm{log}5`$. In the dual basis, $`_mP_{m0}^4`$ is represented as $`_nP_{0n}^4`$. Further applying the operation $``$, $`_nP_{0n}^4\stackrel{}{}P_{00}^3_nP_{0n}`$. So $`\rho _5^{(5)}`$ is quasi-pure and $`E(\rho _5^{(5)})=3\mathrm{log}5`$. For larger $`k`$, the cases above occur similarly.

Ex5 d=6

$`\rho _6^{(2)}`$ $`\stackrel{}{}(P_{00}+P_{30}){\displaystyle \underset{n}{}}P_{0n}+(P_{10}+P_{40}){\displaystyle \underset{n}{}}P_{2n}+(P_{20}+P_{50}){\displaystyle \underset{n}{}}P_{4n},`$ (15a)
$`\rho _6^{(3)}`$ $`\stackrel{}{}(P_{00}^2+P_{20}^2+P_{40}^2){\displaystyle \underset{n}{}}P_{0n}+(P_{10}^2+P_{30}^2+P_{50}^2){\displaystyle \underset{n}{}}P_{3n},`$ (15b)
$`\rho _6^{(4)}`$ $`\stackrel{}{}(P_{00}^3+P_{30}^3){\displaystyle \underset{n}{}}P_{0n}+(P_{10}^3+P_{40}^3){\displaystyle \underset{n}{}}P_{4n}+(P_{30}^3+P_{50}^3){\displaystyle \underset{n}{}}P_{2n},`$ (15c)
$`\rho _6^{(5)}`$ $`\stackrel{}{}P_{00}^4{\displaystyle \underset{n}{}}P_{0n}+P_{10}^4{\displaystyle \underset{n}{}}P_{5n}+P_{20}^4{\displaystyle \underset{n}{}}P_{4n}+P_{30}^4{\displaystyle \underset{n}{}}P_{3n}+P_{40}^4{\displaystyle \underset{n}{}}P_{2n}+P_{50}^4{\displaystyle \underset{n}{}}P_{1n},`$ (15d)
$`\rho _6^{(6)}`$ $`\stackrel{}{}{\displaystyle \underset{m}{}}P_{m0}^5{\displaystyle \underset{n}{}}P_{0n},`$ (15e)
$`\rho _6^{(7)}`$ $`\stackrel{}{}{\displaystyle \underset{m}{}}P_{m0}^6\{{\displaystyle \underset{n}{}}P_{mn}\}.`$ (15f)
It is apparent that $`\rho _6^{(5)}`$ and $`\rho _6^{(7)}`$ are quasi-pure. Explanation is required for $`\rho _6^{(2)}`$, $`\rho _6^{(3)}`$, $`\rho _6^{(4)}`$, and $`\rho _6^{(6)}`$. Now we show that $`P_{00}+P_{30}`$ is quasi-pure.

$`\varphi _{00}`$ $`=`$ $`(|00+|11)(|00+|11+|22),`$
$`\varphi _{30}`$ $`=`$ $`(|01+|10)(|00+|11+|22),`$
$`P_{00}+P_{30}`$ $`=`$ $`(P(|00+|11)+P(|01+|10))P(|00+|11+|22),`$ (16)
which is explicitly quasi-pure. $`P_{10}+P_{40}`$ and $`P_{20}+P_{50}`$ can be proved similarly. So $`\rho _6^{(2)}`$ is quasi-pure and $`E(\rho _6^{(2)})=\mathrm{log}3`$. Now we show that $`P_{00}^2+P_{20}^2+P_{40}^2`$ is quasi-pure.

$`\varphi _{00}`$ $`=`$ $`(|00+|11+|22)(|00+|11),`$
$`\varphi _{20}`$ $`=`$ $`(|01+|12+|20)(|00+|11),`$
$`\varphi _{40}`$ $`=`$ $`(|02+|10+|21)(|00+|11),`$
$`P_{00}^2+P_{20}^2+P_{40}^2`$ $`=`$ $`(P_{00}^2(3)+P_{10}^2(3)+P_{20}^2(3))P_{00}^2(2).`$ (17)
From the $`d=3`$ case, we know that $`_mP_{m0}^2(3)`$ is quasi-pure. $`P_{10}+P_{30}+P_{50}`$ can be proved similarly. So $`\rho _6^{(3)}`$ is quasi-pure and $`E(\rho _6^{(2)})=\mathrm{log}3+2`$. Now we show $`P_{00}^3+P_{30}^3`$ is quasi-pure.

$$P_{00}^3+P_{30}^3=(P_{00}^3(2)+P_{10}^3(2))P_{00}^3(3).$$
(18)
From the $`d=2`$ case, we know that $`_mP_{m0}^3(2)`$ is quasi-pure. The other two terms can be proved simliarly. So $`\rho _6^{(4)}`$ is quasi-pure and $`E(\rho _6^{(4)})=3\mathrm{log}3+2`$. In the dual basis, $`_mP_{m0}^5`$ is represented as $`_nP_{0n}^5`$. Further applying the operation $``$

$$\underset{n}{}P_{0n}^5\stackrel{}{}P_{00}^4\underset{n}{}P_{0n}.$$
(19)
So $`\rho _6^{(6)}`$ is quasi-pure and $`E(\rho _6^{(6)})=4\mathrm{log}6`$. For larger $`k`$, the cases above occur similarly.

Notice that different conclusion is obtained Smolin  for $`k=2`$, where $`\rho _d^{(2)}`$ is demonstrated to be separable. And the relative entropy of entanglement is calculated based on this incorrect conclusion in a few papers. However, We prove that $`\rho _d^{(2)}`$ is distillable when $`d>2`$.

As well-known, irreversibility is generic in the entanglement distillation Vidal . In this paper, we show that the mixture of multiple copies of maximally entangled states is quasi-pure and prove that the entanglement distillation is a reversible process for this nontrivial class at first glance. Whether only in the quasi-pure states is distillation reversible deserves to be considered. Another important question is how to decide the quasi-pure mixed states.

D. Yang thanks S. J. Gu for helpful discussion. The work is supported by the NNSF of China, the NSF of Zhejiang Province (Grant No. 602018) and Guang-Biao Cao Foundation in Zhejiang University.