S T R O N G G R A V I T A T I O N A L L E N S I N G : B L U E P R I N T S F O R G A L A X Y - C L U S T E R C O R E R E C O N S T R U C T I O N by PETER ROBERT NEWBURY B.Sc. (Mathematics) University, of Manitoba j 1990 M.Sc. (Applied Mathematics) University of British Columbia, 1993 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF T H E REQUIREMENTS FOR T H E DEGREE OF - DOCTOR OF PHILOSOPHY in T H E FACULTY OF GRADUATE STUDIES Department of Mathematics Institute of Applied Mathematics We accept this thesis as conforming to the required standard T H E UNIVERSITY OF BRITISH COLUMBIA March 1998 © Peter Robert Newbury, 1998 In presenting this thesis in partial fulfillment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for refer-ence and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of Mathematics The University of British Columbia Vancouver, Canada Date rAoA/JU 13 ) 13 3 & Abstract When rays of light pass by a massive object they are very slightly deflected towards the centre-of-mass of the object. If two or more diverging beams of light re-converge onto an serendipitous observer, this observer may see multiple, magnified images of the source of light. This process is known as gravitational lensing, and has been observed in several dozen spectacular cases. Based on the appearance of the lensed arcs of light, we attempt to "invert" the lens to find the distribution of mass that will produce just such a configuration of lensed objects. In this thesis, we propose a two-stage inversion scheme. First, the distribution of mass on the deflector plane and the geometry of the source-deflector-observer optical system are established. This is done by numerically simulating the lensing of light past a parametric mass model, and interactively adjusting the handful of model parameters to match the positions of the simulated and observed lensed arcs. At the same time, this determines the magnification of the background source induced by the lensing process. The predicted magnification is then removed from the data to reveal the intrinsic, though still distorted, background distribution of light. After tracing each lensed ray back to the source plane, the data are recombined to produce a surface brightness distribution of the source. This two-stage inversion scheme produces a parametric model of the deflector and a pixelised rendering of the background source which together mimic the observed gravitationally lensed features. We test the viability of scheme itself on a well-studied collection of lensed objects in the galaxy-cluster MS 2137. Confident in the algorithm, we apply it second time to predict the distribution of mass in the galaxy-cluster MS 1455 responsible for an observed triplet of lensed arcs. Our predictions about the lens in MS 1455 make it particularly interesting, for a single background source is responsible for both tangential arcs and a radial arc. Table of Contents Abstract ii Table of Contents iii List of Tables v List of Figures vi Acknowledgment xi Foreword xii Chapter 1. Introduction 1 Chapter 2. Review of Gravitational Lensing 7 2.1 Before 1937: Hypothesis 7 2.2 The Period 1937-1979: Analysis 8 2.3 The Period 1979-1991: Discovery 9 2.4 After 1991: Lens Inversion 9 2.4.1 Microlensing 10 2.4.2 Weak Lensing 11 2.4.3 Inversion of Strong Lenses 14 Chapter 3. Gravitational Lens Optics 22 3.1 Deflection of Light 22 3.1.1 Perturbations of Flat Spacetime 23 3.1.2 The Deflection Angle 24 3.1.3 Elliptical Mass Distributions 27 3.2 The Lens Equation 30 3.3 Properties of the Lens Equation 34 3.4 A Toolbox of Lenses 36 3.4.1 Mass Density 39 3.4.2 Cumulative Mass 40 3.4.3 Surface Mass Density 41 3.4.4 Deflection Angle 43 3.4.5 Amplification Factor 47 3.4.6 Line-of-Sight Velocity Dispersion 50 Chapter 4. Numerical Models of Gravitational Lensing 58 4.1 Modeling Scheme 58 4.2 Model Parameters 63 4.2.1 Lens Geometry 63 4.2.2 Deflector Mass Distribution 63 iii iv 4.2.3 Source Light Distribution 65 4.3 Degrees of Freedom in a Parametric Lens Model 68 4.4 Other Features of the Modeling Software 69 Chapter 5. A New Lens Inversion Scheme 74 5.1 Why do we need another lens inversion scheme? 74 5.2 A Two-Stage Inversion Algorithm 76 5.3 Specifying Lensed Image Positions with Polar Moments 77 5.3.1 Applicability of Polar Moments 84 5.4 Optimal Source Reconstruction 85 Chapter 6. Gravitational Lensing in MS 2137 and MS 1455 98 6.1 The Collection of Arcs in MS 2137 98 6.1.1 Existing Models of MS 2137 99 6.2 Inversion of MS 2137 using the New Lens Inversion Scheme 103 6.2.1 Establishing the Lens Geometry and Mass Distribution 103 6.2.2 Other Arcs in MS 2137 109 6.2.3 Reconstructing the Sources 113 6.2.4 Reconstructions of MS 2137 117 6.3 Lensed Arcs in MS 1455 120 6.4 Inversion of the Lens MS 1455 123 6.4.1 Lens Geometry and Mass Distribution 123 6.4.2 Models of MS 1455 123 6.4.3 Source Reconstruction 129 6.4.4 Reconstruction of MS 1455 132 6.4.5 The Extended Appearance of Arc A3 135 6.4.6 Origin of the Spurious Data 143 6.5 The Inversions of MS 2137 and MS 1455 146 Chapter 7. Conclusions 147 7.1 The New Inversion Scheme 147 7.2 Strengths and Weaknesses of the Scheme 151 7.3 Outlook 153 7.3.1 Lensing as an Inverse Problem 153 7.3.2 Astrophysical Applications 155 7.4 Conclusion 159 References 160 List of Tables 3.1 Parameter values for the lens masses. The observed value of oios[cD](2'.'0) is explained in Section 3.4.6 39 4.1 Parameters determining the geometry and resolution of the simulation 64 4.2 Each mass distribution in the lens is described by a small set of model parameters. . 64 4.3 Uniform circular disk source distribution of light. In the example shown here, ys = (0,0), S = 5, r = 0.5 65 4.4 A pseudo-isothermal profile mimics a background elliptical galaxy. This example is de-scribed by ys = (0,0), S = 5, r = 0.5. The core radius is indicated by the circle of radius 0'/5 66 4.5 The shape and orientation of the uniform elliptical disk can be altered to mimic the shape of any background elliptical galaxy. The model shown here is desribed by ys = (0,0), r = 0.5, S = 5, e = 0.7,

0 describes a family of ellipses confocal to the central A = 0 ellipse. The confocal ellipses become more and more circular as the mass lying on the A = 0 shell appears more and more point-like. The value of A defines the distance between a point £ and the mass shell at A = 0 29 3.3 A ray emitted at rj on the source plane pierces the deflector plane at £. The ray is deflected by a (£ ) onto the observer 31 3.4 The angular size distance to objects at redshift z 32 3.5 Three-dimensional mass density p(r) 40 3.6 The 3-dimensional cumulative mass M(r) of the lens models 42 3.7 Two-dimensional surface mass density S(r) of the lens models found by projecting the mass onto the lens plane under the thin lens approximation 44 3.8 The deflection angle a(r) developed about model SIS, PID, and NFW lenses, with and without a central cD galaxy. Curves marked by o are the corresponding dark-matter halos with the cD distribution added at the centre. The dark solid line follows a = r and locates the critical lines developed by the lens (Section 3.4.5) 46 3.9 The Mean Value Theorem guarantees the existence of a point r' between 0 and r where %(r') = 1 48 3.10 The magnification factor \p\(r) around each lens in this collection. Curves marked by o are the corresponding dark-matter halos with the cD distribution added at the centre. The SIS has only one tangential critical line at roughly 17'.'0. The PID and NFW lenses have both outer (tangential) and inner (radial) critical lines 49 3.11a Line-of-sight velocity dispersion profile of a SIS+cD mass distribution. The mass pa-rameter of the SIS adjusted to give o-los[SIS](50'!0) = 1000 kms - 1 , marked by A. The interpolated temperature of the cD component at 2''0 is marked * 54 3.11b Line-of-sight velocity dispersion profile of a PID+cD mass distribution. The mass pa-rameter of the PID adjusted to give aios[PID](50'!0) = 1000 kms - 1 , marked by A. The interpolated temperature of the cD component at 2'/0 is marked * 55 3.11c Line-of-sight velocity dispersion profile of a NFW+cD mass distribution. The mass parameter of the NFW adjusted to give o-los[NFW]{50'.'0) = 1000 kms - 1 , marked by A. The interpolated temperature of the cD component at 2'.'0 is marked * 56 Vll 4.1 The visible array is a linear combination of the arrays generated by each deflector com-ponent. The registration of the field-of-view on each component's layer determines that component's position in the simulation 60 4.2 A simple example with half-width W = 5 and half-height H = 4. The offset (i0,jo)R = (8,6) registers the visible field on the underlying global grid. This registration relates the visible coordinates (3,3)y to the global coordinates (11, 9)G of the Object via Equation (4.1). The global grid is just large enough that the visible "window" can be repositioned with the Object appearing at any location (1, l)y thru (11,9)y 62 4.3 The heavy dashed lines mark the critical lines of the lens where the amplification di-verges. Thin solid lines trace 10x,5x amplication contours (even parity). Thin dashed lines trace contours of —lOx, —5x amplication (odd parity) 70 4.4 Circles, rays, and annular sectors, built from polar moments (Section 5.3) characterise the position and shape of lensed arcs 72 4.5 Lensed arcs in Figure 4.3 convolved with a Moffat point-spread function (4.2) with a 2-pixel FWHM 72 4.6 Gaussian noise mimicking the noise in the observations is added to the simulation. 73 5.1 The two-stage lens inversion algorithm we propose in this thesis 77 5.2 Polar r- and ^-coordinates of the three gravitationally lensed images in the simulation shown in Figure 1.2 are plotted in a Cartesian coordinate system 79 5.3 The sign of the off-diagonal polar moment QTq gives a qualitative description of how the object is rotated inside the annular sector surrounding the polar centroid (6, f). . . . 83 5.4 Neither setting the source plane pixels to the same physical size nor the same angular size as the deflector plane pixels satisfactorily covers the source plane 88 5.5 The "PN" test pattern is lensed through a single PID deflector. Figure 5.6 shows the source distributions reconstructed from the data in this image 90 5.6 Reconstruction of the "PN" test pattern on the source plane at various resolutions from the lensed images in Figure 5.5 91 5.7 The pixelated sources in Figure 5.6 are passed back through the lens. Because larger source pixels spread the information out over larger regions of the source plane, the lensed images contain extraneous lensed light not present in the original simulation shown in Figure 5.5 94 5.8 The software inadvertently reconstructs a source for virtually the entire image plane after the "PN" lensed images are convolved with a point-spread function and noise is added. Dots mark the points where re-traced lensed rays pierce the source plane. The over- and under-sampling of the lens equation is obvious. The astroid and ovoid caustics marked by dotted lines in Figure 5.5 are visible. The limited field-of-view on the deflector planes makes the outer regions of the source plane invisible to the observer 96 6.1 WFPC2 image of MS 2137. The grayscale is linear. The giant arc AO and arclets A2 and A4 are traced to one source. A second source produces the radial arc AR and the arclet A6. Our models suggest that B l , BR, C l , and CR are also lensed arcs. The central cD galaxy G l and the galaxy G7 which obscures part of arc AO are removed in Figure 6.2. North is at an angle of —38.5° relative to the positive rc-axis, based on the HST orientation 100 vm 6.2 Processed WFPC2 image of MS 2137. The grayscale is linear. The central cD galaxy G l and a smaller galaxy G7 obscuring the upper end of arc AO have been removed, revealing more of the structure of the arcs. The irregular under-densities left behind are numerical relics from the removal of the galaxies 101 6.3 Simulation of the arcs in MS 2137 showing the polar moments of each arc. While the surface brightness of the arcs plays no part in the calculation of the moments, the brightness of the two sources is adjusted to mimic the flux in the observed arcs. . . 108 6.4 Two additional pairs of arcs, Bl-BR and Cl-CR, and annular sectors built from their polar moments listed in Table 6.5 110 6.5 Simulation of the Bl-BR pair of arcs. The background elliptical source lies at redshift 0.900 I l l 6.6 Simulation of the Cl-CR pair of arcs. The background elliptical source lies at redshift 1.725 112 6.7 Reconstructed sources SI (top right) and S2 (top left) with grayscales adjusted to bring up the structure of each. The relative brightness of the two sources becomes apparent when they appear together on the source plane. The dashed ellipses mark the positions of the uniform elliptical sources used to find the geometry of the lens and the mass distributions. The dotted lines trace the caustics of the lens 115 6.8 The two-source reconstruction is traced back through the gravitational lens model, and gaussian noise is added to match the WFPC2 data 118 6.9 Cross-section of the surface brightness of source SI parallel to the caustic from point A to point B in the left panel. The two peaks are produced by the double ring structure in the giant arc, AO 119 6.10 An equivalent 4-hour exposure of the X-ray cluster MS 1455 with square root scaling.The candidate tangential arc identified by LeFevre et al. (1994) is labeled A2. The arrow indicates a radial structure in the light of the central cD galaxy. 121 6.11 MS 1455 with the central cD galaxy removed, revealing a candidate radial arc A l . Subsequent modeling of the lens suggests the object A3 is also a gravitationally lensed image. Circles and annular sectors show the polar moments of the arcs 122 6.12 Simulated gravitational lensing in MS 1455 with the PID+cD model. Three lensed features are produced by this configuration. The polar moments of the arcs, represented by the annular sectors surrounding each arc, are tabulated in Table 6.9 127 6.13 Like Figure 6.12 but with a NFW dark-matter halo 128 6.14 The source reconstructed from the data in MS 1455 observed through the PID+cD (left) and NFW+cD (right) lens models. Dotted lines trace the caustics of the lens; the dashed ellipse marks the position of the uniform source used to simulate the lensed features. Both reconstructions suggest the background source is a spiral galaxy. The spurious, isolated pixels in the lower half portions of the plot are discussed in the text 130 6.15 The cleaned reconstructed sources in the PID+cD (left) and NFW+cD (right) lens models. The spurious isolated pixels have been removed, while the genuine faint limb on the right of the source remains 131 6.16 The source reconstructed behind the PID+cD model is passed back through the lens. The simulation is degraded to match the seeing at the CFHT. While the radial arc A l and the arclet A2 are quite well reproduced, tangential arc A3 is much more extended than observed 133 < ix 6.17 The source reconstructed behind the NFW+cD model is passed back through the lens and degraded. Again, the radial arc A l and the arclet A2 are quite well reproduced, while the tangential arc A3 is highly extended 134 6.18a (left) Source reconstructed from the radial arc A l through the PID+cD lens model. Some spurious pixels have been removed, leaving only 245 of the 257 data, (right) The appearance of the source through the lens 136 6.18b (left) Source reconstructed from 287 pixels in giant arc A2 through the PID+cD lens model, (right) Its appearance through the lens 136 6.18c (left) Source reconstructed from 172 pixels in the third arc A3 through the PID+cD lens model, (right) Its appearance through the lens 137 6.19a (left) Source reconstructed from the radial arc A l through the NFW+cD lens model. The source has been cleaned, leaving only 234 of the 257 data selected from arc A l in the observations, (right) Its appearance through the lens 138 6.19b (left) Source reconstructed from 287 pixels in giant arc A2 through the NFW+cD lens model, (right) Its appearance through the lens 138 6.19c (left) Source reconstructed from 172 pixels in the third arc A3 through the NFW+cD lens model. Note the decrease in flux, (right) Its appearance through the lens. . . . 139 6.20 (left) Source reconstructed by combining the A l and A2 data sets. Cleaning removes 53 of 544 data. The pixel size and grayscale is the same as Figure 6.18a. (right) Its appearance through the lens 140 6.21 Simulation shown in Figure 6.20 is convolved with a point-spread function and noise is added. The position and shape of arc A3 rs very similar to the observations 141 6.22 Convolution of the original PID+cD simulation shown in Figure 6.12, in which the source is a featureless, elliptical disk. Differences in intensity of the arcs are due to magnification. Convolution dramatically changes the shape of arc A3, while the noise introduces some sub-structure in the arcs 142 6.23 Each lensed pixel forming the radial arc in the observations is plotted in the PID+cD (left) and NFW+cD (right) lenses. Pixels which are removed from the source recon-struction by the cleaning process on the source plane are shown as filled circles, while those data which contribute to the final model are hollow circles. The discarded points originate entirely from the inner end of the radial arc 143 6.24 (left) A small SIS, marked M3, is added near the centre of the cluster at position (0'.'50, -4'.'05). (right) The original model without the extra SIS. The highly magnified portion of the radial arc significantly changed 145 7.1 Each value of the mass parameter OCD of the cD galaxy in MS 1455 admits a small range of source redshifts at which a radial arc forms while maintaining the position of the tangential arclet. At the same time, OCD determines the mass-to-light ratio T of the cD galaxy through (7.2), plotted as a dashed line. The point marked A is the PID+cD model of MS 1455. The models at points B and C where the radial arc fails to form, are shown in Figure 7.2 157 X 7.2 (left) Model B from Figure 7.1 in which the radial arc fails to form because the source is too near, (right) In model C, the cD galaxy is too massive for a source at redshift zs = 0.825. The two images on either side of the radial critical line are disjoint, and a radial arc does not form. Small squares mark the location of the radial arc, A l ; the small circle marks the location of arc A3. The position of the tangential arclet is maintained in both cases 158 Acknowledgment My thanks are extended first to my supervisor, Greg Fahlman, for humouring my mathematical disposition, for helping me to effectively communicate my ideas and results, and most impor-tantly, for showing me how to use mathematics to do science rather than adapting science to fit the mathematics. I would also like to thank my committee members, Paul Hickson, Brian Seymour, Bill Unruh, and Matt Yedlin, for focusing their diverse backgrounds on the problem of gravitational lensing. This research would not have been possible without the computing facilities in the Applied Mathematics Lab established by Uri Ascher and the Institute of Applied Mathematics. This powerful and friendly computing environment does not simply run jobs faster, but it allows students to ask new questions not even considered in the past. I am also grateful for the assistance of Judy Maxwell who gave me access to all the facilities of the IAM without so much as a second thought. I could not have got this far without the continued interest and enthusiasm of my Mom, and without the support, encouragement, coaching, advice, and insight of my Dad, not just during this thesis, but always. My thanks go to all my friends and collegues in the Math Department, the IAM, and in the AML. In particular, I'd like to thank John Stockie for teaching me the ways of the emacsen, and Ray Spiteri for reminding me how simple this game can be, if you let it. Words are barely sufficient to describe how much love, encouragement, and support I've received from my wife, Margaret. Without her, this thesis would have been infinitely more difficult. I gladly pass onto her the role of graduate student. xi Foreword The Scientific Method Valid, reproducible data taken in context is truth, and the simplest theory that agrees with all the data and permits predictions of reality is the best explanation of truth. Alfred A. Brooks (1997) xii Chapter 1 Introduction A planet, a star, or a galaxy drifts through space like an ocean liner on calm seas. A ray of light, however, is tossed like a leaf floating amongst the ripples. The nature of these ripples is governed by the Theory of General Relativity, which predicts that rays of light are deflected as they pass by massive objects. Through this interaction, light from a distant source may be deflected by an intermediate object onto a serendipitous observer. This behaviour is known as gravitational lensing. In the case of strong gravitational lensing, the subject of this thesis, two or more rays of light diverging from a luminous source may be focused through a gravitational lens onto an observer. This observer sees multiple images of the background source, illustrated in Figure 1.1. Examples of strong gravitational lensing where a single background source appears multiple times through the lens require enormous mass distributions. A galactic nucleus is responsible for the lensing observed in the Einstein Cross 2237+030 (Yee, 1988). The deflecting mass distri-butions we consider in this thesis are on the scale of clusters of galaxies. A significant fraction, 90% or more, of these galaxy-clusters is thought to consist of dark-matter, non-luminous mat-ter which cannot be observed directly. This dark-matter still generates gravitational fields, however, and the light which is bent through a gravitational lens does not distinguish between the luminous matter and the dark-matter. The goal of this thesis is to describe a method of "inverting" such a lens, finding the mass distribution of the lens which produces an observed pattern of lensed images. This total mass distribution can then be compared to the distribu-tion of visible stars, gas, and galaxies to confirm, or refute, the presence of dark-matter in the 1 Chapter 1: Introduction 2 S Figure 1.1: Schematic, two-dimensional gravitational lens. Light from a background galaxy S bends around a massive lens M and focuses on the observer at 0. Multiple, magnified images of S are seen by this observer, appearing at different locations on the sky. galaxy-cluster. While detailed accounts of gravitational lens optics, modeling, and reconstruction are de-scribed in later Chapters, what follows immediately is a simple description of the effects we will be exploring and modeling. The appearance of a gravitational lens on the image plane can be modeled and simulated to arbitrary precision when the mass distribution of the lens on the deflector plane and the intrinsic appearance of the background source on the source plane are known. We generally take the image plane to coincide with the deflector plane, as this is the plane on which the observations are made and, as we shall see below, the plane at which the deflection occurs. It is only a matter of programming a computer to ray-trace the lens to produce a picture of the lensed background source. Unfortunately, just the opposite happens in nature. We observe only the lensed features, and must deduce both the mass distribution of the lens and the intrinsic appearance of the background source of light. It is not surprising that there is Chapter 1: Introduction 3 generally no unique solution of this problem, but rather families of source-deflector-observer optical systems which can reproduce observations of gravitationally lensed features. We need not throw up our hands in despair, however, because we do know the physics of lensing. We seek only the parameters which describe the particular lens at hand, not the physical explanation of this peculiar astronomical behaviour. Based on our knowledge of lens optics, we are able to construct sufficient parametrized models of gravitational lenses. In realisations of lensing systems, beams rather than rays of light pass through the lens. In addition to being deflected, the cross-sections of the beams are distorted. The appearance of the lens is characterised by the locations of the critical lines in the lens plane, the loci of points at which the cross-section of the beam has been squeezed shut. The light which is forced through this 0-area window is extremely (infinitely) amplified. Tracing the critical lines back to the source plane locates the caustics of the lensing system. Any light which is emitted in the vicinity of a caustic will appear to the observer highly magnified and distorted. Because of the great distances between the source, deflector, and observer, often only this extremely amplified light can be observed. A simulated gravitational lens similar to those which appear throughout this thesis is shown in Figure 1.2. The Figure mimics the observations one would obtain with a CCD. The picture is a "negative" of the observations, darker regions representing brighter light and white regions representing blank sky. The grayscale objects in the Figure are three images of a uniform disk of light hidden behind the lens on the source plane. The small circle marked SI near the centre of the Figure indicates how the source would appear in the absence of lensing. The lens in this model is a single large elliptical mass distribution. We do not include in the simulation any light which may come from the lens itself. In fact, we postulate that the major component of the lensing systems described here is a dark-matter halo which emits no radiation. The centre of the distributions, marked M l , is near the centre of the frame, (O'.'0,0''0). There are two critical curves in the simulation, plotted as dashed lines. Near the larger, outer curve, beams originating on the source plane are stretched tangentially. Images which appear on or near this line, like the feature labelled A 3 are known as tangential arcs. The caustic associated with the Chapter 1: Introduction 4 critical line is the "astroid" dotted curve. Because of the great distortion introduced by the lens, the envelope between the observer, the tangential critical line, and the tangential caustic is turned "inside out" producing a region of odd image parity inside the critical line. The origin of arc A3 can be traced to the region of the source disk which lies under the astroid caustic on the source plane. The arc appears very bright because of the extreme magnification. Figure 1.2: A simulated gravitational lens similar to those which appear throughout this thesis. Lensed features are represented by the grayscale objects. The centre of the lensing mass is marked M l . The background source, here a uniform disk, is indicated by the hollow circle marked SI. Dashed lines follow the critical lines of the lens, while dotted lines follow the corresponding caustics on the source plane. Chapter 1: Introduction 5 This simulation contains a second critical line nearer to the centre of the lens. Beams in the vicinity of this critical line are squeezed tangentially while being stretched radially. Lensed images appearing near such a critical line are radially distorted, and are referred to as radial arcs. These are quite rare, because their existence is sensitive to the form of the mass distribution, and the arcs may be obscured by luminous cD galaxies which often inhabit the centres of galaxy-clusters. The caustic associated with this critical line is the ovoid dotted line in the Figure. Radial arc A l is formed by the upper-right portion of the source disk which lies under this ovoid caustic. The great radial distortion of the image is obvious. Because any gravitating mass looks and feels like a point mass at the centre of the dis-tribution to a distant object, there is generally one gravitationally lensed image outside the network of critical lines. This image is, for the most part, the beam of light emitted by the source and bent around the intermediate mass onto the observer. The image is still distorted and amplified, but not to such an extreme as images which appear on the critical lines. Arc A2 in the Figure is such an image. The size of this type of lensed image varies greatly depending on just how far from the lens centre and critical lines the image appears. They are referred to as arclets in this thesis, to distinguish them from tangential and radial arcs. The thesis develops as follows. In Chapter 2, we review the history of gravitational lens analyses, highlighting the significant contributions to the development of the field and also the observational successes. Methods applicable to strong gravitational lensing are discussed in some detail. The optics of gravitational lenses, the formalism we follow, and the parametric lenses involved in this work are introduced in Chapter 3. In Chapter 4, we describe the lens modeling software, lens, which produces the simulations throughout the project and this thesis. A new lens inversion scheme is introduced in Chapter 5. We discuss the motivation for a new scheme, some of the problems encountered, and our strategy for overcoming these problems. It is critical that any work making predictions of astronomical events be tested on actual observations. In Chapter 6 we invert two observed lensing systems. The first, a collection of five arcs observed in the galaxy-cluster MS 2137—23, has already been well studied by Mellier et al. (1993), Miralda-Escude (1995), and Hammer et al. (1996). We use this lens as a test Chapter 1: Introduction 6 of the inversion algorithm itself, for our second lensing system, a collection of three arcs in the cluster MS 1455+22 has not been modeled before. We wish to separate the results into those which are due to lensing and those which may be due to defects in the inversion scheme. Finally, in Chapter 7, we summarize the success of the new inversion scheme, and discuss several directions this research could follow in the future. Chapter 2 Review of Gravitational Lensing The field of gravitational lensing is not yet 100 years old, for it was born in Einstein's Theory of General Relativity (1915). It is only in the last 20 years that lensing has become a familiar topic to astronomers. Strong gravitational lensing of a radio source was detected in MG 1131 + 0456 (Hewitt et al., 1988). With the advent of highly sensitive optical instruments, the arcsecond scale of strong lensing could also be resolved in the optical. The field has grown enormously in the last five years with the introduction of the Hubble Space Telescope (HST). We concentrate on strong gravitational lensing of optical background sources in this work. The development of the field is outlined in this Chapter. Particular attention is paid to the more recent work on strong gravitational lensing. We describe in some detail the lens inversion schemes currently in use throughout the lensing community, to emphasize their often subtle differences, and also to make clearer the new features of the lens inversion scheme described in Chapter 5. 2.1 Before 1937: Hypothesis Even before the beginning of the 20th century, it had long been suspected that light responded to the force of gravity. As most books about General Relativity cite, but only a few like Rindler's Essential Relativity (1979) bother to show, Newtonian mechanics predicts that when an object of mass m grazes a mass M with impact parameter r, the object is deflected in the plane of the orbit by an angle a — - 2 G M m ^ . 7 Chapter 2: Review of Gravitational Lensing 8 By equating mass and energy, this result suggests that an a ray of light will be deflected through an angle _ _2GM r c 2 r - 2 This is Einstein's 1911 prediction, in which he considers only geometric arguments. He im-mediately saw the implications of the deflection of light and proclaimed "[i]t would be a most desirable thing if astronomers would take up the question here raised." (Einstein, 1911) The geometric arguments are revisited in the Theory of General Relativity (1915), where Einstein includes the effects on the ray of light of the potential well around the mass M. The revised prediction is exactly twice the earlier one. In 1919, Eddington's observations of the deflection of star light past the eclipsing Sun (Dyson et al., 1920) supported the new theory and established General Relativity as the preeminent theory of gravitation. 2.2 The Period 1937-1979: Analysis The phenomenon of light bending around massive objects remained yet another Relativistic Gedankenexperiment until Zwicky suggested an astronomical application: gravitational lensing of very distant objects by intermediate galaxies (1937). The arcsecond scale of the effects of lensing were beyond the observational capabilities of the time, however. Despite the inability of astronomers to observe gravitational lenses, analysis proceeded awaiting the first observed lens. Refsdal (1964a; 1964b) showed the cosmological applications of gravitational lenses, in particular, how the Hubble constant could be determined from any measurable time delay in the signals received from different images of the same background source. To provide a more realistic model of the galaxies thought to be responsible for strong gravitational lensing, Bourassa, Kantowski, and Norton (1973) generalized point-mass lenses to ellipsoidal mass distributions. The work is (perhaps) the first to treat the lens as an optical system showing distortion, caustics, critical lines, and other characteristics of laboratory optics. A beautiful analysis of the gravitational deflection about elliptical lenses is given by Schramm Chapter 2: Review of Gravitational Lensing 9 (1990). This work forms the basis of the numerical models produced here, and is discussed further in Section 3.1.3. 2.3 The Period 1979-1991: Discovery It had been suggested as far back as 1965 that quasars were actually gravitationally lensed and magnified galaxies (Barthony, 1965). Astronomers then and now (Kochanek, 1995) concen-trated their search for gravitational lenses on quasars. In 1979, Walsh et al. observed the first gravitational lens 0957+516 A / B , a doubly-imaged quasar. Once lenses were shown to exist in the sky as well as on paper, the race was on to find the next. As more and more multiply imaged systems were discovered, it became apparent (with hindsight) that lensing needs a great amount of mass, and one place to find sufficient mass is in clusters of galaxies. Lens-hunters turned their telescopes to the Abell catalogue of galaxy clusters (Abell, 1958; Abell et al., 1989). Lynds and Petrosian (1986) were rewarded with a new kind of lensed image in the cluster Abell 370: an arc, the highly distorted and magnified image of a background source lying near a caustic of the lens. Since that time, dozens of lensed objects have been observed (Surdej and Soucail, 1993). Some of the most spectacular arcs are seen in Abell 370 (Lynds and Petrosian, 1986; Hammer and Rigaut, 1989; Grossman and Nararyan, 1989), Abell 2390 (Pello et al., 1991; Kassiola et al., 1992), and Cl 2244-02 (Bergmann et al., 1990; Hammer and Rigaut, 1989). The Hubble Space Telescope (HST) has delivered unprecedented images of the lens in Abell 2218 (Kneib et al., 1996), which contains dozens of strongly lensed images, and in MS 2137 (Hammer et al., 1996). High resolution HST images of Abell 317 (Smail et al., 1996) may have revealed a radial arc beyond the reach of current ground based instruments. 2.4 After 1991: Lens Inversion After the excitement of discovery, astronomers paused for a moment to evaluate their good fortune. An excellent review of the state of lensing as of 1992 was produced by Blandford and Narayan (1992), which is responsible for this author's interest in the subject. At the same Chapter 2: Review of Gravitational Lensing 10 time, Schneider, Ehlers, and Falco published the seminal work on the theoretical aspects of lensing (1992). This work collects, and even establishes much of the notation used throughout the lensing community, and has risen to the status of being simply "SEF92," much like the relativist's "MTW" and the mathematician's "Rudin." Other astronomers, intrigued by this newly discovered phenomenon, turned to the question of how gravitational lensing can be used for further exploration. Three branches of lensing analysis have developed, each concentrating on a different scale of lensing. Historically, weak lensing and then microlensing developed from strong lensing, the regime at the heart of this report. We briefly describe these fields below. 2.4.1 Microlensing Microlensing is a vastly different application of the deflection of light than its two stronger cousins. Microlensing is used to search for massive compact halo objects (MACHOs) in our own Galaxy. These non-luminous, Jupiter-sized masses may be the Galaxy's missing dark-matter, if their numbers are large enough. When a MACHO passes in front of a background star, a transient event lasting perhaps a few days, the magnification introduced by lensing produces a well-defined increase followed by a decrease in the brightness of the background star. The size and shape of the star's light curve reveals the mass of the MACHO once the distances to the star and MACHO have been deduced. Extremely coincidental alignment of the star, the MACHO, and the observer is required to produce a measurable effect. The improbability of this alignment occurring is overcome by turning the telescope toward regions of the sky densely populated with stars. The MACHO Collaboration (Alcock et al., 1997) focused first on the Large Magellanic Cloud (LMC) and more recently on the galactic bulge for fields of background stars. This group has identified several dozen microlensing events which very closely match the predictions. Because of the short crossing time, MACHO Project at wwwmacho. anu. edu. au provides an electronic microlensing notification service to observers who can immediately turn their telescopes towards the event. Other groups searching for microlensing within the galactic bulge, the LMC, and the Small Magellanic Cloud include the Optical Gravitational Lensing Chapter 2: Review of Gravitational Lensing 11 Experiment (OGLE) (Paczynnski et al., 1995) and the EROS project (Beaulieu et al., 1995; Renault et al., 1997). 2.4.2 Weak Lensing Weak and strong gravitational lensing occurs on Mpc scales where foreground galaxies act as lenses and background galaxies are lensed. Before one can explore details of the background source or the intermediate space, one needs a highly resolved model of the lens itself. The light from every source behind the lens plane is distorted by some amount as it passes through the lens. Only those background galaxies lying close to the observer-lens axis are strongly lensed. The elliptical appearance of each background galaxy is slightly elongated tangentially about the centre of mass of the lens. While the weak distortion of any single galaxy cannot be measured, by measuring the shapes in a statistically large sample of background galaxies, a coherent shear field can be extracted. This signal can be inverted to produce a mass distribution of the intermediate lens (Kaiser and Squires, 1993). The surface brightness of the source is conserved under gravitational lensing. Therefore, as first suggested by Tyson et al. (1990) and Blandford et al. (1991), the surface brightness I' observed at coordinate x is the intrinsic surface brightness / of the source at position x — V $ , where

