For submission to Monthly Notices
Self-consistent Gravitational Lens Reconstruction

Simon Dye & Andy Taylor
Institute for Astronomy, University of Edinburgh, Royal Observatory, Blackford Hill, Edinburgh EH9 3HJ, U.K.,

We present a new method for directly determining accurate, self-consistent cluster lens mass and shear maps in the strong lensing regime from the magnification bias of background galaxies. The method relies upon pixellisation of the surface mass density distribution which allows us to write down a simple, solvable set of equations. We also show how pixellisation can be applied to methods of mass determination from measurements of shear and present a simplified method of application. The method is demonstrated with cluster models and applied to magnification data from the lensing cluster Abell 1689.

Cosmology: theory – large–scale structure of the Universe, gravitational lensing; Galaxies: clusters

1 Introduction

The possibility of reconstructing cluster lens mass distributions from the magnification bias of background galaxies was first suggested by Broadhurst, Taylor & Peacock (1995) and first demonstrated by Taylor et al. (1998, T98 hereafter). They showed how a direct, local measure of the lens convergence, κ=Σ/Σc𝜅ΣsubscriptΣ𝑐\kappa=\Sigma/\Sigma_{c} , where ΣΣ\Sigma is the mass surface density and ΣcsubscriptΣ𝑐\Sigma_{c} is the critical surface density, could be obtained from knowledge of the lens magnification. In this way, one could measure absolute surface mass densities, thereby breaking the “sheet-mass” degeneracy found in methods based on distortions of background galaxies (Tyson, Valdes & Wenk 1990, Kaiser & Squires 1993, Seitz & Schneider 1995).

van Kampen (1998) and T98 have shown how one can extend magnification analysis into the strong lensing regime. By making reasonable assumptions about γ𝛾\gamma , the lens shear, they showed that one could place quite stringent bounds on κ𝜅\kappa . In addition, T98 found an exact solution for the profile of axisymmetric lenses, although not for more general 2D cases.

Inverse reconstruction methods based on maximum likelihood (Bartelmann et al. 1996) and maximum entropy (Seitz, Schneider & Bartelmann 1998, Bridle et al. 1998) have gone some way in providing a unification of both shear and magnification information. Until now however, no direct method using only magnification has existed.

In this letter, we show how to directly compute an accurate, self-consistent 2D distribution of κ𝜅\kappa and γ𝛾\gamma in the strong lensing regime from magnification. This direct approach has the advantage over indirect alternatives that uncertainties can easily be determined. The method is based on pixellisation of the κ𝜅\kappa distribution suggested by AbdelSalam, Saha & Williams (1998) who used it to estimate the mass of Abell 370 from multiple images. We generalise the method further and also derive a simplified solution to the problem of estimating mass from shear, based on the approach of Kaiser & Squires (1993).

2 Reconstruction of κ𝜅\kappa and γ𝛾\gamma

T98 showed how to estimate cluster surface mass using the magnification measured from the distortion in background galaxy number counts. Here our problem is to find an accurate method for reconstructing the surface mass density, given the magnification by an arbitrary lens. The inverse magnification factor at a given position in the lens plane is

A-1=|(1-κ)2-γ2|,superscript𝐴1superscript1𝜅2superscript𝛾2A^{-1}=|(1-\kappa)^{2}-\gamma^{2}|, (1)

where κ𝜅\kappa is the lens convergence and γ𝛾\gamma is the shear. The shear can be decomposed into two orthogonal polarisation states, γ1subscript𝛾1\gamma_{1} and γ2subscript𝛾2\gamma_{2} , which are related to the lens convergence by

γ1=12-2(12-22)κ,γ2=-212κ.formulae-sequencesubscript𝛾112superscript2subscriptsuperscript21subscriptsuperscript22𝜅subscript𝛾2superscript2subscript1subscript2𝜅\gamma_{1}=\frac{1}{2}\partial^{-2}(\partial^{2}_{1}-\partial^{2}_{2})\kappa,% \quad\gamma_{2}=\partial^{-2}\partial_{1}\partial_{2}\kappa. (2)

