Äîêóìåíò âçÿò èç êýøà ïîèñêîâîé ìàøèíû. Àäðåñ îðèãèíàëüíîãî äîêóìåíòà : http://www.atnf.csiro.au/people/Matthew.Whiting/thesis/honoursThesis.ps.gz Äàòà èçìåíåíèÿ: Wed Mar 29 07:17:03 2006 Äàòà èíäåêñèðîâàíèÿ: Sun Dec 23 11:33:00 2007 Êîäèðîâêà: Ïîèñêîâûå ñëîâà: m 31

MONTE CARLO SIMULATION OF ELECTROMAGNETIC
CASCADES IN THE EXTRAGALACTIC INFRARED
BACKGROUND
Honours Thesis
department of physics and mathematical physics
university of adelaide
By
Matthew Thomas Whiting
Supervisor: Dr. R.J. Protheroe
November 1996

Abstract
The detection by EGRET of 40 active galactic nuclei in high energy gamma rays ahs
sparked interest in high energy extragalactic gamma ray astronomy. However, the presence
of the extragalctic infrared background limits the visibility of distant TeV gamma ray
sources. In this thesis I investigate the infrared background, the interactions by which it
absorbs gamma rays, and I develop a program that uses monte carlo techniques, as well as
suitable approximations to model the e#ect of these interactions on gamma ray spectra,
and discuss the implications for high energy gamma ray astronomy.
ii

Acknowledgements
I would like to thank first of all my supervisor, Dr. Ray Protheroe, for his invaluable
help and patience all through this year. I also want to thank the various members of the
astrophysics group who have given me assistance and support over the year, particularly
Troy and Ken. I also must thank the rest of the '96 honours class, for all their support,
friendship and general entertainment over the year. Finally, I want to thank my family
for putting up with me for the past year.
iii

Contents
Abstract ii
Acknowledgements iii
1 Introduction 1
2 High Energy Gamma Rays 3
3 Infrared Background 6
3.1 Intensity at 60µm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
3.1.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
3.1.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
3.2 Observational results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
3.3 Theoretical Predictions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
3.4 Microwave Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
4 Interactions 15
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
4.2 Cross sections . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
4.2.1 Photon­photon pair production cross sections . . . . . . . . . . . . 15
4.2.2 Inverse compton scattering cross sections . . . . . . . . . . . . . . . 17
4.3 Interaction lengths . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
4.3.1 Interaction lengths for pair production . . . . . . . . . . . . . . . . 19
4.3.2 Interaction lengths for inverse compton scattering . . . . . . . . . . 20
4.4 Geometry of Interactions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
4.4.1 Pair production . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
4.4.2 Inverse compton scattering . . . . . . . . . . . . . . . . . . . . . . . 25
5 Monte Carlo methods 28
5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
5.2 Interaction length . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
5.3 Soft photon energy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
5.4 Interaction angle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
5.5 Scattering Angle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
5.5.1 Pair production . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
5.5.2 Inverse compton scattering . . . . . . . . . . . . . . . . . . . . . . . 33
5.6 Azimuthal angle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
iv

6 Electromagnetic cascades 37
6.1 Program structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
6.2 Mono­energetic cascades . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
6.3 Propagation of a power law spectrum . . . . . . . . . . . . . . . . . . . . . 42
6.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
6.4.1 Cascading in the microwave background . . . . . . . . . . . . . . . 43
6.4.2 Cascading in the infrared background . . . . . . . . . . . . . . . . . 44
7 Conclusions 46
Bibliography 47
A Cosmological parameters 49
B Background intensity program 51
C Cascade program 56
D Mean free path and mean photon energy program 69
E Median scattering angle program 74
v

Contents
1

Chapter 1
Introduction
The discipline of high energy gamma ray astronomy has received a major boost in recent
years following the success of the Compton Gamma Ray Observatory satellite, which dis­
covered gamma ray emission from over 40 active galactic nuclei (AGN). However, it has
long been realised that the gamma ray absorption by the background radiation in infrared
and microwave frequencies would limit the visibility of the universe at these gamma ray
energies, and this has been demonstrated to some extent by the lack of detection at higher
energies of all but one of the AGN detected by the CGRO. In this thesis, I investigate
the processes that lead to attenuation of gamma ray emission in both the cosmic mi­
crowave background (CMB) and the infrared background (IRB), and develop a method of
calculating the expected spectra observed at Earth.
The infrared background is important in determining the visibility of active galaxies
at TeV energies -- energies that are able to be observed by ground­based telescopes -- since
gamma rays of these energies will be absorbed to some extent by the IRB. The degree
of absorption depends on the distance to the source and the intensity of the IRB. There
have been a number of predictions in recent years of the form the IRB spectrum takes,
and in Chapter 3 I will look at several of these, compare them with recent observations
and investigate how one of them a#ects gamma ray propagation.
To model the e#ect of a radiation field on gamma ray emission, one needs to account
for at least two processes. The main absorption process is photon­ photon pair production,
where a gamma ray interacts with a microwave or infrared photon to produce an electron­
positron pair. Gamma rays can be created by the process of inverse compton scattering,
whereby a high energy electron scatters a low energy photon up to gamma ray energies.
These two interactions can combine and feed each other to produce a cascade of electrons,
positrons and photons. The theory of these interactions will be discussed in Chapter 4,
including how the energies of the final particles are calculated.
These processess have an inherent random nature, and in order to model them using
a computer program, one generally needs to sample these random parameters from the
relevant distribution. In Chapter 5, I will detail the methods I used to calculate these
parameters, using both the correct sampling techniques for the required distributions, as
well as using approximations for this sampling in the form of the mean or median value
of the distribution.
In Chapter 6, I will show how a complete program may be constructed to model these
cascades, using the routines mentioned above. I then show how this program may be
used to simulate cascades through either the CMB or the IRB over a given distance, and
calculate the resulting spectrum of gamma rays that would be observed on Earth.
The results generated by this program show that pair production in the IRB has
1

CHAPTER 1. INTRODUCTION 2
important implications for the study of high energy gamma ray emission. In particular,
the IRB places an upper limit of about 30 TeV on the energies that will be observed
from extragalctic sources, in particular the AGN that were observed by EGRET at GeV
energies. This limit applies to the closest EGRET AGN (named Markarian 421), and will
decrease for AGN at greater distances. This is in accord with the lack of detection of all
EGRET AGN except for Markarian 421 at energies greater than the limiting energy of
EGRET.

