UNDERSTANDING NATURE’S PARTICLE ACCELERATORS USING HIGH ENERGY
GAMMA-RAY SURVEY INSTRUMENTS
By
Anushka Udara Abeysekara

A DISSERTATION
Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of
Physics - Doctor of Philosophy
2014

ABSTRACT
UNDERSTANDING NATURE’S PARTICLE ACCELERATORS USING
HIGH ENERGY GAMMA-RAY SURVEY INSTRUMENTS
By
Anushka Udara Abeysekara

Nature’s particle accelerators, such as Pulsars, Pulsar Wind Nebulae, Active Galactic
Nuclei and Supernova Remnants accelerate charged particles to very high energies that then
produce high energy photons. The particle acceleration mechanisms and the high energy
photon emission mechanisms are poorly understood phenomena. These mechanisms can
be understood either by studying individual sources in detail or, alternatively, using the
collective properties of a sample of sources. Recent development of GeV survey instruments,
such as Fermi-LAT, and TeV survey instruments, such as Milagro, provides a large sample
of high energy gamma-ray flux measurements from galactic and extra-galactic sources. In
this thesis I provide constraints on GeV and TeV radiation mechanisms using the X-ray-TeV
correlations and GeV-TeV correlations.
My data sample was obtained from three targeted searches for extragalactic sources
and two targeted search for galactic sources, using the existing Milagro sky maps. The first
extragalactic candidate list consists of Fermi-LAT GeV extragalactic sources, and the second
extragalactic candidate list consists of TeVCat extragalactic sources that have been detected
by Imaging Atmospheric Cerenkov Telescopes (IACTs). In both extragalactic candidate
lists Markarian 421 was the only source detected by Milagro. A comparison between the
Markarian 421 time-averaged flux, measured by Milagro, and the flux measurements of
transient states, measured by IACTs, is discussed. The third extragalactic candidate list
is a list of potential TeV emitting BL Lac candidates that was synthesized using X-ray

observations of BL Lac objects and a Synchrotron Self-Compton model. Milagro’s sensitivity
was not sufficient to detect any of those candidates. However, the 95% confidence flux
upper limits of those sources were above the predicted flux. Therefore, these results provide
evidence to conclude that the Synchrotron Self-Compton model for BL Lac objects is still a
viable model. Targeted searches for galactic candidates were able to measure TeV emission
associated with 14 Fermi-LAT GeV pulsars.
In this thesis I also presented a new multi-wavelength technique that I developed to isolate
the flux correlation factor (fΩ ) of pulsars as a function of pulsar spin down luminosity. The
correlation between fΩ and pulsar spin-down luminosity for a Fermi-LAT GeV pulsar sample
was measured using the measurements obtained in the Milagro targeted search performed
for galactic sources and from the literature. The measured correlation has some features
that favor the Outer Gap model over the Polar Cap, Slot Gap and One Pole Caustic models
for pulsar emission in the energy range of 0.1 to 100 GeV. However, these simulated models
failed to explain many other important pulsar population characteristics. Therefore, further
improvements on the galactic pulsar population simulations are needed to provide tighter
constraints.

I dedicate my dissertation work to my family.

iv

ACKNOWLEDGMENTS

First of all, I like to thank my parents, wife and little son for their continuous caring,
support and encouragement. I cannot express enough thanks to my thesis adviser Prof. Jim
Linnemann for offering me the opportunity to work with him and his continuous support
during my days in graduate school.
I also like to extend my deepest gratitude to all of my thesis committee members: Dr.
Gus Sinnis, Prof. Kirsten Tollefson, Prof. Scott Pratt, and Prof. Mark Voit for their advice
and assistance. Their sincere guidance and comments helped shape this thesis in the best
way possible.
All the members of the Milagro and HAWC collaborations deserve special thanks. There
are too many special people in both collaborations to name here; it has been my good fortune
to be a part of such a supportive and intelligent group of people.
Last but not least I want to thank Michigan State University postdoc Dr. Tilan Ukwatta,
graduate student Samuel Marinelli, my office mates, all of my teachers in the past and
Michigan State University staff members who helped me to successfully finish my PhD.

v

TABLE OF CONTENTS

LIST OF TABLES . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . .

ix

LIST OF FIGURES

. . . . . . . . . . . . . . . . .

xi

. . . . . . . . . . . . .

Chapter 1 Gamma-ray Astronomy . . . .
1.1 Introduction to Gamma-ray Astronomy .
1.2 Modern Gamma-ray Detectors . . . . . .
1.2.1 Space-based Detectors . . . . . .
1.2.2 Ground-based Detectors . . . . .
1.3 Imaging Atmospheric (or Air) Cherenkov
1.4 Extensive Air Shower (EAS) Arrays . . .

. . . . . . . . . . . .
. . . . . . . . . . . .
. . . . . . . . . . . .
. . . . . . . . . . . .
. . . . . . . . . . . .
Telescopes (IACTs)
. . . . . . . . . . . .

Chapter 2 Gamma-ray sources . . . . . . . . . .
2.1 Non-thermal Gamma-ray Production Processes
2.1.1 Synchrotron Radiation . . . . . . . . . .
2.1.2 Curvature Radiation . . . . . . . . . . .
2.1.3 Inverse Compton Radiation . . . . . . .
2.2 Pulsars . . . . . . . . . . . . . . . . . . . . . . .
2.2.1 Properties of Neutron Stars . . . . . . .
2.2.2 Pulsar Spin Evolution . . . . . . . . . .
2.2.3 Spin-down Luminosity Evolution . . . .
2.2.4 Characteristic Age of a Pulsar . . . . . .
2.2.5 Pulsar Magnetosphere . . . . . . . . . .
2.3 Pulsar Wind Nebulae . . . . . . . . . . . . . . .

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

1
1
3
4
4
6
7

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

9
9
10
11
12
13
14
15
16
17
18
19

Chapter 3 The Milagro Observatory . . . . . . . . . .
3.1 Milagro Detector . . . . . . . . . . . . . . . . . . . .
3.1.1 The Main Pond . . . . . . . . . . . . . . . . .
3.1.2 The Outrigger Array . . . . . . . . . . . . . .
3.2 Front-end Electronics and Time to Digital Converter
3.3 Trigger . . . . . . . . . . . . . . . . . . . . . . . . . .
3.4 Reconstruction . . . . . . . . . . . . . . . . . . . . .
3.4.1 Shower Core Location Reconstruction . . . . .
3.4.2 Incidence Angle Reconstruction . . . . . . . .
3.5 Monte Carlo Simulations . . . . . . . . . . . . . . . .

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

24
24
25
29
29
32
35
35
36
38

Chapter 4 Milagro Data Analysis Framework . . . . . . . . . . . . . . . . .
4.1 Gamma/Hadron Separation Parameter . . . . . . . . . . . . . . . . . . . . .
4.2 Background Estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

41
41
42

vi

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.

4.3
4.4
4.5

Significance Maps . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Flux calculation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Energy Reconstruction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

45
46
48

Chapter 5 The Fermi Large Area Telescope (Fermi-LAT) GeV Catalogs
5.1 Fermi Large Area Telescope (Fermi-LAT) . . . . . . . . . . . . . . . . . . . .
5.2 Fermi-LAT Second Source Catalog (2FGL) . . . . . . . . . . . . . . . . . . .
5.3 The Second Fermi-LAT Catalog of Gamma-ray Pulsars (2PC) . . . . . . . .
5.4 Source Naming Convention . . . . . . . . . . . . . . . . . . . . . . . . . . . .

52
53
54
57
58

Chapter 6 The High Altitude Water Cherenkov Observatory
6.1 The HAWC Detector . . . . . . . . . . . . . . . . . . . . . . .
6.2 GPS Timing and Control System (GTC System) . . . . . . . .
6.2.1 GPS Timing System Hardware . . . . . . . . . . . . .
6.2.2 GPS Timing System Firmware . . . . . . . . . . . . . .
6.2.3 Control System Hardware . . . . . . . . . . . . . . . .
6.2.4 Control System Firmware . . . . . . . . . . . . . . . .
6.3 Performance of the GTC system . . . . . . . . . . . . . . . . .
6.4 Additional Capabilities of H Clock cards . . . . . . . . . . . .

(HAWC)
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .

.
.
.
.
.
.
.
.
.

60
60
63
63
65
69
72
74
75

Chapter 7 Data Analysis Methodology . .
7.1 False Discovery Rate method . . . . . . .
7.1.1 Calculation . . . . . . . . . . . . .
7.2 Stack Analysis . . . . . . . . . . . . . . . .

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

77
78
79
79

Chapter 8 Milagro Observations of Extragalactic Sources
8.1 Targeted Search for 2FGL Extragalactic Sources . . . . . . .
8.1.1 Comparison with Previous γ-ray Flux Measurements
8.1.2 Undetected Candidates . . . . . . . . . . . . . . . . .
8.2 Targeted Search for TeVCat Extragalactic Sources . . . . . .
8.3 Targeted Search for Potential TeV Emitting BL Lac Objects

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

81
82
83
84
85
87

Chapter 9 The Milagro Targeted Search for GeV Pulsars.
9.1 Comparison With IACT Flux Measurements . . . . . . . . .
9.2 Population Study of Associated TeV emissions . . . . . . . .
9.3 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . .

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

. 94
. 98
. 102
. 107

Chapter 10 Experimental Constraints on Gamma-Ray Pulsar
10.1 A New Method to Isolate fΩ . . . . . . . . . . . . . . . . . . .
10.2 The Sample of Pulsars and Their Associated
PWNe . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
10.3 Results and Discussion . . . . . . . . . . . . . . . . . . . . . .

Gap Models 115
. . . . . . . . 116
. . . . . . . . 119
. . . . . . . . 124

Chapter 11 Milagro Observations of Fermi-LAT Galactic Sources
vii

. . . . . 129

Chapter 12 Conclusions
APPENDIX

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140

BIBLIOGRAPHY

. . . . . . . . . . . . .

viii

. . . . . . . . . . . . . . . . . 146

LIST OF TABLES

Table 3.1

Table 6.1

Table 8.1

Table 8.2

Table 9.1

Table 9.2

Table 9.3

The reconstructed shower properties that were used in the work presented in this thesis. The term ‘Cal’ stands for calibrated. . . . . . .

40

The error status corresponding to each time-stamp is encoded into
the 4 most significant bits of the time-stamp The meaning of each
error code is shown. . . . . . . . . . . . . . . . . . . . . . . . . . . .

69

The gamma-ray flux of Markarian 421 at 7 TeV in the high flux state,
low flux state and the average flux is summarized . . . . . . . . . .

84

Summary of the search with λ = 0.01 for TeV emission from the 2FGL
list that are identified as candidates off the galactic plane. (Note that
we used the same abbreviations for the source type as the 2FGL, agu
= active galaxy of uncertain type, bzb = BL Lac type of blazar and
bzq = FSRQ type of blazar.) The Milagro flux derived at 7 TeV is
given for candidates that passed the λ = 0.01 FDR cut and the 95%
confidence level flux upper limit is given for the rest. . . . . . . . . .

89

Summary of the search with λ = 0.01 for TeV emission from the
pulsars in the second Fermi-LAT pulsar catalog. The Milagro photon
flux derived in the energy band from 34.5 TeV to 35.5 TeV is given for
the candidates that passed the λ = 0.01 FDR cut and 95% confidence
level flux upper limits are given for the rest. . . . . . . . . . . . . .

97

Properties of the Fermi-LAT pulsars in the Milagro sky coverage.
Column 2 (G100 ) is the the phase-averaged integral energy flux in
the 0.1 to 100 GeV energy band. Column 3 is the Log10 of the
pulsars spin-down luminosity. Column 4 is the pulse period. Column
5 is the first time derivative of the pulse period. Column 6 is the
distance to pulsar. All of these values are taken from the Fermi-LAT
second pulsar catalog. . . . . . . . . . . . . . . . . . . . . . . . . . .

99

The energy flux derived in the energy band from 34.5 TeV to 35.5
TeV is summarized. The first column is the name given in the FermiLAT second pulsar catalog. The second column is the name given
in the TeVCat catalog. The third column is the flux measured by
IACTs. The fourth column is the flux measured by Milagro. (a)
Flux measurements were derived from VERITAS measurements. (b)
Flux measurements were derived from H.E.S.S. measurements. . . . 101
ix

Table 10.1

G100 is the integrated phase-averaged energy flux of the pulsar GeV
emission in the 0.1-100 GeV energy band. FP W N T eV is the integrated energy flux of the PWN TeV emission in the 34.5-35.5 TeV
energy band. Properties of a sample of GeV pulsars cataloged in the
Fermi-LAT Second Pulsar Catalog and the TeV flux of their associated PWNe. (a). VERITAS Measurement. Energy flux derived
by extrapolating the SED, reference [Aliu et al.(2013)] (b). H.E.S.S.
Measurement, reference [Aharonian et al.(2006a)] (c). H.E.S.S. Measurement, reference [Aharonian et al.(2006b)] (d). H.E.S.S. Measurement, Energy flux derived by extrapolating the SED, reference [Aharonian et al.(2006d)] (e). H.E.S.S. Measurement, reference [Aharonian et al.(2005)] (f). H.E.S.S. Measurement, Energy flux derived
by extrapolating the SED, reference [H. E. S. S. Collaboration et
al.(2007)] (g). VERITAS Measurement, Energy flux derived by extrapolating the SED, reference [Aliu et al.(2013)] . . . . . . . . . . 120

Table 11.1

Summary of the search with λ = 0.01 for TeV emission from the
pulsars in the 2FGL list that were not listed in the 0FGL. (Note that
we used the same abbreviations for the source type as the 2FGL:
PSR = Pulsar identified by pulsations and psr = Pulsar identifies by
spatial association.) The Milagro flux derived at 35 TeV is given for
the candidates that passed the λ = 0.01 FDR cut and 95% confidence
level upper limits are given for the rest. . . . . . . . . . . . . . . . 136

Table A.1

Epoch live times . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145

x

LIST OF FIGURES

Figure 1.1

Figure 1.2

Figure 2.1

Figure 2.2

Figure 3.1

Figure 3.2

Figure 3.3

Figure 3.4

Figure 3.5