where i/θisubscript𝑖subscript𝜃𝑖\partial_{i}\equiv\partial/\partial\theta_{i} and -2superscript2\partial^{-2} is the 2D inverse Laplacian. The total shear is given by γ2=γ12+γ22superscript𝛾2subscriptsuperscript𝛾21subscriptsuperscript𝛾22\gamma^{2}=\gamma^{2}_{1}+\gamma^{2}_{2} . One might expect that equation ( 1) could be solved iteratively by first estimating κ𝜅\kappa , using this to calculate γ𝛾\gamma and then updating the estimate of κ𝜅\kappa using equation ( 1) again. This proves to be highly unstable in the strong lensing regime however, rapidly diverging after only a few iterations (Seitz & Schneider 1995).

To find a stable solution to equation ( 1), we first pixellise the image. Following AbdelSalam et al. (1998), we can now write

γin=Dimnκm,i=1,2formulae-sequencesuperscriptsubscript𝛾𝑖𝑛superscriptsubscript𝐷𝑖𝑚𝑛subscript𝜅𝑚𝑖12\gamma_{i}^{n}=D_{i}^{mn}\kappa_{m},\quad i=1,2 (3)

with summation implied over index m𝑚m and where κmsubscript𝜅𝑚\kappa_{m} and γinsuperscriptsubscript𝛾𝑖𝑛\gamma_{i}^{n} are the pixellised convergence and shear distributions respectively. The transformation matrices, Dimnsuperscriptsubscript𝐷𝑖𝑚𝑛D_{i}^{mn} , are

D1mnsuperscriptsubscript𝐷1𝑚𝑛\displaystyle D_{1}^{mn} =\displaystyle= 12(12-22)md2θln|𝜽n-𝜽|12superscriptsubscript12subscriptsuperscript22subscript𝑚superscriptd2superscript𝜃subscript𝜽𝑛superscript𝜽\displaystyle\frac{1}{2}(\partial_{1}^{2}-\partial^{2}_{2})\int_{m}{\rm d}^{2}% \theta^{{}^{\prime}}\ln|\mbox{\boldmath$\theta$}_{n}-\mbox{\boldmath$\theta$}^% {{}^{\prime}}| (4)
=\displaystyle= 1πtan-1[x12-x22(x12+x22)2-1/4],1𝜋superscript1superscriptsubscript𝑥12superscriptsubscript𝑥22superscriptsuperscriptsubscript𝑥12superscriptsubscript𝑥22214\displaystyle\frac{1}{\pi}\tan^{-1}\left[\frac{x_{1}^{2}-x_{2}^{2}}{(x_{1}^{2}% +x_{2}^{2})^{2}-1/4}\right],


D2mnsuperscriptsubscript𝐷2𝑚𝑛\displaystyle D_{2}^{mn} =\displaystyle= 12md2θln|𝜽n-𝜽|subscript1subscript2subscript𝑚superscriptd2superscript𝜃subscript𝜽𝑛superscript𝜽\displaystyle\partial_{1}\partial_{2}\int_{m}{\rm d}^{2}\theta^{{}^{\prime}}% \ln|\mbox{\boldmath$\theta$}_{n}-\mbox{\boldmath$\theta$}^{{}^{\prime}}| (5)
=\displaystyle= 12πln[(x12+x22)2-2x1x2+1/4(1/2+x12+x22)2-(x1-x2)2],12𝜋superscriptsuperscriptsubscript𝑥12superscriptsubscript𝑥2222subscript𝑥1subscript𝑥214superscript12superscriptsubscript𝑥12superscriptsubscript𝑥222superscriptsubscript𝑥1subscript𝑥22\displaystyle\frac{1}{2\pi}\ln\left[\frac{(x_{1}^{2}+x_{2}^{2})^{2}-2x_{1}x_{2% }+1/4}{(1/2+x_{1}^{2}+x_{2}^{2})^{2}-(x_{1}-x_{2})^{2}}\right],

with the integration acting over the mthsuperscript𝑚𝑡m^{th} pixel. 𝒙=𝜽n-𝜽m𝒙subscript𝜽𝑛subscript𝜽𝑚\mbox{\boldmath$x$}=\mbox{\boldmath$\theta$}_{n}-\mbox{\boldmath$\theta$}_{m} is the difference between pixels m𝑚m and n𝑛n which are assumed to be square in calculating these analytic expressions. Equation ( 1) can now be written as the vector equation,