. Looking at a single lensed image reveals nothing about the lens because the intrinsic shape of the background galaxy is unknown. The background population of galaxies, however, come in all shapes and at all position angles, so statistically, where brackets indicate the average over the sample. By measuring many weakly lensed images in a small region of the sky, an approximation for A<&, and hence the surface mass density, in that neighbourhood can be found. In applications of this inversion-via-distortion algorithm, like those used to invert MS 1455 and CL 0016+16 (Smail et al., 1995), MS 1224 (Fahlman et al., 1994), and RXJ 1347.5-1145 (Fischer and Tyson, 1996), the image is subdivided, each cell containing enough background galaxies to extract an estimate of A$ . The values are smoothed and produce a one parameter map of the surface mass density. This free parameter can be regarded as a uniform mass sheet over the entire region, which has no effect on the lensing properties of the images, but obviously contributes to the mass. Wider and wider fields-of-view are being sought in an attempt to map all the way to the "edge" of the cluster to try to remove this degeneracy. The inversion scheme is sensitive to edge effects, which very wide images might suppress as well. More recent reconstruction schemes based on ellipticity and shear try to extend the tech-nique into the non-linear regime of strong lensing (Schneider and Seitz, 1995; Seitz and Schnei-der, 1995), but the results are still preliminary. As the strength of the lensing increases, the constant A $ approximation requires finer and finer resolution. As well, the distorted images (ei) = (e2) = (eie2) = 0 Chapter 2: Review of Gravitational Lensing 13 are less faithfully represented by an equivalent ellipse. The arc A2 in Figure 1.2 is very nearly elliptical, but arcs A l and A3 are not. Kaiser (1995) has made some headway by finding a combination of third derivatives of the potential <&. The integration of these third derivatives to find the underlying mass distribution $ cannot be carried across the critical lines of the lens. This requires determining the parity of each region of the image plane by some other means before the surface mass distribution can be recovered. Another coherent signal that one can use to reconstruct the lens mass distribution comes from the magnification field about the lens. Gravitational lensing preserves the surface bright-ness of the background sources. The magnification factor is the ratio of the area of a lensed image to the area of the image seen in the absence of lensing, where the area is defined by the region on the sky inside a fixed isophotal level (Broadhurst et al., 1995). Of course, the area of the image in the absence of lensing cannot be measured. A coherent signal can be extracted, however, by first dividing the observations into a grid, and then comparing the area of all back-ground sources above a fixed isophote in each cell to the average area of background galaxies observed in fields far outside lensed region. As Broadhurst (1995) recognizes, measuring the areas of very small background galaxies is difficult, even at the O'.'l resolution of the HST. With the same set of observations, one can consider the number of galaxies visible in each region of the sky. Background sources too faint to be seen in the absence of lensing may be magnified into view, changing the number density of galaxies on the sky. This suggests that by averaging over regions of the sky, the magnification in each region can be determined from the number density of galaxies. However, the number of galaxies per magnitude per unit solid angle varies with redshift because of the evolution of the universe and of the galaxies themselves. The number density of the faintest galaxies in the observations may be different than the number density of galaxies just below the observable limit, which become visible because of lensing. This is essentially the argument that resolves Obler's Paradox: the number of galaxies per unit solid angle per unit distance is not constant. These competing effects, fewer background galaxies but deeper views because of magnification, cannot be decoupled without a theory of galaxy evolution. The effect is coherent across field-of-view, though. By comparing the number Chapter 2: Review of Gravitational Lensing of visible galaxies per solid angle per redshift between regions near the lens and far from the lens, changes in the number density can be measured and a coherent magnification signal can be extracted. 2.4.3 Inversion of Strong Lenses The successful strong lens inversion schemes now in use fall into two broad categories: those which iteratively process lens observations to whittle the model down to its simplest form, and those which build up the model and observed arcs based on the geometry and optics of the lens. To each inversion scheme is associated a spectacular gravitational lens system. There are lenses enough to go around, so that different groups and their approaches are rarely compared on the same lensing system. Not surprisingly, the advent of the HST and its high resolution images has allowed for several very successful lens inversions, particularly in CL 0024+1654 (Dell'Antonio et al., 1996; Kochanski et al., 1996), Abell 2218 (Kneib et al., 1996), and MS 2137 (Hammer et al., 1996). LensClean The LensClean algorithm developed by Kochanek and Narayan (1992) is particularly well-suited for inverting extended lensed images, like radio rings. LensClean is a synthesis of the common image-cleaning algorithm Clean (Hogbom, 1974) and the "Ring Cycle" algorithm previously developed by Kochanek et al. (1989). The Clean algorithm iteratively subtracts from an astronomical image the beam-convolution of the brightest peak in the image. The point-spread function used to clean the images is tailored to remove the "side-lobes" often found in radio observations due to constructive interference of the radio signals. The procedure is repeated on the next brightest peak until the residual is comparable with the background or noise. To reconstruct the cleaned image, each peak is replaced by a (^ -function of the same flux, and the collection of ^-functions is convolved with an ideal beam (generally Gaussian in profile). Finally, the residual is added to produce the "cleaned" image. Chapter 2: Review of Gravitational Lensing 15 The Ring Cycle attempts to model the mass distribution of the lens by exploiting the fact that the surface brightness of the source distribution is unaltered by gravitational lensing. Therefore, pixels in the observed lensed features which are traced to the same source point must have the same surface brightness, within the level of the noise in the observations. The Ring Cycle proceeds by tracing a lensed pixel back to a point on source plane under lensing through the current model. Because of the form chosen for the model lens mass distribution, in fact, two-dimensional potentials for singular isothermal spheres with an external shear, all of the solutions of the lens equation stemming from this source point can be predicted analytically. In this way, all the lensed pixels arising from this source point can be found including, of course, the pixel used to find the source point originally. The lens parameters are then adjusted using the simplex method to minimize the differences in the surface brightness between each lensed pixel and the average surface brightness of the collection. At the same time, the best value for the surface brightness of the source point is deduced. Kochanek et al. (1989) do not discuss how the misfit minimisation is carried out or the computational power needed. They apply this algorithm with good success to the radio ring MG 1131+0456. Rings like this one contain many multiply-lensed pixels, providing many independent constraints on the lens model and the source structure. These two algorithms are combined and enhanced in LensClean. The brightest lensed pixels are traced back to the source plane, and then traced forward again through all the solutions to the lens equation describing the current model. Again, finding these solutions is possible because of the choice of the lensing mass potential. This larger collection of rays is convolved with a beam and compared with the observations. The exact position of the parent source point and its surface brightness are set to minimize the difference between the observed and predicted data. The measurement this time takes into account the magnification of the source due to lensing. This point is added to the source surface brightness distribution, and its beam-convolved lensed children are removed from the observations. Continuing this way generally cannot remove the entire lensed feature, precisely because the model may be incorrect. The cleaning loop is halted when the residuals are no longer decreasing, leaving a x2 Chapter 2: Review of Gravitational Lensing 16 statistic measuring the goodness-of-reconstruction based on the current lens model. This cycle is repeated, adjusting the model parameters to find a global minimum in the X 2 statistic produced by each cleaning. The final result is a best set of lens parameters and the corresponding source distribution. The algorithm also returns error estimates on each of the parameters, interpolated from the collection of intermediate calculations towards optimi-sation. Kochanek and Narayan (1992) apply this LensClean algorithm to the radio ring lens PKS 1830—211 to produce a high resolution, well-fitting reconstruction. Chen et al. (1995) "LensClean" the radio ring MG 1131+0456, the original Ring Cycle test, with great success. Neither of these groups discusses the number of iterations or the computation time needed to produce their optimal models. L e n s M E M The Clean algorithm is better suited for beam-convolved observations of point sources. Back-ground sources of gravitational lensed images are often extended galaxies or quasars, however. Wallington et al. (1996) replace the Clean cycle of the LensClean algorithm with a cycle using the maximum entropy method (MEM) (see Narayan and Nityananda (1986) for an excellent review of the MEM.) In an under-constrained inverse problem like lens mass reconstruction, there is generally an infinite family of solutions which reproduce the data to within the errors. The MEM attempts to find the particular solution which encodes the minimum amount of information, the solution which is least biased by our pre-conceived notions of what the solution should look like. This is also, in some sense, the most probable solution. From the observations Dij and for each set of lens model parameters, LensMEM finds the source plane surface brightness distribution Ski and corresponding beam-convolved lensed image Iij which maximizes the function This function can be related directly to Bayes' Formula: The probability P(S\D) of finding a particular source distribution S given the observations D, is proportional to P(S) • P(D\S). (2.1) Chapter 2: Review of Gravitational Lensing 17 The first term is maximized by finding the source S with the maximum entropy, the first term of (2.1). The second term, the probability of predicting the data D given a source distribution S, is measured by the misfit between the observations and the data in the image plane, which is the second term of (2.1). The relative importance of these two terms can be controlled by the Lagrange multiplier A. As Narayan and Nityananda (1986) demonstrate, there is no one correct form for the function J , but the one given here is easily interpreted. Because of the multiple lensing of source pixels, convolution of the image, and the finite resolution of the data, it is unclear how many degrees of freedom are available, although clearly it is less than the number of lensed pixels. Hence defining a target value for a x2 misfit is impossible. Wallington et al. (1996) contend that performing 100 LensMEM iterations will neither over-fit poor/sparse data nor under-fit adequate data. In each of these iterations, the model parameters are varied, using a conjugate gradient approach, to increase the function J. The resulting model and source distribution come with a goodness-of-reconstruction statistic, •Aoo-The key to the LensMEM, and also LensClean algorithms is the ability to find all the image pixels where a given point in the source distribution will appear. To extend the limited choice of lens models demanded by the Ring Cycle, LensMEM uses a seemingly elegant approach, although some obvious limitations are described below. The corners of each image pixel (i,j) are projected back to the source plane, reconnected, and the overlap with each source pixel (k, I) is calculated, producing a four-dimensional array of weights Wijki- From this array, one can calculate the number of images NM of source pixel (k, I) ATfcf _ £ t j ijkl ^ ^ _ a r e a Q £ g o u r c e pixel the magnification in image pixel Pij = \^^M^12M?j ? Aj = area of image pixel and also the lensed image itself T _ T,kl wijklSkl lij — v-> Z f^cZ Wijkl Chapter 2: Review of Gravitational Lensing 18 While this is an elegant approach, it is limited in the resolution it can achieve, unless some even more elegant programming has been done. In tests described by Wallington et al. (1996), the source plane is resolved on a 128 x 128 grid and the image plane on a 64 x 64 grid. The array Wijki, while taking only 0(N2) operations to initialize, appears to occupy some 2 2 6 memory locations or roughly 1 Gbyte of memory! The algorithm was tested on the radio ring lens MG 1654+134, and reproduced the LensClean solutions produced by Kochanek (1995). This approach provides some measure of the errors in the model parameters and can be interpreted as a most probable solution, but once again, there is no discussion of the practical implementation of the algorithm. The Toulouse Approach In the LensClean and LensMEM inversion schemes, ideally one feeds in a data set, processes the data through a programmed set of routines, and produces the inverted lens mass distribution. The Toulouse school of gravitational lens inversion approaches the problem from the opposite direction. A parametric deflector mass distribution is manipulated, based on knowledge about the behaviour of gravitational lenses, to mimic the observed lensed images. This is the path most closely followed by my own work. Mellier, Fort, and Kneib (1993) model the lens in MS 2137—23 with a pseudo-isothermal elliptical mass distribution (PID) centred on the large cD galaxy assumed to lie at the centre of the cluster. A source which reproduces the observed giant arc gives rise to two additional images, which appear to have observed counterparts. Based on this triply-imaged source, the model parameters are adjusted to best reproduce the positions of the lensed images. The optimisation is done by finding three sources, one for each lensed image, and demanding that they coincide on the source plane. An additional cluster galaxy is added to the model to enhance the structure of the giant arc. Again, the model parameters are adjusted. Finally, a second source is added, without adjusting the lens mass itself, to model the observed radial arc. While there is more than enough freedom in the model to reproduce this independent lensed feature, this second source introduces additional lensed features. These features closely match observed Chapter 2: Review of Gravitational Lensing 19 arcs in the lens. As the group contends, there remain small differences between the positions and shapes of all the lensed features, but the consistency of the model is compelling. This educated trial-and-error approach is used by Kneib et al. (1995) to model the lensing mass potential in Abell 2218. From ground-based observations, the multiply-lensed features are used to tune model parameters of PID masses originally placed coincident with the dom-inant luminous cluster galaxies. The best fit is found by minimizing a x 2 statistic measuring differences in background source models. Kneib et al. (1996) were offered the unprecedented opportunity to turn their method loose on the incredible HST images of Abell 2218. Based on the earlier model of four PID masses, this model contains an additional 30 cluster galaxies. The position, orientation, and ellipticity of each is matched to the observations, leaving only the shape of the mass profile and total mass of each galaxy undetermined. A weak-lensing analysis is carried out on 235 lensed arclets of sources located behind the cluster to ascertain the large-scale distribution of mass in the cluster. Parameters describing the 30 cluster galaxies are fit through seven multiply-imaged sources. A refinement to the x 2 statistic is made to take advantage of the large number of lensed arclets. These arclets are used to produce a low-resolution map of the total shear in the neighbourhood of Abell 2218. A new x 2 measures (i) differences in elliptical source distri-butions which are multiply-lensed (as in earlier work), and (ii) differences between the shear produced by the modeled lens mass distribution and the shear interpolated from the appearance of the 235 distorted arclets. This huge number of parameters and equally huge number of con-straints produces a lens mass distribution at the x2-minimum. Because the goodness-of-fit of the strongly lensed (multiply-imaged) features in made on the source plane, the magnification of the gravitational lens is not exploited in any way. The Tyson Group Working with HST images of the lensing cluster CL 0024+1654, the group lead by Tyson has implemented several inversion techniques with the flavour of all the approaches described above. Chapter 2: Review of Gravitational Lensing 20 Colley et al. (1996) "unlens" the cluster by concentrating on four strongly lensed images of a complex background source. Starting from a base mass distribution describing the 48 cluster galaxies, one large pseudo-isothermal dark matter halo is added. A weak lensing analysis of some 1200 observed arclets supplies initial values for the model parameters. Because each of the strongly lensed images shows (supposedly the same) structure of the background source, three prominent features are traced back to the source plane. Colley et al. vary the model parameters of the dark matter halo and 5 significant cluster galaxies until the pre-images of the lensed features coincide on the source plane. The distance (in pixels on the source plane) between the features traced back through the lens through the various lensed images acts as a measure of the goodness-of-fit, the smaller the better. Several different initial values for the dark matter model are used, each converging to roughly the same model. The best model reproduces the morphology of the source in each lensed image. In an attempt to remove the "tuning" factor of this method, Dell'Antonio et al. (1996) and Kochanski et al. (1996) perform a massive x 2 minimisation. Only the results of these models have been released. The method, according to Kochanski, postulates a pseudo-isothermal sphere for each of the 110 cluster galaxies and a large elliptical dark matter halo. The distribution of light on the background source plane is modeled as a (large) number of uniform elliptical disks, with free position, orientation, ellipticity, and surface brightness. Some 400 parameters are needed to describe the model. For each set of parameter values, pixel-by-pixel comparisons between the observed and modeled strongly lensed features are tabulated, and a global minimum is sought. The nature of the minimisation routine is unclear, although Kochanski indicated hundreds of thousands of x 2 tests were performed. The final results of this simulation are among the most detailed lens models constructed to date (Lucent Technologies, 1997). The simulated lensed features reproduce the observations very closely. Reconstruction of MS 2137-2353 High resolution HST images of the cluster MS 2137—2353 show several strongly lensed features. Similar structures within the lensed features suggests that 5 of these giant arcs, including a radial Chapter 2: Review of Gravitational Lensing 21 arc, are images of only two background sources. Hammer et al. (1996) use the multiply-imaged sub-structure in the sources to reconstruct the lensing mass distribution. Bright "knots" within each of the lensed arcs are identified and traced back to source plane assuming a j3 profile lens model, where the three-dimensional density is given by The free model parameters are varied systematically. For each model, a "self-similarity param-eter" measuring the degree of coincidence of the (supposedly) common peak in the pre-images on the source plane is calculated. The set of parameters producing the most "self-similar" collection of sources is chosen. Having fit the bright knots, the entire lensed images are mapped back to the source plane, combined, and then traced forward to test for spurious arcs. The simulations very closely match the observations. The final model is similar to the model of Mellier et al. (1993) described above. With the huge improvement in resolution, models for the surface brightness distribution of the sources are also produced. The interest in gravitational lensing is growing exponentially. Higher quality ground-based observations, better data processing routines, and continued access to the HST is producing more and more cases of gravitational lensing. The collection of available data also continues to grow, allowing researchers to test new lens inversion schemes. We have tried to include the important contributions to the field to date, particularly those focusing on strong gravitational lensing. We may inadvertantly miss important results currently appearing in literature. Chapter 3 Gravitational Lens Optics In this Chapter we explore the mathematics, physics, and optics of gravitational lensing. After deriving an expression for the deflection of light from the Theory of General Relativity, we apply the result to the case of lensing and produce the lens equation which lies at the very heart of the subject. From this simple equation we can deduce the nature of the lens. Armed with these relations, we consider a small collection of lens models which form the basis of our simulations and inversions. By understanding the characteristics of these parametrized lenses, we are able to make the educated guesses needed to model and then invert observed galaxy-cluster lenses. 3.1 Deflection of Light At the centre of General Relativity is the tight relationship between curvature and energy, math-ematics and physics. Einstein postulated that the 4-dimensional manifold describing spacetime curves in response to the presence of matter and energy. Bodies carrying the matter and energy react and evolve according to the equations of motion which, in turn, alter the geometry of the spacetime. In a curved spacetime, freely falling bodies no longer travel in what are classically called straight lines; in respect to Newton's F = ma, deviations from straight paths are traced to a force, the force of gravity. 22 Chapter 3: Gravitational Lens Optics 23 3.1.1 Perturbations of Flat Spacetime The differential geometry describing the mathematical side of Relativity is built upon the spacetime metric defining a metre stick at every point: 3 dr2 = ^ 2 giivdx^dx" = g^dx^dx" . The second expression uses the Einstein summation convention of summing over all values of each repeated index. This relation gives the arclength dr in terms of coordinate differentials, just like ds2 = dx2+dy2 in the Cartesian plane. Einstein's great contribution was in establishing a relationship between differential geometry and physics in the array of equations Rrw - ^Rg^ = ^f-T^ . (3.1) The left-hand side is composed of combinations of the metric g^ and its first and second derivatives wrapped up in the Riemann Curvature tensor. The right-hand side of (3.1) is the stress-energy tensor. Its components encode the energy density, mass density, and momentum flux about each point of spacetime. Equating these two links the energy density at every point and the curvature of its neighbourhood. Another factor Kg^v can be included in the left-hand side to reflect a non-zero cosmological constant A, but we assume A = 0 in this thesis. A solution to the Einstein field equations (3.1) yields a metric g^ to describe the physics of spacetime. There are exact solutions to Einstein's equations, like the Schwarzchild blackhole solution, but all depend on high coordinate symmetry. General solutions are beyond reach. Numerical solutions to the Einstein field equations are very difficult because the coordinate system evolves along with the objects living in that coordinate system. To produce a mathematical description of gravitational lensing, we first consider empty space. The metric is the Minkowski metric, the relativistic counterpart to the Theorem of Pythagoras: dr 2 = (dx0)2 - (dx1)2 - (dx2)2 - (dx3)2 • = t]iiVdxildxv . Chapter 3: Gravitational Lens Optics 24 The 4x4 Minkowski metric rj^u has only diagonal elements {1,-1,—1,-1}. For simplicity, the coordinates are written (x° 2 CC j CC 2 cc 3) ~ (ct,x,y,z). Light, traveling at speed c, lives on the null cone r = 0. This empty space is flat and light travels in classical, straight lines. A stationary point mass M introduced at a fixed 3-dimensional position XM perturbs the local spacetime metric: VIJ.V -> 9fiv = n^v + , |/fy„| < 1 . The small perturbations h^ are determined by Einstein's equations (3.1) when the stress-energy tensor describes a point mass. This stationary point mass M exerts no pressure on its surroundings and T1*" has only one non-zero component T 0 0(!B) = Mc2(53(!E-a!M) . To first order in Einstein's equations reduce to a single relation A l 167TGM.3, , A/i M „ = ^ ° \ x ~ X M > with hfj,v — hpv — \nilvhaa. The simplest of Green's functions gives the solution to this Poisson equation, where h^iii — hSjn/ h(x) = - 5 , r = \x- xM\ • cz r In more physical units, ,rs n n ( M\ / l k m \ For light rays passing within a kpc of a galaxy-sized point mass of 1O 1 O M 0 , h ~ 10 - 6 and the first-order approximation is valid. 3.1.2 The Deflection Angle To calculate the angular deflection experienced by a ray of light as it passes through the local anomaly in the gravitational field produced by the point mass, we compare the incoming and Chapter 3: Gravitational Lens Optics 25 outgoing momenta. A photon entering the region with only an initial ^-component of momen-tum p3 picks up orthogonal components and leaves the region with momentum (Pf,p2,p3)- We assume p3 = p3 in the small \h\ approximation. The ray is deflected through a small angle ip)+3P2f P3f Here, i and j represent unit vectors in the 1- and 2-directions. We reserve the symbol a for the scaled deflection angle (3.10). The outgoing components of the momentum are simply = / ! ( £ ) * * ( 3 - 2 ) and similarly for p"j. The term dp1 /dx3 can be extracted from the geodesic equations, equations of motion describing the photon freely falling through the gravitational field: d2zCT 1 dh^ dx»dxv For the massless m = 0 photon, we cannot naively set pa = m(dxa/dr) as one usually does in switching to a Hamiltonian formulation of this system of second-order differential equations. Instead, we define an afline parameter A = r /m to parameterize the path of the photon, recalling that the arclength r of the photon's path is also 0. Then pa = dxa/dX and the geodesic equations reduce to ^ _ i a v 0 1 2 3 The incoming ray is traveling in the z-direction, so we approximate p in the right-hand side with the initial values p1 = 0, p2 = 0, and p3 = dx3/dX = u>, constant. Conservation of energy-momentum, paPo- = 0, demands p° = u. The term we are seeking in (3.2) can now be written dp1 _ dp^_dX_ _ dp1 1 dx3 dX dx3 dX u so that 5 = = r (—i+—j]dx3 7-00 \dxlt+dx*3) Chapter 3: Gravitational Lens Optics 26 At this point, we make the "thin lens" approximation which assumes the deflecting mass lies completely in the x 3 = 0 plane and that the photon is abruptly deflected as it passes through the plane. This is not an unreasonable approximation. First, the separation between the source, the deflector, and the observer is Mpc in scale while the deflecting mass has scale on the order of kpc. Second, the perturbation dh/dx1 falls off as 1/r3, and is non-negligible only in the neighbourhood of the mass M. Under the thin lens approximation, the deflection angle integrates to where r is the impact parameter of the incoming ray on the x3 = 0 plane as if no deflection occurs. The deflection is parallel to the impact parameter but with the opposite sign, and the ray is bent back towards the deflecting mass. This is the relativistic prediction which is twice the Einstein's original, geometric prediction for the deflection angle. The behaviour can generalized to more complicated distributions of deflecting matter. Un-der the thin lens approximation, the mass distribution is defined only on the 2-dimensional deflector plane. This is achieved with the surface mass density £ (£ ) , found by projecting the 3-dimensional volume mass density down onto the deflector plane. We introduce the 2-dimensional vector £ to represent the physical coordinates on the deflector plane. The deflec-tion angle experienced by a ray striking the lens plane at £ is the superposition of the deflections due to mass elements dM = E(£)d 2 ^ across the deflector plane: The sign of the deflection is now a matter of convention. There is an equally valid development of the deflection angle and subsequent gravitational lensing based on a scalar potential $(£) rather than the surface mass density £ (£ ) . We choose not to follow this route because we wish to base the our models on the form of S. The elegance of this approach, however, warrants at least a brief mention here. Exquisite detail can be found (3.4) in SEF92. The simplicity of the potential formulation of deflection is based upon two mathematical Chapter 3: Gravitational Lens Optics 27 identities. As V l n | £ | = £ / | £ | 2 , (3.4) can be written as a (£) — V $ ( £ ) where 4G Furthermore, as A l n | £ | = 27r<52(£), the surface mass distribution can be extracted from the potential through Poisson's equation: A * « ) =-2 -2 (0 • The lens equation and the network of critical lines and caustics it encodes have potential-based formulations as well, which will be briefly visited below. 3.1.3 Elliptical Mass Distributions The results in this report are based on lenses built by specifying the surface mass distribution and numerically simulating the lensing behaviour. The general surface mass distribution £ (£) in (3.4) is replaced by more symmetric versions, for several reasons. Spherical and elliptical symmetries greatly ease the computation of the deflection angle. More importantly, meaningful parameters, like the mass of the lensing galaxy, can be extracted from the resulting model by choosing physically motivated profiles for the mass distribution. Finally, specifying the parametric form of the lens may remove enough degrees of freedom from the model that it becomes predictive rather than simply adequate. We devote Section 3.4 to spherically symmetric lenses. Here, we consider deflection about mass distributions with elliptical symmetry. The analysis described below expands on the ideas presented by Schramm (1990). In spherically symmetric mass distributions, the deflection angle at radius £ from the lens centre is given by AC 1 Each annulus of mass £(£') 27r£'d£' contributes to the deflection of points at radius £ as if the mass was concentrated at the centre of the lens. This interpretation does not so easily extend to elliptically symmetric distributions, where shells of constant mass density are elliptical. One reason is that the distance between a point at £ and the elliptical shell is not well defined. Chapter 3: Gravitational Lens Optics 28 Consider an elliptical domain with semi-major axis a, semi-minor axis b. Over the semi-major axis, the surface mass density profile is £(£)• Along other rays from the centre to the boundary of the domain, the profile is scaled horizontally to extend the function over the entire region. A family of concentric, similarly oriented ellipses of constant density (Figure 3.1) is defined by + y - 1 , 0 < m < 1 (ma)2 (mb)2 These curves are known as homoeoids. The coordinates x, y are written with respect to the principal axes of the ellipse. The m = 0 homoeoid sits at the origin of the coordinate system; the boundary of the elliptical domain lies at m = 1. Figure 3.1: Parameter m between 0 and 1 describes a family of homoeoidal ellipses which cover the elliptical domain. At a point £ = (re, y) lying inside the elliptical distribution on a homoeoid at m', the contributions to the deflection of all the shells at m > m' exactly balance to zero, just as in the spherically symmetric case. The deflection at £ depends only on the mass of the shells at m < m'. The contribution to the deflection at the point £ due to the mass on any particular shell m depends on the "distance" between the point £ and the shell. To define this distance, we consider the family of ellipses confocal to shell m. Confocal ellipses (Figure 3.2) are a family of ellipses with parameter A > 0 defined through + (am)2 + A (6m)2 + A = 1 (3.5) Chapter 3: Gravitational Lens Optics 29 Figure 3.2: Parameter A > 0 describes a family of ellipses confocal to the central A = 0 ellipse. The confocal ellipses become more and more circular as the mass lying on the A = 0 shell appears more and more point-like. The value of A defines the distance between a point £ and the mass shell at A = 0. The confocal ellipse with A = 0 is the original m-homoeoid. As A increases, the confocal ellipses become larger and more circular as the mass on the the m-homoeoid appears to be more and more point-like. The distance from the homoeoid at m to the point £ is defined through the value of A describing the ellipse passing through £. Combining all the shells lying "inside" £, the x- and y-components of the deflection a are given by av 4G , p p 1 2 = ^ a b J o > S(am) 2-Krndm where a12 - (am)2 + A b12 = (bm)2 + A Chapter 3: Gravitational Lens Optics 30 1 _ x2 y2 p12 = a^+b^ with A defined through (3.5) and m0 = max This limit of integration chooses which homoeoids are "inside" the point £ and contribute to the deflection and those which are "outside" and have no effect on the deflection at £. Equivalently, m0 determines which homoeoids must be included in the numerical calculations and which can be ignored. 3.2 The Lens Equation We have seen how mass distributions alter the path of passing rays of light, deflecting those rays through a prescribed angle. The question remains, of all the rays emitted by a background source, which enter the observer's instrument? Or, more fundamentally, how does this lens focusl As described in this Section, the lens equation selects the rays which leave the source and are deflected through exactly the correct angle to strike the observer's telescope. We derive the lens equation based on geometric arguments. Consider a mass distribution which produces a deflection fi(£) at every point £ on the lens plane (Figure 3.3.) Light emitted at position r) on the source plane which strikes the lens plane at £ will continue on to the observer for every ray satisfying V = ^~Ddsa . (3.6) The distances Ds, Dd, Dds are angular size distances from the observer to the source plane, to the deflector plane, and between the deflector and source planes, respectively. The angular size distance Dd is defined so that an object of physical size £ on the deflector plane subtends an angle £/.E>d- This distance is not simply the line-of-sight (affine) distance. Instead, because of the expansion of the Universe "transverse" to the path of the ray, there is an extra redshift factor in this measure. In a homogeneous universe with density parameter f2, the distance to Chapter 3: Gravitational Lens Optics 31 Figure 3.3: A ray emitted at T/ on the source plane pierces the deflector plane at £. The ray is deflected by a(£) onto the observer. Chapter 3: Gravitational Lens Optics 32 an object at redshift z, based on its angular extent on the sky, is given by D ( z ) = 7 i—^202 (Qz ~ ( 2 - AKv'Jte + l - 1)) 3000/i_1Mpc , (3.7a) where i J 0 = 100/i Mpc/km/s is the Hubble constant (SEF92). Figure 3.4 shows that this measure of distance is linear at small z, but rolls over at larger redshifts. The distance becomes less and less sensitive to the redshift, weakening the redshift-prediction power of lens inversion. Because of the expansion of the universe between the time the ray is emitted, deflected, and observed, Dds =fi Ds — Dd- Instead, the angular size distance between a distant plane at redshift zs and the intermediate plane at redshift Zd, based on the appearance of the further plane from the nearer, is D(zd,zs) = ^ ± . _ . _ v> 3000/i x Mpc Q2 (l + z„)(l + z,) 2 (3.7b) 0.5 1.0 redshift z Figure 3.4: The angular size distance to objects at redshift z. Before we proceed with a detailed analysis of the lens equation (3.6), it is instructive to take a top-down look at the ground-breaking results of Refsdal (1964a; 1964b). A single point mass described by a surface density £ ( £ ) ' = MS2(£) greatly simplifies the lens equation. Because of the spherical symmetry, the problem is reduced to only one dimension with r) > 0: D9t 4GMDdsl A source at 17 will appear at location(s) £ on the deflector plane for all solutions £. The l / £ deflection suggests that rays passing close to the point mass M are too greatly bent, while Chapter 3: Gravitational Lens Optics 33 rays pass far from the deflector are not bent enough. This behaviour is described by the lens equation when we rewrite it as a quadratic in £: 2 Dd 4GMDdDds * = ° • The two solutions have opposite signs indicating the two images will appear on opposite sides of the line-of-sight through the point mass. The solutions merge when 77 = 0: perfect source-deflector-observer alignment. An Einstein Ring of solutions appear surrounding the deflector at the Einstein Radius 2 _ AGMDdDds While this perfect alignment is improbable, the Einstein Ring, or portions of it, is a clear signature of gravitational lensing and is very useful in the lens inversion process. It is convenient to scale the variables in the lens equation (3.6): x = £ / £ 0 i V = ^/ (s^Co)-In much of the results that follow, we set the scale length to *° = 360 • 60 • 60°d ' ( 3 - 8 ) so that the coordinates x and y are measured in arcseconds on the sky, the standard unit for astronomical data. The lens equation simplifies to y = x- a(x) (3.9) where the scaled deflection angle is defined by a(x) = ^^a((0x) . (3.10) In contrast to the geometric formulation of the lens equation, SEF92 shows how to apply Fermat's Principle, that light follows the shortest path, to derive the same relationship. Recall-ing that the deflection angle can be written as the gradient of a scaler field , light will appear at x for each solution (x,y) to V(/)(x,y) = 0, where (/> is the "Fermat potential" Chapter 3: Gravitational Lens Optics 34 The network of critical lines and corresponding caustics can be explored by examining the singularities in the potential (j>. The time delay between multiple images can also be written as a function of this potential. As we are more interested in the appearance of the specific lenses, and less on the characteristics of gravitational lenses in general, we do not follow this approach. 3.3 P r o p e r t i e s of the Lens E q u a t i o n The simplicity of the lens equation (3.9) belittles the amount of information, and complexity, it encodes. In its simplest interpretation, the lens equation shows that any signal emitted from y on the source plane will pierce the deflector plane and deflect onto the observer for each x satisfying the equation. Unfortunately, we observe just the opposite. Images appear at positions x for every pair of unknowns a(x) and y(x) satisfying (3.9). This is why direct inversion is impossible. We need not throw up our hands in despair, though. The lens equation is perfectly tailored for forward modeling; that is, numerically simulating lensing systems. Here we explore how this approach can be exploited. Forward modeling of gravitational lenses is possible because lenses show "optical reci-procity," the property that light travels through the lens the same in either direction. This is valid because of the thin lens approximation and because the wavelength of the light is much smaller than the thickness of the lens. Signals traveling through an inhomogeneous medium do not generally show this behaviour: one cannot "reconstruct" an earthquake based on the signals recorded on a seismograph, for interference and dispersion of the quake-induced shocks remove much of the information about the origin of the event. Rays of light travel in a predictable way going either direction through the lens, so that, for each x and a(x) we specify, a unique source coordinate y is produced by the lens equation. Any source of light that exists at y can be pulled back through the lens, and "observed" at x. Lens modeling is not quite as simple as this, however, as we observe beams of light rather than rays. Distortion and focusing of these beams lead to amplification of the signal. The signal recorded at the telescope, the charge in each pixel of a CCD chip for instance, is a measure of Chapter 3: Gravitational Lens Optics 35 the flux I that strikes that region. This signal is interpreted as having originated from a region of the sky of size (Ax) 2 with surface brightness pj = I/(Ax)2. The signal is traced back to the flux S leaving a region (Ay) 2 of the background source distribution of surface brightness ps — S/(Ay)2. Gravitational lensing preserves this surface brightness, pi = ps, so that the observed signal is related to the emitted signal by r _ c ( A * ) 2 S ( W * The lens equation does not provide y and x(y) needed to evaluate the right-hand side to produce the brightness I, however. Instead, we perform a little mathematical sleight-of-hand. The Jacobian A of the lens equation (3.9), with components A{j = dyi/dxj, i,j = 1,2, relates differential areas on the source plane to corresponding areas on the deflector plane. The ratio of these areas, the determinant |A|, is simply inverted to find the ratio of the differential areas of "pixels" on the deflector and source planes: 1=1- . The amplification factor p,(x) = |A|_1(a;) encodes the structure of the amplification of the lens. The magnitude of the amplification factor p, is the magnification of the intrinsic background signal. The sign of p, describes the parity of the image. Regions of the deflector plane with odd parity, where \A\ < 0, have been turned "inside-out" by the lens, much like the image seen on the inside of a spoon. In regions of even parity, |A| > 0, images are merely distorted, like the reflection on the back of the aforementioned spoon. These two types of regions are separated by the curves where |A| = 0, which we refer to as critical lines of the lens. The amplification diverges along these curves, extremely (infinitely) magnifying any lensed image which forms nearby. Projecting the network of critical lines back to the source plane via the lens equation generates a corresponding network of caustics. A highly magnified (and thus observable) image will form if there is a source of light lying near the caustic on the source plane. The name "caustic" has a slightly different meaning in the gravitational lensing community Chapter 3: Gravitational Lens Optics 36 than it does elsewhere because of the way the lens equation is used backwards to predict source coordinates y from deflector coordinates x. On a sunny day, a bright, "3"-shaped caustic of light is formed on the surface in a cup of coffee. Imagine shining light from, not onto, the surface of the coffee. A sheet of light emitted from the location of the caustic will reflect off the sides of the mug and converge on a point in space-the position of the sun, in fact. In the gravitational lenses described here, the observer sits at this point in space. The caustics, therefore, at not on the image plane, but on the source plane. In spherically symmetric lenses described only by a radial coordinate r = \x\, the lens equation (3.9) reduces to a 1-dimensional relation, y = r — a(r), independent of the position angle. It is easily shown that the Jacobian of the lens equation reduces to Wherever |A| = 0, the amplification factor p, = \A\ 1 diverges, marking the critical lines of the lens. We see, therefore, two conditions which generate critical lines and caustics, Along the critical line tracing the solutions of (3.12a), a = r with no constraint on the position angle. Images which form through this region of the lens are are tangential arcs with a narrow radial extent but wide angular spread. This line is the origin of Einstein Rings which can form in the case of perfect source-lens-observer alignment. The critical line defined by (3.12b) depends on radial variations in the deflection angle. Lensed images appearing near this line have a wider radial extent and are generally refered to as radial arcs. Because of the different types of lensed images formed by these two classes of solutions, we refer to them as the tangential and radial critical lines, and their pre-images, the tangential and radial caustics. (3.11) (3.12a) r (3.12b) 3.4 A T o o l b o x of Lenses While any mass distribution can be used to model the lens, we consider just a small collection of lenses with specific density profiles. The distributions described here are spherically symmetric, Chapter 3: Gravitational Lens Optics 37 so the behaviour of the elliptical versions used in the models will differ slightly. Throughout this section, the scaled coordinate r = £ /£ 0 is measured in arcseconds on the sky. The singular isothermal sphere (SIS) is to gravitational lensing what the simple harmonic oscillator is to mechanics: a simple yet well-motivated model which displays most, if not all, of the interesting behaviour of the system. As shown below, the SIS deflects passing rays of light through a constant angle. This is not merely a model chosen because of its simple analytic form, though. The SIS shows a constant line-of-sight velocity dispersion and models the equilibrium state of cluster of stars and galaxies. The SIS model describes very well the behaviour of large halos of gravitating matter. The simplicity of the model comes at a cost of a singularity in the mass density at the centre of the distribution. To remove this analytic anomaly, the core of the SIS is softened in the pseudo-isothermal (elliptical) mass distribution. We refer to this model as a PID lens avoiding the unwieldy PIEMD acronym of Kassiola & Kovner (1993). The price of removing the singularity of the SIS is the introduction of a model parameter, the core radius rc. The SIS and PID models describe both luminous matter and dark-matter. Recent Af-body simulations of dark-matter by Navarro, Frenk, & White (1995) produced a new density profile which describes the near-equilibrium configuration of cold dark-matter over a wide range of mass scales. We include this model in our toolbox because the major mass component of most strong gravitational lenses is attributed to a dark-matter halo surrounding the galaxy-cluster core. The shape of the NFW profile is characterised by a scale radius rs which is generally a significant fraction of the virial radius of the cluster. The last model we parameterize here was introduced by Miralda-Escude (1995) to model the mass of the luminous cD galaxy often found at the centres of galaxy clusters, under the mass-follows-light hypothesis. This "cD" distribution is characterized by a core radius rc and a larger cut-off radius rv Chapter 3: Gravitational Lens Optics 38 Six functions describe the lensing behaviour of each of these lenses: p - 3-dimensional mass density M - cumulative 3-dimensional mass S - surface mass density a - deflection angle p - amplification factor o~ios - line-of-sight velocity dispersion Because lens models often contain more than a single deflecting mass and because of the large numerical computations that must be made to produce a high resolution model, we strive to find an efficient computation scheme. The density, cumulative mass, surface mass density, and deflection angle of a composite lens are linear superpositions of the contributions of each lens mass. The magnification and line-of-sight velocity dispersion, however, are non-linear com-binations of the components. Thus while the first four of these quantities can be tabulated independently for each lens component, the final two can only be calculated once all the com-ponents are known. By making a judicious choice of the form of the normalisation of the independent functions, though, the computation needed to calculate the non-linear quantities can be greatly reduced. The mass of each lens is controlled by a parameter = 400. The lens plane is assumed to lie at redshift Zd = 0.250, so that 50'.'0 on the sky corresponds to roughly 1 8 0 / i _ I kpc. The source plane lies deep behind the lens at redshift zs = 0.750. The other model parameters used are shown in Table 3.1. In the functions that follow, the leading constants do not carry all the dimensions of the quantities leaving a dimensionless function to encode the shape of each function. Instead, the functions are split to ease the computations. Chapter 3: Gravitational Lens Optics 39 Parameter Values for Sample Lenses Lens mass Scale a aios[cD](2'.'0) (arcsec) without cD with cD (kms - 1) SIS - 999.20 984.24 970.34 PID r c = 4.0 1025.29 1011.97 602.89 NFW rs = 40.0 2441.70 2415.61 642.28 cD rc = 0.569 -rh = 1.428 - -Table 3.1: Parameter values for the lens masses. The observed value of aios[cD](2'.'Q) is explained in Section 3.4.6. 3.4.1 Mass Density The most fundamental description of any mass distribution is its 3-dimensional mass density, p. All other quantities are derived from the density. The density of a lensing mass with mass parameter a and any other (scale) parameters p can be written as p{r;a,p) = p0a2p(r;p) (3.13a) 1 2irGe0 Po = 7TIF^2 (3.13b) where Psisir) = ^ (3.14a) PpiD{r;rc) = J—j (3.14b) e l + fe) PNFw(r',r,) = — T2 (3.14c) PcD(r;rc,rh) = (3.14d) Here and'below, the barred forms like p are functions only of the scale parameters. Note that the density of the SIS diverges as r~2 at the origin, and that the SIS lens model has no intrinsic Chapter 3: Gravitational Lens Optics 40 scale, so that a SIS lens looks the same from all distances. The PID follows the SIS profile decaying as r~2 at large radii, but remains finite at the core. The NFW density diverges only as r _ 1 at the origin and falls off as r - 3 for r » rs. Figure 3.5 compares the densities of these lenses. arcseconds Figure 3.5: Three-dimensional mass density p(r). 3.4.2 Cumulative Mass Critical to the calculation of the line-of-sight velocity dispersion is the 3-dimensional cumulative mass function, simply the total mass within a given sphere: Af(r)= frp(r')Anr'2dr' . J o With the forms of the density p(r) given in (3.13), it is convenient to write the cumulative mass function in the form M(r;a,p) = M0a2M(p) (3.15a) M0 = ^ (3.15b) Chapter 3: Gravitational Lens Optics The different lens masses are described by MSrs(r) r (3.16a) MNFw(r;rs) MPID(r;rc) _r|_ 2r _ rc \ \ r c ) LJ (3.16b) (3.16c) McD(r;rc,rh) —i 1 tan — rc (3.16d) The profiles are shown in Figure 3.6. The cumulative mass of the SIS distribution grows linearly with radius, as does the PID lens. Due to the finite core density, though, the mass of the PID is slightly less than a comparable SIS. The diverging cumulative mass at large radii is clearly unphysical. Models of lens masses remove this problem by specifying an outer truncation radius. The form of the truncation does not generally effect the lensing properties of the lenses, however, as the deflection angle depends only on the mass within a given radius. 3.4.3 Surface Mass Density The 3-dimensional density and cumulative mass describe the local neighbourhood of the lensing mass and are naturally 3-dimensional quantities. The surface mass density and the deflection angle below describe the appearance of the lensing masses in the much larger source-lens-observer system. Through the thin lens approximation, these functions become 2-dimensional on the lens plane. The surface mass density £ is the mass per unit area on the lens plane, found by projecting the full 3-dimensional density onto the lens plane: where z is measured out of the lens plane. This is not a cumulative mass, but merely a mass density. Again, we find it advantageous to write the surface mass density in a form which simplifies computations at a later stage: S(r;cr,p) = £ 0 a 2 E ( r ; p ) (3.17a) Chapter 3: Gravitational Lens Optics 42 Chapter 3: Gravitational Lens Optics 4S E ° " 2cfe ( 3 1 7 b ) The barred functions below describe the shape of each lens' surface mass density. The compli-cated expressions for Y,NFW and the deflection angle (*NFW below are due to Bartelmann (1996). E 5/s(r) = - (3.18a) r Sp/ D (r;r c ) = - . 1 (3.18b) 2 f(r) ^NFw{r;rs) = ——2 , where f(r) = < 2 P ^ T = _ t a n - \ / ^ T I if i > l I V r 3 2 1 - -= = t a n h - 1 y i T f i f f i < 1 (3.18c) 0 if -E- = 1 Despite its daunting appearance, ^NFW{^~ < 1) and S A T F W ( ^ - > 1) evaluate to the same result at all r 7^ r s , the imaginary parts perfectly balancing. At r = rs, UNFW = 5^73- The surface mass densities of these lenses are shown in Figure 3.7. 3.4.4 Deflection Angle The paths that rays follow between the source and the observer are governed by the deflection from the vacuum geodesies due to the intermediate lensing mass. In the thin lens approxima-tion, all the angular deflection is assumed to occur abruptly on the lens plane. In spherically symmetric lenses, the deflection angle at radius r is toward the centre of mass of the lens with a magnitude depending on the total projected mass of the lens within the disk of radius r. Just as in the classical Newtonian theory of gravity, all contributions to the deflection due to mass outside the disk of radius r are exactly balanced, and all mass inside the disk can be treated as a point mass as the centre. Summing the contributions due to constant-density annuli gives a Chapter 3: Gravitational Lens Optics 44 1.8 arcseconds Figure 3.7: Two-dimensional surface mass density E(r) of the lens models found by projecting the mass onto the lens plane under the thin lens approximation. Chapter 3: Gravitational Lens Optics 45 deflection angle We write a(r) = c2 £, 0£>s r j o y ' a(r;a,p) = a0cr2a(r;p) Oin = ZoDsc2 with (*sis(r) ctpiD(r;rc) a N F w i ^ r s ) = 1 • m r 2 r, 1 if £ < 1 if f = 1 acD(r;r c , r f c ) = 3 + 2 f e ) 2 - f e ) ! + (3.19) (3.20a) (3.20b) (3.21a) (3.21b) (3.21c) (3.21d) At large radii where aipjn a l - ^ , the PID and SIS lenses appear almost the same. Figure 3.8 is a semi-logarithmic plot of the deflection angle of the model lenses emphasising the small-radius behaviour of the dark-matter and dark-matter+cD lenses. The deflection is measured in arcseconds, showing that the SIS in this test collection has an Einstein radius of ~ 17'.'0. The PID and NFW lenses are very similar in the core of the lens. The presence of the central cD galaxy only slightly increases the magnitude of the deflection at the centre of the lens. The cD's influence diminishes at larger radii where it contributes less and less to the cumulative projected mass density. Chapter 3: Gravitational Lens Optics 46 arcseconds F i g u r e 3.8: The deflection angle a(r) developed about model SIS, PID, and NFW lenses, with and without a central cD galaxy. Curves marked by o are the cor-responding dark-matter halos with the cD distribution added at the centre. The dark solid line follows a = r and locates the critical lines developed by the lens (Section 3.4.5). Chapter 3: Gravitational Lens Optics 3.4.5 Amplification Factor The background source of a gravitational lens would most likely not be visible in the absence of the lens. It is the extreme magnification that brings the distant source up to an observable level. The observations are therefore greatly biased, as often only those parts of the source which lie on or near the caustics can be observed. Therefore, understanding the nature of a gravitational lens is very closely tied to modeling the network of critical lines and caustics. We take a moment here to dispell a myth which has developed in the lensing community. Mellier et al. (1993) show that for potentials producing PID lenses, the radius of the radial critical line is proportional to the core radius of the distribution. This led to the folklore that without a core radius, no radial critical line will develop. This is true for the PID of course, but unfortunately this "conclusion" was extended to all mass distributions. As we show here, any distribution with mass density proportional to r _ 2 + ^ for 6 > 0 which produces a tangential critical line will also develop a radial critical line. If a lens shows a tangential critical line, there exists some radius r at which a(r) = r. Mass distributions with 3-dimensional mass density p(r) oc r~2+@, 6 > 0 produce no deflection at the lens centre, a(r) oc r@ —r 0 as r —> 0. As is obvious from Figure 3.9, and guaranteed by the Mean Value Theorem, there exists r' between 0 and r at which %p(r') — 1. This is precisely the condition describing the radial critical line. The choice of model parameters in these toy lenses produces tangential critical lines in the SIS, PID, and NFW lenses. We expect, therefore, radial critical lines in the PID lens, following Mellier et al. (1993), and also in the NFW lens, as PNFW OC r _ 1 at the centre of the lens (3.14c.) For lenses composed of a single mass, substitution of the deflection angle (3.19) into (3.11), and an appeal to the Fundamental Theorem of Calculus, yields (3.22) Chapter 3: Gravitational Lens Optics 48 Figure 3.9: The solid line traces the deflection angle a(r); the dashed line follows a — r. The Mean Value Theorem guarantees the existence of a point r' between 0 and r where ~{r') = l. The reason for our choice of function normalisation becomes clear, as we can write this as |A\ = (l - a0a2 ^ ) ( l + Oi0o-2 - s ) ) (3.23) In composite lenses where the deflection angle is a linear superposition of the dark-matter (DM) and luminous cD components, (3.22) expands to ( l + «o (°IM " *DM) + °\„ {^f - %D) ) ) (3.24) While the amplification factor is a non-linear combination of the deflection angles of the lens components, this expression removes the computational difficulty by decoupling the lens com-ponents. Not only can the functions a and £ be tabulated ahead of time, but the mass of the individual components a2 can be easily varied without extensive re-calculation. Figure 3.10 shows the magnification factor around each of the model lenses described here. Because of the huge range in [A]-1, the plot is logarithmically scaled, which excludes showing the sign of \A\ and the parity of the images that form. The parity changes across each divergence in p, with even parity regions at small and large radii, and and odd parity band between. Chapter 3: Gravitational Lens Optics 49 _4i , , ,— • • •— 10° 10 1 10 2 arcseconds Figure 3.10: The magnification factor |/x|(r) around each lens in this collection. Curves marked by o are the corresponding dark-matter halos with the cD distribution added at the centre. The SIS has only one tangential critical line at roughly 17'.'0. The PID and NFW lenses have both outer (tangential) and inner (radial) critical lines. Chapter 3: Gravitational Lens Optics 50 In the SIS lens where a is constant, (3.12b) has no solution and the only critical line lies at the Einstein Radius of the lens: The Einstein Radius is roughly 17'.'0 with the model parameters here, which manifests itself as the constant asis ~ 17 in Figure 3.8 and the divergence in the magnification at ~ 17'.'0 in Figure 3.10. The PID and NFW lenses show very similar critical lines, both an outer tangential line and an inner radial line. Notice that the tangential critical lines appear at radii where apjD and OLNFW cross the a = r curve in Figure 3.8. The radial critical lines appear at radii where apiD and CXNFW are parallel to the a = r curve in Figure 3.8. The presence of the cD galaxy in a composite DM+cD lens does not greatly change the magnification but it does slightly alter the position of critical lines. Because of the extreme bias of observing only strongly lensed light, the cD galaxy can play a significant role in appearance of the lensed images. 3.4.6 Line-of-Sight Velocity Dispersion There are numerous methods employed to estimate the amount of matter in a galaxy by mea-suring the galaxy's influence on it environment. It is difficult to measure the mass of galaxies which orbit each other because of the unknown projection angle of the orbit and the extremely long period of the orbit. Weak gravitational lensing analyses of galaxy-galaxy lensing is be-coming a well-defined method for measuring galaxy masses (Dell'Antonio and Tyson, 1996). The mass of the central cD galaxy in a cluster can also be estimated by its effect on the strong lensing behaviour of the lens (Natarajan and Kneib, 1997; Miralda-Escude, 1995). Since we seek to extract the mass from models of strong lensing, we need a different description of the mass to normalize our results. A common measure of the mass of a galaxy is the line-of-sight velocity dispersion of the stars within the galaxy. This account follows that of Binney and Tremaine (1987). A gas with mass density profile p(r), be it a cloud of molecules, a galaxy of stars, or a cluster of galaxies, develops a line-of-sight velocity dispersion aios depending on its "temperature." In the cases of stars and galaxies, fluctuations in the line-of-sight velocity of the individual members Chapter 3: Gravitational Lens Optics 51 of the distribution lead to a broadening of the spectral lines. The width of the broadening is proportional to the dispersion in radial velocities, which we can interpret as temperature through the following analogy. An isothermal gas of particles with mass m reaches hydrostatic equilibrium when d ( 9 d l n p \ AnGrr, -r"p / 2 lnp\ _ AwGrn 2 dr \ r dr ) ~ kBT Similarly, an isothermal stellar or galactic population with constant speed a reaches self-gravitating equilibrium when d f2d\np\ = torG^ dr \ dr J a2 Comparing these two expressions shows that the astronomical system can be treated as an ideal gas with "temperature" defined through 2 ksT a = m The line-of-sight velocity dispersion of a gas can be determined from the particle number density v through the relations /OO r > , v7?fT=iiT' (3-25) We assume (i) the velocity of the gas is isotropic, 8 = 0, and (ii) the spatial number density v is proportional to the mass density p. The Jeans' Equation balancing kinetic and potential energy states d , 2 x GM(r) (pa2) = - ^ ± p (3.26) where M(r) is the 3-dimensional cumulative mass with radius r. Under the assumptions, the line-of-sight velocity relations (3.25) reduce to roo — r 2 Gl Mfr'WrV „ dr' /oo M(r')p(r')-^ r ) = —^S5 • (3-27) Chapter 3: Gravitational Lens Optics 52 For a system with more than one species of "particle," like the cD galaxy immersed in dark-matter (DM) halo, the Jeans' Equation (3.26) represents a "partial pressure" relationship, where the component with density p contributes only a portion of the total gravitational potential. The cumulative mass M which enters (3.27) is the cumulative sum of both (or all) the masses in the system. We are particularly interested in lenses with more than one mass component, one of lu-minous baryonic matter and one of dark-matter. By measuring the velocity dispersion of the stars in the luminous galaxy, we can estimate its mass. There are no luminous dark-matter "particles" whose line-of-sight velocity can be observed, however. Instead, the velocity disper-sion of the cluster galaxies swirling in the dark-matter halo are observed. It is assumed that these galaxies are test particles in the dark-matter, and their line-of-sight velocity dispersion is a measure of the temperature, or mass, or the dark-matter itself. Thus we independently find the mass of the cD (using stars within the cD itself) and the mass of the dark-matter halo (using cluster galaxies as test particles). Each lens in this collection is described by a mass parameter a related to the observed Gios. Only in the case of an isolated SIS is aios — a. The SIS mass distribution shows a constant, "isothermal," line-of-sight velocity dispersion. In other mass distributions, o~ios varies with radius and scales with a. In the models, we add to the dark-matter halo a cD galaxy with parameter 1 Figure 4.1: The visible array is a linear combination of the arrays generated by each deflector component. The registration of the field-of-view on each component's layer determines that component's position in the simulation. This "initialization and then look-up" routine can quickly display the lensed image once the arrays have been calculated and stored in memory. The goal of lens, though, is to allow for movement of the deflector and source components. Each time a component is shifted, a new image must be calculated and displayed. Re-calculating each array takes far too long. Instead, lens initializes the calculations on grids much larger than are immediately needed. We specify Chapter \: Numerical Models of Gravitational Lensing 61 the half-width W and half-height H for the displayed array. Half-dimensions are used so that the "visible" field-of-view is (2W +1) x (2H +1) and has central pixel. Any object, a deflector or source component, sits in pixel (i,j)v, where the subscript V indicates the coordinates in the visible field. As the window has width 2W + 1, and only one pixel is occupied by the centre of the object, there are 2W alternate values for the i-coordinate of the object. Similarly, there are 2H other possible values for the j-coordinate of the object. To allow for a full range of movement of the object, the calculated arrays must be large enough to contain the the immediately visible field and an additional surrounding buffer of width 2W and height 2H. This underlying, or "global" array has dimension (AW + 1) x (AH + 1). Every object has two sets of coordinates: visible coordinates (i,j)v and global coordinates (i,j)G- The two sets are related by a "registration" factor (ia, J0)R set by the initial geometry of the lens: (i,J)G = (i,J)v + (io,3o)R • (4-1) Figure 4.2 shows an example of the configurations of the global and visible arrays. The components of the deflection angle and terms that make up the amplification are calculated over the entire global array for each mass when the program lens is started. The registration coordinate for each mass is found. Then, for each configuration of objects chosen by the user, only the visible portion of each mass' global grid is used to calculate the deflection and amplification, as indicated in Figure 4.1. Masses which are included in the lens configuration but not visible contribute to a constant background array upon which the visible arrays are built. The amplification factor is a function of position, the deflection angle, and the gradient of the deflection angle. Analytic expressions for the derivatives of the deflection about elliptical lens models are not easily developed, especially in closed form. In lens, the deflection angle is calculated with numerical integration, and the derivative is numerically approximated. The distance A between the centres of the grid cells is set by the user, and typically ranges between 0''1 and 0''5 (for instance, 0'.'207 at CHFT.) To reduce the truncation error, a 5-point centred-Chapter 4'- Numerical Models of Gravitational Lensing 62 global array 17 16 15 14 13 12 ^ 1 1 lio' 't 9 8 g"1 IS 7 £ 6 bjO -5 4 3 2 l " rt o o o .3 i visible field-of-view Object 1 2 3 4 5 6 7 8 9 10 11 i visible coordinate i i 1 ' 2 ' 3 ' 4 ' 5 ' 6 ' 7 ' 8 ' 9 , 1 0 , l l , 1 2 , 1 3 , 1 4 , 1 5 , 1 6 ' l 7 , 1 8 ' l 9 , 2 0 ' 2 1 global coordinate i Figure 4.2: A simple example with half-width W = 5 and half-height i / = 4. The offset (i0,j0)R — (8,6) registers the visible field on the underlying global grid. This registration relates the visible coordinates (3,3)y to the global coordinates (11, 9 ) G of the Object via Equation (4.1). The global grid is just large enough that the visible "window" can be repositioned with the Object appearing at any location ( l , l )y thru ( l l , 9 ) v . differencing scheme with fourth order truncation is used: 8a{i, j) a(i - 2,;') - 8a(i + 8a(i + 1, j) + a(i + 2, j) + 0(A 4 ) dxi 12A To avoid boundary problems, the underlying global grid is surrounded by two more layers of cells used only to calculate the derivatives of the deflection near the edges of the view window. Because of the mass-layers approach to modeling the lens, the computation time to update the image after altering one component's parameters is independent of the number of lens components. This allows for easy manipulation of quite complicated lens models. Of course, the initialization time increases in proportion to the dimensions of the field-of-view in the Chapter 4: Numerical Models of Gravitational Lensing 63 simulations, but this is a one-time calculation. The lens code itself consists of slightly more than 5000 lines of code written in Fortran, with calls to Tim Pearson's graphics library, PGPLOT (version 5.0). This library of routines is "a Fortran-callable, device-independent graphics package for making simple scientific graphs." (Pearson, 1997). A number of system calls are made to the U77 extensions to the standard Fortran library. To implement comparison with the data, lens makes calls to the f i t s io library of FITS data manipulation routines. This Flexible Image Transport System is widely used by the astronomical community. 4.2 Model Parameters Because we restrict our attention to the family of lenses described in Section 3.4, it takes only a small number of parameters to build a lens model. These parameters come in three sets: geometry, deflector, and source. All of the parameters are initialized in a text file (a Fortran namelist) passed to lens on the command line. Those parameters marked by * below can be altered interactively in the current version of lens. 4.2.1 Lens Geometry Before specifying the models for the deflector or source, we establish the geometry of the lens through the parameters listed in Table 4.1. Redshifts zs and z& set the angular size distances (3.7a) to the source and deflector planes, respectively. The size of the numerical computation is set by the half-width W and half-height H of the simulated image, which are limited only by the amount of available memory, and the speed of the computer. Each pixels covers fl = A 2 square arcseconds on the sky. 4.2.2 Deflector Mass Distribution Any number of masses represented by the 2-dimensional parametric models described in Sec-tion 3.4 can be added to the lens. Each mass is described by its centre XM, eccentricity e, and orientation (or position angle) of the major axis M, measured counter-clockwise in de-Chapter 4-' Numerical Models of Gravitational Lensing 64 Lens Geometry Parameters deflector plane redshift* zd source plane redshift* zs half-width of image plane (pixels) W half-height of image plane (pixels) H pixel spacing (arcsec) A centre field-of-view (arcsec) x0 Table 4.1: Parameters determining the geometry and resolution of the simulation. grees from the positive rc-direction. Mass distributions can be truncated at radius R, but this is generally ignored for the strong lensing simulated here: Dark-matter halos extend well be-yond the core of the galaxy-clusters we focus on, and, as discussed in Section 3.1.3, mass lying beyond an outer homoeoidal shell does not contribute to the deflection within. The mass of each lens component is set by the parameter a which can be set via (3.28) or (3.29) to mimic an observed line-of-sight velocity dispersion. Choosing between the SIS, PID, NFW, and cD profiles, together with their model-specific parameters, completes the description of the lensing mass distribution. Table 4.2 list the complete set of deflector parameters. Deflector Component Parameters centre* (arcsec) XM eccentricity e orientation (degrees) truncation radius (arcsec) R mass parameter* a scale parameters (arcsec) SIS PID NFW cD rc rs rc rh Table 4.2: Each mass distribution in the lens is described by a small set of model parameters. Chapter 4'- Numerical Models of Gravitational Lensing 65 4.2.3 Source Light Distribution Solving the lens equation through the centres of each of the (2W + 1) x (2H + 1) image plane pixels generates the same number of corresponding source plane coordinates. These points are the origin of any light which will appear on the image plane, modified by the appropriate amplification. The source coordinates y are not restricted to integer multiples of the scale length A as are the image plane (deflector plane) coordinates x. The source light distribution is specified by an analytic function S(y) which is evaluated at each y to find the flux emitted from this point on the source plane. In Chapter 6, we pull the data in the observations back to the source plane, so it is natural to measure the strength of the source in ADU. Various parameterized distributions are available to model different types of background sources. Again, those parameters marked * can be altered interactively within lens. The uniform circular disk, shown in Table 4.3, is the simplest background source. Any asymmetry in the shape or intensity of the lensed image can be traced to the effects of lensing, rather than a combination of lensing and intrinsic source structure. Uniform Circular Disk centre* (arcsec) flux (ADU) radius (arcsec) Table 4.3: Uniform circular disk source distribution of light. In the example shown here, ys = (0,0), S = 5, r = 0.5. The pseudo-isothermal disk shown in Table 4.4 better approximates the appearance of a background galaxy. This light distribution follows the spherical PID mass distribution, giving the circular disk a more realistic light profile. Chapter 4-' Numerical Models of Gravitational Lensing 66 Pseudo-Isothermal Circular Disk centre* (arcsec) peak flux (ADU) core radius (arcsec) Vs S Table 4.4: A pseudo-isothermal profile mimics a background elliptical galaxy. This example is described by ys = (0,0), S = 5, r = 0.5. The core radius is indicated by the circle of radius 0'.'5. In the lens inversion scheme described in Chapter 5, first the position and shape of the lensed features are established using a uniform elliptical disk for the source. The parameters characterising this source model are shown in Table 4.5. This model was added late in the development of the software and makes redundant the uniform spherical disk. For illustrative purposes, especially to show how a small amount of source structure can produce incredibly complex image structure, a two-armed spiral galaxy model shown in Ta-ble 4.6 is included in the list of source distributions. The profile was adapted from an ingenious model of Bonnet (1995): SM = ^ e - 2 S i n 2 ( ^ ° - T r 2 ) , where r and 6 are polar coordinates of the point y measured with respect to the centre of the distribution at ys. The exponential function favours those points lying near the polar curve r r 2 = 9 — 0o, then fading as a gaussian from the peak. The taux d'enrollement (spiraling rate) r adjusts how tightly the arms spiral about the central bulge. As described below in Section 5.4, the reconstruction of the source distribution from data on the image plane is subject to over- and under-sampling because of the focusing of the lens. To explore the effects and arrive at a strategy, a test pattern in the author's initials "PN"can Chapter 4-' Numerical Models of Gravitational Lensing 67 Uniform Elliptical Disk centre* (arcsec) flux (ADU) major-axis* (arcsec) eccentricity* orientation* (degrees) Vs S a e $ Table 4.5: The shape and orientation of the uniform elliptical disk can be altered to mimic the shape of any background elliptical galaxy. The model shown here is desribed by ys = (0,0), r = 0.5, S = 5, e = 0.7, = 30. Two-Armed Spiral Galaxy centre* (arcsec) ys peak flux (ADU) S0 core radius (arcsec) rc orientation at ys (degrees) d0 spiraling rate r Table 4.6: The two-armed spiral models a background spiral galaxy with a central bulge. The spiral distribution in this Table is modeled with by ys = (0,0), S — 100, rc = 0.5, d0 = 30, r = 1. Chapter 4'- Numerical Models of Gravitational Lensing 68 be viewed through the lens. The pattern, parametrized in Table 4.7, is simply a bitmap on the source plane, except the "P" has twice the intensity of the "N". The pattern takes on a ragged appearance when it is given an orientation and then mapped onto the rectangular array in the image. 'PN" Test Pattern base of "P"* (arcsec) peak flux (ADU) orientation (degrees) Vs S 4> Table 4.7: A test pattern (the author's initials). The simulation shown here is described by parameters ys = (-0.8, -0.8), S = 100, = 15. 4.3 Degrees of Freedom in a Parametric Lens Model Two parameters set the redshifts of the deflector and source planes. Each deflector component has a center (two parameters), orientation, eccentricity, mass, and generally one scale length, for a total of 6 parameters. Following the lens inversion scheme presented in Chapter 5, the background source is reconstructed from the data in the observations. There are no source model parameters in the final model, aside from the source redshift specified elsewhere. To first establish values for the deflector parameters and the redshifts, however, we use the uniform elliptical source. Five parameters, the centre, orientation, eccentrity, and major-axis, charac-terize this source. A single deflector, single source model therefore requires 2 + 6 + 5 = 13 parameters. Each additional mass adds 6 degrees of freedom; each additional elliptical source adds 5. To constrain the lens models following our parametrization scheme, at least 13 statistics must be extracted from the observations. Chapter 4' Numerical Models of Gravitational Lensing 69 4.4 O t h e r F e a t u r e s o f t he M o d e l i n g So f twa re Great effort was taken to make the gravitational lens modeling program lens simple and intuitive to use. The graphical output is a plot on the computer monitor. Each lens component can be moved one pixel (A arcseconds) at a time in any direction by clicking with the pointer in the desired direction. Simulations which require finer adjustment of the deflector components are simply run after specifying a smaller value of A in the initial configuration. The position of any background source is altered by clicking the pointer on the display in the desired direction. Because the source distribution is an analytic function, its centre is not constrained to move in units steps of A. By default, the increments on the source plane are set to A s = (Dd/DS)A, scaled so that A and A s represent the same physical length on the deflector and source planes, respectively. The user can repeatedly double or half step size once the program is running. In addition to the pointer-driven actions which change the positions of the deflector and source distributions, single keystrokes and subsequent keyboard input allow the user to: 1. add symbols or text to the current plot, 2. print the coordinates under the pointer to STDOUT, 3. write the current plot in gif format, primarily for producing images for HTML documents, 4. list the current model parameters of any lens or source component to STDOUT, as shown in Table 4.8, 5. toggle ± 1 0 x and ± 5 x amplification contours, as shown in Figure 4.3, 6. increase or decrease the mass parameter a of any lens component, 7. concatenate a sequence of gif s (using gif merge) to animate a sequence of moves, 8. write the current plot to black-and-white or colour PostScript, 9. increase or decrease the redshift zs of the source plane or the redshift zd of the lens plane, 10. cycle through linear, logarithmic, and square root scaling of the image intensity, Chapter 4-' Numerical Models of Gravitational Lensing 70 11. write the intensity of the pixel under the pointer to STDOUT 12. save the current lens parameters or a user-specified zoom window to file for subsequent runs of lens, 13. alter the grayscale of the image by changing the zero-point and range of the scale, analo-gous to the grayscale manipulation available with SA Olmage, 14. print a list of actions and other help to STDOUT. -1,0 - 5 0 arcseconds 5 10 Figure 4.3: The heavy dashed lines mark the critical lines of the lens where the ampli-fication diverges. Thin solid lines trace 10 x, 5x amplication contours (even parity). Thin dashed lines trace contours of —10x, — 5x amplication (odd parity). Information about deflector mass Ml: Redshift zd= 0.257 Centre ( .00, .00) arcseconds Orientation 30.0 degrees E c c e n t r i c i t y .700 LOS parameter 975.0 km/s Truncation radius 40.0 arcseconds PID — psuedo-isothermal sphere Information about source SI: Redshift zs= 0.825 Centre ( .78, 1.33) arcseconds DSK — uniform disk Radius rc= .40 arcseconds Intensity 5.0 Core radius rc= 4.0 arcseconds Table 4.8: Information describing each lens component (left) or source component (right) is available while lens is running. Chapter 4-' Numerical Models of Gravitational Lensing 71 Several actions tailored to the inversion scheme described in Chapter 5 are included. To find the "shape" of the lensed image in the plot, the polar moments (described in Section 5.3) of any lensed feature can be tabulated and overlayed onto the image, as shown in Figure 4.4. To faithfully compare simulations with the observations, the simulated lensed images can be degraded to mimic the smearing of the signal due to the Earth's atmosphere. The images in Figure 4.3 are convolved against a Moffat point spread function = ( 1 + p ^ 1 ) ( i y < « ) in Figure 4.5. The scale R in the Moffat function is the radius at which the function drops to one-half its central value. The exponent 0 controls the shape of this psf. Both parameters can be estimated matching the function to the profile of stars (point sources) in the observations. Demanding unit volume under the psf sets the normalization constant M0 and conserves the total flux throughout the convolution. Background noise can be added to the simulated images to reproduce the background sky in the observations and the noise introduced at the instrument, Figure 4.6. While data collection is a Poisson process as electrical charge accumulates on the CCD, there are sufficiently large numbers of events to approximate the distribution of noise with a gaussian profile. The mean and standard deviation of an appropriate gaussian can be estimated by looking at blank regions in the observations. One last feature of the lens software is the ability to reconstruct a surface brightness distri-bution of the background source which best reproduces the lensed features based on the lensed pixels in the observations, either real or simulated. Contrary to intuition, it is often difficult to fully recreate even the analytic source models, like the uniform disk or ellipse, because of the non-uniform way in which the lens equation samples the source plane. Details of the source reconstruction scheme are discussed in Section 5.4. Chapter 4- Numerical Models of Gravitational Lensing 72 arcseconds F i g u r e 4.4: Circles, rays, and annular sectors, built from polar moments (Section 5.3) characterise the position and shape of lensed arcs. F i g u r e 4.5: Lensed arcs in Figure 4.3 convolved with a Moffat point-spread function (4.2) with a 2-pixel FWHM. Chapter 4- Numerical Models of Gravitational Lensing IS arcseconds Figure 4.6: Gaussian noise mimicking the noise in the observations is added to the simulation. Chapter 5 A New Lens Inversion Scheme We have set up a formalism for describing gravitational lenses and numerical modeling, and described a software tool to implement this simulation. In this Chapter, we introduce a new lens inversion scheme. The result of this algorithm is a parametric deflector mass distribution and a background source surface brightness distribution, which together reproduce the appearance of observed gravitationally lensed features. After commenting on the need for yet another approach to lens inversion, we describe in detail the scheme proposed. Later, in Chapter 6, the algorithm is implemented on the lensed features observed in the galaxy-clusters MS 2137—2353 and MS 1455+22, members of the EMSS Einstein Medium Sensitivity Survey catalogue (Gioia and Luppino, 1994). 5.1 Why do we need another lens inversion scheme? There is no question that the lens inversion schemes described in Chapter 2 are successful, particularly the outstanding results produced by Tyson et al. (Lucent Technologies, 1997). Why, therefore, are we introducing another approach to lens inversion? A great deal of the success of the schemes described earlier relies on the high resolution HST observations of multiple multiply-lensed background sources, in Abell 2218 (Kneib et al., 1996), in CL 0024+1654 (DelPAntonio and Tyson, 1996; Kochanski et al., 1996; Lucent Technologies, 1997), and in MS 2137 (Hammer et al., 1996). Unfortunately, data of this quality is not always widely available, so a routine that can be implemented on current, ground-based observations, as well as the high-resolution space-based images, is useful. 74 Chapter 5: A New Lens Inversion Scheme 75 It is difficult to treat gravitational lensing as an inverse problem of light traveling through an inhomogeneous medium because the parameters that are sought belong to two, quite distinct, classes. These parameters describe the two objects produced by a successfully inverted lens: the deflecting mass distribution and the surface brightness distribution of the background source. The former is described by a collection of physically motived parameters, like the mass and core radius of each deflector component. Only a handful of these model parameters are needed to completely describe each component. The latter, when the source is reconstructed from observational data, is a pixelised rendering of the surface brightness of the background source. The parameters are an array of values of the surface brightness of each of the source pixels. Based only on the appearance of the lensed images, we cannot easily separate these two families of parameters. Without such a distinction, the surface brightness in any particular source pixel constrains the solution as much as, for instance, the mass parameter describing the dark-matter halo at the centre of the galaxy-cluster. The over-importance of the source pixels or the under-importance of the halo mass in such a scheme does not seem appropriate. A key stage in many of the gravitational lens inversion schemes currently in use is the identification of common structures in two or more lensed images of a single background source. The presence of such structures, often bright "knots" in the images, greatly constrains the geometry of the lens and the deflector mass distribution, for the common knots must trace back to a single source under the lens equation. For instance, the bright central bulge in a background spiral galaxy should give rise to brighter regions in each lensed image that forms. Several of the inversion schemes trace each bright knot back to the source plane, and use the distance between these pre-images as a measure of the goodness-of-fit. The LensClean algorithm (Kochanek and Narayan, 1992) traces the brightest pixel of each lensed image back to the source plane. Hammer et al. (1996) identify larger structures, as do Seitz et al. (1997) in their recent analysis of the lensed features in MS 1512+36. Each of these implementations is prone to the mis-identification of common knots, however. Because of the extreme magnification of images lying near critical lines of the lens, even a very low surface brightness feature in the source which lies near a caustic of the lens may be magnified to appear as a bright knot in the lensed Chapter 5: A New Lens Inversion Scheme 76 image. Adjusting the model parameters to make this region of the source coincide with the bright central bulge of the background galaxy may lead to incorrect parameters and models. Without very high resolution observations to conclusively identify multiply-imaged sources, or multiple features within each multiply-imaged source, these two interpretations of bright knots are indistinguishable. Both of these problems stem from the same source: the appearance of each observed lensed image is the product of two effects, the magnification of the background source due to lensing and the structure of the surface brightness of the source itself. The model parameters describing objects on the deflector plane and on the source plane are mixed. The origin of a bright structure in a lensed image may be due to magnification, or a bright structure in the background source, or even a combination of the two. 5.2 A Two-Stage Inversion Algorithm The key to the inversion scheme introduced here is to decouple the effects of lens magnification and source structure on the observations of a multiply-imaged background source. This is achieved by first deducing the geometry of the lens and the mass distribution of the deflector components by matching the positions of the lensed images in the observations. While the generic problem of distinguishing lensed images from other diffuse objects in the observations relies on the judgement of the observer, the characterization of images by their positions does not depend on variations in brightness, and the problem of mis-identifying knots of light does not arise. Once these "conduits" of light are established, the magnification of the signals passing through them is completely determined, for the magnification is a by-product of the lens geometry. The magnification can then be removed from the observations, revealing the bare, though still distorted, surface brightness of the source. We state the two-stage algorithm in Figure 5.1. This algorithm produces a parametric model for the lens (Step 1) and a reconstructed surface brightness distribution of the source (Step 2). By first establishing the geometry of the lens Chapter 5: A New Lens Inversion Scheme 77 Two-Stage Lens Inversion Algorithm Step 1. Build up a model of the lensing mass distribution and the lens geometry by matching the positions of each lensed image of a common source. This establishes "conduits" through which the solutions of the lens equation pass, and also the magnification these beams undergo. Step 2. Trace the observed data along the solutions to the lens equation back to the source plane, de-magnifying it along the way, and combine the multiple beams in an optimal way. Figure 5.1: The two-stage lens inversion algorithm we propose in this thesis. and deflector mass distribution, the effects of magnification and source structure are decoupled. The two classes of model parameters which are mixed in a one-shot inversion remain separate here, the model parameters produced in Step 1 and the pixelised source distribution in Step 2. 5.3 Specifying Lensed Image Positions with Polar Moments To fit the geometry of the lens and the model parameters describing the deflecting mass dis-tribution, it is necessary to match the position of the observed lensed features. Pixel-by-pixel comparisons between the data and the simulation is less than favourable for several reasons. Errors in the data, for instance hot pixels caused by cosmic-ray impacts on the CCD, could be incorrectly interpreted as lensed pixels of the background source. Ground-based observa-tions are subject to the effects of the seeing at the instrument at the time the observations are collected. The lensed features may be broadened, changing the number and position of pixels containing lensed light. Because of this smearing of the images, the information contained in neighbouring pixels is not independent. The smallest region of the image over which measure-ments are statistically independent, the resolution element, is generally much larger than a single pixel. The number of independent data that can be extracted from a single image is less than total the number of pixels in the array. Matching the observations and simulations pixel-by-pixel predominantly compares the noise in the observations with the numerical detail Chapter 5: A New Lens Inversion Scheme 78 of the simulations. The simulations can be convolved with a suitable point-spread function, but simulations which first match the number of pixels in the instrument may quickly exceed the performance of most desktop workstations. The interactive behaviour of the program lens may be lost. A more general description of the lensed features' positions is needed. We denote by "position" the centroid of the object, as well as any shape parameters that can be derived from the data. As discussed in Section 4.3, parametric models in which a single background source is lensed through a single deflector component require 13 or more parameters, depending on the structure of the source model. The (flux-weighted) centroid of each lensed arc in the observations generates only two constraints per arc. At least 7 images of the background source are required to constrain the model, and the centroids of a single triply-lensed background source are insufficient to fully contrain lens models following our parametrisation scheme. In weak lensing analyses, the second moments of the images are also tabulated. When calculated with respect to the centroid of an image, the ratio of the eigenvalues of the diagonalised matrix of quadrupole moments can be interpreted as the ratio of the major- and minor-axes of an equivalent ellipse lying in the diagonalised coordinate system. This describes quite well the shapes of weakly lensed images which are only small distortions of elliptical background galaxies. Strongly lensed arcs, like the tangential arc A3 and radial arc Al in Figure 1.2, are not elliptical in shape, though, and the constructed equivalent ellipse is not a faithful reproduction of the position of the lensed features. The quadrupole moment approach of weak lensing is tantalizing, however, because of the small number of parameters needed to the describe the position of the arc: a centroid (x~i,x~2) and three quadrupole moments Qn, Qn, Q22- The tilde'd form indicates these moments are measured with respect ,to the centroid of the image. Adapting this idea to strong lensing, we introduce "polar moments" of the lensed images. The polar coordinates (9, r) of the image are interpreted as if they are coordinates in a Cartesian coordinate plane. Figure 5.2 shows the arcs in the simulation in Figure 1.2 in the Cartesian coordinate system. Chapter 5: A New Lens Inversion Scheme 79 T 3 o o u cS 14-r 12-10-8-6-4-2-0-. -180 tangential arc, A3 radial arc, Al tangential arclet, A2 0 -120 120 position angle (degrees) 180 Figure 5.2: Polar r- and -^coordinates of the three gravitationally lensed images in the simulation shown in Figure 1.2 are plotted in a Cartesian coordinate system. The rectangular boxes about each arc are constructed from the quadrupole moments defined in (5.1d)-(5.1f). Chapter 5: A New Lens Inversion Scheme 80 In both the observations and the simulations, the lensed features are represented by a discrete set of values, rather than a continuous distribution of light. We tabulate the moments of the lensed arcs by summing over the individual lensed pixels. In our simulations where the background source is a uniform elliptical disk with a distinct edge, every non-zero pixel in the simulation contains lensed light. Summing over the pixels with intensity I > 0 includes only the lensed pixels. To find the polar moments of the observed lensed arcs in MS 2137 and MS 1455, however, we must introduce a threshold to separate the lensed pixels from the background sky. We write I > I0 to represent the lensed pixels in the data, with IQ = 0 in the simulated images. This threshold is discussed further in Section 5.4. Each of the N pixels in the data containing lensed light has polar coordinates (Oi, rj). The polar moments (5.lb)-(5.If) are defined following the pattern of the Cartesian moments: Q0 = N (number of lensed pixels) (5.1a) r = Qol E n (5-lb) Ii>Io 9 = Qo1 E °i (5-lc) U>Io Qrr = Qo' E f o - ^ ) 2 (5-ld) Ii>Io Qre = Qo1 E in-rM-d) (5.1e) Ii>Io Qee = Qo1 E ^ - ^ ) 2 (5-lf) These moments are not flux-weighted, but depend only on the positions of the lensed features on the image plane. The Oth moment Q0 is the area of the image, in units of fl — A 2 arcsec2. The radial moment f specifies the average radius of lensed image, for there is just as much weight in the image outside the circle of radius f as there is inside this circle. In the case of giant arcs, f should closely approximate the Einstein Radius of the lens. The moment 6 specifies the average position angle of the lensed image: there is just as much weight clockwise from the line at position angle 0 as there is counter-clockwise from this line. As 9 simply measures position Chapter 5: A New Lens Inversion Scheme 81 angle on the sky, the magnitudes of f and 9 are incomparable. The direction 9 = 0 is arbitrary, and care must be taken to use a continuous domain of angles 6 when measuring 9: the average position angle of a pixel at 1° and a pixel at 359° is 0°, not 180°. The second polar moments of an arc characterize its shape: Qrr is a measure of the radial spread, Qgg the angular spread. Analogous to the equivalent ellipse of weak lensing analyses, we construct a representative region based on the values of these components. A uniform rod of length 2L lying along the rc-axis between — L and L has a second moment Qxx = \L2. The length can be recovered from the moment, as L = *J3QXX. We apply this result to the polar moments of the lensed arcs. The region lying between f ± ^/3Q r r and 9 ± \f3Q$g is a rectangle in the Cartesian system, and what we refer to as an "annular sector" in the original polar coordinate system. The rectangular boxes representing the shapes of the the simulated arcs in Figure 1.2 are drawn in Figure 5.2. The annular sectors built from the quadrupole moments are shown in Figure 4.4, where the annular sectors sit at the intersection of a circle of radius f and a radial line at position angle 9. The mathematical description of the annular sector allows for the possibility f — \JzQrr < 0. This gives the annular sector a bow-tie-like shape across the origin, and is surely not a faithful representation of the shape of the lensed arc which generates the moments. A similar patholog-ical case exists in the equivalent ellipse approach used in weak lensing. Formally, the matrix of second (Cartesian) moments can have positive and negative eigenvalues. We would be forced to interpret a negative eigenvalue as a negative semi-major or semi-minor axis. This problem does not arise, however, because the shapes of the weakly lensed background galaxies do not generate such pathological sets of quadrupole moments. In none of the countless simulations of strong gravitational lensing carried out, nor in any of the arcs analysed in MS 2137 and MS 1455, did we find a single instance where f — \J%QTT < 0. We remain aware of the possibility, but we have yet to find or simulate an arc which produces such a result. The polar moments of the three simulated arcs in Figure 1.2 are shown in Table 5.1. The pixel spacing in the simulation is A = 0'/2 so the arcsecond measurements of r j are accurate Chapter 5: A New Lens Inversion Scheme 82 to ~ O'.'l. An error of A/2 arcseconds at a radius of r corresponds to a uncertainty in the position angle SO = degrees . (5.2) nr The uncertainties Sr — A/2 and SO given by (5.2) are included in the Tables of polar moments throughout this thesis. While the absolute error SO is important to note, the relative error SO/0 is meaningless for this angular coordinate. Moment A l A2 A3 arc type radial arc tangential arclet tangential arc Qo (pixels) 95 58 194 r (arcsec) 4.3 ±0.1 12.2 ±0.1 10.3 ±0.1 0 (degrees) -136.9 ±0.7 97.0 ±0.2 -5.7 ±0.3 Qrr 1.7 0.1 0.1 Qr6 -2.2 0.0 -0.4 Qee 17.4 13.2 142.0 Table 5.1: Polar moments of the three gravitationally lensed images in the simulation shown in Figure 1.2. Comparing the polar moments tabulated for the different types of lensed arcs is revealing. In the weak lensing regime, the ratio of the major and minor axes of the equivalent ellipse gives a measure of the ellipticity of the distorted background galaxy. In the strong lensing regime, we define a shape parameter x by measuring the ratio of the dimensions of the annular sector built from the polar quadrupole moments: _ tangential dimension _ nr\JQee ^ radial dimension 180y^ A lensed feature which is circular produces x ~ 1- F° r the three arcs in the simulation, XAI = 0.2, XA2 = 2.4, and XA3 = 6.8. Radially distorted arcs, like arc A l , show x < lj while tangentially distorted arcs give x > lj as in the case of arc A3. The weakly distorted arclet A2 Chapter 5: A New Lens Inversion Scheme 83 is tangentially distorted, x — 2.4, but not nearly as much at A3. The quantity x m a v serve to distinguish between radial and tangential arcs, based only on their polar moments. In weak lensing analyses, the coordinate frame in which the matrix of quadrupole moments is diagonal defines the principal axes of the equivalent ellipse. In the polar moments scheme, the off-diagonal moment Qrg cannot be interpreted so easily. This is due in part to the incomparable dimensions of the polar Q elements. We can extract some information, nevertheless, from the sign of the Qrg moment. In a spherically symmetric lens, arcs have equal weight inside and outside the circle of radius f, and equal weight clockwise and counter-clockwise from the ray at position angle 9. The off-diagonal moment (5.1e) vanishes. Outside (inside) the the centroid circle, r — f is positive (negative); counter-clockwise (clockwise) from the centroid ray, 9 — 9 is positive (negative). As illustrated schematically in Figure 5.3, the sign of QTQ shows any asymmetry of the image inside the annular sector. A lensed image rotated clockwise about the (9, f) centroid more heavily populates the regions where QTQ > 0. Similarly, Qrg < 0 for images rotated counter-clockwise with respect to the polar centroid. Figure 5.3: The sign of the off-diagonal polar moment Qrg gives a qualitative descrip-tion of how the object is rotated inside the annular sector surrounding the polar centroid (9, f). The asymmetry of arcs A l , A2, and A3 shown in Figure 4.4 agrees with this interpretation Chapter 5: A New Lens Inversion Scheme 84 of the off-diagonal quadrupole moment. The radial arc A l is rotated slightly inside the annular sector, and shows a corresponding Qrg < 0 in Table 5.1. The small tangential arclet is very nearly symmetric inside its annular sector, matched by QTQ which vanishes to the number of digits significant to the calculation. The giant tangential arc A3 appears to be rotated counter-clockwise inside the annular sector, and produces a negative r#-quadrupole moment, as expected. 5.3.1 Applicability of Polar Moments The polar moments introduced here depend critically on the location of the origin of the co-ordinate system in which they are measured. This is not the case in the weak lensing analysis using the more familiar Cartesian quadrupole moments. Despite the generally elliptical form of the mass distribution, gravitational lensing shows a remarkable spherical symmetry. This is because the deflection angle depends on the cumulative mass inside a given homoeoidal shell, and the confocal ellipses (Figure 3.2) about this shell quickly become circular at larger radii. The natural circular symmetry of lensing may make it possible to include the origin of the polar coordinate system as a model parameter. In simple mass distributions with a dominant compo-nent, it is clear that the origin of the coordinate system should coincide with the centre of the mass distribution. Clusters of galaxies, however, are more likely to have bi-modal mass distribu-tions, with two (or more) separate mass peaks. In the galaxy cluster A2390 (Pierre et al., 1996), lensed light is squeezed between two mass concentrations, resulting in a long, straight tangen-tial arc which curls around neither mass centre. Without a clearly defined centre-of-lensing, it is inappropriate to characterize the shape of this arc using the polar moments described here. The cluster A2390 also contains a number of tangential arcs and arclets which surround the bi-modal centre of the cluster (Bezecourt and Soucail, 1997). One or the other mass concentra-tion appears to dominate the lensing of several of these arcs. The polar moments of this subset of arcs, measured with respect to the centre of the dominant mass, could be used to constrain the mass distributions independently. The long straight arc, and other arcs which lens equally about both mass centres, should still appear in the subsequent lens model. While these arcs cannot be used to constrain the model, they can be used to check its consistency. Chapter 5: A New Lens Inversion Scheme 85 A strength of the polar moments characterisation of the position of a lensed arc is the manner in which the annular sector smoothly deforms into the equivalent ellipse of weak lensing analyses. The struts extracted from the eigenvalues of the matrix of quadrupole moments of weakly lensed background galaxies, which we interpret as the major- and minor-axes of an equivalent ellipse, are closely matched by the radial bar between r ± y 3Qrr and circular arc between 6± ySQgg in the weak lensing limit. The Cartesian approach breaks down as the strength of the lensing increases and the arcs become less and less elliptical in shape. The polar approach carries on smoothly into the strong lensing regime, breaking down only at the the singular case of an Einstein Ring where the angular moment 6 loses its meaning. The slightest irregularity or eccentricity in the mass distribution breaks this symmetry, however, and a complete Einstein Ring does not form. The viability and applicability of polar moments have yet to be fully explored. Neverthe-less, this small collection of polar moments (5.1) characterizes the position of simulated lensed arcs quite well. Four constraints can be extracted from each image: f, 6, Qrr, and Qgg. Equiv-alently, we can consider f, 6, the shape parameter x> a n d one dimension, \JzQTT. Including the off-diagonal moment Qrg as a constraint is dubious, although matching its sign between the observations and simulations provides an additional check of the model. 5.4 O p t i m a l S o u r c e R e c o n s t r u c t i o n The lens inversion scheme described here first establishes the geometry of the lens and the deflector mass distribution by matching the positions of the lensed arcs. In fixing these model parameters, the amplification factor p and magnification \p\ throughout the lens is completely determined. We can therefore strip away the effect of the magnification of the lens in the observed lensed features. Identifying lensed pixels in the data requires setting a threshold, the quantity I0 in (5.1), at the level of the background sky plus the rms noise. In the application of the inversion scheme to the clusters MS 2137 and MS 1455, we first identify the "skeleton" of the lensed arcs, those pixels which clearly contain lensed light. Then the neighbouring pixels with signals still larger than the background are tagged. In arcs showing significant Chapter 5: A New Lens Inversion Scheme 86 substructure, low surface brightness pixels completely surrounded by lensed pixels are also selected. We postulate that these darker pixels correspond to darker regions of the background source, rather than to cases of extreme twisting and distortion of the background source over very small regions of the deflector plane. The signal recorded in each lensed pixel is divided by the local value of the magnification to "de-magnify" signal. The de-magnified flux of the images is the flux the source as it would appear in the absence of lensing. By tracing the centre of each lensed pixel back to the source plane via the lens equation, we produce a set of points on the source plane, each carrying a surface brightness. The surface brightness distribution of the hidden background source is reconstructed by interpolating between these points. Ideally, we would like to produce a picture of the background source which looks like the observations of a distant galaxy: a pixelised image, with constant surface brightness in each pixel. We adopt a constant source pixel size A s to match the format of the data collected in CCD observations. The interpolation of the data across the source plane could be done over pixels of varying size, but we choose the simpler approach of constant pixel spacing at this time. At some point, it would be interesting to investigate adaptive interpolation of the data, as well as super-resolution of the observations back in the image plane. In CCD observations, the sky in sampled uniformly at a fixed angular resolution, mapping the sky onto the CCD. This generates a contiguous picture of the sky. This is not the case when we try to observe the background source through the gravitational lens. The source plane does not inherit the uniform and contiguous sampling of the image plane for several reasons discussed below. The observations collected by a CCD are an array of numbers which we interpret as the average value of the surface brightness over the region of sky which maps onto each pixel. After de-magnifying and tracing a ray back to the source plane, we find the surface brightness ps at point y given by the lens equation. There are two obvious choices for extrapolating this information over a pixel. Pixels on the source plane can be set to a size As defined by DsAs = DdA . Chapter 5: A New Lens Inversion Scheme 87 Through this relation, pixels of size A s on the source plane and pixels of size A on the deflector plane represent the same physical dimension. Even in the absence of lensing, pixels of this size cannot contiguously cover the source plane (Figure 5.4a) because the physical area of the source plane behind the image is larger. Lensing merely distorts the rays, as shown in Figure 5.4b. Only in the most contrived situation could the distortion produce a perfectly covered region source source plane. A second possibility is choosing the inter-pixel spacing A s to have the same angular size as image plane pixels: A S = A . These pixels cover the source plane (Figure 5.4c) but as illustrated in Figure 5.4d, distortions introduced by gravitational lensing destroy the uniform coverage, multiply sampling some re-gions and perhaps missing others altogether. Furthermore, this approach adds a significant amount of information to the system by assigning a value to the surface brightness of a region of source plane much larger than the corresponding region on the deflector plane where the information is collected. Regions of the source plane near the caustics of the lens are greatly over sampled by the backwards-traced rays. This comes at the expense of other regions of the source plane as the number of rays is fixed in the observations. There is no guarantee that the irregular sampling of the source plane sufficiently covers the background light distribution to produce a contiguous image, no matter what size A s is chosen for the pixels. Uniform sampling of the source plane requires solving the lens equation for rays radiating from a given set of coordinates. This is precisely the gravitational lensing problem which has no solution, and why we must resort to forward modeling. The source reconstruction technique employed in the LensMEM inversion algorithm of Wallington et al. (1996) alleviates the problem of non-contiguous source reconstruction. As described briefly in Chapter 2, each pixel in the image is cut into two triangles. The corners of these triangles are traced back to the source plane and re-connected, forming a triangle on the Chapter 5: A New Lens Inversion Scheme 88 cy Cy \ C~7 : \ Cy Cy cy <=> cy Cy: cy: a y . cy cy CJ; CJ cy' o cy cy cs 0 7 a? cy 5.4a: Source plane pixels have the same physical size as deflector plane pixels. In the absence of lensing, the source plane is not covered. 5.4b: Deflection of rays carrying pixels of equal physical size does not produce a completely cov-ered source plane. 5.4c: Source pixels with the same angular size as those on the deflector plane cover the source plane, but extrapolate the deflector plane infor-mation over a larger physical area. f 5.4d: Deflection of rays carrying pixels of equal angular size destroys the uniform source plane coverage. F i g u r e 5.4: Neither setting the source plane pixels to the same physical size nor the same angular size as the deflector plane pixels satisfactorily covers the source plane. Chapter 5: A New Lens Inversion Scheme 89 source plane. As each vertex is traced several times, once for each triangle to which it belongs, contiguous coverage of the source plane is assured. Because of the extreme distortion induced by lensing, there may be regions of the source plane covered by more than one "sheet" of triangles, so some interpolation or smoothing must be done. The LensMEM routine discretises the source plane into a uniform rectangular grid and calculates the overlap between each projected triangle and each rectangular source pixel. The reconstructed source is found by summing over the rectangles rather than the triangles. Calculating the overlaps is far from trivial and could greatly reduce the speed at which the source can be reconstructed. The interactive nature of the program lens used here would be lost. Furthermore, connecting the vertices of the projected triangle with straight lines is only an approximation, and the overlaps calculated may be spurious. This problem remains no matter what source resolution is chosen. We try to find a compromise which neither ignores nor exploits the amount of information we obtain after de-magnifying the lensed pixels and tracing each to its origin on the source plane. First we establish the pixel size A s which ably reproduces the appearance of the source. A series of source reconstructions based on simulated observations of the "PN" test pattern illustrates our approach. The "PN" background source is lensed through a simple PID lens, producing the simulation shown in Figure 5.5 where the pixel size is A = 0'.'414. The 274 pixels containing lensed light are de-magnified and the centres are traced back to the source plane. We interpolate the data onto the source plane at four different resolutions. Figure 5.6a shows that choosing a high resolution (small A s ) for the source plane very nearly gives each ray its own pixel. The large over-sampling of points near the caustics is evident in this Figure, but it is a poor reproduction of the original source. By decreasing the resolution of the reconstructed source distribution in Figures 5.6b and 5.6c, better and better reproductions of the source are generated. The structure of the source is almost lost in Figure 5.6d with further increases in the source pixel size. Constrained by our choice of uniform source pixel size, it is difficult to decide which re-construction in Figure 5.6 is the most faithful reproduction of the "PN" test pattern. Instead, we are guided by a another constraint, that the amount of information in the source should be Chapter 5: A New Lens Inversion Scheme 90 Figure 5.5: The "PN" test pattern is lensed through a single PID deflector. Figure 5.6 shows the source distributions reconstructed from the data in this image. Chapter 5: A New Lens Inversion Scheme 91 - 4 - 3 a r c s e c o n d s 5.6a: Source reconstruction with 2950 pixels. 5.6b: Source reconstruction with 460 pixels. 5.6c: Reconstruction with 255 pixels, roughly the same number as the 274 incoming rays. 5.6d: Reconstruction with only 110 pixels. The "PN" is almost indistinguishable. F i g u r e 5.6: Reconstruction of the "PN" test pattern on the source plane at various resolutions from the lensed images in Figure 5.5. Dotted lines trace the caustics on the source plane. Small dots mark the points where the backwards-traced rays pierce the source plane. Chapter 5: A New Lens Inversion Scheme 92 comparable to the amount of information available in the observations. We choose the resolu-tion of the source distribution so that the total number of pixels in the result is comparable to the images in Figure 5.5. The reconstructions in Figures 5.6a, 5.6b, 5.6c, and 5.6d are arrays of 2950, 460, 255, and 110 pixels, respectively. We measure the total area of the source reconstruc-tion because pixels without a signal are just as important as those with a signal (the hole in the "P", for instance.) Under this criterion, the reconstruction which neither interpolates too much, nor throws away valuable information is the example containing 255 pixels, comparable to the 274 lensed pixels in the lensed images. The size of each pixel As is determined by dividing the angular extent of the reconstructed source by the width of the image, in pixels. The width of Figure 5.6c is 6'.'70 across 21 pixels, giving each pixels a size A s = 0''32/pixel. This is smaller than the size A = 0'.'414 on the image plane, but larger than the size (Dd/DS)A = 0'.'23 which represents the same physical scale on the source plane as pixels on the deflector plane. Once we have decided on the pixel size at which to display the result, we assign a value to each pixel. In a perfect lens model with error-free data, rays which converge onto a common source point should carry equal surface brightnesses. In practice, though, the model is only a approximation to the deflector mass distribution and the strength of the magnification removed from the observed signal may be incorrect. As well, the data may be contaminated by errors and by extraneous light from luminous galaxies in the cluster. A source pixel pierced by a single backwards-traced ray coming from a lensed pixel with flux I, error a, and magnification \p\, is assigned a flux S = (When a ray coming from this source passes through the lens, it is magnified by |/x|, producing the observed flux I.) Multiply-imaged source pixels require a different strategy. Suppose rays from k lensed pixels, each carrying an observation 7j, error <7j, and magnification strike the same source pixel. We assign this source pixel a flux S which minimizes the error, , in the reconstructed signal, where the number of rays traced from the data. In the examples here, there are 274 lensed pixels in (5.4) I Chapter 5: A New Lens Inversion Scheme 93 This choice of S is better than simply averaging the incoming |/ij|_ 17j, but does not take into account the smearing of the lensed objects in the observations. Kochanek and Narayan (1992) discuss how to include a point-spread function in their LensClean algorithm. Each pixel in the sources in Figure 5.6 is assigned a flux following this procedure. The discretized sources are passed through the lens once more, to check the consistency of the reconstruction scheme. Figure 5.7a shows that the lensed features first modeled in Figure 5.5 can be almost perfectly recovered from the fragmented source containing 2950 pixels. This is to be expected, as the small pixels in the reconstruction closely approximate the points where the collection of rays pierces the source plane. As the size of the source pixels increases, the information contained in each ray is interpolated onto larger regions of the source plane. Larger sources generate larger lensed images, as shown in Figures 5.7b and 5.7c. The number of pixels containing lensed light increases, owing to the larger background sources. In the crude 110 pixel reconstruction of the source, information is lost because of the least-squares averaging (5.4). Furthermore, the inaccurate values for flux found through (5.4) appear over a region of the source plane which is significantly larger than the original analytic source. As Figure 5.7d shows, many extraneous structures are produced in the lensed images. The data travel in a loop from the image plane, to the source plane, and then back to the image plane. This suggests the possibility of iterating the source reconstruction phase of the inversion. Comparing the re-lensed images with the original lensed features reveals which pixels have been introduced by the interpolation on the source plane. It may be possible to eliminate the extraneous pixels in the re-lensed images by removing the region of the source which gives rise to these spurious pixels. Because many parts of the source are multiply-imaged through the lens, however, the removal of one source pixel may result in the elimination of pixels in the re-lensed image in addition to one spurious pixel which lead us to this region of the source. We would then be faced with the decision of how many viable pixels to sacrifice for each spurious pixel. Chapter 5: A New Lens Inversion Scheme 94 arcseconds arcseconds 5.7c: 255 source pixels 5.7d: 110 source pixels F i g u r e 5.7: The pixelated sources in Figure 5.6 are passed back through the lens. Because larger source pixels spread the information out over larger regions of the source plane, the lensed images contain extraneous lensed light not present in the original simulation shown in Figure 5.5. Chapter 5: A New Lens Inversion Scheme 95 A better solution is to iteratively subdivide and refine regions of the source where the data are concentrated. As stated earlier, in this thesis we make only the first approximation of using a constant pixel size across, the entire source plane, over sampling in some regions, and undersampling in others. Any pixel in the source array which is not pierced by an incoming ray is set to 0. There are two explanations for an undetermined pixel in the source reconstruction, however. The first is that no lensed light is detected originating from this region of the source plane. The pixel is empty because it covers a region of the source plane which emits no light. A second explanation is that there is no solution to the lens equation for a ray emitted from this region of the source plane. Even if the pixel contained a non-zero signal, we could not see it through the lens. We have an incomplete view of the source plane. The gaps left by these hidden pixels can be filled by smoothly interpolating between the known data, but we run the risk of assigning a signal to the empty pixels which are visible through the lens. While it is possible we could detect light from an empty pixel, the critical lines and parity inversions produced in highly distorted lenses may be hiding it from our eyes. We can only distinguish the empty pixels from the hidden pixels by forward modeling the source distribution through the lens. This becomes more apparent in the source reconstructed behind the lens in MS 1455 in Chapter 6, where only that part of the source lying under a caustic is visible through the lens. Finally, it is difficult to ignore a serendipitous result, especially one which illustrates so well the complex relationship between the source plane and the image plane. By convolving the "PN" lensed images with a point-spread function and then adding noise to the image, a significant number of the pixels in the image are non-zero. The routine which reconstructs the source re-traces all pixels above the threshold Ia — 0, and inadvertantly interprets all of the non-zero pixels as lensed features. Virtually the entire image plane is reconstructed in Figure 5.8, and we see many of the characteristics of the distortion induced by gravitational lensing. The regions of the source plane near the ovoid and astroid caustics are greatly over-sampled. The caustics are visible in the concentrations of data points. The over density of information near the caustics is matched by empty regions of the source plane which are simply not sampled by Chapter 5: A New Lens Inversion Scheme 96 the image, or equivalently, not visible through the small field-of-view in which we observe the lens. We also see how the "sheet" of light visible at the image plane is twisted and folded on the source plane. - 1 0 - 5 0 5 a r c s e c o n d s F i g u r e 5.8: The software inadvertently reconstructs a source for virtually the entire image plane after the "PN" lensed images are convolved with a point-spread func-tion and noise is added. Dots mark the points where re-traced lensed rays pierce the source plane. The over- and under-sampling of the lens equation is obvious. The astroid and ovoid caustics marked by dotted lines in Figure 5.5 are visible. The limited field-of-view on the deflector planes makes the outer regions of the source plane invisible to the observer. The new lens inversion scheme addresses many problems inherent in gravitational lens inversion. Some of these problems can be traced to the fundamentally unstable nature of inverse Chapter 5: A New Lens Inversion Scheme 97 problems, particularly the large changes in the lensed images resulting from small changes in the position of the source which inadvertently put the source under a caustic of the lens. Other problems arise only when we try to implement the inversion scheme. The non-uniform sampling of the source plane and the disjointed source distribution produced are problems we did not foresee. While care is taken to arrive at well-motivated and robust solutions to these problems, the bulk of the testing has been done on simulated images. In Chapter 6, we test the lens inversion scheme on observations of two galaxy-clusters which show gravitationally lensed features. Chapter 6 Gravitational Lensing in MS 2137 and MS 1455 The strength and usefulness of a lens inversion scheme, or in fact, any mathematical model of natural phenomena, is how well it can predict and reproduce the observations. It is in the observations, after all, that the truth lies, regardless of the elegance of the mathematics which tries to model the behaviour. In this Chapter, we test the inversion scheme described in Chapter 5 on observed lensing systems. First, to test the validity of the scheme, we examine the collection of arcs observed in M S 2137—2353. This lens has already been "inverted" by other groups and allows us to compare our new results against existing ones. Then, confident that the inversion scheme itself is viable, we attempt to invert the gravitational lens observed in M S 1455+22. The strong lensing behaviour of this galaxy cluster has not been studied before. It provides both a second data set on which to test the new inversion scheme, and a new and quite unique example of gravitational lensing. 6.1 The Collection of Arcs in MS 2137 The galaxy-cluster M S 2137 lies at redshift = 0.313. It is visible in the optical, as well as being an X-ray source. Ground-based observations revealed a collection of gravitationally lensed features, including the first recorded radial arc (Fort et al., 1992). Recently, observations from the HST have revealed many more arcs, as well as highly resolved structure within some 98 Chapter 6: Gravitational Lensing in MS 2137 and MS 1455 99 of the previously observed arcs. This data set, available to the public from the STScI,1 is the basis of the work described here. The image in Figure 6.1 results from processing ten WFPC2 800 x 800 pixel frames at a resolution of O'.'100/per pixel with STSDAS routines. The frames, two with exposures of 1900 seconds, and 8 at 2300 seconds, are first aligned and then combined taking the median value the individual frames. This removes the cosmic ray streaks in the data. The light from the large central cD galaxy G l and a smaller galaxy G7 are digitally removed using STSDAS routines which model each galaxy as a sequence of concentric elliptical shells. The resulting image is shown in Figure 6.2. A large collection of arcs become visible through this process. Most obvious is the giant tangential arc, AO. Not showing the signature of gravitational lensing nearly as much, but confirmed by lensing models, are tangential arclets A2, A4, A6, and radial arc AR. The labels follow the notation of Hammer et al. (1996) which itself is based upon the models of Mellier et al. (1993), although the latter could not distinguish arcs B l and A6 in their ground-based data. 6.1.1 Existing Models of MS 2137 The lensed features in MS 2137 have been observed and analysed by different instruments and groups, first by Mellier et al. (1993) (hereafter M93), then by Miralda-Escude (1995), and later by Hammer et al. (1996) (hereafter H96). The models produced by these groups are very similar. The earlier model places a PID mass coincident with the central cD galaxy, perturbed by a mass distribution with density following the light profile of the cD. Miralda-Escude refines the model of the cD galaxy, replacing the observed light profile with a parametric function. The model presented in H96 replaces the pseudo-isothermal profile with a /3-profile mass density, p ~ (1 + (r /r c ) 2 ) - ^ with rc and 3 model parameters (3 = 1 models the PID.) The orientation and eccentricity of the dark-matter halo is allowed to vary slightly from the observed values of 1 Based on observations made with the NASA/ESA Hubble Space Telescope, obtained from the data archive at the Space Telescope Science Institute. STScI is operated by the Association of Universities for Research in Astronomy, Inc. under the NASA contract NAS 5-26555. Chapter 6: Gravitational Lensing in MS 2137 and MS 1455 100 T r~ f. :i at •••^-••-••1""-—.^"^•-vr r——i r r——i —r r i — * p \ i — T T * r j J J i i i i ;f. I 11 it ,i- J L „ „r L 'mi i I i i i i _1 i & - 2 0 - 1 0 0 10 20 arc s e c o n d s Figure 6.1: WFPC2 image of MS 2137. The grayscale is linear. The giant arc AO and arclets A2 and A4 are traced to one source. A second source produces the radial arc AR and the arclet A6. Our models suggest that Bl, BR, Cl , and CR are also lensed arcs. The central cD galaxy Gl and the galaxy G7 which obscures part of arc AO are removed in Figure 6.2. North is at an angle of -38.5° relative to the positive ar-axis, based on the HST orientation. Chapter 6: Gravitational Lensing in MS 2137 and MS 1455 101 arcseconds Figure 6.2: Processed WFPC2 image of MS 2137. The grayscale is linear. The central cD galaxy Gl and a smaller galaxy G7 obscuring the upper end of arc AO have been removed, revealing more of the structure of the arcs. The irregular under-densities left behind are numerical relics from the removal of the galaxies. Chapter 6: Gravitational Lensing in MS 2137 and MS 1455 102 the central cD galaxy. The model in M93 also includes a mass component matching the cluster member G7. All models require two sources at redshift zs ~ 1. The first source SI produces the giant arc AO, with counter images A2 and A4. A second source is necessary as the radial arc AR appears on the same side of the center-of-lensing as the giant arc, and cannot also originate from SI. This second source produces a counter-image, arc A6. The lens geometry is fixed in M93 using arcs AO, A2, and A4, all coming from source SI. The model parameters are systematically adjusted to minimize a x 2 statistic which measures the distance between the centres of the pre-images of each of these three arcs on the source plane, assuming A2 and A4 are only weakly distorted. The inversion scheme of H96 proceeds by tracing the bright knot of light visible in arcs A2 and A4 and twice in arc AO (once on either side of the critical line) back to the source plane. The parameter space is searched for a set of model parameters which minimises the distance between the pre-images. Then each lensed image is pulled back to the source plane and combined to produce a single source ("reconstructed source image", RSI) responsible for all three lensed images. The RSI is rather irregular in shape and surface brightness. The details of the reconstruction, particularly how the problems of source resolution and contiguity are addressed, are not discussed in H96. The models presented in the M93 and H96 analyses of MS 2137 appear to have 16 and 15 parameters, respectively, although the nature of the background source in M93 is not discussed. The values of these parameters are listed in Table 6.2 for comparison with the model we produce below. The redshift of the lens plane and the parameters describing the central cD galaxy are set by the observations. The centre of the dark-matter halo is assumed to coincide with the centre of the cD galaxy. This leave some 7 and 5 parameters free, respectively. The centres of the preimages of the four bright knots provide eight contraints. To check the consistency of the models, a second source S2 is added at the same redshift as SI, which successfully produces arcs AR and A6 in both cases The reconstructed S2 source of H96 is generally elliptical, again with an irregular surface brightness. As Hammer et al. demonstrate, the arc B l and a very small radial arc BR parallel to arc AR are consistent with Chapter 6: Gravitational Lensing in MS 2137 and MS 1455 103 the presence of a third source. Attributing these features to gravitational lensing leads to the prediction of a third source at redshift slightly less than 1. 6.2 Inversion of MS 2137 using the New Lens Inversion Scheme Our lens inversion scheme contains characteristics of the inversion methods of both Mellier et al. and Hammer et al.. We first build a parametric model of the lens with an elliptical background source, similar to M93. The source is later replaced by one generated from the lensed data in the observations, similar to H96. We have the benefit of two existing models as a start for the values of the model parameters. 6.2.1 Establishing the Lens Geometry and Mass Distribution Following the lens inversion scheme outlined in Figure 5.1, we first establish the lens geometry and the deflector mass distribution parameters based on the positions of arcs AO, A2, and A4. We flag the pixels in the data containing lensed light, following the prescription described in Section 5.3, and tabulate the polar moments and shape parameter x> which are listed in Table 6.1. The O'.'lOO resolution of the HST images gives a maximum error of Sr ~ 0'.'050 in the radial positions of the arcs. Errors in position angle are described by (5.2). As expected, X < 1 for the radial arc AR and x 3> 1 for the giant tangential arc AO. Arc r d Qrr QrO Qee X (arcsec) (degrees) AO 15.40 ± 0 . 0 5 -46.7 ± 0 . 2 0.15 -2.1 222.7 10.4 A2 13.07 ± 0 . 0 5 -152.1 ± 0 . 2 0.37 0.2 11.9 1.3 A4 19.45 ± 0 . 0 5 82.6 ± 0 . 1 0.09 -0.5 19.9 5.0 AR 5.43 ± 0.05 -53.6 ± 0 . 5 1.07 1.5 5.7 0.2 A6 24.20 ± 0.05 112.5 ± 0 . 1 0.19 0.4 3.4 1.8 Table 6.1: Polar moments of the observed arcs in MS 2137. The first three arcs in the list are used to build the lens model. The last two are used to check the consistency of the model. Chapter 6: Gravitational Lensing in MS 2137 and MS 1455 104 We postulate a PID lens supplemented by a mass distribution following the cD profile. The orientation and eccentricity of the cD is set to match the observed values and the scale parameters of the cD mass model are fit to the observed profile through an iterative parameter estimation scheme. We allow the centre, orientation, and eccentricity of the PID to vary slightly from the central cD. The core radius rc of the PID is one of the model parameters we fit using the polar moments of the simulations. A single elliptical source SI at redshift zs ~ 1 is observed through the lens, and the polar moments of the simulated arcs are measured. The model parameters, including the source redshift, are adjusted to reproduce the observed polar moments of arcs AO, A2, and A4. Differences in the radial centroid A f = f0(,s — f s j m give a measure of the radial difference between between the observed (obs) and simulated (sim) centroids of the lensed images, in arcseconds. The quantity - - - TT A* = (0obs - 0sim)—fobs (6.1) measures the tangential difference between the polar centroids, also in arcseconds. It is unclear what statistic should be calculated to compare the values of the two sets of moments. A penalty function P which measures the distance between the centroids P = (Af)2 + (At)2 takes no account of the shapes of the arcs. Additional quadrupole moments can be included, perhaps weighted to increase or decrease their importance: P = (Af)2 + (At)2 + \(Qrr,obs ~ Qrr,sim)2 + P(Qee,obs ~ Qoe,sim)2 While the penalty function returns a single statistic describing each configuration, its physical interpretation and the criterion for choosing the weights A and p. remain unclear. Rather than minimizing an arbitrary statistic like P to measure the goodness-of-fit of a model against the data, we adopt a policy of matching the polar centroids to within the errors associated with the measurements. Then we make further small adjustments, particularly in the shape of the background elliptical disk, to match the Qrr and Qgg quadrupole moments. Chapter 6: Gravitational Lensing in MS 2137 and MS 1455 105 Further analysis of the polar moment characterization of lensed arcs, and also what is meant by a "better" model, is needed before a viable statistic that can distinguish subtle improvements of one model over another, can be chosen. The set of model parameter values we chose for the final simulation of MS 2137 are listed in Table 6.2. For comparison, we include the parameter values which are listed or can be estimated from the models found in M93 and H96. The dark-matter halo described in H96 follows a /3-profile mass density, so direct comparison of the mass models is impossible. Also, there is no analytic source model in the lens inversion scheme of Hammer et al., who instead reconstruct the source distribution from the lensed data. The polar moments of the simulated arcs are listed in Table 6.3, and Table 6.4 tabulates the differences between the observed and simulated polar centroids. The centroids of the three arcs AO, A2, and A4 which are used to contrain our model each match to within 1''0. The shape parameter x of the simulated arcs closely matches the shapes of the observed arcs. The sign of the second moment QR$ matches in all cases, providing a further, though qualitative, check on the model. The simulation is shown in Figure 6.3, along with the annular sectors built from the second polar moments. While surface brightness of the arcs plays no part in the fitting process, we adjust the uniform surface brightness of the background source SI to mimic the flux in the weakly magnifified A2 and A4 arcs. Only the portion of the 800 x 800 pixel HST observations which is visible in the field-of-view of the simulations is shown in Figures 6.1 and 6.2, so direct visual comparison between the observations and the simulation can be made. Without changing the model parameters, a second source S2 is added at the same redshift as the first, as in M93 and H96. The second source produces two new lensed images, modeling the observed radial arc AR and its counterimage, A6, which are visible in Figure 6.3. The polar moments of the simulated AR and A6 arcs are included in Table 6.1. As shown in Table 6.4, the polar centroids of the two new simulated arcs closely match the observed arcs. While the magnitudes of the second moments QRR and QQQ of the simulated arcs are different than the Chapter 6: Gravitational Lensing in MS 2137 and MS 1455 106 c p «3 CP u CP CP -a o C O o C O o o o o C O C O C O I o L O o L O o T> CP SH SH O H J o CP c n CP -rt 03 c p SH a> o SH rt o Q o _0J c n o SH OH I Q o Q i—i PH Q Q o o L O 0 0 -H 11 O l L O r- 6 6 I + 1 L O i—l C M O ° „ -H ° ° O L O I I I C O I C M o L O C O oo o II -H L O o >o o O C N 00 H + 1 O i-l o o o o o + 1 o 00 CM 1—1 CO II o o o o o o o L O o o o C O i—( o L O L O o C O L O C O L O C O o rt rt CP SH -s CP CO CP CP feb CP rt o •s s TH o rt CP o o CP CP CP a B B CP o CM C M C O C O C M C O C O T J CP 7§ SH CS OH CP CO C O O -H C M I I I L O o cT C M I C O -00 t--CP X ) o CP CP SH rt o co o CP 03 o CP u •s CP CP o -H O S O 00 O I I I C O C O .71,5. 0.65 o OS 0.74 C M 1 1 L O C O L O C M u a > co o SH ci X cS u o 'c? a i • i-H a CP 03 CP CP fcb CP S-rt o rt S o o fl °(H o Ex 23 tcj 1) and the radial arcs (x < 1)- The corresponding moments of the simulated arcs are listed in Table 6.6. There is very good agreement in the Chapter 6: Gravitational Lensing in MS 2137 and MS 1455 110 a r c s e c o n d s Figure 6.4: Two additional pairs of arcs, Bl-BR and Cl-CR, and annular sectors built from their polar moments listed in Table 6.5. Chapter 6: Gravitational Lensing in MS 2137 and MS 1455 111 i 1 1 1 1 1 r I i . I i i i i I i i i i I , , i i I i i I -20 -10 0 10 20 arcseconds Figure 6.5: Simulation of the Bl-BR pair of arcs. The background elliptical source lies at redshift 0.900. Chapter 6: Gravitational Lensing in MS 2137 and MS 1455 112 - 2 0 -10 0 10 20 a r c s e c o n d s Figure 6.6: Simulation of the C l - C R pair of arcs. The background elliptical source lies at redshift 1.725. Chapter 6: Gravitational Lensing in MS 2137 and MS 1455 113 positions of the polar centroids. The greatest discrepancy between the simulations and the observations in the appearance of radial arc BR in Figure 6.5. This may be due to the fact that we model the source as an elliptical disk. The appearance of B l in the observations, however, clearly suggests a background spiral galaxy. The radial arc BR may originate from one of the spiral arms of this galaxy which lies under the radial caustic, rather than the bright central bulge mimicked by our elliptical disk. Arc r (arcsec) e (degrees) Qrr Qre Qee X B l 22.77 ± 0.05 106.0 ± 0 . 1 0.21 0.0 6.1 2.1 BR 6.00 ± 0.05 -63.9 ± 0 . 5 0.12 0.4 3.1 0.5 C l 27.07 ± 0.05 -124.3 ± 0 . 1 0.09 -0.1 1.4 1.9 CR 7.26 ± 0.05 28.5 ± 0.4 0.11 0.0 0.5 0.3 Table 6.5: Polar moments of the two pairs of arcs, Bl-BR and Cl-CR, visible in Fig-ure 6.4. Arc r (arcsec) d (degrees) Qrr Qre Qee X B l 22.7 ± 0 . 2 105.0 ± 0.5 0.1 -0.2 7.1 3.3 BR 5.7 ± 0 . 2 -63.6 ± 2 . 0 0.9 0.4 28.6 0.6 C l 27.1 ± 0 . 2 -124.5 ± 0 . 4 0.1 -0.2 4.6 3.2 CR 6.7 ± 0 . 2 29.1 ± 1 . 7 1.3 -0.9 3.1 0.2 Table 6.6: Polar moments of the simulated arcs Bl-BR shown in Figure 6.5 and arcs Cl-CR shown in Figure 6.6. 6.2.3 Reconstructing the Sources Having established the geometry of the lens and the deflector mass distribution, we are able to remove the effects of magnification from the lensed arcs. The resulting objects, while still distorted by the lens, have the flux of the background source as it would appear in the absence of lensing. Chapter 6: Gravitational Lensing in MS 2137 and MS 1455 114 With the established lens model, we pull the observed lensed data back through the lens to the source plane, de-magnifying it along the way. The two sources SI and S2 are multiply imaged, and data from different arcs overlap on the source plane. Because of the folding of the caustic surfaces that occurs, data from a single arc, like the giant arc AO, may be folded over. By either means, the uniformly spaced data on the image plane become a highly non-uniform sampling of the source plane. To reconstruct the source responsible for the lensed features in an optimal way requires both a choice of resolution and a way of assigning a surface brightness to each source pixel. As discussed in Chapter 5, we tackle each problem separately. The pixel size on the source plane is set so that the number of pixels in the reconstructed image is comparable to the number of lensed pixels in the observations. To assign a value to the flux in each source pixel, we follow these steps: 1. Source pixels which are not pierced any backwards-traced rays are assigned a value of NaN and appear in the result as pixels of 0 flux. 2. Cells pierced by a single backwards-traced ray are assigned the de-magnified flux. 3. Cells pierced by two or more rays are assigned a flux which minimizes the error in the image, according to the misfit (5.4). The three panels in Figure 6.7 show reconstructions of each source individually and then side-by-side. The reconstruction of SI responsible for the giant arc AO and the two counter-images A2 and A4 contains 513 pixels. Each pixel is 0'.'074 across. There are roughly 2800 lensed pixels identified in the observations of AO, A2, and A4, however. The reason for this discrepancy is that arcs A2 and A4 just barely support of reconstruction at this resolution. The raw data occur at a spacing of 0'.'100 and the A2 and A4 images are not strongly distorted. The data arriving from arc AO sample the source plane at a much higher density, simply because of the great distortion that occurs. The reconstructed SI source is roughly elliptical, and coincides with the uniform elliptical disk used in the simulations. My result is comparable to that shown in H96. Chapter 6: Gravitational Lensing in MS 2137 and MS 1455 115 c o o * . . > • • • • 11 •' - 1 - 0 . 5 0 a r c s e c o n d s arcseconds F i g u r e 6.19c: (left) Source reconstructed from 172 pixels in the third arc A3 through the N F W + c D lens model. Note the decrease in flux, (right) Its appearance through the lens. Chapter 6: Gravitational Lensing in MS 2137 and MS 1455 140 Source reconstruction and relensing shows that gravitational lensing through these para-metric lens models is not entirely responsible for appearance of the object A3. A more likely explanation appears when the source is reconstructed from the data contained only in arcs A l and A2. Cleaning the larger data set removes most of the spurious pixels that remain in Fig-ures 6.18a and 6.19a, which are based on the A l data alone. Without the data from arc A3 to dilute the highly-magnified region of the source near the intersecting caustics, the reconstructed source in Figure 6.20 still contains an improbable bright patch at the upper right, resulting in the bright, highly-distorted appearance of arc A3 in the right panel of the Figure. Source from A1,A2 data Simulated arcs Linear scale 0 — 90 Linear scale 0 — 240 arcseconds orcseconds Figure 6.20: (left) Source reconstructed by combining the A l and A2 data sets. Clean-ing removes 53 of 544 data. The pixel size and grayscale is the same as Figure 6.18a. (right) Its appearance through the lens. The simulated images of the A1-A2 source are convolved with a psf and noise is added, mimicking the effects of seeings on the observations. The result is shown in Figure 6.21. Con-volution spreads the highly concentrated and tangentially distorted arc A3 over a much larger region, matching quite closely the size and shape of the observed A3 arc shown in Figure 6.11. The extreme over-brightness of the simulated arc is due, in part, to the over-brightness of the source in the vicinity of the intersecting caustics. The arc A3 also appears extremely bright Chapter 6: Gravitational Lensing in MS 2137 and MS 1455 141 because of its location on the critical line. In fact, convolution of the original PID+cD simula-tion, in which a featureless elliptical source is used to constrain the model parameters, produces very similar lensed features in Figure 6.22. The relative intensity of arc A3 compared to the tangential arclet A l is due to magnification of the lens, and not structure in the source. n — 1 — 1 — r - ' — 1 — 1 — • — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — ' — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — n L 1 , 1 1 , 1 1 i 1 , , , 1 , 1 , , 1 1 , , 1 1 . . .H -20 -10 0 10 20 a rcseconds Figure 6.21: Simulation shown in Figure 6.20 is convolved with a point-spread func-tion and noise is added. The position and shape of arc A3 rs very similar to the observations. Chapter 6: Gravitational Lensing in MS 2137 and MS 1455 c o o , / \ \ / / / — I I l _ - 2 0 - 1 0 0 arcseconds 1 0 2 0 Figure 6.22: Convolution of the original PID+cD simulation shown in Figure 6.12, in which the source is a featureless, elliptical disk. Differences in intensity of the arcs are due to magnification. Convolution dramatically changes the shape of arc A3, while the noise introduces some sub-structure in the arcs. Reproduction of the extended shape of the arc A3 by convovling a small, bright object with a point-spread function suggests that the appearance of the arc is partly (or maybe wholly) an artifact of the seeing at the time of the observations. This indicates that point-spread function should be deconvolved from ground-based observations fee/ore running the data back to the lens plane. Space-based observations of MS 1455 may resolve this problem. Chapter 6: Gravitational Lensing in MS 2137 and MS 1455 6.4.6 Origin of the Spurious Data The simulated radial arc A l shown in Figures 6.16 and 6.17 does not extend in towards the centre of the lens nearly as much as the observed radial arc. This can be traced directly to the cleaning process which removes the isolated pixels from the sources shown in Figure 6.14 to produce the sources in Figure 6.15. The spurious data on the source plane are removed based only their de-magnified values, not on their origin in one of the three lensed arcs. As the data are discarded, we flag each datum's origin on the deflector plane. These flagged data are not randomly scattered over the three observed arcs, but are concentrated entirely in the inner end of the radial arc. The positions of the lensed data in the observations are plotted in Figure 6.23, with hollow circles representing points included in the source reconstruction, and solid circles marking those data discarded in the cleaning process. /to • o o o o o o o DO#0 O O O O o o o o o o o o o o o o o o o o * ^ o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o OJO c w o O O O O o o o o o o o o o o o o o Oprf) O O O O O O O O O O O O o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o 0 2 a r c s e c o n d s Figure 6.23: Each lensed pixel forming the radial arc in the observations is plotted in the PID+cD (left) and NFW+cD (right) lenses. Pixels which are removed from the source reconstruction by the cleaning process on the source plane are shown as filled circles, while those data which contribute to the final model are hollow circles. The discarded points originate entirely from the inner end of the radial arc. Chapter 6: Gravitational Lensing in MS 2137 and MS 1455 144 It is significant that the discarded data can be traced to a single arc in the observations, for it shows that the pixels are isolated on the source plane because the other lensed arcs do not require light coming from this part of the source. By discarding some of the data on the source plane and failing to produce only (part of) one lensed arc indicates that these discarded data are only singly-imaged through the lens. This is not the case in the bulge of the source where there is redundancy in the data, and the data drawn from arc A2, for instance, reproduces both arc A2 and arc A l (Figure 6.18b). The models we propose for MS 1455 are incapable of producing lensed features close to the cluster core while maintaining the appearance of the other lensed features. The mass distribu-tion near the cluster core is poorly predicted. The parametric models are under-constrained by the 8 statistics drawn from two definitive arcs A l and A2, if we treat the third arc A3 with some skepticism. There may be another point in the multi-dimensional parameter space which better reproduces the three arcs. We are confident that this other model is not in the immediate vicinity of our model in parameter space, for the modeling software allows us to search the local neighbourhood quite carefully. It is also possible that identified lensed features are inconsistent with a single source model, which requires a change in the nature of the model, not merely a better choice of parameter values. Inconsistencies near the centre of the lens are not completely unexpected, however. Grav-itational lensing, the tool we are using to examine the mass distribution of the galaxy-cluster, does not uniformly measure the mass distribution of the deflector, but only the cumulative projected mass distribution. The lensing behaviour far from centre-of-lensing, but still within the strong regime, is insensitive to perturbations in the mass distribution at the lens centre. For example, changing the model parameters of the cD component of the models of MS 1455 has little effect on the position of the arclet A2 which appears far from the centre. As demon-strated by Miralda-Escude (1995) in the case of MS 2137, these same perturbations can remove the radial arc from the lens altogether, because of its proximity to the centre of the lens. We exploit the sensitivity of radial arcs to the mass distribution at the core in Section 7.3.2, where the presence (or absence) of a radial arc constrains the mass of the cD galaxy. In the case of Chapter 6: Gravitational Lensing in MS 2137 and MS 1455 145 MS 1455, the other small galaxies in the vicinity of the cluster core and the radial arc also play a part in the appearance of the radial arc. In the simulation shown in Figure 6.24, a small SIS with oios = 200km s_ 1is added near the core of the lens at position (0'.'50, —4''05), close to one of the cluster galaxies visible in Figure 6.11. This perturbation significantly alters the appearance of the radial arc, for the highly magnified (and hence, observable) region of the arc is reduced in size. Endless refinements to the positions and shapes of the arcs, particularly the radial arc, can be made by including more masses near the centre of the cluster and near the arcs themselves. However, the small number of statistics we extract from the three lensed arcs does not support the inclusion of further masses. arcseconds arcseconds Figure 6.24: (left) A small SIS, marked M3, is added near the centre of the cluster at position (0'.'50, -4'.'05). (right) The original model without the extra SIS. The highly magnified portion of the radial arc significantly changed. We should not judge the simulated arcs as indications of the success of the lens inversion with equal regard, but rather, weight our criticism in proportion to the radius of the arc. In the extreme, what the model predicts at the very centre of the lens is inconsequential. If the data we discarded were traced to the tangential arclet A2, we would have to seriously reconsider the model proposed, and perhaps wonder if gravitational lensing is occurring at all. As this Chapter 6: Gravitational Lensing in MS 2137 and MS 1455 146 spurious data originates in the inner end of the radial arc, where lensing most poorly measures the mass distribution of the deflector, we can remain confident both about the model and the procedure. The spurious data, considerately, highlight for us the weakest part of the model. 6.5 The Inversions of MS 2137 and MS 1455 The observations of MS 2137 and MS 1455 show the characteristic signatures of gravitational lensing, the extended diffuse objects we can attribute to the distorted images of sources behind the lens. The two-stage inversion algorithm by which we arrived at our simulations is quite successful in producing a lensing model of MS 2137. A number of independent source-arc systems are predicted, and appear to have optical counterparts in the observations. Under our parametrization scheme, the lensing behaviour of the galaxy-cluster MS 1455 is not fully constrained by the three identified arcs, because of the effects of seeing on the appearance of the arcs. The two-stage lens inversion scheme presents us with families of lens models with only a small number of degrees of freedom. The under-constrained PID+cD and NFW+cD models we choose do not fully explain the anomalous diffuse objects observed in the cluster. Nevertheless, the ability of our numerical simulations to reproduce the lensed features in both cases is compelling evidence that lensing is occurring around and through these galaxy-clusters. Chapter 7 Conclusions The goal of this thesis is straight-forward: to introduce a new scheme for inverting strong grav-itational lenses. To fully describe the steps of this algorithm, and to appreciate the differences between this and others in use in the astronomical community, requires the development of some gravitational lensing formalism, so that the algorithm itself can be stated simply and clearly. Implementation of the proposed inversion scheme is also an important part of this work. What differentiates this thesis from being merely a mathematical treatise is the application of the method to real data, something absolutely critical to astronomy. 7.1 The New Inversion Scheme The inversion scheme presented in this thesis appears similar to other current approaches that build up a deflector mass distribution and forward model the lensing (Mellier et al., 1993; Kneib et al., 1995; Hammer et al., 1996), as opposed to those methods which distill a discrete representation of the mass distribution on the deflector plane out of the astronomical data (Kochanek and Narayan, 1992; Wallington et al., 1996). The key feature of our routine is decoupling the magnification induced by the lens and the structure of the surface brightness of the background source in the observed lensed features. This is done by establishing the geometry of the lens and finding a set of model parameters describing the deflector components based on the positions in the image plane of the lensed features. We use polar moments, artificial statistics built from the data, to quantify each arc's location on the image plane and its characteristic shape. The model parameters are adjusted, interactively with the lens modeling software lens, to match the polar moments of the simulations and the observations. The values chosen for the 147 Chapter 7: Conclusions 148 parameters fix the geometry of the lens and the mass distribution on the deflector plane, which together determine the magnification of the background source plane visible through the lens. This allows us to remove the magnification from the data, leaving the lensed images distorted, but carrying the flux of the source as it would appear in the absence of lensing. Because light passes predictably through the lens in both directions, we can trace the observed rays of light back to their origins on the source plane. In principle, multiply-imaged arcs should project back to the same neighbourhood on the source plane. However, errors in the data, and more likely, inadequacies of the model generally conspire to upset the simplicity of the scheme. Instead, we combine the beams to minimize the error on the image plane, after deciding on a resolution supported by the data. Our strategy is to choose a source plane resolution so that the total number of pixels in the the neighbourhood of the source, both bright and dark, is comparable to the number of pixels containing lensed light that are available in the observations. This is quite successful, although the distortion of the lens serves to hide some parts of the source plane, while greatly oversampling others. A final step in the inversion scheme is a consistency check. We shine the reconstructed source back through the lens and check that (i) the observed arcs are reproduced, and (ii) no new arcs are introduced, arcs which do not appear in the observations. Because of the non-linear nature of lensing, we cannot definitively predict the result of this re-lensing. Observations of the galaxy-cluster MS 1455 from the CFHT revealed a diffuse, radial object, in addition to the tangential arclet seen previously by LeFevre et al. (1994). This data set offers us the opportunity to demonstrate the existence of a radial arc in the core of the cluster. This radial arc has not been identified before, and inverting the lens with the untested lens inversion scheme is risky proposition. Thus, we turn first to the previously-studied gravitational lens in the galaxy-cluster MS 2137. Mellier et al. (1993) construct a model for MS 2137 using ground-based observations consisting of two sources behind a PID dark-matter halo, perturbed with a cD galaxy. The giant arc visible in the observations is produced by one source, along with several counter-Chapter 7: Conclusions 149 images. A observed radial arc, and its counter-image, are consistent with a second background source in the model. Hammer et al. (1996) refine this model using spectacular HST data. A ^-profile is used for the dark-matter halo, again perturbed by a central cD galaxy. With the higher resolution data, not effected by atmospheric seeing, Hammer et al. reconstruct two sources behind the cluster, consistent with the earlier model of Mellier et al.. The new lens inversion scheme is implemented on the HST data. We have the benefit of the results of Mellier et al. and Hammer et al. to test the validity of the inversion scheme. Following the model of Mellier et al., we propose a PID dark-matter halo, perturbed by a central cD galaxy following a density profile suggested by Miralda-Escude (1995). The geometry of the lens and the deflector mass distribution are fixed by measuring the polar moments of the triply-lensed background source. Without changing the model parameters, an additional source is added. This second source appears twice through the simulated lens, very close to the observed positions of the radial arc and its tangentially distorted counter-arc. The second part of the inversion scheme, in which the background source is reconstructed from the observations, agrees with the work of Hammer et al., both in the relative position of the two sources and in their structure and surface brightness. The final consistency check of viewing the reconstructed source through the lens model quite nicely reproduces the lensed features. In particular, the double ring structure in the giant arc is visible, and no additional spurious arcs are formed. Further sources are subsequently added at various redshifts, again producing features which appear to have optical counterparts. Since the nature of the gravitational lensing in MS 2137 had been previously established, the success of the inversion validates the two-stage inversion scheme at the centre of this thesis. We apply the algorithm to MS 1455, in the belief that our results will yield a well-defined model consistent with the observational data. We model the dark-matter halo of MS 1455 with a NFW profile and separately with a PID profile. Both models contain a central cD galaxy fit to the observations. These models suggest that single background source is responsible for both the tangential arclet identified by LeFevre and the radial arc. The models predict a third tangential arc, which apears to have an optical counterpart in the observations, and Chapter 7: Conclusions 150 is subsequently included in the modeling procedure. The peculiar configuration of the source lying over the intersection of tangential and radial caustics on the source plane greatly restricts the model parameter space. Polar moments of the three lensed arcs sufficiently quantify the position and shape features to admit a close fit between the data and the simulated arcs. Upon reconstructing the source from the data, a collection of spurious pixels appears on the source plane. By cleaning the data, we can exclude these points, and find that the source has a generally elliptical shape with a brighter central component, in both the PID+cD and NFW+cD models. When the reconstructed source is passed back through the models, the third arc is greatly extended. By processing the data from the three arcs separately, we discover that the extended appearance of arc A3 is due to misidentifying lensed light in the observations. The small, bright image of A3 predicted by the simulations is spread over a much wider region of the image, and appears quite similar to the observations, when the simulations are convolved against a point-spread function. Assuming gravitational lensing is entirely responsible for the appearance of arc A3 requires a much larger source than is consistent with the source of the other arcs. The error is magnified when the source is relensed, reproducing the extended source, but also a large amount of additional light. We conclude that the seeing at the instrument plays a significant role in the appearance of arc A3. Deconvolution of ground-based images, or space-based observations, may remove this problem. Alternatively, each arc could be "cleaned" with the LensClean algorithm, using a beam matching the point-spread function. The radial arc formed by passing the reconstructed source back through the lens does not extend towards the centre of the lens as far as the observed arc. We trace this to the data-cleaning that removed the spurious pixels on the source plane. As the pixels are removed, we flag their origin in the image, and find that 100% of the discarded points come from the inner end of the radial arc. This indicates that our model incorrectly describes the lensing behaviour near the core of the cluster. This is not unexpected, as gravitational lensing measures only the cumulative mass about the centre of the lens. Chapter 7: Conclusions 151 The models we choose to produce the results show inconsistencies characteristic of gravita-tional lensing: the highly non-linear response of the appearance of the lens to small changes in the source, and the inability to determine the mass distribution at the centre of the lens, where small perturbations have a significant impact on the cumulative mass. 7.2 Strengths and Weaknesses of the Scheme The lens inversion algorithm is based on forward modeling the gravitational lensing behaviour. This is weaker than a true inversion scheme, and perhaps even the "cleaning" approaches of LensClean. There is always the possibility that another, undiscovered, model might produce better results. We are confident that the models we have presented are close to optimal in this neighbourhood of parameter space, but we can say nothing about the existence of another model that lies far from this neighbourhood. A choice was made at the outset, however, to pursue the forward modeling approach. The weakest component of the scheme implemented here is the lack of an objective, quanti-tative measure of the closeness of fit between the observations and the simulations. The statistic one should employ to make this quantitative measurement, however, is unclear. Pixel-by-pixel comparisons on the image plane are unsuitable because of the effects of seeing in ground-based observations and presence of noise in the data. We might consider finding a suitable statistic on the source plane. We have noted the limitations of tracing bright knots or peaks of light back to the source plane, measuring the distance between the supposedly coincident source points, and using this distance as a quantitative measurement of the closeness of fit: bright knots can be mis-identified because of the great distortion induced in strongly lensing systems. In the case in the triply-imaged source in MS 1455, the source of the bright tangential arc can be modeled as the dim extension of what appears to be a background spiral galaxy, not the bright central bulge. Polar moments are introduced as an alternative measurement of the lensed features. Al-though artificial, the first moments f, 0 and second moments Qrr, QTg, Qgg characterise the Chapter 7: Conclusions 152 position and shape of the lensed arcs remarkably well. The shape parameter x, which mea-sures the ratio of the tangential and radial "lengths" of an arc, distinguishes between radial and tangential arcs. In Section 5.3.1 we remark upon some improvements that could be made, in particular introducing the origin of the polar coordinate system as another model param-eter. We make no attempt to combine the moments into a single quantitative statistic, but independently fit the polar centroid (9, f) and the diagonal quadrupole moments, QTT and Qee-One of the greatest strengths of the routines introduced here is the interactive modeling component. Because of the highly non-linear nature of lensing, it is very difficult to make quantitative predictions about the position and shape of each lensed feature. We can make general statements about how many images should occur, and roughly where the arc/counter-arc pairs (or triplets or more) should lie, but the detailed structure of the arcs, particularly the effect of amplification, can only be discovered through models. One need only "play" with a lens for short time to discover the critical features of the lens. A real-time, interactive program like lens quickly gives an intuitive understanding of the optics of gravitational lenses. We have noted the weaknesses of parametric models in inverse problems like gravitational lensing: small changes in the parameters may result in large changes in the outcome. There are strengths to the parametric models used here, as well. As outlined in Section 4.3, a single mass, single source model is described by 13 parameters. Often, some of these parameters are set by the observations, like the redshift of the deflector plane. While we use 5-parameter elliptical sources to fit the deflector mass distribution and the lens geometry, these parameters are removed from the final model when the source is reconstructed from the observations. From each lensed image we tabulate the centroid (9, f) and three second moments Qrr, QRQ, Qgg. We treat the mixed quadrupole moment Qrg only as a test of the consistency of the simulated arc. Thus four parameters can be extracted from each lensed feature, except in exceptional cases like arc A3 in MS 1455 where seeing dominates the appearance of the arc. Even in these cases, the centroid (9, f) is reliable statistic. Therefore, three images of a single source should be sufficient to constrain the model. Two arcs are virtually guaranteed by the definition of strong gravitational lensing. Asymmetry of the dominant mass distribution may introduce a third arc, Chapter 7: Conclusions 153 as occurs in MS 1455 and with one of the sources in MS 2137. If there is an additional background source (5 parameters), two or three more arcs may form with 8-12 measurable statistics, and the model becomes predictive, an essential feature of astronomical models. Through our choice of physically-motivated model parameters (mass, orientation, redshift, etc.) these predictions can be easily tested, confirmed, or refuted. It is the nature of research that we cannot list every direction this work may point. Two paths immediately come to mind, however. 7.3.1 Lensing as an Inverse Problem Gravitational lensing is a classic example of an inverse problem: a signal is emitted, passes through a black-box, and the outcome is recorded. The goal is to discover the contents of the black-box. We know from experience that subtle changes in the deflector mass distribution and the background source can cause great changes in the lensed outcome. This fairly defines an unstable inverse problem. The inversion of the lens, therefore, is a delicate process. It is very tempting to plunge into an inverse problem approach to lensing based on the appearance of the deflection angle: The deflection angle is simply the convolution of the surface mass density E(£) against a lensing "kernel" at point £, a Fredholm equation of the first kind, which is described by Parker (1994) in the context of geophysics. It may be possible to employ the techniques of geophysical inverse modeling to gravitational lensing to recover the surface mass density. A solution comes in the form of a discretised array, much like the solutions generated by the LensClean routines of Kochanek et al.. 7.3 Outlook (7.1) Inversions of this type generally proceed by minimizing an objective function, subject to reproducing the data to within the limits of the errors. The objective function gives us great Chapter 7: Conclusions 154 freedom to persuade the surface mass distribution to show specific behaviour. By minimizing S 2 we keep the model small; minimizing | V S | 2 favours smoother models. To extract meaningful parameters from the solution, it may be favourable to minimize a functional of SpjD — S, for instance, which finds the discrete distribution most closely matching the PID mass profile. By including model parameters describing the mass distribution in the inversion, the solution can predict the best SIS, or PID, or NFW model. In the geophysical setting, the data which must be reasonably reproduced at all times are the direct result of the convolution (7.1). In the case of lensing, however, we do not observe the deflection angle a. Rather, we see only the light (if any) coming from the point on the source plane r) = (Ds/D,j)^ — -D^fi, the lens equation (3.6). It is as if there is a black-box inside a black-box. Solving for the surface mass distribution £ requires opening both of these boxes, one after the other. The non-linearity and instabilities are multiplied. Perhaps we can proceed by first finding the deflection field a(£), based on the positions of the lensed features. We could follow the philosophy of the Maximum Entropy Method to find the surface brightness distribution on the source plane which can reproduce a collection of arcs. The objective function of the M E M can be chosen to give the smoothest, or most localized source, for instance. Then, having established the deflection field, and also the magnification, the second black-box could be opened to find a surface mass distribution which generates such a deflection field. We note that this approach is very nearly the reverse of the two-stage lens inversion scheme presented in this thesis. If this approach is successful, the solution will have many useful features not available in the forward modeling solution. A quantitative measure of the goodness-of-fit is available, the minimum value of the objective function, and different models can be objectively compared. Analysis of the kernel functions produces quantitative measures of the sensitivity of the solution to each kernel, highlighting the pixels in the array, and the regions in the lens, which play a dominant role in forming the solution. At the same time, the resolution of each kernel can be measured, which in turn tells us how confident we can be about the value of the surface mass density at any given location on the deflector plane. These are all features that would be Chapter 7: Conclusions 155 invaluable to the interpretation of the result. Because of non-linear and unstable nature of the problem, the minimum in the objective function must be slowly and carefully approached. We run the risk of finding a local minimum rather than the global minimum. How the minimum is approached is a problem addressed by the specific implementation of the iterative process. Locating the global minimum of the objective function is a separate question, however. To find the optimal solution, it may be necessary to make an excellent initial guess. Perhaps this is where the forward modeling lens inversion scheme described in this report would be most useful. After producing a model which is 70 or 80 or even 90 percent of the way to optimal, the forward modeling solution could be passed to the inverse problem routines for the final step to the global minimum. 7.3.2 Astrophysical Applications While the author approached this work from the field of applied mathematics, one cannot spend any amount of time working with astronomical data without being intrigued by the astrophysical processes at play. One area which can be explored by the lens inversion schemes presented here is the mass-to-light ratio of the central cD galaxy. Among the many characteristics of a galaxy which are are linked to the mass-to-light ratio is the amount of dark matter associated with the collection of luminous stars making up the galaxy (for example, (Miralda-Escude, 1995; Natarajan and Kneib, 1997; Deh"Antonio and Tyson, 1996)). Adjusting the mass parameter ac£> of the central cD galaxy in MS 1455 is equivalent to adjusting the mass-to-light ratio of cD, for the amount of light collected in the observations remains fixed. Because we have chosen to model the surface mass density of the cD component under a mass-follows-light hypothesis, the shape of the mass distribution and the shape of the light distribution are identical, only an over all scaling constant distinguishing the two. This constant is precisely the mass-to-light ratio, V. By design, V does not vary as a function of radius or position angle, only a first approximation of the true situation. It is a simple matter to find the surface mass density E c o of the central cD galaxy for each Chapter 7: Conclusions 156 value of the mass parameter a and hence for any line-of-sight velocity dispersion o~ios[cD]. The amount of light received at the instrument from this cD galaxy is a much more challenging question, for we must deduce the intrinsic surface brightness of the cD galaxy. The (equivalent to) 4-hour 7-band exposure of MS 1455 at the CFHT gives calibrated 7-band magnitudes for the data. Two 30 minute V-band exposures were also made, unfortunately with much poorer seeing. Analysis of the observations produces (V — I) = 1.46, slightly smaller than the measurement of 1.8 of Small et al. (1994). The observed /-band is shifted into the rest frame V over the distance Zd = 0.257 to MS 1455, and converted to rest frame luminosity. The scale factor between the projected mass distribution and the luminosity is the mass-to-light ratio, Ty-measured in solar masses M Q per solar luminosity L Q . As shown by Miralda-Escude (1995) in spherically symmetric lenses, the presence of a central cD galaxy can inhibit the formation of a radial arc. This is due, again, to the signif-While maintaining the large radii lensing characteristics of MS 1455, namely the position of the tangential arclet A2, we can constrain the range of the parameter aco by the appearance (or disappearance) of the radial arc. Recall the delicate interplay of the radial and tangential caustics on the source plane in our models of MS 1455: a single source must lie on the inter-section of these two curves to produce the three observed images. Altering the mass parameter of the cD significantly effects the location of the radial caustic, while having little effect on the tangential caustic. Of course, these changes can often be countered by changing one or more of the other free parameters at our disposal. For the sake of simplicity, we fix all the parameters on the lens plane except the mass of the cD galaxy, and alter the source plane parameters, the redshift z3 and the source position y. The radial arc forms when even and odd parity images merge across the radial critical line. While it is difficult to determine precisely when the radial arc fails to form, it is clear when there is no radial arc whatsoever. In Figure 7.1, we show a range of source redshifts which permit the formation of a radial arc close to the position of arc A l in MS 1455 for a variety of mass parameters OCD, all the while maintaining the position (7.2) icant portion of the cumulative mass near the centre of the lens coming from the cD galaxy. Chapter 7: Conclusions 157 of the tangential arclet A2. Outside the redshift band at each value of cr, we can confidently conclude that a radial arc fails to form. To each value of acn there corresponds a mass-to-light ratio T via (7.2) which is also shown in the Figure. 1.000 ' V 12.0 0.900 A t 0.800 o S H s o 0.700 A