The Fermi-LAT sky map derived from gamma-rays of energies between 100 MeV and 10 GeV. Colored circles indicate the TeV sources
listed in TeVCat online catalog (http://tevcat.uchicago.edu/) as of
April 16th 2013. This sky map was created using the online tool
available in the TeVCat web site. . . . . . . . . . . . . . . . . . . . .

3

A schematic diagram of the technique used in Imaging Atmospheric
(or Air) Cherenkov Telescopes (IACTs). The original image was
downloaded from http:// veritas. sao. arizona. edu /about - veritas
- mainmenu - 81/ atmospheric- cherenkov - technique- and- veritastechnologies-mainmenu-87 . . . . . . . . . . . . . . . . . . . . . . .

7

A simplistic toy model for a neutron star and its magnetic sphere.
Reference: Handbook of Pulsar Astronomy by Lorimer and Kramer.

20

Pulsar Gap regions in the pulsar magnetosphere are marked in different colors. Reference: Alice Harding’s presentation “Gamma-Ray
Pulsars Theory and Modeling” at Fermi Summer School 2011. . . .

21

An aerial photograph of the Milagro detector. The pond is at the
center of this photograph and red circles mark the locations of the
outrigger tanks. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

25

This picture was taken during the cleaning process of the Milagro site
a few years after the Milagro shutdown. The Milagro water pond is
visible on the left side of this picture and several outrigger tanks are
visible around the pond. . . . . . . . . . . . . . . . . . . . . . . . .

26

This is an inside view of the Milagro pond. Two layers of PMTs are
visible in this picture. . . . . . . . . . . . . . . . . . . . . . . . . . .

28

This schematic diagram shows the PMT placement inside the Milagro
pond. The first layer, the Air Shower Layer, was placed 1.5 m below
the water surface. The second layer, the Muon layer, was placed 6 m
below the water surface. . . . . . . . . . . . . . . . . . . . . . . . . .

29

Functional block diagram of the Milagro front-end analog electronic
board. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

30

xi

Figure 3.6

Figure 3.7
Figure 3.8

Figure 3.9

Figure 4.1

Figure 4.2

Figure 4.3

The two analog signals shown at the top of this figure represent the
shape of a typical PMT output signal. The one on the left has a
high amplitude and the one on the right has a lower amplitude. The
two digital signals shown at the bottom are the outputs of the two
discriminators. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

31

Functional block diagram of the Milagro front-end digital electronic
board. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

33

A conceptual diagram of the shower core location reconstruction. As
an EAS lands on the Milagro detector, PMTs at the core of the
shower detect larger numbers of photo electrons (PEs) compared to
the PMTs at the edges of the shower. The histogram at the bottom
shows the conceptual amplitude distribution of PEs. . . . . . . . . .

37

A conceptual diagram of the shower incidence angle reconstruction in
a 2 dimensional space. As an EAS lands on the Milagro detector with
an angle θ, one end of the shower lands on the detector earlier than
the other end. The histogram at the bottom shows the conceptual
distribution of the arrival times of a shower on each PMT relative to
the shower arrival time at the first PMT. . . . . . . . . . . . . . . .

39

R (τ ) shown in the top histogram is the all sky background event
rate. ε (h, δ) shown in the second histogram from the top is the unitnormalized local-coordinate distribution of background event rate.
B (α, δ) shown in the third histogram from the top is the convolution
of the R (τ ) and ε (h, δ). B (α, δ) is the background estimation of the
two hour data set as a function of right ascension α and declination
δ. The fourth histogram from the top is the raw signal collected
within a 4 hour period, S (α, δ). The histogram at the bottom is the
background subtracted signal, S (α, δ) − B (α, δ). . . . . . . . . . . .

44

Simulated dN as a function of energy for a pure power law spectrum
dE
is shown, dN − I◦ (E/TeV)−α . I◦ ’s are selected to give the same
dE
excess for three sources with three different spectra, green α = 2.2,
red α = 2.6 and blue α = 3.0. Flux for all three sources are calculated
using spectral optimization of α = 2.6. The region where the lines
intersect indicates the energy range least dependent on the assumed
spectrum. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

48

Two significance contour plots in the α−Ec and α−I◦ space obtained
from forward folding are shown[Abdo et al.(2012)]. The spectral assumption is a power law with an exponential cutoff, as shown in
equation 4.10. The figure on the left shows the 1σ and 2σ contours
in the α − Ec space. The figure on the right shows the 1σ and 2σ
contours in the α − I◦ space. . . . . . . . . . . . . . . . . . . . . . .

50

xii

Figure 5.1

Schematic diagram of Fermi-LAT. A precision converter-tracker and
calorimeter is shown at the center of the figure. Here the calorimeter
is taken apart from the precision converter-tracker. The red doted line
indicates the track of a gamma-ray and the blue dotted lines indicate
the tracks of the cascade electron and positron. The physical size of
this detector is 1.8 m × 1.8 m × 0.72 m. It weighs 2789 kg and in
normal operation mode consumes 650 W. This figure was obtained
from the publication [Atwood et al.(2009)]. . . . . . . . . . . . . . .

55

The Aitoff projection of the significance sky map derived from the
data set used to make the Fermi Large Area Telescope Second Source
catalog. This map includes the gamma-ray energy flux measured in
the range between 100 MeV and 10 GeV. Refer to [Nolan et al.(2012)]
for the original figure and a more detailed description of data. Color
scale shows the gamma-ray significance of each sky location. . . . .

56

Figure 6.1

A drawing of a HAWC tank with 4 PMTs at the bottom. . . . . . .

61

Figure 6.2

A photograph of a fully assembled H Clock card is shown. The Clock
version of the H Clock card with 2 general purpose input ports and
8 general purpose output ports are used in the Clock System. The
Control version of the H Clock card used in the Control System has
4 general purpose input ports and 6 general purpose output ports. .

64

A block diagram of the Clock System is shown. The NAVSYNC
CW46 GPS receiver and the Clock Card are the two main components
of the Clock system. The GPS receiver is used to obtain the GPS
time and a low jitter 10 MHz signal. The Clock Card produces a 40
MHz global clock signal and two time-stamps for the scaler system
and for the TDC system. The communication between the Clock
Card and the control computer is accomplished through a A24D16
VME interface. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

65

A simplified functional block diagram of the Clock firmware is shown.
This firmware maintains an internal clock synchronized with the GPS
time and produces TDC readable time-stamps. The communication
between this module and the control computer was done by a A24D16
VME interface. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

66

Encoding of the time-stamp 12.34567 seconds with no errors is shown.
Each digit is encoded into a 4 bit binary number, and each binary
number is then encoded into a pulse. A 1µs wide pulse is used to
indicate logic 0 and, a 2µs wide pulse is used to indicate logic 1. . .

68

Figure 5.2

Figure 6.3

Figure 6.4

Figure 6.5

xiii

Figure 6.6

Figure 6.7

Figure 6.8

Figure 8.1

Figure 8.2

Figure 9.1

A block diagram of the Control System. The Control system consists
of two VME cards: the Control Card and CB Fan card. The Control
card does all the logical operations and the CB Fan card does the
appropriate level conversion to interface TDCs to the Control Card.

70

A photograph of the fully assembled CB Fan card is shown. A
CB Fan card is able to provide the level conversions and Fan-outs
required to handle 6 TDCs. . . . . . . . . . . . . . . . . . . . . . .

71

A simplified functional block diagram of the Control firmware is
shown. The Trigger, CLR, and CRST modules generate the TDC
control signals, and LNE, Pause Pulses, Busy Pulses and 10 MHz
Reference produces signals to the scaler system. Similar to the Control firmware, communication between the Control system and the
control computer was done by a A24D16 VME interface. . . . . . .

73

This map shows the 5◦ ×5◦ region around Markarian 421. The FermiLAT source position is marked by a white dot. This map is made with
− E
the spectral optimization dN/dE ∝ E −2.0 e 5T eV and the data
have been smoothed by a Gaussian point spread function. The color
of a bin shows the statistical significance (in standard deviations)
of that bin. The horizontal axis is right ascension in hours and the
vertical axis is declination in degrees. . . . . . . . . . . . . . . . . .

92

The horizontal axis is the declination in the Milagro sky coverage.
The vertical axis is the energy flux in erg cm−2 s−1 . The red triangle
data points are Milagro upper limits of BL Lac objects at 7 TeV, the
blue square data paints are the theoretically predicted flux at 7 TeV
and the black solid line is the 95% confidence level expected upper
limit of the flux for extragalactic sources corresponding to zero excess
derived at 7 TeV for each declination band of the Milagro sky maps
made with spectral assumption E −2.0 exp −E/5TeV. . . . . . . . .

93

These maps show the 5◦ × 5◦ region around the pulsars that passed
the λ = 0.01 FDR cut. The LAT source positions are marked by
white dots. These maps are made with the spectral optimization
dN/dE ∝ E −2.6 and the data have been smoothed by a Gaussian
point spread function. The color of a bin shows the statistical significance (in standard deviations) of that bin. The horizontal axis is
right ascension in hours and the vertical axis is declination in degrees. 111
xiv

Figure 9.2

PWN TeV luminosity vs. Pulsar GeV luminosity normalized with
respect to the beaming factor fΩ . Blue squares are the luminosity
derived in the energy band from 34.5 TeV to 35.5 TeV of the candidates that passed the λ = 0.01 FDR cut. Red triangles are the upper
limit of the luminosity derived in the energy band from 34.5 TeV to
35.5 TeV of the candidates that failed the λ = 0.01 FDR cut. The
linear correlation coefficient calculated only using the candidates that
passed the λ = 0.01 FDR cut is 0.73, and the best fit line for these
data points have a slope of 0.79 ± 0.05 and intercept of 3.5 ± 1.8 and
χ2 /N DF = 129.531/9. . . . . . . . . . . . . . . . . . . . . . . . . . 112

Figure 9.3

The PWN gamma-ray flux in the 34.5 to 35.5 TeV energy band is
shown (FT eV ). The blue squares are the flux of the candidates that
passed the λ = 0.01 FDR cut. The red triangles are the upper limits
of the candidates that failed the λ = 0.01 FDR cut. . . . . . . . . . 113

Figure 9.4

Fraction of the PWNe that passed the λ = 0.01 FDR cut in one˙ 2 . . . . . . . . . . . . . . . . . . . . . . . . . 113
decade bins of the E/d

Figure 9.5

Distributions of PWNe luminosity in the 34.5 to 35.5 TeV energy
band LP W N T eV vs. the spin-down luminosity E˙ of the associated
pulsars (a), and LP W N T eV vs. the characteristic age (τc ) of the
associated pulsars (b) are shown. Blue squares are the LP W N T eV
of the candidates that passed the λ = 0.01 FDR cut. Red triangles
are the upper limit of the LP W N T eV of the candidates that failed
the FDR cut. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114

Figure 9.6

TeV efficiency of PWNe, εγ , vs. E˙ (a), and εγ vs. τc (b) distributions
are shown. εγ is defined as shown in Equation 10.16. Blue squares
are the εγ of the candidates that passed the λ = 0.01 FDR cut. Red
triangles are the εγ upper limit of the candidates that failed the FDR
cut. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114

Figure 10.1

PWN TeV luminosity vs. Pulsar GeV luminosity normalized with
respect to the flux correlation factor fΩ . Blue squares are measured
by Milagro. Red stars are measured by H.E.S.S. Green circles are
measured by VERITAS. The dotted line is the best fit line. The linear
correlation coefficient is R = 0.8. The error bars are dominated by
the distance measurement uncertainties. . . . . . . . . . . . . . . . . 121

Figure 10.2

PWN TeV luminosity vs. spin-down luminosity E˙ of the associated
pulsar. Blue squares are for PWNe measured by Milagro, red stars by
H.E.S.S., and green circles by VERITAS. The distance uncertainty
contributes significantly to the error bars. The data are consistent
˙ . . . . . . . . . . . . . . 122
with no dependence of LP W N T eV on E.
xv

Figure 10.3

˙ Blue
Experimentally obtained fΩ vs. pulsar spin-down luminosity E.
squares are PWNe measured by Milagro, red stars by H.E.S.S, and
green circles by VERITAS. The best fit line for the pulsars with
E˙ > 1035.1 ergs s−1 is shown as a dotted line. . . . . . . . . . . . . 124

Figure 10.4

The distribution of the beaming factor fΩ as a function of spindown luminosity E˙ for four models derived by [Pierbattista et
al.(2012)]. This plot was reproduced by M. Pierbattista with linear fits for pulsars above E˙ = 1035 ergs s−1 . For the original figure
refer to [Pierbattista et al.(2012)]. Purple and green markers refer to
the radio-loud and radio-quiet pulsars, respectively. Black lines refer
to the best linear fits for the pulsars with E˙ > 1035 ergs s−1 . . . . 126

Figure 11.1

The fraction of Fermi-LAT sources seen by Milagro as a function of
half-decade bins of the integrated time-averaged Fermi flux ( photons
cm−2 s−1 ) in the energy range from 100 MeV to 100 GeV. . . . . 134

Figure 11.2

The fraction of Fermi-LAT sources seen by Milagro as a function of
half-decade bins of the integrated phase-averaged Fermi flux ( photons
cm−2 s−1 ) in the energy range from 100 MeV to 100 GeV. . . . . 135

Figure A.1

Comparison between Monte Carlo simulated data and real data for
cosmic-rays. In the top figure, the red data points show the distribution of the energy estimator of the Monte Carlo simulated data.
In the top figure, the blue data points show the distribution of the
energy estimator of the real data. In the bottom figure, the blue data
points show the ratio data/monte Carlo simulations. . . . . . . . . . 143

Figure A.2

Correlation between the energy estimators and true energy. Histogram to the left is the correlation between EBNN and log(Energy).
Histogram to the right is the correlation between FRASOR, the old
energy estimator, and log(Energy). Both of these plots were made
using a simulated data sample. . . . . . . . . . . . . . . . . . . . . . 143

Figure A.3

Unit-normalized distribution of true energies for events in EBNN bins,
0.2-0.4, 0.4-0.6, 0.6-0.8, 0.8-1.0, 1.0-1.2, 1.2-1.4, 1.4-1.6, 1.6-1.8 and
1.8-2.0. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144

Figure A.4

Unit-normalized distribution of true energies for events in FRASOR
bins, 0.2-0.4, 0.4-0.6, 0.6-0.8, 0.8-1.0, 1.0-1.2, 1.2-1.4, 1.4-1.6, 1.6-1.8
and 1.8-2.0. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144

xvi

Chapter 1
Gamma-ray Astronomy
Astronomy is one of the oldest natural sciences dating back to the prehistoric period. Historical records show that almost every civilization such as the Babylonians, Greeks, Chinese,
Indians, Iranians and Mayans observed the night sky and recorded them in their own way.
In the early stages, astronomers observed the night sky from the visible light emission or
reflection from cosmic objects. With the development of new technology astronomers were
able to observe the universe using electromagnetic radiation all the way from radio waves to
TeV gamma-rays.

1.1

Introduction to Gamma-ray Astronomy

The first cosmic gamma-rays were observed by Vela satellites, which were part of a US
military program. A study performed in 1973 was able to find sixteen gamma-ray bursts
in the energy rage 0.2 to 1.5 MeV1 [Klebesadel et al.(1973)]. The first dedicated satellite
missions designed to detect gamma-rays were the Small Astronomy Satellite (SAS-2) in
1 1 MeV = 106 eV
1

1972 and the COS-B mission in 1975. Both missions were able to observe the Vela and
the Crab pulsar wind nebulae. In 1981 the COS-B satellite produced a source catalog with
25 detections above 100 MeV [Swanenburg et al.(1981)]. Since then gamma-ray emission
from many different types of cosmic objects have been detected, such as pulsars, pulsar wind
nebulae (PWNe) and Active Galactic Nuclei. The latest Fermi Large Area Telescope (FermiLAT) catalog known as the 2FGL [Nolan et al.(2012)] contains 1873 sources detected in the
100 MeV to 100 GeV energy range. Most of these sources do not have a spatial association
with the sources listed in the earlier catalogs such as 1st AGILE catalog, the 3rd EGRET
Gamma-ray catalog or the EGRET Revised catalog.
The first TeV gamma-ray emission from a cosmic object, the Carb Pulsar Wind Nebula,
was detected in 1989 by a ground-based instrument called the Whipple Imaging Atmospheric
Cherenkov Telescope [Weekes et al.(1989)]. Since then TeV gamma-ray emission from many
different types of cosmic objects have been identified. As of April 16th 2013 the online
TeV source catalog called TeVCat(http://tevcat.uchicago.edu/) had 153 sources. Figure 1.1
shows the TeVCat sources marked on the Fermi-LAT sky map derived from γ-rays in the
energies between 100 MeV and 10 GeV.
At present, the GeV catalogs such as the Fermi-LAT catalogs and the TeV catalogs such
as TeVCat can be considered as two categories by the energy range in which the sources
were detected. However, with the current advancement of instrumentation the energy gap
between the GeV catalogs and the TeV catalogs is filling with new measurements. When
this gap is filled, catalogs in these two energy bands will become a single catalog with
measurements all the way from 10s GeVs to 10s of TeVs. The expansion of the energy
range is necessary to understand the characteristics of the sources, and the measurements of
a broad energy range can be turned into a coherent picture of the underlying acceleration
2

Figure 1.1: The Fermi-LAT sky map derived from gamma-rays of energies between 100
MeV and 10 GeV. Colored circles indicate the TeV sources listed in TeVCat online catalog
(http://tevcat.uchicago.edu/) as of April 16th 2013. This sky map was created using the
online tool available in the TeVCat web site.

mechanisms. Even many decades after the first discovery of GeV and TeV sources, the exact
acceleration mechanisms that produce gamma-rays observed from individual sources are
poorly understood. Most of the generally accepted mechanisms require the acceleration of
charged particles to relativistic speeds. These relativistic charged particles produce gammarays through leptonic or hadronic interactions, such as curvature radiation, synchrotron
radiation, inverse Compton scattering and Bremsstrahlung.

1.2

Modern Gamma-ray Detectors

Today’s gamma-ray detectors can be divided into two main categories: space-based and
ground-based detectors. Fermi-LAT is an example of a space-based GeV gamma-ray telescope. VERITAS and Milagro are examples of ground-based TeV gamma-ray observatories.
3

1.2.1

Space-based Detectors

Space-based detectors are carried by satellites that orbit the earth or sometimes are carried
by high altitude balloons. They are capable of detecting gamma-rays directly using the
interactions of gamma-rays with the detector material when gamma-rays pass through the
detector. The main advantages of this technique are excellent background rejection capability, very good energy resolution and long exposure. A larger duty cycle is also an advantage
of space-based detectors over the ground-based observatories that use the Air Cherenkov
technique. However, water Cherenkov gamma-ray observatories such as Milagro have comparable duty cycles. The main disadvantage of space-based detectors is the limited effective
area due to the limited physical size of the detector. The gamma-ray flux from sources
typically decreases as a power law with increasing energy. Therefore the probability of a
very high energy photon going through a detector of size ∼ 1 m2 becomes very small. This
limits the ability to use space-based detectors as gamma-ray observatories above a few 100s
of GeVs. At present the most sensitive space-based GeV gamma-ray telescope is the Fermi
gamma-ray observatory, which is discussed in Chapter 5.

1.2.2

Ground-based Detectors

Ground-based TeV observatories use the by-products of gamma-ray-earth-atmosphere interactions to determine the properties of the primary gamma-ray. When a gamma-ray enters
the earth’s atmosphere, it produces a relativistic electron-positron pair by pair production.
As these particles travel through the atmosphere, they lose energy by producing secondary
photons through bremsstrahlung; these photons also produce new electron-positron pairs.
This electromagnetic cascade continues until the photon energy drops below the threshold
4

for pair production. The result of this process is a shower of relativistic charged particles,
which is known as an Extensive Air Shower (EAS). When an EAS front propagate through
the atmosphere its lateral extent gets broader and the total energy of the shower front dissipates through the atmosphere. The altitude at which the number of particles in the shower
reaches its maximum is called the shower maximum altitude. The shower maximum altitude
depends on the energy of the initial photon and the altitude of its first pair production. Below
the shower maximum altitude the number of particles that can be detected becomes lower.
Therefore the altitude of a ground-based detector is important for detecting shower particles.

Typically TeV gamma-rays can create EAS shower fronts with a diameter of lateral extent
of ∼ 100 m. It is cheaper to built ground-based detectors to cover a large area, compared
to space-based detectors. Therefore, it is possible to build ground-based detectors with an
effective area of ∼ 105 m2 in contrast to the ∼1m2 effective area of space-based detectors.
The large effective area of a ground-based detectors increases the probability of detecting
gamma-rays with low flux at high energies.

Currently existing ground-based TeV gamma-ray observatories are divided into two categories based on the technique used : Imaging Atmospheric (or Air) Cherenkov Telescopes
and Extensive Air Shower Arrays.

5

1.3

Imaging Atmospheric (or Air) Cherenkov Telescopes
(IACTs)

Relativistic electrons and positrons in an EAS produce Cherenkov light along the air shower
direction. The maximum Cherenkov light emission occurs when the number of particles in
the cascade process is the largest. This maximum occurs at an altitude of ∼ 10 km for a
shower generated from a primary photon with energy in the 100 GeV - 1 TeV range, and
this creates a light pool on the ground. After absorption and scattering in the atmosphere,
this Cherenkov light arrives on earth as brief flashes of a few nanoseconds duration. These
Cherenkov light pulses can be collected to a focal point by placing a reflector telescope inside
the light pool. The light collected on the focal point can be detected using a fast photon
detector such as Photo Multiplier Tubes(PMTs). Figure 1.2 shows a schematic diagram of
this concept.
Not only high energy gamma-rays but also hadronic cosmic rays can create EASs, and the
Cherenkov light generated by hadronic cosmic rays are the largest background for IACTs.
Cherenkov images that result from gamma-rays appear as narrow elongated ellipses in the
focal plane, while hadronic cosmic ray initiated showers produce wider and less regular patterns. These unique features of the Cherenkov images are used to discriminate gamma-ray
photon events from hadronic cosmic ray events, and to reconstruct the photon energy and the
arrival direction of the initial gamma-ray photon. Examples of modern IACTs are High Energy Stereoscopic System (H.E.S.S.) [Hinton(2004)] and VERITAS [Gibbs et al.(2003)], and
an example of a future IACT is the Cherenkov Telescope Array (CTA) [Pareschi et al.(2013)].
The advantages of the IACT technique are high angular resolution < 0.2◦ , energy resolution (∼ 15%) and point source sensitivity. The sensitivity of modern IACTs is sufficient to
6

Figure 1.2: A schematic diagram of the technique used in Imaging Atmospheric (or Air)
Cherenkov Telescopes (IACTs). The original image was downloaded from http:// veritas.
sao. arizona. edu /about - veritas - mainmenu - 81/ atmospheric- cherenkov - techniqueand- veritas- technologies-mainmenu-87
detect the Crab pulsar wind nebula within order of hours. This sensitivity allowed IACTs to
successfully perform multi-wavelength campaigns. On the other hand the main disadvantage
of this technique is the limited duty cycle (∼ 10%) and small field-of-view. The limited duty
cycle in combination with a small field-of-view makes IACTs harder to use as survey instruments or for observing diffuse sources with large angular extents. Ground-based Extensive
Air Shower (EAS) Arrays remove these limitations.

1.4

Extensive Air Shower (EAS) Arrays

Examples of different types of modern EAS arrays are the Auger observatory, the Tibet ASgamma Experiment and the Milagro observatory. The Auger observatory consists mainly of
1,600 water tanks, separated by 1.5 km. These tanks are dark inside except when relativistic
particles from an EAS pass through the water in the tank and produce Cherenkov light.
7

These Cherenkov photons are detected by PMTs placed inside the tanks. Properties of
the EAS can be measured by studying the PMT signals. Apart from these water tanks
Auger is also equipped with fluorescence detectors to detect the atmospheric Cherenkov
light. However, the primary goal of Auger is not to detect EAS initiated by gamma-ray
photons, but to detect EAS generated by hadronic cosmic ray particles.
The Tibet AS-gamma experiment consists of scintillation counters made of plastic scintillator plates. These scintillator plates are placed on a lattice. When relativistic particles
from an EAS pass through these plates they produce photons that can be detected by the
attached PMTs. The properties of the PMT signals are used to reconstruct the properties
of the EAS, which can be related to the properties of the gamma-ray that initiated the EAS.
The major disadvantages of this technique are the high cost of scintillation plates, the poor
sensitivity to low energies, and no background rejection. The Milagro gamma-ray detector
overcame these disadvantages by using a water filled pool as the detector medium. Details
of the Milagro detector are discussed in Chapter 3.

8

Chapter 2
Gamma-ray sources
TeV gamma-ray emission has been identified from many different types of gamma-ray sources,
including Active Galactic Nuclei (AGN), Pulsar Wind Nebulae (PWNe), and Super Nova
Remnants (SNRs). Many of the detected sources are AGNs and PWNe. As of June 24th ,
2013 about 40% of the reported TeV sources in TeVCat were AGNs and about 25% were
PWNe. In this chapter I discuss the basics of non-thermal gamma-ray production processes
and briefly discuss properties of pulsars and PWNe.

2.1

Non-thermal Gamma-ray Production Processes

Non-thermal gamma-rays are produced by several processes. Sometimes those processes
are divided into two categories based on the types of high energy particles which produced
gamma-rays: leptonic mechanisms and hadronic mechanisms. In leptonic mechanisms high
energy leptons produce the gamma-rays, and in hadronic mechanisms high energy hadrons
produce the gamma-rays. Sometimes these processes are divided based on their interactions: particle-field interaction mechanisms and particle-matter interaction mechanisms. In
9

the particle-field interaction mechanisms, the gamma-rays are produced by the interactions
between particles and fields, such as synchrotron radiation of electrons in the presence of
magnetic fields. In the particle-matter interaction mechanisms, the gamma-rays are produced by the interactions between particles and the matter in the ambient field, such as
bremsstrahlung. In the very high energy gamma-ray band the most dominant radiation
mechanisms are synchrotron radiation, curvature radiation, inverse compton radiation, relativistic bremsstrahlung, and π 0 decays.

2.1.1

Synchrotron Radiation

Synchrotron radiation is created when electrons accelerate radially1 . An electron can be accelerated radially and move in a helical path when it moves with a non zero pitch angle2 and
a non-zero velocity component along the magnetic field. In such motions, accelerating electrons emit photons within an angle of θ ∼ me c2 /E [Gustavo E.R. and K.S. Cheng (2004)]
of its velocity, where me is the electron mass, c is the speed of light and E is the energy of the
electron. The total energy rate loss by synchrotron radiation generated by a helically moving
electron in a magnetic field, B, can be written as [Gustavo E.R. and K.S. Cheng (2004)],

dEe
dt

2
=− c
3

2
e2
2 γ 2,
B⊥
2
me c

(2.1)

where γ is the Lorentz factor of the electron, B⊥ = B sin θ, and θ is the pitch angle. By
averaging this over an isotropic pitch angle distribution, the average energy rate for all pitch
1 Other charged particles also produce radiation when they accelerate radially, but the
mass dependence of equation 2.1 shows that radiation by e+ and e− are much stronger than
more massive particles.
2 Pitch angle is the angle between the velocity vector of the particle and the magnetic
field.
10

angles can be written as:
dEe
dt

= 0.66 × 103 B 2 γ 2 eV/s,

(2.2)

where B is measured in Gauss. If the electron distribution is homogeneous and follows an
isotropic power-law distribution with power-law index p,

−p
Ne (Ee ) dEe = Ke Ee dEe ,

(2.3)

in a random magnetic field, then the resulting spectrum of the synchrotron radiation can be
written as:
e3
I Eγ = a(p)
me c2

3e
4πm3e c5

(p−1)/2

−(p−1)/2
B (p+1)/2 Ke LEγ
,

(2.4)

where L is the characteristic size of the emitting region, a(p) is 0.147, 0.103, 0.0852, and
0.0742 for p =1.5, 2, 2.5, and 3, respectively. The spectrum of the synchrotron radiation is
a power-law with the power-law index α = (p − 1)/2.

2.1.2

Curvature Radiation

Curvature radiation is created when charged particles move along curved field lines. In the
presence of strong magnetic fields, charged particles tend to move along the field lines. If
the curvature radius (Rc ) of the magnetic field lines is small, particles moving along the
field lines will radiate photons. The energy loss rate by curvature radiation from an electron
moving on a magnetic field line with a curvature radius Rc and zero pitch angle, can be
11

written as [Gustavo E.R. and K.S. Cheng (2004)],

dE
dt

=−

2 ce2 γ 4
,
3 Rc2

(2.5)

where γ is the Lorentz factor of the electron. When the pitch angle is zero synchrotron
radiation losses become zero. However, electrons can move in curved magnetic fields with
non zero pitch angles. The critical pitch angle is defined by equating the synchrotron and
the curvature radiation losses, which yields:

sin θcr =

γme c2
,
eBRc

(2.6)

where θcr is the critical pitch angle.
Curvature radiation becomes the dominant radiation process when the pitch angle is
much smaller than the critical pitch angle. However, the synchrotron radiation losses become
significant when the pitch angle is of order of the critical pitch angle. Therefore, in this case
the total radiation losses can not be explained by either synchrotron radiation or curvature
radiation. The radiation losses when the pitch angle is on order of the critical pitch angle is
explained by the synchro-curvature radiation, which is a mixture of the two processes (eg.
[Cheng & Zhang(1996), Zhang & Yuan(1998), Harko & Cheng(2002)]).

2.1.3

Inverse Compton Radiation

Relativistic electrons can be scattered off low energy photons by transferring the energy
of the relativistic electrons to the low energy photons. This process goes in the inverse
direction of Compton scattering (where high energy photons scatter off low energy electrons
12

by transferring energy to the electrons). Therefore, this process is called Inverse Compton
scattering (IC scattering). When Ee Eph

2
me c2
the electron can transfer most of its

energy to the photon, Eγ ∼ Ee ,where Eph is the average energy of the target photons and
Eγ is the average energy of the photons after the IC scattering..
Target photons for IC scattering can come from multiple sources, such as microwave
photons from the cosmic microwave background, ambient photons from thermal emissions, or
high energy photons from synchrotron radiation. One can consider a special case where target
electrons for IC scattering come from synchrotron radiation of the same electron population
that produces the IC scattering. This special case is called Synchrotron Self-Compton (SSC)
scattering. In this scenario the radiation from the synchrotron radiation and from the IC
scattering are highly correlated. The spectrum of the inverse compton component is expected
to be similar to the synchrotron component [Costamante & Ghisellini(2002)].

2.2

Pulsars

Pulsars are highly magnetized, rotating neutron stars that emit two beams of electromagnetic
photons in opposite directions. A pulsar can only be seen when a beam is pointing towards
the observer. Therefore, the observer sees the pulsar as a pulsating star. The first pulsar was
accidentally discovered by Jocelyn Bell Burnell and Antony Hewish on November 28, 1967
[Hewish et al.(1968)]. Since then more than 2000 radio pulsars3 has been discovered, and
the recently published second Fermi-LAT pulsar catalog observed 117 pulsars in the GeV
energy band [Abdo et al.(2013)].
3 According to the ATNF online catalog http://www.atnf.csiro.au/people/pulsar/psrcat/
13

2.2.1

Properties of Neutron Stars

The composition of neutron stars was predicted in 1939 by Oppenheimer and Volkov, 28
years before the discovery of the first pulsar [Oppenheimer & Volkoff(1939)]. A neutron star
is created as a result of the core collapse of a massive star during a supernova. The resulting
neutron star retains most of its angular momentum. However, the radius of the star is
sharply reduced. Therefore, the angular velocity is sharply increased and the density of the
neutron star becomes high. The equation of state for neutron stars, which consist of highly
compressed matter, remains uncertain, and the predicted maximum neutron star mass is
∼ 2M

[Lattimer & Prakash(2001)]. Measurements of neutron star masses have been found

to typically be ∼ 1.4M

[Lorimer & Kramer(2012)].

The possible range of radii for neutron stars can be derived using theoretical arguments.
The theoretically possible minimum radius for a neutron star is derived based on arguments that the speed of sound should be less than the speed of light in a neutron star,
and that the equation of state should have a smooth transition from low to high densities
[Lattimer et al.(1990), Glendenning(1992)],

Rmin

1.5Rs =

3GM
= 6.2km
c2

M
1.4M

,

(2.7)

where Rs is the Schwarzschild radius. The theoretically possible maximum radius is derived
based on the argument that the neutron star has to be stable against break-up due to
centrifugal forces leading to,

Rmax

GM P 2
4π

1/3
= 16.8km
14

M
1.4M

1/3

P 2/3
,
ms

(2.8)

where P is the period of rotation.
Neutron stars are nearly spherical objects. Therefore, the momentum of inertia can be
calculated using I = kM R2 . For a sphere of uniform density k = 0.40. However, neutron
stars are not perfect spheres with uniform density. Therefore, the exact value of k depends on
the density profile, which depends on the equation of state. The typical range of k predicted
by most models is 0.3 < k < 0.45 for a mass/radius range of 0.1M km−1 < M/R <
0.2M km−1 [Lorimer & Kramer(2012)]. Most practical calculations adopt M = 1.4M ,
R = 10 km and k = 0.4 to calculate I, so I = 1038 kg m2 .

2.2.2

Pulsar Spin Evolution

It has been observed that pulsar periods increase with time, P˙ = dP > 0. Pulsar periods
dt
are increasing because of the loss of rotational kinetic energy. Therefore, P˙ can be related
to the rate of loss of rotational kinetic energy,
d IΩ2 /2
dE
rot
=−
= −IΩΩ˙ = 4π 2 I P˙ P −3 ,
E˙ ≡ −
dt
dt

(2.9)

where Ω = 2π/P . The E˙ of a pulsar is called the spin-down luminosity, which represents
the total rate of energy output of the neutron star, by neglecting changes of the energy due
to changes in the magnetic field. Equation 2.9 can be rewritten using the typical value for
the momentum of inertia, I = 1038 kg m2 ,

E˙

3.95 × 1031 erg s−1
15

P˙
10−15

P −3
.
s

(2.10)

A small fraction of the E˙ is converted into electromagnetic radiation, and to a wind of
electrons and positrons flowing out from the pulsar. The bulk of the E˙ converts into magnetic
dipole radiation.

2.2.3

Spin-down Luminosity Evolution

A pulsar is a rotating magnetic dipole4 , which radiates an electromagnetic wave at its rotational frequency,
2
E˙ dipole =
|m|2 Ω4 sin2 α,
3
3c

(2.11)

where m is the magnetic moment and α is the angle between the magnetic moment and
the spin axis. Since the bulk of the E˙ converts into the magnetic dipole radiation, E˙ is
approximately equal to E˙ dipole . This equality yields:

Ω˙ = −

2|m|2 sin2 α
3Ic3

Ω3 .

(2.12)

This equation can be more generally expressed as follows:

Ω˙ = −KΩn ,

(2.13)

where K is a constant and n is the braking index. The braking index can be determined by
taking the second derivative of Ω and eliminating K,

n=

¨
ΩΩ
.
Ω˙ 2

4 Assuming the magnetic field is constant
16

(2.14)

The braking index of several pulsars were measured and found to have a typical range of
1.4–2.9. These measurements indicate that the assumption of n = 3 is usually not correct.
The theoretical expectation of the rotational energy loss rate of a neutron star is a
function of the braking index [Mattana et al.(2009)],

˙
E(t)
=

E˙ ◦
,
1 + t/tdec n

(2.15)

where E˙ ◦ is the E˙ at the time of the pulsar birth, tdec is a characteristic decay time, t is
the time elapsed since the pulsar’s birth, and n is the braking index.

2.2.4

Characteristic Age of a Pulsar

The generalized model for angular velocity evolution that is shown in Equation 2.13 can be
expressed in terms of the pulsar period (P) and braking index (n),

P˙ = kP 2−n

(2.16)

This is a first-order differential equation, which can be integrated assuming k is a constant
and n = 1,

T =

P
1−
(n − 1)P˙

P◦ n−1
,
P

(2.17)

where T is the pulsar age and P◦ is the pulsar period at birth. This equation can be
rewritten under the assumptions that the current pulsar spin period is much larger than
that at the birth of the pulsar (P◦ << P ) and the dominant spin-down is due to magnetic
17

dipole radiation (n = 3),
τc =

P
,
2P˙

(2.18)

where τc is the characteristic age. Due to these assumptions τc does not necessarily estimate
the correct age of a pulsar. For example the characteristic age of the Crab pulsar is 1240
years, compared to its known age of 950 years. Another example is the pulsar J0205+6449m
which was born in the historical supernova of AD 1181, but has the characteristic age of 5370
years. Therefore, characteristic age cannot be considered a reliable parameter to estimate
the real age of a pulsar, and it should be interpreted with some care.

2.2.5

Pulsar Magnetosphere

The structure of a pulsar magnetosphere and its acceleration processes can not be solved
analytically in a closed form [Lorimer & Kramer(2012)]. Therefore, numerical simulations
are used to study the structure and the acceleration processes of pulsar magnetospheres.
A simplistic toy model for a neutron star and its magnetosphere is shown in Figure 2.1,
where the magnetosphere is the region surrounding the neutron star in which its magnetic
field is the predominant effective magnetic field. In a simplistic scenario a neutron star
can be considered as a superconducting object. In the case of a pulsar neutron star it is a
spinning magnetic dipole. This conducting and spinning magnetic dipole acts as a uni-polar
generator, which creates an electric quadrupole field around the neutron star. The electric
field component parallel to the magnetic field strips the electrons from the pulsar surface
and accelerates them along the magnetic field lines. These accelerated electrons produce
curvature radiation. Due to the geometry of the magnetic field, the photons from curvature
radiation create two photon beams along the magnetic field axis. The photons from curvature
18

radiation produce secondary electron-positron pairs (by interacting with the strong magnetic
field), which are also accelerated along the magnetic field lines by the electric field.
Co-rotating magnetic field lines are closed loops near the rotating neutron star. However, the co-rotating magnetic field lines become open as they cross the light cylinder. The
light cylinder is a co-rotating imaginary cylinder, centered on the pulsar, aligned with the
rotational axis, and with a surface speed equal to the speed of light. Electrons and positrons
that accelerate along the open field lines create an electron-positron wind flowing out from
the pulsar. This wind is called the pulsar wind.
Different acceleration mechanisms require different conditions to accelerate electrons and
positrons to very high energies. Some models expect large plasma densities, which are not
possible to generate on the superconducting neutron star. These models predict regions in
the magnetosphere called gap regions, such as: Polar cap, Outer gap and Slot gap. The
Polar cap is the region where open field lines meet the neutron star surface, the Outer gap
is the region where Ω · B = 0, and the Slot Gap is the edge of the open field region from the
neutron star surface to near the light cylinder. These gap regions are shown in the schematic
diagram in Figure: 2.2.

2.3

Pulsar Wind Nebulae

A Pulsar Wind Nebula (PWN) is a nebula that is powered by the ultra-relativistic pulsar
wind of a pulsar. The pulsar wind flows out ward from the pulsar along its magnetic field
lines. This wind is composed from ultra relativistic electrons and positrons. The total
injection rate of electrons and positrons to the PWN by its associated pulsar can be written
as [Mattana et al.(2009)]:
19

Figure 2.1: A simplistic toy model for a neutron star and its magnetic sphere. Reference:
Handbook of Pulsar Astronomy by Lorimer and Kramer.

20

Figure 2.2: Pulsar Gap regions in the pulsar magnetosphere are marked in different colors.
Reference: Alice Harding’s presentation “Gamma-Ray Pulsars Theory and Modeling” at
Fermi Summer School 2011.

21

N˙ =

E˙
,
Γw me c2 (1 + σ)

(2.19)

where Γw is the bulk Lorentz factor of the wind and σ sets the fraction of the spin-down
˙ converted into the kinetic energy of the wind. Note that only a very small
luminosity (E)
fraction of the E˙ is converted into the kinetic energy of the wind, therefore σ

1. To

simplify the math, in this discussion Γw was assumed to be a constant. As I described in
Section 2.2.3, E˙ of a pulsar decreases with time. Therefore, the value of N˙ is proportional
˙
to the value of the spin-down luminosity, E(t),
at a given time.

When this wind meets the nebula it creates a standing shock, which makes non-zero pitch
angles (pitch angle is the angle between the magnetic field and the velocity of the particle).
These particles with non-zero pitch angles lose their energy by synchrotron radiation. The
electrons also create very high energy gamma-ray photons by Inverse Compton scattering.
The cooling time of electrons that cool by Inverse Compton scattering is typically larger than
the age of young pulsars [Mattana et al.(2009)]. Therefore, the electrons that accumulated
from the birth of a PWN could still produce TeV photons by IC scattering. The total number
of electrons that accumulated in the PWN can be written as:
τ
nu ∝

˙
E(t)dt,

(2.20)

0

˙
where nu is the number of accumulated electrons and τ is the pulsar’s age. E(t)
can be
substituted from Equation 2.15, giving:
τ
nu ∝

0

1 + t/tdec
22

E˙ ◦
dt.
(n+1)/(n−1)

(2.21)

In the case where the dominant spin-down losses are due to magnetic dipole radiation (n = 3),

nu ∝ E˙ ◦ tdec

Typically τ

τ
τ + tdec

.

(2.22)

tdec [Mattana et al.(2009)], which yields:

nu (E, t) ∝ E˙ ◦ tdec .

(2.23)

Since the luminosity is roughly proportional to the population of the radiating particles,
luminosity is also proportional to E˙ ◦ tdec . Therefore, one can expect that for PWNe LT eV ∝
E˙ 0 (PWNe TeV luminosity is proportional to the spin-down luminosity to the power of zero).

23

Chapter 3
The Milagro Observatory
The Milagro gamma-ray observatory was a ground-based water Cherenkov detector located
near Los Alamos, New Mexico, USA at latitude 35.9◦ north, longitude 106.7◦ west and
altitude 2630 m. Milagro was operated from 2001 to 2008 and was sensitive to gamma-ray
initiated EASs from a few hundred GeV to ∼100 TeV. Its high duty cycle (> 90%) and wide
field of view (∼ 2sr) made it an excellent instrument to survey the northern hemisphere
−7◦ < DEC < 80◦ . The Milagro archived data set is an ideal place to study phenomena
of gamma-ray sources as well as hadronic cosmic rays. This chapter describes the working principles, electronics, event reconstruction and background subtraction algorithm of
Milagro.

3.1

Milagro Detector

The Milagro detector consisted of a light tight pond of 80 m (l) × 60 m (w) × 8 m (h) filled
with 24 million liters of ultra purified water surrounded by 175 small outrigger tanks. The
area of the Milagro detector including the outrigger tanks was about 40, 000 m2 . Figure 3.1
24

Figure 3.1: An aerial photograph of the Milagro detector. The pond is at the center of this
photograph and red circles mark the locations of the outrigger tanks.

shows an aerial photograph of Milagro and the distribution of the outrigger tanks are marked
with red circles. Figure 3.2 shows some outrigger tanks. The main pond was instrumented
with 723 Photo Multiplier Tubes (PMTs) arranged in two layers and each outrigger tank
was instrumented with a single PMT.

3.1.1

The Main Pond

A photograph of the inside of the Milagro pond is shown in Figure 3.3. This photograph
shows the two layers of the PMTs. The top layer, called the Air Shower (AS) layer is placed
1.5 m below the water surface and consists of 450 PMTs. The bottom later, called the
Muon (MU) layer is placed 6 m below the water surface and consists of 273 PMTs. The
PMTs in each of these layers were arranged on a 2.8 m × 2.8 m grid and the MU layer is
25

Figure 3.2: This picture was taken during the cleaning process of the Milagro site a few
years after the Milagro shutdown. The Milagro water pond is visible on the left side of this
picture and several outrigger tanks are visible around the pond.

26

horizontally shifted by a half grid space from the AS layer. A schematic side view of this
PMT arrangement is shown in Figure 3.4.
When an EAS front lands on the Milagro pond, the EAS’s relativistic electrons and
positrons travel through the water making Cherenkov light in the water that can be detected by PMTs. In the meantime, the secondary photons generated in an EAS also enter
the pond. Since the AS layer is placed 1.5 m below the water surface (equivalent to 4 radiation lengths) these gamma-rays can cascade into electron-positron pairs before reaching
this layer. These tertiary electron-positron pairs are also capable of producing Cherenkov
light. The secondary showering increases the density of Cherenkov photons inside the water,
which increased the sensitivity of Milagro. Since the MU layer is placed 6 m under the water
surface (16 radiation lengths of water) the energy of electrons and positrons are not enough
to reach that depth. However, muons and hadrons that mostly constitute hadronic showers
(showers that initiated from a highly energetic hadronic cosmic particle) can easily penetrate
to this depth. Therefore, the ratio between the number of AS layer PMT hits to the number
of MU layer PMT hits gives a hint about the origin of the shower. This property was used
to reject the background of hadronic generated showers. Aous Ahmad Abdo [Abdo(2007)]
was able to identify a more efficient method to reject the hadronically generated showers by
combining the AS layer, MU layer and OR tanks.

The PMTs used both in the main water pond and in the outrigger tanks were 20 cm
Hamamatsu #R5912SEL with a custom water-proof base made from polyvinyl chloride
(PVC). The base protects the electronic connections of the PMT, which carry the high
voltage to the PMT and the PMT pulses back to the front-end electronics. As can be seen
in Figure 3.3, each PMT is surrounded by a conical collar. The inside of this collar is white
27

Figure 3.3: This is an inside view of the Milagro pond. Two layers of PMTs are visible in
this picture.

and outside is black. The white interior increased the light collection area of each PMT and
the black exterior blocked the Cherenkov light generated from horizontally traveling muons.
In order to detect water Cherenkov light using PMTs, the pond had to be dark and the
water ultra-pure. Darkness inside the pond was accomplished by covering the pond using
a 1 mm thick polypropylene cover. The top surface of the cover is painted with a highly
reflective paint to reduce the inside temperature during the day time. The quality of water
inside the pond was maintained by constantly recirculating water at a rate of 200 gallons per
minute. During the re-circulation process, the water went through a set of filters : a charcoal
filter, a 10 µm filter, a 1 µm filter, a carbon filter, a 0.2 µm filter1 and finally through UV
lamps to prevent any biological growth.

1 This filter was removed after the second year.
28

Figure 3.4: This schematic diagram shows the PMT placement inside the Milagro pond.
The first layer, the Air Shower Layer, was placed 1.5 m below the water surface. The second
layer, the Muon layer, was placed 6 m below the water surface.

3.1.2

The Outrigger Array

The outrigger array consisted of 175 tanks of 2.4 m in diameter and 1 m in height. Each tank
was filled with ultra-purified water and instrumented with a PMT pointing down. These
outrigger tanks were unevenly spread across an area of 40, 000 m2 . This array allowed for
more accurate determination of the shower core location and its arrival direction.

3.2

Front-end Electronics and Time to Digital Converter

The Milagro front-end electronics were designed to digitize the analog signal outputs coming
from the PMTs. This task was accomplished using two types of custom electronic cards,
called analog and digital front-end Boards, and a commercial Time to Digital Converter
(TDC). Each analog and digital card pair was connected to 16 PMTs.
The functional block diagram of one channel of the analog cards is shown in Figure 3.5.
The input of each channel is the analog signal output coming from a PMT. Then the input
signal splits into two and goes to two transconductance amplifiers (where output current is
proportional to the amplitude of the input signal) with two different gains. The output of
29

Figure 3.5: Functional block diagram of the Milagro front-end analog electronic board.
each transconductance amplifier is then sent to an RC(Resistor-Capacitor) integrator. The
voltage across the RC integrator’s capacitor can be written as:
t
Vc ∝

0

Vin dt.

(3.1)

Each Vc is then sent to one of two discriminators with a preset threshold voltage. One discriminator has a high threshold voltage (of 80 mV) and the other one has a lower threshold
voltage (of 30 mV). Whenever the PMT pulse crosses either of these thresholds the discriminator outputs produce an edge as shown in Figure 3.6. Smaller pulses create an edge pair
only at the low threshold discriminator output and pulses with larger amplitudes create edge
pairs at both discriminator outputs. The width of each digital pulse is equal to the time
spent over the preset threshold. Therefore these signals are called Time Over Threshold
(TOT) signals.
30

Figure 3.6: The two analog signals shown at the top of this figure represent the shape of a
typical PMT output signal. The one on the left has a high amplitude and the one on the
right has a lower amplitude. The two digital signals shown at the bottom are the outputs
of the two discriminators.

31

The TOT signals are then sent to digital cards, where two discriminator output signals
are combined to make one signal. Figure 3.7 shows the functional block diagram of the
digital cards. The basic functionality shown in this figure can be summarized in several
steps.
1. Delay the high TOT signal by 20 ns.
2. Invert the high TOT signal.
3. Add the inverted signal to the low TOT signal.
The result of this process is a four edge event with the edge 1 to edge 4 distance equal to
the low TOT pulse width, and the edge 2 to edge 3 distance equal to the high TOT pulse
width. If the signal is too small to cross the high threshold, it would produce a two edge
event.
The next step is to convert these digital signals to a number that can be read by a
computer. This was done using LeCroy FASTBUS Time to Digital Convertors (TDCs), that
produce an integer proportional to the TOT values. Since these TOT values are a function of
the PMT pulse size, these digital numbers can then be calibrated to get the charge generated
by the PMTs, under the assumption that all PMT signals have a similar shape.

3.3

Trigger

A better way to have operated the Milagro data acquisition system would have been to
record the data continuously and throw out the background events after a high-level data
analysis. However, this was not possible to do at the time due to technical limitations at
that time. Instead the Milagro data acquisition system only accepted events that fulfilled a
32

Figure 3.7: Functional block diagram of the Milagro front-end digital electronic board.

33

trigger condition, then the triggered raw data was saved to Digital Linear Tapes (DLTs). At
the beginning of the Milagro operation, the trigger condition was the number of PMTs fired
during a given time period, which is called the simple multiplicity trigger (SMT). The trigger
conditions of the SMT were changed from time to time but in general the trigger condition
was ∼ 60 AS layer PMT hits within a 200 ns time window. Lowering this threshold would
have increased the sensitivity to low energy gamma-rays, which is physically motivated because of the low flux at very high energies. However, lowering the threshold would have
increased the data rate (typical trigger rate was ∼ 2kHz), which could not be handled by
the data acquisition system.
The dominant background events that create low multiplicities of PMT hits are high energy
muons that travel nearly horizontal with respect to the Milagro pond. An alternative way
to identify such events is by their time profile. When an air shower lands on the detector
the shower front sweeps along the pond faster than a muons travel through water. Therefore
muon events can be identified by using their rise time parameter, which is the time interval
within which 10% to 90% of the hits occur for a given time. The rise time parameter is
larger for muon events compared to events generated by air showers.
The calculation of the rise time parameter was not as easy as the calculation of the multiplicity. Therefore, Milagro designed a customized piece of electronics with programmable logic
to perform this task. This device was used for triggers for a while, but had to be removed
after it failed.

34

3.4

Reconstruction

The triggered raw data consists of a collection of edge times. The next step is to reconstruct
the properties of each shower, such as its core location and incidence angle, from these edge
times.

3.4.1

Shower Core Location Reconstruction

The shower core is the location where the primary particle would hit the Milagro detector if
there were no interactions with the atmosphere. The particle density of a shower is highest
at the core location and the core contains the highest energy particles. Therefore the number
of photoelectrons (PEs) recorded by the PMTs are highest at the shower core. This property
was used in Milagro to identify shower cores, and a conceptual diagram is shown in Figure 3.8.
Over Milagro’s life time several algorithms were used to identify the shower core location.
The simplest algorithm used in Milagro was the center of mass fitter, which calculates
the core location using the weighted sum of the PMTs.

xcore =

P Ei × x i

i
i

ycore =

P Ei
P Ei × y i

i
i

P Ei

(3.2)

(3.3)

In Equation 3.2 and 3.3, ( xcore , ycore ) are the Cartesian coordinates of the shower
core location,

xi , yi

are the Cartesian coordinates of the PMTs and P Ei is the number

of PEs detected by the PMT at

xi , yi . If the shower core location was likely to be on

the pond only the PMTs in the pond were used, but if it was likely to be off the pond only
the PMTs in the outrigger tanks were used. The likelihood of the core being on or off the
35

pond was calculated using the ratio of the number of outrigger hits to the number of AS
layer PMT hits.
The core fitter algorithm used later in Milagro was based on a 2-D Gaussian least-squares
fit. The underlying concept of this algorithm is the same as the previous algorithm. It fits
a 2-D Gaussian function by considering xi and yi as the two dimensions and the number of
PEs detected by each PMT as the amplitude. The peak of the Gaussian distribution is the
core location.
These simplified versions of the algorithms explained here assume that EAS fronts are flat
and have a uniform thickness. However, EAS fronts are not flat and are thin near the shower
core, compared to the edge. Therefore the angle reconstruction process makes corrections to
remove these effects.

3.4.2

Incidence Angle Reconstruction

The incidence angle of an EAS is reconstructed using the timing profile of the PMT hits.
This concept is easy to understand using the two dimensional example shown in Figure 3.9.
When an EAS front lands on the ground with a non-zero zenith angle, one end of the shower
lands on the detector earlier than the other, and then the shower sweeps through the detector.
The time that a shower takes to travel from one PMT to the other can be written as,

δt =

L sin θ
,
V

(3.4)

where δt is the relative time that the shower takes to travel from one PMT to the other,
L is the gap between PMTs, θ is the incidence angle and V is the velocity of the shower
front. Since shower fronts travel with a speed close to the speed of light, this equation can
36

Figure 3.8: A conceptual diagram of the shower core location reconstruction. As an EAS
lands on the Milagro detector, PMTs at the core of the shower detect larger numbers of
photo electrons (PEs) compared to the PMTs at the edges of the shower. The histogram at
the bottom shows the conceptual amplitude distribution of PEs.

37

be written as follows.
θ = arcsin

cδt
L

(3.5)

Hence, θ can be calculated using δt and L as measurable parameters. In order to get a good
fit, the angle fitter algorithm uses all the PMT hits in a shower plane to reconstruct the
shower angle.
The angle fitter algorithm used in Milagro was an iterative algorithm that fit a shower
plane in 3D space. The iteration starts by assigning a weight to each PMT and a cut based
on the number of PEs in each PMT. The PMTs that detected fewer PEs get a lower weight.
This weighted timing distribution is then fit to a shower plane using a 3D version of the 2D
algorithm explained in the above paragraph. At the end of an iteration, PMTs with poor
residuals are removed from the fit and the process is repeated with a looser PE cut. This
process is iterated five times and counts the number of PMTs used in the final fit (nFit).

3.5

Monte Carlo Simulations

Milagro Monte Carlo (MC) simulations were done in two stages. In the first stage EASs were
generated and propagated through the atmosphere using CORSIKA (COsmic-Ray SImulations for KAskade) v6.5021. CORSIKA simulations start with the first interaction of the
primary particle in the atmosphere. Then the shower propagates through the atmosphere
with the interactions of the secondary particles until the shower front reaches the ground.
In the Milagro MC simulations, both gamma-ray and hadron initiated showers were thrown
with a power law spectrum of E −2.0 . Later in the MC simulation chain showers were properly re-weighted to match the hadronic cosmic ray spectrum and to match the ratio of hadron
generated showers to gamma-ray generated showers.
38

Figure 3.9: A conceptual diagram of the shower incidence angle reconstruction in a 2 dimensional space. As an EAS lands on the Milagro detector with an angle θ, one end of the
shower lands on the detector earlier than the other end. The histogram at the bottom shows
the conceptual distribution of the arrival times of a shower on each PMT relative to the
shower arrival time at the first PMT.

39

The second stage of the MC simulation chain is to simulate the response of the Milagro
detector to EASs, which was done using GEANT4(GEometry ANd Tracking). GEANT4
takes the CORSIKA output at ground level and propagates the shower particles through a
model of the Milagro detector. The output of GEANT4 includes the number of PEs and
the arrival time of each PMT hit, which can be turned into a TOT signal as described in
Sction 3.2. Then this MC simulated data sample is sent through the reconstruction software
to reconstruct the shower properties. The reconstructed shower properties that were used in
the work presented in this thesis are shown in Table 3.5.
Data field name Content
nCalAS
Number of PMT hits in the Air Shower Layer.
nLiveASTubes
Number of live PMTs in the Air Shower Layer.
nCalOR
Number of PMT hits in the Outrigger array.
nLiveORTubes
Number of live PMTs in the Outrigger array.
nCalMU
Number of PMT hits in the Muon Layer.
Theta
Zenith angle of the shower.
XCore
X coordinate of the shower core.
YCore
Y coordinate of the shower core.
nFit
Number of PMTs used in the shower plane fit,
.
Refer to Section 3.4.2 for more details.
A5
This is the gamma-hadron separation parameter,
Refer to Section 4.1 for more details.
Table 3.1: The reconstructed shower properties that were used in the work presented in this
thesis. The term ‘Cal’ stands for calibrated.

40

Chapter 4
Milagro Data Analysis Framework
The Milagro reconstructed data set is a collection of gamma-ray and hadron initiated EASs.
In order to search for gamma-rays coming from celestial gamma-ray sources, this data set
needs to be processed into a useful form. This chapter briefly explains the separation of
gamma-ray and hadron initiated showers, estimation of the background, calculation of the
significance and finally how to calculate the TeV flux of a sky location.

4.1

Gamma/Hadron Separation Parameter

The Milagro reconstructed data set consists of both gamma-ray initiated EASs, which are
signal events, and hadron initiated EASs, which are background events. Distinguishing the
signal from the background with a 100% efficiency is an impossible task. Therefore, the data
has been statistically discriminated using the parameter A5, which was defined as:

A5 = 400 ·

N
NAS
+ OR · ζ (t) · Ff it
live
live
NAS
NOR
,
M axP EM U
41

(4.1)

live is the fraction of live Air Shower layer PMTs that detected Cerenkov light
where NAS /NAS
live is the fraction of live Outrigger PMTs that detected Cerenkov
in the event, NOR /NOR
light in the event, Ff it is the fraction of live PMTs used in the shower fitting algorithm explained in Section 3.4.2, ζ (t) is a run-dependent correction which corrects for the systematic
variations of each run such as changes in calibration, the denominator M axP EM U is the
maximum number of photo-electrons recorded by a PMT in the muon layer and the constant
400 is an arbitrary scaling factor that gives A5 typical values between 0 and 10.
The numerator of A5 is a function of the shower size, which becomes larger for showers
with larger size. The denominator of A5 is expected to be larger for hadron initiated EASs
relative to the gamma-ray initiated EASs, because of the high muon content in hadron
initiated EASs. Therefore, A5 is expected to typically be larger for gamma-ray initiated
EASs, and rejecting events with lower A5 values decreases the background contamination in
the Milagro data set. Typical data analyses use a cut of A5 > 1 for background rejection.
Even after this cut was applied, the Milagro data set was ∼ 99% background events.

4.2

Background Estimation

The search for gamma-ray events in the Milagro data set is obscured by the ∼ 99% background contamination. Therefore, the correct estimation of the background is a very important step. Underestimation of the background could cause background fluctuations to
appear as real signals, and overestimation could cause real signals to be hiddan in the data.
The work presented in this thesis used the direct integration method [Abdo et al.(2012)] to
estimate the background.
For a short time interval (typically 2 hours) the local angular distribution of the back42

ground events can be written as a function F (h, δ) depending on the local hour angle (h)
and declination (δ). This function can be turned into an efficiency distribution for a given
declination band by normalizing the function into a unit area, ε (h, δ) = F (h, δ)/ F (h, δ)dh,
which gives the efficiency of a hadronic initiated EAS that comes from the local coordinates
of (h, δ). The convolution of the efficiency distribution with the all-sky background rate
R (τ ) gives the expected background coming from the local coordinates of (h, δ). Therefore,
the background coming from a sky location with coordinates (α, δ) can be written as:

B (α, δ) =

ε (h, δ) · R (t) · Γ (h, R.A., t) dt dΩ,

(4.2)

where Γ (h, R.A., t) is 1 if the hour angle, right ascension, and sidereal time are such that
the event falls within the right ascension, declination bin of interest, and zero otherwise. In
Milagro data analysis, this algorithm was applied to the full Milagro data set with a two hour
integration time, because ε is a constant for short time intervals but we also have to chose
a large enough time interval to get enough statistics. An example of this analysis is shown
in Figure 4.1 using a four hour long Milagro data set. After computing B (α, δ), the excess
distribution can be obtained by subtracting B (α, δ) from signal S (α, δ). This final step is
shown in the bottom plot of the Figure 4.1. The example shown in Figure 4.1 extracted
the excess distribution of a narrow-declination band as a function of right ascension using
a four hour long data sample, with integration time also equal to four hours. The result
of this process is an excess map of the celestial sphere in the equatorial coordinate system.
Real Milagro maps are two dimensional maps of right ascension and declination. The excess
distribution for each declination band was extracted using this algorithm.

43

R (τ)

2500
2000
1500
1000
500
0

∈(h,δ )

0

2

4

6

8

10

12

14

16

18

Sidereal Time (hrs)
20
22
24

0.08
0.06
0.04
0.02

B(α ,δ )

0

-10

-5

0

Hour Angle (hrs)
10

5

2000
1500
1000
500
0

S(α ,δ )

0

Right Ascension (hrs)
20
22

2

4

6

8

10

12

14

16

18

24

2

4

6

8

10

12

14

16

18

Right Ascension (hrs)
20
22
24

2

4

6

8

10

12

14

16

18

Right Ascension (hrs)
20
22
24

2000
1500
1000
500
0

S(α ,δ ) - B(α ,δ )

0

100
0
-100

0

Figure 4.1: R (τ ) shown in the top histogram is the all sky background event rate. ε (h, δ)
shown in the second histogram from the top is the unit-normalized local-coordinate distribution of background event rate. B (α, δ) shown in the third histogram from the top is the
convolution of the R (τ ) and ε (h, δ). B (α, δ) is the background estimation of the two hour
data set as a function of right ascension α and declination δ. The fourth histogram from the
top is the raw signal collected within a 4 hour period, S (α, δ). The histogram at the bottom
is the background subtracted signal, S (α, δ) − B (α, δ).

44

4.3

Significance Maps

The next step after obtaining the excess map is to calculate the statistical significance of the
excess. The significance calculation has been done using an extended version of the statistical
method explained in Li & Ma(1983) with weighted events. The weights are assigned to
optimize the statistical significance of a gamma-ray source in the sky. The Milagro analysis
framework applied an event weight to each event based on two factors: the A5 value and the
expected signal-to-background ratio for an assumed spectrum as a function of the parameter
call FRASOR (F ):
F =

N
NAS
+ OR .
live N live
NAS
OR

(4.3)

FRASOR is an energy dependent parameter, therefore the expected number of events in
each F bin1 depends on the energy spectrum. The expected number of gamma-ray events
for a given spectrum can be estimated using the Monte Carlo simulation. In this algorithm,
an event that falls in an F bin with a higher gamma-ray event expectancy gets a higher
weight compared to an event that falls in an F bin with a lower gamma-ray event expectancy.
Similarly events with a larger A5 value get a higher weight compared to events with a smaller
A5. Therefore, an event with a larger A5 that falls in an F bin with a higher gamma-ray
event expectancy gets the highest event weight.
In addition to this weight, events are given another weight to account for the angular
resolution of Milagro. For a given event, the incident angle can only be reconstructed with a
finite angular resolution. Therefore, each recorded event is spread around the reconstructed
angle assuming the angular resolution function is a 2-D Gaussian distribution. The width
of the Gaussian distribution depends on the F bin, which ranges from 1.2◦ for small F bins
1 Typically F is binned into 10 equal size bins.
45

to 0.35◦ for large F bins.
Since the weights depend on the assumed spectrum, the sky map obtained with this
algorithm has a bias towards the assumed spectrum. Therefore, the Milagro analysis framework produced multiple Milagro sky maps with different spectral assumptions. The analysis
presented in this thesis used two Milagro sky maps: a sky map made with a pure power
law spectrum of E −2.5 was used to search for galactic sources, and a sky map made with
a power law with an exponential cutoff, E −2.0 e−E/5T eV , was used to search for extragalactic gamma-ray sources.

4.4

Flux calculation

The differential photon flux for a given source is defined as the number of photons received
on the earth from the source per unit energy per unit area per unit time. For a power law
(of index α) spectrum without a cutoff the differential photon flux can be written as:

dN
= I◦
dE

E
1 T eV

−α

Photons TeV−1 cm−2 s−1 ,

(4.4)

where I◦ is the flux amplitude (also known as the flux normalization) and α is the power
law index. For a power law spectrum with an exponential cutoff this can be written as:

dN
= I◦
dE

E
1 T eV

−α − E
e Ec Photons TeV−1 cm−2 s−1 ,

(4.5)

where Ec is the cutoff energy. A plot of differential photon flux versus energy (sometimes
frequency or wavelength) is called the Spectral Energy Distribution (SED).
The Milagro detector does not observe the photons with a constant efficiency across a
46

broad energy band. Therefore, the number of electrons measured by Milagro per second is
the convolution between the Milagro effective area,Aef f (E, θ), and the true SED:

Nobserved =

dN
dθdE
Aef f (E, θ)
dE

(4.6)

The integral of Equation 4.6 over energy and the Milagro detector’s lifetime gives the expected excess that can be expected from the Milagro detector for a given SED:

NE =

dN
Aef f (E, θ(t))
dE dt,
dE

(4.7)

where NE is the expected excess for a given energy spectrum.
In the Milagro data analysis framework, the differential photon flux of an observed source
is derived from the ratio of the observed excess NO to the expected excess NE . The NO
is obtained from a sky map made with an assumed SED, and the NE is obtained from
Monte Carlo simulations with the same spectral assumption. These two parameters can be
combined to derive the flux amplitude of the observed source:

I◦ =

NO
× I◦ Simulated ,
NE

(4.8)

where I◦ Simulated is the flux amplitude used in Monte Carlo simulations to obtain NE .
Once I◦ is derived, the differential flux at any energy can be derived using equation 4.4 or
4.5. Since this algorithm assumed an SED to obtain NO and NE , the calculated flux has a
bias towards the assumed spectrum. This bias for three simulated sources, at a declination
of 22, is shown in Figure 4.2. In this figure, I derived the flux in the range of 1 to 100 TeV
of three simulated sources with α = −2.2, −2.6 and -3.0. As it can be seen in this figure,
47

dN photons TeV-1cm-2s-1
dE

α = -2.2
α = -2.6
α = -3.0

10-12

-13

10

10-14

10-15

-16

10

1

10

E TeV

Figure 4.2: Simulated dN as a function of energy for a pure power law spectrum is shown,
dE
dN − I (E/TeV)−α . I ’s are selected to give the same excess for three sources with three
◦
◦
dE
different spectra, green α = 2.2, red α = 2.6 and blue α = 3.0. Flux for all three sources
are calculated using spectral optimization of α = 2.6. The region where the lines intersect
indicates the energy range least dependent on the assumed spectrum.
the dependency of the calculated flux on the true spectrum is minimal between 20 and 30
TeV. The energy region with the minimum energy dependence is a function of the source
declination, and varies between ∼ 20 TeV and ∼ 50 TeV for galactic sources. Therefore, the
Milagro collaboration chose to publish the flux at 35 TeV for galactic sources and 7 TeV for
extragalactic sources.
The work presented in this thesis also used the integral energy flux over an energy range:

F =

E2
NO
× I◦ Simulated ×
E · E α dE TeV.
NE
E1

(4.9)

which can be determined numerically once the flux amplitude is derived.

4.5

Energy Reconstruction

The algorithm presented in the previous section can be used to derive the differential photon
flux at a given energy. However, that algorithm cannot be used to derive the SED in a given
48

energy range. Milagro used a parameter called FRASOR (defined in equation 4.3) as an
energy estimator to obtained the SED of the Crab pulsar wind nebula using the forward
folding technique[Abdo et al.(2012)]. The forward folding technique is only capable of identifying the best SED, that describes the data, among a set of assumed SEDs. Therefore, the
forward folding technique depends on how well the true SED can be guessed. In Abdo et
al.(2012) a power law with an exponential cutoff was assumed,

dN
=
dE

I◦
1 TeV

−E
E −α e Ec .

(4.10)

Then a sample of expected FRASOR distributions was produced using the Milagro Monte
Carlo simulations, by changing the α, (I◦ ) and Ec . In Milagro technical jargon, the table
that summarizes the results from these simulations is called the Gamma-table. More details
about the Gamma-table are discussed in Appendix . In the next step each expected FRASOR
distribution was compared with the observed FRASOR distribution using the χ2 test. The
significance contours of the χ2 test in α − Ec space and α − I◦ space is shown in Figure 4.3.
These two significance contour plots led the authors of Abdo et al.(2012) to conclude that
the SED of the Crab Pulsar Wind Nebula can be best described by:





+0.7
−12
2.5
× 10
E
dN
 E 2.5±0.4 exp −
 cm−2 s−1 TeV−1
=  −0.4
+39
dE
1 T eV
32−18 TeV

(4.11)

The accuracy of this spectrum was limited for three reasons: poor energy resolution,
poor linear correlation between FRASOR and log10 (true energy), and the limitations of the
forward folding technique. FRASOR is derived from two parameters: the number of Air
Shower layer PMT hits and the number of Outrigger tank hits. These two parameters are
49

Figure 4.3: Two significance contour plots in the α − Ec and α − I◦ space obtained from
forward folding are shown[Abdo et al.(2012)]. The spectral assumption is a power law with
an exponential cutoff, as shown in equation 4.10. The figure on the left shows the 1σ and
2σ contours in the α − Ec space. The figure on the right shows the 1σ and 2σ contours in
the α − I◦ space.

highly sensitive to the size of the EAS. However, the size of the EAS depends not only on the
energy of the initial gamma-ray photon but also on its zenith angle and shower core location.
FRASOR is not sensitive to the zenith angle and the shower core location, these are limiting
Milagro’s energy resolution. I developed an improved energy estimator by combining the
zenith angle, shower core location and some other variables available in Milagro reconstructed
data. The performance of my energy estimator is discussed in Appendix .
The main limitation of the forward folding technique is that it picks the SED that best
describes the observed data from a set of known SEDs, instead of deriving the true SED from
the data. Therefore, picking the correct sample of SEDs is crucial, because it is not possible
to generate the expected FRASOR distribution using all possible SEDs. The expected
FRASOR distributions are generated from a set of discrete parameters for a limited number
of SED functional forms. An alternative method to overcome this issue is unfolding. In my
thesis work I applied the unfolding technique to get the true SED shape of a source observed
by Milagro. However, it turned out that the energy resolution of the Milagro detector, even
with an improved energy estimator, was insufficient to apply unfolding algorithms to the
50

Milagro data set.

51

Chapter 5

The Fermi Large Area Telescope
(Fermi-LAT) GeV Catalogs

The Fermi Gamma-ray Space Telescope (formerly known as the Gamma-ray Large Area
Space Telescope, GLAST) was launched by NASA on June 11, 2008 on a Delta II heavy
launch vehicle. The Fermi space telescope consists of two on-board instruments: the Gammaray Burst Monitor (GBM) and the Large Area Telescope (LAT). Fermi-GBM monitors
Gamma-ray bursts in the energy range of 8 keV to 40 MeV, and the Fermi-LAT surveys
the sky in the energy range of ∼ 20 MeV to ∼ 300 GeV. The work presented in this thesis
used two catalogs produced by Fermi-LAT: the Fermi Large Area Telescope Second Source
Catalog, also known as 2FGL [Nolan et al.(2012)] and the Second Fermi Large Area Telescope Catalog of Gamma-ray Pulsars, also known as 2PC [Abdo et al.(2013)]. This chapter
gives a brief introduction to Fermi-LAT, 2FGL and 2PC.
52

5.1

Fermi Large Area Telescope (Fermi-LAT)

The Fermi-LAT consists of 16 precision converter-trackers and 16 calorimeters arranged in a
4 × 4 array. A schematic diagram of Fermi-LAT is shown in Figure 5.1. An internal view of a
precision converter-tracker and an incoming gamma-ray is shown at the center of this figure.
Each precision converter-tracker consists of alternating layers of tungsten conversion foil and
silicon strip detectors, which have two layers of strips oriented perpendicular to each other.
At the bottom of each precision converter-tracker is a cesium iodide (CsI) calorimeter.
When a gamma-ray enters a precision converter-tracker it interacts with one of the conversion foils and produces an electron-positron pair. The resulting electron-positron pair is
traced by the Silicon Strip Detectors in successive layers. The incidence angle of the initial
gamma-ray can be constructed using the traces of the electron and positron. The cascading
of a gamma-ray is illustrated in Figure 5.1 using a red-dotted line and two blue-dotted lines.
After passing through the precision converter-tracker the electron-positron pair enters
the calorimeter, which consists of 96 optically isolated CsI crystals of size 2.7 cm × 2.0 cm
× 32.6 cm. These crystals are arranged horizontally in eight layers of 12 crystals each with
each layer rotated 90◦ with respect to its adjacent layers. When the electron or positron
passes through a crystal, it converts its the electron’s or positron’s energy into scintillation
light, which is measured by two PIN diodes mounted on each end of the crystal. The sum
of the light deposited in the PIN diodes provides a determination of the energy deposited
on the crystal. The difference in light levels deposited in the two PIN diodes provides a
determination of the position of the energy deposition along the crystal. In addition to these
two measurements, the physical location of the crystal in the calorimeter provides two more
coordinates of the location of the energy deposit. These measurements can be combined
53

to measure the energy of the electron-positron pair produced by a gamma-ray, which can
be used to reconstruct the energy of the primary photon, and to reconstruct the shower
development profile.
The dominant background for the Fermi-LAT is high-energy charged cosmic-ray particles.
In order to identify these particles Fermi-LAT was covered with an Anti-Coincident Detector
(ACD). The ACD is made from plastic scintillation tiles, and the light from each tile is
collected by wavelength shifting fibers that are coupled to two photo multiplier tubes. When
a charged cosmic-ray passes through the tiles it produces scintillation light, while a gammaray can pass through the tiles without interaction. Therefore, charged cosmic-ray particles
produce coincident events both at the ACD and in the precision converter-trackers or in the
calorimeters, while gamma-rays produce a signal only in the precision converter-trackers and
in the calorimeter. This feature has been used to identify the charged cosmic-ray events.

5.2

Fermi-LAT Second Source Catalog (2FGL)

The Fermi-LAT Second source catalog, which is known as 2FGL, is a catalog of highenergy gamma-ray sources detected by Fermi-LAT in the first two years of its operation
[Nolan et al.(2012)]. This catalog includes 1873 sources that were observed in the 100 MeV
to 100 GeV energy range during the period from August 4, 2008 (15:43 UTC) to August 1,
2010 (01:17 UTC). During this period, Fermi mostly operated in sky-scanning survey mode,
which evenly scans the full sky by rocking north and south of the zenith on alternate orbits.
With this rocking motion the corresponding exposure at all sky locations are relatively uniform. However, the exposure is minimum at the celestial equator and it is maximum at the
north celestial pole. The ratio of maximum exposure to minimum exposure depends on the
54

Figure 5.1: Schematic diagram of Fermi-LAT. A precision converter-tracker and calorimeter
is shown at the center of the figure. Here the calorimeter is taken apart from the precision
converter-tracker. The red doted line indicates the track of a gamma-ray and the blue
dotted lines indicate the tracks of the cascade electron and positron. The physical size of
this detector is 1.8 m × 1.8 m × 0.72 m. It weighs 2789 kg and in normal operation mode
consumes 650 W. This figure was obtained from the publication [Atwood et al.(2009)].

55

Figure 5.2: The Aitoff projection of the significance sky map derived from the data set
used to make the Fermi Large Area Telescope Second Source catalog. This map includes
the gamma-ray energy flux measured in the range between 100 MeV and 10 GeV. Refer to
[Nolan et al.(2012)] for the original figure and a more detailed description of data. Color
scale shows the gamma-ray significance of each sky location.

rocking angle, which was 35◦ at the beginning of the data collection and later changed to
50◦ . With a 35◦ rocking angle the ratio of maximum exposure to minimum exposure was
1.33 in contrast to 1.75 for rocking angle 50◦ . The sky map of the energy flux derived from
this data set is shown in Figure 5.2.
Among the sources reported in this catalog, ∼68% are sources off the Galactic plane (
galactic longitude > |10◦ |) that have no association with known galactic sources. From these
extra galactic sources 709 are located within Milagro’s sky coverage (7◦ < DEC < 80◦ ). This
provided a targeted list of extragalactic candidates for Milagro to search for TeV gammaray emission associated with extragalctic GeV sources. The results from a Milagro targeted
search on this source list is presented in Chapter 8.
56

5.3

The Second Fermi-LAT Catalog of Gamma-ray Pulsars (2PC)

The second Fermi-LAT Catalog of Gamma-ray Pulsars (2PC) is a catalog of 117 GeV-pulsars
measured by Fermi-LAT in the energy range of 0.1 GeV - 100 GeV, using data collected
from August 14, 2008 to August 14, 2011[Abdo et al.(2013)]. Similar to the data collection
period for 2FGL, Fermi mostly operated in the sky-scanning mode during this period. They
excluded the data when Fermi-LAT was not operating in the science operations mode, the
rocking angle exceeded 52◦ , or the gamma-rays were detected with a zenith angle greater
than 100◦ . This data set was searched for GeV pulsars using these different methods: a
targeted search using known pulsars measured in radio or X-ray energies, a targeted search
with a candidate list made from known GeV point sources (e.g. SNR), and a blind search.
The first method searched for GeV pulsation in the LAT data using 2286 known pulsar
ephemerides and identified 60 GeV pulsars. This is the most sensitive method for observing
pulsars with low gamma-ray fluxes. However, this method can not identify radio and X-ray
quiet pulsars. The other two methods were able to discover 57 new GeV pulsars.
From these 117 GeV pulsars, 32 were located within Milagro’s sky coverage (7◦ < DEC
< 80◦ ). I performed a Milagro targeted search on this list of 32 pulsars and the results are
presented in Chapter 9. The analysis presented in Chapter 9 used the pulsar period (P ), the
first derivative of the pulsar period P˙ , the spin-down luminosity E˙ , the phase averaged
integrated energy flux in the 0.1 to 100 GeV energy band (G100 ) and the distance to the
pulsars, all of which are readily available in the 2PC.
The 2PC obtained the P and P˙ of a pulsar from observed data, and E˙ was derived from
57

the P and P˙ using equation 5.1.
E˙ = 4πI◦ P˙ /P 3 ,

(5.1)

where I◦ = 1045 g cm2 . Full derivation of this equation is given in Section 2.2.2. The
integrated energy flux in the 0.1 to 100 GeV energy band was derived by integrating the
Spectral Energy Density (SED). The SEDs were obtained by forward folding a power law
with an exponential cutoff,

dN
= I◦
dE

E
E −α
exp −
E◦
Ecut

.

(5.2)

The pivot energy E◦ was arbitrarily defined as 1 GeV for pulsars unassociated with 2FGL
sources. For the pulsars associated with 2FGL the E◦ listed in 2FGL was used. The I◦ , α
and Ecut are the input parameters, and the best fit values of these input parameters were
obtained using a likelihood analysis method. Once the best SED was identified, the G100
was derived by integrating the SED,
100 GeV
G100 =

5.4

E
100 M eV

dN
dE
dE

(5.3)

Source Naming Convention

Both the 2FGL catalog and the second pulsar catalog use the same naming convention,
XXXX JHHMM.m+DDMM[c]. The first set of characters, XXXX denotes the name of the
catalog. In the 2FGL catalog XXXX is 2FGL, and in the second pulsar catalog XXXX is
PSR. In the 2FGL catalog ‘2’ refers to the second-year catalog and ‘FGL’ represents Fermi
Gamma-ray LAT. HHMM.m is the right ascension written in hours and truncated to 0.1
58

decimal minutes. DDMM is the declination written in degrees and truncated to 1 degree.
The suffix ‘c’ is optional, which indicates that the source is considered to be potentially
confused with Galactic diffuse emission.

59

Chapter 6
The High Altitude Water Cherenkov
Observatory (HAWC)
The High Altitude Water Cherenkov observatory (HAWC) is the successor to the Milagro
observatory. HAWC is currently under construction near the Volcano Sierra Negra, Mexico
at latitude 180 59 48

North, longitude 970 18 34

West and altitude 4100m. As a part

of my thesis work I developed the firmware for HAWC’s GPS Timing and Control (GTC)
system. In this chapter I briefly present the design of the HAWC observatory and the working
principles of the GTC system.

6.1

The HAWC Detector

The HAWC detector is designed on the same principles as Milagro. The primary difference
between Milagro and HAWC is that instead of a large water pond HAWC consists of densely
packed steel water tanks of 7.3 m in diameter and 4.5 m in height. Each of these tanks holds
a light-tight bladder filled with about 200,000 liters of ultra purified water, and 4 PMTs
60

Figure 6.1: A drawing of a HAWC tank with 4 PMTs at the bottom.
are placed near the bottom of the bladder. A drawing of the tank and the 4 PMTs at the
bottom is shown in Figure 6.1. The PMT at the center is a Hamamatsu 10 high quantum
efficiency, 7081HE PMT, and the other three PMTs are 8 Hamamatsu R5912 PMTs reused
from Milagro. HAWC also reused the front-end electronics boards from Milagro. CAEN
VX 1190A TDCs are used to digitize the TOT signal outputs of the front-end electronics
boards. Each of these TDCs can digitize 128 input signals. Currently HAWC has 111 tanks
operational, which uses 4 TDCs. When HAWC is completed in 2014, it will have 300 tanks
and use 9 full TDCs and 48 channels from a 10th TDC to digitize the PMT signals. Other
than the 10 TDCs used to record PMT signals, an 11th TDC will be used to record 32 signals
coming from the GPS Timing System. These signals are similar to the TOT signals, but
they are encoded with the current GPS time, which is the time-stamp of an event. Further
details of the time-stamps are discussed in section 6.2.2.
CAEN Vx 1190A TDCs are designed to record TOT measurements within a given time
window around a trigger signal, and are equipped with an output buffer to store the recorded
61

data until they are read out. As each TDC buffer fills with data, a GE XVB602 Intel Core
i7 based VME Single Board Computer (SBC) reads each individual TDC in parallel and
delivers the data to the online reconstruction farm. However, each SBC cannot perform the
readout process at the same rate. Therefore, the online reconstruction farm receives different
fragments of a single event at different times. The HAWC online reconstruction software
identifies the event fragments belonging to a given trigger using the event identification
number (Event ID), which is a 12 bit number that is stamped on each event by the TDCs.
After identifying the fragments of a single event, the online reconstruction process combines
the fragments into a single event and decodes the time-stamp into a computer readable data
type. This defragmentation is possible only if all the TDCs are working synchronously and
maintain an unique Event ID for a given trigger. This was achieved by the control part of
the GTC system that controlled the TDCs through its control bus.

The controlling of the TDCs is done using three signals of the TDC control bus: TRG,
CLR and CRST. The TRG is the trigger signal to the TDCs. In HAWC, the trigger signal
is a periodic signal that is provided by the GTC system. In a typical data run the periodic
trigger frequency is 40 kHz (period = 25 µs) and the TDCs record the data in a 25.2 µs
window around each trigger. The data saved in a given time window is called an event. The
CLR is the clear command, which clears the data from the output buffer, resets the event
counter, and performs a TDC global reset. The CRST is the reset command, which performs
the same actions as the CLR, in addition to resetting some of the other internal counters of
the TDC. Further details of the interface between the control bus and the GTC system are
discussed in section 6.2.3.
62

6.2

GPS Timing and Control System (GTC System)

Because of its central role, the GTC system is considered the central nervous system of
HAWC. As its name suggests the GTC system combines two sub-systems: the GPS Timing
system and the Control system. The primary tasks of the GPS Timing system is to provide
a time-stamp for each recorded event, which is the absolute time of the trigger, and to derive
a low jitter 40 MHz signal to use as the global clock signal for HAWC. The primary task
of the Control system is to provide the clock and the control signals to the TDCs and to
provide trigger and detector status information to the scaler system.

6.2.1

GPS Timing System Hardware

The GPS Timing system consists of a GPS receiver and a customized piece of electronics
called the Clock type H Clock card. A functional block diagram of the GPS Timing system
is shown in Figure 6.3. The NAVSYNC CW46 GPS receiver is a commercial GPS receiver
that carries three signals to the Clock type H Clock Card. H Clock cards are general-purpose
custom cards, and Clock type refers to a version of the H Clock cards. A photograph of a
fully assembled H Clock card is shown in Figure 6.2.
An H Clock card is a 2 unit wide 6U VME-64X module that is equipped with a Phase
Lock Loop (PLL), 10 LVDS type General Purpose Input Output (GPIO) connectors, a 16
pin connector to the GPS receiver, a A24D16 VME interface and a Virtex II FPGA on a
mezzanine board. In Figure 6.2 the ten GPIO connectors are on the left. Each of these ports
has 17 signal pairs; 16 of them are LVDS GPIO signals to/from the FPGA, and the 17th
pair is a 40 MHz clock signal which is also a LVDS type signal. The direction of the GPIO
ports are reversible by replacing the IO driver chips. Clock type H Clock cards are created
63

Figure 6.2: A photograph of a fully assembled H Clock card is shown. The Clock version
of the H Clock card with 2 general purpose input ports and 8 general purpose output ports
are used in the Clock System. The Control version of the H Clock card used in the Control
System has 4 general purpose input ports and 6 general purpose output ports.
by configuring H Clock cards with two input ports and eight output ports. The FPGA is
mounted in a mezzanine card (on the right in the picture). Since the performance and the
resources of the Virtex II family FPGAs are adequate for the requirements of HAWC, a
Virtex II xc2v1000-4fg456 FPGA is used. A stock of this FPGAs was also available at MSU.
However, if a future upgrade needs to change the FPGA, it can be done by making a new
mezzanine card.
The GPS receiver is used to obtain the GPS time and a low jitter 10 MHz sine wave
signal. The internal PLL of the Clock type H Clock Card uses this 10 MHz sine wave signal
to derive a low jitter 40 MHz digital clock signal and makes several exact copies that are
delivered to the Control type H Clock Card, to the FPGA and to the 17th signal pair of all
the GPIO connectors. This 40 MHz signal is used as the global clock signal for HAWC. The
GPS receiver also transmits a one pulse per second (1PPS) pulse stream and a set of data
strings via the RS232 protocol. The rising edges of these 1PPS pulses mark the top of each
64

second. The firmware running inside the FPGA uses this 1PPS signal and the data strings
to replicate the current GPS time.

GPS Receiver
NAVSYNC
CW46

40 MHz

H_Clock 01: Clock Card

1 pps

(All LVDS I/O)

10 MHz
1 pps
ASCII Time
Control
Power

H_Clock: Control Card

10 MHz
40 MHz
PLL

Time Code Pulses (32 digit BCD)

Xilinx Virtex FPGA

Scaler System
(WBS 4.7)

32

Time Code Pulses (32 digit BCD)

TDC System
(WBS 4.3)

32

VME Interface

VME Bus
NTP Time
Server

VME to PCIe
Bridge

Run Control
Computer

NTP Time Service Signal

Figure 6.3: A block diagram of the Clock System is shown. The NAVSYNC CW46 GPS
receiver and the Clock Card are the two main components of the Clock system. The GPS
receiver is used to obtain the GPS time and a low jitter 10 MHz signal. The Clock Card
produces a 40 MHz global clock signal and two time-stamps for the scaler system and for
the TDC system. The communication between the Clock Card and the control computer is
accomplished through a A24D16 VME interface.

6.2.2

GPS Timing System Firmware

The logic operations of the GPS Timing system are done by the firmware running in the
FPGA. Part of my PhD thesis work was to develop, test and implement this firmware.
A simplified functional block diagram of this firmware is shown in Figure 6.4. This is a
65

sequential logic design with several state machines implemented using VHDL. In this section
I briefly discuss the functionality of this firmware.
1PPS
GPS Time

RS232 String

Serial to Parallel
Converter

Receiver's
Health

Encoded Timestamp

Time

Internal Clock

28

Error Code
4

Timestamp
Encoder

28

Encoded Error Code
4

VME Module
A24 D16
RS232 String

GPS Configuration
Strings

Figure 6.4: A simplified functional block diagram of the Clock firmware is shown. This
firmware maintains an internal clock synchronized with the GPS time and produces TDC
readable time-stamps. The communication between this module and the control computer
was done by a A24D16 VME interface.

The first module of the firmware, the Serial to Parallel Converter, reads the serial data
strings coming from the GPS receiver and extracts the current GPS time and the GPS
receiver’s health information from the GPS signals. The GPS receiver installed at the HAWC
site is configured to send three data strings followed by a 1PPS pulse. These three strings
(POLYT, GPGSA and POLYP) are standard NMEA 0183 strings, which carry the current
GPS time, GPS receiver operating mode, satellites used to derive the time, and dilution of
precision (DOP) values. The Serial to Parallel Converter module extracts the current GPS
time from the POLYT string, GPS Fix status from the GPGSA string1 , and dilution of
precision measurements from the POLYP string. Then the GPS time and health information
goes to the Internal Clock module, which is a continuously running 8 digit binary-coded
1 GPS Fix is an integer between 1 and 3. If the GPS receiver is not getting enough GPS
signals and it is unable to fix, GPS Fix is 1. If GPS receiver is being able to become to a
2D or 3D fix this GPS Fix becomes 2 or 3 respectively. Refer to the CW25 GPS Receiver
User Manual for more information.
66

decimal (BCD) clock that runs using the 40 MHz clock signal as the reference frequency.
The least significant digit of the internal clock is microseconds and the most significant
digit is tens of seconds. Apart from the signals coming from the Serial to Parallel Converter
module, the Internal Clock module also receives the 1PPS signal that comes from the GPS
receiver, which is used to identify the top of each second. At the top of each second, the
Internal Clock module compares its clock time with the GPS clock time, and overwrites the
internal clock if the time does not match (and the GPS receiver is in good health). This
allows the GPS Timing System to have an internal clock that runs synchronously with the
GPS clock.
The second stage is the Timestamp Encoder module, where the time of the internal
clock is encoded into a 32 bit TDC readable time-stamp and then transmitted to the TDCs
every ∆T µs interval (∆T µs can be configured to 10µs or 20µs). The first 28 bits of the
time-stamps are 7 BCD-digits that carry the time in the format: 10’s of s, 1’s of s, 100’s
of milliseconds, 10’s of milliseconds, 1’s of milliseconds, 100’s of microseconds and 10’s of
microseconds. The last four digits of a time-stamp are the error code of the time-stamp,
which is described in Table 6.1. The encoding of this time-stamp to a TDC readable format
is done using a simple algorithm. Each bit is denoted by a pulse; a 1 µs wide pulse denotes
a logic zero bit, a 2 µs wide pulse denotes a logic one bit. As an example, the timing
diagram shown in Figure 6.5 is the encoding for the time 12.34567 seconds with the error
code 0000. This TOT-like encoding allows HAWC to use the same software to read both the
FEB outputs and the time-stamps.
Besides these major modules, the Clock Firmware has a VME module that provides an
A24D16 VME interface and several FIFOs. These FIFOs are filled with the information that
users need to monitor the health of the GPS Timing System.
67

Figure 6.5: Encoding of the time-stamp 12.34567 seconds with no errors is shown. Each
digit is encoded into a 4 bit binary number, and each binary number is then encoded into a
pulse. A 1µs wide pulse is used to indicate logic 0 and, a 2µs wide pulse is used to indicate
logic 1.

68

Error Code
xxx1
xx1x
x0xx
1xxx

Error
Internal Clock was over written by the current GPS clock.
GTC system lost communication with the GPS receiver.
GPS receiver do not have enough satellite reception to make a GPS Fix.
NMEA string coming from the GPS receiver has errors.

Table 6.1: The error status corresponding to each time-stamp is encoded into the 4 most
significant bits of the time-stamp The meaning of each error code is shown.

6.2.3

Control System Hardware

The Control System is made from two custom VME boards to provide several services to
HAWC: 1) keep all the TDCs working synchronously , 2) issue a synchronous trigger signal
to the TDCs, and 3) issue a trigger signal and send the status of the detector to the scaler
system. The functional block diagram of the Control system is shown in Figure 6.6.
The H Clock: Control card and CB Fan card shown in Figure 6.6 are the two custom
VME boards. The H Clock card used in the Control system (Control type H Clock) is a
version of the H Clock card with 6 input ports and 4 output ports. The CB Fan card is
a 2 unit wide 6U VME-64X module that is designed primarily to provide appropriate level
conversions and fan-outs for the interface between the Control type H Clock Card and the
TDCs. A photograph of a fully assembled CB Fan card is shown in Figure 6.7.
In contrast to the Clock type H Clock card, the Control type H Clock card gets the 40
MHz global clock signal from the Clock card. Then the Control type H Clock Card makes
several copies of this 40 MHz clock signal, and distributes them to the 17th signal pair of
the 10 GPIO connectors and to the FPGA. One input GPIO port and one output GPIO
port of the Control type H Clock card is connected to the scaler system. The other GPIO
ports are connected to two CB Fan cards.
The interface of the Control type H Clock Card to the Scaler system consists of 4 output
signals and one input signal: 10 MHz reference, Pause Pulses, Busy Pulses, LNE, and LNE
69