𝟏-2𝜿+𝜿𝑮𝜿t-𝒫𝑨-1=012𝜿superscript𝜿𝑮𝜿t𝒫superscript𝑨10\mbox{\boldmath$1$}-2\mbox{\boldmath$\kappa$}+\mbox{\boldmath$\kappa$}\mbox{% \boldmath$G$}\mbox{\boldmath$\kappa$}^{\rm t}-{\cal P}\mbox{\boldmath$A$}^{-1}=0 (6)

where 𝑨-1superscript𝑨1\mbox{\boldmath$A$}^{-1} is the N𝑁N -dimensional vector of pixellised inverse magnification values, 𝜿tsuperscript𝜿t\mbox{\boldmath$\kappa$}^{\rm t} is the transpose of the vector 𝜿𝜿\kappa of pixellised convergence values and 𝟏11 is the vector (1,1,1,)111(1,1,1,\cdots) . The matrix 𝑮𝑮G is the N×N×N𝑁𝑁𝑁N\times N\times N matrix

Gpqn=δpnδqn-D1pnD1qn-D2pnD2qnsubscript𝐺𝑝𝑞𝑛subscript𝛿𝑝𝑛subscript𝛿𝑞𝑛superscriptsubscript𝐷1𝑝𝑛superscriptsubscript𝐷1𝑞𝑛superscriptsubscript𝐷2𝑝𝑛superscriptsubscript𝐷2𝑞𝑛G_{pqn}=\delta_{pn}\delta_{qn}-D_{1}^{pn}D_{1}^{qn}-D_{2}^{pn}D_{2}^{qn} (7)

where δijsubscript𝛿𝑖𝑗\delta_{ij} is the Krönecker delta, and summation is only over indices p𝑝p and q𝑞q . The parity of the measured inverse amplification A-1(𝜽)superscript𝐴1𝜽A^{-1}(\mbox{\boldmath$\theta$}) is handled by 𝒫𝒫{\cal P} which flips from being +11+1 outside regions bounded by critical lines to -11-1 within such regions.

The amplification equation in the form of equation ( 6) is the first main result of this letter. We can now solve for κ𝜅\kappa given a measured inverse amplification. Having solved for κ𝜅\kappa , the corresponding shear distribution can then be calculated from equation ( 3).

3 Application to Cluster Models

We apply the method to two types of idealised cluster models. Starting with a predetermined cluster mass density distribution, the corresponding shear distribution is derived using Fourier methods (see for example, Bartelmann & Weiss 1994). From these, the resulting magnification is calculated from equation ( 1) and then windowed to remove boundary effects. Using equation ( 6), we solve for κ𝜅\kappa . γ𝛾\gamma is then solved using equation( 3). A grid of 32 by 32 pixels is used in both models.

3.1 Truncated Isothermal Sphere Model

We first test the method with a simple truncated isothermal lens model. The pixellated mass distribution is laid down using κ(r+r0)-1proportional-to𝜅superscript𝑟subscript𝑟01\kappa\propto(r+r_{0})^{-1} , where r𝑟r is the radial distance from the centre of the sphere and r0subscript𝑟0r_{0} is a constant.

Figure 1shows the κ𝜅\kappa and γ𝛾\gamma distribution from which the magnification distribution was calculated, the solved κ𝜅\kappa and γ𝛾\gamma distribution and the difference between them. The plotted distributions are smoothed from the underlying grid and the white dashes highlight the lens’ critical line. The residuals are shown as percentage deviations from the true distribution. These are less than one percent for κ𝜅\kappa over most of the grid which is negligible in comparison to the errors typically found in practice from background clustering, shot noise (T98) and the uncertainties resulting from use of local κ𝜅\kappa estimators (see van Kampen 1998). The recovered shear distribution is more affected although still fares better than γ𝛾\gamma calculated from uncorrected Fourier techniques. The main contribution to these residuals is from boundary effects arising from trying to recover a nonlocal shear in a finite area. Since much work has been carried out into the removal of such effects (see Squires & Kaiser 1996 and Seitz & Schneider 1995 for example) which have little impact on the recovered κ𝜅\kappa , we shall address the problem elsewhere.

