Inverting Gravitational Lenses

Simulating gravitational lenses with `lens`

This page describes how to simulate gravitational lenses using the matlab routine lens. This routine accompanies the article Inverting Gravitational Lenses by Peter R. Newbury and Raymond J. Spiteri.

The Figure below is a 128-by-128 pixel simulation showing two gravitationally lensed arcs of light. The arcs are highly distorted images of a uniformly bright, elliptical disks of light hidden behind the lens. The position of this background source is outlined by the magenta ellipse in the Figure. The colour of the arcs is a measure of the magnification, with high magnification at the red end of the spectrum and low magnification at the blue end of the spectrum. The deflecting mass is a singular isothermal sphere (SIS) represented by its Einstein Ring (drawn in red) which has a radius of about 6 arcseconds in this example.

This simulation is a parametric model of lensing, in which the model parameters are fed into a "general-purpose" modeling routine. You can generate the simulation above, and other like it, with the matlab routine lens and the parameter file SIS1.

Downloading the software Running the simulations Changing the simulations

Downloading the software

These simulations require the matlab routine lens and a text file of parameters. These parameters can later be changed to generate different simulations.

1. Download the matlab routine lens.m.: Be sure to save this text file as "lens.m" (the default "Save as..." when using Netscape.) There are comments throughout the code lens.m which may be useful for those who are interested in altering the routine.
2. Download the parameter file SIS1.: Save the parameter file as a text file. As described below, you can later use this parameter file as a template for generating different simulations.

Running the simulations

The routine lens is a matlab function, so you must be able to run matlab. The routine was originally written for "Version 5.3.0.14912a (R11) Jun 19 1999". The commands below are entered into the matlab console:

1. Run the routine (be sure to enclose the parameter file in single quotes.)

matlab>> lens('SIS1')

2. The simulation above should appear momentarily on the desktop.

If the simulation does not appear and the matlab console indicates an error has occured, check that the parameter file is valid, following the instructions on altering the parameter file.

The position of the figure window has been coded into the routine, and you may wish to adapt it to better fit your screen. Look for the line

set(gcf,'Position',[300 100 400 400]);

which generates a 400-by-400 pixel figure window starting 300 pixels over and 100 pixels up from the bottom left of the screen.

The console prompts you for further action:

matlab>> lens('SIS1') Gravitational lensing ========================= Options q - quit w - write image to file Select an action:

3a. Quit

Press 'q' to halt the routine. Be sure the matlab console is the active window (on some platforms, the figure window automatically becomes the active window.)

Once the routine is stopped, you can print the figure, manipulate the graphic, etc., as usual.

3b. Write image to file

Press 'w' to save the image to file. A data file SIS1.mat using matlab's internal storage format will be created in the current directory. (In general, the parameter file foo will produce a data file foo.mat.)

We include this action because the forward modeling routine is used to create "observations" which will later be inverted following the inversion scheme described Inverting Gravitational Lenses and implemented with the invert routine.

Changing the simulations

There are basically two ways to change the simulation shown in the Figure above.

The interested student of lensing will alter the values of the parameters to explore the complex and unpredictable nature of lensing, even around this simple SIS mass. We encourage users to spend some time experimenting with the simulations. The details of the model parameters and parameter file are given below.
Hard-core lensing and matlab enthusiasts will alter the lens.m code itself. We offer some interesting and useful modifications below.

Changing the simulations - Altering the parameters

The whole idea behind a parametric model is that you can produce different simulations without having to re-write the program. You simply change the values of a collection of model parameters. Here we use the parameter file SIS1 as a template.

The parameter file is a text file with one parameter given on each line, in the form parametername=parametervalue

Lines which begin with % are comments and are ignored by the routine.

Parameter File 'SIS1'

Line/Parameter Description

% Parameter file 'SIS1' Comments

% SIS with elliptical background source

%

% observation parameters The parameters come in three sets, one describing the nature of the "observation", one describing the lensing mass, and one describing the background source.

rho=0.2 Each pixel is rho arcseconds in width (and height)