Analog Trigger
Scaler System (WBS 4.7)

Trigger System
(WBS 4.4)

Paused Pulses x 10 MHz
Busy Pulses x 10 MHz

2

H_Clock: Control Card

LNE

2
(All LVDS I/O)

40 MHz, LNE,
Busy (ECL)

H_Clock: Clock Card

TDC

LNE Enable

3
ECL to LVDS
40 MHz
CLR
CRST
TRIG
Almost Full 01
Almost Full 02

Busy

40 MHz

CB_Fan x 2

LNE

1 pps
1

40 MHz
CLR

Adjustable
Periodic
Trigger

......................

2
Fanout

CRST

LVDS to Diff. ECL 110 Ω
TRIG

6 (8)
6 (8)

4

16 Control
Buses

Xilinx Virtex FPGA

12

Almost Full 01
Distribute or
Aggregate

16
Trigger Signal (LVDS)

15
16

6 (8)
6 (8)

Almost Full 02
16

VME Interface

TDC System
(WBS 4.3)

16
Trigger System
(WBS 4.4)

4

Readout

12

Power Only

VME Bus

VME to PCIe
Bridge

Run Control
Computer

Figure 6.6: A block diagram of the Control System. The Control system consists of two
VME cards: the Control Card and CB Fan card. The Control card does all the logical
operations and the CB Fan card does the appropriate level conversion to interface TDCs to
the Control Card.
Enable input. The 10 MHz reference is a continuous 10 MHz square wave signal output.
The Pause Pulses produce a 10 MHz signal in-phase with the 10 MHz reference when the
Control system is in the pause state. The scaler system counts both of these signals. The
ratio of the number of pause pulses to the number of 10 MHz reference pulses gives the dead
time of the detector. The functionality of the Busy Pulse output is similar to the Pause
Pulses except Busy Pulses produce a 10 MHz signal when at least one TDC is almost full.
The LNE (Load Next Event) is the trigger signal for the scaler system, which is a periodic
signal.
The interface of the Control type H Clock Card to the CB Fan cards consist of four
70