Figure 1: Truncated Isothermal Sphere Model: The initial κ𝜅\kappa and γ𝛾\gamma used to form the magnification distribution from which the solved κ𝜅\kappa and γ𝛾\gamma are derived. Underlying grid dimensions are 32 by 32. White dashes show the position of the critical line. Contours are linearly spaced and set at the same levels in both κ𝜅\kappa plots and in both γ𝛾\gamma plots. Residuals are expressed as percentages of (κinit-κsolved)/κinitsubscript𝜅initsubscript𝜅solvedsubscript𝜅init(\kappa_{\rm init}-\kappa_{\rm solved})/\kappa_{\rm init}

3.2 Dumb-bell Mass Model

The method was also tested with a more general dumb-bell model. Magnification was determined in the same fashion as for the isothermal model, setting a negative parity inside the critical lines, shown by the white dashes in figure 2. Once again the residuals between the initial and solved κ𝜅\kappa are typically less than one percent, while those for γ𝛾\gamma are typically 10%percent1010\% and again come mainly from boundary effects.

Figure 2: Dumb-bell Model: Critical lines are shown as white dashes. Underlying grid dimensions are 32 by 32. Linearly spaced contours are set at the same levels for κ𝜅\kappa and at the same levels for γ𝛾\gamma .

4 Practical Considerations

We solve equation ( 6) with the hybrid Powell method (NAG routine C05PCF). The number of equations needed to solve for κ𝜅\kappa is equal to the total number of grid pixels which can prove computationally intensive for especially fine grids. We find that this is not a problem for grid resolutions used to measure magnification bias in practice however. The 32 by 32 grid of pixels used for the models in section 3was solved in approximately one minute on an average workstation. The residuals exhibit no noticeable dependence on grid size.

The Powell algorithm is an iterative process and therefore requires an initial estimate of the solution to start from. The choice of the initial estimate turns out to be irrelevant. We have tried a wide range of initial distributions and even starting from a uniform distribution arrived at the same final solution.

We have found that correct choice of pixel parity (especially for low grid resolutions) is essential in order to achieve a sensible result. Inappropriate assignment of parities to pixels manifest themselves, as one would expect, by κ𝜅\kappa being overestimated when a pixel is wrongly assumed to lie inside a critical line and underestimated in the reverse situation. This gives a means of checking whether critical line positions have been properly defined by looking for large discontinuities in the κ𝜅\kappa distribution. Models with dual critical lines requiring dual parity flips have also been tested and we find that κ𝜅\kappa can be recovered just as well.

Finally, to ensure that the method does not break down with noisy data, we introduced a random noise term to the amplification. Errors in κ𝜅\kappa resulting from noise in the inverse amplification propagate as one would expect from equation ( 6). For an isothermal lens we recovered the expected result, δκ=δA/2A2𝛿𝜅𝛿𝐴2superscript𝐴2\delta\kappa=\delta A/2A^{2} , indicating that pixellisation does not lead to spurious noise properties.

5 Application to Abell 1689

We apply the method to the magnification data presented in T98 for the lensing cluster A1689. A 12 by 12 grid is used as the best compromise between shot noise in galaxy counts per bin and the resolution of the derived κ𝜅\kappa map. Identification of the critical line was achieved by locating giant arc positions in the observed image.

Figure 3shows the solved mass density and shear distribution. Comparison with the mass density map illustrated in T98 (figure 6) which was produced with the sheet κ𝜅\kappa estimator shows very similar structure. We find that the value of κ𝜅\kappa at the peak calculated here is approximately 10% lower than the peak value in T98 since the sheet estimator over-estimates κ𝜅\kappa inside critical line regions. This has little effect on the total integrated mass of A1689 found in T98. The γ𝛾\gamma distribution is shown for completeness although undoubtedly suffers from boundary effects typically found in the models.

Figure 3: A1689 solved convergence and shear distributions. Darker areas represent a higher distribution density. White dashes show the observed critical line. The plots are smoothed from a 12 by 12 grid with north up and east to the left.

6 Shear Analysis