N=128 The simulation is N-by-N pixels (so that the simulation covers a N*rho-by-N*rho arcsecond window on the sky.)

x10=0.0
x20=0.0 Coordinates, in arcseconds, of the center of the field-of-view. In this example, the "telescope" is "pointed" at the center of the SIS mass, but you can "look" at other parts of the simulation, too. This can be useful for zooming in on interesting features by changing (x10,x20) and increasing the resolution by decreasing rho.

% lens parameters The next four parameters describe the SIS mass distribution.

x1SIS=0.0 x2SIS=0.0 Coordinates, in arcseconds, of the center of the SIS mass distribution.

sigma=600 Mass parameter of the SIS mass distribution, typically in the range 100-1000. See the article for a description of sigma.

zd=0.3 Redshift of the deflector plane. Most of the lenses which have been discovered lie at redshift z_d <1. This redshift must be less than the redshift of the source plane z_s given below.
See the CfA-Arizona Space Telescope Lens (CASTLE) Survey for a list of observed lenses and their properties.

% source parameters The next seven parameters describe the elliptical background source of light which is observed through the gravitational lens. In the modeling routine lens, we have chosen to "observe" very simple sources. In reality, the background source is often a distant galaxy, which introduces new structures in the lensed images.

zs=0.7 Redshift of the source plane. This must be greater than the redshift z_d of the deflector plane.

y1SRC=1.0 y2SRC=1.5 Coordinates, in arcseconds, of the center of the background source.

radiusSRC=1.5 Semi-major axis of the ellipse.

ellipticitySRC=0.4 Ellipticity 1- b/a of the source ellipse, a variation on the eccentricity of the ellipse.

orientationSRC=40 Orientation of the ellipse: the semi-major axis is rotated this number of degrees counter-clockwise from the positive x-axis.

fluxSRC=1.0 The sources in these simulations are assumed to be uniformly bright, with this signal. With fluxSRC=1.0, the colorbar and shading of the lensed images are direct measures of the magnification factor.

Parameter File '`SIS1`'
Line/Parameter	Description
`% Parameter file 'SIS1'`	Comments
`% SIS with elliptical background source`
`%`
`% observation parameters`	The parameters come in three sets, one describing the nature of the "observation", one describing the lensing mass, and one describing the background source.
`rho=0.2`	Each pixel is `rho` arcseconds in width (and height)
`N=128`	The simulation is `N`-by-`N` pixels (so that the simulation covers a `Nrho`-by-`Nrho` arcsecond window on the sky.)
`x10=0.0` `x20=0.0`	Coordinates, in arcseconds, of the center of the field-of-view. In this example, the "telescope" is "pointed" at the center of the SIS mass, but you can "look" at other parts of the simulation, too. This can be useful for zooming in on interesting features by changing `(x10,x20)` and increasing the resolution by decreasing `rho`.
`% lens parameters`	The next four parameters describe the SIS mass distribution.
`x1SIS=0.0 x2SIS=0.0`	Coordinates, in arcseconds, of the center of the SIS mass distribution.
`sigma=600`	Mass parameter of the SIS mass distribution, typically in the range 100-1000. See the article for a description of `sigma`.
`zd=0.3`	Redshift of the deflector plane. Most of the lenses which have been discovered lie at redshift z_d <1. This redshift must be less than the redshift of the source plane z_s given below. See the CfA-Arizona Space Telescope Lens (CASTLE) Survey for a list of observed lenses and their properties.
`% source parameters`	The next seven parameters describe the elliptical background source of light which is observed through the gravitational lens. In the modeling routine `lens`, we have chosen to "observe" very simple sources. In reality, the background source is often a distant galaxy, which introduces new structures in the lensed images.
`zs=0.7`	Redshift of the source plane. This must be greater than the redshift z_d of the deflector plane.
`y1SRC=1.0 y2SRC=1.5`	Coordinates, in arcseconds, of the center of the background source.
`radiusSRC=1.5`	Semi-major axis of the ellipse.
`ellipticitySRC=0.4`	Ellipticity 1- b/a of the source ellipse, a variation on the eccentricity of the ellipse.
`orientationSRC=40`	Orientation of the ellipse: the semi-major axis is rotated this number of degrees counter-clockwise from the positive x-axis.
`fluxSRC=1.0`	The sources in these simulations are assumed to be uniformly bright, with this signal. With `fluxSRC=1.0`, the colorbar and shading of the lensed images are direct measures of the magnification factor.