Figure 6.7: A photograph of the fully assembled CB Fan card is shown. A CB Fan card is
able to provide the level conversions and Fan-outs required to handle 6 TDCs.

output signals and 16 input signals per card: 40 MHz, CLR, CRST, TRIG and 16 Almost
Full. The Control card makes four identical copies of all its output signals and can accept up
to 64 inputs. Therefore, one Control type H Clock Card can be connected to four CB Fan
cards. The four output signals (40 MHz, CLR, CRST, TRIG) and the input signals (Almost
Full) of the Control type H Clock Card are the control bus signals of the TDC. But the
TDCs can not directly connect to the Control type H Clock Card, because these I/Os are
LVDS signals but the TDC control bus is compatible with only ECL signals. Therefore, the
CB Fan card is designed to provide the level shifting from Control type H Clock Card LVDS
outputs to ECL and TDC ECL outputs to LVDS. Apart from the level shifting, the CB Fan
card makes 6 identical copies of the 40 MHz,CLR, CRST and TRIG signals. Therefore, one
CB Fan card can be used to control up to 6 TDCs. Since one Control type H Clock Card
can interface with 4 CB Fan cards, the GTC system is capable of controlling up to 24 TDCs.
Each CB Fan card also fan-out three copies of the LVDS signal type control buses. These
71

control buses are not used in the current design.

6.2.4

Control System Firmware

A simplified functional block diagram of the Control firmware is shown in Figure 6.8. Similarly to the Clock firmware, this firmware is a sequential logic design implemented in VHDL.
However, unlike the Clock firmware, individual modules of the Control firmware are not
connected sequentially. Coordination of these modules are done through the VME module.

The first module shown in Figure 6.8 is the Trigger module, which coordinates the trigger
signals that go to the TDCs. The trigger module can work in three modes: pause, periodic
trigger and external trigger. In the pause mode, the trigger module is not issuing any
triggers. In the periodic trigger mode, the trigger module issues a periodic trigger signal
with a known frequency that is set by the VME module. In the external trigger mode, the
trigger module issues a trigger signal upon a request coming from the external trigger. This
trigger mode so far not has been used in HAWC. In a typical HAWC data taking run, the
trigger module runs in the periodic trigger mode with a trigger frequency of 40 kHz. At the
end of each run, the HAWC experiment control system sends a request to the VME module
to switch the trigger module to the pause mode. The 40 kHz periodic trigger frequency was
chosen because it is the optimum trigger frequency that could be handled by the HAWC
DAQ system. However, this periodic trigger frequency can be changed by the run control
system sending a request to the VME module. The CLR and CRST modules issue the clear
and reset signals to the TDCs upon a request coming from the VME module. These requests
come to the VME module from the run control system at the beginning of each run.
The next three modules in Figure 6.8 provide signals to the scaler system. The func72

Figure 6.8: A simplified functional block diagram of the Control firmware is shown. The
Trigger, CLR, and CRST modules generate the TDC control signals, and LNE, Pause Pulses,
Busy Pulses and 10 MHz Reference produces signals to the scaler system. Similar to the
Control firmware, communication between the Control system and the control computer was
done by a A24D16 VME interface.

tionality of the Pause Pulse module is equivalent to a multiplexer with two inputs and one
output: Logic Lo input, 10 MHz square wave input and Pause Pulse output. When the
Trigger module is in the Pause state the Pause Pulse output switches to the 10 MHz square
wave. When the Trigger module is not in the Pause state, the output switches to the Logic
Lo level. The Busy Pulse module has a similar functionality, except selection between Logic
Lo and 10 MHz square is done using the OR of the Almost Full signals. If any of the Almost
73

Full inputs are Logic Hi, the Busy Pulse output gets connected with the 10 MHz square
wave, otherwise the out put stays in the Logic Lo level. The 10 MHz Reference module is a
10 MHz square wave signal generator that generates a reference pulse stream to the scaler
system.

6.3

Performance of the GTC system

The fully functional GTC system was integrated with the HAWC 30 array (30 tanks) in
October, 2012. The operation of HAWC 30 with 4 PMTs in each tank needed only 120 TDC
channels (30 × 4) to record PMT signals. Therefore, 1 TDC was enough to record the PMT
signals, and 32 channels from a second TDC was used to record the time-stamps.
At the beginning of HAWC 30 operations we were concerned about getting missing edges
of the time stamps. Therefore, we configured the GTC system to issue a periodic timestamp every 10 µs, and the TDCs to record the data in a 25.2 µs time window. In this
configuration two time-stamps are guaranteed for each trigger window. After analyzing data
from 8 months we did not find any time-stamps with missing edges. Therefore, the frequency
of time-stamps was reduced to one time-stamp per 20 µs. With these settings, only one timestamp is guaranteed for each trigger window. This new configuration reduces the band width
of the TDC readout by 50%.
The health of the GTC system has been continuously monitored since early 2013, and it
has revealed that the 1PPs signal has a jitter of less than 50 ns with respect to the 10 MHz
output. Each time the 1PPs signal jitters more than 25 ns the GPS clock over writes the
GTC internal clock, and the firmware produces an error flag in the time-stamp. The average
rate of this error flag was 11 per hour. This jitter introduced an upper limit of 50 ns accuracy
74

to the GTC generated time-stamp. The 25 ns accuracy is well below the required accuracy
of 1µs for HAWC. The Clock firmware was upgraded in October 2013 to be immune to this
jitter.
The GTC health monitoring system is also continuously monitoring the GPS Fix status
and Time Dilution of Precision. During the period from October 2012 to November 2013
the GPS Fix was lost only once. The GPS receiver was able to recover the lost GPS Fix
within 5 minutes. The internal clock was not over written by the GPS time after the GPS
receiver recovered from the lost GPS Fix. This observation yields the conclusion that the
internal clock of the GTC system was able to maintain its internal clock during that 5
minuets without the assistance of the GPS receiver. During the period from October 2012
to November 2013 the Time Dilution of Precision (TDOP) was in the excellent range (TDOP
< 2) for 95% of the time and was in the good range (2 < TDOP < 5) 4.9% of the time. There
were few occurrences, 0.1% of the time, with moderate (5 < TDOP < 10) Time Dilution of
Precision.

6.4

Additional Capabilities of H Clock cards

Since H Clock cards are made with a re-programmable FPGA, H Clock cards can be used
for different applications. As a part of my service work projects, I developed the HAWC
pattern generator and Scaler Time Stamp reader using spare H Clock cards. The HAWC
pattern generator was made from a Clock type H Clock card. This pattern generator is able
to produce TOT-like signals that are similar to the output signals of the front-end electronic
boards. A key feature of this pattern generator is that the patterns are not hard-coded. The
patterns can be configured through special software.
75

The Scaler Time Stamp reader was made from a Control type H Clock card. The purpose
of the Scaler Time Stamp reader is to read time-stamps produce by the GPS Timing system
and store them in a FIFO to be read through the VME back-plane. The Scaler Time Stamp
reader was integrated with the scaler system in May 2013.

76

Chapter 7

Data Analysis Methodology

In the Milagro skymaps, the expected significance at a sky location with no true emission is
a Gaussian random variable with mean 0 and unit standard deviation [Atkins et al.(2004)].
A common treatment of N candidate searches is to use a trials correction technique. Here
one can choose a significance threshold, calculate the tail probability (p-value) λ, and adjust
λ . The purpose of the trials correction is to maintain, at the value
the p-value threshold to N
λ, the probability of a background fluctuation producing one or more false discoveries among
the N searches. The False Discovery Rate (FDR) technique discussed in [Miller et al.(2001)]
offers some advantages over the trials correction technique. In Section 7.1 I briefly discuss
the FDR technique. The FDR technique can be used to search for individual candidates
with a TeV association. A stacking analysis can be used to search for evidence of collective
TeV emission among the undetected candidates by studying their mean flux. In Section 7.2
I briefly discuss the stacking analysis used in this thesis.
77

7.1

False Discovery Rate method

Instead of controlling the expected probability of having even one false detection, FDR
controls the expected fraction of false discoveries among a set of detections; that is, it controls
the contamination fraction of the lists of associations, rather than the probability of a random
individual association being accepted. The key input parameter is again a probability λ, but
now λ represents the expected fractional contamination of any announced set of detections.
Based on this input parameter, the method dynamically adjusts the detection threshold but
in a way that depends on the properties of the entire list of search significances (converted
into p-values). This dynamic adjustment is sensitive to whether the distribution of p-values
is flat (as would be expected if there were no detectable sources) or skewed to small p-values
(i.e. large significances). This adjustment lowers the significance threshold for detection
if a list is a “target-rich environment” in such a way that the expected fraction of false
discoveries among the announced detections remains at the fraction λ. In particular, the
most significant candidate is required to have a p-value of λ/N just as in the trials-correction
method, but the n-th most significant candidate need only have a p-value less than λ × n/N .
As a result, this technique has a higher efficiency for finding real detections, while producing
the same results as a trials-correction method in target-poor environments where the only
decision is whether to report zero or one detection. The method adjusts for both the length
of the search list and the distribution of the significances found within the search list. Note
that as a result, a given candidate location might pass the FDR criteria on one search list,
but fail in another. Also note that λ controls the expected contamination, i.e. averaged over
potential lists of associations, not the contamination fraction on a specific list.
78