Having shown that pixellisation allows us to accurately reconstruct surface mass densities from magnification data, we now apply it to shear analysis. Shear analysis exploits the idea that a given distribution of images of galaxies lying behind a lensing cluster will, in the statistical mean, have regions of lens-induced correlations in image orientation and ellipticity. Measuring the quadrupole moments of individual galaxy images enables the construction of a map of the ellipticity parameters, eijsubscript𝑒𝑖𝑗e_{ij} (Valdes, Tyson & Jarvis 1983). The ellipticity parameters relate to the surface mass density and shear via (Kaiser 1995, K95 hereafter):

eij=γij1-κ,γij=(γ1γ2γ2-γ1)formulae-sequencesubscript𝑒𝑖𝑗subscript𝛾𝑖𝑗1𝜅subscript𝛾𝑖𝑗subscript𝛾1subscript𝛾2subscript𝛾2subscript𝛾1e_{ij}=\frac{\gamma_{ij}}{1-\kappa},\quad\gamma_{ij}=\left(\begin{array}[]{cc}% \gamma_{1}&\gamma_{2}\\ \gamma_{2}&-\gamma_{1}\\ \end{array}\right) (8)

One way of solving this for κ𝜅\kappa in the weak lensing regime is to follow the approach of Kaiser & Squires (1993). Generalisations of this to the strong regime have been made by K95. One would have hoped that an alternative to such approaches would be to pixellise equation ( 8) and use equation ( 3) to solve it by matrix inversion. However, the resulting matrix equation is ill-conditioned, since the matrix D1mnsubscriptsuperscript𝐷𝑚𝑛1D^{mn}_{1} is singular and D2mnsubscriptsuperscript𝐷𝑚𝑛2D^{mn}_{2} is itself ill-conditioned. Instead, we show a new, simplified expression for the solution to Kaisers’ ellipticity equation and then pixellise it.

Starting with the equation (K95),

iκ=jγijsubscript𝑖𝜅subscript𝑗subscript𝛾𝑖𝑗\partial_{i}\kappa=\partial_{j}\gamma_{ij} (9)

and using equation ( 8), one can show that

iln(1-κ)=-jln(δij+eij).subscript𝑖1𝜅subscript𝑗subscript𝛿𝑖𝑗subscript𝑒𝑖𝑗\partial_{i}\ln(1-\kappa)=-\partial_{j}\ln(\delta_{ij}+e_{ij}). (10)

The term on the right hand side is obtained from the definition,

ln(I+𝑩)=𝑩+12𝑩2+13𝑩3+𝐼𝑩𝑩12superscript𝑩213superscript𝑩3\ln(I+\mbox{\boldmath$B$})=\mbox{\boldmath$B$}+\frac{1}{2}\mbox{\boldmath$B$}^% {2}+\frac{1}{3}\mbox{\boldmath$B$}^{3}+\cdots (11)

where I𝐼I is the identity matrix and 𝑩𝑩B is an arbitrary square matrix. Using this expansion and collecting even and odd terms we find,

ln(δij+eij)=-12ln(1-e2)δij+12ln(1+e1-e)eije,subscript𝛿𝑖𝑗subscript𝑒𝑖𝑗121superscript𝑒2subscript𝛿𝑖𝑗121𝑒1𝑒subscript𝑒𝑖𝑗𝑒\ln(\delta_{ij}+e_{ij})=-\frac{1}{2}\ln(1-e^{2})\delta_{ij}+\frac{1}{2}\ln% \left(\frac{1+e}{1-e}\right)\frac{e_{ij}}{e}, (12)

where e2=e12+e22superscript𝑒2superscriptsubscript𝑒12superscriptsubscript𝑒22e^{2}=e_{1}^{2}+e_{2}^{2} and ei=γi/(1-κ)subscript𝑒𝑖subscript𝛾𝑖1𝜅e_{i}=\gamma_{i}/(1-\kappa) . This result requires that e<1𝑒1e<1 . Inserting equation ( 8) into the magnification equation ( 1) we find

A-1=|(1-κ)2(1-e2)|.superscript𝐴1superscript1𝜅21superscript𝑒2A^{-1}=|(1-\kappa)^{2}(1-e^{2})|. (13)

Hence the parity changes when e>1𝑒1e>1 . Since eijsubscript𝑒𝑖𝑗e_{ij} and eij-1subscriptsuperscript𝑒1𝑖𝑗e^{-1}_{ij} are observationally indistinguishable and flip from one to another whenever there is a parity change, we can satisfy the criterion e<1𝑒1e<1 just by noting the critical line positions and inverting the ellipticity matrix when one is crossed.