Notes about the simulations

It is really only the displacement or off-set between the center of the SIS lens and the center of the elliptical source that matters. All you really observe, after all, is the relative position of the foreground SIS and background source.
Perfect source-lens-observer alignment, when (x1SIS,x2SIS)=(y1SRC,y2SRC) generates an Einstein Ring. If the background source is not simply a point-source, but has some size like the ellipses used here, an Einstein Ring will form if any part of the background source lies directly behind the center of the SIS. This parameter file is almost identical to SIS1 except the background ellipse is moved to create an Einstein Ring. Notice the extreme magnification that occurs, compared to the SIS1 example, when an Einstein Ring forms.
The deflection angle depends only on the ratio sigma² (z_s - z_d)/ z_s so these three parameters are not independent. A distant, weak lens and a nearby, strong lens can produce the same set of observations.
One arc or two? In the special case of spherically symmetric SIS lenses explored here, the lens equation reduces to a one-dimensional (radial) equation y = x - R_E x / | x | where R_E is the Einstein Radius of the lens.
If y < R_E (that is, the ellipse lies within the Einstein Ring, like SIS1 in the Figure above) then there are 2 solutions to the lens equation
x₁ = y + R_E > R_E
x₂ = y - R_E < 0 and two arcs are formed, one outside the Einstein Ring at x₁ and one inside the Einstein Ring at x₂ on the other side of the center of the lens.
However, if y > R_E (that is, the ellipse lies outside the Einstein Ring) then there is only one solution to the lens equation
x₁ = y + R_E and only one arc is formed, the same arc that formed outside the Einstein Ring in the earlier case. It is very interesting to look at the family of models in which the source moves from inside the Einstein Ring to outside the Einstein Ring, because you can watch the second (interior) arc disappear. And running the family in the opposite direction, you can watch the interior arc form as the background source begins to pass behind the Einstein Ring. Animating this beautiful behaviour is one of the modifications discussed below.

Changing the simulations - Changing the routine

Those familiar with matlab's capabilites will soon feel the need for a modified version of the lens routine. Lensing is a non-linear, ill-posed problem, and small changes in the model parameters can cause disproportionately large changes in the results. To see these large changes, it is very helpful to see the before and after observations next to one other. Or better yet, to animate the change.

The simplest way to turn this one-time simulation routine into an interactive routine is to add more cases to the routine. In the form supplied here, the routine lens has only two actions, quit and write-to-file, but more can be added. For example, you can use the arrows on the numpad (2=down, 6=right, 8=up, 4=left) to move the background source by one pixel:

... switch lower(button) case '4' y1SRC=y1SRC-rho ...

Then you'll have to re-draw the figure: cut and paste all the graphics command and insert them inside the while button~='q' loop. With double-buffering turned on, the figure will respond to your keystrokes seamlessly.

Other possible changes include

'm' and 'M' keystrokes to decrease or increase the mass of the SIS
'z' and 'Z' keystrokes to decrease or increase the redshift of the source
If you want to get carried away, think about using ginput(1) to read the coordinates of the pointer, and have the background source move in the direction of your mouse click.
Create an animation by using matlab animation routines, upgrading this old MPEG simulation.
Micro-lensing simulation: In a sequence of simulations, pass the source from one side of the lens to the other. Arcs will form and the magnification will peak. Then the magnification will drop again as the source passes out the other side of the lens. Animate this sequence and watch the "movie." The rapid rise and then fall in magnification is the signature of gravitational lensing that micro-lensing studies search for. The mass of the intermediate lens can be extracted by measuring the height of the magnification peak and the rate at which the magnification rises and falls.