Hands-on Workshop:
How to calculate Quasar Microlensing
Joachim Wambsganss
Zentrum für Astronomie der Universität Heidelberg (ZAH/ARI)
Journal of Computational and Applied Mathematics 109 (1999) 353–372
www.elsevier.nl/locate/cam
Gravitational lensing: numerical simulations with a hierarchical
tree code
Joachim Wambsganss
Astrophysikalisches Institut Potsdam, An der Sternwarte 16, 14482 Potsdam, Germany
Received 12 May 1998; received in revised form 6 January 1999
Abstract
The mathematical formulation of gravitational lensing — the lens equation — is a very simple mapping R2 → R2 ,
between the lens (or sky) plane and the source plane. This approximation assumes that all the de!ecting matter is in one
plane. In this case the de!ection angle ! is just the sum over all mass elements in the lens plane. For certain problems
— like the determination of the magni"cation of sources over a large number of source positions (up to 1010 ) for very
many lenses (up to 106 ) — straightforward techniques for the determination of the de!ection angle are far too slow.
We implemented an algorithm that includes a two-dimensional tree-code plus a multipole expansion in order to make
such microlensing simulations “inexpensive”. Subsequently we modi"ed this algorithm such that it could be applied to a
three-dimensional mass distribution that "lls the universe (approximated by many lens planes), in order to determine the
imaging properties of cosmological lens simulations. Here we describe the techniques and the numerical methods, and we
c 1999 Elsevier Science B.V. All rights reserved.
mention a few astrophysical results obtained with these methods. "
September
18, 2012; XI-thGravitational
School of Cosmology,
IESC
Cargese;
Keywords:
Astrophysics;
lensing;
Tree
code Joachim Wambsganss: “Workshop: How to calculate Quasar Microlensing”
Efficient
Fig. 2.criterion
Example of is
the something
use of the hierarchical
tree code
in a microlensing
scenario
ran
like the
opening
angle of
the with
cell31 as
(a, left) Cell structure of the point lenses; (b, middle) Illustration of the cell structure as “t
this Ray
ratio
is smaller
than
a chosen
value
of thedot).“accuracy pa
Inverse
A
Tree-Code
Approach
con!guration
to beShooting:
used
for a certain
ray
position
(marked
with a circled
angles of all individual stars acting as
just
a the
summation
of int
th
(Wambsganss
1990,is
1999)
one lens,lenses
this cellhence
is used. is
Otherwise,
cell
resolved
lengths
again are
with
criterion
is something
likecompared
the opening
anglethe
of distances
the celln as of
seentheir
fromcenter-o
the pos
n
!and
and issosmaller
on. Typically,
!!
ranges
0.4
0.9.ij If! aorparti
this ratio
than a chosen
value ofbetween
the “accuracy
parameter”
if th
Deflection angle
lenses:
oneforlens,
this de"ection
cell is used. angle,
cell is resolved
into its
iOtherwise,
j (up
ofn the
it isthejiconsidered
w
2 to four)
2 as a pseudo-lens
lengths
again are
withofthe
distances
of their center-of-mass
t
ij positions
located
at compared
the center
mass
determined
by all these
particle
j=1
j=1
and so on. Typically, ! ranges between 0.4 and 0.9. If a particular cell is used f
of the de"ection
angle, it is considered as a pseudo-lens with the total mass
Number of computational
operations:
2.3.3.
jpseudo-lenses
located
at theLenses
center ofand
mass
determined by all these particles.
4G
!˜ =
c
!˜ =
r
M :
r
Here M is the mass of point lens j, r
In thelight
approach
justi described,
we use
the j,
hierarchical
tree
ray
and
point
lens
and
r
is
ij
2.3.3.ofLenses
and pseudo-lenses
de"ection
!˜N*by
two split
parts
according
to the directly include
Calculation of deflection
angle
for
lenses
into
two
parts:
In the approach
we useof
the hierarchical
tree code
in; order
toisa
y
)
is Njust
thedescribed,
position
ray
i,
and
(x
j
j
Nto the directly included lenses !˜ and th
Naccording!
of de"ection !˜ by!
two parts!
far
most
the!˜L computing
time
!˜ = By
!˜ ≈
!˜ + of
!˜k =:
+ !˜C :
!
! i
! j
!˜ =
!˜ ≈i=1 !˜ + j=1!˜ =: !˜ k=1
+ !˜ :
of these de!ection
angles. As seen abo
The N ’s denote the following:
The N ’s denoteone
the following:
light ray takes about ten mathem
number
of all lenses,
∗ is
athelarge
number
of pixels Npix is requir
• N •is N
the
number
of all lenses,
the number
to be
included directly,
of lenses of
to lenses
be included
directly,
• N •theNLnumber
a (=
high
density
rays
per
• N •the
of cells like
(=
to
be included.of
NCnumber
thewould
number
of pseudo-lenses)
cells
pseudo-lenses)
to be
included.
number
ofupindividua
" N , the calculation
is speeded
considerably.
For cases in whichThe
N + N(minimal)
N∗
i=1
j
j=1
L
NC
NL
i
C
L
∗
k
L
C
k=1
∗
L
C
For cases in which NL + NC " N∗ , the calculation is speeded
C
∗
September 18, 2012; XI-th School of Cosmology, IESCLCargese; Joachim
Wambsganss:
“Workshop: How to calculate Quasar Microlensing”
and so
! ranges
between 0.4
If a p
for each ray.
In on.
theTypically,
general case,
a smoothed
outand
or 0.9.
continuo
of to
thethe
de"ection
angle,
is considered
as a pseudo-lens
contributes
de"ection
as itwell
(as an additional
constan
Efficient Inverse
Ray Shooting:
A Tree-Code
Approach
mass determined by all these part
equation located
then is at the center of (Wambsganss
1990, 1999)
Lens Equation:
!
"
N∗
#
mi (x − xi )
2.3.3.
Lenses
and
pseudo-lenses
1−"
0
y = In the approach justxdescribed,
− !c x − we use the hierarchical
;
t
2
0
1+"
(x
−
x
)
i
i=1
of de"ection !˜ by two parts according
to the directly incl
where the macromodel
of !
the
lens has to specify
NC
N∗
NL gravitational
!
!
Tree code approach:
density !c (cf.!˜ [26]);
=
!˜x
˜j + positions
!˜k =: !˜Lof+ the
!˜C : lenses. Tha
i are !the
i ≈
j=1 withk=1
are mapped onto ai=1rectangle
a side ratio (1 − ! − ")=(
360
J. Wambsganss / Journal of Computational and Applied Mathematics 109 (1999) 353–372
The Nmass
’s denote
the following:
is the surface
density
in compact objects). As we wan
choose the shooting #eld, i.e., the area in the lens plane in w
• N is the number of all lenses,
with the side∗ ratio (1 − ! + ")=(1 − ! − ").
• NL the number of lenses to be included directly,
We produced
maps
in the
plane with
• NC thesuch
number
of cells
(= source
pseudo-lenses)
to be very
includd
purpose and the desired resolution: the smallest receiving sq
NC " N∗ ,#0the
spee
For cases of
in which
NL 0+ (where
had a sidelength
L = 0:8#
= calculation
$E × DS isis the
the light
shown
in Fig.
cells with more th
the source For
plane).
Thisraywas
useful
for2c,
thesixhigh-resolution
s
lenses),
whereas all
other lenses
areLtreated
individual
[36]. The16largest
sidelength
considered
was
= 2000#
0 in an
= 6. This is not very remarkable yet. But in a calc
objects ofNCplanetary
mass [24]. In practice, we have to use
" = 1:17; # = 0:83; L = 150$0 ; ! = 0:6), the average number o
than L=(1 tree
−lenses
!code±in")aused
because
due
to 31the
grainyness of the mat
Fig. 2. Example of the use of the hierarchicalof
microlensing
scenario
with40.
randomly distributed lenses.
directly
was
(a, left) Cell structure of the the
point receiving
lenses; (b, middle)
Illustrationand
of theothers
cell structure
as “tree”;
(c, right)
Cell/lens scattere
from
farther
out
Thesquare,
computing
time for this
hierarchical
treeare
method (in
con!guration to be used for a certain ray position (marked with a circled dot).
whereas a direct summation would increase as O(N∗ ). In
September 18, 2012; XI-th School of Cosmology, IESC Cargese; Joachim Wambsganss: “Workshop: How to calculate Quasar Microlensing”
Efficient Inverse Ray Shooting: A Tree-Code Approach
(Wambsganss 1990, 1999)
J. Wambsganss / Journal of Computational and Applied Mathematics 109 (1999) 353–372
363
Fig. 4. Illustration of the way we “shoot” light rays: For each (test-)ray marked with a cross the cell/lens con!guration
is determined. For n2!x (test-)rays (marked with circles) the de"ection angle due to this !xed cells is determined (here
n!x = 4; in the calculations we usually took n!x = 10). On the lowest level n2int rays (marked with dots) are “really” shot,
with the de"ection angle of the cells determined by interpolation and that of the NL close lenses calculated directly (here
September 18, 2012;
School
Cosmology,
nint = 4; in the calculations
we XI-th
took
nintof =
10). IESC Cargese; Joachim Wambsganss: “Workshop: How to calculate Quasar Microlensing”
71
Efficient
Inverse
Ray Shooting: A Tree-Code Approach
362
J. Wambsganss / Journal of Computational and Applied Mathematics 109 (1999) 353–372
(Wambsganss 1990, 1999)
Fig. 3. The accuracy increase of the de!ection angle obtained with the inclusion of higher-order multipole moments (see
September 18, 2012; XI-th School of Cosmology, IESC Cargese; Joachim Wambsganss: “Workshop: How to calculate Quasar Microlensing”
Section 2.3.4) is shown here. In each panel the di"erence between the directly determined de!ection angle (!˜direct ) and
Efficient Inverse Ray Shooting: A Tree-Code Approach
(Wambsganss
1999)
J. Wambsganss / Journal of Computational and Applied Mathematics
109 (1999)1990,
353–372
365
Fig. 6. Left: Magni!cation pattern for surface mass density of ! = 0:32 and external shear of " = 0:18 (parameters for
image A of the double quasar Q0957+561). The sidelength corresponds to 40 Einstein radii. The “lenses” here are objects
of planetary mass (mLens = 10−5 M" ), and the white bar indicates location and length of the track along which the light
curve (displayed in the right panel) is evaluated. The length of the light curve (for a quasar with a Gaussian brightness
pro!le of size 3 × 1014 cm) corresponds to 160 days (for an assumed transverse velocity of 600 km/s). More details on
this project can be found in [24].
September 18, 2012; XI-th School of Cosmology, IESC Cargese; Joachim Wambsganss: “Workshop: How to calculate Quasar Microlensing”
Efficient Inverse Ray Shooting: A Tree-Code Approach
364
J. Wambsganss / Journal of Computational and Applied Mathematics
109 (1999) 353–372
(Wambsganss
1990, 1999)
Fig. 5. “Wrong” magni!cation pattern for ! = 0:2 (left panel): here the accuracy parameter is too large: " = 1:2 and the
shooting region is too small. Compare this frame with the correct one (right panel). A magni!cation pattern indicates
September
2012; XI-th School
of Cosmology,
IESCsource
Cargese; Joachim
“Workshop:
How
to calculate
Quasar Microlensing”
the magni!cation
as a18,function
of position
in the
plane. Wambsganss:
Darker gray
means
higher
magni!cation.
On top of the
74
Quasar Microlensing: How to do simulations!
1) copy file
Wambsganss-MicrolensingCode-Cargese-2012.tar
to your disk
2) untar this file ... should produce directory:
Wambsganss-MicrolensingCode-Cargese-2012
September 18, 2012; XI-th School of Cosmology, IESC Cargese; Joachim Wambsganss: “Workshop: How to calculate Quasar Microlensing”
Quasar Microlensing: How to do simulations!
1) copy file
Wambsganss-MicrolensingCode-Cargese-2012.tar
to your disk
2) untar this file ... should produce directory:
Wambsganss-MicrolensingCode-Cargese-2012
3) cd cfitsio
4) ./configure
5) make
(still in directory cfitsio)
6) ..
(now in directory Wambsganss-MicrolensingCode-Cargese-2012)
7) make
(should produce executable “microlens”
8) run the program by typing:
)
./microlens
September 18, 2012; XI-th School of Cosmology, IESC Cargese; Joachim Wambsganss: “Workshop: How to calculate Quasar Microlensing”
Quasar Microlensing: How to do simulations!
1) copy file
Wambsganss-MicrolensingCode-Cargese-2012.tar
to your disk
2) untar this file ... should produce directory:
Wambsganss-MicrolensingCode-Cargese-2012
3) cd cfitsio
4) ./configure
5) make
(still in directory cfitsio)
6) ..
(now in directory Wambsganss-MicrolensingCode-Cargese-2012)
7) make
(should produce executable “microlens”
8) run the program by typing:
9)
)
./microlens
newly produced files:
dat.401
log-file
IRIS401
magnification pattern (unformatted)
IRIS401.fits
magnification pattern (FITS format)
10) display magnification pattern with IDL:
./rnew dis_1000
September 18, 2012; XI-th School of Cosmology, IESC Cargese; Joachim Wambsganss: “Workshop: How to calculate Quasar Microlensing”
IRIS401
September 18, 2012; XI-th School of Cosmology, IESC Cargese; Joachim Wambsganss: “Workshop: How to calculate Quasar Microlensing”
IRIS401
dat.401
September 18, 2012; XI-th School of Cosmology, IESC Cargese; Joachim Wambsganss: “Workshop: How to calculate Quasar Microlensing”
Quasar Microlensing: How to do simulations!
11) in order to extract a lightcurve: compile
lightcurve.f
(I use: gfortran lightcurve -o lightcurve)
12) run lightcurve routine:
./lightcurve
13) output produced:
out_line
(lightcurve data, pixels convolved with source profile)
IRIS401-track (magnification pattern WITH track marked)
14) display magnification pattern with track AND lightcurve:
dis_light
September 18, 2012; XI-th School of Cosmology, IESC Cargese; Joachim Wambsganss: “Workshop: How to calculate Quasar Microlensing”
14) display magnification pattern with track AND lightcurve:
dis_light
Source size: sigma = 3
September 18, 2012; XI-th School of Cosmology, IESC Cargese; Joachim Wambsganss: “Workshop: How to calculate Quasar Microlensing”
14) display magnification pattern with track AND lightcurve:
dis_light
Source size: sigma = 1
Source size: sigma = 3
September 18, 2012; XI-th School of Cosmology, IESC Cargese; Joachim Wambsganss: “Workshop: How to calculate Quasar Microlensing”
14) display magnification pattern with track AND lightcurve:
dis_light
Source size: sigma = 10
Source size: sigma = 3
September 18, 2012; XI-th School of Cosmology, IESC Cargese; Joachim Wambsganss: “Workshop: How to calculate Quasar Microlensing”
Quasar Microlensing: How to do simulations!
11) in order to extract a lightcurve: compile
lightcurve.f
(I use: gfortran lightcurve -o lightcurve)
12) run lightcurve routine:
./lightcurve
13) output produced:
out_line
(lightcurve data, pixels convolved with source profile)
IRIS401-track (magnification pattern WITH track marked)
14) display magnification pattern with track AND lightcurve:
dis_light
15) modify input file for microlens:
September 18, 2012; XI-th School of Cosmology, IESC Cargese; Joachim Wambsganss: “Workshop: How to calculate Quasar Microlensing”
Quasar Microlensing: How to do simulations!
surface mass density kappas (or sigma)
15) modify input file for microlens:
replace nray = “20” by “100” ... and run again!
September 18, 2012; XI-th School of Cosmology, IESC Cargese; Joachim Wambsganss: “Workshop: How to calculate Quasar Microlensing”
Quasar Microlensing: Now YOU do simulations!
Deal:
1) You can use the code “microlens” freely
2) On first scientific paper using “microlens”: J.W offered co-authorship
3) This (and subsequent) papers cite:
Wambsganss, J.: 1999, Journ. Comp. Appl. Math. 109, 353
Wambsganss, J.: 1990, PhD Thesis, Ludwig-Maximilians-University
Munich (also available as MPA report 550)
September 18, 2012; XI-th School of Cosmology, IESC Cargese; Joachim Wambsganss: “Workshop: How to calculate Quasar Microlensing”