Finally, inserting equation ( 12) into equation ( 10), and solving for κ𝜅\kappa we find the pixellised solution is

κn=1-(1-en2)1/2exp[-12(D1nms1m+D2nms2m)]subscript𝜅𝑛1superscript1superscriptsubscript𝑒𝑛21212subscriptsuperscript𝐷𝑛𝑚1subscriptsuperscript𝑠𝑚1subscriptsuperscript𝐷𝑛𝑚2superscriptsubscript𝑠2𝑚\kappa_{n}=1-(1-e_{n}^{2})^{1/2}\exp\left[-\frac{1}{2}(D^{nm}_{1}s^{m}_{1}+D^{% nm}_{2}s_{2}^{m})\right] (14)


si=eieln(1+e1-e),i=1,2formulae-sequencesubscript𝑠𝑖subscript𝑒𝑖𝑒1𝑒1𝑒𝑖12s_{i}=\frac{e_{i}}{e}\ln\left(\frac{1+e}{1-e}\right),\quad i=1,2 (15)

Equation ( 14) is the second main result of this letter. We can directly calculate κ𝜅\kappa given a measured ellipticity field. Figure 4shows the results of reconstructing κ𝜅\kappa using equation ( 14) for the dumb-bell model. The ellipticity parameters are calculated from equation ( 8) using the κ𝜅\kappa and γ𝛾\gamma distribution. We normalise the reconstructed κ𝜅\kappa to both peaks in the initial κ𝜅\kappa distribution. The residuals, again being dominated by boundary effects, show that reconstruction is possible to within approximately 10 percent across the field of view.

Figure 4: Reconstruction of κ𝜅\kappa from the ellipticity parameters. Contours are at the same levels in both κ𝜅\kappa plots. The distortion field is illustrated by plotting the apparent shape of an intrinsically circular background object.

7 Summary

We have outlined a method for directly calculating accurate, self-consistent surface mass density and shear distributions from the lens amplification and critical line positions. The method has been demonstrated with the isothermal sphere and dumb-bell cluster models. We find it reconstructs the surface density to within a percent over most of the field of view. The reconstruction of the shear pattern only has fractional accuracy of a few tenths due to boundary effects. We have applied the method to magnification data from Abell 1689, and reconstructed its surface mass and shear distribution.

We have also found a simplified solution to the problem of estimating surface mass density from galaxy ellipticities. This approach puts the calculation of surface mass from shear and magnification on an equal footing, and we shall investigate the combined analysis elsewhere.


SD is supported by a PPARC studentship. ANT is a PPARC research associate. Thanks to our collaborators Tom Broadhurst, Txitxo Benitez, Eelco van Kampen for allowing us to use data from Abell 1689. We also thank Hanadi AbdelSalam for bringing the benefits of pixellisation to our attention and for useful discussion.


AbdelSalam H.M., Saha P., Williams L.L.R., 1998, MNRAS, 294, 734

Bartelmann M., Narayan R., Seitz S., Schneider P., 1996, ApJ, 464, 115

Bartelmann M., Weiss A., 1994, A&A, 287, 1

Bridle S.L., Hobson M.P., Lasenby A.N., Saunders R., 1998, astro-ph/9802159

Broadhurst T.J., Taylor A.N., Peacock J.A., 1995, ApJ, 438, 49

Kaiser N., (K95), 1995, ApJ, 439, L1

Kaiser N., Squires G., 1993, ApJ, 404, 441

Seitz C., Schneider P., 1995, A&A, 297, 287

Seitz C., Schneider P., Bartelmann M., 1998, astro-ph/9803038

Squires G., Kaiser N., 1996, ApJ, 473, 85

Taylor A.N., Dye S., Broadhurst T.J., Benitez N., van Kampen E., (T98), 1998, ApJ, 501, 539

Tyson J.A., Valdes F., Wenk R.A., 1990, ApJ, 349, L1

Valdes F., Tyson J., Jarvis J., 1983, ApJ, 271, 431

van Kampen, 1998, submitted to MNRAS