7.1.1

Calculation

Calculation of the FDR threshold have 4 main steps. First the parameter λ has to be
selected, such that 0 < λ < 1. The FDR procedure guarantees that the average rate of false
discoveries is less than or equal to λ. The second step is to calculate the correspoinding P
values for each candidate. Let P1 ,..,Pj ,..,PN denote the P-values correspond to N candidates
that are listed from smallest to largest P-value. The thired step is to find the maximum j
value that satisfies,
jλ
},
jmax = max{j : Pj <
cN N

(7.1)

where cN = 1 when the P-values are based on statistically independent tests. If the P-values
are dependent,
N
cN =

i−1 .

(7.2)

i=1
The last step is to reject all hypotheses with Pj < Pjmax .

7.2

Stack Analysis

A stacking analysis can be used to search for evidence of collective TeV emission among the
undetected candidates by studying their mean flux. This thesis used the stacking methodology of Section 3 in [Kurczynski & Gawiser(2010)]. The significance of the stacked flux is
given by Equation 7.3 below.

Significance =
79

I
∨( I )

,

(7.3)

where I is the weighted average flux as defined in Equation 7.4 below and ∨ I is its
variance, defined in Equation 7.5.

I =

Ii
σi2
1
σi2
1
σi2

∨I =

1 2
σi

(7.4)

(7.5)

Here Ii is the flux of each candidate and σi is the standard deviation of flux of each candidate.

80

Chapter 8
Milagro Observations of Extragalactic
Sources
The Milagro collaboration performed blind searches of the Milagro sky maps and found a
number of TeV sources ([Atkins et al.(2004)] and [Abdo et al. (2007)]). Blind searches for
excess events over the full sky have a high probability of picking up random fluctuations.
Therefore after trials correction, a full sky blind search is less sensitive than searches using a
predefined list of potential TeV emitting candidates. Therefore, I performed three targeted
searches for extragalactic sources using the existing Milagro sky maps. This chapter discusses
the results from these three targeted searches.
The first extragalactic candidate list is compiled from the extragalactic sources in the
Fermi-LAT second source catalog (2FGL) [Nolan et al.(2012)]. The second extragalactic candidate list is made from the extragalactic sources in the TeVCat catalog [Wakely & Horan(2008)].
All of the TeVCat extragalactic sources were detected by IACTs. Measurements available
in TeVCat include the transient states of variable extragalactic sources. The analysis pre81

sented in Section 8.2 searched for the time average flux of those candidates by averaging their
flux over approximately 8 years. The third and last extragalactic candidate list is compiled
from a list of potential TeV emitting BL Lac objects that were theoretically predicted by L.
Costamante and G. Ghisellini 2002. The analysis presented in Section 8.3 searched for the
TeV emission associated with these BL Lac objects.
In order to optimize the sensitivity to extragalactic sources, these three searches used a
− E
Milagro sky map made with a power law spectrum with an exponential cutoff, E −2.0 e 5TeV .
This choice reflects the fact that when TeV gamma-rays travel cosmological distances they
are attenuated due to interactions with photons from the extragalactic background light
[Dwek & Krennrich(2005)]. This means that the energy spectrum of extragalactic sources
is cut off at high energies. This spectral assumption is similar to the power law spectral
index and the cut-off energy measured for Markarian 421 and Markarian 501 by VERITAS
[Krennrich et al.(2001)]. The detection thresholds for these searches were obtained using the
False Discovery Rate (FDR) technique, discussed in Section 7.1, with a 1% contamination
fraction.

8.1

Targeted Search for 2FGL Extragalactic Sources

The 2FGL catalog has 1808 sources off the Galactic plane

|l| > 100

that have no as-

sociations with known galactic sources. From these extragalactic sources 709 were within
Milagro’s sky coverage −7◦ < DEC < 80◦ , of which 72% have spatial associations with
blazars. Among these blazars 4 are firmly identified as BL Lac blazars and 12 are firmly
identified as FSRQ1 type blazars.
1 Flat Spectrum Radio Quasar
82

In this list only 2FGL J1104.4+3812 (also known as Markarian 421) is classified as a
source by the FDR procedure with a 1% contamination. The Milagro sky map of a 5◦ × 5◦
region around the Fermi-LAT position of the Markarian 421 is shown in Figure 8.1. This is
a 9.57σ detection and the measured differential photon flux is 3.89 ± 0.4 × 10−15 photons
TeV−1 s−1 cm−2 at 7 TeV. Table 8.2 gives the 95% confidence level flux upper limit for the
brightest 20% in the energy band of 3 to 10 GeV in this list.

8.1.1

Comparison with Previous γ-ray Flux Measurements

Markarian 421 is one of the closest (z=0.030) known, best studied TeV photon emitting
blazars. Markarian 421 was discovered to emit TeV photons by Whipple in 1992 [Punch et al.(1992)].
Since then, Markarian 421 has been observed by several gamma-ray observatories and identified as the source with the fastest observed TeV flux variation [Albert et al.(2007)]. The
gamma-ray emission from Markarian 421 is known to have multiple states: high flux states
and a low flux states. The gamma-ray flux at a high flux state was measured by the H.E.S.S.
observatory between April and May of 2004. The measured gamma-ray flux between 1.5
and 10 TeV was found to vary by a factor of ∼4 on a time scale of days, and the averaged
differential energy distribution was consistent with a power law with an exponential cutoff
of:
dN
= 1.55 ± 0.08stat ± 0.4sys × 10−10 (E TeV)−2.1±0.1stat ±0.3sys
dE

(8.1)

exp −E/(3.1 ± 0.4stat ± 0.9sys TeV
The MAGIC telescope observed Markarian 421 between November, 2004 and April, 2005
[Albert et al.(2007)], when it was considered to be in a low flux state. The measured gamma83

State
Milagro Average Flux
High Flux State
Low Flux State

Differential flux photons ×10−15 TeV−1 s−1 cm−2 at 7 TeV
3.89 ± 0.4
270+170
−110 (H.E.S.S.)
4.98+9.45
−3.69 (MAGIC)

Table 8.1: The gamma-ray flux of Markarian 421 at 7 TeV in the high flux state, low flux
state and the average flux is summarized
ray flux between 100 GeV and 3 TeV was also found to vary by a factor of 4 on a scale of
days, and the averaged differential energy distribution was consistent with a power law with
an exponential cutoff of:

dN
= 0.45 ± 0.01 × 10−10 (E)−2.20±0.08 exp (−E/(1.44 ± 0.27 TeV))
dE

(8.2)

Contrasted with these two measurements, the Milagro flux measurement is the average flux
over approximately eight years. The gamma-ray flux of Markarian 421 at 7 TeV in the
high flux state, low flux state and the average flux is summarized in Table 8.1.1. The
Milagro flux measurement of the average flux over approximately eight years is less than
the flux measured by MAGIC at a low flux state. Therefore, the Milagro time average
flux measurement suggests that Markarian 421 has lower flux states than the low flux state
measured by MAGIC, and Markarian 421 stays in a lower flux state for most of the time.

8.1.2

Undetected Candidates

Milagro also observed a signal excess at the Fermi-LAT sky locations of Markarian 501 (2.9σ),
TXS 1720+102 (2.8σ) and 1ES 0502+675 (2.5σ). Their significances are 2.93, 2.84 and 2.53
respectively, which is insufficient to pass the FDR cut of 1% contamination (λ = 0.01).
Among these three candidates, Markarian 501 and 1ES 0502+675 have been already reported
84

as TeV sources in TevCat. However, TXS 1720+102 has not yet been identified as a source
with TeV emission. This is a radio quasar type blazar identified at a red shift of 0.732
[Afanas’Ev et al.(2005)]. The lowest FDR cut that TXS 1720+102 passes is λ = 0.32. With
this looser FDR cut, three candidates become TeV associations: Markarian 421, Markarian
501 and TXS 1720+102. However, the expected contamination of the resulting candidates
list is 32% so it is likely that TXS 1720+102 is a background fluctuation. While it is hard to
advocate a dedicated IACT observation of TXS 1720+102, a clearer answer to whether this
is a background fluctuation will be provided by HAWC, which is a survey type instrument.
The FDR method is capable of identifying individual candidates which have a TeV association. However, there might be some collective TeV emission coming from the candidates
that fail the λ = 0.01 FDR cut. I searched for evidence of collective TeV emission on the
candidates that fail the λ = 0.01 FDR cut by using the stacking method described in Section 7.2 . I stacked the 2FGL Extragalactic candidates in two different ways: first all FDR
False 2FGL Extragalactic sources and then the FDR False sources among the brightest 20%
in the Fermi energy band 3-10 GeV. These two lists had 0.7σ and 0.6σ significance, respectively. Neither of these stacking results indicate statistically significant collective gamma-ray
emission from the rejected candidates.

8.2

Targeted Search for TeVCat Extragalactic Sources

As of February 8th, 2012 the TeVCat catalog ([Wakely & Horan(2008)], http://tevcat.uchicago.edu)
had 135 sources, of which 31 sky locations were off the Galactic plane and within Milagro’s
sky coverage. These 31 candidates were IACT detections and 23 are identified as BL Lac objects. Among these 31 TeVCat sources only Markarian 421 is classified as a source in Milagro
85

maps with a 1% contamination. Stacking the TeVCat candidates other than Markarian 421,
gives a 0.9σ upward fluctuation. This is a slightly higher upward fluctuation than the 2FGL
extragalactic sources stacking. However, this still does not indicate a significant collective
gamma-ray emission from the rejected candidates.

Among the candidates that fail the λ = 0.01 FDR cut, TeVCat reports Markarian 501
as a well studied, nearby (z=0.034) TeV emitting blazar. Milagro observed a 2.93 σ positive
excess for Markarian 501, and the measured 95% confidence upper limit is 1.8 × 10−15
photons cm−2 s−1 TeV−1 at 7 TeV. Similar to Markarian 421, Markarian 501 is also known
as a blazar with extreme spectral variability [Acciari et al.(2011)]. Most of the studies of
Markarian 501 were done during its high flux states. The HEGRA observatory observed that
during the high flux states (during the strong outburst in 1997) Markarian 501’s spectral
energy distribution has an exponential cut off at ≈ 6 − 8 TeV [Aharonian et al.(1999)] and
[Aharonian et al.(2001)]. However, only a few gamma-ray observations have been performed
when Markarian 501 is in the low flux state (also known as the quiescent-state). In 2001, a
group of physicists performed a combined multiwavelength observation campaign to measure
the Markarian 501 flux at its quiescent state[Acciari et al.(2011)]. In this campaign both
VERITAS and MAGIC were utilized to measure the TeV flux in the energy range from ∼0.2
to ∼ 3 TeV. The TeV gamma-ray spectrum measured by both instruments was compatible
with a simple power law spectrum.

dN
= I0 × 10−12 ×
dE

−α
E
1 TeV

(8.3)

For the VERITAS measurements, the best-fit parameters are α = 2.72 ± 0.15stat ± 0.1sys
and I0 = 5.78±0.83stat ±1.16sys , and for the MAGIC measurement, the best-fit parameters
86

are α = 2.67 ± 0.21stat ± 0.2sys and I0 = 8.34 ± 1.53stat ± 2.50sys . The extrapolation of
−14 and 4.6+3.5 × 10−14 photons
these spectra gives the differential flux of 2.9+1.5
−1.0 × 10
−2.1
cm−2 s−1 TeV−1 at 7 TeV. Both of these extrapolated fluxes without a cutoff are above the
Milagro upper limit within their error bars. Therefore, I can conclude that the flux derived
at 7 TeV by extrapolating a pure power law SED of Markarian 501 is an over-estimate.
This suggests a power law spectrum with an exponential cut between 3 and 7 TeV for the
Markarian 501 spectrum.

8.3

Targeted Search for Potential TeV Emitting BL
Lac Objects

The observations of blazars by the Energetic Gamma-Ray Experiment Telescope (EGRET)
on NASA’s Compton Gamma-Ray Observatory satellite led to the discovery that blazars
emit most of their power in the gamma-ray band [Costamante & Ghisellini(2002)]. The
EGRET measured spectral energy distributions (SEDs) of the blazars were best characterized by two broad peaks, which are explained as synchrotron and Inverse Compton radiation. L. Costamante and G. Ghisellini studied these two-peak SEDs from the EGRET
and three other blazar samples available at that time (before year 2002), and developed a
model to predict a list of TeV emitting blazars from a subclass of blazars known as BL
Lacs [Costamante & Ghisellini(2002)]. Of this candidate list, 27 fall within the Milagro sky
coverage −7◦ < DEC < 80◦ . One method to constrain this model is to search for TeV
emission from these 27 BL Lac candidates and compare with the predictions. Therefore,
Elizabeth A. Hays performed a targeted search for this list of BL Lac candidates in the 2.7
year long Milagro archived data set [Hays(2004)]. In her search, she did not find significant
87

TeV emission from any of these TeV BL Lac candidates.
Since then the Milagro archived data set has grown by about a factor of three, and Milagro’s background rejection capability has substantially improved [Abdo(2007)]. Therefore,
I re-performed a targeted search for this list of TeV BL Lac candidates in the eight year
long Milagro archived data set. In addition to the larger data set and the new background
rejection capability, the Milagro map I used in my analysis was optimized for a spectral distribution of E−2.0 exp (−E/5 TeV), which increases the sensitivity to extragalactic sources.
Even though the new search was performed in a more sensitive Milagro data set, it was

still not sensitive enough to detect these TeV BL Lac candidates [Abeysekara & Milagro Collaboration(2013
Stacking of these candidates shows only a 0.28σ deviation from the background, which does
not indicate a significant collective gamma-ray emission. A graphical comparison between
the theoretically predicted TeV flux at 7 TeV, Milagro’s upper limit of each sky location
at 7 TeV, and the 95% expected upper limit of the flux for extragalactic sources in each
declination band corresponding to zero excess derived at 7 TeV is shown in Figure 8.2.
These results show that Milagro doesn’t observe collective or individual gamma-ray emission from these candidates. Therefore, the Milagro measurements were not able to constrain
L. Costamante and G. Ghisellini’s proposed model. However, these Milagro measurements
do provide upper limits to estimate the TeV cutoffs for new models. Since Figure 8.2 provides the expected upper limit at 95% confidence level for all the declination bands between
−7◦ and 80◦ , a theoretician can use these limits to estimate the TeV cutoffs of any source
in Milagro sky coverage.

88

Fermi Name
2FGL

RA
(deg)

DEC
(deg)

l
(deg)

b
(deg)

J0007.8+4713
J0009.1+5030
J0022.5+0607
J0045.3+2127
J0100.2+0746

1.97
2.29
5.64
11.34
15.06

47.23
50.51
6.12
21.45
7.78

115.3
116.09
110.02
121.04
126.74

-15
-11.8
-56.02
-41.4
-55.03

Flux/Flux Limit
(×10−17 TeV−1
s−1 cm−2 )
< 65.06
< 85.8
< 279.67
< 120.71
< 378.43

Source Type

Significance
σ

Associated source

bzb
agu
bzb
bzb
bzb

-0.78
-0.26
-0.1
0.48
1.69

MG4 J000800+4712
NVSS J000922+503028
PKS 0019+058
GB6 J0045+2127
GB6 J0100+0745

J0106.5+4854
J0108.6+0135
J0112.1+2245
J0112.8+3208
J0115.4+0358

16.65
17.17
18.03
18.21
18.87

48.91
1.59
22.76
32.14
3.97

125.49
131.85
129.15
128.19
134.43

-13.88
-60.98
-39.86
-30.51
-58.37

< 95.99
< 566.7
< 63.03
< 48.2
< 296.62

bzq
bzb
bzq
bzb

0.21
0.64
-1.29
-1.71
-0.59

4C 1.02
S2 0109+22
4C 31.03
PMN J0115+0356

J0136.5+3905
J0136.9+4751
J0144.6+2704
J0153.9+0823
J0211.2+1050

24.14
24.24
26.16
28.49
32.81

39.09
47.86
27.08
8.4
10.84

132.42
130.78
137.29
148.21
152.59

-22.95
-14.32
-34.31
-51.38
-47.39

< 70.3
< 117.84
< 85.04
< 209.19
< 106.04

bzb
bzq
bzb
bzb
bzb

-0.46
0.89
-0.03
-0.2
-1.48

B3 0133+388
OC 457
TXS 0141+268
GB6 J0154+0823
MG1 J021114+1051

J0217.4+0836
J0217.9+0143
J0221.0+3555
J0222.6+4302
J0237.8+2846

34.35
34.48
35.27
35.66
39.47

8.61
1.73
35.93
43.04
28.78

156.17
162.2
142.6
140.14
149.48

-48.63
-54.41
-23.49
-16.77
-28.55

< 138.64
< 347.2
< 85.82
< 85.18
< 101.74

bzb
bzq
bzq
BZB
bzq

-1.14
-0.84
0.21
0.13
0.58

ZS 0214+083
PKS 0215+015
S4 0218+35
3C 66A
4C 28.07

J0238.7+1637
J0316.1+0904
J0319.8+4130
J0326.1+0224
J0333.7+2918

39.68
49.05
49.97
51.55
53.43

16.62
9.08
41.51
2.41
29.31

156.78
172.1
150.58
180.74
160.49

-39.1
-39.59
-13.25
-42.45
-21.49

< 131.21
< 339.0
< 74.54
< 802.03
< 103.39

BZB
bzb
rdg
bzb
agu

0.24
1.59
-0.15
1.66
0.62

AO 0235+164
GB6 J0316+0904
NGC 1275
1H 0323+022
TXS 0330+291

J0423.2-0120
J0433.5+2905
J0442.7-0017
J0448.9+1121
J0508.0+6737

65.81
68.39
70.69
72.24
77.01

-1.34
29.09
0.29
11.36
67.63

195.28
170.52
197.21
187.4
143.8

-33.15
-12.62
-28.44
-20.77
15.9

< 1039.11
< 114.07
< 863.75
< 198.53
< 648.35

BZQ
bzb
bzq
bzq
bzb

1.35
0.93
1.04
0.41
2.53

PKS 0420-01
MG2 J043337+2905
PKS 0440-00
PKS 0446+11
1ES 0502+675

J0509.4+0542
J0532.7+0733
J0534.8-0548c
J0541.8-0203c
J0547.1+0020c

77.37
83.19
83.72
85.45
86.8

5.7
7.56
-5.81
-2.06
0.34

195.4
196.84
209.36
206.69
205.15

-19.62
-13.71
-19.66
-16.41
-14.1

< 320.94
< 275.2
< 942.7
< 1007.43
< 666.58

bzb
bzq

0.34
0.68
-1.51
0.66
0.24

TXS 0506+056
OG 50

J0607.4+4739
J0612.8+4122
J0650.7+2505
J0654.2+4514
J0654.5+5043

91.87
93.21
102.7
103.57
103.65

47.66
41.37
25.1
45.24
50.72

165.64
171.83
190.24
171.2
165.68

12.87
10.92
11.02
19.36
21.14

< 76.14
< 86.2
< 45.99
< 77.21
< 164.55

bzb
bzb
bzb
bzq
bzq

-0.4
0.24
-2.17
-0.14
1.79

TXS 0603+476
B3 0609+413
1ES 0647+250
B3 0650+453
GB6 J0654+5042

J0714.0+1933
J0719.3+3306
J0721.9+7120
J0725.3+1426
J0738.0+1742

108.51
109.83
110.48
111.33
114.52

19.57
33.11
71.35
14.44
17.7

197.68
185.06
143.97
203.63
201.85

13.61
19.85
28.02
13.93
18.06

< 103.44
< 115.37
< 512.3
< 169.64
< 152.56

bzq
bzq
bzb
BZQ
bzb

-0.03
1.02
0.49
0.49
0.83

MG2 J071354+1934
B2 0716+33
S5 0716+71
4C 14.23
PKS 0735+17

Table 8.2: Summary of the search with λ = 0.01 for TeV emission from the 2FGL list that
are identified as candidates off the galactic plane. (Note that we used the same abbreviations
for the source type as the 2FGL, agu = active galaxy of uncertain type, bzb = BL Lac type
of blazar and bzq = FSRQ type of blazar.) The Milagro flux derived at 7 TeV is given for
candidates that passed the λ = 0.01 FDR cut and the 95% confidence level flux upper limit
is given for the rest.

89

Table 8.2: (cont’d)
Fermi Name
2FGL

RA
(deg)

DEC
(deg)

l
(deg)

b
(deg)

Source Type

Significance
σ

Associated source

11.39
30.79
13.99
32.91
33.41

Flux/Flux Limit
(×10−17 TeV−1
s−1 cm−2 )
< 293.61
< 1032.12
< 1883.54
< 136.9
< 113.39

J0739.2+0138
J0805.3+7535
J0807.1-0543
J0809.8+5218
J0818.2+4223

114.82
121.34
121.78
122.46
124.57

1.65
75.59
-5.72
52.31
42.4

216.96
138.88
227
166.26
178.21

bzq
bzb
bzb
bzb
bzb

-1.35
-0.11
0.61
1.07
1.15

PKS 0736+01
RX J0805.4+7534
PKS 0804-05
0806+524
S4 0814+42

J0831.9+0429
J0854.8+2005
J0905.6+1357
J0909.1+0121
J0909.7-0229

127.99
133.71
136.4
137.29
137.43

4.49
20.1
13.96
1.37
-2.5

220.72
206.83
215.04
228.93
232.8

24.36
35.83
35.96
30.92
28.99

< 417.64
< 93.21
< 190.04
< 576.86
< 476.09

bzb
BZB
bzb
bzq
bzq

0.4
-0.35
0.88
0.19
-1.54

PKS 0829+046
OJ 287
MG1 J090534+1358
PKS 0906+01
PKS 0907-023

J0915.8+2932
J0920.9+4441
J0957.7+5522
J1012.6+2440
J1015.1+4925

138.96
140.24
149.43
153.17
153.79

29.54
44.7
55.38
24.68
49.43

196.67
175.7
158.59
207.74
165.53

42.93
44.81
47.94
54.36
52.73

< 134.29
< 107.06
< 86.48
< 77.99
< 124.22

bzb
bzq
bzq
bzq
bzb

1.58
0.84
-0.63
-0.37
1.04

B2 0912+29
S4 0917+44
4C 55.17
MG2 J101241+2439
1H 1013+498

J1016.0+0513
J1033.9+6050
J1037.6+5712
J1058.4+0133
J1058.6+5628

154.01
158.48
159.42
164.61
164.67

5.23
60.84
57.21
1.57
56.48

236.51
147.8
151.77
251.5
149.57

47.04
49.13
51.77
52.77
54.42

< 412.37
< 136.0
< 89.82
< 451.8
< 110.53

bzq
BZQ
bzb
bzb
bzb

0.74
-0.18
-0.82
0.01
0.02

TXS 1013+054
S4 1030+61
GB6 J1037+5711
4C 1.28
TXS 1055+567

J1104.4+3812
J1117.2+2013
J1121.5-0554
J1132.9+0033
J1150.5+4154

166.12
169.31
170.39
173.23
177.63

38.21
20.23
-5.91
0.56
41.91

179.82
225.63
266.27
264.33
159.14

65.03
67.39
50.45
57.42
70.67

389.74±40.7
< 90.1
< 2362.34
< 798.62
< 80.64

bzb
bzb
bzq
bzb
bzb

9.57
-0.45
1.39
1.3
0.12

Mkn 421
RBS 958
PKS 1118-05
PKS B1130+008
RBS 1040

J1159.5+2914
J1217.8+3006
J1221.3+3010
J1221.4+2814
J1224.9+2122

179.88
184.47
185.35
185.37
186.23

29.25
30.11
30.18
28.24
21.38

199.41
188.93
186.33
201.69
255.07

78.37
82.06
82.74
83.28
81.66

< 79.56
< 74.31
< 77.42
< 113.67
< 168.32

bzq
bzb
bzb
bzb
BZQ

-0.11
-0.25
-0.12
0.9
1.59

Ton 599
1ES 1215+303
PG 1218+304
W Comae
4C 21.35

J1226.0+2953
J1229.1+0202
J1231.7+2848
J1239.5+0443
J1243.1+3627

186.52
187.28
187.94
189.88
190.78

29.9
2.04
28.81
4.73
36.45

185.02
289.95
190.66
295.18
133.13

83.78
64.35
85.34
67.42
80.51

< 60.44
< 459.81
< 54.27
< 405.6
< 67.63

BZQ
bzb
bzq
bzb

-0.87
0.01
-1.28
0.73
-0.45

3C 273
B2 1229+29
MG1 J123931+0443
Ton 116

J1248.2+5820
J1253.1+5302
J1256.1-0547
J1303.1+2435
J1309.4+4304

192.06
193.28
194.04
195.78
197.37

58.35
53.05
-5.79
24.6
43.08

123.74
122.36
305.1
349.62
111.17

58.77
64.08
57.06
86.35
73.64

< 107.89
< 72.69
< 1072.78
< 117.62
< 87.93

bzb
bzb
BZQ
bzb
bzb

-0.44
-0.99
-1.07
0.85
0.26

PG 1246+586
S4 1250+53
3C 279
MG2 J130304+2434
B3 1307+433

J1310.6+3222
J1312.8+4828
J1418.4-0234
J1427.0+2347
J1438.7+3712

197.67
198.21
214.6
216.76
219.68

32.38
48.47
-2.57
23.8
37.21

85.59
113.32
341.56
29.48
63.72

83.29
68.25
53.64
68.2
65.27

< 65.29
< 71.38
< 1091.11
< 96.03
< 50.77

bzq
bzq
bzb
bzb
bzq

-0.68
-0.57
1.01
0.13
-1.33

OP 313
GB 1310+487
BZB J1418-0233
PKS 1424+240
B2 1436+37B

90

Table 8.2: (cont’d)
Fermi Name
2FGL

RA
(deg)

DEC
(deg)

l
(deg)

b
(deg)

Source Type

Significance
σ

Associated source

56.46
60.34
54.58
42.48
57.02

Flux/Flux Limit
(×10−17 TeV−1
s−1 cm−2 )
< 109.12
< 114.2
< 185.18
< 820.68
< 101.94

J1454.4+5123
J1501.0+2238
J1504.3+1029
J1520.8-0349
J1522.1+3144

223.62
225.28
226.1
230.22
230.54

51.4
22.64
10.49
-3.83
31.74

87.66
31.46
11.37
358.11
50.18

bzb
bzb
BZQ
bzb
bzq

0.49
0.55
0.0
-0.48
0.64

TXS 1452+516
MS 1458.8+2249
PKS 1502+106
NVSS J152048-034850
B2 1520+31

J1542.9+6129
J1553.5+1255
J1555.7+1111
J1607.0+1552
J1625.2-0020

235.73
238.39
238.94
241.77
246.3

61.49
12.93
11.19
15.88
0.33

95.38
23.77
21.92
29.4
13.92

45.4
45.21
43.95
43.42
31.83

< 93.21
< 109.28
< 160.64
< 85.05
< 762.2

bzb
bzq
bzb
bzb

-1.34
-0.97
-0.17
-1.16
0.66

GB6 J1542+6129
PKS 1551+130
PG 1553+113
4C 15.54

J1635.2+3810
J1637.7+4714
J1640.7+3945
J1642.9+3949
J1640.7+3945

248.81
249.43
250.18
250.18
250.75

38.17
47.24
39.76
39.76
39.83

61.13
73.38
63.35
63.48
63.35

42.34
41.88
41.38
40.95
41.38

< 58.66
< 68.76
< 122.41
< 122.41
< 115.23

bzq
bzq
BZQ
BZQ
BZQ

-0.95
-0.6
1.31
1.31
1.11

4C 38.41
4C 47.44
NRAO 512
3C 345
NRAO 512

J1642.9+3949
J1653.6-0159
J1653.9+3945
J1700.2+6831
J1709.7+4319

250.75
253.4
253.48
255.06
257.45

39.83
-2
39.76
68.52
43.32

63.48
16.59
63.61
99.58
68.41

40.95
24.93
38.85
35.19
36.21

<
<
<
<
<

115.23
847.78
186.65
316.95
113.14

BZQ

3C 345

BZB
bzq
bzq

1.11
0.15
2.93
-0.1
1.01

Mkn 501
TXS 1700+685
B3 1708+433

J1719.3+1744
J1722.7+1013
J1725.0+1151
J1734.3+3858
J1748.8+7006

259.83
260.68
261.27
263.58
267.22

17.74
10.23
11.87
38.98
70.11

39.53
32.22
34.11
64.04
100.54

28.07
24.3
24.47
31.02
30.69

< 183.16
< 425.7
< 288.43
< 95.96
< 275.01

bzb
bzq
bzb
bzq
bzb

1.41
2.84
1.82
0.48
-1.01

PKS 1717+177
TXS 1720+102
1H 1720+117
B2 1732+38A
1749+70

J1751.5+0938
J1754.3+3212
J1800.5+7829
J1806.7+6948
J1811.3+0339

267.88
268.58
270.15
271.68
272.83

9.64
32.2
78.48
69.8
3.66

34.91
57.75
110.06
100.1
31.62

17.65
25.38
29.07
29.18
10.59

< 183.43
< 60.01
< 1134.67
< 497.66
< 279.71

bzb
bzb
bzb
bzb
bzb

0.0
-0.97
-0.91
0.75
-0.73

OT 81
RX J1754.1+3212
S5 1803+784
3C 371
NVSS J181118+034114

J1824.0+5650
J1838.7+4759
J1849.4+6706
J1852.5+4856
J1903.3+5539

276
279.7
282.35
283.13
285.84

56.84
47.99
67.1
48.94
55.67

85.72
76.9
97.5
78.6
85.96

26.09
21.82
25.03
19.94
20.51

< 65.02
< 107.68
< 287.58
< 41.8
< 77.36

bzb
bzb
bzq
bzq
bzb

-1.93
0.63
0.2
-2.49
-1.1

4C 56.27
GB6 J1838+4802
S4 1849+67
S4 1851+48
TXS 1902+556

J1927.0+6153
J2000.0+6509
J2116.2+3339
J2121.0+1901
J2133.9+6645

291.77
300.02
319.05
320.26
323.49

61.9
65.16
33.66
19.03
66.75

93.31
98.02
79.82
69.25
105.17

19.71
17.67
-10.64
-21.25
10.96

<
<
<
<
<

123.93
208.54
124.07
168.15
316.94

bzb
bzb
bzb
bzq

-0.79
-0.07
1.32
1.26
0.47

1RXS J192649.5+615445
1ES 1959+650
B2 2114+33
OX 131

J2143.5+1743
J2147.3+0930
J2202.8+4216
J2203.4+1726
J2236.4+2828

325.88
326.84
330.71
330.87
339.1

17.72
9.51
42.27
17.44
28.48

72.09
65.85
92.6
75.68
90.12

-26.08
-32.28
-10.46
-29.63
-25.66

< 108.01
< 162.46
< 47.82
< 164.9
< 86.69

bzq
bzq
bzb
bzq
bzb

-0.21
-0.35
-1.53
0.9
0.05

OX 169
PKS 2144+092
BL Lacertae
PKS 2201+171
B2 2234+28A

91

Table 8.2: (cont’d)
Fermi Name
2FGL

RA
(deg)

DEC
(deg)

l
(deg)

b
(deg)

Source Type

Significance
σ

Associated source

-33.37
-15.77
-38.18
-24.02
-58.23

Flux/Flux Limit
(×10−17 TeV−1
s−1 cm−2 )
< 84.11
< 100.26
< 135.72
< 62.49
< 552.24

J2243.9+2021
J2244.1+4059
J2253.9+1609
J2311.0+3425
J2323.6-0316

341
341.03
343.5
347.77
350.91

20.36
40.99
16.15
34.43
-3.28

86.59
98.5
86.12
100.42
77.78

bzb
bzb
BZQ
bzq
bzq

-0.7
0.69
0.23
-0.78
-1.9

RGB J2243+203
TXS 2241+406
3C 454.3
B2 2308+34
PKS 2320-035

J2323.8+4212
J2325.3+3957
J2334.8+1431
J2339.6-0532

350.95
351.33
353.72
354.91

42.2
39.96
14.53
-5.54

106.06
105.52
96.56
81.36

-17.78
-19.98
-44.39
-62.47

< 81.89
< 47.82
< 167.75
< 2061.09

bzb
bzb
bzb

0.13
-1.62
0.7
0.84

1ES 2321+419
B3 2322+396
BZB J2334+1408

J1104.4+3812
14
40

12
10
8
6

38

4
2
0
-2

36

-4
11h10m

11h00m

Figure 8.1: This map shows the 5◦ × 5◦ region around Markarian 421. The Fermi-LAT
source position is marked by a white dot. This map is made with the spectral optimization
− E
dN/dE ∝ E −2.0 e 5T eV and the data have been smoothed by a Gaussian point spread
function. The color of a bin shows the statistical significance (in standard deviations) of
that bin. The horizontal axis is right ascension in hours and the vertical axis is declination
in degrees.

92

Figure 8.2: The horizontal axis is the declination in the Milagro sky coverage. The vertical
axis is the energy flux in erg cm−2 s−1 . The red triangle data points are Milagro upper limits
of BL Lac objects at 7 TeV, the blue square data paints are the theoretically predicted flux
at 7 TeV and the black solid line is the 95% confidence level expected upper limit of the flux
for extragalactic sources corresponding to zero excess derived at 7 TeV for each declination
band of the Milagro sky maps made with spectral assumption E −2.0 exp −E/5TeV.

93

Chapter 9
The Milagro Targeted Search for GeV
Pulsars.
In 2009, the Milagro collaboration performed a targeted search [Abdo et al.(2009a)] of the
galactic sources in the Fermi Bright Source List (also known as Fermi-BSL) [Abdo et al.(2009b)],
and 14 of these sources were observed by Milagro above 3σ. Among these 14 sources, 9 have
spatial associations with pulsars. However, in the Fermi-BSL these 9 sources are not confirmed as GeV pulsars, and the GeV pulsed emissions were not isolated from the unpulsed
GeV emission. Therefore, these 9 sources can not be used to study the correlation between the pulsar GeV emission and the associated TeV emission. This issue was resolved
when the Fermi-LAT collaboration published their second catalog of gamma-ray pulsars
[Abdo et al.(2013)] with 117 high-confidence

0.1 GeV gamma-ray pulsar detections, using

a three year long data set acquired by the Fermi-LAT. Of these 117 GeV pulsars 32 are
located in the Milagro sky coverage.
I performed a targeted search for the Fermi-LAT GeV pulsars in the Milagro sky coverage
94

using a Milagro sky map made with the spectral assumption of a pure power law with spectral index of -2.6 and no cutoff. The detection threshold for this search was obtained using
the False Discovery Rate (FDR) technique, discussed in Section 7.1, with a 1% contamination fraction. In this list of 32 GeV pulsars, 14 passed the FDR cut. Milagro sky maps of a
5◦ × 5◦ region around the Fermi-LAT position of these 14 pulsars are shown in Figure 9.1.
The differential photon flux of the sources that were identified with TeV associations were
calculated at 35 TeV to minimize the dependency of the derived flux on the spectral assumption (refer to Section 4.4). The 95% confidence upper limit was derived for the Fermi GeV
pulsars that did not pass the FDR cut. The results of this search is summarized in Table 9.1,
and some other useful parameters of the Fermi-LAT pulsars are summarized in Table 9.2.
The authors of [Abdo et al.(2009a)] mention that the Milagro measured TeV emission associated with pulsars might come from the pulsar and/or from the associated PWN. However,
it is very unlikely to get a significant contribution from the pulsar to the TeV flux measured
by Milagro. The best example is the Crab pulsar, which is the brightest TeV object measured
by Milagro. The VERITAS Collaboration [VERITAS Collaboration et al.(2011)] observed
pulsed gamma-rays from the Crab pulsar in the energy range of ∼ 100 GeV to ∼ 200 GeV.
The measured energy spectrum for pulsed gamma-rays is well described by a simple power
law, without a cut-off:

dN
−11
= 4.2 ± 0.6stat ±2.4
1.4 syst ×10
dE

−3.8±0.5stat ±0.2syst
E
TeV−1 cm−2 s−1 .
150 GeV
(9.1)

An extrapolation of this energy spectrum gives a differential photon flux of 4.2 × 10−20
photons TeV−1 cm−2 s−1 at 35 TeV. This is 0.003% of the TeV flux observed by Milagro coincident with the Crab pulsar. In addition, the theoretical model proposed in
95

[Aharonian et al.(2012)] predicts a sharp cut-off below ∼500 GeV, so the extrapolated flux
at 35 TeV from the pulsar might be even lower. These considerations led me to conclude
that the TeV emissions observed coincident with pulsars come predominantly from their
associated PWNe.

96

Pulsar Name
PSR

R.A.
(deg)

DEC
(deg)

Significance
(σ’s)

FP W N
−14
×10
s−1 cm−2

J0007+7303
J0106+4855
J0205+6449
J0248+6021
J0357+3205

1.75
16.5
31.25
42
59.25

73.05
48.92
64.82
60.35
32.08

2.66
0.09
-0.95
1.4
-0.23

2.22 ± 0.83
< 0.81
< 0.67
< 1.55
< 0.7

J0534+2200
J0622+3749
J0631+1036
J0633+0632
J0633+1746

83.5
95.5
97.75
98.25
98.25

22
37.82
10.6
6.53
17.77

17.59
1.57
3.45
1.25
3.43

6.45 ± 0.36
< 1.28
1.67 ± 0.48
< 2.07
1.38 ± 0.4

J0659+1414
J1838-0537
J1836+5925
J1846+0919
J1907+0602

104.75 14.23
279.5 -5.62
279
59.42
281.5 9.32
286.75 6.03

0.84
1.27
-0.92
-1.06
6.86

<
<
<
<
4.22

1.43
6.56
0.54
0.64
± 0.61

J1952+3252
J1954+2836
J1957+5033
J1958+2846
J2021+3651

298.5 28.6
298
32.87
299.25 50.55
299.5 28.77
305.25 36.85

4.41
-0.61
1.86
3.98
12.38

1.45
<
<
1.31
4.18

± 0.32
0.57
1.35
± 0.32
± 0.33

J2021+4026
J2030+4415
J2028+3332
J2030+3641
J2032+4127

305.25 40.43
307.5 44.25
307
33.53
307.5 36.68
308
41.45

4.15
1.51
-0.29
4.43
7.81

1.37
<
<
1.48
2.53

± 0.33
1.21
0.68
± 0.33
± 0.32

Table 9.1: Summary of the search with λ = 0.01 for TeV emission from the pulsars in the
second Fermi-LAT pulsar catalog. The Milagro photon flux derived in the energy band from
34.5 TeV to 35.5 TeV is given for the candidates that passed the λ = 0.01 FDR cut and 95%
confidence level flux upper limits are given for the rest.

97

Table 9.1: (cont’d)

9.1

Pulsar Name
PSR

R.A.
(deg)

DEC
(deg)

Significance
(σ’s)

FP W N
−14
×10
s−1 cm−2

J2043+2740
J2055+2539
J2111+4606
J2139+4716
J2229+6114

310.75
313.75
317.75
324.75
337.25

27.67
25.65
46.1
47.27
61.23

-0.88
-0.12
1.11
-0.61
6.59

<
<
<
<
2.76

J2238+5903
J2240+5832

339.5
340

59.05
58.53

4.77
3.22

1.85 ± 0.38
1.25 ± 0.38

0.47
0.76
1.13
0.56
± 0.41

Comparison With IACT Flux Measurements

From this list of GeV pulsars J0534+2200, the Crab pulsar, was found to be the brightest
TeV association. It is the brightest TeV PWN in the northern hemisphere and is one of the
best measured TeV point sources by gamma-ray observatories. According to TeVCat, the
Crab was seen by Whipple, HEGRA, CAT, H.E.S.S., MAGIC, Milagro, Telescope Array,
CANGAROO, TACTIC, VERITAS, ARGO-YBJ and FACT observatories. Among these
observations the Milagro observatory measured the Crab spectral energy distribution at
the highest TeV energies[Abdo et al.(2012)], while H.E.S.S. measured the spectral energy
distribution in the second highest energy band. The spectral energy distribution measured
by H.E.S.S.[Aharonian et al.(2006a)] between 400 GeV and 40 TeV was found to follow a
power law with an exponential cutoff,

dN
= (3.76 ± 0.07) × 10−11 × E −2.39±0.03 exp(−E/14.3 ± 2.1) TeV−1 s−1 cm−2 . (9.2)
dE
98

Pulsar Name
PSR

G1 00
×10−11

log10

E˙
1 erg s−1

p

p˙

Distance

(ms)

×10−15

kpc

erg cm−2 s−1
J0007+7303
J0106+4855
J0205+6449
J0248+6021
J0357+3205

40.1 ± 0.4
1.9 ± 0.2
5.4 ± 0.2
5.2 ± 0.4
6.4 ± 0.2

35.6
34.4
37.4
35.3
33.7

315.9
83.2
65.7
217.1
444.1

357
0.43
190
55
13.1

1.4±0.3
3+1.1
−0.7
1.95±0.04
2±0.2
< 8.2

J0534+2200
J0622+3749
J0631+1036
J0633+0632
J0633+1746

129.3 ± 0.8
1.4 ± 0.1
4.7 ± 0.3
9.4 ± 0.5
423.3 ± 1.2

38.6
34.4
35.2
35
34.5

33.6
333.2
287.8
297.4
237.1

420
25.4
104
79.6
11

2±0.5
< 8.3
1±0.2
< 0.87
0.2+0.2
−0.1

J0659+1414
J1838-0537
J1836+5925
J1846+0919
J1907+0602

2.5 ± 0.2
18.8 ± 0.9
60.6 ± 0.4
2.4 ± 0.2
25.4 ± 0.6

34.5
36.7
34
34.5
36.4

384.9
145.7
173.3
225.6
106.6

55
465
1.5
9.9
86.7

0.28±0.03
< 24.1
0.5±0.3
< 22
3.2±0.3

J1952+3252
J1954+2836
J1957+5033
J1958+2846
J2021+3651

10.3 ± 0.4
13.8 ± 0.3
2.6 ± 0.1
9.1 ± 0.4
49.4 ± 0.8

36
36.5
33.6
35.5
36.5

39.5
92.7
374.8
290.4
103.7

5.8
21.2
6.8
212
95.6

< 18.6
2±0.5
< 14.5
< 18.5
10+2
−4

J2021+4026
J2030+4415
J2028+3332
J2030+3641
J2032+4127

95.5 ± 0.9
5.8 ± 0.4
5.8 ± 0.4
3.1 ± 0.3
10.6 ± 0.6

35
34.3
34.5
34.5
35.4

265.3
227.1
176.7
200.1
143.2

54.2
6.5
4.9
6.5
20.4

1.5±0.4
< 15.7
< 17.2
3±1
3.7±0.6

Table 9.2: Properties of the Fermi-LAT pulsars in the Milagro sky coverage. Column 2
(G100 ) is the the phase-averaged integral energy flux in the 0.1 to 100 GeV energy band.
Column 3 is the Log10 of the pulsars spin-down luminosity. Column 4 is the pulse period.
Column 5 is the first time derivative of the pulse period. Column 6 is the distance to pulsar.
All of these values are taken from the Fermi-LAT second pulsar catalog.

99

Table 9.2: (cont’d)
Pulsar Name
PSR

G1 00
×10−11

log10

E˙
1 erg s−1

p

p˙

Distance

(ms)

×10−15

kpc

erg cm−2 s−1
J2043+2740
J2055+2539
J2111+4606
J2139+4716
J2229+6114

1 ± 0.1
5.4 ± 0.2
4.4 ± 0.3
2.3 ± 0.2
25.3 ± 0.4

34.7
33.6
36.1
33.4
37.3

96.1
319.6
157.8
282.8
51.6

1.2
4.1
143
1.8
77.9

1.8±0.3
< 15.3
< 14.8
< 14.1
0.8+0.15
−0.2

J2238+5903
J2240+5832

6.4 ± 0.3
1.1 ± 0.3

35.9
35.3

162.7
139.9

97
15.2

< 12.4
7.7±0.7

Using this energy spectrum the flux derived in the energy band from 34.5 TeV to 35.5 TeV is
−14 TeV−1 s−1 cm−2 . Within error bars the H.E.S.S. measurement is consistent
2.4+6
−2 ×10
with the TeV flux measured by Milagro of: 5.4 ± 0.3 × 10−14 TeV−1 s−1 cm−2 .
Among the other candidates that passed our standard FDR cut, IACTs have measured the TeV flux coincident with three of these sources: J0007+7303, J2032+4127 and
J2229+6114. The names given to their spatially associated IACT TeV sources are CTA1,
TeV J2032+4130 and Boomerang, respectively. The most precise spectral distributions of
CTA1 and Boomerang were measured by VERITAS [Aliu et al.(2013), Acciari et al.(2009)]
and for TeV J2032+4130 the best fit spectral distribution is measured by H.E.S.S. [Aharonian et al.(2005)]. Among the candidates that failed our standard FDR cut, IACTs have
measured the TeV flux coincident with two sources: J0633+0632 and J1838-0537. The names
given to their spatially associated IACT TeV sources are HESS J0632+057 and HESS J1837069. The parameters of the best fit spectral distributions of these two sources are measured
by H.E.S.S. [Aharonian et al.(2007), Aharonian et al.(2006c)]. The integral photon flux de100

rived in the energy band from 34.5 TeV to 35.5 TeV of these five PWNe using the spectral
parameters measured by IACTs are summarized in column 2 of Table 9.3, and the integral photon flux in the same energy band measured by Milagro is summarized in column
3. As shown in this table the Milagro flux measurements and the IACT flux measurements
of the candidates that passed the FDR cut are consistent within their error bars. For the
two candidates that failed the FDR cut the Milagro upper limits are above the IACT flux
measurements. Therefore, I can conclude that the Milagro measurements are in agreement
with the IACT measurements.

Fermi Name

TeVCat Name

Milagro Flux
−14
×10
TeV s−1 cm−2

CTA1

IACT Flux
−14
×10
TeV s−1 cm−2 a
b
1.4+5
−1.1

J0007+7303
J2032+4127

TeV J2032+4130

a
3.5+2.7
−1.6

2.1 ± 0.3

J2229+6114

Boomerang

a
1.3+15
−1.2

2.76 ± 0.41

J0633+0632

HESS J0632+057

b
0.4+0.7
−0.3

< 2.07

J1838-0537

HESS J1837-069

b
5.4+1.6
−1.2

< 24.1

2.2 ± 0.83

Table 9.3: The energy flux derived in the energy band from 34.5 TeV to 35.5 TeV is summarized. The first column is the name given in the Fermi-LAT second pulsar catalog. The
second column is the name given in the TeVCat catalog. The third column is the flux
measured by IACTs. The fourth column is the flux measured by Milagro. (a) Flux measurements were derived from VERITAS measurements. (b) Flux measurements were derived
from H.E.S.S. measurements.

101

9.2

Population Study of Associated TeV emissions

In this section, I discuss the correlations between pulsar properties and the TeV flux of the
associated PWN. Table 9.2 lists the phase-averaged integral energy flux in the 0.1 to 100
˙ pulse period (P),
GeV energy band (G100 ), Log10 of the pulsar’s spin-down luminosityE,
the first time derivative of the pulse period P˙
luminosity of a PWN LP W N T eV

and the distance (d) to each pulsar. The

is related to its flux FP W N by:

LP W N T eV = 4πd2 FP W N .

(9.3)

However, the luminosity of pulsars LP SR GeV depends on the beaming factor fΩ (refer
to Section 2.2) in addition to the flux and distance:

LP SR GeV = 4πd2 fΩ G100 .

(9.4)

It is not possible to derive the fΩ that corresponds to each pulsar using the pulsar’s properties
measured by Fermi-LAT. Therefore, it is not possible to calculate the pulsar luminosity in
the 0.1 to 100 GeV energy band. However, the Pulsar GeV luminosity normalized with
respect to the beaming factor can be calculated by:

LP SR GeV /fΩ = 4πd2 G100 .

(9.5)

First, I examined the correlation between Pulsar GeV luminosity normalized with respect to
the beaming factor (LP SR GeV /fΩ ) and TeV luminosity of the associated PWNe (LP W N T eV ).
Figure 9.2 shows the correlation between LP SR GeV /fΩ and LP W N T eV . The blue
102

squares are the luminosity of the PWNe that passed the λ = 0.01 FDR cut and the red
triangles are the upper limit of the luminosity of the other sources. The large error bars on
this plot are due to the distance uncertainties. The candidates that passed the λ = 0.01
FDR cut yielded a linear correlation coefficient of 0.73. Furthermore, the best fit of the
candidates that passed the λ = 0.01 FDR cut, assuming a linear function, yielded

Log10 LP W N T eV = (0.79 ± 0.05) × Log10 LP SR GeV + 3.5 ± 1.8,

(9.6)

with χ2 /N DF = 129.531/9. The best fit of the same candidates, assuming a constant (slope
of the best fit line is zero), yielded:

Log10 LP W N T eV = 31.3 ± 0.03,

(9.7)

with χ2 /N DF = 326.25/10. This means that the improvement of the best fit to a linear
function is equivalent to 13.9σ. The linear correlation coefficient and the comparison between
the two fits demands a common underlying cause for the two emission mechanisms.

Next, in Figure 9.3 and Figure 9.4, I examined the correlation between the observed
˙ 2 of each PWNe.
flux on earth in the 34.5 to 35.5 TeV energy band FP W N and the E/d
˙ 2 are: E˙ is a good estimator for the rate of the
Motivations to choose the parameter E/d
energy that is deposited in the PWN by the associated pulsar, making a PWN associated
with a high E˙ pulsar have more energetic leptons to create TeV photons, and E˙ was divided
by d2 to consider the effect that the flux density decreases quadratically with the distance to
a source. The authors of [Carrigan et al.(2008)] argue that the probability of detecting very
˙ 2 , and they showed that for the
high energy gamma-rays from a PWN should vary with E/d
103

radio loud pulsars the fraction of PWNe detected with H.E.S.S. is about 70% for pulsars with
˙ 2 > 1035 ergs−1 kpc−2 while for pulsars with E/d
˙ 2 < 1033 ergs−1 kpc−2 it is 5%.
E/d
Before my thesis work started, the Milagro collaboration also investigated the association be˙ 2 of the associated GeV pulsars [Abdo et al.(2009a)].
tween the detected PWNe and the E/d
They noted that “the four high-confidence Milagro detections associated with pulsars of
˙ 2 > 1035 .
known periodicity and distance” have E/d
˙ 2 and F
In Figure 9.3 I examined the distribution between E/d
P W N for the GeV pulsars.
˙ 2
The fraction of the PWNe that passed the λ = 0.01 FDR cut in one-decade bins of the E/d
is shown in Figure 9.4. The slope of the best fit line for these data points is 0.02 ± 0.1, which
is consistent with zero. Therefore, these results do not demonstrate a connection between
˙ 2 of the associated
the probability of detecting TeV gamma-rays from a PWN and the E/d
GeV pulsar, as H.E.S.S. observed for radio loud pulsars. (Note that H.E.S.S. studies radio
loud pulsars and I studied GeV loud pulsars. So their are different selection effects between
the two samples.).
Next, I examined the distribution of the PWNe luminosity in the 34.5 to 35.5 TeV energy
˙ and characteristic age
band (LP W N T eV ) as a function of the spin-down luminosity (E)
(τc ) of the associated pulsars. In 2009, Mattana et al. investigated the LP W N T eV vs.
E˙ correlation, and the LP W N T eV vs. τc correlation using the PWNe observed by the
H.E.S.S. observatory, and they discussed the theoretical expectations of these two distributions in the leptonic scenario [Mattana et al.(2009)], where the PWNe are energized by
the relativistic electrons injected into the PWNe by the associated pulsars. In Chapter 2
I summarized their argument, which leads one to theoretically expect that the PWN luminosity in TeV energies should not be a function of the E˙ or τc of the associated pulsar,
LP W N T eV ∝ E˙ 0 and LP W N T eV ∝ τc0 .
104

In Figure 9.5(a) and Figure 9.5(b), I examined the distribution of PWNe luminosity in the
34.5 to 35.5 TeV energy band LP W N T eV vs. the spin-down luminosity of the associated
pulsars, and the LP W N T eV vs. the characteristic age (τc ) of the associated pulsars. The
linear correlation coefficient between the Log10 E˙

and Log10 LP W N T eV , calculated

only using the candidates that passed the λ = 0.01 FDR cut was, 0.3. Furthermore, the best
fit of the candidates that passed the λ = 0.01 FDR cut, assuming a linear function, yielded

Log10 LP W N T eV = (0.06 ± 0.03) × Log10 E˙ + 29 ± 1,

(9.8)

with χ2 /N DF = 327/9. The best fit of the same candidates assuming a constant (slope of
the best fit is zero) yielded

Log10 LP W N T eV = 31.3 ± 0.03,

(9.9)

with χ2 /N DF = 326.25/10. This means the linear assumption is consistent with zero with a
2σ error bar and the improvement of the best fit to a linear function is only a 1.1σ equivalent.
Therefore, the LP W N T eV vs. E˙ distribution derived from Milagro observations is not
consistent with the theoretical expectations of [Mattana et al.(2009)] within 1σ error bar,
but they become consistent with 2σ error bars. Note that these error bars are estimated
only using the statistical errors. Therefore, these error bars become larger with inclusion of
the systematic errors. I re-examined this correlation in Chapter 10, adding 7 data points
measured by VERITAS and H.E.S.S. With those new data points the slope of the best fit
for the LP W N T eV vs. E˙ distribution became consistent with zero, within 1σ error bar.

The linear correlation coefficient between Log10 LP W N T eV
105

and Log10 (τc ), calcu-

lated using only the candidates that passed the λ = 0.01 FDR cut was -0.28. Furthermore,
the best fit of the candidates that passed the λ = 0.01 FDR cut, assuming a linear function,
yielded

Log10 LP W N T eV = (−0.036 ± 0.056) × Log10 E˙ + 31.42 ± 0.09,

(9.10)

with χ2 /N DF = 327/9. The slope of the best fit with a linear assumption is consistent with
zero. Therefore, I can conclude that the LP W N T eV vs. τc distribution derived from Milagro data set is consistent with the theoretical expectation derived by [Mattana et al.(2009)].

Lastly, I examined the TeV efficiency of the PWNe at 35 TeV (εγ ) as a function of E˙
and τc . The efficiency is defined as:

εγ =

LP W N T eV
,
E˙

(9.11)

where εγ is the fraction of the spin-down luminosity that becomes the TeV flux in the 34.5 to
35.5 TeV energy band. The distributions of Log10 εγ vs. Log10 E˙ and Log10 εγ vs.
Log10 (τc ) are shown in Figure 9.6(a) and Figure 9.6(b), respectively. The linear correlation
coefficient between Log10 E˙

and Log10 εγ , calculated using only the candidates that

passed the λ = 0.01 FDR cut was -0.6. However, this correlation was mostly driven by the
two pulsars with E˙ > 1037 erg s−1 . The linear correlation coefficient of this distribution
excluding those pulsars is 0.08. Therefore, these results were not favorable to concluding
that there is a dependency between εγ and the E˙ of the associated pulsars. However, the
εγ of the two PWNe that have associated pulsars with E˙ > 1037 erg s−1 are more than an
order of magnitude less than the εγ of those associated with low E˙ pulsars. The distribution
106

of Log10 (τc ) vs. Log10 εγ also had similar features. The linear correlation coefficient
between Log10 (τc ) and Log10 εγ , calculated using only the candidates that passed the
λ = 0.01 FDR cut was 0.6. However, this correlation was mostly driven by the two pulsars
that have τc < 101.1 kyr. The linear correlation coefficient of this distribution excluding
those pulsars become 0.08. Therefore, these results were not favorable to conclude that there
is a dependency between εγ and τc of the associated pulsars.

9.3

Discussion

The most widely discussed gamma-ray emission mechanisms are either hadronic or leptonic.
In the hadronic scenario, gamma-rays are produced by secondary π 0 decays, which are
produced by the collisions between the highly energetic hadrons and nuclei in the ambient
medium [Gabici et al.(2009)]. In the leptonic scenario, gamma-rays are produced by relativistic electrons and positrons. Most of the pulsar gamma-ray emission models predict that
the GeV photons in the 0.1 to 100 GeV energy range should be dominated by curvature radiation of the relativistic leptons emitted by the rotating neutron star [Muslimov & Harding(2004)].
Most of the PWN TeV gamma-ray emission models predict that TeV photons are produced
from inverse Compton(IC) up-scattering of ambient photon fields, such as cosmic microwave
background, stellar radiation, infrared emission from dust, or synchrotron radiation by the
ultra-relativistic electrons/positrons [de Jager & Djannati-Ata¨ı(2009), de Jager et al.(2009)].
This suggests that in the leptonic scenario both the pulsar GeV luminosity in the 0.1 to 100
GeV energy range LP SR GeV

and the PWN TeV luminosity LP W N T eV

are tightly

related to the properties of the electron-positron wind coming from the rotating neutron star.
Therefore, the studies of correlations between LP W N T eV and pulsar properties (such as
107

˙ τc and L
E,
P SR GeV ) are important to constrain the leptonic scenario.
The correlation between LP SR GeV and LP W N T eV shown in Figure 9.2 leads one
to suspect a common underlying cause for the two emission mechanisms. In the leptonic
scenario, one property relevant to both emissions is the electron-positron current emitted
by the rotating neutron star, N˙e . The GeV energy flux from pulsars is thought to be
directly related to the instantaneous value of N˙e , because the GeV pulsed emission from
the magnetosphere is often thought to be produced from curvature emission by the most
recently produced electron-positron population from the rotating neutron star.
Inside PWNe, TeV photons are often thought to be produced by IC. Therefore, LP W N T eV
should depend both on the relativistic electron-positron population and the ambient photon
population in the PWN. By considering the cooling time [Mattana et al.(2009)] suggested
that the electron-positron population that produces TeV photons by IC becomes proportional
to the integral of N˙e over the pulsar lifetime, instead of being proportional to the instantaneous value of N˙e . Hence, [Mattana et al.(2009)] argued that LP W N T eV should be
proportional to the integral of N˙e over the pulsar lifetime. However, [Mattana et al.(2009)]
did not account for the ambient photon field density ρph when deriving LP W N T eV . I
suggest that the proportionality between LP W N T eV and LP SR GeV might occur due to
a proportionality between the ambient photon field density ρph and N˙e .
There are two different ambient photon fields which could be relevant to the production
of TeV gamma-rays and directly related to N˙ e : photons from synchrotron radiation and farinfrared photons [Atoyan & Aharonian(1996)]. The density of synchrotron radiation photons
in the x-ray energy band is roughly proportional to the density of the freshly injected pulsar
wind [Mattana et al.(2009)]. In addition, far-infrared seed photons can be made by heating
the pulsar wind, as described in Section 2.2 of [Arons(1996)]. Therefore ρph may be roughly
108

proportional to N˙e , instead of the integral of N˙e . While these considerations are suggestive,
a more detailed theoretical study is clearly needed to fully understand this correlation, which
is beyond the scope of the work presented in this dissertation.
Using Figure 9.3 and Figure 9.4, I investigated the detectability of a TeV PWN as a
˙ The lack of correlation between the detectability and E/d
˙ 2 suggests that the
function of E.
rate of energy that a pulsar feeds to the associated PWN is not the only parameter relevant
for the PWN TeV flux. The PWN TeV flux could also result from some other mechanisms
that are correlated with the pulsar or PWN, such as the age of the pulsar and adiabatic
cooling. For example, if the age of a pulsar is shorter than the cooling time of electrons that
produce TeV photons, then the relic electrons from earlier epochs might also contribute to
the PWN TeV emission.
Using Figure 9.6(a) and Figure 9.6(b), I investigated the efficiency (εγ ) of PWN TeV
emission at 35 TeV as a function of E˙ and τc . Neither, the εγ vs. E˙ nor the εγ vs.τc distribution demonstrate strong evidence for correlations. However, the two PWNe associated
with the youngest and highest E˙ pulsars have an εγ more than an order of magnitude less
than the other PWNe. These results are consistent with the IC model discussed by Mattana
et al (2009). When a PWN gets older it accumulates relativistic electrons that produce TeV
photons, and the current injection rate of electrons becomes negligible. Therefore, relatively
older PWNe might have a larger population of relativistic electrons to produce TeV photons.
However, the population of the relativistic electrons in a PWN cannot continue to increase
forever, because of cooling effects and that the magnetic field strength decays with time.
Mayer et al.(2012) have shown that the IC peak is expected to have a negative correlation
with the pulsar’s characteristic age, because the strong magnetic field strength decays with
time, and the cooling effects. My data set was not able to constrain the Mayer et al.(2012)’s
109

argument, because of the lack of older pulsars.

110

J0007+7303

J0534+2200

J0631+1036

14

J0633+1746

14

14

12

12

12

10

10

8

8

8

6

6

6

4

4

14

20

24

74

12

12

10

10
8
6

18

22

4
72

10

2

2

0

0

0

-2

-2

-2

4

2

2
0
-2

16

20

-4
00h40m

00h20m

00h00m

-4

23h40m

05h40m

(a)

06h40m

J1907+0602

(c)

J1952+3252

6

J1958+2846
14

12

12

10

10

30

8

8

8

6

6

6

4

4

4

28

2

2

28

0

0

0

-2

-2

-2

-4

-4

-4

19h00m

20h00m

(e)

19h50m

20h00m

(f)

J2021+4026

12
10

38

8
6
4
2

36

0
-2
-4

19h50m

20h30m

J2031+4127

20h10m

J2229+6114

14

14

12

12

12

10

10

38

20h20m

(h)

14

10

14

(g)

J2030+3641

42

J2021+3651

14

30

06h30m

(d)

12

2

4

06h40m

14

10

-4

06h30m

(b)

8

19h10m

-4

05h30m

8

8

6

6

8
6

4

4

4

2

2

0

0

14
12
10
62

8

42

40

2

36

0

6
4
2
60

0

40

38

20h30m

20h20m

-2

-2

-2

-4

-4

-4

20h10m

20h40m

(i)

20h30m

20h20m

20h40m

20h30m

20h20m

(j)

(k)

J2238+5903

J2240+5832

-2
-4
22h45m

22h30m

22h15m

(l)

14

14

12

12
60

60

10

10

8

8

6

6

4

4
58

58

2

2

0

0

-2

-2

-4
22h50m

22h40m

22h30m

22h20m

(m)

-4
22h50m

22h40m

22h30m

(n)

Figure 9.1: These maps show the 5◦ × 5◦ region around the pulsars that passed the λ = 0.01
FDR cut. The LAT source positions are marked by white dots. These maps are made with
the spectral optimization dN/dE ∝ E −2.6 and the data have been smoothed by a Gaussian
point spread function. The color of a bin shows the statistical significance (in standard
deviations) of that bin. The horizontal axis is right ascension in hours and the vertical axis
is declination in degrees.

111

TeV s-1)
PWN TeV

32

10

Log (L

33

31

30

29

33.5

34

34.5

35

35.5

36
Log (L
10

36.5
37
/f Ω TeV s-1)

PSR GeV

Figure 9.2: PWN TeV luminosity vs. Pulsar GeV luminosity normalized with respect to
the beaming factor fΩ . Blue squares are the luminosity derived in the energy band from
34.5 TeV to 35.5 TeV of the candidates that passed the λ = 0.01 FDR cut. Red triangles
are the upper limit of the luminosity derived in the energy band from 34.5 TeV to 35.5
TeV of the candidates that failed the λ = 0.01 FDR cut. The linear correlation coefficient
calculated only using the candidates that passed the λ = 0.01 FDR cut is 0.73, and the
best fit line for these data points have a slope of 0.79 ± 0.05 and intercept of 3.5 ± 1.8 and
χ2 /N DF = 129.531/9.

112

TeV s-1)
TeV
10

Log (F

-13
-13.2
-13.4
-13.6
-13.8
-14
-14.2
-14.4
33

34

35

36

37

38
Log (E/d2 erg s-1 kpc-2)
10

Detected Fraction of PWNe

Figure 9.3: The PWN gamma-ray flux in the 34.5 to 35.5 TeV energy band is shown (FT eV ).
The blue squares are the flux of the candidates that passed the λ = 0.01 FDR cut. The red
triangles are the upper limits of the candidates that failed the λ = 0.01 FDR cut.

1
0.8
0.6
0.4
0.2
0
33

34

35

36

37

38

39
Log (E/d2)
10

Figure 9.4: Fraction of the PWNe that passed the λ = 0.01 FDR cut in one-decade bins of
˙ 2.
the E/d

113

TeV s-1)

TeV s-1)
Log (L

PWN TeV

32

33

32

10

10

Log (L

PWN TeV

33

31

31

30
30

29

28

29

34

35

36

37

28
0

38
39
Log (E erg s-1)

0.5

1

1.5

2

2.5

3

3.5
4
Log (τc kyr)
10

10

(a)

(b)

Figure 9.5: Distributions of PWNe luminosity in the 34.5 to 35.5 TeV energy band
LP W N T eV vs. the spin-down luminosity E˙ of the associated pulsars (a), and LP W N T eV
vs. the characteristic age (τc ) of the associated pulsars (b) are shown. Blue squares are the
LP W N T eV of the candidates that passed the λ = 0.01 FDR cut. Red triangles are the
upper limit of the LP W N T eV of the candidates that failed the FDR cut.

10

Log (εγ )

10

Log (εγ )

-0.5
-1

-1.5

-1

-2

-2
-2.5

-3
-3
-3.5

-4

-4
-4.5

-5

-5
-5.5

34

35

36

37

-6
0

38
39
Log10(E erg s-1)

0.5

1

1.5

2

2.5

3
3.5
Log (τc kyr)
10

(a)

(b)

Figure 9.6: TeV efficiency of PWNe, εγ , vs. E˙ (a), and εγ vs. τc (b) distributions are
shown. εγ is defined as shown in Equation 10.16. Blue squares are the εγ of the candidates
that passed the λ = 0.01 FDR cut. Red triangles are the εγ upper limit of the candidates
that failed the FDR cut.

114

Chapter 10
Experimental Constraints on
Gamma-Ray Pulsar Gap Models
In this chapter I present a new method to isolate the flux correlation factor fΩ of GeV
pulsars as a function of the pulsar spin down luminosity E˙ . I applied this method to a
sample of PWNe measured in TeV energies, and obtained the fΩ vs. E˙ distribution of the
pulsars. I then compared my measured fΩ vs. E˙ distribution with the theoretically expected
fΩ vs. E˙ distributions of 4 pulsar gap models.
Currently, there are several pulsar gap models being discussed in the pulsar community.
These models can be constrained either by studying individual pulsars in detail or alternatively using the collective properties of a sample of pulsars. The large sample of pulsars
discovered by the Fermi Large Area Telescope (Fermi-LAT) provides a good place to study
the collective properties of GeV pulsars.
The GeV luminosity as a function of other pulsar parameters is a fundamental quantity
which the models must predict. However, the potential utility of this measurement is limited
115

by two factors: inability to measure fΩ and imprecise distance measurements. The factor
fΩ provides the correction to extrapolate the observed phase-averaged flux from the earth
line-of-sight to the full sky flux for a given beam shape. It is an essential factor needed to
convert observed fluxes to luminosity measurements:

LP SR GeV = 4πd2 fΩ F,

(10.1)

where LP SR GeV is the luminosity, d is the distance to the pulsar from earth and F is the
phase averaged flux measured at earth. Since fΩ is a model dependent parameter, luminosity calculations are also model dependent. Therefore, one option to constrain GeV pulsar
emission models is to use the collective properties of the luminosity distribution, but the
uncertainty of distance measurements degrades the accuracy of the luminosity distribution.
This issue can be resolved by studying the ratio of the flux from pulsars (G100 ) to that
of their associated PWNe FP W N , which is a distance independent parameter. In the
˙ using this distance
next section I present a new technique to obtain fΩ as a function of E,
independent parameter: G100 /FP W N .

10.1

A New Method to Isolate fΩ

The ratio of a pulsar’s GeV flux (G100 ) to PWN flux FP W N

can be turned into a

luminosity ratio by multiplying both the numerator and denominator by 4πd2 , yeilding:
G100
4πd2 G100
=
FP W N
4πd2 FP W N
116

(10.2)

The PWN luminosity LP W N is related to flux by:

LP W N = 4πd2 FP W N .

However, pulsar’s GeV luminosity LP SR GeV

(10.3)

is also dependent on the flux correlation

factor fΩ :
LP SR GeV = 4πd2 fΩ G100 .

(10.4)

Substituting Equation 10.3 and 10.4 into Equation 10.2 gives:

L
/f
G100
= P SR GeV Ω .
FP W N
LP W N

(10.5)

˙
This result can be combined with the theoretical correlation between LP SR GeV and E.
1
Most of the models expect LP SR GeV ∝ E˙ 2 above some threshold of open field-line voltage
(eg. [Arons(1996)], [Abdo et al.(2013)] , [Arons, J. 2006], [Harding(1981)],
[Harding & Muslimov(2002)], [Pierbattista et al.(2012)], [Muslimov & Harding(2003)],
[Harding & Muslimov(2003)], [Romani & Watters(2010)], [Takata et al.(2010)],
[Thompson, D.J. (2000)] , and [Thompson(2001)]). Using this relation Equation 10.5 becomes:

1
E˙ 2 /fΩ
G100
∝
.
FP W N
LP W N

(10.6)

˙ by using L
This equation can be used to obtain fΩ as a function of E,
P W N measured in
an appropriate energy band. The choice of the energy band to measure LP W N must fulfill
two conditions: have a tight correlation with LP SR GeV /fΩ and be a known function of
˙
E.

117

Two good possibilities are the PWN luminosity in the TeV energy band LP W N T eV
and the X-ray energy band LP W N X . As I discussed in Chapter 9 Mattana et al. (2009)
I argued that the expected distribution of LP W N T eV is LP W N T eV ∝ E˙ 0 . They also
validated their argument with a sample of PWNe measured in the TeV energy range. By
substituting their theoretical expectation, LP W N T eV ∝ E˙ 0 , in Equation 10.5 I can rewrite
Equation 10.5 as:
1
G100
E˙ 2
∝
FP W N T eV
fΩ
1 F
E˙ 2 × P W N T eV
G100

fΩ ∝

(10.7)

.

(10.8)

Mattana et al. (2009) also argued that the expected distribution of LP W N X is, LP W N X ∝
˙ and they also validated this argument with X-ray measurements. Therefore, I can rewrite
E,
˙ as:
Equation 10.5 with their theoretical expectation, LP W N X ∝ E,
−1
G100
E˙ 2
∝
FP W N X
fΩ

fΩ ∝

−1 F
E˙ 2 × P W N X
G100

(10.9)

.

(10.10)

In the work presented in this thesis I derived the fΩ as a function of E˙ for a Fermi GeV
pulsar sample using the pulsar GeV flux and the TeV flux of their associated PWNe.

118

10.2

The Sample of Pulsars and Their Associated
PWNe

In Chapter 9 I presented the PWNe TeV flux measurements made with Milagro for 14
Fermi GeV pulsars. From a literature survey I found the PWN TeV fluxes for 7 additional
PWNe observed with IACTs. Among these 7 PWNe the SEDs of 5 PWNe are measured
by H.E.S.S.: Crab, Vela, K3 in Kookabura, MSH 15-52, and G 21.5-0.9. In order to be
consistent with the Milagro measurements, I integrated the energy flux within a 1 TeV
energy band around 35 TeV using the H.E.S.S. SEDs. The derived energy flux of these
sources are shown in Table: 10.1. The SEDs of the Crab, Vela and MSH 15-52 PWNe were
measured in the energy ranges of 440 GeV - 40 TeV, 550 GeV - 65 TeV and 250 GeV - 40
TeV, respectively. Therefore, their integrated energy flux around 35 TeV can be obtained
without any extrapolation. However, the SEDs of K3 in Kookabura and G 21.5-0.9 PWNe
were measured in the energy ranges of 200 GeV - 25 TeV and 150 GeV - 5 TeV, respectively.
Therefore, their integrated energy fluxes obtained around 35 TeV by extrapolation of their
SEDs might be an overestimate, if there is a cutoff below 35 TeV.
VERITAS has published SEDs of the Boomerang and CTA 1 PWNe. In both cases the
SEDs of these sources were measured in the energy range of 1-15 TeV. Therefore, as for K3
and G 21.5-0.9, the integrated energy flux obtained around 35 TeV by extrapolating the SEDs
might be an overestimate, if there is a cutoff before 35 TeV. For all H.E.S.S. and VERITAS
measured PWNe, errors of the integrated flux are estimated by a standard Gaussian Monte
Carlo propagation of the uncertainties of the SED fit, with the 16th percentile as the lower
1σ error bar and the 84th percentile as the upper.
Before applying these measurements to obtain fΩ , I examined whether these measure119

Pulsar Name

Association

log10 E˙

J0007+7303

CTA 1

35.65

40.1±0.4

1.4±0.3

J0534+2200

Crab

38.64

129.3±0.8

2.0±0.5

J0835-4510

Vela

36.84

906±2

0.29 ± 0.02

J1420-6048

37.01

17.0±1.4

5.6±0.9

J1509-5850

K3
in Kookabura
MSH 15-52

FP W N T eV
10−14
TeV s−1 cm−2
a
1.4+4.9
−1.1
b
2.3+6.0
−2.0
c
16.4+9
−7
d
5.4±1.7
1.3

35.71

12.7±0.7

2.6±0.5

e
6.2±0.9
0.8

J1833-1034

G21.5-0.9

37.53

5.9±0.5

J2229+6114

Boomrang

37.35

25.3±0.4

4.7±0.4
0.8±+0.15
−0.20

f
0.99±10.6
g
1.3+15
−1.2

PSR

Distance

erg s−1

G100
10−11
erg cm−2 s−1

(kpc)

Table 10.1: G100 is the integrated phase-averaged energy flux of the pulsar GeV emission
in the 0.1-100 GeV energy band. FP W N T eV is the integrated energy flux of the PWN
TeV emission in the 34.5-35.5 TeV energy band. Properties of a sample of GeV pulsars
cataloged in the Fermi-LAT Second Pulsar Catalog and the TeV flux of their associated
PWNe. (a). VERITAS Measurement. Energy flux derived by extrapolating the SED, reference [Aliu et al.(2013)] (b). H.E.S.S. Measurement, reference [Aharonian et al.(2006a)]
(c). H.E.S.S. Measurement, reference [Aharonian et al.(2006b)] (d). H.E.S.S. Measurement, Energy flux derived by extrapolating the SED, reference [Aharonian et al.(2006d)] (e).
H.E.S.S. Measurement, reference [Aharonian et al.(2005)] (f). H.E.S.S. Measurement, Energy flux derived by extrapolating the SED, reference [H. E. S. S. Collaboration et al.(2007)]
(g). VERITAS Measurement, Energy flux derived by extrapolating the SED, reference
[Aliu et al.(2013)]

120

TeV s-1)
PWN TeV
10

Log (L

35
34
33
32
31
30
29
28
27
32.5

33

33.5

34

34.5

35

35.5

36
36.5
37
37.5
Log (L
/f Ω erg s-1)
10

PSR GeV

Figure 10.1: PWN TeV luminosity vs. Pulsar GeV luminosity normalized with respect to the
flux correlation factor fΩ . Blue squares are measured by Milagro. Red stars are measured
by H.E.S.S. Green circles are measured by VERITAS. The dotted line is the best fit line.
The linear correlation coefficient is R = 0.8. The error bars are dominated by the distance
measurement uncertainties.

121

TeV s-1)
PWN TeV

34
33

10

Log (L

35

32
31
30
29
28
34.5

35

35.5

36

36.5

37

37.5

38

38.5
39
Log (E erg s-1)
10

Figure 10.2: PWN TeV luminosity vs. spin-down luminosity E˙ of the associated pulsar.
Blue squares are for PWNe measured by Milagro, red stars by H.E.S.S., and green circles
by VERITAS. The distance uncertainty contributes significantly to the error bars. The data
˙
are consistent with no dependence of LP W N T eV on E.

122

ments fulfill the two conditions required to apply the method: have a tight correlation with
˙ Using Figure 10.1 I examined the correlaLP SR GeV /fΩ and be a known function of E.
tion between LP SR GeV /fΩ and LP W N T eV . The blue squares are PWNe measured by
Milgro, the red stars by H.E.S.S., and the green circles by VERITAS. It appears that the
PWNe measurements of the three experiments are generally consistent. Four of the TeV
measurements (LP SR = 34.3, 35.0, 35.2, 35.8) have been extrapolated from lower energy
bands, but do not appear to be outliers to the general scatter of the distribution. The large
error bars on this plot are due to the distance uncertainties. Even though the error bars are
large a reasonable correlation is obtained with a linear correlation coefficient of R = 0.82,
and the best fit line for these data points yielded:

log10 (LP W N T eV ) = (0.84 ± 0.1) log10 (LP SR GeV /fΩ ) + (1.7 ± 3.8),

(10.11)

with χ2 /N DF = 35.653/19. This provides evidence to conclude that this data fulfills the
first condition listed above. Second, I examined the correlation between LP W N T eV and
E˙ , which was previously studied by Mattana et al. (2009) using H.E.S.S. measurements.
The correlation between LP W N T eV and E˙ for my PWN sample is shown in Figure 10.2.
The best fit of these data points assuming a linear function yielded

˙ + (30.6 ± 2.5),
log10 (LP W N T eV ) = (0.02 ± 0.06) × log10 (E)

(10.12)

which has a slope consistent with zero and also with the slope obtained from the PWN
sample in Mattana et al. (2009) of 0.07 ± 0.10. Those two observations led me to conclude
that my data sample fulfills both required conditions.

123

10 Ω

Log (f )

16
15.5
15
14.5
14
13.5
13
12.5
12
11.5
11

34.5

35

35.5

36

36.5

37

37.5

38

38.5
39
Log (E erg s-1)
10

˙ Blue squares
Figure 10.3: Experimentally obtained fΩ vs. pulsar spin-down luminosity E.
are PWNe measured by Milagro, red stars by H.E.S.S, and green circles by VERITAS. The
best fit line for the pulsars with E˙ > 1035.1 ergs s−1 is shown as a dotted line.

10.3

Results and Discussion

˙ using Equation 10.8,
The distribution of the experimentally obtained fΩ as a function of E,
is shown in Figure 10.3. Note that the fΩ values obtained by this technique are not the
˙ because the correct normalization between
actual values of fΩ but only its dependence on E,
the two energy bands is not predicted. It appears that above E˙ ≈ 1035 erg s−1 , the fΩ vs.
E˙ distribution has a linear correlation with a slope of 0.28 ± 0.03 (χ2 /N DF = 33/18), and
below E˙ ≈ 1035 ergs s−1 the slope is steeper and fΩ distribution for a given E˙ has a tail
towards smaller fΩ values.

This work allowed me to compare the experimentally obetained fΩ vs. E˙ distribution
with the theoretically expected fΩ vs. E˙ distributions made by Pierbattista et al.(2012)
124

under four pulsar gap models. Pierbattista et al.(2012) studied four gamma-ray pulsar
acceleration models. They synthesized a pulsar population based on a radio emission model
and four gamma-ray pulsar gap models: Outer Gap (OG), Polar Cap (PC), Slot Gap (SG)
and One Pole Caustic (OPC). Their model simulations of the correlation between fΩ and
E˙ are shown in Figure 10.4. In all four model predictions, for pulsars with E˙ > 1035 ergs
s−1 , fΩ can be reasonably fit by a straight line in log − log space and the best fits give
the following slopes: mP C = 0.05 ± 0.23, mSG = 0.003 ± 0.004, mOG = 0.12 ± 0.01 and
mOP C = 0.026 ± 0.007.1 The slopes of the PC and SG models are consistent with zero, and
the slope of the OPC is very small. The OG model is the only model that predicts a positive
slope that is not consistent with zero. However, the small number of data points for the PC
model have a large scatter, and, while the slope is consistent with zero, the uncertainty on
the slope is much larger than for the other models.
I can now compare my results from Figure 10.3 with the expectations from the four pulsar
gap models. The experimental fΩ vs. E˙ distribution has a non-zero slope of 0.28 ± 0.03 for
E˙ > 1035 erg s−1 . This would tend to disfavor the PC, SG and OPC models, which have
slopes consistent with zero. However, the slope in our data is over twice that expected by
the OG model for E˙ > 35 ergs s−1 . Below E˙ = 1035 ergs s−1 the expected correlation
between fΩ and E˙ for the OG model becomes more dispersed and the fΩ distribution for
a given E˙ has a tail towards smaller fΩ values, especially for radio quiet pulsars2 . This
feature is also consistent with the experimentally obtained fΩ vs. E˙ distribution. If the two
data points that have lower fΩ ’s (PSR J0633+1746 and PSR J2021+4026) are radio quiet
pulsars I would expect smaller fΩ values compared to radio loud pulsars. However, the radio
1 These best fit parameters were provided in a private conversation by the authors of
[Pierbattista et al.(2012)].
2 Not detected by radio telescopes.
125

Figure 10.4: The distribution of the beaming factor fΩ as a function of spin-down luminosity E˙ for four models derived by [Pierbattista et al.(2012)]. This plot was reproduced
by M. Pierbattista with linear fits for pulsars above E˙ = 1035 ergs s−1 . For the original
figure refer to [Pierbattista et al.(2012)]. Purple and green markers refer to the radio-loud
and radio-quiet pulsars, respectively. Black lines refer to the best linear fits for the pulsars
with E˙ > 1035 ergs s−1 .

126

measurements at 1400 MHz (S1400 ) for these two pulsars are only upper limits: for PSR
J0633+1746 S1400 < 0.507 mJy and for PSR J2021+4026 S1400 < 0.02 mJy. Therefore, it
is unclear whether to categorize them as radio loud or radio quiet pulsars. By considering
both features above and below E˙ ≈ 1035 ergs s−1 , I can conclude that my data sample has
some features which favor the OG model for pulsar emission in the energy range 0.1-100 GeV
over the SG and OPC models, even though the OG model does not quantitatively match
our measurements. However, I cannot reach any conclusions about the PC model, because
of the lack of data-points in the simulated fΩ vs. E˙ distribution.
The discrepancies between my experimental results and the model expectations might be
due to the systematic limitations or inadequacies of Pierbattista et al’s simulations or biases
in my data sample. As Pierbattista et al. mention, results of all four simulated models are
lacking pulsars with E˙ > 3 × 1035 ergs s−1 and characteristic age < 100 kyr, and overpredict the number of low E˙ pulsars. Furthermore they also mention that the simulated
OG model fails to explain many of the most important pulsar population characteristics.
This could certainly effect their fΩ vs. E˙ distributions. Therefore, we can not provide tight
constraints using this synthetic pulsar population. However, a comparison of my results with
an improved model simulation could provide tighter constraints.
My estimates of fΩ may also be biased if there is an unexpected residual dependence
˙ which is hidden in the TeV data scatter, or depending on
of the TeV luminosity on E,
the selected TeV energy range. The energy range of a 1 TeV band around 35 TeV was
chosen because of the Milagro dataset. Ideally, I would prefer to do this study with a more
uniform sample of PWN TeV energy fluxes obtained around the inverse Compton peak of
1
each individual PWN. This result also depends on the relation LP SR GeV ∝ E˙ 2 . In
˙ I can introduce an systematic
order to account for systematic errors of the exponent of E,
127

error of ±0.1 to the exponent of this relation,LP SR GeV ∝ E˙ 0.5±0.1 . After propagating
this systematic error to the slope of log10 (fΩ ) vs. log E˙ distribution, the slope becomes
0.28 ± 0.03statistical ± 0.1systematic . This shows that even after introducing a systematic
error of ±0.1, the slope is still not consistent with zero.

128

Chapter 11

Milagro Observations of Fermi-LAT
Galactic Sources

In 2009 the Milagro collaboration performed a targeted search [Abdo et al.(2009a)] for FermiLAT bright sources (also known as the 0FGL catalog) [Abdo et al.(2009b)] using the Milagro
sky maps and found that Fermi bright sources, that were measured at or above 3 standard
deviations in significance (3σ) by Milagro were dominated by Fermi-LAT bright sources that
have associations with known pulsars. The Fermi-LAT bright source list was made using
the first 3 months of Fermi-LAT data. In 2012 the Fermi-LAT collaboration extended their
source list using the bright sources in the first 2 years of Fermi-LAT data. The new source list
is called the Fermi-LAT 2FGL catalog [Nolan et al.(2012)]. This new catalog increased the
number of Fermi-LAT sources that have associations with known pulsars and in the Milagro
skycoverage by 32. Therefore, I extended the previous Milagro targeted search by making
a candidate list from these 32 new Fermi-LAT sources. Results from the search source in
this new source list is summarized in Table 11.1. In this list, the FDR procedure with a 1%
129

contamination classified 3 Fermi-LAT sources (2FGL J2238.4+5902, 2FGL J2030.0+3640
and 2FGL J1928.8+1740c) as having coincident TeV emission.
The Milagro targeted search on 0FGL sources measured the Milagro flux/flux limit at the
locations of 16 Fermi-LAT sources that were associated with known pulsars [Abdo et al.(2009b)].
Among these 16 sources, 9 passed the standard FDR cut with a 1% contamination. 0FGL
J2055.5+2540, 0FGL J2214.8+3002 and 0FGL J2302.9+443 were categorized as sources of
unknown source type and 0FGL J1954.4+2838 was identified as a source with a spatial association with a known supernova remnant. In the 2FGL catalog these four sources have
been identified with associations to known pulsars, and only 0FGL J1954.4+2838 passed
the standard FDR cut of a 1% contamination. Therefore altogether 52 Fermi-LAT sources
detected by Fermi-LAT have been observed by Milagro and 13 Fermi-LAT sources that have
associations with known pulsars were identified with TeV associations.
One can use this sample to study how the fraction of pulsars FT with a TeV counterpart
changes as a function of the time-averaged GeV flux. FT can be defined as the fraction of
pulsars that passed the standard FDR cut in a given bin of GeV flux,

FT =

Number of FDR true candiates in a given flux bin
.
Total number of candidates in a given flux bin

(11.1)

As shown in Figure 11.1 FT increases with the Fermi flux. The best fit of the FT vs. Fermi
flux assuming a liner function yield:

FT = (0.49 ± 0.09) × (Time-average Fermi flux) + (3.9 ± 0.7),

(11.2)

with χ/N DF = 0.68/3. The fit of the FT vs. Fermi flux assuming a constant has a
130

χ/N DF = 14.22/4. This means that the improvement of the best fit to a linear function
is 3.4σ. This gives evidence that the brighter GeV sources in the 2FGL catalog that have
associations with known pulsars are more likely to have a detectable TeV counterpart.
However, one has to interpret this result with care, because of the known caveats in
the 2FGL catalog. The GeV flux measurements reported in the 2FGL catalog are timeaveraged fluxes of regions with GeV gamma-ray excesses. If there is a pulsar associated with
this region, then the GeV emissions are coming from the pulsar, its associated PWN and
from the Galactic diffuse GeV emission. Therefore, the time-averaged flux is the average
of these 3 components. The ratio of pulsar time-averaged flux to PWN flux is not always
the same. For example, the GeV flux only from the Crab PWN in the energy range of 10
- 316 GeV is 486 ± 188 × 10−12 erg cm−2 s−1 [Acero et al.(2013)] and the time-average
flux from the 2FGL catalog in the same energy range is 9208.41 ± 226 × 10−12 . For 2FGL
J1959.5+2047 the unpulsed component of the time-average GeV flux is not visible after
removing the pulsed component [Acero et al.(2013)]. Thus for the region around the Crab
75% of the time-average GeV flux comes from the time-averaged pulsar emission, but for
2FGL J1959.5+2047 ∼ 100% of the time-average GeV flux is coming from the time-average
pulsar emission. Another disadvantage of the 2FGL flux measurements is the contamination
of Galactic diffuse GeV emission. The best example is 2FGL J1928.8+1740c, which was
measured with 4.03 σ by Milagro. This source was not detected either as a GeV pulsar
[Abdo et al.(2013)] or as a PWN [Acero et al.(2013)] after removing the Galactic diffuse
emission.
One may argue that the dominant GeV flux emission for this set of 2FGL sources is
coming from the time-average pulsar flux, because for the Crab 75% of the flux is coming
from the pulsar and for 2FGL J1959.5+2047 ∼ 100% of the time-average GeV flux is coming
131

from the pulsar. However, the time-averaged flux of a pulsar is a parameter that has to be
interpreted with extra care. The time-averaged flux of a pulsar is the average flux that can
be observed from the earth, but not the true flux emitted by the pulsar. The true flux of the
pulsar is the phase-averaged flux. For some pulsars the width of the pulse is very narrow, for
other pulsars the pulse width is very wide. When the pulse is narrow the time-averaged flux
is considerably less than the phase-averaged flux, but for wide pulses the time-averaged flux
can become roughly equal to the phase-averaged flux. This can be led to misidentify faint
pulsars as bright pulsars and vice versa. The time-averaged flux of J2030+3641 is greater
than the time-averaged flux of J0631+1036. Therefore, one can misidentify J2030+3641 as
a pulsar brighter than J0631+1036. However, J0631+1036 is the brighter pulsar, because
the phase-averaged flux, which is the true flux of the pulsar, is smaller for J2030+3641
than the phase-average flux of J0631+1036. This observation leads me to conclude that the
time-averaged flux is a misleading parameter to determine the brightness of a pulsar.

Since phase-averaged flux is a better estimator for pulsar brightness, one can propose to
remake Figure 11.1 with the phase-averaged fluxes. As shown in Figure 11.2 FT also increases
with the phase-averaged Fermi flux. The best fit of the FT vs. Fermi flux assuming a linear
function yields:

FT = (0.28 ± 0.15) × (Time-average Fermi flux) + (3.3 ± 1.5),

(11.3)

with χ/N DF = 0.95/4. The fit of the FT vs. Fermi flux assuming a constant has a
χ/N DF = 2.16/5. This means that the improvement of the best fit to a linear function is
1σ. This gives a 1.3σ evidence to say that the brighter GeV pulsars are more likely to have
a detectable TeV counterpart. One also should account for the fact that phase-average flux
132

still does not gives the correct estimate for the pulsar brightness. In order to get the correct
brightness, the phase-averaged flux has to be corrected with the beaming factor, fΩ . As
shown in Figure 10.4 the theoretically expected fΩ can be varied within about an order of
magnitude. Therefore, the correction can be quite substantial, which can alter equation 11.3
Altogether, the GeV flux measurements reported in the 2FGL catalog are coming from
3 different GeV sources: pulsars, PWNe and Galactic diffuse emission. Therefore, in the
2FGL catalog one can not identify a single physical source that produces the measured GeV
fluxes. However, this data set gives a 3σ evidence to say that brighter GeV sources in the
2FGL catalog that have associations with known pulsars are more likely to have a detectable
TeV counterpart. As it is shown in this chapter and in chapter 9 there is no evidence to say
that brighter GeV pulsars are more likely to have a detectable TeV counterpart. Therefore,
the correlation between 2FGL GeV sources and their associated TeV emission might be due
to the PWN component or due to the Galactic diffuse emission component. It also possible
that this correlation is a selection effect.

133

FT
1

0.8

0.6

0.4

0.2

0

-8

-7.5
-7
-6.5
-6
Log (Time-averaged Fermi Flux)
10

Figure 11.1: The fraction of Fermi-LAT sources seen by Milagro as a function of half-decade
bins of the integrated time-averaged Fermi flux ( photons cm−2 s−1 ) in the energy range
from 100 MeV to 100 GeV.

134

FT
1

0.8

0.6

0.4

0.2

-11

-10.5
-10
-9.5
-9
-8.5
Log (Phase-average Fermi Flux)
10

Figure 11.2: The fraction of Fermi-LAT sources seen by Milagro as a function of half-decade
bins of the integrated phase-averaged Fermi flux ( photons cm−2 s−1 ) in the energy range
from 100 MeV to 100 GeV.

135

Fermi Name
2FGL

RA
(deg)

DEC
(deg)

l
(deg)

b
(deg)

J0023.5+0924
J0034.4-0534
J0102.9+4838
J0205.8+6448
J0218.1+4233

5.89
8.61
15.74
31.45
34.53

9.41
-5.58
48.65
64.81
42.55

111.5
111.55
124.9
130.74
139.5

-52.85
-68.08
-14.18
3.07
-17.51

J0248.1+6021
J0308.3+7442
J0340.4+4131
J0659.7+1417
J0751.1+1809

42.04
47.08
55.1
104.93
117.78

60.36
74.71
41.53
14.29
18.15

136.89
131.73
153.78
201.05
202.7

0.69
14.23
-11.01
8.27
21.09

<
<
<
<
<

J1023.6+0040
J1142.9+0121
J1301.5+0835
J1312.7+0051
J1549.7-0657

155.92
175.74
195.39
198.18
237.43

0.68
1.35
8.58
0.85
-6.96

243.43
267.56
310.76
314.82
1.23

45.78
59.44
71.3
63.23
35.03

<
<
<
<
<

J1714.0+0751
J1745.6+1015
J1810.7+1742
J1846.4+0920
J1928.8+1740c

258.5
266.4
272.69
281.61
292.22

7.86
10.27
17.7
9.34
17.68

28.84
34.84
44.62
40.7
52.87

25.21
19.23
16.76
5.34
0.03

J1957.9+5033
J1959.5+2047
J2030.0+3640
J2017.3+0603
J2043.2+1711

299.48
299.9
307.51
304.35
310.81

50.56
20.79
36.68
6.05
17.18

84.61
59.18
76.12
48.63
61.9

J2043.7+2743
J2046.7+1055
J2129.8-0428
J2215.7+5135
J2234.7+0945

310.95
311.69
322.47
333.94
338.69

27.72
10.93
-4.48
51.59
9.75

J2238.4+5902
J2239.8+5825 2

339.61
339.97

59.05
58.43

1

Flux/Flux Limit
(×10−17 TeV−1
s−1 cm−2 )
< 22.42
< 54.58
< 22.31
< 20.75
< 41.80

Source Type

Significance
σ

Associated source

psr
PSR
psr
PSR
PSR

-0.73
-2.06
0.51
-0.88
2.94

PSR
PSR
PSR
PSR
PSR

J0023+09
J0034-0534
J0103+48
J0205+6449
J0218+4232

38.32
53.95
14.87
37.31
19.51

PSR
psr
PSR
PSR
PSR

1.54
-0.15
-0.52
1.18
-0.43

PSR
PSR
PSR
PSR
PSR

J0248+6021
J0308+7442
J0340+4130
J0659+1414
J0751+1807

43.66
37.41
21.62
58.01
88.75

psr
psr
psr
psr
psr

-0.33
-0.9
-0.82
0.5
-0.8

PSR
PSR
PSR
PSR
PSR

J1023+0038
J1142+01
J1301+08
J1312+00
J1549-06

< 18.13
< 32.77
< 18.39
< 23.81
46.41±11.50

PSR
psr
psr
PSR
psr

-1.58
0.48
-0.6
-0.54
4.03

PSR J1713+0747
PSR J1745+10
PSR J1810+17
PSR J1846+0919
PSR J1928+1746

10.98
-4.7
-1.45
-16.02
-15.3

< 28.56
< 15.32
42.68±9.55
< 27.2
< 17.82

PSR
PSR
PSR
PSR
PSR

1.41
-0.85
4.46
-0.71
-0.76

PSR J1957+5033
PSR J1959+2048
PSR J2030+3641
PSR J2017+0603
PSR J2043+1710

70.65
57.02
48.93
99.89
76.29

-9.14
-19.57
-36.96
-4.18
-40.43

< 14.62
< 37.49
< 104.35
< 13.28
< 31.44

PSR
psr
psr
psr
psr

-0.72
0.99
0.55
-1.1
0.38

PSR
PSR
PSR
PSR
PSR

106.55
106.41

0.47
-0.16

50.41±11.10
< 51.39

PSR
PSR

4.53
3.01

PSR J2238+5903
PSR J2240+5832

J2043+2740
J2047+10
J2129-04
J2215+51
J2234+09

Table 11.1: Summary of the search with λ = 0.01 for TeV emission from the pulsars in the
2FGL list that were not listed in the 0FGL. (Note that we used the same abbreviations for the
source type as the 2FGL: PSR = Pulsar identified by pulsations and psr = Pulsar identifies
by spatial association.) The Milagro flux derived at 35 TeV is given for the candidates that
passed the λ = 0.01 FDR cut and 95% confidence level upper limits are given for the rest.

136

Chapter 12
Conclusions
In this thesis I presented the GeV-TeV correlations of point-like sources, upper limits to the
Synchrotron Self-Compton model for BL Lac type AGNs, and constraints for the pulsar GeV
emission models. I also discussed the GPS Timing and Control system for the High Altitude
Water Cherenkove (HAWC) observatory, for which I developed the FPGA firmware. HAWC
is a next generation TeV gamma-ray observatory, which will begin to survey the northern
hemisphere of the sky at its full capacity in 2015.
My data sample was obtained from three targeted searches for extragalactic sources and
one targeted search for pulsars, using the existing Milagro sky maps. The detection threshold for these targeted searches was derived using the False Discovery Rate procedure with a
1% contamination fraction. The first extragalactic candidate list was made from extragalactic candidates in the Fermi-LAT 2FGL catalog and the second list form the extragalactic
candidates in the TeVCat catalog. From these two extragalactic lists, Markarian 421 was
the only extragalactic source that Milagro was able to detect. The time averaged flux measurements of Markarian 421 measured by Milagro compared to the IACT measurements of
the transient states of Markarian 421 lead me to conclude that the Markarian 421 stays in
137

the low flux state most of the time. A research paper that summarizes the above results,
of which I was the primary author, was accepted to the Astroparticle Physics Journal. The
third extragalactic candidate list is a list of potential TeV emitting BL Lac candidates that
was compiled using x-ray observations of BL Lac objects and a Synchrotron Self-Compton
model. Milagro’s sensitivity was not enough to detect any of those candidates. However,
the 95% confidence flux upper limits of those sources were above the predicted flux, so the
Synchrotron Self-Compton model used to predict TeV emission is still viable. This work was
published in [Abeysekara & Milagro Collaboration(2013)].
The pulsar candidate list was made from GeV pulsars in the Fermi-LAT second pulsar
catalog. This search was able to measure TeV emission from associated PWNe for 14 GeV
pulsars. The TeV PWN luminosity LP W N T eV vs. E˙ distribution and the LP W N T eV
vs. τc distribution of this sample were both consistent with the leptonic model proposed
by [Mattana et al.(2009)]. For this sample, a linear correlation between pulsar GeV emission and PWN TeV emission was observed, with a linear correlation coefficient of R = 0.82.
This observed GeV –TeV correlation suggests the possibility of a linear relationship between gamma-ray emission mechanisms in pulsars and TeV emission mechanisms in PWNe.
However, it is not possible to explain this linear relationship using the electron-positron populations of curvature radiation inside the magnetosphere and synchrotron radiation in the
PWN. An alternative possibility is a linear relationship between the ambient photon density
in the PWN and the pulsar wind current. A more detailed theoretical study will be needed
to fully understand this correlation.
˙ 2 , TeV efficiency εγ
The correlations between PWN detectability and E/d

˙
and E,

and TeV efficiency and τc were studied for this sample of 14 GeV pulsars. There was no
˙ 2 , which suggests that the current rate of
correlation between PWN detectability and E/d
138

energy a pulsar feeds its associated PWN is not the only parameter relevant for the PWN
˙ and εγ and τc . However, the two
TeV flux. There were no correlations between εγ and E,
PWNe associated with the youngest pulsars have an εγ more than an order of magnitude
less than the other PWNe. This observation supports the argument that TeV photons are
produced by the accumulated electrons in the PWN.
In this thesis I also presented a new multiwavelength technique that I developed to
isolate the flux correlation factor (fΩ ) of GeV pulsars as a function of pulsar spin down
luminosity. This technique requires three input parameters: pulsar spin down luminosity,
pulsar phase-averaged GeV flux, and TeV flux from the associated PWN. The measured
correlation has some features that favor the Outer Gap model over the Polar Cap, Slot
Gap and One Pole Caustic models for pulsar emission in the energy range of 0.1 to 100
GeV, though these simulated models failed to explain many of the most important pulsar
population characteristics. A research paper that summarizes this new technique, of which
I was the first author, was submitted to the Astrophysical Journal.
Altogether my findings support that the TeV emission from PWNe are produced by the
Inverse Compton scattering of the accumulated electrons. Using this property, I developed
a new multiwavelength technique to study GeV pulsar acceleration mechanisms. In this
thesis work I successfully applied this technique to constrain the existing GeV pulsar gap
models. My attempt to constraint a Synchrotron Self-Compton model for BL Lac objects
using the Milagro data set motivate the significance of building more sensitive gamma-ray
observatories, such as HAWC.

139

APPENDIX

140

Appendix
Energy Estimators

Neural Network Based Energy Estimator
An improved energy estimator for Milagro was developed using the Bayesian Neural Network
(BNN) in the ROOT framework [Zhong et al.(2011)]. The BNN was trained using five input
parameters: incident angle, core radius, nFit, nCalMU and A5 (refer to Chapter 4 for the
definitions of these parameters). Both test and training samples were also selected using
two cuts: A5 > 1 and FRASOR > 0.2. The A5 cut removes the hadron like showers in
the data sample, and the FRASOR cut removes the low energy events which barely passed
the Milagro threshold. The Milagro simulated data sample has only about 30000 effective
events. I divided the sample into two and trained the BNN using one sample and used the
other sample for testing. I will call the new energy estimator EBNN.
Figure A.1 shows the comparison between the EBNN obtained using MC simulated data
and are obtained using a sample of real Milagro data. As it can be seen in the top histogram
in Figure A.1, the simulated EBNN follows the general trend of the EBNN obtained with real
data. However, they are not consistent with the one sigma error bars. The bottom histogram
141

in Figure A.1 shows the ratio between the EBNN obtained with a real data sample and the
EBNN obtained with a simulated data sample. If the simulations are ideal this ratio should
equal 1. Data points are closer to 1 in the range from 3.25 to 4.25, but deviated from 1 at
both ends.
Figure A.2 shows the correlations between the energy estimators and the expected true
energy. Note that these plots were made using a simulated data sample, therefore the
expected energy is known. The histogram to the left shows the correlation between EBNN
and expected energy, which has a linear correlation coefficient of 0.84. The histogram to the
right shows the correlation between FRASOR, the old energy estimator, and the expected
energy, which has a liner correlation coefficient of 0.54. The improvement of the linearity is
substantial for the EBNN compared to the existing energy estimator, FRASOR.
Figure A.3 shows the typical dependence of EBNN on the expected energy. Note that
the width of a single EBNN bin covers a wide range of energies and these energies overlap
significantly. The corresponding figure with FRASOR is shown in Figure A.4. Distributions
are narrower for the EBNN compared to the FRASOR.

The Gamma-Table
The Gamma-Tables are a set of tables that summarize the Milagro expected excess of a
Milagro map made with the spectral optimization of E −2.6 for a source with a known
spectrum and at a known location in the sky. Each table summarizes the Milagro expected
excess for a given Milagro epoch and for a source in a given declination. The naming
convention of files is, E<N> dec<ddd>.0.dat, where <N> is the Milagro epoch and <ddd>
is the declination of a source in degrees. For example, “E8 dec42.0.dat” summarizes the
142

Real data
1200

1000

800

600

400

200

0

2

2.5

3

3.5

4

4.5

5

5.5

6
Log(Energy)

4.5

5

5.5

6
Log(Energy)

Data / MC
3

2.5

2

1.5

1

0.5

0

2

2.5

3

3.5

4

Figure A.1: Comparison between Monte Carlo simulated data and real data for cosmic-rays.
In the top figure, the red data points show the distribution of the energy estimator of the
Monte Carlo simulated data. In the top figure, the blue data points show the distribution of
the energy estimator of the real data. In the bottom figure, the blue data points show the
ratio data/monte Carlo simulations.

Correlation of EBNN
0.35

5.5

FRASOR

EBNN

Correlation of FRASOR
2

0.4

1.8

0.3

0.35

1.6

5
0.25
4.5

0.3

1.4

0.25

1.2
0.2
1

0.2

4
0.15
3.5

0.1

0.8

0.15

0.6
0.1
0.4

3

0.05

0.05

0.2
2.5
1

1.5

2

2.5

3

3.5

4

4.5

5
5.5
6
Log(Energy)

0

0
1

1.5

2

2.5

3

3.5

4

4.5

5
5.5
6
Log(Energy)

0

Figure A.2: Correlation between the energy estimators and true energy. Histogram to the left
is the correlation between EBNN and log(Energy). Histogram to the right is the correlation
between FRASOR, the old energy estimator, and log(Energy). Both of these plots were
made using a simulated data sample.

143

EBNN
0.1
0.09
0.08
0.07
0.06
0.05
0.04
0.03
0.02
0.01
0

2

2.5

3

3.5

4

4.5

5

5.5

6
Log(Energy)

Figure A.3: Unit-normalized distribution of true energies for events in EBNN bins, 0.2-0.4,
0.4-0.6, 0.6-0.8, 0.8-1.0, 1.0-1.2, 1.2-1.4, 1.4-1.6, 1.6-1.8 and 1.8-2.0.

FRASOR
0.09
0.08
0.07
0.06
0.05
0.04
0.03
0.02
0.01
0

2

2.5

3

3.5

4

4.5

5

5.5

6
Log(Energy)

Figure A.4: Unit-normalized distribution of true energies for events in FRASOR bins, 0.2-0.4,
0.4-0.6, 0.6-0.8, 0.8-1.0, 1.0-1.2, 1.2-1.4, 1.4-1.6, 1.6-1.8 and 1.8-2.0.

144

Epoch
1
2
3
4
5
6
7
8
9

live time in days
0.
322.2
387.2
443.9
160.2
156.2
181.3
257.9
386.1

Table A.1: Epoch live times

Milagro expected excess per day of Milagro epoch 8 for a source at declination 42.0. Each
table has 4860 rows and 42 columns. The first column is the spectral index and the second
column is the log10 (Energy cutoff) of the simulated source. Every 5th column starting from
column 3 is the expected excess from FRASOR bin 1 through 10. Every 5th column starting
from column 4 is the error of expected excess from FRASOR bin 1 through 10. Every 5th
column starting from column 5 is the mean energy of events that landed in each FRASOR
bin. Every 5th column starting from column 6 is the standard deviation of energy for the
events that landed in each FRASOR bin. For example, column 3, 4, 5 and 6 are the expected
excess, error of the expected excess, mean energy of the events and the standard deviation of
energy for the events landed in FRASOR bin 1, 0 < FRASOR < 0.2. Note that all of these
excess measurements are per day, in order to get the total excess from the full epoch this
number has to be multiplied by the number of live days. Table A.1 summarizes the number
of live days in each epoch.

145

BIBLIOGRAPHY

146

BIBLIOGRAPHY
[Abdo(2007)] Abdo, A. A. 2007, Ph.D. Thesis,
[Abdo et al. (2007)] Abdo, A. A., et al. 2007, ApJ. Lett., 664, L91
[Abdo et al.(2009a)] Abdo, A. A., Allen, B. T., Aune, T., et al. 2009, ApJ. Lett., 700, L127
[Abdo et al.(2009b)] Abdo, A. A., Ackermann, M., Ajello, M., et al. 2009, ApJS, 183, 46
[Abdo et al.(2012)] Abdo, A. A., Allen, B. T., Atkins, R., et al. 2012, ApJ., 750, 63
[Abdo et al.(2013)] Abdo, A. A., Ajello, M., Allafort, A., et al. 2013, ApJS, 208, 17
[Abeysekara & Milagro Collaboration(2013)] Abeysekara, A. U., & Milagro Collaboration
2013, AAS/High Energy Astrophysics Division, 13, #123.01
[Acciari et al.(2009)] Acciari, V. A., Aliu, E., Arlen, T., et al. 2009, ApJ. Lett., 703, L6
[Acciari et al.(2011)] Acciari, V. A., Arlen, T., Aune, T., et al. 2011, ApJ., 729, 2
[Acero et al.(2013)] Acero, F., Ackermann, M., Ajello, M., et al. 2013, ApJ., 773, 77
[Afanas’Ev et al.(2005)] Afanas’Ev, V. L., Dodonov, S. N., Moiseev, A. V., et al. 2005,
Astronomy Reports, 49, 374
[Aharonian et al.(2005)] Aharonian, F., Akhperjanian, A. G., Aye, K.-M., et al. 2005, A.&
A., 435, L17
[Aharonian et al.(1999)] Aharonian, F. A., Akhperjanian, A. G., Barrio, J. A., et al. 1999,
A.& A., 349, 11
[Aharonian et al.(2001)] Aharonian, F. A., Akhperjanian, A. G., Barrio, J. A., et al. 2001,
A.& A., 366, 62
[Aharonian et al.(2005)] Aharonian, F., Akhperjanian, A., Beilicke, M., et al. 2005, A.& A.,
431, 197
147

[Aharonian et al.(2006a)] Aharonian, F., Akhperjanian, A. G., Bazer-Bachi, A. R., et al.
2006, A.& A., 457, 899
[Aharonian et al.(2006b)] Aharonian, F., Akhperjanian, A. G., Bazer-Bachi, A. R., et al.
2006, A.& A., 448, L43
[Aharonian et al.(2006c)] Aharonian, F., Akhperjanian, A. G., Bazer-Bachi, A. R., et al.
2006, ApJ., 636, 777
[Aharonian et al.(2006d)] Aharonian, F., Akhperjanian, A. G., Bazer-Bachi, A. R., et al.
2006, A.& A., 456, 245
[Aharonian et al.(2007)] Aharonian, F. A., Akhperjanian, A. G., Bazer-Bachi, A. R., et al.
2007, A.& A., 469, L1
[Aharonian et al.(2012)] Aharonian, F. A., Bogovalov, S. V., & Khangulyan, D. 2012, Nature, 482, 507
[H. E. S. S. Collaboration et al.(2007)] H. E. S. S. Collaboration, :, Djannati-Atai, A., et al.
2007, arXiv:0710.2247
[Albert et al.(2007)] Albert, J., Aliu, E., Anderhub, H., et al. 2007, ApJ., 663, 125
[Aliu et al.(2013)] Aliu, E., Archambault, S., Arlen, T., et al. 2013, ApJ., 764, 38
[Arons(1996)] Arons, J. 1996, A&AS, 120, 49
[Arons, J. 2006] Arons, J. 2006, On the Present and Future of Pulsar Astronomy, 26th
Meeting of the IAU, Join Discussion 2, JD02, 39, 2
[Atkins et al.(2004)] Atkins, R., Benbow, W., Berley, D., et al. 2004, ApJ., 608, 680
[Atoyan & Aharonian(1996)] Atoyan, A. M., & Aharonian, F. A. 1996, MNRAS, 278, 525
[Atwood et al.(2009)] Atwood, W. B., Abdo, A. A., Ackermann, M., et al. 2009, ApJ., 697,
1071
[Carrigan et al.(2008)] Carrigan, S., Hinton, J. A., Hofmann, W., & et al. 2008, International
Cosmic Ray Conference, 2, 659
148

[Cheng & Zhang(1996)] Cheng, K. S., & Zhang, J. L. 1996, ApJ., 463, 271
[Costamante & Ghisellini(2002)] Costamante, L., & Ghisellini, G. 2002, A.& A., 384, 56
[Dwek & Krennrich(2005)] Dwek, E., & Krennrich, F. 2005, ApJ., 618, 657
[de Jager & Djannati-Ata¨ı(2009)] de Jager, O. C., & Djannati-Ata¨ı, A. 2009, Astrophysics
and Space Science Library, 357, 451
[de Jager et al.(2009)] de Jager, O. C., Ferreira, S. E. S., Djannati-Ata¨ı, A., et al. 2009,
arXiv:0906.2644
[Gabici et al.(2009)] Gabici, S., Aharonian, F. A., & Casanova, S. 2009, MNRAS, 396, 1629
[Gibbs et al.(2003)] Gibbs, K., Criswell, S., Falcone, A., et al. 2003, International Cosmic
Ray Conference, 5, 2823
[Glendenning(1992)] Glendenning, N. K. 1992, Phys. Rev. D, 46, 4161
[Gustavo E.R. and K.S. Cheng (2004)] in Cosmic Gamma-ray Sources, ed. K.S. Cheng 2004
Astrophysics and Space Science Library, 2004
[Harding(1981)] Harding, A. K. 1981, ApJ., 245, 267
[Harding & Muslimov(2002)] Harding, A. K., & Muslimov, A. G. 2002, ApJ., 568, 862
[Harko & Cheng(2002)] Harko, T., & Cheng, K. S. 2002, MNRAS, 335, 99
[Hays(2004)] Hays, E. A. 2004, Ph.D. Thesis
[Hewish et al.(1968)] Hewish, A., Bell, S. J., Pilkington, J. D. H., Scott, P. F., & Collins,
R. A. 1968, Nature, 217, 709
[Hinton(2004)] Hinton, J. A. 2004, New Astron. Rev., 48, 331
[Klebesadel et al.(1973)] Klebesadel, R. W.,Strong, I. B., & Olson, R. A. 1973, ApJ. Lett.,
182, L85
[Konopelko et al.(2003)] Konopelko, A., Mastichiadis, A., Kirk, J., de Jager, O. C., &
Stecker, F. W. 2003, ApJ., 597, 851
149

[Krennrich et al.(2001)] Krennrich, F., et al. 2001, ApJ. Lett., 560, L45
[Kurczynski & Gawiser(2010)] Kurczynski, P., & Gawiser, E. 2010, AJ., 139, 1592
[Lattimer et al.(1990)] Lattimer, J. M., Prakash, M., Masak, D., & Yahil, A. 1990, ApJ.,
355, 241
[Lattimer & Prakash(2001)] Lattimer, J. M., & Prakash, M. 2001, ApJ., 550, 426
[Li & Ma(1983)] Li, T.-P., & Ma, Y.-Q. 1983, ApJ., 272, 317
[Lorimer & Kramer(2012)] Lorimer, D. R., & Kramer, M. 2012, Handbook of Pulsar Astronomy, by D. R. Lorimer , M. Kramer, Cambridge, UK: Cambridge University Press,
2012
[Mattana et al.(2009)] Mattana, F., Falanga, M., G¨otz, D., et al. 2009, ApJ., 694, 12
[Mayer et al.(2012)] Mayer, M., Brucker, J., Holler, M., et al. 2012, arXiv:1202.1455
[Miller et al.(2001)] Miller, C. J., et al. 2001, AJ., 122, 3492
[Muslimov & Harding(2003)] Muslimov, A. G., & Harding, A. K. 2003, ApJ., 588, 430
[Harding & Muslimov(2003)] Harding, A. K., & Muslimov, A. G. 2003, Pulsars, AXPs and
SGRs Observed with BeppoSAX and Other Observatories, 121
[Muslimov & Harding(2004)] Muslimov, A. G., & Harding, A. K. 2004, ApJ., 606, 1143
[Nolan et al.(2012)] Nolan, P. L., Abdo, A. A., Ackermann, M., et al. 2012, VizieR Online
Data Catalog, 219, 90031
[Oppenheimer & Volkoff(1939)] Oppenheimer, J. R., & Volkoff, G. M. 1939, Physical Review, 55, 374
[Pacini & Salvati(1973)] Pacini, F., & Salvati, M. 1973, ApJ., 186, 249
[Pareschi et al.(2013)] Pareschi, G., Armstrong, T., Baba, H., et al. 2013, Proceedings of the
SPIE, 8861,
150

[Pierbattista et al.(2012)] Pierbattista, M., Grenier, I. A., Harding, A. K., & Gonthier, P. L.
2012, A.& A., 545, A42
[Punch et al.(1992)] Punch, M., Akerlof, C. W., Cawley, M. F., et al. 1992, Nature, 358, 477
[Romani & Watters(2010)] Romani, R. W., & Watters, K. P. 2010, ApJ., 714, 810
[Swanenburg et al.(1981)] Swanenburg, B. N., Bennett, K., Bignami, G. F., et al. 1981, ApJ.
Lett., 243, L69
[Takata et al.(2010)] Takata, J., Wang, Y., & Cheng, K. S. 2010, ApJ., 715, 1318
[Thompson, D.J. (2000)] in High Energy Gamma-Ray Astronomy, ed. F.A. Aharonian & H.
J. Volk (AIP:New Yor), 103
[Thompson(2001)] Thompson, D. J. 2001, American Institute of Physics Conference Series,
558, 103
[VERITAS Collaboration et al.(2011)] VERITAS Collaboration, Aliu, E., Arlen, T., et al.
2011, Science, 334, 69
[Wakely & Horan(2008)] Wakely, S. P., & Horan, D. 2008, International Cosmic Ray Conference, 3, 1341
[Weekes et al.(1989)] Weekes, T. C., Cawley, M. F., Fegan, D. J., et al. 1989, ApJ., 342, 379
[Zhang & Yuan(1998)] Zhang, J. L., & Yuan, Y. F. 1998, ApJ., 493, 826
[Zhong et al.(2011)] Zhong, J., Huang, R.-S., & Lee, S.-C. 2011, Computer Physics Communications, 182, 2655

151