Chapter 2
High Energy Gamma Rays
The study of high energy (E > 100 MeV) gamma ray emission has risen in importance in
the past five years, since the launch of the Compton Gamma Ray Observatory (CGRO)
in 1991, along with the development of sensitive atmospheric Cerenkov detectors such as
the Whipple observatory in the USA.
Our knowledge of high energy gamma ray sources has been vastly increased by the
Energetic Gamma Ray Experiment Telescope (EGRET) on board the CGRO. EGRET is
a multi­level spark chamber that detects gamma rays by converting them into an electron­
positron pair and following the paths of these particles through the spark chamber. It has
a number of levels to determine the direction the gamma ray came from, and this allows
it to image strong sources to within 5 arc minutes (1 arc minute = 1/60 degree). This
means it is able to pinpoint fairly accurately any gamma ray emission and allow it to be
identified with known objects.
EGRET has detected more than 40 active galactic nuclei (AGN) in high energy gamma
rays, ranging in energy from 30 MeV to EGRET's limit of 30 GeV (Thompson et al. [1]).
The locations of these are shown in Fig. 2.1, along with a small number of Galactic gamma
ray pulsars, as well as many sources that are yet to be identified with any other known
objects. Also shown is the location of the Large Magellanic Cloud (LMC). Note that
most of the unidentified sources are situated in or near the Galactic plane, indicating that
they are quite likely in the Galaxy. The sources identified with AGN are located fairly
uniformly about the sky, as you would expect for extragalactic sources.
AGN is a broad name for a class of extragalactic sources that are characterised by
bright, variable central nuclei. AGN are observed in most regions of the spectrum, notably
in radio frequencies, where many AGN are particularly bright radio galaxies. A key
feature of many AGN are large jets of relativistic matter that shoot out from either side
of the galaxy. These jets originate near the nucleus, which is generally thought to be a
supermassive black hole accreting matter from the galaxy. These jets are mostly seen in
radio frequencies, where they are very visible due to synchrotron emission from electrons
caught up in the magnetic fields carried by the jet.
Most AGN associated with gamma ray emission are a type of AGN termed blazars.
These are sources characterised by a strongly polarised and often variable optical emission
and a flat radio spectrum. This means that the radio spectrum is represented by a power
law: S(#) # # -# , where # is the spectral index and has a value # # 0.5 over a given range
of frequencies (usually 100 MHz to 30 GHz). In blazars, the jets are aligned very close to
the line of sight, and it is thought that the gamma ray emission originates in or near these
jets, due to the relativistic beaming of particles caused by the jet motion. A summary of
the properties of gamma ray AGN is given by Dermer and Schlickheiser [3].
3

CHAPTER 2. HIGH ENERGY GAMMA RAYS 4
Figure 2.1: EGRET's observations, showing locations of identified AGN, as well as gamma ray
pulsars in the Galaxy, and a number of sources yet to be identified. The map is in Galactic
coordinates, with the Galactic centre in the middle of the map, and the plane of the Galaxy
extending along the centre horizontal line.
Most of the AGN observed by EGRET have spectra that are well represented by a power
law, where the number of photons per unit energy is given by dN/dE # E -# , where #, the
spectral index, is typically between 1.4 and 3.0, with most from 2.0 to 2.4 (von Montigny
et al. [2]). Note that for gamma rays, the spectrum is represented by dN/dE, the photon
count per unit of energy, since the flux of gamma rays is small enough and the energy of
the individual photons is high enough to talk about the actual photon count. In the radio
part of the spectrum, however, the energies of the individual photons are far smaller, so it
is more meaningful to talk about the emission in terms of electromagnetic waves, leading
to the use of the intensity S(#) in describing the spectrum.
The spectra of blazars are relatively flat (i.e. not decreasing too rapidly with energy),
indicating that emission could extend well into TeV energies and higher (1 TeV = 10 12 eV).
In fact, one of the EGRET sources, the BL Lac object Markarian (Mrk) 421, has been
detected at TeV energies by both the Whipple observatory (Punch et al. [4]) and the
HEGRA experiment in the Canary Islands (Petry et al. [5]), as has another BL Lac, Mrk
501 (Quinn et al. [6] ­ Whipple observatory), although this was not detected at GeV
energies by EGRET.
The Whipple observatory is located in Arizona, USA. It uses atmospheric Cerenkov

CHAPTER 2. HIGH ENERGY GAMMA RAYS 5
imaging to detect gamma rays with energies from 100 GeV to 10 TeV. When a gamma
ray enters the atmosphere, it creates a cascade of photons and particles (such as elec­
trons and positrons) by pair production and bremmstrahlung. The particles are generally
highly relativistic, and move at a speed greater than the local velocity of light. This
causes the emission of Cerenkov light, resulting in a flash of optical light that can be seen
by ground­based telescopes. The Whipple observatory is one such telescope, although
there are Cerenkov observatories located around the world, including the CANGAROO
collaboration at Woomera that the University of Adelaide is involved in.
Mrk 421 remains the only EGRET AGN to have been detected at TeV energies. The
reason for this lies in the fact that gamma rays travelling through a radiation field will
interact with it and get absorbed. The radiation field relevant for TeV gamma rays is
the infrared background (IRB). The average distance gamma rays of these energies travel
through the IRB before interacting (known as the mean interaction length) is of the order
of 100 Mpc (1 Mpc = 3.26 million light years). This is about the distance to Mrk 421,
which is the closest of the EGRET sources. Since the other AGN are much further away,
more absorption takes place, making them less likely to be seen at TeV energies.
At higher gamma ray energies (greater than about 0.1 PeV or 10 14 eV), the cosmic
microwave background (CMB) becomes the dominant source of absorption. The number
density of photons in the CMB is much greater than that of the IRB, so an interaction is
much more likely to occur. This means that the mean interaction length in the CMB is
much shorter (by up to several orders of magnitude) than in the IRB. (See Chapter 4 for
details on the mean interaction length.)
Thus, at PeV energies, the distances over which absorption becomes important are
much shorter. This has the e#ect of ruling out the possibility of detecting any PeV
emission from even the closest AGN (such as Mrk 421). It also means that relatively close
objects, such as some of those in our Galaxy or in the Magellanic Clouds, will be a#ected
by absorption due to the CMB. It is unclear, however, what sources, if any, there are of
PeV gamma rays in these regions.
Notable among the Galactic sources is Cygnus X­3, an X­ray binary system. It is
believed that this system consists of a neutron star (a small, very dense star) orbiting
a normal stellar companion. Accretion from the companion onto the neutron star pro­
duces the X­rays that were first detected from this system (hence the name: X­3 meaning
the third X­ray source found in Cygnus). The neutron star is su#ciently close to the
companion that high energy particles emitted by the neutron star may interact with the
companion's outer atmosphere, and produce gamma ray emission. This system is far
enough away (# 10 kpc) for absorption to have a noticeable e#ect. There have been
detections reported of Cygnus X­3 at TeV and greater energies (e.g. Lloyd­Evans et al.
[8] and Samorsky and Stamm [9]), although recent observations at TeV and PeV energies
(Aglietta et al. [10]) have cast some doubt on the significance of these detections.
Objects of interest in the Magellanic Clouds (at a distance of about 50 kpc) include the
X­ray binary systems LMC X­4 and SMC X­1 (which are of a similar type to Cygnus X­3
­ a small, dense object orbiting a normal stellar companion), as well as the X­ray pulsar
PSR 0540­693 (a neutron star, spinning rapidly and giving o# pulsed X­ray emission).
Searches have been made attempting to detect TeV or greater energy emission from these
objects (Allen et al. [11]), but no significant emission was found.

Chapter 3
Infrared Background
3.1 Intensity at 60µm
3.1.1 Background
It is possible to calculate the background intensity of the infrared background (IRB) at
a particular frequency or wavelength fairly easily. In this section I'll briefly describe the
background to the calculation of the intensity at a wavelength of 60µm.
There are three quantities used to descibe the amount of radiation from a source ­
luminosity, flux, and intensity. The luminosity is the rate of energy output from a source.
The total luminosity of a source, L is in SI units of W. One can also define the specific,
or spectral luminosity L # by L = # L # d#, with SI units of W Hz -1 .
The flux is the energy flow through an area. Here I use the definitions from Rybicki and
Lightman [28], where flux is the energy carried by all rays through a given area per unit
time. As for luminosity, one can define the total flux, S (sometimes F is used instead),
with SI units W m -2 , and the specific flux S # , defined by S = # S # d#, with SI units
W m -2 Hz -1 .
Similar to flux is the intensity, I. This is defined as the energy carried through a given
area by rays within an infinitesimal solid angle d# The SI units for I are W m -2 sr -1 , and
the SI units for I # (which is defined similarly to L # and S # ) are W m -2 Hz -1 sr -1 . Another
unit often used for intensity is the Jansky, Jy, defined as 1 Jy = 10 26 W m -2 Hz -1 sr -1 .
Thus, the units of flux are Jy sr. Throughout this thesis, I # , S # and L # will be refered to
as intensity, flux and luminosity respectively, that is, the term ``specific'' will be dropped.
The key ingredient in calculating the intensity is the luminosity function, #(L # , z). This
gives the number of sources of luminosity L # at a redshift z, per comoving volume, per
unit of luminosity. Its units are Mpc -3 (W Hz -1 ) -1 . Comoving volume is a volume in
comoving coordinates, defined such that the volume in these coordinates remains constant
as the universe expands (i.e. the coordinates move with the expansion). This is di#erent
to the normal physical volume, which increases with the expansion. A definition of the
comoving volume is given in Appendix A.
A general form of a luminosity function covering sources at all redshifts has three
components: the local luminosity function # 0 (z), dealing with sources at redshift z = 0;
and functions that deal with the evolution with redshift (i.e. evolution over time). This
evolution can described in two di#erent ways. Firstly, the density of the sources can
change with time. Since the density is defined as number per comoving volume, the only
way for the density to change is for the number of sources to change. This will happen
6

CHAPTER 3. INFRARED BACKGROUND 7
predominantly by the creation of sources (i.e. galaxy formation) and by the merging of
galaxies. This form of evolution is described by the function g(z), defined by:
#(L, z) = g(z) â #(L, z = 0)
The second form of evolution is how the actual luminosity of the sources changes with
time. This is due, for instance, to changes in the physical processes involved in the source.
This is described by the function f(z), defined by:
L(z) = f(z) â L(z = 0)
So, we can combine these two relations, to give an expression for the luminosity function
in terms of the local luminosity function at a luminosity L # , at frequency #, and the two
evolution functions:
#(L # , z) =
g(z)
f(z)
# 0
# L #
f(z)
# (3.1)
As the universe expands, any radiation that has been emitted will be redshifted, mean­
ing its frequency will decrease. Thus, a photon emitted at redshift z, and observed with
frequency # will have been emitted at a frequency # # = #(1 + z). Also, if you observe a
flux S # , then the luminosity of the source at # # will be given by
L # # d# # = 4#d 2
L S # d#
L # # = S #
4#d 2
L
(1 + z)
(3.2)
where d L is the luminosity distance (see Appendix A for full expression of d L ).
However, since we will want to use L # , we need to make use of the spectrum of lu­
minosities, to relate L # to L # # . This is done by assuming a simple power law L # # -# ,
giving
L # = L # #
# # #
#
# #
or, using the relation between # # and #,
L # = L # # (1 + z) # (3.3)
where # is called the spectral index.
Using the luminosity function, we are able to evaluate the di#erential source counts.
The number of sources with observed flux greater than S # at frequency # is denoted
N(> S # ). This is given by the integral
N(> S # ) = # dz
dV c
dz
# #
L # #
dL # # #(L # # , z) (3.4)
where dV c /dz is the rate of change of comoving volume with redshift. The di#erential
source count, n(S # ) is defined as the number of sources with observed flux in the region
S # # S # + dS # per unit flux, per unit solid angle, with units of Jy -1 sr -2 . This is related
to the derivative of N(> S # ) with respect to flux:
n(S # ) =
1
4#
dN
dS #
= # dz
dV c
dz
d 2
L
(1 + z)
#(L # , z)(1 + z) # (3.5)

CHAPTER 3. INFRARED BACKGROUND 8
where L # = L # # (1 + z) #
The extra factor of (1 + z) # comes from the definition of #, as being the number of sources
per comoving volume, per unit of luminosity. Since L # = L # # (1 + z) # , we must have
d
dL # #
= (1 + z) # d
dL #
.
To get the intensity, I # (in units of Jy), you simply integrate the sourcecounts over
observed flux, to give
I # = # n(S # ) S # dS # (3.6)
3.1.2 Results
The calculations that I made were to calculate the intensity at 60µm for di#erent lumi­
nosity functions and evolutionary models.
The first luminosity function that I used was the 60µm local luminosity function from
Hacking et al. [16]:
# 0 (L 60 ) = 2.94 â 10 28 10 Q (3.7)
Q = Y - # B 2 + # log L 60 -X
W
# 2
# 1/2
- 2.5 log L 60
Here I use the values quoted in Protheroe and Biermann [17], of B = 1.51, W = 0.85,
X = 23.96, and Y = 5.93. Note that these values assume that H 0 = 100 km s -1 Mpc -1 ,
and q 0 = 0.5. This is the best fit to the data of Soifer et al. [18], as computed by [16].
The data in question are from the Infrared Astronomical Satellite (IRAS) all­sky survey.
There were three types of evolution that I considered with this model, using those
considered in [17]. The first is no evolution. This means you set f(z) = 1 and g(z) = 1,
so that:
#(L # , z) = # 0 (L # )
Secondly, I considered pure luminosity evolution, where g(z) = 1 and
f(z) = # (1 + z 0 ) 4 z # z 0
(1 + z) 4 z < z 0
(3.8)
where I used z 0 = 0.8, the best fit given by [17]. Finally, I used the evolution given by
Table 1. in Condon [19], which gives values of log(f) and log(g) for values of log(1 + z)
between 0 and 1.For values of z not listed, the log(f) and log(g) values are obtained by
interpolating or extrapolating as necessary. These functions are shown in Fig 3.1.
The other luminosity function that I used was the 60µm luminosity function from
Saunders et al. [20], with density evolution, given by:
#(L, z) = (1 + z) 6.7 #(L) (3.9)
where
#(L) =
C
L log 10
# L
L #
# 1-#
exp # -
1
2# 2
log 2
10
# 1 +
L
L #
##
with: C = (2.6 ± 0.8) â 10 -2 h -3 Mpc -3 , # = 1.09 ± 0.120, # = 0.724 ± 0.031, and
L # = 10 (8.47±0.23) h -2 L# . Also note that this requires H 0 = 66 km s -1 Mpc -1 , instead

CHAPTER 3. INFRARED BACKGROUND 9
Figure 3.1: Evolution functions. Solid and dotted lines are log(f(z)) and log(g(z)) respectively
from Condon [19]. Dashed line is log(f(z)) for the pure luminosity evolution (Eq. (3.8)) of
Protheroe and Biermann [17].
of the value of 100 that I was using for the first luminosity function. This will alter the
calculations of dL and dV c /dz.
The extra (L log 10) -1 factor compared to [20] is because their luminosity function is
defined as the di#erential luminosity function per decade of luminosity, i.e.
d(log 10 L) = d # log L
log 10
#
=
dL
L log 10
so
dN
d(log 10 L)
=
1
L log 10
dN
dL
.
Given a luminosity function and the relevant evolution functions, it is fairly straightfor­
ward to calculate the 60µm intensity, using Equations (3.2), (3.5), and (3.6). Care has to
be taken, however, when considering the spectral index # (Eq.( 3.3)). A source generally
does not have a single spectral index, but instead some distribution of indices about some
mean. To consider this, I used the following gaussian distribution approximation from
Hacking, Condon and Houck [16].

CHAPTER 3. INFRARED BACKGROUND 10
The fraction of sources with given luminosity L and redshift z having spectral index #
is well approximated by
u(#|L, z) =
1
# 2##(L, z)
exp # -[# - ##(L, z)#] 2
2# 2 (L, z)
# (3.10)
where #(L, z) = 0.5 is the dispersion and ##(L, z)# is the mean spectral index for a given
luminosity and redshift, parametrised by
##(L, z)# = # # 25µm for z > 0.41
# 25µm # log(1+z)
log(1.41) # + # 60µm (L, z = 0) # log[1.41/(1+z)]
log(1.41) # for z # 0.41
where
# 60µm (L, z = 0) =
# # # # # # #
+2.7 for L < 10 22
+2.7 - 0.5 log # L
10 22 W Hz -1 # for 10 22 < L < 10 25.4
+1.0 for L > 10 25.4
and # 25µm = +2.4. Here, L is measured in units of W Hz -1
The function u(#|L, z) is e#ectively a weighting function. We can replace the parts
of Eq.( 3.5) that depend on # by their weighted mean. Since L # = L # # (1 + z) # , #(L # , z)
depends on #, so we replace #(L # , z)(1 + z) # by
# #
#(L # # , z)(1 + z) # u(#|L # # , z).
Note that L # # is the luminosity at emission.
From here it is fairly straightforward to calculate the di#erential sourecounts n(S # ) as
a function of S # , using Eq. (3.5), for the di#erent evolutionary models discussed earlier.
These are shown in Fig. 3.2.
Similarly, the intensity can be calculated by integrating the sourcecounts according to
Eq.( 3.6). I integrated this over the range of S # values shown in Fig. 3.2, i.e. S # = 10 -2 to
S # = 10 2 (in the case of the luminosity function from Saunders et al. [20], the minimum
value of S # is 0.6). This gives the following results for the intensity at 60µm (5 â10 12 Hz):
Model I # (W m -2 Hz -1 sr -1 )
No evolution 1.55 â 10 -21
Pure luminosity evolution (Eq. (3.8)) 3.69 â 10 -21
Evolution from Condon [19] 4.27 â 10 -21
Luminosity function from Saunders et al. [20] 1.51 â 10 -22
The program used to calculate the pure luminosity evolution case is included in Ap­
pendix B.
3.2 Observational results
The extragalactic IRB is very hard to measure directly, as there is a lot of contamination
from foreground sources, in particular interplanetary dust in the solar system, as well as
interstellar dust and other emissions from the Galaxy. To get an accurate estimate of
the IRB, one needs to find models for each of these emissions and remove them from the
observed spectrum.

CHAPTER 3. INFRARED BACKGROUND 11
Figure 3.2: Di#erential sourcecounts. Solid line is no evolution, dotted line is pure luminos­
ity evolution (Eq. (3.8)), dashed line is evolution from Condon [19], and dot­dash line is the
luminosity function from Saunders et al. [20].
An important recent result comes from Puget et al. [25], using data obtained from the
Far Infrared Absolute Spectrometer (FIRAS) on board COBE. Their technique involved
removing model contributions from these foreground sources, as well as the microwave
background (since these measurement were made in a region of the spectrum dominated
by the CMB). The result is a spectrum that roughly follows a power law of
#I # # 3.45 â 10 -9 # #
400µm
# -3
W m -2 sr -1
for 400µm # # # 1000µm. For I # in terms of frequency, this becomes
I # # 4.5 â 10 -21 # #
7.5 â 10 11 Hz
# 2
W m -2 Hz -1 sr -1
for 7.5 â 10 11 Hz # # # 3 â 10 11 Hz.
Another recent result from COBE, this time from the Di#use Infrared Background
Explorer (DIRBE), is that from Hauser [27], quoted in Kashlinsky et al. [26]. The residual
results from [27], that is, the measurements after some model of foreground emission has
been removed, are summarised in the following table, where the values given are the
upper limits of the residuals (i.e. the values with the least amount of foreground emission
removed).

CHAPTER 3. INFRARED BACKGROUND 12
# (µm) # (Hz) I # (W m -2 Hz -1 sr -1 )
1.25 2.4 â 10 14 4.3 â 10 -22
2.2 1.4 â 10 14 1.9 â 10 -22
3.5 8.6 â 10 13 3.1 â 10 -22
These points are shown on each of the plots in Fig. 3.3, as is the power law from [25].
3.3 Theoretical Predictions
Because of the di#culty in observing the IRB directly, a lot of work has been done to
find suitable models for it. In recent years, there have been a number of predictions made
regarding the IRB over a range of frequencies. In this section, I review several of these
and compare their predictions with each other.
A quite useful result is that of Macminn and Primack [12], who give a comprehensive
prediction covering wavelengths from 0.1µm to 1000µm (3â10 15 Hz to 3â10 11 Hz). Their
prediction is based on di#erent models of galaxy formation: either with cold dark matter
(CDM) (consisting of a universe of 90% cold dark matter) or a mixture of cold and hot
dark matter (CHDM) (60% CDM and 30% massive neutrinos). The other 10% in each
case is the normal, baryonic matter.
Besides the composition, the main di#erence between the two models is when the
era of star formation (redshift z f ) takes place, as this is the main source of the infra­
red background. In the CDM models, galaxies form earlier, at formation redshifts of
1 < z f < 3, while in the CHDM models, they form later, 0.2 < z f < 1. Both sets of models
represent average properties of starburst and other infrared luminous galaxies, as well as
normal galaxies. Starburst galaxies are those whose output is dominated by large numbers
of young, massive, hot stars. Dust in the galaxies contributes to the infrared emission by
absorbing radiation and re­emitting in the infrared. (See, for example, Calzetti [21], or
the recent review by Morwood [22].) Sources not included in the modelling would either
be at high redshift (in which case there is no known population to identify them with), or
would be intrinsically dim (which would produce little e#ect to the overall background).
They do caution, however, that they have not considered more ``exotic'' mechanisms for
producing the IRB, so the lowest of their curves should be taken as the lower limit to the
IRB. In Fig 3.3, I have shown the upper and lower envelopes of the set of curves.
De Zotti et al. [13] and Franceschini et al. [14] both give predictions for the IRB in the
far infrared (long wavelengths or low frequencies). The predictions from [13] bound what
they term a ``conservative allowed region'', that is, a region where the IRB is most likely
to be found, since the lower curve assumes no evolution and is inconsistent with counts at
60µm, and the lower curve assumes a rate of evolution that is inconsistent with the low
number of faint high redshift galaxies observed by IRAS. The curve plotted on Fig 3.3
from [14] is the expected minimal contribution from galaxies, assuming no evolution, and
is therefore to be taken as a lower limit. Both of these papers also give predictions for the
integrated emission of starlight from galaxies in the near infrared (1 - 10µm).
Another useful prediction is from V˜ais˜anen [15], this time in near infrared to optical
wavelengths. He models a number of di#erent classes of galaxies, taking into account
di#erent evolutionary models. I have included the upper bound of his predictions (taken
from Fig.7 of [15]).
Finally, a recent paper by Stecker and De Jager [23] includes two theoretical predictions,
based on the models of Franceschini et al. (1994) [24] and normalised to fit the observational

CHAPTER 3. INFRARED BACKGROUND 13
data from Puget et al. [25].
These various models are shown on Fig. 3.3, in comparison with the observational
constraints from Puget et al. [25] and Hauser [27]. Also shown are the results of my
calculation of the IRB at 60µm. These points correspond to, from the top: evolution
functions from Condon [19]; pure luminosity evolution according to Eq. (3.8); no evolution
(f(z) = g(z) = 1); and the luminosity function and evolution from Saunders et al. [20].
The dashed line in each plot is the CMB, discussed briefly in the next section.
3.4 Microwave Background
As well as the IRB, the other background radiation field of interest for gamma ray astron­
omy is the cosmic microwave background (CMB). This is the relic photon field that is left
over from the big bang. The photons have been travelling through space since the surface
of last scattering at redshift z # 1100.
It has been shown by the FIRAS instrument on board COBE (Fixsen et al. [30]) that
the microwave spectrum is almost an exact blackbody, with deviations of less than 50
parts per million, with a temperature of T = 2.728 ± 0.004 K. Since it is a blackbody
spectrum, its intensity is given by the Planck formula
I # =
2h# 3
c 2
1
exp( h#
kT ) - 1
where h is Planck's constant, and k is the Boltzmann constant. The blackbody spectrum
is shown on each of the plots in Fig. 3.3, demonstrating how the intensity of the CMB is
far greater than that of the IRB.

CHAPTER 3. INFRARED BACKGROUND 14
Figure 3.3: Plot of infrared background intensity predictions. The vertical axes in each case
are log(I # ), where I # is in units of W m -2 Hz -1 sr -1 . Each plot shows the power law from the
observations of Puget et al. [25] (thick line), the upper limits of the residuals from the observations
from Hauser [27] (triangles), and the results from my program, calculating the background at
60µm (+ signs). Also shown is the microwave background (dashed line). a. Upper and lower
envelopes of the curves from Macminn and Primack [12]. b. Upper and lower curves are curves
1. and 2. repssectively from Stecker and de Jager [23]. The curves are shown continuing into
the region dominated by the CMB. c. The upper and lower solid curves are curves c) and a)
respectively from DeZotti et al. [13]. The dotted curve is curve d) from the same paper. d. The
solid curve is curve G2 from Franceschini et al. [14], while the dotted curve is curve G3 from the
same paper. The dot­dash curve is the upper envelope from V˜ais˜anen [15].

Chapter 4
Interactions
4.1 Introduction
When considering the propagation of gamma­rays through an extragalactic background
radiation field (either infra­red or microwave), there are two main interactions that will
alter the spectrum of the gamma­rays: photon­photon pair production, where the gamma­
ray interacts with a soft (infrared or microwave) photon to produce an electron­positron
pair; and inverse compton scattering, where an electron produced via pair production
collides with a soft photon, increasing its energy.
These processes will alter the shape of the gamma­ray spectrum emitted from such
sources discused in the first chapter. Pair­production in the soft photon field will decrease
the spectrum, whereas the inverse compton scattering by relativistic electrons can increase
the energy of the soft photons up to gamma­ray energies.
4.2 Cross sections
The cross section, #, is a key quantity in the analysis of such interactions. It is a measure
of the probability of the interaction occuring, for given values of the initial energies and
momenta. The units of # are those of area, e.g. cm 2 , so # gives the e#ective cross sectional
area of the interaction. Also of use is the di#erential cross section d#/d# which is the
cross section for the interaction given that the final momenta lie in a di#erential range of
solid angle d#
4.2.1 Photon­photon pair production cross sections
Photon­photon pair production occurs when two photons interact to produce an electron­
positron pair: # + # # e + + e - . Throughout this thesis, there will be no distinction made
between the e + and e - particles, and the term ``electron'' will be used to refer to both
particles. The interaction is considered in the centre of momentum (CM) frame, where
the two photons have equal energy. By energy conservation, the energies of the electrons
in the CM frame are the same as those of the photons. The di#erential cross section for
such a process can be shown (a full derivation is given in Jauch and Rohrlich [29]) to be
d# ##
d cos #
=
#r 2
0
2
# # m e c 2
E #
# 1 - # 4 cos 4 # + 2(m e c 2 /E # ) 2 # 2 sin 2 #
(1 - # 2 cos 2 #) 2
. (4.1)
15

CHAPTER 4. INTERACTIONS 16
Figure 4.1: Total cross section for photon­photon pair production vs. s.
Here, m e c 2 is the electron rest mass energy, E # is the energy of one of the photons in the
CM frame, r 0 is the classical electron radius with the value r 0 = 2.82 â 10 -13 cm, # is the
scattering angle, that is, the angle the electron directions make with respect to the photon
directions, and # is defined by
# =
# # # # 1 - # m e c 2
E #
# 2
=
# # # # 1 - # 4m 2
e c 4
s
#
where s = (2E # ) 2 , the total CM frame energy squared, and m e c 2 is the rest mass energy
of the electron, with the value m e c 2 = 5.11 â 10 -4 GeV. Note that since Eq. (4.1) is
independent of azimuthal angle #, the solid angle is given by
d# = d# d cos # = 2# d cos #,
and so
d# ##
d#
= 2#
d# ##
d cos #
This can be integrated to give the total cross section for photon­photon pair production:
# ## =
#r 2
0
2
(1 - # 2 ) # (3 - # 4 ) ln
1 + #
1 - # - 2#(1 - # 2 ) # (4.2)
As seen in Fig 4.1, # ## increases with s from the threshold value of s = (2mc 2 ) 2 to a
maximum, then decreases as the energy increases.

CHAPTER 4. INTERACTIONS 17
4.2.2 Inverse compton scattering cross sections
Inverse compton scattering occurs when a relativistic electron collides with a photon,
increasing its energy. This is di#erent to the usual case of compton scattering, where
energy is transfered from the photon to the electron. In this case, the final photon energy
# 1 , in terms of #, the initial energy, and #, the angle through which the photon is scattered,
is
# 1 =
#
1 + #
mec 2 (1 - cos #)
.
This equation applies to the case of an initially stationary electron. Thus, inverse compton
scattering can be considered by boosting from the LAB frame to the electron's rest frame
(ERF), calculating the change in energy of the photon by the above equation, then boosting
back to the LAB frame.
The di#erential cross section can be calculated (e.g. in [29]) to be (in the ERF)
d# #e
d#
=
r 2
0
2
# # 1
#
# 2
# # 1
#
+
#
# 1
- sin 2 # # (4.3)
where # 1 is given by the equation above. This expression can be integrated to give the
total cross section (which is independent of the choice of frame) as follows
# #e =
3
4
#T # 1 + x
x 3
# 2x(1 + x)
1 + 2x - ln(1 + 2x) # +
ln(1 + 2x)
2x -
1 + 3x
(1 + 2x) 2
# (4.4)
where x # #/mc 2 , and #T is the Thomson cross section, #T = 0.655 â 10 -24 cm 2 . This
is the classical value for scattering of low energy photons o# electrons, as first done by
J.J.Thomson.
For small energies, # # m e c 2 , there is negligible change in the photon energy, # 1 # # or
# 1 /# # 1, so the di#erential cross section becomes
d# #e
d# # r 2
0
2
(2 - sin 2 #) =
r 2
0
2
(1 + cos 2 #).
This can be integrated, using
d# = 2#d cos #,
# = #r 2
0 # 1
-1
(1 + cos 2 #)d cos #
= #r 2
0 # cos # +
1
3
cos 3 # # cos #=1
cos #=-1
=
8#
3
r 2
0 = #T
As the energy of the photon decreases, it moves into this classical regime, and the cross
section becomes constant, as shown in Fig 4.2 (here the cross section is plotted against s,
the total CM frame energy squared, which is directly related to the photon energy in the
ERF by a Lorentz boost). As the energy increases, the cross section decreases from its
classical value.
4.3 Interaction lengths
An important quantity to know when dealing with these interactions over large (e.g.
intergalactic) distances is the mean interaction length, or mean free path. This is the

CHAPTER 4. INTERACTIONS 18
Figure 4.2: Total cross section for inverse compton scattering vs. s.
average distance an incident particle (photon, electron, or otherwise) will travel before
undergoing an interaction. The mean free path of a particle of energy E through a radiation
field is denoted by #(E), and is given in general by the following integral
1
#(E)
=
1
8#E 2
# #
# min
n(#)
# 2
# smax
s min
#(s)(s -m 2 c 4 ) ds d# (4.5)
where: the incident particle has mass m and velocity #c; n(#) is the di#erential number
density of photons of energy # in the radiation field; s is the CM frame energy squared,
given by
s = m 2
c 4
+ 2#E(1 - cos #) (4.6)
with # the angle between the directions of the incident particle and the soft photon (i.e.
the photon from the radiation field); and #(s) is the total cross section for the interaction
in question, as a function of s. This equation applies for all types of particles incident on
a radiation field, be they photons, electrons, protons, or otherwise. Of course, I shall only
consider gamma­ray photons and electrons.

CHAPTER 4. INTERACTIONS 19
Figure 4.3: Mean interaction lengths for gamma rays against pair production in the microwave
and infra­red backgrounds. The dashed line gives the mean interaction length in the cosmic
microwave background. The lower and upper dashed lines give the mean interaction lengths
using the upper and lower envelopes respectively of the IRB predictions from Macminn and
Primack [12].
4.3.1 Interaction lengths for pair production
For gamma­rays undergoing photon­photon pair production, m = 0 and # = 1, so the
integral becomes
1
# ## (E)
= 1
8E 2
# #
# min
n(#)
# 2
# smax
s min
# ## (s)s ds d# (4.7)
where s = 2#E(1 - cos #), s max (#, E) = 4#E and s min = (2m e c 2 ) 2 , since there are two
electrons produced in the CM frame, and they each must have an energy of at least their
rest mass, m e c 2 . The cross section, # ## , is the total cross section for photon­photon pair
production, and is given by Eq. (4.2).
Since s min = (2m e c 2 ) 2 , a value for # min can be obtained. The lowest energy collision
will be a head on one (# = #), involving the lowest energy photon, of energy # min . Thus
(2m e c 2 ) 2 = 4E# min
# min =
(m e c 2 ) 2
E
.

CHAPTER 4. INTERACTIONS 20
The di#erential number density n(#) is given by
n(#) =
4#
hc
I #
h#
where I # is the specific intensity of the radiation field at frequency # = #/h. Thus, given
an intensity spectrum of a radiation field, such as the blackbody spectrum of the CMB or
one of the IRB spectra discussed in Section 3.3, one can calculate the mean free path of a
gamma ray through the field against photon­photon pair production interactions.
Firstly, consider the CMB, which is a blackbody spectrum of temperature 2.728 K [30].
The specific intensity of a blackbody spectrum of temperature T is given in terms of # by
I # =
2h# 3
c 2
1
exp( h#
kT ) - 1
where k is the Boltzmann constant. Since # = h#, this can be rewritten in terms of #.
I # =
2# 3
h 2 c 2
1
exp( #
kT ) - 1
Hence, one gets the expression for the di#erential number density of CMB photons as a
function of energy:
n(#) =
8#
h 3 c 3
# 2
exp( #
kT ) - 1
(4.8)
This expression can then be used in Eq. (4.7) to calculate exactly the mean free path of
gamma rays through the CMB.
For the IRB, such as that calculated in [12], n(#) is calculated directly from the pre­
diction of I # (#). In my case, this involves digitising the data from the plot in the relevant
paper, to form a series of discrete data points. Then, in calculating the integral numeri­
cally, I interpolate between the data points to get the necessary values of I # . In Fig 4.3, I
have included the CMB along with the results from the envelope curves of [12]. Obviously,
below a gamma­ray energy of about 10 5 GeV, the mean interaction length in the IRB is
much less than that in the CMB, meaning the gamma­ray will interact with an infra­red
photon, which is of higher energy than a microwave photon. Above about 10 5 GeV, the
opposite is true, with the gamma­ray interacting primarily with the CMB.
4.3.2 Interaction lengths for inverse compton scattering
For inverse compton scattering, the incident particle is an electron, which is massive
(m e #= 0) and has # < 1. Thus, we can use Eq. (4.5) without alteration, except to define
the integration limits. Since s = m 2
e c 4 + 2#E(1 - # cos #), the maximum and minimum
values of s for given values of E and # are
s max = m 2
e c 4 + 2#E(1 + #)
and s min = m 2
e c 4 + 2#E(1 - #)
respectively. In practice, since the electrons being considered have energies # 100 GeV,
# is very close to 1, so s min # m 2
e c 4 . Also, there is no threshold value for # for inverse
compton scattering (unlike pair production, where there needs to be enough energy to

CHAPTER 4. INTERACTIONS 21
Figure 4.4: Mean interaction lengths for electrons against inverse compton scattering in the
microwave and infra­red backgrounds. The dashed line gives the mean interaction length in
the cosmic microwave background. The lower and upper solid lines give the mean interaction
lengths using the upper and lower envelopes respectively of the IRB predictions from Macminn
and Primack [12].
create at least two electron rest masses), so # min = 0. So for an incident electron with
energy E, travelling at velocity #c, the mean free path # #e (E) is given by
1
# #e (E)
=
1
8#E 2
# #
0
n(#)
# 2
# smax
s min
# #e (s)(s -m 2
e c 4 ) ds d# (4.9)
where # #e is the total cross section for inverse compton scattering, given by Eq (4.4). The
calculations for n(#) are as for pair production, so the mean free path can be calculated
for both the IRB and the CMB in the same way. However, one notices from Fig 4.4 that
the mean free path in the CMB is around three orders of magnitude less than in the IRB,
so inverse compton scattering in the IRB can be neglected. That is, a relativistic electron
will interact with a microwave photon far earlier than with an infra­red photon.
Note from Fig 4.4 that below about 1000 GeV, # is constant. This is because, for low
energies, the compton scattering cross section # #e is given by the Thomson cross section,
#T = 0.655 â 10 -24 cm 2 . The mean free path equation, Eq. (4.9), can be written as
1
# #e (E)
= # #
0
n(#) # 1
-1
# #e (s)
2
(1 - # cos #) d(cos #) d#,

CHAPTER 4. INTERACTIONS 22
and if # = #T , then
1
# #e (E)
=
#T
2
# #
0
n(#) # 1
-1
(1 - # cos #) d(cos #) d#
= #T # #
0
n(#) d#.
Hence the value of # at low energies is completely determined by the di#erential photon
density. If the exact solution for the blackbody photon density, Eq. (4.8), is used, a value
for the mean free path at low energies can be obtained:
1
# #e (E)
= 16##T
h 3 c 3
(kT ) 3 #(3),
where #(3) # 1.202057, and k is the Boltzmann constant, with the value k = 8.617 â
10 -14 GeV K -1 . Evaluating this expression gives the value of # #e (E) = 3.624 â 10 21 cm =
1.174 kpc
4.4 Geometry of Interactions
This section is concerned with the actual interaction itself. It looks at the directions of the
particles involved, the lorentz transformations needed to get from one frame to another,
and how the energies of the final particles are computed. However, the techniques used to
calculate quantities such as the initial soft photon energy, the angle of incidence and the
scattering angle will be discussed in the following chapter.
4.4.1 Pair production
First consider the photon­photon pair production interaction. A high energy photon (i.e. a
gamma­ray) of energy E interacts with a lower­energy (or soft) photon of energy #, which
is travelling at an angle # to the direction of propagation of the gamma­ray. Such a setup
is shown in Fig 4.5, where I follow the method used by Protheroe [31]. This frame is called
the LAB frame, and this is the frame in which the gamma­rays will be observed, so the
final energies of all particles need to be calculated in this frame. However, the interaction
needs to be considered in the CM frame, where the two photons have equal energy and
travel in opposite directions.
To get to the CM frame, you need to apply a boost at an angle # to the gamma­ray
direction (as shown in Fig 4.5). This angle can be calculated by equating the momenta in
the LAB frame perpendicular to the boost direction:
E sin # = # sin(# - #) = # sin # cos # - # cos # sin #
(E + # cos #) sin # = # sin # cos #
# = tan -1 # # sin #
(E + # cos #)
# (4.10)
The velocity of the boost is # # c, where # # can be found by setting the momenta in the
CM frame to be equal and opposite. The momenta in the CM frame are given in terms of
the LAB frame energies and angles by: (setting c = 1)
p # 1 = # # (E cos # - # # E)

CHAPTER 4. INTERACTIONS 23
y
f
q e
E
Boost direction
(a)
(b)
E'
E'
Figure 4.5: (a) LAB frame: Gamma ray of energy E and soft photon of energy #, showing
direction of boost needed to get to CM frame. (b) CM frame: Photon directions, as well as
boost direction. Diagram taken from Protheroe [31].
p # 2 = # # (# cos(# - #) - # # #)
where p # 1 and p # 2 are the momenta of the gamma ray and the soft photon respectively in
the boost direction in the CM frame, and
# # = (1 - # #2 ) -1/2 .
By setting p # 1 = -p # 2 , as must be the case in the CM frame, we get
# # =
E cos # + # cos(# - #)
E + #
. (4.11)
In the CM frame, the photons travel in opposite directions, so that they undergo a
head­on collision. The angle they make with the boost direction is # (see Fig 4.5), given
by the aberration formulae:
sin # =
sin #
# # (1 - # # cos #)
, cos # =
cos # - # #
1 - # # cos #
So that we have
# = tan -1 # sin#
# # (cos # - # # )
# (4.12)
The energies of the photons in the CM frame are calculated using conservation of 4­
momentum. Let k 1 and k 2 denote the 4­momenta of the gamma­ray and the soft photon

CHAPTER 4. INTERACTIONS 24
respectively in the LAB frame, and k #
1 and k #
2 the respective 4­momenta in the CM frame.
Then, by conservation of 4­momentum,
(k 1 + k 2 ) 2 = (k # 1 + k # 2 ) 2 .
Since k · k = 0 for photons, we have
k 1 · k 2 = k # 1 · k #
2
E#(1 - cos #) = E # # # (1 - cos # # )
where a # means the quantity is measured in the CM frame. This frame is defined by
requiring E # = # # and # # = #. So, the energy of each photon in the CM frame is given by
E # = # 1
2
E#(1 - cos #) (4.13)
or, in terms of s, the total CM frame energy squared,
E # =
# s
2
(4.14)
Again by conservation of 4­momentum, the electrons will be emitted each with an
energy equal to E # . However, to calculate their energies in the LAB frame, the direction
that the electrons are travelling needs to be known, in particular, momentum of each
electron in the direction of the boost.
To do this, consider the components of the momenta of the particles. Let the plane
that the two photons lie in be the x­y plane, and choose the x­axis to be in the direction
of the boost. Then the momenta of the gamma­ray and soft photon in the LAB frame are
given respectively by
“ k 1 = (cos #, sin #, 0)
“ k 2 = (cos(# - #), sin(# - #), 0)
and in the CM frame they are
“ k # 1 = (cos #, sin #, 0)
“
k #
2 = (- cos #, - sin #, 0)
where # is given by Eq. (4.12).
Now define new axes “ i, “ j, “ k defined in terms of the usual cartesian axes “
x, “
y, “
z by
“ i = “
k #
1 = (cos #, sin #, 0)
“ j =
“
y â “ i
|“y â “ i|
= (0, 0, -1)
“ k = “ i â “ j = (- sin #, cos #, 0)
This way, the photon trajectories lie on the “ i axis. The electron, when emitted, is scattered
by an angle # from the direction of the photon, so that the momentum vector of the electron
lies on a cone centred on the “ i axis, with half­angle #. To determine where on this cone the
vector lies, another parameter is needed. This is the so­called azimuthal angle, #, measured
around the ``base'' of the cone, starting with # = 0 when the momentum perpendicular to

CHAPTER 4. INTERACTIONS 25
the “ i axis is in the “ j direction. Both of these angles ­ the scattering angle and the azimuthal
angle ­ need to be calculated and/or sampled. The techniques for this are discussed in the
next chapter.
Thus, “
p # 1 can be written as
“
p # 1 = cos # “ i + sin # cos # “ j + sin # sin # “ k
or, in terms of the “ x, “ y, “
z coordinates
“
p # 1 = (cos # cos #- sin # cos # sin #, cos # sin #+ sin # cos # cos #, - sin # sin #).
The other electron has momentum “
p # 2 = -“p # 1 , since it is emitted in the opposite direction to
conserve momentum. The CM frame energy of each of the electrons is given by Eq. (4.13).
To get the energies in the LAB frame, a boost is applied in the opposite direction to
the one that changed frames from the LAB to the CM. If p 1x denotes the momentum of
electron 1 in the “ x direction (i.e. the direction of the boost), then the LAB frame energy
of that electron is
E 1 = # # (E # + # # p 1x c),
so that the energies of the two electrons are given by
E 1,2 = # # # E # ± # # (##mc 2
)(cos # cos #- sin # cos # sin #) # (4.15)
where E 1 has a + sign and E 2 has a - sign. The factor of ##mc 2 comes from the magnitude
of the momentum in the CM frame. Since “
p 1 is a unit vector, p 1 “
p 1 gives the vector #p 1 ,
with magnitude
p 1 = #mv = #mc#
# p 1 c = ##mc 2
where v = #c is the velocity of the electrons in the CM frame, and # = (1 - # 2 ) -1/2 .
4.4.2 Inverse compton scattering
The analysis of the inverse compton scattering interaction is similar to the above. In the
LAB frame, a relavtivistic electron with energy E interacts with a soft photon of energy
# 1 , which is travelling at an angle # to the direction of the electron, as shown in Fig 4.6.
The interaction itself is considered in the electron's rest frame (ERF), where the electron
is stationary, and so compton scattering is used. A lorentz boost is applied in the direction
of the electron's motion to move from the LAB frame to the ERF, where the speed of the
boost is that of the electron in the LAB frame.
To calculate the energy of the scattered photon in the LAB frame, a similar procedure
is used to that for pair production. First, take the plane of the electron and photon
directions to be the x­y plane, as shown in Fig 4.6. The unit vector of momentum for the
soft photon in the LAB frame is thus given by
“ k 1 = (cos #, sin #, 0)
and in the ERF by
“
k #
1 = (cos # # , sin # # , 0),

CHAPTER 4. INTERACTIONS 26
1
(a)
q
e
E
y
x
(b)
x
y
a
incident
scattered
Figure 4.6: (a) LAB frame: Electron of energy E and soft photon of energy # 1 , with angle of
incidence #. (b) ERF frame: Photon directions, before and after scattering, showing scattering
angle #.
where, the # means measured in the ERF. The expressions involving # # come from the
aberration formulae
sin # # =
sin #
#(1 - # cos #)
, cos # # =
cos # - #
1 - # cos #
where #c is the boost velocity, and # = (1 - # 2 ) -1/2 .
Now, as before, define axes “ i, “ j, “ k in terms of the cartesian axes by
“ i = “ k # 1 = (cos # # , sin # # , 0)
“ j =
“
y â “ i
|“y â “ i|
= (0, 0, -1)
“ k = “ i â “ j = (- sin # # , cos # # , 0)
Then, when the photon is scattered with scattering angle # and azimuthal angle # (see
next chapter for the method of obtaining these angles), it has momentum unit vector in

CHAPTER 4. INTERACTIONS 27
these coordinates of
“ k # 2 = cos # “ i + sin # cos # “ j + sin # sin # “ j,
or, in the original “
x, “
y, “
z coordinates,
“
k #
2 = (cos # cos # # - sin # cos # sin # # , cos # sin # # + sin # cos # cos # # , - sin # sin #).
This gives the direction of the momentum vector of the scattered photon in the ERF.
The magnitude of the momentum is found from the energy of the photons. The energy of
the incident photon in the ERF is obtained by a simple lorentz boost
# # 1 = ## 1 (1 - # cos #),
and the energy of the scattered photon is given by the equation for compton scattering:
# # 2 =
# # 1
1 + # # 1
mc 2 (1 - cos #)
The magnitude of the momentum of the scattered photon in the ERF is then k = # # 2 /c.
Given the energy of the scattered photon, and its momentum in the “
x­direction (i.e. in
the direction of the boost), the energy of the scattered photon is given by
# 2 = ## # 2 # 1 + #(cos # cos # # - sin # cos # sin # # ) # (4.16)

Chapter 5
Monte Carlo methods
5.1 Introduction
In the previous chapter, expressions were derived for the energies of the final particles
resulting from pair production and inverse compton scattering. This derivation, however,
assumed the knowledge of quantities that were not derived, namely the soft photon energy
(#), the angle of interaction (#), the scattering angle (#), and the azimuthal angle from the
scattering (#). These quantities are each distributed according to their own distribution,
from which they need to be sampled in some way to be used in the earlier expressions.
This chapter will discuss these sampling techniques, as well as the approximations I
used to speed up the computational process. The sampling techniques discussed, which
use random numbers, are known generally as monte carlo methods. These are named after
the Monte Carlo casino in Monaco, since the roulette wheel was often used as an early
random number generator.
5.2 Interaction length
Before the quantities mentioned above are calculated, one needs to know the interaction
length ­ the distance that the incident particle has travelled before interacting. This has
to be sampled from an appropriate distribution.
The interaction length, x, is distributed according to an exponential distribution
p(x) =
1
#
exp(-x/#) (5.1)
where # is the mean interaction length, as given by Eq. (4.5). To sample from this
distribution, the transformation method is used. Let w(#) be a uniform distribution of
random numbers:
w(#) = # 1 0 # # # 1
0 otherwise
where # is a random number (or, to use the correct term, a pseudo­random deviate),
generated by a suitable pseudo­random generator. In the case of my programs, I used the
IMSL library function RNUNF().
Then, x is sampled according to
# x
0
p(x # ) dx # = # #
0
w(# # ) d# # = #.
28

CHAPTER 5. MONTE CARLO METHODS 29
In the case of the above exponential distribution,
# x
0
1
#
exp(-x # /#) = #
# - exp(-x # /#) # x
0
= #
1 - exp(-x/#) = #
x = -# ln(1 - #) = -# ln(# # )
where # # is another random number. (Since # is distributed uniformly between 0 and 1,
so is (1 - #).)
5.3 Soft photon energy
Before you are able to perform the interaction and calculate the energies of the products,
you need to know the energy # of the soft photon that the incident particle (either a
gamma­ray or an electron) is interacting with, and the angle # between the directions of
the two particles.
Firstly, the energy of the soft photon. By inspection of Eq. (4.5), the energy of the soft
photon is distributed according to
p(#) =
#(E)
8#E 2
n(#)
# 2
# smax
s min
#(s)(s -m 2 c 4 ) ds
where #(E) is the mean free path of the relevant incident particle with energy E, and # is
the total cross­section for the process in question. This distribution has been normalised
so that # #
# min
p(#) d# = 1.
However, instead of sampling from this distribution, I instead calculate the mean soft
photon energy from this distribution, by the integral
•# = # #
# min
# p(#) d#.
=
#(E)
8#E 2
# #
# min
n(#)
#
# smax
s min
#(s)(s -m 2 c 4 ) ds.
The advantage of doing this is that it is much easier and quicker than sampling from the
distribution, as •# is able to be calculated at the same time as the mean interaction length
#.
In this way, •# can be calculated for a given soft photon spectrum and interaction cross
section. These are shown for pair production and inverse compton scattering in Fig. 5.1.
Note that the range of energies for •# in the case of pair production is much greater than
that for inverse compton.
For low energies in the inverse compton case, •# is constant. This is the so­called
Thomson regime, where the compton scattering cross section is given by #T , the Thomson
cross section. Then, the analytic expression for n(#) in the CMB (Eq. (4.8)) can be
used to find the exact value for • # (following the method used for the mean free path in
Section 4.3.1).The result is that for low energies (below about 1000 GeV), the value of the
mean soft photon energy is •# = 6.366 â 10 -13 GeV.

CHAPTER 5. MONTE CARLO METHODS 30
Figure 5.1: Top: • # for photon­photon pair production. Solid upper and lower curves respectively
are for the lower and upper envelopes of the IRB prediction from Macminn and Primack [12].
The dashed curve is for the CMB. Bottom: •# for inverse compton scattering in the CMB. Note
the di#erent scales on the two plots.
5.4 Interaction angle
Once the energies of the two interacting particles are known, a third parameter is needed
­ the angle # between their respective directions. Given that the energies of the incident
particle and the soft photon are E and •
# respectively, # is defined in terms of s, the total
CM frame energy squared (defined in Eq. (4.6)), by
# = cos -1 # 1 -
(s -m 2 c 4 )
2E•#
# (5.2)
where m is the mass of the incident particle.
Thus, # can be sampled by first sampling a value of s, and then calculating # by
Eq (5.2). The advantage of this is that the distribution of s is already known from the
calculation of the mean interaction length. Properly normalised, this distribution is
p(s) =
(s -m 2 c 4 )#(s)
# smax
s min
(s -m 2 c 4 )#(s) ds
(5.3)
so that # smax
s min
p(s) ds = 1.

CHAPTER 5. MONTE CARLO METHODS 31
For pair production, where m = 0, this becomes
p(s) =
s#(s)
# 4E•#
2m 2
e c 4 s#(s) ds
(5.4)
To sample s, the transformation method is used, as detailed in Section 5.2, in which s
is sampled according to
# s
s min
p(s # ) ds # = # #
0
w(# # ) d# # = #.
So, for pair production,
# s
2m 2
e c 4
s # #(s # ) ds # = # # 4E•#
2m 2
e c 4
s # #(s # ) ds # .
Now, the integral on the left (which is a function of s) can not be inverted analytically to
give an expression for s, so it has to be done numerically.
The first step in doing this is to evaluate the integral
# smax
2m 2
e c 4
s # #(s # ) ds #
for an array of s max values, up to s max = 4E•#, storing the integrals in an array. A value
of the integral is then sampled, by sampling #, and then a value of s is obtained by
interpolating in the s max array.
The interpolation can be done quite accurately with a relatively small number of array
points if the arrays used are log( # s#(s) ds) and log(s max ), since the integral forms a power
law of s.
For inverse compton scattering, a similar process is used. However, since in this case
s = m 2
e c 4 + 2E#(1 - # cos #), a change of variable can be used:
s # x = s -m 2
e c 4 = 2E#(1 - # cos #)
so that a value of x is sampled such that
# x
2E#(1-#)
x # #(x # +m 2
e c 4 ) dx # = # # 2E#(1+#)
2E#(1-#)
x # #(x # +m 2
e c 4 ) dx # .
To see why this change of variable is useful, consider how small values of x are represented
on the computer. If s was being represented for su#ciently small values of x, it would
take the value of m 2
e c 4 , since decimal numbers can only have a limited number of decimal
places. However, if such numbers are represented by x, an exponential representation can
be utilised, with the result that far smaller values of x can be expressed. This is useful,
since in practice (1 - #) is very small for the electrons under consideration, and so s min
will be small.
5.5 Scattering Angle
Once the initial parameters E, #, and # have been determined, calculation of the interac­
tion, and in particular the energies of the final particles, can proceed. As seen in Chapter 4,
to calculate the final energy of a particle, one needs to know the angle through which it
has been scattered, that is, the angle between its direction before the interaction, and its
direction after the interaction.

CHAPTER 5. MONTE CARLO METHODS 32
5.5.1 Pair production
The scattering angle # is used in the di#erential cross section d# ## /d cos # (Eq. (4.1))
(taken in the CM frame) and this gives adistribution for cos # of
p(cos #) =
1
# ##
d# ##
d cos #
(5.5)
=
1
# ##
#r 2
0
2
# # m e c 2
E #
# 1 - # 4 cos 4 # + 2(m e c 2 /E # ) 2 # 2 sin 2 #
(1 - # 2 cos 2 #) 2
(5.6)
where # = # 1 - (4m 2
e c 4 )/s, and E # is the photon energy in the CM frame. Again, p(cos #)
is normalised so that # 1
0 p(cos #) d cos # = 1. Note that this integral is from cos # = 0
to cos # = 1, since there are two identical particles involved. In determining the total
cross section, one only integrates over those angles that correspond to physically di#erent
events, so if there are two identical (and hence indistinguishable) particles, the integration
is only made over one hemisphere, hence the integral is made from cos # = 0 to cos # = 1
instead of -1 to 1. (Note that electrons and positrons can be considered indistinguishable
particles for this problem, as they are exactly the same except for their charge, and the
cross section, and hence the distribution, is independent of charge.)
This distribution can not be sampled from easily, using any simple rule, so I instead
calculate cos #m , the median value of cos #. This is defined by the relation
# cos #m
0
d# ##
d cos #
d cos # =
1
2
# 1
0
d# ##
d cos #
d cos # (5.7)
=
1
2
# ## ,
i.e. cos #m is the value of cos # that divides the area under d# ## /d cos # in half.
This value is found by a simple search method. The value of cos # is slowly increased
from some initial value, and d# ## /d cos # is evaluated at that point and compared with
1
2 # ## . If it is smaller, increase cos # and continue; whereas if it is larger, move back a step,
decrease the step size, and continue until the step size is less than some fixed value.
For values of cos # close to 1, d# ## /d cos # increases and becomes very sharply peaked.
This means that cos #m is very close to 1 (i.e. #m is close to 0). For higher energies, this
peak becomes sharper, and cos #m becomes closer to 1. To calculate cos #m accurately, an
appproximation needs to be made for these high energies, as the di#erence between cos #m
and 1 becomes similar in size to the step size.
Since at high energies, E # # m e c 2 (i.e. photon energy # electron rest mass energy),
we make the approximation of neglecting the 2(m e c 2 /E # ) 2 # 2 sin 2 # term in d# ## /d cos #,
so that
d# ##
d cos #
=
#r 2
0
2
# # m e c 2
E #
# 1 - # 4 cos 4 #
(1 - # 2 cos 2 #) 2
.
Now, let x = cos # (and xm = cos #m ), so that
d# ##
dx
=
#r 2
0
2
# # m e c 2
E #
# 1 - # 4 x 4
(1 - # 2 x 2 ) 2
=
#r 2
0
2
# # m e c 2
E #
# (1 - # 2 x 2 )(1 + # 2 x 2 )
(1 - # 2 x 2 ) 2
=
#r 2
0
2
# # m e c 2
E #
# 1 + # 2 x 2
1 - # 2 x 2

CHAPTER 5. MONTE CARLO METHODS 33
Integrate this:
# xm
0
d# ##
dx
dx =
#r 2
0
2
# # m e c 2
E #
# # xm
0
1 + # 2 x 2
1 - # 2 x 2
dx
=
#r 2
0
2
# m e c 2
E #
# # log(1 - # 2 x 2
m ) - #xm #
So, the criterion for cos #m , Eq. (5.7), reads
log(1 - # 2 x 2
m ) - #xm =
1
2
log(1 - # 2 ) -
1
2
#
log # 1 - # 2 x 2
m
# 1 - # 2
# = #(xm -
1
2
). (5.8)
So, for high energies, where E # # m e c 2 , cos #m can be found by solving this transcen­
dental equation. This is generally faster and more accurate than solving by the regular
method indicated above. A simple root finding algorithm is used, to find the root of the
function
f(xm ) = log # 1 - # 2 x 2
m
# 1 - # 2
# - #(xm -
1
2
),
that is, find xm such that f(xm ) = 0.
Fig. 5.2 shows #m plotted as a function of energy. It shows the approximation that
I made (i.e. Eq. (5.8)), comparing it with the exact formula (i.e. using Eq. (5.7)), and it
demonstrates how the approximation is used at high energies.
5.5.2 Inverse compton scattering
A similar approach is used for inverse compton scattering. The distribution of cos # is
given by the di#erential cross section d# #e /d cos #:
p(cos #) =
1
# #e
d# #e
d cos #
(5.9)
=
1
# #e
#r 2
0 # # 1
#
# 2
# # 1
#
+
#
# 1
- (1 - cos 2 #) # (5.10)
where
# 1
#
=
1
1 + #
mec 2 (1 - cos #)
.
This is normalised so that # 1
-1 p(cos #) d cos # = 1. Note that this integral is from ­1
to 1, since the two particles under consideration are distinguishable, that is, they are an
electron and a photon, as opposed to two electrons or two photons in the pair production
case.
Again, the median value of cos # is found, using the following relation
# cos #m
-1
d# #e
d cos #
d cos # =
1
2
# 1
-1
d# #e
d cos #
d cos # (5.11)
The method for finding cos # is essentially the same as that for pair production ­ a
simple step search. However, again the search breaks down at high energies, where the
di#erential cross section becomes very sharply peaked at cos # = 1, so an approximation

CHAPTER 5. MONTE CARLO METHODS 34
Figure 5.2: Median scattering angle for pair production against photon energy. The solid line is
the curve used in the evaluation of #m . The dotted line shows what happens at high energies if
cos #m is calculated using Eq. (5.7). The dashed line shows the continuation of the approximation
(Eq. (5.8)) into low energies, where it is not applicable. The plot cuts o# at E # = m e c 2 , the
minimum photon energy in the CM frame.
needs to be made. For # # m e c 2 , the di#erential cross section becomes (letting x = cos #
for brevity)
d# #e
dx #
#r 2
0
1 + #
mec 2 (1 - x)
.
Integrate this:
# xm
-1
d# #e
dx
dx = #r 2
0 # xm
-1
dx
1 + #
mec 2 (1 - x)
= #r 2
0
# -m e c 2
#
log # 1 +
#
m e c 2 (1 - x) # # xm
-1
= #r 2
0
m e c 2
#
# log # 1 + 2
#
m e c 2
# - log # 1 +
#
m e c 2
(1 - x) ##
So, Eq. (5.11) reads
log # 1 + 2
#
m e c 2
# -log # 1 +
#
m e c 2
(1 - xm ) # =
1
2
log # 1 + 2
#
m e c 2
# -log # 1 +
#
m e c 2
(1 - 1) #

CHAPTER 5. MONTE CARLO METHODS 35
Figure 5.3: Median scattering angle for inverse compton scattering as a functionof photon energy
in the ERF. The solid line is the curve used to evaluate #m . The dotted line shows what happens
at high energies if cos #m is calculated using Eq. (5.11). The dashed line shows the continuation
of the approximation (Eq. (5.12)) into low energies, where it is not applicable.
log # 1 +
#
m e c 2
(1 - xm ) # =
1
2
log # 1 + 2
#
m e c 2
#
1 +
#
m e c 2
(1 - xm ) = # 1 + 2
#
m e c 2
xm = 1 +
m e c 2
# -
m e c 2
#
# 1 + 2
#
m e c 2
(5.12)
Thus, at high energies, there is an analytic expression for #m , which speeds up the
process enormously. Fig. 5.3 shows how this approximation is needed at high energies,
and how it is not valid at low energies. As you can see, the agreement between the
exact result and the approximation is never exact, although I have chosen the valuefor
the change­over to minimise the di#erence between the two. The values that I used are
indicated by the solid line.
Note that below about 10 -6 GeV, #m is constant at 90 # . This is because, for low
energies the photon is in the Thomson regime, where # 1 # #, so the di#erential cross
section becomes
d# #e
dx
= #r 2
0 (1 + x 2 ).

CHAPTER 5. MONTE CARLO METHODS 36
Thus,
# xm
-1
d# #e
dx
dx = #r 2
0 # 4xm
3
+
4
3
#
and
# 1
-1
d# #e
dx
=
8#r 2
0
3
= #T
Then Eq. (5.11) simply gives xm = 0, or #m = 90 # .
5.6 Azimuthal angle
Recalling the analysis of the interactions from Section 4.4, a set of coordinates can be
defined such that, in the frame of interest (either the CM frame for pair production or
the electron rest frame (ERF) for inverse compton scattering), the photon before the
interaction moves along a given axis. Then, the resulting particle (photon or electron) has
its momentum vector lying on a cone, centred on the axis, of half angle #m , calculated
above.
To find where on this cone the vector lies, we sample uniformly an angle between 0
and 2#, and measure this from some axis perpendicular to the particle's initial direction.
Thus, the azimuthal angle # is defined by
# = 2##
where # is a random number, discussed in Section 5.2.
Thus, I have sampled two parameters from their respective distributions, namely # and
# (ignoring the distance for the moment), and found either the mean or the median values
of the other two, # and #. The next step is to combine these routines together to model a
cascade of electrons and photons in a background radiation field.

Chapter 6
Electromagnetic cascades
Electron photon cascades in the microwave and infrared backgrounds, and their e#ects on
the gamma ray spectra of distant sources have been considered in recent years by a number
of authors. Protheroe [31] used monte carlo methods to model cascades in the CMB, from
sources such as Cygnus X­3, the Large Magellanic Cloud and Centaurus A. This was
done by using matrices to describe propagation of the particles over a small distance,
and then repeating the calculations to describe larger distances. This procedure was
modified by Protheroe and Stanev [33], who considered cascading in the IRB for distances
corresponding to several AGN (Mrk 421, 3C 273, 3C 279, and NGC 4151 (an AGN not
detected by EGRET)). A more recent result comes from Biller [34]. This considers the
e#ect of the IRB on the gamma ray spectra of AGN, using a di#erent approach to the
computing by parametrising the interaction length and using pre­computed libraries to
find the energy­momenta of the secondary particles.
6.1 Program structure
In the previous two chapters, methods for calculating the results of photon­photon pair
production and inverse compton scattering have been discussed. These have been consid­
ered as individual interactions, with an incident particle interacting with a soft photon to
produce two final particles.
However, it is obvious that the electrons produced in the pair production interaction
may undergo interact themselves, producing a gamma ray which then pair produces, giving
more electrons, and so on, producing an electromagnetic cascade. The extent to which the
electrons produce gamma rays by inverse compton scattering is limited by the presence of
magnetic fields, in which the electrons will both spiral (thereby changing their direction so
that they are not travelling in the same direction as the initial particle), and lose energy
by synchrotron emission. In this model, it is assumed that the extragalactic magnetic field
is su#ciently low that these e#ects can be neglected. Thus, the electrons lose energy by
inverse compton scattering alone, and travel in straight lines.
The program shown schematically in Fig. 6.1 models a cascade consisting of photon­
photon pair production and inverse compton scattering, for an intial injection of a single
gamma ray of given energy. This was then run several hundred times (typically 250) to
obtain good statistics for the random components of the program (sampling the interaction
angle # and the azimuthal angle #). This simulates 250 photons of the same energy being
injected in the photon field.
For the variables that were not sampled from a distribution (such as #, •#, and #), I
37

CHAPTER 6. ELECTROMAGNETIC CASCADES 38
Calulate mean free path
and soft photon energy
Sample distance
Past Earth?
Is smax greater
than threshold?
Calculate
electron energies
Store
electrons
Calulate mean free path
and soft photon energy
Sample distance
Past Earth?
Calculate
photon energy
Initial gamma ray
Store
Output
photon
photon
Output
Store
photon
Is electron energy
less than minimum?
Get particle
from store
Output
photon
spectrum
Store
Empty
No
No
Yes
Yes
Yes
Yes
No
Pair Production Inverse Compton
Photon Electron
Figure 6.1: Flow diagram for the cascade program with a single gamma ray input.

CHAPTER 6. ELECTROMAGNETIC CASCADES 39
used arrays to store values of these variable for a range of energies and then interpolate in
these arrays when a value of the variable is needed. This removes the need for unnecessary
calculations (since the arrays can be calculated beforehand, and read in at the start of the
program) and so speeds up the program.
The program only deals with one particle at a time, so when particles have been
produced (either in a pair production in the case of electrons or scattered up to gamma
ray energies in an inverse compton interaction), their energy, type (i.e. electron or photon)
and distance travelled are stored in an array. Each time an interaction is complete, the
particle at the top of this array is retrieved.
Since the cascade is taken over some fixed distance, the distance travelled by each
particle is stored, and incremented each time the interaction length is sampled. The
particles are assumed to travel in the same direction as the initial gamma ray ­ that is,
in a straight line from the source to Earth. This is an approximation, but not one that
will a#ect the results greatly, since the directions of the emitted particles will be highly
collimated in the direction of the high energy particle.
When the photons have either travelled past Earth (have travelled a greater distance
than that from the source to Earth), or are not energetic enough to undergo pair production
(i.e. s max < (2m e c 2 ) 2 , where s is the total CM frame energy squared), their energy is stored
in an array ready for output. (If s max < (2m e c 2 ) 2 , the photons will not interact again, so
they travel straight through to Earth.) When the cascade has finished, that is, there are
no more particles in the store, the spectrum of photons at Earth is stored in a file.
The spectrum of electrons at Earth is not under consideration, so if an electron travels
further than Earth, it is disregarded, and a new particle is retrieved from the store. Since
the electron is not destroyed in the interaction, it is available to scatter again, provided it
is above a minimum energy that is imposed on electrons. If that is the case, the electron
is used straight away, without being stored.
The electrons have a minimum energy imposed upon them to reduce the amount of
particles produced and hence the computing time. Low energy photons have a short mean
free path, and so produce a lot of gamma rays (of low energy) in a given distance, all
of which need to be considered before the program terminates. If the minimum electron
energy is set at the lowest photon energy of interest, then any photons produced by
electrons below the minimum will have less energy, and so will not a#ect the region of
interest. This also reduces the computational time, as there will be much less photons
produced, meaning less electrons as well.
6.2 Mono­energetic cascades
Mono­energetic cascades are cascades initiated by gamma­rays of a single energy. This is
modelled by the basic program structure shown in Fig. 6.1, where a single photon of the
required energy is injected, and the program is run a large number of times to simulate the
injection of a large number of photons of the same energy (i.e. a large number of cascades).
When a photon reaches Earth, its energy is recorded. If the photon has energy between,
say, (E-#E/2) and (E+#E/2), then the number of photons received in the energy bin of
width #E, centred on energy E, is incremented. These bins cover the full range of energies
of interest. The width of the bins is a constant value of log E , that is, #(log E) = const,
so #E # E. The value of the constant used in the program was 0.1, so that there were
ten bins per decade.
When the program is finished, it takes the number of photons received in each bin

CHAPTER 6. ELECTROMAGNETIC CASCADES 40
Figure 6.2: Monoenergetic spectra resulting from cascading over a distance of 50 kpc. The
initial photon energy, E inj , is shown for each plot. The photon spectrum times E 2 (units of GeV)
is plotted against photon energy in GeV.
and divides it by the width of the bin and the number of cascades (N casc ). This gives the
di#erential photon spectrum dN/dE, defined by
dN
dE
=
# photons
N casc #E
.
This then gives the spectrum of photons produced by a cascade initiated by a single
photon energy. Examples of these are shown in Fig. 6.2, where log(E 2 dN/dE) is plotted
against log(E). The reason I multiply the photon spectra by E 2 is that dN/dE is generally
decreasing quite steeply, so multiplying it by E 2 smoothes it out and allows deviations
from a dN/dE # E -2 spectrum to be seen more clearly (many gamma ray sources have a
spectrum that is approximately E -2 ).
The spectra in Fig. 6.2 are for di#erent initial gamma ray energies (E inj ), all propagated
over the same distance ­ in this case, 50 kpc, which is about the distance to the Large
Magellanic Cloud ­ with the photons interacting with the CMB. The di#erence in the
spectra is due to the di#erence in the mean free paths of the photons through the CMB.
The values of the mean free paths are shown in the following table.

CHAPTER 6. ELECTROMAGNETIC CASCADES 41
log(E inj (GeV )) # ## (kpc)
9.0 343.2
7.5 23.2
6.0 8.2
4.5 5.5 â 10 14
Since the interaction lengths are distributed according to an exponential distribution:
p(x) =
1
#
exp(-x/#),
the number of photons remaining after travelling a distance x is
N = N casc exp(-x/#),
where N casc is the number of photons initially. Thus, the photon spectrum for the energy
bin containing the initial photon energy is given by
dN
dE
=
exp(-x/#)
#E
.
For example, for photons of energy 10 9 GeV and x = 50 kpc, this becomes
dN
dE
=
exp(-50/343.2)
#E
=
0.864
(10 9.05
- 10 8.95 )
= 3.746 â 10 -9
and so the value of log(E 2 dN/dE) is log(3.746 â 10 9 ) = 9.57. This is approximately
the value of the peak at E = 10 9 GeV on the plot in Fig. 6.2, indicating that the correct
amount of absorption is occuring.
The photons present at energies from about 10 7.5 GeV down are the result of the
cascading process. The electrons, created by those initial photons that did interact, scatter
microwave photons up to gamma ray energies. These energies will be less than the electron
energies, which in turn are less than the energies of the injected photons. As the cascade
proceeds, the energies of the particles decrease slowly, which is why there are no photons
with energies between 10 8 and 10 9 GeV.
At around 10 5 GeV, the mean free path increases rapidly, so any photons produced
at these energies are most likely to travel stright through without interacting. Thus, you
tend to get a lot of photons around these energies. This is particularly obvious for longer
distances, when the photons have more time to lose their energy through pair production
and move down to the low energy part of the spectrum.
Since the mean free path for 10 4.5 GeV photons is so large, none of them interact, so
the photon spectrum for the energy bin containing 10 4.5 GeV should be given by
dN
dE
=
N casc
N casc #E
=
1
#E
which gives, for E inj = 10 4.5 , a value of log(E 2 dN/dE) = 5.14, which is approximately
the value of the peak on the plot in Fig. 6.2.

CHAPTER 6. ELECTROMAGNETIC CASCADES 42
6.3 Propagation of a power law spectrum
In reality, however, you do not have mono­energetic cascades occuring. Instead, emission
from typical astronomical sources takes the form of a power law, dN/dE # E -# , where #
is the spectral index. Typically, the value of # is between 1.5 and 3, and mostly close to
2. This is particularly true for the AGN discovered by EGRET (von Montigny et al. [2]).
I will be assuming for my work a spectral index of # = 2, since this is the spectral index
of Mrk 421 (from [2]).
To model the injection of a power law spectrum, one simply performs mono­ energetic
cascades for the range of energies of interest, spacing the injection energies with the same
spacing as the energy bins. Once you have each of the mono­energetic spectra, you combine
them in a way that weights each spectra according to the relative weighting of the initial
energy in the initial power law.
Assuming an E -2 power law, the injected photon spectrum is set to be
dN inj
dE
= E -2 .
Then, the number of injected photons in the energy bin centered on energy E i of width
#E is given by
N inj
i = # E i +#E/2
E i -#E/2
E -2 dE
=
1
E i -#E/2 -
1
E i +#E/2
=
#E
(E i -#E/2)(E i +#E/2)
.
Then,
N inj
i
#E
=
1
(E i -#E/2)(E i +#E/2)
#
1
E 2
i
(provided #E is small), as expected for an E -2 power law.
Now, let dN i /dE j denote the mono­energetic photon spectrum for the energy bin cen­
tered on E j , due to the mono­energetic injection of photons of energy E i . Then, the output
photon spectrum in the j th energy bin is given by
dN
dE j
= # i#j
N inj
i
dN i
dE j
(6.1)
The sum is made over all bins greater than or equal to the energy in question, since any
photons that are in a given bin at Earth must have either originated in that bin, in which
case they did not interact, or they have been produced by a cascade that was originated
by a higher energy photon. This is because there is no process that enables any photons
produced to have a greater energy than the primary photon.
In this way, the final photon spectrum can be calculated for each energy bin in the
region of interest, to obtain a complete photon spectrum resulting from the injection of a
power law.

CHAPTER 6. ELECTROMAGNETIC CASCADES 43
Figure 6.3: An E -2 spectrum after cascading over a distance of 50 kpc. The histogram is the
resulting spectrum from the program described in the text. The dotted line shows the shape of
the initial spectrum, i.e. what the spectrum would look like if no interactions took place. The
dashed line is the spectrum from Protheroe [31].
6.4 Results
6.4.1 Cascading in the microwave background
The first case to consider is that shown in Fig. 6.2, cascading through the microwave
background over a distance of 50 kpc, which is approximately the distance to the Large
Magellanic Cloud. The combined output spectrum is shown in Fig. 6.3, over energies
ranging from 10 3 GeV up to 10 10 GeV.
The histogram in Fig. 6.3 is a histogram of results of mono­energetic injection produced
by the program detailed in Section 6.1, and combined to represent the evolution of an E -2
power law spectrum by the method explained above.
The qualitative features of the spectrum are as follows. At high energies, the spectrum
su#ers little absorption, because the mean interaction length is much greater than the
distance of propagation. For energies between 10 5.5 GeV to 10 8 GeV, the amount of
absorption is considerable, due to the mean interaction length being less than 50 kpc.
The amount of absorption, however, is less than would be expected from pair production
alone, since the inverse compton scattering by higher energy electrons produces photons
in this energy range to replace some of those absorbed in the pair production. At lower

CHAPTER 6. ELECTROMAGNETIC CASCADES 44
energies, the spectrum increases above the value it had upon injection, due to the addition
of photons from the inverse compton scattering of higher energy electrons, and because
the mean interaction length is much greater than the distance travelled.
The dashed line on Fig. 6.3 is the spectrum after cascading over 50 kpc as calculated
by Protheroe [31]. This was calculated using a program that sampled each parameter
correctly from its distribution, with no approximations made, unlike my program, where
I used the mean soft photon energy and the median scattering angle.
As is evident in Fig. 6.3, my approximate program gives results that almost match the
more accurately sampled results from the program by Protheroe. This indicates that the
distributions from which I took the mean (in the case of • #) or the median (in the case
of #m ) would produce values not too di#erent from • # and #m had the full sampling been
made instead. Thus, my program seems to provide a very good approximation to the more
accurately sampled results, at least up to a distance of 50 kpc.
One test to make on this program would be to run it over a larger distance. I had
originally intended to run the program over a distance of 6.4 Mpc, which is roughly the
distance to Centaurus A, a well known active galaxy, which has been suspected of being
a gamma ray emitter, although no definite evidence of this has been found. However,
such a program would take a lot longer to run, since the mean interaction lengths are a
much smaller proportion of the total distance travelled, meaning that you get a lot more
interactions occuring. The problem here is that the electrons, once they get down to low
energies, travel on average about 1 kpc, so they undergo a lot of interactions in the given
distance before they either reach Earth or drop below the cuto# energy. This large increase
in particle numbers requires much more processing time, as well as more storage space, to
store the particles before they are interacted.
Thus, time constraints prevented me from completing these calculations. This would
be a good direction to extend this work, as there are already previous calculations of
spectra from this distance (Protheroe [31]) that used the full sampling routines, and so
the e#ectiveness of the program could be evaluated for longer distances.
6.4.2 Cascading in the infrared background
The second case to consider is that of cascading in the infrared background. For this case,
I consider gamma rays emitted from a distance of 100 Mpc, comparable to the distance
to Mrk 421, a known TeV gamma ray emitter. The models of the IRB that I considered
were the upper and lower envelopes from the predictions of Macminn and Primack [12]
The resulting spectra from these cascades are shown in Fig. 6.4. The solid­line his­
togram is the spectrum resulting from cascading through the IRB given by the lower
envelope of the set of curves given by Macminn and Primack [12]. The dashed­line his­
togram uses the upper envelope of the set of curves from [12].
Since the upper envelope has the greater intensity, the mean free path through it for
photon­photon pair production is shorter, so more absorption is going to take place. This is
seen in Fig. 6.4 for energies greater than about 10 3 GeV, where the photon spectrum is less
than that for the lower envelope. However, below about 10 2.5 GeV, the upper envelope's
spectrum is greater than the lower envelope's. This is because the extra absorption via
pair production produces more electrons, which produce more lower energy photons.
An important prediction resulting from these results is that very little emission will
be received at energies above about 10 4.5 GeV (or about 30 TeV) from sources at this
distance. Thus, at these high energies, any gamma rays detected will be far more likely to

CHAPTER 6. ELECTROMAGNETIC CASCADES 45
Figure 6.4: An E -2 spectrum after cascading over a distance of 100 Mpc. The histogram is the
resulting spectrum from the program described in the text. The dotted line shows the shape of
the initial spectrum, i.e. what the spectrum would look like if no interactions took place. The
solid histogram is from cascading through the lower envelope from Macminn and Primack [12],
while the dashed histogram is from cascading through the upper envelop of [12].
have come from more local sources. Also, the further away a source is, the more absorption
will take place, and so the cuto# energy will be reduced. This means that, since these
spectra correspond to Mrk 421 ­ the closest of the EGRET AGN ­ one can just about
rule out the possibility of detecting emission at TeV energies from any of the more distant
EGRET AGN.

Chapter 7
Conclusions
Motivated by the detection by EGRET of more than 40 AGN at GeV gamma ray energies,
and the detection of only one of these at TeV energies, I have investigated the e#ect of
gamma ray­initiated cascades in the infrared and microwave backgrounds on the gamma
ray spectra from extragalactic sources.
The infrared background (IRB) has recieved a lot of attention in recent years, both
from a theoretical and observational perspective. The predictions of the shape of the IRB
spectrum that I have reviewed have been in broad agreement over the range of frequencies
of interest, and in general have been consistent with the constraints from observations from
COBE. I also made calculations of the background intensity at a wavelength of 60µm.
The spectrum of gamma rays from a distant source is modified in a background photon
field by pair production and inverse compton scattering interactions. This produces a
cascade of photons and electrons. In modelling such a cascade, it is necessary to take
into account the random nature of these interactions by using monte carlo simulations.
I constructed a program that made several approximations, including taking the mean
energy of the soft photon, and the median scattering angle, instead of sampling from the
appropriate distributions. I showed that these approximations had little e#ect on the
resulting gamma ray spectrum.
I have also demonstrated that the IRB causes considerable attenuation of the gamma
ray spectrum from Mrk 421 at energies greater than about 30 TeV. This has implications
in the area of AGN gamma ray emission, since the amount of absorption will increase
for more distant sources, which implies that the cuto# energy for the observed spectra
will decrease. Since Mrk 421 is the nearest of all the AGN observed by EGRET at GeV
energies, this result indicates that none of the other EGRET AGN will be visible at TeV
energies or greater. This is in agreement with the lack of detection of all the EGRET
AGN (except Mrk 421) at energies above the EGRET limit of 30 GeV.
46

Bibliography
[1] Thompson.D.J. et al., Ap. J. Supp., 101, 259, (1995)
[2] von Montigny.C. et al., Ap. J., 440, 525, (1995)
[3] Dermer.C.D. and Schlickheiser.R., Science, 257, 1642, (1992)
[4] Punch.M. et al., Nature, 358, 477, (1992)
[5] Petry.D. et al., Astron. Astrophys., in press, preprint astro­ph/9606159, (1996)
[6] Quinn.J. et al., Ap. J., 456, L83, (1996)
[7] Lin.Y.C. et al., Ap. J., 401, L61, (1992)
[8] Lloyd­Evans.J. et al., Nature, 305, 784, (1983)
[9] Samorsky.M. and Stamm.W., Ap. J., 268, L17, (1983)
[10] Aglietta.M. et al., Astroparticle Phys., 3, 1, (1995)
[11] Allen.W.H. et al., Ap. J., 403, 239, (1993)
[12] Macminn.D. and Primack.J.R., Space Science Reviews, 75, 413, (1996)
[13] De Zotti.G. et al., preprint astro­ph/9503007, (1995)
[14] Franceschini.A. et al., Astron. Astrophys Suppl., 89, 285, (1991)
[15] V˜ais˜anen.P, Astron. Astrophys., in press (1996)
[16] Hacking.P., Condon.J.J., and Houck.J.R., Ap. J., 316, L15, (1987)
[17] Protheroe.R.J. and Biermann.P.L., Astroparticle Phys., in press, (1996)
[18] Soifer.B.T. et al., Ap. J. (Letters), 303, L41, (1986)
[19] Condon.J.J., Ap. J., 287, 461, (1989)
[20] Saunders.W. et al., M.N.R.A.S, 242, 318, (1990)
[21] Calzetti.D., A. J., in press, preprint astro­ph/9610184, (1996)
[22] Morwood.A.F.M., Space Science Reviews, 77, 303, (1996)
[23] Stecker.F.W. and De Jager.O.C., Ap. J. (submitted), preprint astro­ph/9608072,
(1996)
47

BIBLIOGRAPHY 48
[24] Franceschini.A. et al., Ap. J., 427, 140, (1994)
[25] Puget.J.­L. et al., Astron. Astrophys., 308, L5, (1995)
[26] Kashlinsky.A. et al., Ap. J., in press, preprint astro­ph/9604182, (1996)
[27] Hauser.M.G., in Unveiling the cosmic infrared background, AIP Conf. Proc 348, pp
11­24, ed. Dwek.E. (1996)
[28] Rybicki.G.B. and Lightman.A.P., Radiative Processes in Astrophysics, John Wiley
and Sons, 1st Ed., (1979)
[29] Jauch.J.M and Rohrlich.F., The Theory of Photons and Electrons, Springer­Verlag,
2nd Ed., (1980)
[30] Fixsen.D.J, et al., Ap.J. (submitted), preprint astro­ph/9605054 (1996)
[31] Protheroe.R.J., M.N.R.A.S., 221, 769, (1986)
[32] Mandl.F. and Shaw.G., Quantum Field Theory, John Wiley and Sons, Revised Edi­
tion, (1994)
[33] Protheroe.R.J. and Stanev.T., M.N.R.A.S., 264, 191, (1993)
[34] Biller.S.D., Astroparticle Phys., 3, 385, (1995)
[35] Kolb.E.W. and Turner.M.S., The Early Universe, Addison­Wesley, 1st Ed., (1990)

Appendix A
Cosmological parameters
Two cosmological parameters were used in the calculation of the background intensity
at 60µm: the luminosity distance d L , and the rate of change of comoving volume with
redshift dV c /dz. Here I give the formulae for these, taken from Kolb and Turner [35]. The
universe is assumed to be a matter­dominated Freidmann­Robertson­Walker cosmology,
with a scale factor a(t).
The general expression for the luminosity distance is
d 2
L = a 2 (t 0 )r 2
1 (1 + z) 2
where a(t 0 ) is the scale factor of the universe at the present epoch (i.e. at time t 0 ). This
is the luminosity distance to a source that emitted light at redshift z, with comoving
coordinate r = r 1 . The light is recieved now (t = t 0 ) at coordinte r = 0.
The expression for the comoving coordinate r 1 is given by, for a matter dominated
universe:
r 1
=2#
0 z
+(2# 0 - 4) #
## 0 z + 1 - 1 #
a 0 H
0# 2
0 (1 + z)
where
# 0 =
# 0
# c
=
present density of universe
critical density
.
Then, the expression for d L becomes
d L =
zq 0 + (q 0 - 1) # # 2q 0 z + 1 - 1 #
q 2
0 H 0
where
H 0 =
—
a(t 0 )
a(t 0 )
is the Hubble parameter, and
q 0 = -˜aa
—
a 2
= -˜a
aH 2
0
is the deceleration parameter.
The unit of comoving volume is given by
dV c =
r 2
(1 - kr 2 ) 1/2
dr
d#
49

APPENDIX A. COSMOLOGICAL PARAMETERS 50
which means that the rate of change of comoving volume with redshift is, using the full
expression for r from above:
dV c
dz
=
4# # zq 0 + (q 0 - 1) # # 2q 0 z + 1 - 1 ## 2
a 3
0 H 3
0 (1 + z) 3 q 4
0 [1 - 2q 0 + 2q 0 (1 + z)] 1/2
=
4#d 2
L
a 3
0 H 0 (1 + z) 3 [1 - 2q 0 + 2q 0 (1 + z)] 1/2

Appendix B
Background intensity program
program intensity
c background_60micron_lumevol.f
c calculates the intensity for pure luminosity evolution
implicit real*8 (a­h,o­z)
c
data c/299792458.d0/ !value of the speed of light in m/s
pi=4.d0*datan(1.d0)
h0=1.d2 !value of the hubble constant in km/s/Mpc
q0=0.5d0
zmin=0.d0
zmax=9.d0
smin=1.d­2
smax=1.d2
call background(zmin,zmax,q0,h0,smin,smax,bg)
freq60=c/6.d­5
write(*,*) dlog10(freq60)
write(*,*) dlog10(bg)
end
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
function dlum(z,q0,h0)
c returns luminosity distance for redshift z in units of Mpc
implicit real*8 (a­h,o­z)
c convert Hubble const. to units of Mpc
t0=1.d0/(h0*1.d5/3.d24)
ct0=3.d10*t0/3.d24
a=(z*q0 + (q0­1.d0)*(dsqrt(2.d0*q0*z + 1.d0)­1.d0))
dlum=a*ct0/(q0*q0)
return
end
51

APPENDIX B. BACKGROUND INTENSITY PROGRAM 52
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
function dvcdz(z,q0,h0)
c returns rate of change of comoving volume with redshift,
c at redshift z, in units of Mpc^3.
implicit real*8 (a­h,o­z)
pi=4.d0*datan(1.d0)
c convert Hubble const. to units of Mpc
t0=1.d0/(h0*1.d5/3.d24)
ct0=3.d10*t0/3.d24
a=(z*q0 + (q0­1d0)*(dsqrt(2.d0*q0*z + 1.d0)­1.d0))
dvcdz=4.d0*pi*a*a*ct0**3/(dsqrt(1.d0­2.d0*q0+2.d0*q0*(1+z))
+ *q0**4*(1+z)**3)
return
end
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
function f(z)
c luminosity evolution function
implicit real*8 (a­h,o­z)
z0 = 0.8
if (z.lt.z0) then
f = (1+z)**4
else
f = (1+z0)**4
endif
return
end
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
function g(z)
c density evolution function
implicit real*8 (a­h,o­z)
g = 1.d0

APPENDIX B. BACKGROUND INTENSITY PROGRAM 53
return
end
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
function rho(alum, z)
c returns luminosity function for luminosity alum
c (in W/Hz) and redshift z
c rho is given in units of Mpc^(­3)*(W/Hz)^(­1)
implicit real*8 (a­h,o­z)
c parameter values, given that H0=100 km/s/Mpc, and q0=0.5
c Using function quoted in Protheroe & Biermann
Y = 5.93d0
B = 1.51d0
W = 0.85d0
X = 23.96d0
dummy = dlog10(alum/f(z))
Q = Y ­ dsqrt ( B**2 + ((dummy ­ X)/W)**2) ­ 2.5d0 * dummy
rho = 2.94d28 * 10**Q * g(z)/f(z)
return
end
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
subroutine sourcecounts(zmin,zmax,q0,h0,snu,sum)
c returns sourcecounts n(snu) in the variable sum
c (units = #/Jy/sr) , where snu is the measured flux in Jy,
c and zmin & zmax are the limits of the integration over redshift.
implicit real*8 (a­h,o­z)
pi=4.d0*datan(1.d0)
n=1000
dz=(zmax­zmin)/n
sum=0.0d0
do i=1,n
z=zmin + (i­0.5)*dz
alum=4.d0*pi*snu*dlum(z,q0,h0)**2*9.d18/(1+z)
c this produces the luminosity in W/Hz
c alum60 = alum * (1+z)**alpham(alum,z)

APPENDIX B. BACKGROUND INTENSITY PROGRAM 54
c this adjusts the luminosity to account for the change in
c frequency, by using the mean spectral index
c note that 1+z = freq/freq60
sigma=0.5d0
unorm=0.d0
alphasum=0.d0
alphamean=alpham(alum,z)
do dalpha = ­2.d0,2.d0,0.2d0
ualpha = dexp(­dalpha**2/(2.d0*sigma**2))
alumalpha = alum * (1+z)**(dalpha+alphamean)
if ((alumalpha.gt.1.d21).and.
+ (alumalpha.lt.1.d27)) then
alphasum = alphasum + ualpha * rho(alumalpha,z)
+ *(1+z)**(dalpha+alphamean)
endif
unorm = unorm + ualpha
enddo
alphasum = alphasum/unorm
c now perform the integral
sum = sum + dz*dvcdz(z,q0,h0)*dlum(z,q0,h0)**2
+ *9.d44*alphasum/(1+z)
c sum is now in units of m^2*Hz/W, since 1Mpc = 3.d22 m
enddo
c now convert sum into #/Jy/sr, by using the fact that
c 1Mpc^2*Hz/W = 1.d­26 Jy^(­1)sr^(­1)
sum = sum * 1.d­26
return
end
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
function alpham(alum,z)
c calculates the mean spectral index at a given
c luminosity alum (in W/Hz) and redshift z.
implicit real*8 (a­h, o­z)
data alpha25/2.4d0/
if (z.gt.0.41d0) then

APPENDIX B. BACKGROUND INTENSITY PROGRAM 55
alpham = alpha25
else
if (alum.lt.1.d22) alpha60 = 2.7d0
if (alum.gt.(10**25.4)) alpha60 = 1.d0
if ((alum.lt.(10**25.4)).and.(alum.gt.1.d22))
+ alpha60 = 2.7d0 ­ 0.5d0*dlog10(alum/1.d22)
alpham = alpha25*(dlog10(1.d0+z)/dlog10(1.41d0)) +
+ alpha60 * (dlog(1.41d0/(1.d0+z))/dlog10(1.41d0))
endif
return
end
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
subroutine background(zmin,zmax,q0,h0,smin,smax,bg)
c Calculates the background intensity (Iv) at 60microns
c for observed fluxes between smin and smax (in Jy).
c Returns the intensity in units of W/Hz/m^2/sr. (SI units)
implicit real*8 (a­h,o­z)
pi = 4.d0 * datan(1.d0)
n = 1000
ds=dlog(smax/smin)/dfloat(n)
sum=0.d0
do i=1,n
s=smin * dexp( (dfloat(i)­0.5d0)*ds)
call sourcecounts(zmin,zmax,q0,h0,s,c)
sum = sum + ds*s*s*c
enddo
c convert from Jy/sr to W/Hz/m^2/sr
bg = sum * 4.d0 * pi * 1.d­26
return
end

Appendix C
Cascade program
This is the program used to calculate the mono­energetic spectrum for the injection energy
given (in this case 10 9 GeV).
PROGRAM CASCADE_MBR
c 6/11/96
c Modelling an electromagnetic cascade through the microwave background
IMPLICIT REAL*8 (A­H,O­Z)
PARAMETER(NCASC=250)
DATA EMASS/0.5110034D­3/ ! ELECTRON MASS (GEV)
DATA BK/8.617E­14/ ! BOLTZMANN CONST (GEV/K)
DATA PC6/3.086D24/ ! 1 MPC IN CM
DATA SIGMAT/0.665D­24/ ! THOMSON'S CROSS­SECTION (CM^2)
T = 2.728D0 ! TEMP OF CMB BLACKBODY
PI = DACOS(­1.D0)
STHR_PP = (2.D0*EMASS)**2 ! THRESHOLD FOR PAIR PRODUCTION
EMIN_EL = 1.D4 ! MINIMUM ELECTRON ENERGY IS 10^5 GEV
DISTANCE = 50.d­3 ! DISTANCE TO SOURCE IN MPC
E0 = 10.d0**(9.0) ! INJECTION GAMMA RAY ENERGY
DO NUMBER=1,NCASC
E = E0
X=0.d0
ITP=0
NSTORE=0
CALL STORE(X,ITP,E,NSTORE,1)
1 IF(NSTORE.LE.0)GOTO 5
CALL STORE(X,ITP,E,NSTORE,­1)
2 IF(ITP.EQ.0)THEN
CALL PPMFP(E,GGMFP,EBAR)
DX = SAMPLE_INTLENGTH(ggmfp)
X=X+DX
IF (X.LT.DISTANCE) THEN
56

APPENDIX C. CASCADE PROGRAM 57
SMAX = 4.D0*E*EBAR
IF (SMAX.GT.STHR_PP) THEN
S = SSAMP_PP(SMAX)
if (s.gt.smax) then
print *,'s > smax pp',s,smax,S/(2.D0*E*EBAR)
endif
THETA = DACOS(1.D0 ­ S/(2.D0*E*EBAR))
CALL TRANSFORM_PP(E,EBAR,THETA,E1,E2)
if (e1.gt.e0) print *,'e1>e0',e1,e,ebar,theta
if (e2.gt.e0) print *,'e2>e0',e2,e,ebar,theta
CALL STORE(X,1,E1,NSTORE,1)
CALL STORE(X,1,E2,NSTORE,1)
ELSE
if (e.gt.e0) print *,e,x
CALL RECORD(E)
ENDIF
ELSE
if (e.gt.e0) print *,e,x
CALL RECORD(E)
ENDIF
GOTO 1
ELSE
CALL ELMFP(E,EMFP,EBAR)
DX = SAMPLE_INTLENGTH(EMFP)
X=X+DX
IF (X.LT.DISTANCE) THEN
ONEMBETA = 0.5D0*(EMASS/E)**2 + 0.125*(EMASS/E)**4
YMAX = 2.D0*E*EBAR*(2­ONEMBETA)
YMIN = 2.D0*E*EBAR*ONEMBETA
Y = XSAMP_IC(YMIN,YMAX)
if (Y.gt.Ymax) print *,'s>smax ic',Y,Ymax
if (Y.lt.0) print *,'s COSTHETA = (1.D0 ­ Y/(2.D0*E*EBAR))
if (abs(costheta).gt.1.d0) print *,'costheta',costheta,s
THETA = DACOS(COSTHETA)
CALL TRANSFORM_IC(E,EBAR,THETA,EPH)
if (eph.gt.e0) then
print *,'eph>e0',eph,theta,ebar
endif
E = E + EBAR ­ EPH
CALL STORE(X,0,EPH,NSTORE,1)
IF (E.LT.EMIN_EL) THEN
GOTO 1
ENDIF
ELSE
GOTO 1
ENDIF
ENDIF

APPENDIX C. CASCADE PROGRAM 58
GOTO 2
5 ENDDO
CALL OUTPUT(NCASC)
STOP
END
C=================================================================
SUBROUTINE PPMFP(E,GGMFP,EBAR)
IMPLICIT REAL*8(A­H,O­Z)
DIMENSION ENERGY(200),XPP(200),EBARPP(200)
DATA IFIRST/1/
IF (IFIRST.EQ.1) THEN
ifirst = 0
open (unit=12,file='mbr_pp.dat',status='old')
read(12,*) nsize
read(12,*) (energy(i),i=1,nsize)
read(12,*) (xpp(i),i=1,nsize)
read(12,*) (ebarpp(i),i=1,nsize)
close(unit=12)
endif
GGMFP = 0.D0
EBAR = 0.D0
DELTA_EL = ENERGY(2)­ENERGY(1)
EL = DLOG10(E)
IE = (EL­ENERGY(1))/DELTA_EL + 1
IF (IE.GT.NSIZE) IE=NSIZE­1
IF (IE.LT.1) IE=1
c interpolate in arrays
SLOPE1 = (XPP(IE+1)­XPP(IE))/DELTA_EL
GGMFP = 10.D0**(XPP(IE) + SLOPE1*(EL­ENERGY(IE)))
SLOPE2 = (EBARPP(IE+1)­EBARPP(IE))/DELTA_EL
EBAR = 10.D0**(EBARPP(IE) + SLOPE2*(EL­ENERGY(IE)))
RETURN
END
C=================================================================
SUBROUTINE ELMFP(E,EMFP,EBAR)
IMPLICIT REAL*8(A­H,O­Z)
DIMENSION ENERGY(200),XIC(200),EBARIC(200)
DATA IFIRST/1/

APPENDIX C. CASCADE PROGRAM 59
IF (IFIRST.EQ.1) THEN
ifirst = 0
open (unit=12,file='mbr_ic.dat',status='old')
read(12,*) nsize
read(12,*) (energy(i),i=1,nsize)
read(12,*) (xic(i),i=1,nsize)
read(12,*) (ebaric(i),i=1,nsize)
close(unit=12)
endif
EMFP = 0.D0
EBAR = 0.D0
DELTA_EL = ENERGY(2)­ENERGY(1)
EL = DLOG10(E)
IE = (EL­ENERGY(1))/DELTA_EL + 1
IF (IE.GT.NSIZE) IE=NSIZE­1
IF (IE.LT.1) IE=1
c interpolate in arrays
SLOPE1 = (XIC(IE+1)­XIC(IE))/DELTA_EL
EMFP = 10.D0**(XIC(IE) + SLOPE1*(EL­ENERGY(IE)))
SLOPE2 = (EBARIC(IE+1)­EBARIC(IE))/DELTA_EL
EBAR = 10.D0**(EBARIC(IE) + SLOPE2*(EL­ENERGY(IE)))
RETURN
END
C==================================================================
SUBROUTINE STORE(X,ITP,E,NSTORE,IDEPOSIT)
c from R.J.Protheroe's electron_photon_cascade.f
PARAMETER (NMAX=1000)
C NMAX IS THE MAXIMUM NUMBER OF PARTICLES THAT CAN BE STORED AT
C ANY ONE TIME IN THE STORE
IMPLICIT REAL*8(A­H,O­Z)
DIMENSION XS(NMAX),ITPS(NMAX),ES(NMAX)
IF(NSTORE.EQ.0.AND.IDEPOSIT.NE.1)
& WRITE(6,*)'CANNOT WITHDRAW FROM STORE'
IF(NSTORE.EQ.NMAX.AND.IDEPOSIT.EQ.1)
& WRITE(6,*)'STORE FULL'
IF(IDEPOSIT.EQ.1)THEN
NSTORE=NSTORE+1
XS(NSTORE)=X
ITPS(NSTORE)=ITP
ES(NSTORE)=E

APPENDIX C. CASCADE PROGRAM 60
ELSE
X=XS(NSTORE)
ITP=ITPS(NSTORE)
E=ES(NSTORE)
NSTORE=NSTORE­1
ENDIF
RETURN
END
*******************************************************************
SUBROUTINE RECORD(E,ifirst)
IMPLICIT REAL*8 (A­H,O­Z)
PARAMETER(NBIN=62)
DIMENSION NUMBER(NBIN),Esqdnl(NBIN),EPT(NBIN)
DATA NUMBER/NBIN*0/
EMINL = 3.95D0
DELTAE = 0.1D0
IE = (DLOG10(E)­EMINL)/DELTAE + 1
IF ((IE.LE.NBIN).AND.(IE.GT.0)) THEN
NUMBER(IE) = NUMBER(IE) + 1
ELSE
PRINT *,'PHOTON OUTSIDE BINNED AREA.', E
ENDIF
RETURN
ENTRY OUTPUT(NCASC)
DO K=1,NBIN
EPT(K) = EMINL + (K­0.5d0)*DELTAE
WIDTH=10.D0**EPT(K)*(10.D0**0.05 ­ 10.D0**(­0.05))
EDIST = DFLOAT(NUMBER(K))/(WIDTH*DFLOAT(NCASC))
IF (EDIST.NE.0.D0) THEN
ESQDNL(K)=2.D0*EPT(K)+DLOG10(EDIST)
ELSE
ESQDNL(K)=­100.D0
ENDIF
ENDDO
OPEN (UNIT=9,FILE='spec­9.0­50k­4.dat',STATUS='unknown')
write(9,*) NBIN
WRITE(9,*) (EPT(K),K=1,NBIN)
WRITE(9,*) (ESQDNL(K),K=1,NBIN)
CLOSE(UNIT=9)

APPENDIX C. CASCADE PROGRAM 61
RETURN
END
*******************************************************************
FUNCTION SAMPLE_INTLENGTH(AMFP)
IMPLICIT NONE
c Samples an interaction length from an exponential distribution
c given a mean interaction length of AMFP.
c Both GGMFP and GAMINTLENGTH are in Mpc.
REAL *8 SAMPLE_INTLENGTH,AMFP
REAL R,RNUNF
SAMPLE_INTLENGTH = 0.D0
R = RNUNF()
R = RNUNF()
SAMPLE_INTLENGTH = ­AMFP*DLOG(1.d0­R)
RETURN
END
C==================================================================
FUNCTION SSAMP_PP(SMAX)
C
C CALCULATES INTEGRAL OF S*SIGMA(S) OVER S FROM STHR TO SMAX
C FOR PAIR PHOTOPRODUCTION.
C SMAX IN GEV^2
C UNIT: GEV^4 CM^2
C
IMPLICIT REAL*8 (A­H,O­Q,S­Z)
IMPLICIT REAL (R)
DIMENSION sintl(400),sarrayl(400)
DATA IFIRST/1/
DATA EMASS/0.5110034D­3/
DATA DELTA/0.0001D0/
data de/0.025d0/
STHR=(2*EMASS)**2
IF(IFIRST.EQ.1)THEN
IFIRST=0
SUM=0.D0
VAL=0.D0
DLOGS=DELTA*2.3025851D0
DO ISMAX=1,100000
OLDVAL=VAL
S=STHR*10.D0**(ISMAX*DELTA)

APPENDIX C. CASCADE PROGRAM 62
VAL=S*(GGEETOT_CS(S))*S*DLOGS
SUM=SUM+0.5D0*(OLDVAL+VAL)
if (mod(ismax,250).eq.0) then
sarrayl(ismax/250) = dlog10(s)
sintl(ismax/250) = dlog10(sum)
endif
ENDDO
ENDIF
SSAMP_PP=0.D0
GGEE_SIGINT=0.D0
IF(SMAX.LE.STHR)RETURN
SMAXL=DLOG10(SMAX/STHR)
ISMAX=SMAXL/DE+1.D0
IF(ISMAX.GT.100000) THEN
WRITE(6,*)'GGEE_SIGINT: (ISMAX.GT.100000). ',SMAX,ISMAX
ELSEIF(ISMAX.EQ.1)THEN
GGEE_SIGINT=(SMAXL/DE+1.D0­ISMAX)*SINTl(ISMAX)
ELSE
GGEE_SIGINT=SINTl(ISMAX­1)+(SMAXL/DE+1.D0­ISMAX)*
& (SINTl(ISMAX)­Sintl(ISMAX­1))
ENDIF
GGEE_SIGINT=10.d0**GGEE_SIGINT
20 R = RNUNF()
R = RNUNF()
VAL = dlog10(GGEE_SIGINT*R)
SSAMP_PP=0.D0
k=0
10 k=k+1
IF (SINTl(k).LT.VAL) GOTO 10
if (k.eq.1) k=2
SLOPE = (SINTl(k) ­ SINTl(k­1))/DE
if (slope.eq.0) print *,'GGEE_SINT: SLOPE=0',smax,k
SSAMPLOG = sarrayl(k­1) + (VAL­SINTl(k­1))/SLOPE
SSAMP_PP = 10.d0**SSAMPLOG
IF (SSAMP_PP.LE.STHR)then
c write(6,*) 'SSAMP_PP: (S.LT.STHR).',SSAMP_PP,STHR
GOTO 20
ENDIF
RETURN
END
C =================================================================
FUNCTION GGEETOT_CS(S)
IMPLICIT REAL*8 (A­H,O­Z)

APPENDIX C. CASCADE PROGRAM 63
DATA EMASS/0.5110034d­3/
GGEETOT_CS=0.D0
IF(S.LE.(2.D0*EMASS)**2)RETURN
c B=DSQRT(1.D0­(2.D0*EMASS)**2/S)
x = (2.D0*EMASS)**2/S !1 ­ beta^2
y = 0.5d0*x + 0.125*x**2 !1 ­ beta
SIG=­2.D0*(1.d0­y)*(x+1.D0)+(2.d0+x*(2.d0­x))*DLOG((2.d0­y)/y)
GGEETOT_CS=1.2474D­25*x*SIG
return
end
*******************************************************************
SUBROUTINE TRANSFORM_PP(EPH1,EPH2,THETA,EE1,EE2)
c transforms to electron rest frame, performs pair production
c and transforms back to lab frame.
c EPH1,EPH2 = energies of photons (hard&soft) in lab frame
c THETA = angle between photons' directions in lab frame.
c EE1, EE2 = energies of electrons in lab frame after interaction
IMPLICIT REAL*8 (A­H,O­Q,S­Z)
IMPLICIT REAL (R)
DATA EMASS/0.51100337D­3/ ! ELECTRON MASS (GEV)
PI = DACOS(­1.D0)
c Calculate energy of each photon in CM frame
S = 2.D0*EPH1*EPH2*(1­COS(THETA))
EPRIME = DSQRT(S)/2.D0
IF (EPRIME.LT.EMASS) THEN
EE1 = 0.D0
EE2 = 0.D0
PRINT *,'EPRIME.LT.EMASS'
RETURN
ELSE
c
c Now calculate geometry of boosted frame
c First angle of boost to hard­photon direction
tanPHI = EPH2*SIN(THETA)/(EPH1+EPH2*COS(THETA))
phi = datan(tanphi)
cosphi = cos(phi)
sinphi = sin(phi)
c Next the beta factor of the boost (velocity)
c ONEMINUSBB = (EPH1*(1­COS(PHI)) + EPH2*(1­COS(THETA­PHI)))
c & /(EPH1 + EPH2)
c BB = 1.D0 ­ ONEMINUSBB

APPENDIX C. CASCADE PROGRAM 64
oneminusbb = (eph1*(1.d0­cosphi)+eph2*(1.d0­cosphi*cos(theta)
& ­sinphi*sin(theta)))/(eph1+eph2)
c
c Gamma factor of the boost
GB = 1.D0/DSQRT(2.D0*ONEMINUSBB ­ ONEMINUSBB**2)
c
c Angle between boost direction and photon direction in CM frame
ALPHA = DATAN(SIN(PHI)/(GB*(COS(PHI)­1.d0+oneminusBB)))
c Now beta and gamma for the electrons
GE = EPRIME/EMASS
onembe = 0.5/ge**2 + 0.125/ge**4
c BE = DSQRT(1.D0 ­ 1.D0/GE**2)
c
c now find momentum (in CM frame) and energy (LAB frame) of electrons
c PC = GE*EMASS*BE !GeV
pc = ge*emass*(1.d0­onembe)
SCAT_ANG = THETAM_PP(EPRIME)
c scat_ang = 1.d0
R=RNUNF()
R=RNUNF()
AZIMUTH = 2.D0*PI*R
PXPRIME = COS(SCAT_ANG)*COS(ALPHA) ­ SIN(SCAT_ANG)*
& COS(AZIMUTH)*SIN(ALPHA)
EE1 = GB*(EPRIME + (1.d0­oneminusBB)*PC*PXPRIME)
EE2 = GB*(EPRIME ­ (1.d0­oneminusBB)*PC*PXPRIME)
RETURN
ENDIF
END
********************************************************************
FUNCTION THETAM_PP(EPH)
IMPLICIT REAL*8(A­H,O­Z)
DIMENSION ENERGY(200),THETAM(200)
DATA IFIRST/1/
IF (IFIRST.EQ.1) THEN
ifirst = 0
open (unit=12,file='thetam_pp.dat',status='old')
read(12,*) nsize
read(12,*) (energy(i),i=1,nsize)
read(12,*) (thetam(i),i=1,nsize)
close(unit=12)
endif
THETAM_PP=0.D0

APPENDIX C. CASCADE PROGRAM 65
DELTA_EL = ENERGY(2)­ENERGY(1)
EL = DLOG10(EPH)
IE = (EL­ENERGY(1))/DELTA_EL + 1
IF (IE.GT.NSIZE) IE=NSIZE­1
IF (IE.LT.1) IE=1
c interpolate in array
SLOPE = (THETAM(IE+1)­THETAM(IE))/DELTA_EL
THETAM_PP = THETAM(IE) + SLOPE*(EL­ENERGY(IE))
RETURN
END
C===================================================================
FUNCTION XSAMP_IC(XMIN,XMAX)
C
C CALCULATES INTEGRAL OF S*SIGMA(S) OVER S FROM XMIN TO XMAX
C FOR INVERSE COMPTON SCATTERING
C X = S ­ EMASS**2
C XMAX,XMIN IN GEV^2
C UNIT: GEV^4 CM^2
C
IMPLICIT REAL*8 (A­H,O­Q,S­Z)
IMPLICIT REAL (R)
PARAMETER(IMAX=500000)
DIMENSION xarrayl(500),xintl(500)
DATA IFIRST/1/
DATA EMASS/0.5110034D­3/
DATA DELTA/0.0001D0/
DATA DE/0.1D0/
IF(IFIRST.EQ.1)THEN
IFIRST=0
XMINLOG = ­45.D0
SUM=0.D0
VAL=0.D0
DLOGX=DELTA*2.3025851D0
DO IXMAX=1,IMAX
OLDVAL=VAL
X=10.D0**(XMINLOG+IXMAX*DELTA)
VAL=X*COMPTON_CS(X)*X*DLOGX
SUM=SUM+0.5D0*(OLDVAL+VAL)
if (mod(ixmax,1000).eq.0) then
xarrayl(ixmax/1000) = dlog10(x)
xintl(ixmax/1000) = dlog10(sum)
endif
ENDDO

APPENDIX C. CASCADE PROGRAM 66
ENDIF
XSAMP_IC = 0.D0
COMPT_SIGINT=0.D0
IF(XMAX.LE.0.D0)RETURN
XMAXL= DLOG10(XMAX) ­ XMINLOG
IXMAX=XMAXL/DE+1.D0
IF (IXMAX.EQ.0) PRINT *,'COMPTON_SINT: (IXMAX.EQ.0)',xmaxl,xmax
IF(IXMAX.GT.IMAX) THEN
WRITE(6,*)'COMPTON_SINT: (IXMAX.GT.IMAX). ',SMAX,IXMAX
ELSEIF(IXMAX.EQ.1)THEN
UPPER_INT=(XMAXL/DE+1.D0­IXMAX)*XINTL(IXMAX)
ELSE
UPPER_INT=XINTL(IXMAX­1)+(XMAXL/DE+1.D0­IXMAX)*
& (XINTL(IXMAX)­XINTL(IXMAX­1))
ENDIF
IF(XMIN.LE.0.D0)RETURN
XMINL= DLOG10(XMIN) ­ XMINLOG
IXMIN=XMINL/DE+1.D0
IF(IXMIN.GT.IMAX) THEN
WRITE(6,*)'COMPTON_SINT: (IXMIN.GT.IMAX). ',XMIN,IXMIN
ELSEIF (IXMIN.LT.0) THEN
WRITE(6,*)'COMPTON_SINT: (IXMIN.LT.0). ',XMIN,IXMIN
ELSEIF(IXMIN.EQ.1)THEN
BOTTOM_INT=(XMINL/DE+1.D0­IXMIN)*XINTL(IXMIN)
ELSE
BOTTOM_INT=XINTL(IXMIN­1)+(XMINL/DE+1.D0­IXMIN)*
& (XINTL(IXMIN)­XINTL(IXMIN­1))
ENDIF
COMPT_SIGINT = 10.D0**UPPER_INT ­ 10.D0**BOTTOM_INT
XSAMP_IC = 0.D0
20 R = RNUNF()
R = RNUNF()
VAL = DLOG10(COMPT_SIGINT*R)
J=0
10 J=J+1
IF (XINTL(J).LT.VAL) GOTO 10
if (j.eq.1) j=2
SLOPE = (XINTL(J)­XINTL(J­1))/DE
XSAMPLOG = XARRAYL(J­1) + (VAL­XINTL(J­1))/SLOPE
XSAMP_IC = 10.D0**XSAMPLOG
IF (XSAMP_IC.LT.0.D0) THEN
PRINT *,'COMPTON_SINT: (XSAMPLE.LT.0.D0)',R,XSAMP_IC
GOTO 20
ENDIF
RETURN
END

APPENDIX C. CASCADE PROGRAM 67
C ==================================================================
FUNCTION COMPTON_CS(X)
c cross section for inverse compton scattering, where S is the
c CM frame energy squared.
c returns COMPTON_CS in cm^2
IMPLICIT REAL*8 (A­H,O­Z)
DATA R/2.818d­13/ !classical radius of electron (cm)
DATA EMASS/0.51100337d­3/ !mass of electron (GeV)
PI = DACOS(­1.D0)
c calculate photon energy in electron's rest frame
W = X/(2.D0*EMASS)
W = W/EMASS !need W in units of EMASS
if (w.le.0.d0) print *, 'compton_cs:w<0',w
if (w.gt.1.d­3) then
SIGMA = (1.D0+W)*(2.D0*W*(1.D0+W)/(1.D0+2.D0*W) ­
& DLOG(1.D0+2.D0*W))/W**3
SIGMA = SIGMA + DLOG(1.D0+2.D0*W)/(2.D0*W) ­ (1.D0+3.D0*W)/
& (1.D0+2.D0*W)**2
COMPTON_CS = SIGMA * R**2 * 2.D0 * PI
else
compton_cs = 8.d0*pi*r**2/3.d0
endif
RETURN
END
*******************************************************************
SUBROUTINE TRANSFORM_IC(EE,EPH,THETA,EPHNEW)
c transforms to electron rest frame, performs compton scattering
c interaction and transforms back to lab frame.
c EE = electron energy (GeV), EPH = photon energy (GeV)
c THETA = angle between particles' directions ­ all in lab frame.
c EPHNEW = energy of photon in lab frame after interaction
IMPLICIT REAL*8 (A­H,O­Q,S­Z)
IMPLICIT REAL (R)
DATA EMASS/0.51100337D­3/ ! ELECTRON MASS (GEV)
PI = DACOS(­1.D0)
G = EE/EMASS !gamma factor of electron
B = 1.D0 ­ 0.5D0/G**2 ­ 0.125D0/G**2
EPH_ERF = G*EPH*(1.D0­B*COS(THETA))
if(eph_erf.le.0) print *,'transform_ic: eph_erf.le.0',eph_erf
COSTHETA_ERF = (COS(THETA)­B)/(1.D0­B*COS(THETA))
SINTHETA_ERF = SIN(THETA)/(G*(1.D0­B*COS(THETA)))
SCAT_ANG = THETAM_IC(EPH_ERF) !sample the scattering angle
c scat_ang = 1.d0
EPRIME_ERF = EPH_ERF/(1.D0+EPH_ERF*(1.D0­COS(SCAT_ANG))/EMASS)

APPENDIX C. CASCADE PROGRAM 68
R=RNUNF()
R=RNUNF()
PHI = R*2.D0*PI
PXPRIME = COSTHETA_ERF*COS(SCAT_ANG)­SINTHETA_ERF*SIN(SCAT_ANG)
& *SIN(PHI)
EPHNEW = G*EPRIME_ERF*(1.D0+B*PXPRIME)
RETURN
END
********************************************************************
FUNCTION THETAM_IC(EPH)
IMPLICIT REAL*8(A­H,O­Z)
DIMENSION ENERGY(500),THETAM(500)
DATA IFIRST/1/
IF (IFIRST.EQ.1) THEN
ifirst = 0
open (unit=12,file='thetam_ic.dat',status='old')
read(12,*) nsize
read(12,*) (energy(i),i=1,nsize)
read(12,*) (thetam(i),i=1,nsize)
close(unit=12)
endif
THETAM_IC=0.D0
c interpolate in array
DELTA_EL = ENERGY(2)­ENERGY(1)
EL = DLOG10(EPH)
IE = (EL­ENERGY(1))/DELTA_EL + 1
IF (IE.GT.NSIZE) IE=NSIZE­1
IF (IE.LT.1) IE=1
SLOPE = (THETAM(IE+1)­THETAM(IE))/DELTA_EL
THETAM_IC = THETAM(IE) + SLOPE*(EL­ENERGY(IE))
RETURN
END

Appendix D
Mean free path and mean photon
energy program
PROGRAM GGEE_BB
C MEAN FREE PATH (MPC) FOR PAIR PRODUCTION ON BLACK BODY
C RADIATION OF TEMPERATURE T, AND MEAN SOFT PHOTON ENERGY (GEV)
C ADAPTED FROM A PROGRAM BY R.J. PROTHEROE
IMPLICIT REAL*8 (A­H,O­Z)
PARAMETER (M=150)
DIMENSION EMID(M),XINT(M),EBAR(M)
DATA EMASS/0.5110034D­3/ ! ELECTRON MASS (GEV)
DATA BK/8.617E­14/ ! BOLTZMANN CONST (GEV/K)
DATA PC6/3.086D24/ ! 1 MPC IN CM
T=2.728D0
BKT=BK*T
DO IP=1,M
EMID(IP)=1.D0*10.D0**((IP­0.5D0)*0.1D0)
ENDDO
DO IP=1,M
write(6,*)ip
EM=EMID(IP)
CALL PPMFP(EM,T,GGMFP,E)
XINT(IP)=GGMFP/PC6
ebar(ip)=e
ENDDO
k=0
10 k=k+1
if (ebar(k).eq.0.d0) goto 10
k=k­1
OPEN (UNIT=1, FILE='mbr_MFP.dat', STATUS= 'unknown')
write(1,*) M
WRITE(1,*) (DLOG10(EMID(I)), I=1,M)
WRITE(1,*) (DLOG10(XINT(I)), I=1,M)
close(unit=1)
open (unit=2,file='mbr_ebarpp.dat',status='unknown')
69

APPENDIX D. MEAN FREE PATH AND MEAN PHOTON ENERGY PROGRAM70
write(2,*) (m­k)
write(2,*) (dlog10(emid(i)), i=k+1,m)
WRITE(2,*) (dlog10(EBAR(I)),I=k+1,M)
CLOSE (UNIT=2)
STOP
END
C==================================================================
SUBROUTINE GGEE_SINT(SMAX,GGEE_SIGINT)
C
C CALCULATES INTEGRAL OF S*SIGMA(S) OVER S FROM STHR TO SMAX
C FOR PAIR PHOTOPRODUCTION.
C SMAX IN GEV^2
C UNIT: GEV^4 CM^2
C
IMPLICIT REAL*8 (A­H,O­Q,S­Z)
IMPLICIT REAL (R)
DIMENSION SIGSINT(1000000)
DATA IFIRST/1/
DATA EMASS/0.5110034D­3/
DATA DELTA/0.0001D0/
STHR=(2*EMASS)**2
IF(IFIRST.EQ.1)THEN
IFIRST=0
SUM=0.D0
VAL=0.D0
DLOGS=DELTA*2.3025851D0
DO ISMAX=1,260000
OLDVAL=VAL
S=STHR*10.D0**(ISMAX*DELTA)
VAL=S*(GGEETOT_CS(S))*S*DLOGS
SUM=SUM+0.5D0*(OLDVAL+VAL)
SIGSINT(ISMAX)=SUM
ENDDO
ENDIF
GGEE_SIGINT=0.D0
IF(SMAX.LE.STHR)RETURN
SMAXL=DLOG10(SMAX/STHR)
ISMAX=SMAXL/DELTA+1.D0
IF(ISMAX.GT.260000) THEN
WRITE(6,*)'GGEE_SIGINT(SMAX):
& (ISMAX.GT.150000). ',SMAX,ISMAX
ELSEIF(ISMAX.EQ.1)THEN

APPENDIX D. MEAN FREE PATH AND MEAN PHOTON ENERGY PROGRAM71
GGEE_SIGINT=(SMAXL/DELTA+1.D0­ISMAX)*SIGSINT(ISMAX)
ELSE
GGEE_SIGINT=SIGSINT(ISMAX­1)+(SMAXL/DELTA+1.D0­ISMAX)*
& (SIGSINT(ISMAX)­SIGSINT(ISMAX­1))
ENDIF
GGEE_SIGINT=GGEE_SIGINT
RETURN
END
C =================================================================
SUBROUTINE PPMFP(EP,T,GGMFP,EBAR)
C
C INTERACTION LENGTH (CM) OF PHOTONS IN BLACK BODY RADIATION
C EP (GEV), T (K)
IMPLICIT REAL*8 (A­H,O­Z)
DATA BK/8.617E­14/,PI/3.14159D0/
DATA EMASS/0.5110034d­3/
BKT=BK*T
DLOGEPH=0.01D0*2.3025851D0
BETA=1.D0
STHR=(2.d0*EMASS)**2
EPHTHR=STHR/(2.D0*EP*(1.D0+BETA))
SUM=0.D0
AVSUM=0.D0
DO IEPH=1,1000
EPH=EPHTHR*10.D0**(0.01D0*IEPH)
SMAX=2.D0*EPH*EP*(1.D0+BETA)
CALL GGEE_SINT(SMAX,SINT,0,DUMMY)
if (sint.le.0.d0) print *,'GGEE_SINT.LE.0.D0',sint,eph
VAL = BB_DENS(EPH,BKT)*SINT*DLOGEPH/EPH
SUM=SUM+VAL
AVSUM=AVSUM+VAL*EPH
ENDDO
GGMFP=(8.D0*EP**2*BETA)*1.D80
EBAR = AVSUM * 1.D80
if (sum.lt.0.d0) print *,'GGEE_SINT:sum<0',sum,ep
IF(SUM.GT.1.D­80) THEN
GGMFP=(8.D0*EP**2*BETA)/SUM
EBAR = AVSUM/SUM
ENDIF

APPENDIX D. MEAN FREE PATH AND MEAN PHOTON ENERGY PROGRAM72
RETURN
END
C ==================================================================
FUNCTION GGEETOT_CS(S)
IMPLICIT REAL*8 (A­H,O­Z)
DATA EMASS/0.5110034d­3/
GGEETOT_CS=0.D0
IF(S.LE.(2.D0*EMASS)**2)RETURN
x = (2.D0*EMASS)**2/S !1 ­ beta^2
y = 0.5d0*x + 0.125*x**2 !1 ­ beta
SIG=­2.D0*(1.d0­y)*(x+1.D0)+(2.d0+x*(2.d0­x))*DLOG((2.d0­y)/y)
GGEETOT_CS=1.2474D­25*x*SIG
return
end
C ==============================================================
DOUBLE PRECISION FUNCTION BB_DENS(E,KT)
C
C GIVES DIFFERENTIAL NUMBER DENSITY OF PHOTONS IN A
C BLACK BODY FIELD OF TEMPERATURE T (KELVIN)
C
C UNITS PHOTONS CM^­3 GEV^­1
C
C E: PHOTON ENERGY (GEV)
C KT: (BOLTZMANN CONST)*T (GEV)
C
IMPLICIT NONE
DOUBLE PRECISION EKT,KT,E,FACTOR
EKT=E/KT
IF(EKT.LT.1.D­7)THEN
FACTOR=1/(EKT+EKT**2/2+EKT**3/6)
ELSE
IF(EKT.LT.1.D2)THEN
FACTOR=1.D0/(DEXP(EKT)­1.D0)
ELSE
FACTOR=0.D0
ENDIF
ENDIF

APPENDIX D. MEAN FREE PATH AND MEAN PHOTON ENERGY PROGRAM73
BB_DENS=1.31D40*E**2*FACTOR
RETURN
END

Appendix E
Median scattering angle program
PROGRAM THETA_SAMPLE_PP2
C sample the mean scattering angle theta (called alpha in the
c thesis) in the centre of momentum frame for a
C photon­photon pair production interaction
IMPLICIT REAL*8 (A­H,O­Z)
PARAMETER(N=140)
DIMENSION ENERGY(N),THETAM(N)
DATA EMASS/0.51100337D­3/ ! ELECTRON MASS (GEV)
DE = 0.05
EMINLOG = DLOG10(EMASS)
DO I=1,N
print *,i
ENERGY(I) = 10**(EMINLOG+DFLOAT(I)*DE)
THETAM(I) = THETAM_PP(ENERGY(I))
ENDDO
OPEN(UNIT=9,FILE='thetam_pp.dat',STATUS='unknown')
WRITE(9,*) N
WRITE(9,*) (DLOG10(ENERGY(I)), I=1,N)
WRITE(9,*) (THETAM(I), I=1,N)
CLOSE(UNIT=9)
END
********************************************************************
FUNCTION THETAM_PP(EPH)
IMPLICIT REAL*8 (A­H,O­Z)
PARAMETER(N=100000)
DATA EMASS/0.51100337D­3/ ! ELECTRON MASS (GEV)
DATA R/2.818d­13/ !classical radius of electron (cm)
PI = DACOS(­1.D0)
WM = EPH/EMASS
IF (WM.LT.1) THEN
THETAM_PP = 0.D0
74

APPENDIX E. MEDIAN SCATTERING ANGLE PROGRAM 75
PRINT *,'THETAM_PP : EPH < EMASS'
RETURN
ELSE
B = 1.D0 ­ 0.5D0/WM**2 !beta
bsqd = 1.d0 ­ 1.d0/wm**2 !beta ^2
IF (WM.LT.620) THEN
SUM = 0.D0
VAL = 0.D0
DX = 1.D0/N
X = 0.D0 ­ DX
DO I=1,N+1
X = X + DX
VAL = 1.D0­Bsqd**2*X**4+2.D0*Bsqd/WM**2*(1­X**2)
VAL = 0.5D0*PI*R**2 * B/(WM**2)*VAL/
& (1.D0­Bsqd**2*X**2)**2
SUM = SUM + VAL*DX
ENDDO
if (sum.le.0.d0) print *,'THETAM_PP: SUM<0'
HALFSUM = dlog10(sum/2.D0)
THETAM_PP = 0.D0
SUM = 0.D0
VAL = 0.D0
X = 0.D0 ­ DX
30 CONTINUE
if ((x+dx).gt.1.d0) then
print *,'THETAM_PP:reached limit of x'
goto 20
endif
X = X + DX
VAL = 1.D0­Bsqd**2*X**4+2.D0*Bsqd/WM**2*(1­X**2)
VAL = 0.5D0*PI*R**2 * B/(WM**2)*VAL/
& (1.D0­Bsqd**2*X**2)**2
if (val.le.0.d0) print *,'THETAM_PP: VAL<0',val,x
SUM = SUM + VAL*DX
TOTALLOG = DLOG10(SUM)
IF (TOTALLOG.GE.HALFSUM)THEN
IF (ABS(TOTALLOG­HALFSUM).LT.0.0001) THEN
GOTO 20
ELSE
SUM = SUM ­ VAL*DX
X=X­DX
DX = DX/10.D0
GOTO 30
ENDIF
ENDIF
GOTO 30
20 THETAM_PP = DACOS(X)
ELSE

APPENDIX E. MEDIAN SCATTERING ANGLE PROGRAM 76
TOLX = 1.D­6
X = 0.5D0
FOLD = DLOG(WM*(1.D0­Bsqd*X*X)) ­ B*(X­0.5D0)
DX = 0.01D0
40 CONTINUE
X = X+DX
FNEW = DLOG(WM*(1.D0­Bsqd*X*X)) ­ B*(X­0.5D0)
IF ((FOLD*FNEW).LT.0.D0) THEN
X=X­DX
DX=DX/2.D0
ENDIF
IF (ABS(DX).GT.TOLX) GOTO 40
THETAM_PP = DACOS(X)
ENDIF
RETURN
ENDIF
END