CONSTRAINTS ON LORENTZ-INVARIANCE VIOLATION WITH THE

HAWC OBSERVATORY

By

Samuel Stephens Marinelli

A DISSERTATION

Submitted to

Michigan State University

in partial fulfillment of the requirements

for the degree of

Physics — Doctor of Philosophy

2019

ABSTRACT

CONSTRAINTS ON LORENTZ-INVARIANCE VIOLATION WITH THE

HAWC OBSERVATORY

By

Samuel Stephens Marinelli

Lorentz invariance is one of the fundamental symmetries assumed by the Standard Model

(SM) of particle physics. However, many extensions to the SM allow or predict Lorentz-

invariance violation (LIV), in which this symmetry is either spontaneously broken or is only

approximately valid at low energies. Under models in which LIV occurs in the photon

sector, photon decay to electron/positron pair becomes allowed at sufficiently high energies.

Astrophysical γ (gamma) rays propagate over sufficient distances before reaching Earth that

very nearly all of them should decay above the energy at which the decay becomes allowed.

The detection of high-energy photons can therefore be used to constrain such models.

The High-Altitude Water-Cherenkov (HAWC) observatory is a γ-ray telescope located on

the Sierra Negra volcano in the state of Puebla, Mexico. The detector consists of 300 water-

filled tanks instrumented with four photomultiplier tubes each, which detect the Cherenkov

light produced by charged particles in extensive air showers initiated by TeV cosmic γ rays.

This thesis presents an analysis of high-energy emission from several HAWC γ-ray sources

and sets new constraints on LIV-induced photon decay. This is done using a new algorithm

for estimating the energies of γ rays detected by HAWC, utilizing a neural network (NN)

trained on Monte Carlo (MC) simulations of the detector. The NN uses a multilayer per-

ceptron architecture and maps several geometric features of the air shower to an estimate of

the primary photon’s energy.

HAWC’s high-energy sensitivity combined with its improved energy resolution using the

NN allows an LIV-induced cutoff in γ-ray spectra below 209 TeV to be ruled out with 95%
confidence. This corresponds to limits on the quantum-gravity energy scale of 8.74 × 1030
or 4.27 × 1022 eV for LIV perturbing the photon dispersion relation at first or second order

respectively. These limits constitute a new best constraint on photon-sector LIV.

Copyright by
SAMUEL STEPHENS MARINELLI
2019

ACKNOWLEDGMENTS

Thanks to my thesis advisor and mentor, Jim Linnemann. You taught me about all of

the technical facets of the work: astrophysics, particle physics, HAWC, statistics, machine

learning—the list goes on. But you also taught me how to conduct myself as a scientist and

how to communicate my ideas and results. Thanks for all your patience and guidance.

Thanks to my thesis committee, Kirsten Tollefson, Gus Sinnis, Ed Brown, and Scott

Pratt. You helped me plan out a thesis about SNRs and didn’t bat an eye when I came back

with a thesis about LIV.

Thanks to Tilan Ukwatta and Tolga Yapıcı, who were postdocs in our group while I was

here. You set an example for me as a scientist that I still strive to emulate to this day.

Thanks to my predecessor, Udara Abeysekara. You were my first guide to HAWC, and

you left some seriously big shoes to fill. And thanks to my successors, Joe Lundeen, Alison

Peisker, and Dan Salazar-Gallegos. You tolerated all my jargon and acronyms. Now go fill

my shoes (size 81⁄2).

You’re welcome to all the young people in the collaboration who got to hang out with

me at collaboration meetings and read my musings on Slack. I don’t know what you’ll do

without me.

Thanks to everyone in the HAWC collaboration. You were the giants whose shoulders

I got to stand on. Thanks to John Pretz and Andy Smith for all the wisdom and advice

on data, algorithms, and everything else. Thanks to Humberto Mart´ınez-Huerta and Pat

Harding for letting your LIV theory knowledge diffuse into my head (albeit slowly).

Thanks to my parents, George and Kathy Marinelli. You encouraged me literally since

birth, and it’s no exaggeration to say none of this work would have been done without you.

v

And thanks to Xueying Huyan, who put up with the travel, the nights and weekends

spent working late, and the endless rants about software. You were a constant source of

comfort and support while I was writing this document. Now go write your own!

vi

TABLE OF CONTENTS

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

LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Chapter 1

Lorentz-invariance violation . . . . . . . . . . . . . . . . . . . . . .
1.1 The Lorentz transformation . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2 Violation phenomonology . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3 Photon decay . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Chapter 2

γ-ray production at PWNe and SNRs . . . . . . . . . . . . . . .
2.1 Pulsar properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2 Pulsar evolution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3 Mechanisms of hadronic acceleration . . . . . . . . . . . . . . . . . . . . . .
2.4 Models of leptonic acceleration . . . . . . . . . . . . . . . . . . . . . . . . . .

Chapter 3 The HAWC observatory . . . . . . . . . . . . . . . . . . . . . . . .
3.1 Atmospheric air showers . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2 Water-Cherenkov detectors . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3 Data-acquisition system . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.4 Science goals
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Chapter 4 Air-shower reconstruction . . . . . . . . . . . . . . . . . . . . . . .
4.1 Calibration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2 Edge finding . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.3 Core reconstruction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.4 Angle reconstruction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.5 Background rejection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.6 MC simulations of air showers . . . . . . . . . . . . . . . . . . . . . . . . . .

5.1.1

Chapter 5 Energy estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.1 Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Input Variables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.1.1.1 Multiplicity . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.1.1.2 Containment
. . . . . . . . . . . . . . . . . . . . . . . . . .
5.1.1.3 Atmospheric Attenuation . . . . . . . . . . . . . . . . . . .
5.1.2 Training . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.1.3 Resulting Functional Form . . . . . . . . . . . . . . . . . . . . . . . .
5.2 Performance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2.1 Other Energy Variables . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2.1.1 Fraction Hit/Fraction-Hit NN . . . . . . . . . . . . . . . . .
5.2.1.2
. . . . . . . . . . . . . . . . .

Likelihood Energy Estimator

ix

xi

1
1
4
5

9
9
12
14
19

22
22
24
28
31

35
35
39
40
42
44
46

48
48
49
49
50
50
52
53
56
56
57
57

vii

5.2.1.3 Ground Parameter Energy Estimator . . . . . . . . . . . . .
5.2.1.4 Ground Energy Estimator . . . . . . . . . . . . . . . . . . .
5.2.2 Linearity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2.3 Bin Overlap . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2.4 Bias
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2.5 Resolution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2.6 RMS Error
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2.7 Bin Contamination . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2.8 Agreement of Estimators . . . . . . . . . . . . . . . . . . . . . . . . .

Chapter 6 Crab Nebula spectral reconstruction . . . . . . . . . . . . . . . .
6.1 Skymaps . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.2 Binned-likelihood formalism . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.3 Flux-points calculation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.4 Fit result
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Chapter 7

LIV limits . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.1 Analysis technique
7.2 Results and interpretation . . . . . . . . . . . . . . . . . . . . . . . . . . . .

57
57
58
59
63
69
69
71
71

75
75
77
85
86

92
92
99

Chapter 8 Conclusions and future prospects . . . . . . . . . . . . . . . . . . 105

APPENDIX . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108

BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125

viii

LIST OF TABLES

Table 5.1: NN input variables that characterize energy. The annuli are about the
shower axis. Annuli 0 – 8 are each 10 m thick. “Annulus” 9 extends from
90 m to infinity.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Table 6.1: Best-fit values of log-parabola spectral parameters with statistical and sys-
tematic uncertainties. Systematic uncertainties are the sum in quadrature
of the deviations from the nominal result produced by performing the anal-
. . . . . . . .
ysis using the systematics MCs as described in Section 6.3.
Table 6.2: Energy (in TeV) and energy squared times flux (in TeV/(cm2 · s)) for each
point in Figure 6.7. Statistical errors are propagated from the uncertainties
on the flux normalization as refit in each energy bin as described in Section
6.3. Systematic errors are sums in quadrature of deviations of the flux
points computed using systematics MCs from the nominal flux points.
. .

49

87

88

Table 7.1: Top-hat radii for the seven sources. . . . . . . . . . . . . . . . . . . . . . .

94

Table 7.2: Best-fit values of the spectral parameters in the LIV model. The σ param-

eter is fixed for each source.

. . . . . . . . . . . . . . . . . . . . . . . . . 102

Table 7.3: p-values and 95%-confidence limits on ELIV. Limits are expressed in TeV.
The limit using the nominal MC is provided along with a systematic error
obtained by adding in quadrature the deviations of the limits calculated
using the systematics MCs. The systematics MCs are those having a lower
than nominal charge smearing, higher charge smearing, lower QE, higher
QE, lower PMT threshold, and no broad-pulse effect. The systematics-
degraded limit is then computed as the nominal limit minus its lower error.
See Section 4.6 for more details on the systematics MCs.

. . . . . . . . . 103

Table 7.4: 95%-confidence lower limits on the QG energy EQG for n values of 1 and 2.
All limits are expressed in eV. These limits can be computed from the ones
given in the nominal degraded column of Table 7.3 using Equation (1.18).

103

Table 7.5: Systematics-degraded limits on ELIV at various confidences. Limits are
expressed in TeV. A Zσ limit is defined to have confidence F (Z) where F
is the standard-normal cumulative distribution. As a result, the 3σ, 4σ, and
5σ limits correspond to confidences of 99.87%, 99.9968%, and 99.999971%
respectively.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104

ix

Table 7.6: Observed and expected bin contents from the 2HWC J1825−134 analysis,
as plotted in Figure 7.4. Columns are central reconstructed bin energy, ob-
served count, background expectation, signal-plus-background expectation
with ELIV set to its 95% lower limit of 259 TeV, and signal-plus-background
expectation with ELIV set to infinity.

. . . . . . . . . . . . . . . . . . . . 104

x

LIST OF FIGURES

Figure 3.1: An atmospheric air shower. Dots indicate shower particles. Image taken
. . . . . . . . . . . . . . . . . . . . . . . . . .

from [1] with permission.

25

Figure 3.2: Left: a HAWC tank with a charged particle passing through (the thick
line coming from above) and the resulting Cherenkov cone (the thinner
lines). Right: the layout of HAWC’s 300 tanks (large circles) and their
PMTs (smaller filled circles). The gap in the center is the location of the
. . .
counting house, which holds the electronics. Image taken from [2].

Figure 3.3: The output voltage from the PMT (top) and the resulting FEB output.
From this output the time over threshold (ToT) for each threshold can be
extracted and used to estimate the number of PEs produced in the PMT.
Image taken with permission from [1].
. . . . . . . . . . . . . . . . . . .

27

29

Figure 3.4: Diagram of the HAWC data-acquisition system [3].

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

32

Figure 4.1: Slewing times shown as functions of low and high ToTs. The curves show
. . . . . . . . . . . . . . . . . . . .

the fit function from Equation (4.1).

37

Figure 4.2: The effective charge in each PMT in a γ-ray event. Image taken from [2].

41

Figure 4.3: The arrival time of each hit used in reconstructing the direction from
which the example event in Figure 4.2 originated. Image taken from [2].

Figure 4.4: The lateral distribution of charge within a background event and a signal
. . . . . . . . . . . . . . . . . . . . . . . .

event. Image taken from [2].

Figure 5.1: Dependence of NN energy estimate on each input variable, with the others
. . . . . . . . . . . . . . . . . . . . . . .

frozen at their median values.

43

45

54

Figure 5.2: Correlation between energy variables and MC-truth energy.

. . . . . . .

60

Figure 5.3: Gaussians representing distributions of MC-truth energies in reconstructed-
energy or energy-proxy bins. The resulting bin purity (fraction of events
not migrating into adjacent bins) is shown in Figure 5.8.
. . . . . . . . .

64

Figure 5.4: Distribution of MC-truth energies in energy bins.

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

66

Figure 5.5: Bias of energy variables.

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

68

Figure 5.6: Resolution of energy variables. . . . . . . . . . . . . . . . . . . . . . . . .

70

xi

Figure 5.7: RMS error of energy variables. . . . . . . . . . . . . . . . . . . . . . . . .

72

Figure 5.8: Fraction of MC events reconstructed in each NN bin whose MC-truth
. . . . . . . . . . . . . . . . . . . . . .

energies actually lie in that bin.

73

Figure 5.9: Correlations of three energy variables.

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

74

Figure 6.1: The skymap for fraction-hit bin 3 with NN energies between 3.16 and
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5.62 TeV.

Figure 6.2: The background skymap for fraction-hit bin 3 with NN energies between
. . . .

3.16 and 5.62 TeV. This bin uses the DI background calculation.

Figure 6.3: The skymap for fraction-hit bin 9 with NN energies between 56.2 and
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

100.0 TeV.

Figure 6.4: The background skymap for fraction-hit bin 9 with NN energies between
. . . . . .

56.2 and 100.0 TeV. This bin uses background randomization.

Figure 6.5: The 68% quantile of the distribution of the angular distance between the
. . . . . . .

true and reconstructed directions of events in data and MC.

Figure 6.6: The Crab log-parabola spectrum fit using HAWC data binned in NN
energy, along with similar measurements made by the H.E.S.S. [4] and
HEGRA [5] experiments. Errors on the HAWC spectrum are statistical
. . . . . . . . . . . . . . . . .
only. See Figure 6.7 for systematic errors.

Figure 6.7: The Crab energy spectrum produced by the flux-points analysis described
in Section 6.3, along with measurements by H.E.S.S. [4] and HEGRA [5].
. . . . . . . .
Errors are statistical and systematic added in quadrature.

Figure 6.8: Pulls (residuals divided by their expected fluctuations) for each of the 2D
bins included in the Crab analysis. The x axis indicates NN bins; the
color scale indicates fraction-hit bins.
. . . . . . . . . . . . . . . . . . .

Figure 7.1: A log-parabola energy spectrum subject to a hard cutoff. The blue curve
shows the true-energy spectrum; the orange curve shows the reconstructed-
energy distribution.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 7.2: 95%-confidence lower limit on ELIV using 2HWC J1825-134 data as a
. . . . . . . . . . . . . . . . . . . . . . . . .

function of top-hat radius.

Figure 7.3: Two times the natural logarithm of the profile likelihood defined in Equa-
. . . . . .

tion (7.4), plotted as a function of ELIV for the Crab Nebula.

78

79

80

81

84

88

89

90

93

95

98

xii

Figure 7.4: Expected and observed contents of the reconstructed-energy bins from the
2HWC J1825−134 limit-setting procedure. The orange points show the
expected distribution under the signal-plus-background model in which
the LIV cutoff energy is set to infinity and the other parameters are opti-
mized. The green points show the expected distribution under the signal-
plus-background model with the hard cutoff set to its 95%-confidence
lower limit and the parameters optimized. The red points indicate the
expectations from background only.

. . . . . . . . . . . . . . . . . . . . 100

Figure 8.1: Mean free path of photons subject to the CMB [6].

. . . . . . . . . . . . 107

Figure .1: Data/MC comparison of fraction hit. p = 0.0997.
Figure .2: Data/MC comparison of fraction of tanks hit. p = 1.10 × 10−5.

. . . . . . . . . . . . . 110

. . . . . 111

Figure .3: Data/MC comparison of log of core-fit amplitude. p = 0.0538.

. . . . . . 112

Figure .4: Data/MC comparison of core distance from array center. p = 0.0538.

. . 113

Figure .5: Data/MC comparison of cosine of zenith angle. p = 0.715.

. . . . . . . . 114

Figure .6: Data/MC comparison of fraction of charge contained in annulus 0 . p =

0.0762. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115

Figure .7: Data/MC comparison of fraction of charge contained in annulus 1 . p =

0.116.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116

Figure .8: Data/MC comparison of fraction of charge contained in annulus 2 . p =

0.550.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117

Figure .9: Data/MC comparison of fraction of charge contained in annulus 3 . p =

0.929.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118

Figure .10: Data/MC comparison of fraction of charge contained in annulus 4 . p =

0.165.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119

Figure .11: Data/MC comparison of fraction of charge contained in annulus 5 . p =

0.0118. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120

Figure .12: Data/MC comparison of fraction of charge contained in annulus 6 . p = 0.11.121

Figure .13: Data/MC comparison of fraction of charge contained in annulus 7 . p =

0.0196. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122

xiii

Figure .14: Data/MC comparison of fraction of charge contained in annulus 8 . p =

1.34 × 10−7. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123

Figure .15: Data/MC comparison of fraction of charge contained in annulus 9. p =

4.54 × 10−15.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124

xiv

Chapter 1

Lorentz-invariance violation

Symmetry principles are fundamental to the modern understanding of the laws of physics.

Many diverse properties of the Universe can be seen as arising out of the symmetries that

characterize it. Momentum conservation is guaranteed by space-translation invariance; en-

ergy conservation is guaranteed by time-translation invariance; the gauge bosons’ existence is

necessitated by their corresponding gauge symmetries. And one such fundamental symmetry

assumption is that the Universe obeys Lorentz symmetry, a symmetry group corresponding

to invariance of the laws of physics under rotations and changes of reference frame. Lorentz

invariance is an assumption of the Standard Model (SM) of particle physics, but there are

reasons to believe it may be violated, as described in Section 1.2. The subject of this thesis

is a search for evidence of the violation of the Lorentz symmetry group, carried out using

the High-Altitude Water-Cherenkov (HAWC) observatory, described in Chapter 3.

1.1 The Lorentz transformation

The Lorentz group is the group of linear transformations that can be used to transform

quantities in one inertial reference frame into the values they assume in another. Let the

spacetime interval between two events, or locations in spacetime, be defined as

s ≡

(c∆t)2 − |∆r|2

(1.1)

(cid:113)

1

where c is the speed of light, t is time, r is the position 3-vector, and ∆ represents a difference

taken between the two events. It follows from the axioms of special relativity that the interval

between two events is invariant under a change of reference frame [7]. The transformations

of the Lorentz group must therefore have the property that they leave the interval between

any two events unchanged. The only such transformations are translations and rotations

(where rotations involving the time direction must be hyperbolic as described below) [7].

Since translations are affine transformations rather than linear transformations, they do not

fall into the Lorentz group, and so the Lorentz group consists of rotations in 4-dimensional

spacetime1.

The following derivation of the functional form of the Lorentz transformation follows

[7]. An arbitrary rotation in 4 dimensions can be expressed as a composition of (cid:0)4

(cid:1) = 6

2

rotations such that each rotation occurs in a plane containing a different pair of the axes

of the coordinate system (i.e. the tx-, ty-, tz-, xy-, xz-, and yz-planes). The xy-, xz-, and

yz-rotations are ordinary circular rotations and can be interpreted as reorientations of the

spatial coordinates of the reference frame. For example, a rotation by an angle θ in the

xy-plane can be represented as

r(cid:48) =

1

 r

0

0

0
0 cos θ − sin θ 0
0
0 sin θ
0
1

cos θ

where r is the 4-position (cid:2)ct x y z(cid:3)T and the prime indicates a transformed quantity.

0

0

(1.2)

However, the minus sign in Equation (1.1) implies that, in order to leave the interval un-

changed, rotations involving the time direction must be hyperbolic rather than circular. For

1The full group of transformations which leave the interval invariant and therefore represent changes of
inertial reference frame is called the Poincar´e group or the inhomogeneous Lorentz group and consists of
translations and rotations (and compositions thereof). These transformations are all affine.

2

example, the transformation in the tx-plane is

 cosh ψ − sinh ψ 0 0

− sinh ψ

cosh ψ

 r

(1.3)

r(cid:48) =

0
0

0
0

0 0
1 0
0 1

for some real parameter ψ, called the rapidity of the transformation. Such a transformation

is known as a Lorentz boost and corresponds to a change from one inertial frame to another

which has some (constant) velocity with respect to the first. Let the velocity in question

be V , directed in the x-direction. We can attempt to find the relation between ψ and V

by considering the movement of the spatial origin of the boosted frame in both coordinate

systems. Substituting r =(cid:2)ct V t 0 0(cid:3)T and r(cid:48) =(cid:2)ct(cid:48) 0 0 0(cid:3)T into Equation (1.3) and

considering the x-component of the resulting vector equation produces

0 = −ct sinh ψ + V t cosh ψ,

which implies that the boost speed and rapidity are related by

tanh ψ =

≡ β,

V
c

and so the boost transformation, in terms of the scaled boost velocity β, is

 Γ −βΓ 0 0
 r

−βΓ
0
0

0 0
1 0
0 1

Γ
0
0

r(cid:48) =

where Γ is the Lorentz factor, defined by

1(cid:112)

1 − β2

.

Γ ≡

3

(1.4)

(1.5)

(1.6)

(1.7)

1.2 Violation phenomonology

In general, Lorentz-invariance violation (LIV) can be introduced by the addition of a not

Lorentz-invariant term into the free-particle Lagrangian of a theory [8]. This results in a

modification to the theory’s free-particle dispersion relation of the following form:

E2 − p2c2 + f (q) q2 = m2c4

(1.8)

where q is either pc or E and f is an arbitrary dimensionless function. In the ultrarelativistic
limit, where p (cid:29) mc, the choice of either pc or E for q is equivalent. f should be small,

at least at low energies, so that the low-energy physics remains unperturbed, reproducing

the Lorentz-invariant Universe we observe. In order to consider the behavior in the energy

regime in which LIV is just becoming perceptible, it is desirable to simplify Equation (1.8)

by Taylor-expanding f about q = 0 and keeping only the lowest-order nonzero term in q,

resulting in

E2 − p2c2 ± αnqn+2 = m2c4

(1.9)

where n is the order of the lowest-order nonzero term in the Taylor expansion, henceforth

called the order of the violation, and αn is a positive real number called the LIV coefficient

of order n. Models of quantum gravity (QG) that induce LIV are expected to produce a

coefficient

αn ∼ E−n
QG,

(1.10)

where EQG is the energy at which quantum-gravitational effects become nonnegligible. EQG

is in turn expected to be of order 1028 eV, corresponding to the Planck energy scale [9].

For photons, m = 0, and the sign of the LIV term in Equation (1.9) therefore determines

4

whether the particle’s energy or momentum is greater. This has a qualitative effect on the

allowed physics; in particular, choosing the negative sign means that the photon’s energy

exceeds its momentum, and when they differ sufficiently, the extra energy can be used to

create a electron/positron pair as described in Section 1.3.

Thus far we have characterized LIV generically. However, a theoretical framework called

the minimal Standard Model Extension (mSME) has been developed that specifically exam-

ines the functional form of the Lagrangian density associated with such Lorentz-symmetry-

breaking phenomena [10]. The mSME Lagrangian adds to the SM Lagrangian the minimal

set of non-Lorentz-invariant terms expressible as products of SM and gravitational fields that

maintain gauge invariance and renormalizability. Some of these break CPT symmetry. Since

we will be concerned with LIV in the photon sector, we will assume only the photon-sector

contributions to the Lagrangian, which, in addition to the SM terms, are the following:

LLIV = −1
4

(kF )ρλµν F ρλF µν +

1
2

(kAF )ρ 
ρλµνAλF µν

(1.11)

where kF and kAF are constant coefficients which are free parameters of the model, A is
the photon field, F µν ≡ ∂µAν − ∂νAµ, and 
 is the Levi-Civita symbol. The left term of

Equation (1.11) is CPT-even, and the right is CPT-odd [9].

1.3 Photon decay

Under LIV models in which Lorentz-invariance-violating term is preceeded by a minus sign

in Equation (1.9), photon decay to an electron-positron pair becomes possible. This is only

allowed at very high energies, and so this statement is compatible with the fact that photon

5

decay is apparently forbidden at all energies thus far observed. Here we will derive the energy

at which such a decay becomes allowed. The derivation follows [9].

The dispersion relation for a photon with energy ω and wavenumber k under the generic

LIV framework described in Section 1.2 is

ω2 = k2 (1 + αnkn) .

(1.12)

Let p be the momentum of the positron produced in the decay, and let θ be half the angle

between the emitted positron and electron. Then energy conservation implies

k

1 + αnkn =

m2

e + p2 +

=

m2

e + p2 +

e + |k − p|2
m2
e + p2 + k2 − 2pk cos θ.
m2

(cid:113)
(cid:113)

Squaring this equation, rearranging, and squaring again gives the condition

αnkn + 2 sin2 θ

p2 − 2αnkn+1p cos θ + 2m2

e + αnknm2

e = 0,

(cid:112)

(cid:16)

(cid:113)
(cid:113)

(cid:17)

(1.13)

(1.14)

(1.15)

.

(1.16)

which admits solutions

αnkn+1 cos θ ±(cid:113)

nk2n+2 cos2 θ − 4(cid:0)sin2 θ + αnkn(cid:1) (1 + αnkn) m2

e

α2

p± =

2αnkn + 2 sin2 θ

In order for the decay to be allowed, the positron momentum must be real, and so the

discriminant must be nonnegative. For a given value of the LIV coefficient αn and a given

photon momentum k, this may be possible at all θ in [0, π], at some values of θ, or at no

values. In the last case, the decay is forbidden. The critical relationship between the LIV

6

coefficient and photon momentum in order for the decay to be allowed is

αn ≥

kn(cid:0)k2 − 4m2

4m2
e

e

(cid:1) .

(1.17)

When the left- and right-hand sides of Equation (1.17) are exactly equal, the discriminant
of Equation (1.16) assumes a value of −4m2

e)2, which is negative for all

ek4 sin2 θ/(k2 − 4m2

opening angles except 0 and π, where it is 0, and so the decay is allowed only for particles

emitted in exactly the same direction or in opposite directions. When the greater-than

relationship is satisfied, more angles are allowed around 0 and π.

This behavior, in which photon decay is prohibited below a particular αn-dependent

momentum, allows the value of αn to be constrained by observations of high-energy photons

because the decay probability of photons from relevant astrophysical sources is precisely zero

below the critical momentum and, as will be shown in Chapter 7, very nearly one above it.

Such a limit in turn implies a constraint on the QG energy scale via Equation (1.10):

(cid:115)
EQG ≥ k n

− 1.

k2
4m2
e

(1.18)

The following chapter, Chapter 2, describes pulsar wind nebulae (PWNe) and supernova

remnants (SNRs), two related classes of astrophysical sources capable of producing γ rays of

sufficient energies to set tight constraints on αn. Chapter 3 discusses the HAWC observatory,

a detector with unprecedented sensitivity to γ rays of such energies, and Chapters 4 and 5

detail the algorithms by which γ-ray event data are analyzed and properties of the particles

are reconstructed. Chapter 6 describes a spectral analysis of the Crab Nebula, the brightest

source of TeV γ rays, and Chapter 7 shows that the highest-energy measurements of photons

7

derived from the Crab imply a tighter constraint on photon-sector LIV than has been set

by any past experiment. Chapter 8 contains some closing remarks about the results and

possible future extensions of the analysis.

8

Chapter 2

γ-ray production at PWNe and SNRs

Pulsars are a class of fast rotating neutron stars. Some pulsars accelerate electrons and

positrons to form PWNe, nebulae that surround them and can in turn radiate TeV γ rays.

PWNe are some of the most common sources of TeV γ rays, and one PWN, the Crab Nebula,

is the brightest such source and will be the subject of the analysis carried out in Chapter 6

and one of the subjects of the analysis discussed in Chapter 7. Pulsars are formed during

supernova explosions, and these explosions often produce SNRs, which consist of supernova

ejecta and interstellar matter accumulated by these ejecta as they propagate.

2.1 Pulsar properties

The neutron was discovered by James Chadwick in 1932 [11]. Lev Landau predicted based

on the existence of the neutron that there would exist very dense stars composed primarily

of neutrons. Walter Baade and Fritz Zwicky published a similar prediction in 1934 [12].

Pulsars were later discovered in 1967 via their radio emission by Antony Hewish and Jocelyn

Bell [13].

Pulsars are products of supernova explosions. They are minute in size compared to their

progenitors. The radius of a neutron-star progenitor is greater than 106 km whereas that of

a pulsar is approximately 10 km. This shrinking results in a pulsar with a core density in
excess of 2× 1014 g/cm3, making pulsars approximately as dense as atomic nuclei. The high

9

density allows for production of neutrons and neutrinos via the process

p + e− → n + νe.

(2.1)

The resulting neutrons are stable in the pulsar.

If they β-decayed, they would produce

electrons with energy no greater than 0.78 MeV, but the Fermi energy of the surrounding

electron gas in the pulsar is hundreds of MeV, and so the β-decay process is disallowed by

the Pauli exclusion principle because there are no electron states unoccupied with a suitable

energy [14].

The high rotational frequencies of pulsars can be understood in terms of angular-momentum

conservation. Prior to the supernova, the star’s radius is around six orders of magnitude

larger than it is after collapsing into a pulsar, and so, assuming that the star possesses some

nonzero angular momentum in its initial state, the frequency of rotation will be magnified

by orders of magnitude by the collapse to compensate for the shrinking radius. Since the
angular momentum is L = Iω ∝ mR2ω (where I is the star’s moment of inertia, ω is its ro-

tational frequency, m is its mass, and R is its radius), the pulsar frequency is approximately

ωpulsar =

star

R2
R2

pulsar

ωstar,

corresponding to a period of

Tpulsar =

R2
pulsar
R2

star

Tstar,

(2.2)

(2.3)

assuming that a negligible fraction of the stellar mass and angular momentum are lost during

collapse. A representative pulsar period of order 1 ms is then obtained by assuming a stellar

radius of 106 km, a pulsar radius of 20 km, and a stellar period of 1 month [14].

10

An order-of-magnitude estimate of the magnitude of the pulsar’s magnetic fields can be

obtained in a similar fashion. It is reasonable to assume that the magnetic flux through the

surface of the northern or southern hemisphere of the star is approximately conserved as the

star collapses into a pulsar:

(cid:90)

pulsar

(cid:90)

star

Bstar · dA.

Bpulsar · dA =

(2.4)

Then the magnetic field scales like

Bpulsar =

star

R2
R2

pulsar

Bstar.

(2.5)

For a representative stellar magnetic field of 1000 G, a pulsar magnetic field of 2 × 1012 G is

computed.

If the axis of symmetry of the pulsar’s magnetic field differs from its axis of rotation,

the magnetic field will be time-dependent as a result of the rotation. This time-dependent

magnetic field then begets an electric field. The speed with which the magnetic field rotates

at the pulsar surface is

v =

2πR

T

,

(2.6)

which gives a value of 4 × 106 m/s for a pulsar with radius 20 km and period 30 ms. The

magnetic field at the pulsar surface should be approximately perpendicular to the velocity

of rotation, so the electric-field magnitude is simply

E = vB,

(2.7)

11

which produces a representative value of 1015 V/m. An electron or positron subject to such

a field would gain energy at a rate of 1 PeV/m [14].

A pulsar’s rotational kinetic energy is

E =

Iω2,

1
2

(2.8)

which gives a value of 4×1061 eV for a pulsar with period 30 ms, mass 2×1030 kg, and radius

20 km. A pulsar that converts its kinetic energy into the energy of accelerated particles with

an efficiency of 1% over a lifetime of 1010 years must convert that energy at a time-averaged

rate of 1042 eV/s [14].

2.2 Pulsar evolution

Since the pulsar’s rotational kinetic energy is

E =

Iω2 =

1
2

2π2I
T 2 ,

(2.9)

its spin-down luminosity, defined as

˙E ≡ −dE/dt, is related to its period and the time

derivative of its period by

˙E =

4π2I ˙T

T 3

.

(2.10)

This quantity represents the rate at which the pulsar’s rotational energy is transformed into

other forms of energy. Measured spin-down luminosities vary by several orders of magnitude,
from 3×1028 erg/s for PSR J2144-3933 to around 5×1038 erg/s in the case of the Crab pulsar
and PSR J0537-6910 [15]. Pulsars with spin-down luminosities exceeding about 4×1036 erg/s

12

normally produce perceptible PWNe [16, 17].

If the assumption is made that the time derivative of the pulsar’s frequency of rotation

has a power-law relationship with the frequency itself, i.e.

˙ω = −kωn,

(2.11)

where n is called the braking index, then solving this differential equation for ω(τ ) (where τ

is the pulsar age) and inverting the relationship to find the age as a function of the period

gives

τ =

T

(n − 1) ˙T

(cid:34)

1 −

(cid:19)n−1(cid:35)

(cid:18) T0

T

(2.12)

for n (cid:54)= 1 [18]. Observations are consistent with the statement that the braking index assumes

a value between 2 and 3 [19]. A value of 3 would be expected assuming the rotational energy

is converted entirely into magnetic-dipole radiation. For the case where n = 3 and the pulsar
has spun down enough that T (cid:29) T0, Equation (2.12) simplifies to

τc =

T
2 ˙T

(2.13)

˙T , and τ have shown
where τc is called the pulsar’s characteristic age. Measurements of T ,
that the true age often exceeds the characteristic age, indicating that the assumption T (cid:29) T0

is violated for those pulsars [20].

For a pulsar producing a purely dipole magnetic field (i.e. neglecting the contribution

from accelerated particles), the constant of proportionality k in Equation (2.11) is related

to the dipole moment by

k =

2M 2⊥
3Ic3

13

(2.14)

where M⊥ is the component of the magnetic dipole moment perpendicular to the pulsar’s

axis of rotation. This implies that the magnitude of the magnetic field at the pulsar’s equator

is

(cid:16)

3.2 × 1019 G/s1/2(cid:17)(cid:112)

T ˙T .

B =

(2.15)

Millisecond pulsars have magnetic fields of order 108 G while the fields of slower rotating pul-

sars called magnetars can exceed 1015 G. For pulsars with perceptible PWNe, the magnetic
field estimated via Equation (2.15) is usually in the range 1 × 1012 G – 5 × 1013 G [17].

If n is time-independent, the pulsar’s spin-down luminosity evolves in time according to

(cid:18)

(cid:19)−(n+1)/(n−1)

˙E = ˙E0

1 +

t
τ0

(2.16)

where ˙E0 and ˙T0 are the spin-down luminosity and period time derivative at the pulsar’s

birth and

τ0 ≡

T0

(n − 1) ˙T0

.

(2.17)

The spin-down luminosity is approximately constant for t (cid:28) τ0; for t (cid:29) τ0, it decays
according to a power law with index (n + 1) / (n − 1).

2.3 Mechanisms of hadronic acceleration

A charged particle moving in the nonuniform magnetic field of the pulsar experiences a

time-dependent magnetic field. The particle’s motion in the direction perpendicular to the

field is approximately circular. According to Faraday’s equation, the time dependence of the

14

magnetic field induces an electric field perceived by the particle:

∇ × E = −1
c

∂B
∂t

.

(2.18)

This electric field varies the perpendicular component v⊥ of the particle’s velocity. Specifi-
cally E⊥ ≡ mv2⊥/2 changes by
(cid:90)

(cid:90)

(2.19)

∆E⊥ =

qE · dl = q

(∇ × E) · dS = −q
c

(cid:90) ∂B

∂t

· dS.

The time derivative of the magnetic field can be approximated in terms of the Larmor period

TL:

∆E⊥ =

q
c

∆B
TL

πR2
L.

Since TL = 2π/ωL, RL = v⊥/ωL, and ωL = qB/(mc),

∆E⊥ =

1
2

mv2⊥

∆B
B

= E⊥ ∆B
B

.

(2.20)

(2.21)

The fractional change in E⊥ is equal to the fractional change in B; in other words, E⊥ ∝ B.

In the observer reference frame, no electric field is present, only a magnetic field, which does

no work on the particle. Since the particle’s total energy is conserved, when it moves into a

region of greater magnetic field and E⊥ increases, E(cid:107) decreases so that E⊥ + E(cid:107) is constant,

and vice versa. As a result, a particle encountering an increasing magnetic field will slow

down and eventually reverse direction. Such a region of heightened magnetic field is thus

known as a magnetic mirror [6].

A particle reflected by a magnetic mirror that is approaching the particle will be reflected

15

with a speed greater than its initial speed (because the initial and final speeds are equal in

the rest frame of the mirror). In the following derivation, the particle’s energy is E, the

component of its momentum in the direction towards (or away from) the magnetic mirror is

p(cid:107), the magnetic mirror approaches with speed U , unprimed and primed quantities refer to

the observer and mirror rest frames respectively, and pre- and post-reflection quantities are

indicated with i and f respectively. The transformed initial energy and momentum are given

by

E(cid:48)
i = Γ
p(cid:48)
(cid:107)i = Γ

(cid:16)
(cid:18)

(cid:17)
(cid:19)

;

Ei + U p(cid:107)i

p(cid:107)i +

U
c2 Ei

.

In the mirror rest frame, the reflection simply flips the sign of the momentum:

E(cid:48)
f = E(cid:48)
i;
p(cid:48)
(cid:107)f = −p(cid:48)
(cid:107)i,

and so the final energy and momentum in that frame are given by

(cid:16)
E(cid:48)
f = Γ
p(cid:48)
(cid:107)f = −Γ

(cid:17)

(cid:18)

;

Ei + U p(cid:107)i
U
c2 Ei

p(cid:107)i +

(cid:19)

.

(2.22)

(2.23)

(2.24)

(2.25)

(2.26)

(2.27)

Performing the inverse Lorentz transformation back into the observer frame gives the final

16

particle energy:

Ef = Γ2

(cid:20)(cid:18)

1 +

(cid:19)

U 2
c2

(cid:21)

(cid:18)

= Γ2Ei

Ei + 2U p(cid:107)i

1 +

2U v cos θ

c2

+

U 2
c2

(cid:19)

(2.28)

where v is the particle’s starting speed and θ is the angle between its initial velocity and the

direction toward the mirror [6].

Taylor-expanding to first order in U/c, which will be small in astrophysical environments,

produces

Ef ≈ Ei

(cid:18)

1 +

2U v cos θ

c2

.

(2.29)

(cid:19)

(cid:18)

(cid:19)

4U
3c

When Equation (2.29) is averaged with respect to θ from 0 to π/2, the result is

Ef = Ei

1 +

,

(2.30)

which gives the factor by which the particle’s energy increases when it interacts with the mir-

ror. If the particle traverses a distance d between two approaching mirrors and is repeatedly

reflected between them, its energy grows in time according to the equation

dE
dt

=

∆E
∆t

=

4U E

3c
d
c

=

4U E
3d

(2.31)

assuming v ≈ c. So the energy increases exponentially with time with a time constant of

τF ≡ 3d
4U

.

(2.32)

This acceleration process is known as first-order Fermi acceleration, in reference to the fact

that τF is proportional to the first power of 1/U . There is also an analogous process in which

17

particles scatter off of a single moving region of increased magnetic field (presumably a gas

cloud); this is called second-order Fermi acceleration due to the proportionality of its time

scale to 1/U 2 [6].

The CR diffusion equation, in the absence of sources, assuming a steady-state solution,

and with an energy-gain rate of E/τF, reads

(cid:21)

(cid:20) E

τF

d
dE

n(E)

+

n(E)
τesc

= 0

(2.33)

where n is the number density of particles per unit energy and τesc is the timescale on which

particles escape from the system in which they are being accelerated. Performing the energy

derivative and rearranging gives

dn(E)

dE

= −

which admits solutions of the form

(cid:18)

1 +

(cid:19) n(E)

E

,

τF
τesc

n(E) ∝ E−(1+τF/τesc).

(2.34)

(2.35)

Fermi acceleration produces a power-law energy spectrum [6].

Because τF ∝ 1/d, first-order Fermi acceleration is efficient in regions in which magnetic-

field inhomogeneities occur on short distance scales. This condition is satisfied in astro-

physical shock waves, and first-order Fermi acceleration occurring at such sites is known as

diffusive shock acceleration (DSA) [6]. The DSA hypothesis is relevant to SNRs because a

PWN’s wind expands supersonically and drives a shock wave into the surrounding SNR [17].

As the SNR initially expands after the supernova, it accumulates mass from the inter-

18

stellar and circumstellar media and eventually enters the Sedov-Taylor phase, in which the

energy of the SNR is approximately conserved. The outer edge of the SNR is a shock, called

the forward shock, that compresses and heats these media; there is also a reverse shock at

the inner boundary of the SNR, which initially moves outward but eventually begins moving

back in towards the pulsar and PWN. At a time several thousand years after the supernova,

still prior to the beginning of the Sedov phase, the SNR reverse shock and the PWN forward

shock meet. Until this happens, the two shocks constitute a pair of converging magnetic

mirrors capable of accelerating the hadronic CRs [17].

In the presence of molecular clouds, the accelerated CRs can interact with hadronic

matter to produce very-high-energy γ rays that are detectable by instruments such as HAWC.

2.4 Models of leptonic acceleration

In the synchrotron self-Compton (SSC) model of γ-ray emission at sites of leptonic accelera-

tion such as PWNe, synchrotron radiation produced by accelerating electrons interacts with

the same electron population again via inverse Compton (IC) scattering to produce γ rays.

The synchrotron photons have energies in the X-ray band. The parent electrons lose energy

to synchrotron radiation at a rate of

−dEe
dt

=

σTE2
6πm2

e B2
ec3

(2.36)

19

where σT is the Thompson cross section. Assuming a power-law electron energy spectrum

with index αe, the spectral energy distribution (SED) of synchrotron photons is then

E2
γ

dNγ
dEγ

= −dEe
dt

dNe
dEe

dEe
dνγ

(2.37)

where N represents a particle flux (particles per unit area per unit time) and νγ is the linear
frequency of a photon. Using Equation (2.36), dNe/dEe ∝ E

, and the relation

−αe
e

(cid:114)2πmecνγ

eB

Ee = mec2

(2.38)

between an electron’s energy and the energy of a resultant synchrotron photon, the SED is

found to be

[6].

E2
γ

dNγ
dEγ

∝ E

−(αe−1)/2
γ

(2.39)

These synchrotron photons are then subject to IC scattering with the same PWN elec-

trons. The energy density of the synchrotron photons in the rest frame of an electron with

Lorentz factor Γ is

u(cid:48) = n(cid:48)E(cid:48)

γ = ΓnΓEγ = Γ2u

(2.40)

where u is the photon energy density, n is the photon number density, and primes indicate

the electron rest frame. One factor of Γ comes from the effect of length contraction on

the number density; the other comes from the Lorentz transformation of the photon’s 4-

momentum. The rate at which an electron loses energy to IC scattering is then

−dEe
dt

=

4
3

σTcu(cid:48) =

4
3

σTcΓ2u.

20

(2.41)

The factor of 4/3 comes from averaging over the angle of collision of the particles.
In
the ultrarelativistic Klein-Nishina regime, where ΓEγ (cid:29) mec2, the average energy of the

produced IC photon is simply determined by the electron energy according to

EIC

γ =

1
2

Ee.

(2.42)

Combining this result with the per electron energy-loss rate (Equation (2.41)) via a derivation

similar to the one leading to Equation (2.39) yields

E2
γ

dNγ
dEγ

∝ E−αe

γ

ln Eγ.

(2.43)

Notably, in the Klein-Nishina regime the photon spectral index approaches the electron index

[6].

IC scattering can also occur via other varieties of target photons. These include the

photons that comprise the cosmic microwave and infrared backgrounds.

The mechanisms of particle acceleration and γ-ray production described in this chapter

can yield the very-high-energy photons needed to test for the signature of Lorentz-invariance

violation. The following three chapters discuss the HAWC observatory and the algorithms

by means of which it reconstructs the properties of these photons.

21

Chapter 3

The HAWC observatory

The HAWC observatory is located at 4100 m above sea level on the side of the Sierra Negra

volcano in the state of Puebla, Mexico. The observatory detects cosmic γ rays, with the goal

of investigating their sources and the mechanisms by which they acquire such high energies.

3.1 Atmospheric air showers

When γ rays and charged CRs strike the Earth’s atmosphere, they produce showers of

lower-energy particles known as extensive air showers. The HAWC observatory detects some

of the shower particles that reach the ground and uses information about their locations,

energies, and arrival times to reconstruct information about the initial cosmic particle, which

is referred to as the primary particle. This reconstruction process is described in detail in

Chapters 4 and 5.

The air showers produced by γ-ray, electron, or positron primaries consist almost exclu-

sively of photons, electrons, and positrons and are called electromagnetic showers. The five

varieties of particle interaction most relevant to the propagation of these showers are pair

production, brehmsstrahlung, ionization, Cherenkov radiation, and trident (three-lepton)

production. Pair production occurs when a photon strikes the nucleus of an atom of

an atmospheric molecule such as nitrogen and is converted to an electron/positron pair.

Brehmsstrahlung occurs when an electron or positron is slowed down by the electromagnetic

22

field of the nucleus of an atmospheric molecule and emits a photon due to this acceleration.

Ionization occurs when an electron or positron strikes a molecule and knocks off an electron.

Cherenkov radiation occurs when an electron or positron moving faster than the speed of

light in the medium in which the shower is propagating radiates photons by analogy with

a sonic boom. Finally, trident production occurs when an electromagnetic interaction be-

tween a shower electron or positron and the nucleus of an atmospheric molecule causes the

production of an electron/positron pair. Above 80 MeV (an energy known as the critical

energy), brehmsstrahlung is the dominant interaction of shower electrons and positrons with

the atmosphere, such that the number of particles participating in the shower grows more

or less exponentially since energy is approximately conserved within the shower. Below the

critical energy, the electrons’ and positrons’ energy is quickly dissipated into the atmosphere

via ionization [21].

An approximate form for the longitudinal evolution of the total number of particles in

the shower is given by Greisen [22]:

N (t) ≈

(cid:112)ln(E0/Ec)

0.31

(cid:26)

exp

t

(cid:20)

1 − 3
2

ln

(cid:21)(cid:27)

.

3t

t + 2 ln(E0/Ec)

(3.1)

Here t is the shower depth (defined as the distance the shower has propagated in units of

37.1 g/cm2, the radiation length of electrons and positrons in air [21]), E0 is the primary

energy, and Ec is the critical energy. According to this formula, the depth at which the

number of particles is maximized, referred to as shower maximum, is

t = ln(E0/Ec) .

(3.2)

23

In order to determine the lateral width of the shower, it is necessary to consider also

the Coulomb scattering of electrons and positrons off atmospheric molecules. This process

dominates the transverse spreading out of the shower particles as the shower propagates.

The characteristic length scale for this Coulomb-induced spreading is the Moli`ere radius,

defined as

r1 ≡ 2mec2X0

Ec

(cid:114) π

α

= 9.3 g/cm2,

(3.3)

where me is the electron mass, X0 is the radiation length, and α is the fine-structure constant.

At HAWC’s altitiude this corresponds to 124 m. Greisen [22] and Kamata and Nishimura [23]

provide the following expression, called the NKG function, for the lateral shower distribution:

(cid:19)s−2(cid:18)

(cid:18) r

r1

ρ(r) ∝

(cid:19)s−4.5

1 +

r
r1

(3.4)

where s is the shower age, a parameter appearing in the Mellin-transform solution to the

differential equations governing shower evolution. The age has a value of 0 at the start of

the shower and increases as it propagates.

The NKG function is sharply peaked at r = 0. This peak is called the shower core.

The core propagates along a line called the shower axis. The region of space containing the

shower particles is approximately an elliptic paraboloid but is often referred to as the shower

plane. The shower core, axis, and plane are depicted in Figure 3.1.

3.2 Water-Cherenkov detectors

The HAWC observatory is made up of 300 water-filled tanks that utilize the water-Cherenkov

technique to detect the charged particles produced in air showers. These tanks occupy 57% of

24

Figure 3.1: An atmospheric air shower. Dots indicate shower particles. Image taken from
[1] with permission.

25

an area of 22000 m2 of a flat region constructed on the side of the Sierra Negra volcano. Each

steel tank is cylindrical in shape with a diameter of 7.3 m and a height of 5 m and contains

a light-tight plastic bladder filled to a depth of 4.5 m with water. The water is purified

to the point that its attenuation length is around 10 m, such that Cherenkov photons can

propagate unattenuated for distances of order the tank size [24]. The bottom of the tank

is instrumented with three Hamamatsu R5912 photomultiplier tubes (PMTs), referred to

as the A, B, and D PMTs, and one high-quantum-efficiency (high-QE) Hamamatsu R7081

PMT, referred to as the C PMT. The R5912 PMTs are 8 in in diameter, and the R7081 is

10 in [2]. The former were repurposed from the Milagro experiment to reduce the cost of

constructing HAWC. The bladders contained by the tanks are black to minimize reflection

of Cherenkov photons so that the PMTs only detect photons emitted directly by the shower

particles [25].

When a cosmic γ ray interacts with the atmosphere above the detector, it induces an air

shower. Once the shower propagates down to the altitude of the detector, shower electrons

and positrons continue to shower in HAWC’s tanks and also produce Cherenkov light if they
are moving faster than the speed of light in water (2.26 × 108 m/s). Shower photons also

enter the tank and continue to shower into electrons and positrons in the water, which can

produce Cherenkov light as well. The Cherenkov photons are then detected by the PMTs.

This process and the layout of the tanks and PMTS are shown in Figure 3.2. Since the

Cherenkov threshold for positrons and electrons in water is 1 MeV and the ionization energy

loss of positrons and electrons in water is 2 MeV/cm [26], positrons and electrons entering

the water with energies below 900 MeV will drop below their Cherenkov threshold before

they exit the tanks if they pass through the full 4.5 m of water. However, only a small

fraction of this energy loss is via Cherenkov radiation; most is via ionization (as with a

26

Figure 3.2: Left: a HAWC tank with a charged particle passing through (the thick line
coming from above) and the resulting Cherenkov cone (the thinner lines). Right: the layout
of HAWC’s 300 tanks (large circles) and their PMTs (smaller filled circles). The gap in the
center is the location of the counting house, which holds the electronics. Image taken from
[2].

shower propagating in air). Conversely, the median energy of ground-level muons, which

are relatively rare in γ-induced showers but common in hadron-induced showers, is such

that they will not go below their Cherenkov threshold of 213 MeV before exiting the tanks.

Therefore such a muon will usually pass closer to one PMT than to the rest, producing an

excess signal in one PMT, which can be used to discriminate hadronic showers from photonic

ones [24]. Such algorithms are described in more detail in Section 4.5.

The 10-in PMTs located at the centers of HAWC’s tanks have a quantum efficiency

(QE, the probability of producing a photoelectron (PE) that reaches the first dynode of the

PMT) about twice that of the 8-in PMTs. This fact is taken into account by the various

27

reconstruction algorithms described in Chapter 4. The PMTs are provided with a voltage

of approximately 1700 V, but the exact values are chosen per PMT in order to match their
gains, which are approximately 1.6 × 107 [24]. The 8-in and 10-in PMTs both contain 10

dynodes each.

If a photon shower with primary energy 1 TeV originates from directly above HAWC,

on average 7% of that primary energy will reach the ground in the form of the energy of

the shower particles. For 100-TeV photons, this number is 28%. Overhead γ showers with

primary energy above 1 TeV are detected with nearly 100% efficiency, where a shower is

considered to trigger the detector if it produces hits in at least 6.7% of the active PMTs in

the event. Those of lower energies can also be detected if their first interactions occur at

lower than average altitudes [2].

3.3 Data-acquisition system

Each PMT is connected by an RG-59 cable to one of the front-end boards (FEBs), which

prepare the signals to be fed into time-to-digital converts (TDCs). The cables are all of

the same length, 610 ft, so that signals from all PMTs require the same time to propagate

to the counting house, the building located in the center of the HAWC array containing

the data-acquisition system (DAQ). Each FEB is actually two boards connected to a Versa

Module Europa (VME) crate (Wiener 6023x610) by which they can interact. A FEB has

16 input channels.

It amplifies and integrates each input signal. The signal is then sent

to two comparators with two different voltage thresholds. The comparators’ outputs are

multiplexed by the digital half of the FEB, producing an ECL signal with two or four

transitions depending on whether the voltage ever reached the higher of the two thresholds.

28

Figure 3.3: The output voltage from the PMT (top) and the resulting FEB output. From
this output the time over threshold (ToT) for each threshold can be extracted and used to
estimate the number of PEs produced in the PMT. Image taken with permission from [1].

See Figure 3.3. The low threshold corresponds to a charge of about 1/4 PE deposited in the

PMT; the high threshold corresponds to about 4 PEs. The edge times contain information

sufficient to reconstruct the time and deposited charge of the original PMT pulse. The

process by which this information is extracted is described in Section 4.1. The analog

portion of the FEB also distributes the high voltage (HV) required by the PMTs; each FEB

provides HV to the 16 PMTs from which it receives data [3].

The times are measured by a CAEN V1190A-2eSST TDC, which can handle up to 128

PMTs’ worth of FEB signals. The complete detector uses 10 TDC units. The TDCs receive

an external 40-MHz clock signal, with 8 output bits used to represent 256 subdivisions of

the 25-ns clock period, for a timing resolution of 98 ps. TDC output is triggered every 25 µs;

at each trigger, each TDC writes out 26 µs of edge-time data (the overlap avoids dropping

29

VoltageToTt0t1t2t3timelow thresholdhigh thresholdlow ToThigh ToTedges on the time-window boundary) along with other peripheral data such as time stamps

and error information. The resulting data rate for each TDC is 45 MB/s, for a per PMT

data rate of 350 kB/s [3].

Each TDC module is controlled by a General Electric XVB602-13240010 single-board

computer (SBC). The SBC runs a custom program called the Readout Process on the CentOS

operating system that receives the TDC output and sends it via Gigabit Ethernet using

ZeroMQ [27] to the online-reconstruction farm. The readout rate of TDC data is about

60 MB/s, limited by the clock speed of the VME chip on the TDC. When a TDC has taken

25 TDC events’ worth of data, corresponding to about 625 µs, it sets a data-ready bit, which

is read by the SBC, initiating the data transfer to the online-reconstruction clients. The

SBC recombines 42 such blocks of data, the last two of which must be redundant with the

following set of 42 due to occasional TDC firmware errors, and pushes them to the Data

Queue, which feeds data to the Reconstruction Clients [3].

The Reconstruction Client is a program that runs on several dedicated computers, col-

lectively known as the online-reconstruction farm, and reconstructs various parameters char-

acterizing the air showers detected by HAWC. See Chapter 4 for the details of the recon-

struction process and its outputs. A Reconstruction Client polls the Data Queue for a block

of 25 ms of data and runs the reconstruction algorithms on them [3].

The various processes running on the DAQ and online-reconstruction farm are coordi-

nated by a system called Experiment Control. Experiment Control monitors and adjusts

several settings of the DAQ and starts and stops data-taking runs. It also reports and logs

significant errors in the operation of the detector and will attempt to restart the experiment

when a failure causes it to crash.

In the modern configuration of the experiment, these

attempts are almost always successful although crashes and subsequent restarts often occur

30

frequently during lightning storms. The restart process takes several minutes. Experiment

Control communicates with the processes that it oversees by serializing data to the JSON

format via the cJSON library [28]. Experiment Control is also responsible for sending alerts

to external systems which provide alerts regarding possible astrophysical transients such as

γ-ray bursts (GRBs), among them the GRB Coordinates Network (GCN) [3].

Monitoring of the experiment is performed by a dedicated Python application. This

package obtains data directly from some monitored components and indirectly from others

via Experiment Control. These data are stored via MySQL and copied to servers at the

University of Maryland (UMd) and the National Autonomous University of Mexico (UNAM).

A web page provides plots generated via the Google Charts API to shifters, the collaboration

members assigned to monitor and control the operations of the experiment [3].

Figure 3.4 shows the flow of data among all components of the DAQ and online-reconstruction

system.

3.4 Science goals

The HAWC observatory was constructed with the goals of studying TeV γ-ray emission from

Galactic and extra-Galactic astrophysical sources as well as various possible new particle-

physics processes occurring in the Universe. The detector is also sensitive to charged CRs,

including their energy spectrum and anisotropy.

PWNe, described in Chapter 2, are the most common TeV sources in the Galaxy [17].

These include the Crab Nebula, the brightest source of TeV γ rays. The AGILE and Fermi-

LAT experiments have observed flares from the Crab between 100 MeV and 100 GeV [29, 30].

The ARGO-YBJ experiment has found evidence of possible TeV time variability [31], but

31

Figure 3.4: Diagram of the HAWC data-acquisition system [3].

no imaging air-Cherenkov telescopes (IACTS) have seen indication of Crab flare events in

the TeV. According to the Fermi-LAT analysis, the observed Crab flare represented, at its

peak, a 50-fold increase in the synchrotron emission from the Crab. The suggestion has been

made, in order to explain the relatively short time scale on which such large flares can occur,

that the underlying particle-acceleration mechanism is magnetic reconnection and not shock

acceleration [32]. Among the motivations of HAWC’s construction was its sensitivity to a

possible TeV inverse-Compton component to these flares: HAWC would be able to detect a

10-fold Crab flare in 10 min. The experiment also monitors other PWNe for flares [33]. Crab

time dependence on the time scale of its rotational period has also been observed [34], but

HAWC is unlikely to be sensitive to this.

SNRs, also described in Chapter 2, are another class of sources of great interest in the

study of TeV γ rays. They are a leading candidate for the source of Galactic CRs. The

32

observation of TeV γ rays from SNRs could offer evidence of hadronic acceleration at such

sources. This is because an SNR accelerating hadrons near a molecular cloud would cause

the production of γ rays via a reaction such as

2p → 2p + π0 → 2p + 2γ

(3.5)

The spectral shape of such γ-emission would be distinct from that of inverse-Compton emis-

sion.

IC443 is an interesting SNR from which TeV emission has been observed by the

VERITAS [35] and Fermi telescopes [36]. Specifically Fermi claims to have detected the

distinct shape of π0-decay in the IC443 energy spectrum. The possibility of observing the

source with HAWC at higher energies was a particularly strong motivation for HAWC’s con-

struction because of Milagro’s 3σ detection of IC443 at 35 TeV [37], suggesting the spectrum

was harder than measured by VERITAS [33].

Another motivation for the construction of HAWC was the prospect of observing active

galactic nuclei (AGNs). AGNs are prone to large flares, sometimes up to a factor of 10 in

flux. Because AGNs are extra-Galactic sources, γ rays detected from them must propagate

over large distances, and so they can also be utilized for the study of axions [38]. HAWC’s

status as an all-sky TeV detector is of particular significance because it implies sensitivity to

“orphan” flares, inverse-Compton AGN flares which occur in the absence of contemporaneous

synchrotron flares [39]. These suggest hadronic acceleration [40, 41]. As the only current

all-sky TeV γ-ray observatory, HAWC has the unique ability to detect these flares despite

the lack of lower-energy counterparts on which to trigger [33].

HAWC is also sensitive to dark-matter-induced TeV γ-ray emission. Objects with high

mass-to-light ratios are the primary targets for such an analysis since relatively little γ-

33

emission induced by CR acceleration is expected from these sources [42, 43, 44]. The peak γ-

production rate from dark-matter decay or annihilation typically occurs at an energy roughly

an order of magnitude below the dark-matter mass, so HAWC’s dark-matter sensitivity is

competitive for a mass above 10 TeV [33].

Lastly, HAWC was constructed with the intent of searching for GRBs in the TeV. The

Fermi Large-Area Telescope (LAT) observed a GRB emitting a photon with energy 95 GeV

[45]. Milagro’s predecessor Milagrito found 3σ evidence for TeV γ rays from GRB 970417A

[46]. Milagro continued to search for TeV emission from GRBs but did not find any [47].

HAWC was designed to have sufficient TeV GRB sensitivity to determine whether such

sources radiate in the TeV. The experiment’s all-sky nature implies that triggered searches

can be performed retroactively simply by looking back at existing HAWC data from the

time and location of a known burst seen at lower energies. Determining whether the GRB

spectrum cuts off at TeV energies has implications regarding the propagation of γ rays

through the interstellar medium. A cutoff is expected based on models of the interaction

of γ rays with the Extragalactic background light (EBL) [48], and so a measurement of a

higher than predicted cutoff (or none) would suggest physics beyond the SM, such as axions

[49].

34

Chapter 4

Air-shower reconstruction

HAWC uses dedicated software to perform event reconstruction, the process by which raw

data are analyzed and quantities such as shower-core location, shower direction, primary-

particle species, and primary energy are extracted. These shower features can then be used

to do high-level data analysis of astrophysical sources. The event-reconstruction software is

part of HAWC’s Analysis and Event-Reconstruction Integrated Environment (AERIE) C++

framework.

4.1 Calibration

Two varieties of calibration must be performed on HAWC data in order for it to be usable

for event reconstruction. These are the timing calibration, which corrects for PMT-to-PMT

timing variations and for the slewing time needed for a PMT’s voltage to reach the low

threshold, and the charge calibration, which produces an estimate of the number of PEs

produced in a PMT from the time its output voltage spends over each of the two thresholds.

The timing calibration begins with the taking of laser calibration data. A laser is fired

into each HAWC tank through each of 68 neutral density filters with transmittances ranging

from 22% to 87%. The resulting data are stored in a ROOT [50] file. A function of the form

35

S(T ) = exp

−T − p0

p1

− exp

T − p2

p3

+ p4 − p5T

(4.1)

is fit to the data, where T is the time over threshold (ToT), S is the average slewing time,

and the quantities pi are real numbers that float in the fit. Figure 4.1 shows the slewing

time as a function of ToT in the laser calibration data. The relationship is approximately

linear. The upturn at small ToT, modeled by the decaying exponential in Equation (4.1),

occurs for pulses that barely exceed the threshold because of the curvature near the top of

these pulses. The growing exponential accomodates the very small downward curvature at

large ToT.

Several PMT-specific timing corrections are applied. Since the diffuser that emits laser

light into the tank is located at the center of the tank’s roof over the C PMT, the laser

light reaches the C PMT earlier than the A, B, and D PMTs by on average 2.2 ns, so the

slewing times for C PMTs are corrected by this amount. A per tank correction must also

be made for the lengths of the optical cables that lead from the laser in the counting house

to the tanks. These corrections can be as large as 10 ns. Finally a per PMT time-residual
correction is applied. This is determined by iteratively reconstructing batches of 2 × 105

events and computing a correction for each PMT equal to its time residual averaged over the

sample of events. Three iterations of this procedure are performed. The result is an additive

correction to p4 for each PMT.

Charge calibration is also carried out using the laser data. The laser is fired several times

on each filter-wheel setting, and the fraction of laser pulses resulting in at least one hit in

the PMT is computed. A PMT is considered to have triggered on a laser pulse if it registers

a hit within 2 µs of the pulse. The fraction of pulses causing a PMT to trigger is an estimate

36

Figure 4.1: Slewing times shown as functions of low and high ToTs. The curves show the fit
function from Equation (4.1).

37

of the probability that laser light and/or a CR causes a hit in that PMT during the 2-µs

time window. This can then be related to the probability of a laser trigger:

Ptrig = Plaser + PCR − PlaserPCR

(4.2)

where Plaser is the occupancy (i.e. the probability of triggering on the laser by production

of a PE passing the low threshold), PCR is the probability of triggering on a CR, and Ptrig

is the probability of triggering on either. PCR can be measured from data taken while the

laser is not firing. Assuming that the number of PEs produced in the PMT is a Poisson

variate, the Poisson mean for a given laser intensity can be computed using

Plaser = P (n > 0) = 1 − P (n = 0) = 1 − e−µµ0

0!

= 1 − e−µ

(4.3)

where n is the number of PEs produced and µ is the mean. So, combining Equations (4.2)

and (4.3), the measured mean number of PEs from a laser pulse of a given intensity is

computed as

[51].

µ = ln

1 − PCR
1 − Ptrig

(4.4)

This technique is only useful when the occupancy Plaser is signifantly less than 1; for an

arbitrarily high laser intensity, the occupancy saturates at 1, meaning that the occupancy

no longer contains useful information about the laser intensity. However, since µ should be

directly proportional to the laser intensity, this relationship is fit for occupancies between

0.05 and 0.90 and then extrapolated [51].

For each filter setting applied to the laser, a simple Monte Carlo (MC) simulation is

38

performed in which a Poisson distribution with the appropriate mean µ for that laser con-

figuration is sampled and Gaussian noise is added. This noise represents the PMTs’ charge

resolution. The resulting PE-count distribution is compared with the ToT distribution from

the laser data, and nine ordered pairs of PE count and ToT are generated by matching the

tenth through nintieth (in increments of 10%) percentiles of the two distributions. These

pairs are aggregated across different laser intensities to produce a curve relating ToT with

PE count [51].

4.2 Edge finding

As discussed in Section 3.3, larger (“four-edge”) hits will produce four TDC edges as the

voltage exceeds first the low threshold and then the high threshold and then falls below the

high threshold and then the low threshold. Smaller (“two-edge”) hits never exceed the high

threshold. There is the potential for ambiguity between a single four-edge hit and a pair of

two-edge hits because in each case four edges are present. These cases can be distinguished

using the time between the first two edges. The FEB forces the low ToT for any hit to be

at least 53 ns by delaying the trailing low-threshold edge if necessary. But the time between

the rising low- and high-threshold edges is normally much smaller than this. So four-edge

hits can be distinguished from two-edge hits based on whether the time between the first

two edges is greater than 53 ns [51].

The edge finder also identifies and flags hits that need to be removed prior to recon-

struction. Two very small pulses can create a spurious two-edge hit, which the edge finder

marks with a “bad” flag. PMT afterpulsing can also produce unusable data. An afterpulse

that occurs during a four-edge hit and is strong enough to cause the voltage to rise back

39

above the high threshold after having already fallen below it results in a “six-edge” hit,

which is marked using an “ambiguous” flag. Bad and ambiguous hits are not used by the

reconstruction algorithms described below [51].

4.3 Core reconstruction

The location of the shower core is used by several other reconstruction algorithms, including

the angle fitter, background-rejection variables, and energy estimators. In order to fit the

location of the shower core, a model of the lateral distribution of particles within the shower

must be assumed.

In principle the correct model is the NKG function (Equation (3.4)).

However, properly normalizing the NKG distribution requires computing Γ functions, which

are relatively expensive to evaluate and also require numerical differentiation [2].

Instead a more practical functional form has been developed:

(cid:18)

(cid:19)

ρ(r) ∝ 1

2πσ2 exp

− r2
2σ2

+

N

(0.5 + r/r1)3

(4.5)

where σ ≡ 10 m and N ≡ 5 × 10−5. This function is a linear combination of a two-

dimensional Gaussian and a power-law-like function meant to approximate the behavior of

the NKG function. However, there is no free shower-age parameter; only the normalization

and the x- and y-coordinates of the core float in the fit. Equation (4.5) is also analytically

differentiable and does not diverge at the core like the NKG function [2].

This function is fit to the charge observed in each PMT in order to estimate the location

of the shower core. Figure 4.2 shows the effective charge in each PMT for an example event.

The core fitting algorithm is known as the Super Fast Core Fit (SFCF). The median error of

40

Figure 4.2: The effective charge in each PMT in a γ-ray event. Image taken from [2].

the SFCF core estimate is approximately 4 m for smaller air showers and 2 m for larger ones,

among showers whose cores land on the array. Core reconstruction is less precise for off-array

showers; a shower with a core 50 m from the edge of the detector has a core resolution of

approximately 35 m [2].

Prior to the execution of the SFCF, a simple estimate of the core location is performed.

This estimate is the average of the PMT locations weighted with their estimated charges

(the “center of charge”) and is used as the initial guess for the core location in the SFCF

minimization [2].

41

4.4 Angle reconstruction

As shown in Figure 3.1, shower particles in a shower originating not directly from zenith

arrive at the detector at different times because some must propagate a greater distance than

others. The arrival times of the hits can therefore be used to estimate the direction from

which the primary particle came. These times are shown for an example event in Figure

4.3. The simplest model of the propagating shower consists of a plane of particles. Two

considerations complicate this model: the differences in distances travelled by particles at

different distances from the shower axis and the fact that each PMT measures the time of

the first Cherenkov photon to strike it [2].

Because the shower begins with the interaction of the primary γ ray at one point in the

atmosphere, it propagates as a cone, with particles far from the shower axis therefore falling

behind those closer to it. This curvature effect results in the nonplanar surface shown in

Figure 3.1 [2].

Furthermore, the stochastic nature of the particle interactions between shower particles

and the atmosphere means that even particles at the same distance from the axis travel

different distances since they scatter multiple times and take nonlinear paths to the ground.

The resulting time delay between the earliest and latest particles in the shower is of order

several ns near the axis and greater far from the core. When several shower particles reach

the water of a tank at different times, the resulting Cherenkov light can produce overlapping

hits in the PMTs. The recorded time is that of the first light. The difference in shower thick-

ness between particles near the axis and those far from it necessitates a position-dependent

sampling correction since the problem is more severe far from the axis [2].

These two effects are corrected for together. MC simulations of showers and their in-

42

Figure 4.3: The arrival time of each hit used in reconstructing the direction from which the
example event in Figure 4.2 originated. Image taken from [2].

43

teractions with HAWC were used to deduce a suitable functional form for the correction.

This function was then fit to a sample of actual HAWC data events and verified on another

disjoint population. Fitting the function to data events is necessary because the MC simula-

tion is not sophisticated enough to adequately capture the subtleties of the effects discussed

above. A χ2 fit is then used to fit a plane to the corrected times of the hits [2].

4.5 Background rejection

The primary background to the γ rays detected by HAWC consists of various hadronic CR

species, mostly protons. In order to distinguish between CR events and γ events, two vari-

ables have been developed that exploit the differences in the signatures of these classes of

events in HAWC data. As discussed in Section 3.1, γ-induced air showers consist mostly of

electrons, positrons, and photons. Muons occur uncommonly in these showers. By contrast,

hadronic showers contain muons and mesons such as pions, and the subsequent electro-

magnetic sub-showers produced by hadronic secondaries appear as inhomogeneities in the

shower’s lateral distribution [2].

The background-rejection variable known as compactness is defined by

C ≡

Nhit

maxi∈F qi

(4.6)

where Nhit is the number of hits within 20 ns of the reconstructed shower plane, qi is the

effective charge in the ith PMT, and F is the set of indices of PMTs greater than 40 m from

the reconstructed core of the shower. Because hadronic showers contain muons and sub-

showers, they are likely to have large hits beyond the 40-m radius. This is not the case for

44

(a) A hadron event.

(b) A photon event.

Figure 4.4: The lateral distribution of charge within a background event and a signal event.
Image taken from [2].

γ events, wherein the lateral distribution is fairly smooth. Compactness is therefore larger

on average for γ events than for hadronic ones [2].

The second variable used to discriminate between hadronic and photonic showers, known

as the Parameter for Identifying Nuclear Cosmic rays (PINCness), is defined as

(cid:32)

Nhit(cid:88)

i=1

P ≡ 1
Nhit

(cid:33)2

ζi − (cid:104)ζi(cid:105)

σζi

(4.7)

where ζi ≡ log10 qi, (cid:104)ζi(cid:105) is the average log charge among PMTs with radial distance within

5 m of the ith PMT’s radial distance, and σζi

is the expected fluctuation of the log charge

in the ith PMT as measured on a sample of events near the Crab Nebula. Figure 4.4 shows
the log charges and the local average (cid:104)ζi(cid:105) used in the PINCness computation. In hadronic

showers, such as Figure 4.4a, the charges tend to differ more from their local averages than in

γ showers, such as Figure 4.4b. This is again due to a combination of muons and secondary

sub-showers and results in a higher PINCness value for hadrons than for photons.

45

020406080100120140PMT Dista ce to Reco structed Core [m]−101234log10(Qeff)SFCF FitPINC Movi g Average <ζ>Qeff020406080100120140PMT Dista ce to Reco structed Core [m]−101234log10(Qeff)SFCF FitPINC Movi g Average <ζ>Qeff4.6 MC simulations of air showers

A detailed MC simulation of the HAWC experiment is used to optimize and evaluate the

reconstruction algorithms described above and in Chapter 5. This same simulation is also

used in spectral likelihood analyses such as the ones described in Chapters 6 and 7. The

simulation of a single event consists of propagating the air shower through the atmosphere,

simulating the interactions of the shower particles with the detector on the ground, and

running the reconstruction algorithms on the resulting pseudodata.

The shower propagation is performed using version 7.4000 of the CORSIKA package [52].

GEANT version 4.10.00 [53] is used to simulate the propagation of the shower particles and

Cherenkov photons through the water in the tanks. The simulation of the PMTs is performed

using custom software. This software models each PMT as having an efficiency chosen to

match vertical-muon rates in data. The PMT signal is then fluctuated according to a charge-

smearing parameter, and single-PE noise is added. The width of the charge smearing and the

distribution of single-PE noise are chosen such that distributions of reconstructed quantities

in the simulation match the distributions of the same quantities in data. The simulated

PMT signals are then reconstructed using the same reconstruction software used with data

[2].

Any inaccuracies in the assumptions made in modeling the detector in these simulations

will produce systematic errors in analysis results that depend on them.

It is therefore

necessary to estimate uncertainties on several quantities in the simulation, perform analyses

using several variants of the simulation, and cite the spread in the estimates of the relevant

physics parameters as a systematic error. The quantities in the simulation for which this

procedure is performed are the charge-smearing parameter described above, the QE of the

46

PMTs, the PE threshold in the PMT, and the presence or absence of an effect called “broad

pulse” in the simulation, described below.

The width of the charge smearing is expressed as a fraction of the recorded charge. The

nominal value of this quantity in HAWC’s simulations is 0.1100, with 0.0935 – 0.1265 being

the range used for evaluating systematic errors.

A QE of 55% has been chosen as the best value for the A, B, and D PMTs; 80% is used

as the nominal QE for the high-QE C PMTs. The A/B/D QE is varied between 44% and

66% when studying systematics, and the C QE is varied between 64% and 96%.

The PMT threshold is assigned a nominal value of 0.55 PE. A lower PMT threshold of

0.30 PE has been shown to match data event rates better than the 0.55 PE value but shows

poorer agreement with respect to the distributions of hadron-rejection variables, and so the

threshold is varied to 0.30 PE to study systematics.

Because calibration data are produced using laser pulses, HAWC’s charge calibration

assumes photons arrive at PMTs promptly, with arrival times differing by less than 10 ns.

However, in actual air showers photons do not arrive simultaneously at the PMTs and as

a result produce larger than expected ToTs and reconstructed charges. To compensate for

this, a broad-pulse effect is introduced into the simulation, in which PE estimates in PMTs

are increased in order to resemble the values observed in data. This effect is part of the

nominal simulation; to study systematics, it is removed.

The uncertainty ranges for these MC parameters are felt to be somewhat conservative.

Studies of the HAWC detector and simulation are ongoing and are expected to result in more

precise measurements of these quantities in the near future.

47

Chapter 5

Energy estimation

The Neural-Net (NN) Energy Estimator is a γ-ray energy-estimation algorithm developed

for HAWC with the dual goals of calculating an accurate and precise energy estimate while

still being efficient enough to be run on all HAWC data as part of event reconstruction. It

uses the Toolkit for Multivariate Analysis (TMVA) [54] NN implementation.

The MC-based plots presented throughout this thesis (with the exception of the Ap-

pendix) were made using the HAWC MC simulation internally referred to as DAQSim 4.

HAWC collaborators can find more details about this simulation in HAWC internal memo

2397.

5.1 Method

An NN is a function mapping several quantities associated with an event (input variables) to

some regression target or output variable, in this case log10 E where E is the primary energy

of the shower. A detailed description of the functional form of the NN can be found in [54].

This function is characterized by many (479 in this implementation) free parameters called

weights. The optimal values of the weights are determined in a process called training.

48

Shower Characteristic Described

Multiplicity

Fraction of shower contained in the detector

Atmospheric attenuation of shower

Name
Fraction of channels hit
Fraction of tanks hit
Log core amplitude
Core distance
Cosine zenith angle
Fraction of charge in annulus 0
Fraction of charge in annulus 1
Fraction of charge in annulus 2
Fraction of charge in annulus 3
Fraction of charge in annulus 4
Fraction of charge in annulus 5
Fraction of charge in annulus 6
Fraction of charge in annulus 7
Fraction of charge in annulus 8
Fraction of charge in annulus 9

Table 5.1: NN input variables that characterize energy. The annuli are about the shower
axis. Annuli 0 – 8 are each 10 m thick. “Annulus” 9 extends from 90 m to infinity.

5.1.1 Input Variables

Fifteen input variables are used to characterize the shower energy. These variables have been

chosen in order to capture three qualities of the shower: the multiplicity in the detector,

which is a crude measure of how much energy reaches the ground, the containment of the

shower within the detector, which indicates by how much the multiplicity in the detector

needs to be scaled up, and the atmospheric attenuation of the shower, which indicates how

much energy was lost on the way to the ground. The input variables are summarized in

Table 5.1.

5.1.1.1 Multiplicity

Three variables are used to quantify the overall multiplicity of the shower. The first, the

fraction of channels hit, is defined here as the number of channels hit in the event divided

49

by the total number of functioning channels in the detector. Channels with bad hits are

conventionally excluded from the denominator in this quantity in other contexts in HAWC

analysis. However, in practice both choices have been found to yield comparable estimator

performances, and the calculation is faster when bad channels are not excluded.

The second variable relating to the shower multiplicity is the fraction of tanks hit. A

tank is considered hit if at least one of its PMTs receives a hit.

The third multiplicity variable is the logarithm of the overall normalization parameter

from the SFCF core fit. This quantity should in principle be more robust against off-

array showers than the fraction of channels or tanks hit since, with an ideal fitter, the

fit normalization would be independent of the containment of the shower.

In practice,

this parameter is one of the least important NN inputs, perhaps indicating that the fit

normalization parameter is a poor estimate of the integral of the shower’s lateral distribution.

5.1.1.2 Containment

One input variable is used to infer the fraction of the shower contained within the array: the

distance of the reconstructed core from the array center. This is somewhat crude because

the detector is not rotationally symmetric. The “core fiducial scale” C, defined as 100 times

the factor by which the detector area would have to be scaled in order to barely contain the

reconstructed core, has been evaluated as an alternative measure of the containment of the

shower but does not yield a better energy estimate.

5.1.1.3 Atmospheric Attenuation

Eleven variables are used to account for the attenuation of the shower by the atmosphere on

its way to the ground: the zenith angle (in cosine) and ten annular variables representing

50

the lateral distribution of charge in the shower front. The zenith angle θ is included because

showers at higher zenith must pass through more atmosphere. Training was attempted using

θ, cos θ, or sec θ as an input variable; cos θ proved to yield the best performance.

The annular variables are defined as follows. For i between 0 and 8 inclusively, the ith

annular variable is the summed charge of all hits, excluding those in bad PMTs, landing a

distance between (10 m)i and (10 m)(i + 1) from the shower axis, divided by the total charge

in the event. The distance is measured in the reconstructed shower plane, not the plane

of the detector. The ninth is defined as the fraction of the charge landing more than 90 m

from the shower axis. These ten variables capture the lateral energy distribution within the

shower.

The purpose of the annular variables is to estimate the first-interaction height of the

shower, which is also relevant in determining how much atmosphere the shower has passed

through. One way of estimating the first-interaction height is via the shower-age parameter

of the NKG function, which is higher for “older” showers. However, introducing the age

parameter from an NKG fit as an NN input does not improve the energy reconstruction.

It may be a poor estimate of the age. Since higher-age showers are wider, there should be

some information about the first-interaction height in the lateral distribution of the shower,

and this is why the annular variables are included: as a more empirical way of mapping the

lateral distribution to the first-interaction height.

Using the fraction of charge deposited in an annulus instead of the total charge desen-

sitizes the NN algorithm to systematic errors involving the charge scale as modeled by the

MC. Three of the known systematic uncertainties affecting HAWC’s MC are uncertainties in

the PMT threshold, charge fluctuations, and QE, all of which impact estimates of PE counts

in the MC. The NN is made more robust against these sources of uncertainty by eliminating

51

the dependency on the total PE count in an event in favor of utilizing only the shape of the

lateral distribution.

5.1.2 Training

The goal of the training process is to choose the vector of NN weights w that minimizes the

error function, defined as

D(w) ≡ 1
2

(cid:104)

ui

Nevents(cid:88)

i=1

(cid:105)2

log10

ˆE(xi; w) − log10 Ei

(5.1)

where ui is the weight of the ith event (the event’s relative importance in the training,

unrelated to w, which are also called weights), Ei is the MC-truth energy of the ith event,

xi is the ith event’s vector of input variables, and ˆE is the energy estimate for the event.

This minimization is performed by repeatedly iterating over the sample of training events,

with a small update to w performed at each event. The details of the Broyden-Fletcher-

Goldfarb-Shanno (BFGS) algorithm for updating the weights can be found in [54]. The MC

events are weighted (the ui are set) to resemble an E−2.63 (Crab Nebula-like) power-law

spectrum during the training. This weighting has been observed to produce a relatively flat

RMS error between 1 and 100 TeV (see Section 5.2.6).

The event weighting manipulates the thrown E−2 spectrum into the desired E−2.63

spectrum but also corrects the distributions of two other quantities that are thrown with

unphysical distributions in the MC: the zenith angle and distance from array center. The

zenith distribution, like the energy distribution, will be different for different varieties of

sources (e.g., an isotropic angular distribution for cosmic rays or a delta function for a

point-like γ source), and so its effective distribution in modeling a specific source must be

52

determined using weights. The core distance from the array center is importance-sampled,

with events near the array center being oversampled because they more often pass event-

selection cuts. The weighting removes this bias in the distribution.

5.1.3 Resulting Functional Form

The function mapping inputs to outputs that is produced by the training process is difficult

to visualize because its domain is 15-dimensional. One way of doing this is to plot the

output as a function of each individual input while freezing the others at their respective

median values. This is done in Figure 5.1. The medians are chosen using MC, and the input

variables are plotted over the ranges they span in this sample.

The dependence of the output on fraction of channels hit is surprising: the slope is

mostly negative. This is because the other inputs, including the fraction of tanks hit, have

been frozen: among events with a given fraction of tanks hit, those with a higher number

of channels hit have a higher density (number of hits per unit area) in the detector. These

events probably have their cores on the array, and such events have a lower energy on average.

With other inputs fixed, the energy estimate is an increasing function of core distance.

This is because an off-array event must have a greater primary energy in order to produce

the same number of hits in the detector.

The inner annular variables produce low energy estimates when their values are very

high. This occurs because only a small, low-energy event can fit almost entirely within the

central region of the detector.

53

Figure 5.1: Dependence of NN energy estimate on each input variable, with the others frozen
at their median values.

54

0.10.20.30.40.50.60.70.80.9Fraction of channels hit11.51212.51313.5 (eV)]E[10log0.20.30.40.50.60.70.80.91Fraction of tanks hit11.51212.51313.5 (eV)]E[10log3.544.555.56(core amplitude)10log11.51212.51313.5 (eV)]E[10log20406080100Core distance (m)11.51212.51313.5 (eV)]E[10log0.750.80.850.90.951cos(zenith angle)11.51212.51313.5 (eV)]E[10log00.10.20.30.40.50.60.70.80.9Fraction of charge between 0 and 10 m of core11.51212.51313.5 (eV)]E[10log00.10.20.30.40.50.60.70.8Fraction of charge between 10 and 20 m of core11.51212.51313.5 (eV)]E[10log00.10.20.30.40.50.60.7Fraction of charge between 20 and 30 m of core11.51212.51313.5 (eV)]E[10logFigure 5.1 (cont’d)

55

00.10.20.30.40.50.6Fraction of charge between 30 and 40 m of core11.51212.51313.5 (eV)]E[10log00.050.10.150.20.250.30.350.40.45Fraction of charge between 40 and 50 m of core11.51212.51313.5 (eV)]E[10log00.050.10.150.20.250.30.350.4Fraction of charge between 50 and 60 m of core11.51212.51313.5 (eV)]E[10log00.050.10.150.20.250.30.350.4Fraction of charge between 60 and 70 m of core11.51212.51313.5 (eV)]E[10log00.050.10.150.20.250.30.350.40.450.5Fraction of charge between 70 and 80 m of core11.51212.51313.5 (eV)]E[10log00.050.10.150.20.250.30.35Fraction of charge between 80 and 90 m of core11.51212.51313.5 (eV)]E[10log00.10.20.30.40.5Fraction of charge more than 90 m from core11.51212.51313.5 (eV)]E[10log5.2 Performance

In this section the performance of several energy estimators is evaluated using photon MC.
Events are weighted to resemble an E−2.63 spectrum as in training, and the cuts applied are

the same as those used in training as well. If the NN becomes “overtrained” on the training

sample, it will perform anomalously well on that sample, better than on independent MC

samples or on data, and so the evaluation MC sample is chosen to be a disjoint set of events

from the training sample in order to avoid overestimating the performance of the estimator.

This analysis is performed only using events that pass background-rejection cuts. Al-

though the MC sample consists only of γ-ray events, standard hadron-rejection cuts are still

imposed to ensure that only events that are liable to be used in an actual analysis are con-

sidered. These cuts (on compactness and PINCness) were chosen to optimize the expected

statistical significance of the Crab Nebula using γ MC and data background events. Each

energy estimator is evaluated without and with a fiducial cut designed to select events well

contained within the HAWC array:

C < 100.

(5.2)

5.2.1 Other Energy Variables

Several energy or energy-proxy variables are evaluated in this section. In addition to the NN

energy, the following will be considered.

56

5.2.1.1 Fraction Hit/Fraction-Hit NN

Here fraction hit is defined as the fraction of PMTs that are hit within an event, excluding bad

hits. This variable represents the shower multiplicity and is therefore positively correlated

with the event’s primary energy. In order to determine what an explicit energy estimate

based on fraction hit would look like, an NN was trained using only fraction hit as an input

variable. This is the fraction-hit NN.

5.2.1.2 Likelihood Energy Estimator

The Likelihood (LH) Energy Estimator uses the maximum-likelihood (ML) formalism to es-

timate the shower energy. The energy estimate is chosen to maximize the likelihood function

associated with the observed event, which is the product over all hits of single-hit likelihood

functions. More detail about this reconstruction algorithm can be found at [1].

5.2.1.3 Ground Parameter Energy Estimator

The Ground Parameter (GP) Energy Estimator uses the value of the shower energy density

estimated at some fraction-hit-dependent distance from the shower core to derive an energy

estimate. This is corrected for zenith as well.

5.2.1.4 Ground Energy Estimator

It is worth considering whether any energy-estimation algorithm could perform better than

a mapping from energy deposited on the ground to primary energy, corrected for zenith

angle. To address this question, a neural net has been trained using two input variables:

MC-truth ground energy and cosine of reconstructed zenith angle. This will be referred to as

the Ground Energy Estimator. This estimator cannot actually be implemented since it uses

57

MC-truth information. It simply represents the best achievable performance by an estimator

seeking to use only the zenith angle and some measure of the energy deposited on the ground

by a shower.

5.2.2 Linearity

Figure 5.2 contains the joint distributions of the energy variables with MC-truth energy. The

identity line is drawn in black for variables with a dimension of energy. The NN Estimator

(Figure 5.2a) displays good linearity above 1 TeV with a bias at the lowest energies due to

the threshold set by the fraction-hit cut applied in the MC sample. Applying the fiducial cut

(Equation (5.2)) removes some poorly reconstructed events at high energies, but otherwise

the NN performs quite well even on off-array events. This is because the NN input variables

already model the containment of the shower within the array (see Table 5.1), allowing the

NN to handle off-array events reasonably well. The most noticeable change is an improvement

in resolution at high energies.

Fraction hit is much less linear than the NN for the following reason. All low-energy

events are inherently low-multiplicity events and therefore have low fraction hit. However,

the converse is not true: a high-energy event, although its shower contains a large number

of particles, may have an arbitrarily low multiplicity in the detector depending on where it

lands. As a result, Figures 5.2c and 5.2e show that at all values of fraction hit, high-energy

events are present. This effect is only somewhat lessened when the fiducial cut is applied

(Figures 5.2d and 5.2f). Because the fraction of PMTs hit cannot exceed one, this variable

saturates at an energy around 30 TeV and cannot distinguish among higher energies.

The LH Estimator (Figure 5.2g) is generally linear above 10 TeV but with a greater

thickness (resolution) to the distribution than the NN displays. This improves somewhat

58

with the addition of the fiducial cut (Figure 5.2h).

Figure 5.2i shows that without the fiducial cut, the GP Estimator suffers from high-

energy contamination at low values of the energy variable due to events reconstructed off

the array. Adding the fiducial cut removes most of these events (5.2j).

The Ground Estimator displays a linearity comparable to that of the NN (Figure 5.2k).

Its performance is not changed much by imposing the fiducial cut (Figure 5.2l) because the

MC-truth ground energy is not sensitive to the fraction of the event contained within the

array. The only change is due to the better zenith-angle estimate for well contained events.

5.2.3 Bin Overlap

Figure 5.3 shows the extent of the overlap in the MC-truth energy content of bins in the

various energy variables. Each bin is represented by a Gaussian with a mean and standard

deviation equal to the mean and standard deviation respectively (in log space) of the MC-

truth energies of the events reconstructed in that bin.

The NN energy bins (Figure 5.3a) show some overlap. The extent of the overlap is

relatively unchanged by the fiducial cut (Figure 5.3b) for the reason described in Section

5.2.2, namely that the NN explicitly handles off array events using the core distance from the

array center. The fraction-hit bins (Figure 5.3c) show a much greater overlap especially at low

energies due to the high-energy contamination of all bins. This also decreases the dynamic

range of the variable. The LH variable (Figure 5.3g) possesses a greater dynamic range than

fraction hit, but there is still significant bin overlap at low energies. Its dynamic range is also

worsened slightly by the fiducial cut (Figure 5.3h). There is substantial overlap in the GP bin

contents without the fiducial cut due to the presence of high-energy events whose energies

are substantially underestimated at high energy (Figure 5.3i). This is improved significantly

59

(a) NN Estimator correlation without fiducial cut.
Correlation coefficient = 83.6%.

(b) NN Estimator correlation with fiducial cut.
Correlation coefficient = 79.9%.

(c) Fraction hit correlation without fiducial cut.
Correlation coefficient = 59.3%.

(d) Fraction hit correlation with fiducial cut.
Correlation coefficient = 66.6%.

Figure 5.2: Correlation between energy variables and MC-truth energy.

60

11-1010-109-108-107-106-105-104-103-10Fraction of Events1010.51111.51212.51313.51414.515[True energy (eV)]10log1010.51111.51212.51313.51414.515[Reconstructed energy (eV)]10log11-1010-109-108-107-106-105-104-103-10Fraction of Events1010.51111.51212.51313.51414.515[True energy (eV)]10log1010.51111.51212.51313.51414.515[Reconstructed energy (eV)]10log11-1010-109-108-107-106-105-104-103-10Fraction of Events1010.51111.51212.51313.51414.515[True energy (eV)]10log1.2-1-0.8-0.6-0.4-0.2-0hitf 10log11-1010-109-108-107-106-105-104-103-10Fraction of Events1010.51111.51212.51313.51414.515[True energy (eV)]10log1.2-1-0.8-0.6-0.4-0.2-0hitf 10logFigure 5.2 (cont’d)

(e) Fraction-hit NN correlation without fiducial
cut. Correlation coefficient = 58.9%.

(f) Fraction-hit NN correlation with fiducial cut.
Correlation coefficient = 66.7%.

(g) LH Estimator correlation without fiducial cut.
Correlation coefficient = 69.1%.

(h) LH Estimator correlation with fiducial cut.
Correlation coefficient = 64.7%.

61

11-1010-109-108-107-106-105-104-103-10Fraction of Events1010.51111.51212.51313.51414.515[True energy (eV)]10log1010.51111.51212.51313.51414.515[Reconstructed energy (eV)]10log11-1010-109-108-107-106-105-104-103-10Fraction of Events1010.51111.51212.51313.51414.515[True energy (eV)]10log1010.51111.51212.51313.51414.515[Reconstructed energy (eV)]10log11-1010-109-108-107-106-105-104-103-10Fraction of Events1010.51111.51212.51313.51414.515[True energy (eV)]10log1010.51111.51212.51313.51414.515[Reconstructed energy (eV)]10log11-1010-109-108-107-106-105-104-103-10Fraction of Events1010.51111.51212.51313.51414.515[True energy (eV)]10log1010.51111.51212.51313.51414.515[Reconstructed energy (eV)]10logFigure 5.2 (cont’d)

(i) GP Estimator correlation without fiducial cut.
Correlation coefficient = 54.7%.

(j) GP Estimator correlation with fiducial cut.
Correlation coefficient = 63.0%.

(k) Ground Estimator correlation without fiducial
cut. Correlation coefficient = 76.0%.

(l) Ground Estimator correlation with fiducial cut.
Correlation coefficient = 70.2%.

62

11-1010-109-108-107-106-105-104-103-10Fraction of Events1010.51111.51212.51313.51414.515[True energy (eV)]10log1010.51111.51212.51313.51414.515[Reconstructed energy (eV)]10log11-1010-109-108-107-106-105-104-103-10Fraction of Events1010.51111.51212.51313.51414.515[True energy (eV)]10log1010.51111.51212.51313.51414.515[Reconstructed energy (eV)]10log11-1010-109-108-107-106-105-104-103-10Fraction of Events1010.51111.51212.51313.51414.515[True energy (eV)]10log1010.51111.51212.51313.51414.515[Reconstructed energy (eV)]10log11-1010-109-108-107-106-105-104-103-10Fraction of Events1010.51111.51212.51313.51414.515[True energy (eV)]10log1010.51111.51212.51313.51414.515[Reconstructed energy (eV)]10logby the fiducial cut although the problem persists (Figure 5.3j). Even the Ground Estimator

(Figures 5.3k and 5.3l) has a smaller dynamic range than the NN, especially at low energies,

because it tends to overestimate the energies of low-energy events more than the NN does.

The raw histograms used to create these distributions can be seen in Figure 5.4.

5.2.4 Bias

Figure 5.5 shows the log-energy bias of the energy variables (except for fraction hit, which

does not have a dimension of energy). Here bias is defined as the mean deviation of the

estimated log energy from its MC-truth value:

(cid:69)

(5.3)

b ≡(cid:68)
(cid:42)

(cid:43)
ˆE − log10 E
ˆE
E

log10

=

log10

where E is the MC-truth energy and ˆE is the estimate thereof. The estimator bias is a way

of quantifying the average offset of the estimator as a function of MC-truth energy. The

ideal bias of 0 is shown as a black line on the plot.

All of the estimators examined here suffer from a positive bias at low energies and a

negative bias at high energies (Figure 5.5a). The former is due to the threshold cut used in

selecting the training and evaluation MC samples. The latter is caused by saturation of the

detector at high energies. The LH estimator displays the least bias of all estimators above

1.5 TeV.

Imposing the fiducial cut has a small effect on the NN performance, worsening the bias

between 1.5 and 10 TeV and improving it elsewhere (Figure 5.5b).

63

(a) NN bin overlap without fiducial cut.

(b) NN bin overlap with fiducial cut.

(c) Fraction-hit bin overlap without fiducial cut.

(d) Fraction-hit bin overlap with fiducial cut.

(e) Fraction-hit-NN bin overlap without fiducial
cut.

(f) Fraction-hit-NN bin overlap with fiducial cut.

Figure 5.3: Gaussians representing distributions of MC-truth energies in reconstructed-
energy or energy-proxy bins. The resulting bin purity (fraction of events not migrating
into adjacent bins) is shown in Figure 5.8.

64

1010.51111.51212.51313.51414.515[True energy (eV)]10log00.20.40.60.81Event Density (Normalized)[Energy Estimate (eV)]10log11.30 - 11.6011.60 - 11.9011.90 - 12.2012.20 - 12.5012.50 - 12.8012.80 - 13.1013.10 - 13.4013.40 - 13.7013.70 - 14.001010.51111.51212.51313.51414.515[True energy (eV)]10log00.20.40.60.81Event Density (Normalized)[Energy Estimate (eV)]10log11.30 - 11.6011.60 - 11.9011.90 - 12.2012.20 - 12.5012.50 - 12.8012.80 - 13.1013.10 - 13.4013.40 - 13.7013.70 - 14.001010.51111.51212.51313.51414.515[True energy (eV)]10log00.20.40.60.81Event Density (Normalized)Fraction Hit0.0655 - 0.10300.1030 - 0.16100.1610 - 0.24600.2460 - 0.35500.3550 - 0.48400.4840 - 0.61700.6170 - 0.73900.7390 - 0.83900.8390 - 1.01001010.51111.51212.51313.51414.515[True energy (eV)]10log00.20.40.60.81Event Density (Normalized)Fraction Hit0.0655 - 0.10300.1030 - 0.16100.1610 - 0.24600.2460 - 0.35500.3550 - 0.48400.4840 - 0.61700.6170 - 0.73900.7390 - 0.83900.8390 - 1.01001010.51111.51212.51313.51414.515[True energy (eV)]10log00.20.40.60.81Event Density (Normalized)[Energy Estimate (eV)]10log12.20 - 12.5012.50 - 12.8012.80 - 13.1013.10 - 13.4013.40 - 13.7013.70 - 14.001010.51111.51212.51313.51414.515[True energy (eV)]10log00.20.40.60.81Event Density (Normalized)[Energy Estimate (eV)]10log12.20 - 12.5012.50 - 12.8012.80 - 13.1013.10 - 13.4013.40 - 13.7013.70 - 14.00Figure 5.3 (cont’d)

(g) LH bin overlap without fiducial cut.

(h) LH bin overlap with fiducial cut.

(i) GP bin overlap without fiducial cut.

(j) GP bin overlap with fiducial cut.

(k) Ground-energy bin overlap without fiducial
cut.

(l) Ground-energy bin overlap with fiducial cut.

65

1010.51111.51212.51313.51414.515[True energy (eV)]10log00.20.40.60.81Event Density (Normalized)[Energy Estimate (eV)]10log11.60 - 11.9011.90 - 12.2012.20 - 12.5012.50 - 12.8012.80 - 13.1013.10 - 13.4013.40 - 13.7013.70 - 14.001010.51111.51212.51313.51414.515[True energy (eV)]10log00.20.40.60.81Event Density (Normalized)[Energy Estimate (eV)]10log11.60 - 11.9011.90 - 12.2012.20 - 12.5012.50 - 12.8012.80 - 13.1013.10 - 13.4013.40 - 13.7013.70 - 14.001010.51111.51212.51313.51414.515[True energy (eV)]10log00.20.40.60.81Event Density (Normalized)[Energy Estimate (eV)]10log11.30 - 11.6011.60 - 11.9011.90 - 12.2012.20 - 12.5012.50 - 12.8012.80 - 13.1013.10 - 13.4013.40 - 13.7013.70 - 14.001010.51111.51212.51313.51414.515[True energy (eV)]10log00.20.40.60.81Event Density (Normalized)[Energy Estimate (eV)]10log11.30 - 11.6011.60 - 11.9011.90 - 12.2012.20 - 12.5012.50 - 12.8012.80 - 13.1013.10 - 13.4013.40 - 13.7013.70 - 14.001010.51111.51212.51313.51414.515[True energy (eV)]10log00.20.40.60.81Event Density (Normalized)[Energy Estimate (eV)]10log11.30 - 11.6011.60 - 11.9011.90 - 12.2012.20 - 12.5012.50 - 12.8012.80 - 13.1013.10 - 13.4013.40 - 13.7013.70 - 14.001010.51111.51212.51313.51414.515[True energy (eV)]10log00.20.40.60.81Event Density (Normalized)[Energy Estimate (eV)]10log11.30 - 11.6011.60 - 11.9011.90 - 12.2012.20 - 12.5012.50 - 12.8012.80 - 13.1013.10 - 13.4013.40 - 13.7013.70 - 14.00(a) NN bin content without fiducial cut.

(b) NN bin content with fiducial cut.

(c) Fraction-hit bin content without fiducial cut.

(d) Fraction-hit bin content with fiducial cut.

(e) Fraction-hit-NN bin content without fiducial
cut.

(f) Fraction-hit-NN bin content with fiducial cut.

Figure 5.4: Distribution of MC-truth energies in energy bins.

66

1010.51111.51212.51313.51414.515[True energy (eV)]10log7-106-105-104-103-102-101-10Fraction of Events[Energy Estimate (eV)]10log11.00 - 11.3011.30 - 11.6011.60 - 11.9011.90 - 12.2012.20 - 12.5012.50 - 12.8012.80 - 13.1013.10 - 13.4013.40 - 13.7013.70 - 14.001010.51111.51212.51313.51414.515[True energy (eV)]10log7-106-105-104-103-102-101-10Fraction of Events[Energy Estimate (eV)]10log11.00 - 11.3011.30 - 11.6011.60 - 11.9011.90 - 12.2012.20 - 12.5012.50 - 12.8012.80 - 13.1013.10 - 13.4013.40 - 13.7013.70 - 14.001010.51111.51212.51313.51414.515[True energy (eV)]10log7-106-105-104-103-102-101-10Fraction of EventsFraction Hit0.0437 - 0.06550.0655 - 0.10300.1030 - 0.16100.1610 - 0.24600.2460 - 0.35500.3550 - 0.48400.4840 - 0.61700.6170 - 0.73900.7390 - 0.83900.8390 - 1.01001010.51111.51212.51313.51414.515[True energy (eV)]10log7-106-105-104-103-102-101-10Fraction of EventsFraction Hit0.0437 - 0.06550.0655 - 0.10300.1030 - 0.16100.1610 - 0.24600.2460 - 0.35500.3550 - 0.48400.4840 - 0.61700.6170 - 0.73900.7390 - 0.83900.8390 - 1.01001010.51111.51212.51313.51414.515[True energy (eV)]10log7-106-105-104-103-102-101-10Fraction of Events[Energy Estimate (eV)]10log11.00 - 11.3011.30 - 11.6011.60 - 11.9011.90 - 12.2012.20 - 12.5012.50 - 12.8012.80 - 13.1013.10 - 13.4013.40 - 13.7013.70 - 14.001010.51111.51212.51313.51414.515[True energy (eV)]10log7-106-105-104-103-102-101-10Fraction of Events[Energy Estimate (eV)]10log11.00 - 11.3011.30 - 11.6011.60 - 11.9011.90 - 12.2012.20 - 12.5012.50 - 12.8012.80 - 13.1013.10 - 13.4013.40 - 13.7013.70 - 14.00Figure 5.4 (cont’d)

(g) LH bin content without fiducial cut.

(h) LH bin content with fiducial cut.

(i) GP bin content without fiducial cut.

(j) GP bin content with fiducial cut.

(k) Ground-energy bin content without fiducial
cut.

(l) Ground-energy bin content with fiducial cut.

67

1010.51111.51212.51313.51414.515[True energy (eV)]10log7-106-105-104-103-102-101-10Fraction of Events[Energy Estimate (eV)]10log11.00 - 11.3011.30 - 11.6011.60 - 11.9011.90 - 12.2012.20 - 12.5012.50 - 12.8012.80 - 13.1013.10 - 13.4013.40 - 13.7013.70 - 14.001010.51111.51212.51313.51414.515[True energy (eV)]10log7-106-105-104-103-102-101-10Fraction of Events[Energy Estimate (eV)]10log11.00 - 11.3011.30 - 11.6011.60 - 11.9011.90 - 12.2012.20 - 12.5012.50 - 12.8012.80 - 13.1013.10 - 13.4013.40 - 13.7013.70 - 14.001010.51111.51212.51313.51414.515[True energy (eV)]10log7-106-105-104-103-102-101-10Fraction of Events[Energy Estimate (eV)]10log11.00 - 11.3011.30 - 11.6011.60 - 11.9011.90 - 12.2012.20 - 12.5012.50 - 12.8012.80 - 13.1013.10 - 13.4013.40 - 13.7013.70 - 14.001010.51111.51212.51313.51414.515[True energy (eV)]10log7-106-105-104-103-102-101-10Fraction of Events[Energy Estimate (eV)]10log11.00 - 11.3011.30 - 11.6011.60 - 11.9011.90 - 12.2012.20 - 12.5012.50 - 12.8012.80 - 13.1013.10 - 13.4013.40 - 13.7013.70 - 14.001010.51111.51212.51313.51414.515[True energy (eV)]10log7-106-105-104-103-102-101-10Fraction of Events[Energy Estimate (eV)]10log11.00 - 11.3011.30 - 11.6011.60 - 11.9011.90 - 12.2012.20 - 12.5012.50 - 12.8012.80 - 13.1013.10 - 13.4013.40 - 13.7013.70 - 14.001010.51111.51212.51313.51414.515[True energy (eV)]10log7-106-105-104-103-102-101-10Fraction of Events[Energy Estimate (eV)]10log11.00 - 11.3011.30 - 11.6011.60 - 11.9011.90 - 12.2012.20 - 12.5012.50 - 12.8012.80 - 13.1013.10 - 13.4013.40 - 13.7013.70 - 14.00(a) All estimators with fiducial cut.

(b) NN estimator with and without fiducial cut.

Figure 5.5: Bias of energy variables.

68

1010.51111.51212.51313.51414.515[True energy (eV)]10log0.5-00.511.5Log-energy biasNeural NetFraction-Hit Neural NetLikelihoodGround ParameterGround1010.51111.51212.51313.51414.515[True energy (eV)]10log1-0.8-0.6-0.4-0.2-00.20.40.60.8Log-energy biasNeural Net with fiducial cutNeural Net without fiducial cut5.2.5 Resolution

Figure 5.6 shows the energy variables’ log-energy resolution, defined as the standard deviation

of the log energy estimate:

(cid:115)(cid:28)(cid:16)

r ≡

ˆE −(cid:68)

log10

(cid:69)(cid:17)2(cid:29)

.

log10

ˆE

(5.4)

This quantity is a measure of the width of the estimator’s distribution. Figure 5.6b shows

that the NN energy resolution is improved by the addition of the fiducial cut, particularly at

high energies. With the cut, the NN displays a significantly smaller resolution than the LH

Estimator, GP Estimator, or Ground Estimator between 100 GeV and 68 TeV (Figure 5.6a).

All the estimators’ resolutions become small at low energies. This should not be understood

to mean that the energy reconstruction is best at low energies, for all estimators are also

characterized by substantial biases in this regime (see Figure 5.5).

5.2.6 RMS Error

Figure 5.7 shows the log-energy RMS error of the different energy variables. This is defined

as the RMS deviation of the estimated log energy from the MC-truth log energy:

(cid:17)2(cid:29)

(5.5)

(cid:115)(cid:28)(cid:16)
(cid:118)(cid:117)(cid:117)(cid:116)(cid:42)

ρ ≡

=

log10

ˆE − log10 E
(cid:43)

log2
10

ˆE
E

.

69

(a) All estimators with fiducial cut.

(b) NN estimator with and without fiducial cut.

Figure 5.6: Resolution of energy variables.

70

1010.51111.51212.51313.51414.515[True energy (eV)]10log00.050.10.150.20.250.3Log-energy resolutionNeural NetFraction-Hit Neural NetLikelihoodGround ParameterGround1010.51111.51212.51313.51414.515[True energy (eV)]10log00.050.10.150.20.250.30.350.4Log-energy resolutionNeural Net with fiducial cutNeural Net without fiducial cutThis quantity is also equal to the sum in quadrature of the bias and resolution:

(cid:112)

ρ =

b2 + r2.

(5.6)

Without the fiducial cut, the NN displays an RMS error that is roughly constant, with a

value around 0.17, in the interval from 1 to 32 TeV (Figure 5.7b). The RMS error is smaller

for the NN than for all other estimators at all energies.

5.2.7 Bin Contamination

In a spectral analysis, there is the concern that the highest-energy bins may have a low purity

(i.e. they may contain mostly misreconstructed lower-energy events). This can happen if the

spectrum is sharply decreasing, meaning that there are far more lower-energy events, which

may fluctuate upwards, than there are higher-energy events. Figure 5.8 shows, for the NN,

LH, and GP Estimators, the fraction of events reconstructed in each half-decade energy bin
whose MC-truth energies are also within the bin, assuming an E−2.63 spectrum. With the

fiducial cut (Figure 5.8b), the NN bin purity is highest at all energies and increases with

energy. The purity is 86% in the highest-energy NN bin with the fiducial cut.

5.2.8 Agreement of Estimators

As a check on the consistency of the various energy variables, it is interesting to plot them

against each other. This is done in Figure 5.9. A fiducial cut is applied. The estimators

are all fairly correlated with each other. There are classes of events that lie away from the

identity line, and these can be understood in terms of the similar classes of events in Figures

5.2h and 5.2j.

71

(a) All estimators with fiducial cut.

(b) NN estimator with and without fiducial cut.

Figure 5.7: RMS error of energy variables.

72

1010.51111.51212.51313.51414.515[True energy (eV)]10log00.20.40.60.811.21.4Log-energy RMS errorNeural NetFraction-Hit Neural NetLikelihoodGround ParameterGround1010.51111.51212.51313.51414.515[True energy (eV)]10log00.10.20.30.40.50.60.70.80.9Log-energy RMS errorNeural Net with fiducial cutNeural Net without fiducial cut(a) NN bin purity without fiducial cut.

(b) NN bin purity with fiducial cut.

(c) LH bin purity without fiducial cut.

(d) LH bin purity with fiducial cut.

(e) GP bin purity without fiducial cut.

(f) GP bin purity with fiducial cut.

Figure 5.8: Fraction of MC events reconstructed in each NN bin whose MC-truth energies
actually lie in that bin.

73

1212.51313.51414.5[Reconstructed energy (eV)]10log0.50.550.60.650.70.75Bin purity1212.51313.51414.5[Reconstructed energy (eV)]10log0.660.680.70.720.740.760.780.80.820.840.86Bin purity1212.51313.51414.5[Reconstructed energy (eV)]10log0.540.560.580.60.620.64Bin purity1212.51313.51414.5[Reconstructed energy (eV)]10log0.50.550.60.650.7Bin purity1212.51313.51414.5[Reconstructed energy (eV)]10log0.550.60.650.70.75Bin purity1212.51313.51414.5[Reconstructed energy (eV)]10log0.550.60.650.70.75Bin purity(a) Correlation coefficient = 90.9%.

(b) Correlation coefficient = 91.7%.

(c) Correlation coefficient = 93.3%.

Figure 5.9: Correlations of three energy variables.

74

11-1010-109-108-107-106-105-104-103-10Fraction of events1010.51111.51212.51313.51414.515[Neural Net energy (eV)]10log1010.51111.51212.51313.51414.515[Likelihood energy (eV)]10log11-1010-109-108-107-106-105-104-103-10Fraction of events1010.51111.51212.51313.51414.515[Neural Net energy (eV)]10log1010.51111.51212.51313.51414.515[Ground Parameter energy (eV)]10log11-1010-109-108-107-106-105-104-103-10Fraction of events1010.51111.51212.51313.51414.515[Likelihood energy (eV)]10log1010.51111.51212.51313.51414.515[Ground Parameter energy (eV)]10logChapter 6

Crab Nebula spectral reconstruction

The Crab Nebula is a PWN located at 5 hours, 34 minutes, and 32 seconds right ascension

(R.A.) and 22 degrees, 0 arcminutes, and 52 arcseconds declination. It is the brightest γ-ray

source seen by HAWC, and a spectral analysis of this source using the NN energy variable

is therefore an excellent proof of principle for the NN technique. Specifically the application

of the NN energy-estimation method to the Crab allows the Crab’s energy spectrum to be

resolved at unprecedentedly high energies.

6.1 Skymaps

The first step in HAWC’s high-level data analysis is the production of skymaps. A skymap

is a kind of two-dimensional histogram representing the distribution of the reconstructed

directions of events detected by HAWC. The bins in this histogram are referred to as pixels,

and the pixelization of the sky is determined using the HEALPix [55] software. Each pixel

subtends a solid angle of 0.00328 square degree. Skymap files are stored using the FITS [56]

data format.

The skymaps are produced using 836 sidereal days of HAWC data. These data are

subjected to cuts on PINCness and compactness that are chosen to optimize the expected

statistical significance of the Crab Nebula, where this expectation is computed using γ MC

and background data.

75

Events are binned using a two-dimensional binning scheme with rectangular bins in frac-

tion hit and NN energy, and one skymap is produced for each analysis bin. The fraction-hit

bin edges are 0.067, 0.105, 0.162, 0.247, 0.356, 0.485, 0.618, 0.740, 0.840, and 1.000. The
NN-energy bins are each one quarter decade wide and span the range from 10−1/2 TeV to

105/2 TeV (316 GeV to 316 TeV). These are sometimes referred to by letter names ranging

from a to l. Since the analysis technique described below relies on being able to characterize

the distributions of photon and hadron events using MC and data respectively in each bin,

it is necessary to exclude from the analysis bins that are very poorly populated, which are

generally bins with high fraction hit and low NN energy or vice versa. Fortunately these

low-statistics bins contribute only negligibly to the statistical power of the analysis. Only

37 of the 108 two-dimensional bins are expected to contribute nonnegligibly to the Crab’s

statistical significance. The other bins are not included.

Some hadron events survive the compactness and PINCness cuts, and so it is necessary

to estimate the amount of background in the maps in order to correctly model the expected

event count in each pixel during the spectral analysis. For bins with reasonably good statis-

tics, this is done using a technique called direct integration (DI) [57, 58]. The background

rate is modeled as a function of hour angle (H.A.) and declination, estimated for each two-

hour interval by integrating the data over that period while excluding known bright γ-ray

sources such as the Crab, Markarians 421 and 501, and the Galactic plane.

In bins with lower statistics, DI becomes imprecise. For bins with total counts below
5 × 105, an approach called background randomization is used instead. For each bin, a

two-dimensional histogram of H.A. and declination is populated using the entire data set.

This distribution is then sampled 10000 times, and each resulting (H.A., declination) pair is

combined with the time from one randomly chosen event from the data set, and the pixel

76

in which this pseudoevent would have fallen is incremented. The background map is then

scaled to have the appropriate normalization.

Following DI or background randomization, an additional smoothing is performed on the
background maps using a flat, circular kernel with radius 0.5◦. This further damps out local

fluctuations in the background maps.

Figures 6.1 and 6.2 show the near-Crab regions of the count and background maps re-

spectively for the 2D analysis bin consisting of events with reconstructed energies between

3.16 and 5.62 TeV from fraction-hit bin 3. This bin’s background is estimated using DI.

Figures 6.3 and 6.4 show the count and background maps for events with energies recon-

structed between 56.2 and 100 TeV coming from fraction-hit bin 9. This bin’s background is

computed using background randomization. The background map shows that the per pixel

background expectation is a fraction of an event in all pixels in this lower-statistics bin.

6.2 Binned-likelihood formalism

The Crab energy spectrum is modeled as a log parabola:

= Φ0 (E/E0)−α−β ln(E/E0)

dN
dE

(6.1)

where Φ0, α, and β are free parameters. E0, the pivot energy, is fixed. Its value is chosen

to decorrelate Φ0 and α in the spectral fit. 7 TeV has been empirically determined to be a

suitable value for E0 for the Crab Nebula.

The log parabola’s name refers to the fact that the logarithm of the flux depends quadrat-

77

Figure 6.1: The skymap for fraction-hit bin 3 with NN energies between 3.16 and 5.62 TeV.

78

8384 []2223 []0369121518212427CountFigure 6.2: The background skymap for fraction-hit bin 3 with NN energies between 3.16
and 5.62 TeV. This bin uses the DI background calculation.

79

8384 []2223 []0369121518212427BackgroundFigure 6.3: The skymap for fraction-hit bin 9 with NN energies between 56.2 and 100.0 TeV.

80

8384 []2223 []01234CountFigure 6.4: The background skymap for fraction-hit bin 9 with NN energies between 56.2
and 100.0 TeV. This bin uses background randomization.

81

8384 []2223 []01234Backgroundically on the logarithm of the energy:

(cid:18) dN/dE

(cid:19)

(cid:19)2

(cid:18)

ln

E
E0

= −α ln

− β

E
E0

(6.2)

ln

Φ0

The three free parameters in Equation (6.1) are determined using a maximum-likelihood

analysis of the skymaps. The number of events falling into a given pixel is modeled as a

Poisson variate with a mean equal to the background expectation in that pixel plus the signal

expectation under the log-parabola model of the Crab spectrum:

(cid:89)

e−(bi+si) (bi + si)ci

i

ci!

L(Φ0, α, β) =

(6.3)

where L is the likelihood and bi, si, and ci are the background mean, signal mean, and

observed count in the ith combination of skymap and pixel.

Calculating the expected signal si in a given pixel requires knowledge of the detector

response. This consists of the distribution of reconstructed energies produced by a given

true energy and the distribution of reconstructed shower directions produced by a given true

direction. The latter is known as the point-spread function (PSF). These distributions are

built using MC and are stored in a ROOT file called a detector-response file. For each analysis

bin, the file contains the distribution of reconstructed energies produced by a nominal true-
energy distribution of E−2.63, which can then be reweighted for an arbitrary spectrum in

order to evaluate its likelihood. The file also contains for each bin a histogram of the angular

distances between MC events’ true and reconstructed directions and a double-Gaussian fit

82

to this distribution:

ρ(r) =

A
2π

(cid:40)

f
σ2
1

(cid:34)

exp

(cid:19)2(cid:35)

(cid:18) r

σ1

−1
2

+

1 − f
σ2
2

(cid:34)

exp

(cid:19)2(cid:35)(cid:41)

(cid:18) r

σ2

−1
2

(6.4)

where ρ is the number of events per unit reconstructed solid angle at a distance r from the

true direction and the normalization A, the radii σ1 and σ2, and the coefficient f float in

the fit. The double-Gaussian function is chosen to accomodate two classes of events: well

reconstructed and poorly reconstructed. Figure 6.5 shows the 68% quantiles of the double-

Gaussian PSFs fit to data and MC. In most bins, the MC underestimates the PSF width.

This is treated as a systematic error, described in Section 6.3.

It is useful to define a log-likelihood-ratio (LLR) test statistic (TS) comparing the signal-

plus-background hypothesis described above to a background-only hypothesis with no signal

(Φ0 = 0). The likelihood of the latter model is

(cid:89)

e−bibi
ci!

ci

.

Lbkg =

The TS is then defined as

i

(6.5)

(6.6)

(6.7)

.

D ≡ 2 ln

maxΦ0,α,β L(Φ0, α, β)

= 2

ci ln

(cid:20)

(cid:88)

i

(cid:19)

(cid:18)

Lbkg
si
bi

1 +

(cid:21)

− si

This TS indicates the extent to which the signal-plus-background hypothesis better mod-

els the data than the background-only hypothesis. According to Wilks’ theorem [59], the

TS is distributed under the background-only hypothesis according to a χ2 distribution with

83

Figure 6.5: The 68% quantile of the distribution of the angular distance between the true
and reconstructed directions of events in data and MC.

84

1c2c3c4c2d3d4d5d2e3e4e5e6e4f5f6f7f8f5g6g7g8g9g6h7h8h9h7i8i9i9j9kBin name0.20.40.60.81.068% Containment Radius (deg)SimulationMeasured Crabthree degrees of freedom in the high-statistics limit. The TS can therefore be employed to

compute the statistical significance with which HAWC detects an excess of events from the

source above background.

Statistical uncertainties of the spectral parameters can also be computed by means of

Wilks’ theorem. Considering for instance Φ0, the LLR

0,α(cid:48),β(cid:48) L(cid:0)Φ(cid:48)

0, α(cid:48), β(cid:48)(cid:1)

2 ln

maxΦ(cid:48)
maxα(cid:48),β(cid:48) L(Φ0, α(cid:48), β(cid:48))

(6.8)

should be asymptotically distributed according to a χ2 distribution with one degree of free-

dom if Φ0 assumes its true value. This LLR can therefore be evaluated for a given value of

Φ0 to test whether that value is compatible with the observed data, and the set of values

of Φ0 at which the p-value derived from this χ2 is not less than 0.32 represents the 68%-

containment (1σ) error band for Φ0. The same procedure can be carried out for α and β as

well.

6.3 Flux-points calculation

It is also desirable to be able to compute a flux estimate in a way that is less reliant on the

assumption that the true Crab emission spectrum is exactly a log parabola. A second style

of analysis is carried out in addition to the one described in Section 6.2. In this analysis,

an additional fit is done for each NN-energy bin using only the two-dimensional bins falling

into that energy bin and allowing only Φ0 to float in the fit. The result is an estimate of

the flux in each energy bin that is obtained assuming that the shape of the spectrum in that

bin is still approximately log-parabolic but allowing the event rate in that bin to deviate

85

from that predicted by the log-parabola fit from Section 6.2. The computed flux for each

reconstructed-energy bin is plotted at the mean true log energy contributing to that bin

assuming the best-fit log-parabola shape.

This choice places the point at an energy that is representative of the true energy at

which the flux is being measured only if the functional form of the source’s true energy spec-

trum is approximately log-parabolic in shape. It is not sufficient for detecting substantial,

qualitative deviations from the log-parabola functional form, such as an LIV-induced hard

cutoff. Chapter 7 will discuss a technique for doing this.

The calculation of expected counts in skymap pixels is reliant upon the MC simulation of

HAWC described in Section 4.6. The computed spectrum is therefore subject to systematic

uncertainties deriving from the uncertainties in the MC simulation described in Section 4.6.

To compute these uncertainties, the spectral analysis is performed using dedicated system-

atics MC simulations in which the charge uncertainty, quantum efficiency, PMT threshold,

and broad-pulse effects have been assigned values at the extremes of the ranges described in

Section 4.6. In order to account for the data/MC discrepancy in the angular resolution, an

additional iteration of the spectral analysis is performed in which the PSF is modeled using

the double-Gaussian parameters measured using data. Each of these effects is varied one at

a time, and the resulting uncertainties on the flux points are added in quadrature with each

other and with the statistical errors derived from the likelihood fitting procedure.

6.4 Fit result

Figures 6.6 and 6.7 show the results of the log-parabola and flux-points analyses of HAWC

Crab data respectively. Table 6.1 contains the best-fit values of the log-parabola spectral

86

Parameter

Φ0[10−13/(cm2 · s · TeV)]

α
β

Value Stat. error
3.27
2.623
0.08

0.05
0.014
0.01

Syst. error
+1.10/−0.76
+0.081/−0.109
+0.06/−0.05

Table 6.1: Best-fit values of log-parabola spectral parameters with statistical and systematic
uncertainties. Systematic uncertainties are the sum in quadrature of the deviations from the
nominal result produced by performing the analysis using the systematics MCs as described
in Section 6.3.

parameters along with their uncertainties, and Table 6.2 contains the values of the quantity

energy squared times flux at each energy plotted in Figure 6.7 along with uncertainties.

The error band in Figure 6.6 shows the statistical uncertainty on the flux propagated from

the uncertainties on the fit parameters. The error bars in Figure 6.7 are statistical errors

added in quadrature with the systematic errors described in Section 6.3. Figure 6.7 shows

that the HAWC measurement of the Crab flux at energies between 10 and 50 TeV exceeds

that measured by H.E.S.S. and HEGRA. The result also constitutes the highest-energy

observation of the Crab γ-ray spectrum by any experiment. The production of γ rays by the

Crab at these energies has been predicted by multiple theoretical models [60, 61]; however,

this measurement is the first empirical measurement.

Figure 6.8 shows the “pulls” for the log-parabola fit to the Crab. The pull for a bin is

defined as the residual divided by the expected fluctuation of that residual under the best-fit

model:

t ≡ c − s − b

√
s + b

(6.9)

where c, s, and b are the observed count, expected signal, and expected background in the

bin added over pixels that lie within a certain angular distance of the Crab. The distance

is chosen for each bin to be the radius that would optimize the sensitivity of a hypothetical

√
aperture-based analysis (i.e. it maximizes the quantity s/

b). The resulting quantity t

87

Figure 6.6: The Crab log-parabola spectrum fit using HAWC data binned in NN energy,
along with similar measurements made by the H.E.S.S. [4] and HEGRA [5] experiments.
Errors on the HAWC spectrum are statistical only. See Figure 6.7 for systematic errors.

Energy Energy squared times flux Stat. error Syst. error
0.8 +2.7/−1.4
37.8
0.6 +5.1/−3.9
29.6
0.5 +7.1/−4.6
22.8
0.4 +6.5/−4.1
15.8
0.4 +5.5/−2.9
11.0
0.4 +2.9/−2.0
7.4
0.4 +2.3/−1.3
4.7
1.8 +0.4/−0.3 +1.9/−0.5
0.4 +0.6/−0.3
1.1

1.2
2.2
3.8
6.9
12.6
22.7
40.4
72.8
126.0

Table 6.2: Energy (in TeV) and energy squared times flux (in TeV/(cm2 · s)) for each point
in Figure 6.7. Statistical errors are propagated from the uncertainties on the flux normal-
ization as refit in each energy bin as described in Section 6.3. Systematic errors are sums
in quadrature of deviations of the flux points computed using systematics MCs from the
nominal flux points.

88

100101102E (TeV)1013101210111010E2dN/dE (TeV/cm2/s)Neural NetH.E.S.S.HEGRAFigure 6.7: The Crab energy spectrum produced by the flux-points analysis described in
Section 6.3, along with measurements by H.E.S.S. [4] and HEGRA [5]. Errors are statistical
and systematic added in quadrature.

89

100101102E (TeV)1013101210111010E2dN/dE [TeV/(cm2 s)]Neural NetH.E.S.S.HEGRAFigure 6.8: Pulls (residuals divided by their expected fluctuations) for each of the 2D bins
included in the Crab analysis. The x axis indicates NN bins; the color scale indicates
fraction-hit bins.

should have a mean of zero and a standard deviation of one if the data are compatible with

the model within statistical errors. However, a substantial negative bias can be seen in

the pulls. This results from the systematic underestimation of HAWC’s angular resolution

by the MC (Figure 6.5). Because a fixed radius is predicted by the MC to contain a larger

fraction of events than it actually does in data, the MC-based prediction exceeds the observed

count. This source of systematic error is accounted for via the angular-resolution systematic

described in Section 6.3.

The results presented in this chapter demonstrate that the HAWC experiment can make

90

100101102103E (TeV)42024PullBin 1Bin 2Bin 3Bin 4Bin 5Bin 6Bin 7Bin 8Bin 9unprecedentedly high-energy measurements of γ-ray flux using the NN energy variable. How-

ever, measurements of the kind shown in Figure 6.7 are not in and of themselves sufficient

to constrain LIV. The flux-points technique relies on the assumption that true emission

spectrum of the source is at least approximately log-parabolic in shape. In particular, the

kind of hard cutoff that would be expected under LIV would qualitatively violate this as-

sumption. A bin containing an excess of events reconstructed at high energies would still be

assigned a positive flux and a high energy even if all of these events were lower-energy events

whose energies were overestimated by the NN. Chapter 7 will present a technique by which

HAWC data reconstructed using the NN algorithm can be used to set constraining limits

on photon-sector Lorentz-invariance violation and demonstrate the presence of high-energy

flux.

91

Chapter 7

LIV limits

HAWC’s observations of unprecedentedly high-energy γ rays provide a vantage point from

which to search for and constrain the sort of hard spectral cutoff that might occur in the

presence of photon-sector LIV as described in Section 1.3.

7.1 Analysis technique

According to Section 1.3, a photon dispersion relation like Equation (1.12) (in which the

LIV perturbation contributes with the same sign as the k2 term) allows photon decay to

electron/positron pair for a photon momentum k satisfying Equation (1.17). For photons

propagating on astrophysical distance scales, the decay probability above the energy cutoff

is very nearly one [9]. As a result, if such a decay is allowed, the observed energy spectrum

from an astrophysical source of very-high-energy γ rays should have the distinct feature of

a hard cutoff, i.e. an energy above which there is absolutely no emission. The blue curve in

Figure 7.1 shows the effect of such a cutoff on a log-parabola spectrum.

However, because HAWC’s energy resolution is imperfect, an event with true energy

below the cutoff energy still has a nonzero probability of being assigned a reconstructed

energy above the cutoff. This is shown in the orange curve in Figure 7.1. As a result, it

cannot simply be concluded that there is no LIV-induced cutoff below the highest energy

assigned to any photon in HAWC data.

Instead a profile-likelihood technique is used to

92

Figure 7.1: A log-parabola energy spectrum subject to a hard cutoff. The blue curve shows
the true-energy spectrum; the orange curve shows the reconstructed-energy distribution.

determine a 95%-confidence lower limit on any LIV cutoff that may be present.

In this analysis the data are binned in a manner similar to that described in Section

6.1 with two modifications. First, each of the top two reconstructed-energy bins, the one

spanning 100 TeV to 178 TeV and the one spanning 178 TeV to 316 TeV, is subdivided into

three bins of equal size in log space. This finer binning is done to increase the precision with

which the high-energy spectral shape can be resolved. Second, the pixel-by-pixel analysis is

replaced with an aperture- or “top-hat”-based analysis: instead of treating each pixel as a

separate bin in the likelihood, we add the contents of all pixels within a certain radius of

the center of the source, and the likelihood is constructed using only these summed counts.

93

log Elog(E2dN/dE)True spectrumConvolved with resolutionSource name

Crab
2HWC J1825−134
J1839−057
2HWC J1844−032
2HWC J1908+063
2HWC J2019+367
2HWC J2031+415

Top-hat radius (◦)
0.6
0.4
3.2
1.1
0.9
0.8
0.4

Table 7.1: Top-hat radii for the seven sources.

This reduces the systematic dependency of the result on the modeling of the PSF, which is

known to be poor in the MC (as shown in Figure 6.5).

The top-hat radius is chosen independently for each source and used in all bins in the

analysis of that source. The limit-setting procedure is carried out using successively larger

values of the top-hat radius until the value levels off; the limit improves as the radius is

increased until most or all of the high-energy signal photons are within the top hat. The

chosen top-hat radii are shown in Table 7.1. Figure 7.2 shows the top-hat optimization

procedure for 2HWC J1825-134.

Observations of seven sources are used. These have been chosen for the presence of
high-energy emission in HAWC data. The sources are the Crab Nebula, 2HWC J1825−134,
J1839−057, 2HWC J1844−032, 2HWC J1908+063, 2HWC J2019+367, and 2HWC J2031+415.

The 2HWC sources are all in the vicinities of known PWNs or SNRs and are described in the

2HWC catalog [62]. J1839-057 is an unidentified source that is near the SNRs Kes 73 and

G26.6-0.1 [62]. The Crab Nebula is modeled as a point source with a log-parabola spectrum

as in Chapter 6. The remaining sources are treated as extended; the morphology is modeled

as a two-dimensional Gaussian distribution centered at the location of the source (known

from other experiments and shown in Table 7.3). The width σ of this Gaussian is fixed to

a value determined by a prior HAWC analysis of these sources. This analysis used a pixel-

94

Figure 7.2: 95%-confidence lower limit on ELIV using 2HWC J1825-134 data as a function
of top-hat radius.

95

0.10.20.30.40.50.6Top-hat radius (deg)1024×1016×1012×10295% lower cutoff limit (TeV)Chosen top hatby-pixel likelihood technique and was therefore able to estimate the sources’ morphologies.

These sources’ energy spectra are modeled as exponentially cut off power laws:

(cid:18) E

(cid:19)−α

(cid:18)

(cid:19)

dN
dE

= Φ0

E0

exp

− E
Eexp

(7.1)

where Eexp is the energy at which the spectrum cuts off.

Since in this analysis we are concerned only with the LIV cutoff energy, the parameters

of the sources’ emission spectra will be treated as nuisance parameters. These are Φ0, α, and

in the case of the Crab β or in the case of the other sources Eexp. The nuisance parameters

are dealt with using the profile-likelihood formalism [63]. Let θ be the vector of nuisance

parameters ((cid:2)Φ0 α β(cid:3)T or(cid:2)Φ0 α Eexp

(cid:3)T). Again using Wilks’ theorem, we can consider

an LLR of the form

D(ELIV) ≡ 2 ln

maxθ(cid:48),E(cid:48)

L(θ(cid:48), E(cid:48)
maxθ(cid:48) L(θ(cid:48), ELIV)

LIV

LIV)

(7.2)

where ELIV is the LIV-induced hard-cutoff energy. This will be asymptotically distributed

as a χ2 with one degree of freedom. This LLR can be rewritten as

D(ELIV) = 2 ln

maxE(cid:48)

Lp(E(cid:48)
LIV
Lp(ELIV)

LIV)

where Lp is the so called profile likelihood, defined as

Lp(ELIV) ≡ max

θ(cid:48) L(cid:0)θ(cid:48), ELIV

(cid:1) .

(7.3)

(7.4)

The equivalency of Equations (7.2) and (7.3) implies that we can treat the profile likelihood

as if it were the ordinary likelihood for a one-parameter model.

96

Under the hypothesis that there is no LIV cutoff (ELIV = ∞), the TS D(∞) can be used

to perform a likelihood-ratio test to determine the significance with which the data deviate

from their expectations under the null (no-cutoff) hypothesis. Under the null hypothesis,
Wilks’ theorem states that D(∞) should be distributed as a χ2 with one degree of freedom.

A p-value quantifying the statistical significance of the deviation from the null hypothesis is

computed by integrating the χ2 distribution with one degree of freedom from the value of
D(∞) to ∞. These values are reported in Table 7.3.

Since a 95%-confidence lower limit is simply the low edge of a 90%-confidence uncertainty

interval, such a lower limit on ELIV can then be computed by solving

D(ELIV) = F−1(0.9) = 2.71

(7.5)

where F is the cumulative distribution for a χ2 variable with one degree of freedom. Such

limits are computed for the nominal and systematics MCs for each source and listed in Table

7.3. A systematic uncertainty for each nominal limit is calculated by adding in quadrature

the deviations of the systematics limits from the nominal limits. The final result, a degraded

limit, is then computed as the nominal limit minus its lower systematic error.

The 95%-confidence limit is computed in two parts, shown in Figure 7.3. First the

likelihood is maximized over all parameters (ELIV and the nuisance parameters). The profile

likelihood is then computed at successive lower values of ELIV until the profile likelihood

is less than the maximum likelihood minus 2.71. The value of ELIV at which the profile

likelihood equals the maximum likelihood minus 2.71 is then found by repeatedly bisecting

the resulting interval.

Because the hard cutoff parameter has a negligible effect on the low-energy bins in the

97

Figure 7.3: Two times the natural logarithm of the profile likelihood defined in Equation
(7.4), plotted as a function of ELIV for the Crab Nebula.

98

101102103Cutoff energy (TeV)87868584832 ln LBest fit95.0% lower limitanalysis, their expectations are the effectively the same under the null and alternative hy-

potheses. As a result, they nearly cancel in the likelihood ratio. The limit is therefore driven

by the content of the high-energy bins, whose expected contents can change considerably

between hypotheses.

7.2 Results and interpretation

Figure 7.4 shows the contents of the reconstructed-energy bins using 2HWC J1825−134 data.

It also shows the expected background distribution as well as the expected counts under the

signal-plus-background model in which the hard cutoff is set to its 95%-confidence lower

limit. The limit derived from these data, prior to degrading for systematics, is 259 TeV. This

value lies in the second to last bin in the data. These bin contents and expectations are also

shown in Table 7.6.

Table 7.3 contains the p-values and 95%-confidence lower limits on ELIV for each of the
seven sources. p = 1 when the likelihood is maximized at ELIV = ∞. The lowest p-value of
0.082 occurs for 2HWC J1844−032; after correcting for seven independent statistical tests,

this corresponds to a p-value of 0.450. We therefore conclude that we have not observed

statistically significant evidence of an LIV-induced cutoff.

The strongest systematics-degraded limit of 209 TeV is given by 2HWC J1825−134. By

far the dominant source of systematic uncertainty is the QE, which allows the limit to vary

between 209 and 327 TeV. The QE systematic dominates for the other sources as well. In-

creasing the PMT QE in the MC increases the number of PMT hits expected per event,

thereby increasing the expected event rate at high estimated energy and making the fixed

data resemble a deficit compared to the null-hypothesis model. Conversely, a lower QE cor-

99

Figure 7.4: Expected and observed contents of the reconstructed-energy bins from the 2HWC
J1825−134 limit-setting procedure. The orange points show the expected distribution under
the signal-plus-background model in which the LIV cutoff energy is set to infinity and the
other parameters are optimized. The green points show the expected distribution under the
signal-plus-background model with the hard cutoff set to its 95%-confidence lower limit and
the parameters optimized. The red points indicate the expectations from background only.

100

100101102Reconstructed energy (TeV)102101100101102103104105106Number of eventsDataSignal + background (no LIV)Signal + background (95%-limit LIV)Background onlyresponds to a lower expected high-reconstructed-energy event rate, making the data counts

high compared to their model expectations and reducing the possibility of a cutoff. This

is exacerbated by the large QE uncertainties: as described in Section 4.6, the A, B, and D

PMTs’ QE is allowed to vary between 44 and 66%, and the C PMTs’ QE is varied between

64 and 96%.

As an additional systematic check, the Crab limit was computed modeling the Crab

spectrum instead as an exponentially cut off power law. The result was 137 TeV, somewhat

lower than the log-parabola-derived limit. Incorporating this value in the systematics for

the Crab limit inflates the overall systematic by 4 TeV.

Equation (1.18) can be used to compute 95%-confidence lower limits on EQG from the

ELIV limits given in Table 7.3. Table 7.4 shows the resulting limits. These substantially

exceed the previous best limits on EQG. An analysis [9] of HEGRA data set EQG lower
limits of 1.5 × 1029 eV for n = 1 and 2.8 × 1021 eV for n = 2. The limits from all seven

HAWC sources are more constraining than the HEGRA limits. The strongest limits are set
by 2HWC J1825−134: 8.74 × 1030 eV for n = 1 and 4.27 × 1022 eV for n = 2.

These limits exceed the previous best photon-decay-based LIV limits by more than an

order of magnitude. Limits set by searching for vacuum birefringence (such as [64]) exceed

these limits by several orders of magnitude; however, they are not applicable to all possible

photon-sector LIV scenarios [9]. With respect to models in which birefringence is disallowed,

the limits presented here are the most constraining.

Table 7.5 shows higher-confidence limits on the hard-cutoff energy. These have been
degraded to reflect systematic uncertainties as well. 2HWC J1825−134 provides 3σ evidence

for γ emission continuing up to 141 TeV. A cutoff below 55 TeV can be ruled out at the 5σ

level. The highest existing 5σ measurement of γ-ray flux was performed by HEGRA at an

101

Source

Crab
2HWC J1825−134
J1839−057
2HWC J1844−032
2HWC J1908+063
2HWC J2019+367
2HWC J2031+415

Φ0[(cm2 · s · TeV)−1]
3.6 × 10−13
3.6 × 10−13
4.0 × 10−13
2.6 × 10−13
1.2 × 10−13
0.7 × 10−13
0.7 × 10−13

α

2.63
1.95
2.29
2.34
2.22
1.65
1.89

β

0.10

Eexp(TeV) ELIV(TeV) σ(◦)

134
131
9904
185
48
86

5160
792
104
105
1166
214
251

0.55
1.01
1.29
0.64
0.29
0.42

Table 7.2: Best-fit values of the spectral parameters in the LIV model. The σ parameter is
fixed for each source.

energy of 42 TeV [5]. The HAWC measurement of 2HWC J1825−134 therefore constitutes

the highest-energy 5σ measurement of cosmic γ rays and a discovery of γ-ray flux above

55 TeV.

102

Source

p

Crab
2HWC J1825−134
J1839−057
2HWC J1844−032
2HWC J1908+063
2HWC J2019+367
2HWC J2031+415

0.412
0.998
0.997
0.082
1.000
0.791
0.709

Nom. Syst. err. Low CS High CS Low QE High QE Low
151 +44/−24
259 +68/−50
79 +23/−12
81 +23/−20
214 +60/−46
125 +30/−17
150 +36/−24

thresh.
139
245
71
75
208
124
148

131
211
70
63
168
108
126

154
263
82
81
217
128
153

145
259
83
88
212
123
152

194
327
101
103
273
155
186

No BP Nom.

degraded

162
266
82
75
222
129
153

127
209
67
61
168
108
126

Table 7.3: p-values and 95%-confidence limits on ELIV. Limits are expressed in TeV. The limit using the nominal MC is provided
along with a systematic error obtained by adding in quadrature the deviations of the limits calculated using the systematics
MCs. The systematics MCs are those having a lower than nominal charge smearing, higher charge smearing, lower QE, higher
QE, lower PMT threshold, and no broad-pulse effect. The systematics-degraded limit is then computed as the nominal limit
minus its lower error. See Section 4.6 for more details on the systematics MCs.

Source

Crab
2HWC J1825−134
J1839−057
2HWC J1844−032
2HWC J1908+063
2HWC J2019+367
2HWC J2031+415

n = 1

1.96 × 1030
8.74 × 1030
2.88 × 1029
2.17 × 1029
4.54 × 1030
1.21 × 1030
1.92 × 1030

n = 2

1.58 × 1022
4.27 × 1022
4.39 × 1021
3.64 × 1021
2.76 × 1022
1.14 × 1022
1.55 × 1022

Table 7.4: 95%-confidence lower limits on the QG energy EQG for n values of 1 and 2. All limits are expressed in eV. These
limits can be computed from the ones given in the nominal degraded column of Table 7.3 using Equation (1.18).

103

Source
Crab
2HWC J1825−134
J1839−057
2HWC J1844−032
2HWC J1908+063
2HWC J2019+367
2HWC J2031+415

95% 3σ 4σ 5σ
127
48
55
209
45
67
62
41
53
168
46
108
126
26

64
141
58
54
118
71
42

54
75
51
48
78
55
32

Table 7.5: Systematics-degraded limits on ELIV at various confidences. Limits are expressed
in TeV. A Zσ limit is defined to have confidence F (Z) where F is the standard-normal
cumulative distribution. As a result, the 3σ, 4σ, and 5σ limits correspond to confidences of
99.87%, 99.9968%, and 99.999971% respectively.

E (TeV) Obs.
225124
76029
19175
6271
1733
633
212
43
6
5
5
2
0
1

1.33
2.37
4.22
7.50
13.34
23.71
42.17
74.99
110.07
133.35
161.56
195.73
237.14
287.30

B

209268.32
66813.93
14233.04
3780.80
759.90
196.44
57.06
13.31
2.49
1.92
1.28
0.64
0.36
0.20

S + B (limit) S + B (no LIV)
225211.97
76502.31
19154.73
6306.56
1734.13
591.83
199.08
54.02
8.10
5.68
3.48
1.72
0.86
0.37

225468.61
76423.32
19082.14
6298.96
1761.70
623.18
204.23
25.53
2.53
1.92
1.28
0.64
0.36
0.20

Table 7.6: Observed and expected bin contents from the 2HWC J1825−134 analysis, as
plotted in Figure 7.4. Columns are central reconstructed bin energy, observed count, back-
ground expectation, signal-plus-background expectation with ELIV set to its 95% lower limit
of 259 TeV, and signal-plus-background expectation with ELIV set to infinity.

104

Chapter 8

Conclusions and future prospects

In addition to constraining models of LIV, the results presented in this thesis constitute the

first measurement of γ-ray flux above 100 TeV. This measurement suggests the existence of a

wealth of new science that can be carried out in this new energy window. It also demonstrates

the currently unparalleled ability of the HAWC experiment to carry out this science. As the

only γ-ray observatory able to resolve these energies, HAWC will be uniquely poised to carry

out measurements of γ-ray sources in this energy band. Several such analyses are already

underway and will be published within a year of this writing. Furthermore, the lower limit

on ELIV of 209 TeV derived from 2HWC J1825-134 implies that with 95% confidence, there

exist photons of this energy in the Universe and they can propagate on Galactic distance

scales.

The HAWC detector is currently undergoing a major enhancement. Smaller water tanks

called outriggers, each housing a single PMT, have been deployed and are currently taking

data. Once the development of software for calibration, reconstruction, and analysis of

outrigger data is complete, the background-rejection and reconstruction capabilities of the

experiment will improve further. A future iteration of the analysis presented here might gain

substantial sensitivity from the incorporation of outrigger data.

Another avenue for improvement is the reduction of systematic errors. The systematic

uncertainties on the MC described in Section 4.6 are deliberately conservative and reflect

105

the relatively wide ranges considered plausible for the parameters in question at the time

when those ranges were determined. A new iteration of HAWC’s MC and accompanying

systematics simulations are currently underway, and ongoing studies of these parameters

generally tend to suggest that some of their uncertainty ranges can be significantly reduced

although this work is not yet complete.

If a future version of this analysis is able to rule out an LIV-induced cutoff at even higher

energies, the attenuating effect of the cosmic microwave background (CMB) may need to

be taken into account. Figure 8.1 shows the mean free path of photons subject to CMB

interactions as a function of energy. At 259 TeV, the highest limit set in this thesis, the
mean free path is approximately 60 kpc. The distance to 2HWC J1825−134, the source

from which the most constraining limit is derived, is 4 kpc, and so only 6% of photons ought

to be lost to CMB interactions. However, Figure 8.1 shows that the mean free path drops

steeply at these energies. This suggests that a future analysis might need to take into account

the CMB interaction in modeling the sources’ energy spectra.

While the work presented in this thesis on the topic of LIV has ruled out certain physics

models, the work on energy reconstruction will also provide a new lens through which HAWC

can explore the high-energy universe. For this reason and others, the near future will be a

time of great scientific productivity for the HAWC experiment.

106

Figure 8.1: Mean free path of photons subject to the CMB [6].

107

101510161017Energy (eV)102101100Mean free path (Mpc)APPENDIX

108

Appendix

Data/MC agreement of neural-net

input variables

In Section 5.2, the NN Energy Estimator was shown to outperform the other available

HAWC γ-ray energy-reconstruction algorithms in MC. In order for the same performance to

be achieved in data, the 16-dimensional joint distribution of the 15 inputs and the output

must be correctly modeled in the MC. This is difficult to verify because of the curse of

dimensionality. A first-order check and necessary condition is to see that the one-dimensional

distributions of input variables are consistent between data and MC.

Of course the data also contain an admixture of cut-passing hadrons. Instead of modeling

these with MC, their distributions can be built from data and subtracted off. To this end,

two different locations in the sky must be considered: the location of the Crab Nebula
(the γ-source used for this comparison) and a location 4◦ away from the Crab at the same

declination. The cuts used for this analysis are the hadron-rejection cuts described in Chapter
5, the fiducial cut (Equation (5.2)), and a zenith cut at 40◦. The MC is weighted to resemble

a Crab-like spectrum with spectral index 2.63 and normalized to the data. The latter is

done to deconvolve the question of the data/MC agreement of the overall event rate, which

is not a function of the energy-estimation algorithm, from the question of the probability

distributions obeyed by the NN inputs. This analysis was performed using the HAWC-250

109

Figure .1: Data/MC comparison of fraction hit. p = 0.0997.

MC simulation that was current as of March, 2016, and used in the first publication of

HAWC’s measurement of the Crab Nebula energy spectrum, published in 2017.

The resulting distributions are shown in the figures in this appendix. To quantify any

disagreement in the distributions, a Kolmogorov-Smirnov (KS) test is performed for each

variable. The p-values from these tests are displayed in the captions of the respective figures.

Three input variables show a discrepancy at the 3σ level: the fraction of tanks hit and

the outermost two annuli. The uncertainty in the modeling of these variables and the

resulting uncertainty in the the mixing matrix of the NN energy variable are accounted

for by performing spectral analyses using a variety of different simulations as described in

Section 4.6.

110

Number of events10210DataMonte CarloFraction of channels hit0.30.40.50.60.70.80.9Data/MC00.511.52Figure .2: Data/MC comparison of fraction of tanks hit. p = 1.10 × 10−5.

111

Number of events110210DataMonte CarloFraction of tanks hit0.40.50.60.70.80.91Data/MC00.511.52Figure .3: Data/MC comparison of log of core-fit amplitude. p = 0.0538.

112

Number of events1-10110210310DataMonte Carlo(core amplitude)10log44.555.56Data/MC00.511.52Figure .4: Data/MC comparison of core distance from array center. p = 0.0538.

113

Number of events110210DataMonte CarloCore distance (m)20406080100Data/MC00.511.52Figure .5: Data/MC comparison of cosine of zenith angle. p = 0.715.

114

Number of events10210310DataMonte Carlocos(zenith angle)0.80.850.90.951Data/MC00.511.52Figure .6: Data/MC comparison of

fraction of charge contained in annulus 0 . p = 0.0762.

115

Number of events1-10110210DataMonte CarloFraction of charge between 0 and 10 m of core00.10.20.30.40.50.60.70.80.9Data/MC00.511.52Figure .7: Data/MC comparison of

fraction of charge contained in annulus 1 . p = 0.116.

116

Number of events1-10110210310DataMonte CarloFraction of charge between 10 and 20 m of core00.10.20.30.40.50.60.70.8Data/MC00.511.52Figure .8: Data/MC comparison of

fraction of charge contained in annulus 2 . p = 0.550.

117

Number of events1-10110210310DataMonte CarloFraction of charge between 20 and 30 m of core0.10.20.30.40.50.60.7Data/MC00.511.52Figure .9: Data/MC comparison of

fraction of charge contained in annulus 3 . p = 0.929.

118

Number of events2-101-10110210310DataMonte CarloFraction of charge between 30 and 40 m of core0.10.20.30.40.50.6Data/MC00.511.52Figure .10: Data/MC comparison of

fraction of charge contained in annulus 4 . p = 0.165.

119

Number of events110210310DataMonte CarloFraction of charge between 40 and 50 m of core00.10.20.30.40.50.6Data/MC00.511.52Figure .11: Data/MC comparison of fraction of charge contained in annulus 5 . p = 0.0118.

120

Number of events2-101-10110210310DataMonte CarloFraction of charge between 50 and 60 m of core00.050.10.150.20.250.30.350.40.45Data/MC00.511.52Figure .12: Data/MC comparison of

fraction of charge contained in annulus 6 . p = 0.11.

121

Number of events1-10110210310DataMonte CarloFraction of charge between 60 and 70 m of core00.050.10.150.20.250.30.350.40.45Data/MC00.511.52Figure .13: Data/MC comparison of fraction of charge contained in annulus 7 . p = 0.0196.

122

Number of events1-10110210310DataMonte CarloFraction of charge between 70 and 80 m of core00.050.10.150.20.250.30.350.4Data/MC00.511.52Figure .14: Data/MC comparison of
1.34 × 10−7.

fraction of charge contained in annulus 8 . p =

123

Number of events2-101-10110210310DataMonte CarloFraction of charge between 80 and 90 m of core00.050.10.150.20.250.3Data/MC00.511.52Figure .15: Data/MC comparison of fraction of charge contained in annulus 9. p = 4.54 ×
10−15.

124

Number of events2-101-10110210310DataMonte CarloFraction of charge more than 90 m from core00.10.20.30.40.50.6Data/MC00.511.52BIBLIOGRAPHY

125

BIBLIOGRAPHY

[1] Z. Hampel-Arias. Cosmic Ray Observations at the TeV Scale with the HAWC Observa-

tory. PhD thesis, The University of Wisconsin - Madison, 2017.

[2] A. U. Abeysekara et al. Observation of the Crab Nebula with the HAWC Gamma-Ray

Observatory. The Astrophysical Journal, 843(1):39, 2017.

[3] A. U. Abeysekara et al. Data Acquisition Architecture and Online Processing System

for the HAWC gamma-ray observatory. Nucl. Instrum. Meth., A888:138–146, 2018.

[4] M. Holler, D. Berge, C. van Eldik, J.-P. Lenain, V. Marandon, T. Murach, M. de
Naurois, R. D. Parsons, H. Prokoph, D. Zaborov, and for the H. E. S. S. collaboration.
Observations of the Crab Nebula with H.E.S.S. Phase II. arXiv e-prints, September
2015.

[5] F. Aharonian et al. The Crab nebula and pulsar between 500-GeV and 80-TeV. Observa-
tions with the HEGRA stereoscopic air Cerenkov telescopes. Astrophys. J., 614:897–913,
2004.

[6] M. Spurio. Particles and Astrophysics: A Multi-Messenger Approach. Astronomy and

Astrophysics Library. Springer International Publishing, 2014.

[7] L.D. Landau and E.F. Lifshitz. The Classical Theory of Fields. Course of theoretical

physics. Butterworth-Heinemann, 1975.

[8] Sidney Coleman and Sheldon L. Glashow. High-energy tests of Lorentz invariance. Phys.

Rev. D, 59:116008, Apr 1999.

[9] H. Mart´ınez-Huerta and A. P´erez-Lorenzana. Restrictions from Lorentz invariance vi-

olation on cosmic ray propagation. Phys. Rev. D, 95:063001, Mar 2017.

[10] D. Colladay and V. Alan Kosteleck´y. Lorentz-violating extension of the standard model.

Phys. Rev. D, 58:116002, Oct 1998.

[11] James Chadwick. The existence of a neutron. Proceedings of the Royal Society of London

Series A, 136, Jun 1932.

[12] W. Baade and F. Zwicky. Cosmic Rays from Super-Novae. Proceedings of the National

Academy of Sciences, 20(5):259–263, 1934.

[13] A. HEWISH, S. J. BELL, J. D. H. PILKINGTON, P. F. SCOTT, and R. A. COLLINS.

Observation of a Rapidly Pulsating Radio Source. Nature, 217(5130):709–713, 1968.

126

[14] C. Grupen. Astroparticle Physics. SpringerLink: Springer e-Books. Springer, 2005.

[15] R. N. Manchester, G. B. Hobbs, A. Teoh, and M. Hobbs. The Australia Telescope
National Facility Pulsar Catalogue. The Astronomical Journal, 129:1993–2006, April
2005.

[16] International Astronomical Union. Symposium, F. Camilo, and B.M. Gaensler. Young
Neutron Stars and Their Environments: Proceedings of the 218th Symposium of the
International Astronomical Union Held During the Iau General Assembly XXV, Sydney,
Australia, 14-17 July 2003. IAU Symposium Series. Astronmical Society of the Pacific,
2004.

[17] Bryan M. Gaensler and Patrick O. Slane. The evolution and structure of pulsar wind

nebulae. Ann. Rev. Astron. Astrophys., 44:17–47, 2006.

[18] R.N. Manchester and J.H. Taylor. Pulsars. A Series of books in astronomy and astro-

physics. W. H. Freeman, 1977.

[19] Margaret A. Livingstone, Victoria M. Kaspi, and Fotis P. Gavriil. Long-Term Phase-
coherent X-Ray Timing of PSR B0540–69. The Astrophysical Journal, 633(2):1095,
2005.

[20] J. M. Migliazzo, B. M. Gaensler, D. C. Backer, B. W. Stappers, E. van der Swaluw,
and R. G. Strom. Proper-Motion Measurements of Pulsar B1951+32 in the Supernova
Remnant CTB 80. The Astrophysical Journal Letters, 567(2):L141, 2002.

[21] T.K. Gaisser. Cosmic Rays and Particle Physics. Cambridge University Press, 1990.

[22] Progress in Elementary Particle and Cosmic Ray Physics. Number v. 3 in Series in

physics. North Holland Publishing Company, 1956.

[23] Koichi Kamata and Jun Nishimura. The Lateral and the Angular Structure Functions

of Electron Showers. Progress of Theoretical Physics Supplement, 6:93–155, 1958.

[24] Joshua Wood. An All-Sky Search for Bursts of Very High Energy Gamma Rays with

HAWC. PhD thesis, Maryland U., 2016.

[25] I. G. Wisher. Real-time Transient Monitoring With the HAWC Detector: Design and

Performance. PhD thesis, The University of Wisconsin - Madison, 2016.

[26] M.J. Berger, J.S. Coursey, M.A. Zucker, and J. Chang. Stopping-Power & Range
https://www.nist.gov/pml/

Tables for Electrons, Protons, and Helium Ions.
stopping-power-range-tables-electrons-protons-and-helium-ions.

[27] ZeroMQ home page. http://zeromq.org/.

127

[28] cJSON GitHub repository. https://github.com/DaveGamble/cJSON.

[29] M. Tavani et al. Discovery of Powerful Gamma-Ray Flares from the Crab Nebula.

Science, 331:736, February 2011.

[30] A. A. Abdo et al. Gamma-Ray Flares from the Crab Nebula. Science, 331(6018):739–

742, 2011.

[31] G. Aielli et al. Enhanced TeV gamma ray flux from the Crab Nebula observed, 09 2010.

[32] Benoit Cerutti et al. Extreme Particle Acceleration in Magnetic Reconnection Layers:
Application to the Gamma-Ray Flares in the Crab Nebula. The Astrophysical Journal,
746(2):148, 2012.

[33] A. U. Abeysekara et al. Sensitivity of the high altitude water Cherenkov detector to

sources of multi-TeV gamma rays. Astroparticle Physics, 50:26–32, December 2013.

[34] S. Ansoldi et al. Teraelectronvolt pulsed emission from the Crab pulsar detected by

MAGIC. Astron. Astrophys., 585:A133, 2016.

[35] Thomas Brian Humensky. VERITAS Studies of the Supernova Remnants Cas A and

IC 443. AIP Conf. Proc., 1085:357–360, 2009.

[36] M. Ackermann et al. Detection of the Characteristic Pion-Decay Signature in Supernova

Remnants. Science, 339(6121):807–811, 2013.

[37] A. A. Abdo et al. Milagro Observations of Multi-TeV Emission from Galactic Sources
in the Fermi Bright Source List. The Astrophysical Journal Letters, 700(2):L127, 2009.

[38] Dieter Horns et al. Hardening of TeV gamma spectrum of AGNs in galaxy clusters by

conversions of photons into axion-like particles. Phys. Rev., D86:075024, 2012.

[39] H. Krawczynski et al. Multiwavelength Observations of Strong Flares from the TeV

Blazar 1ES 1959+650. The Astrophysical Journal, 601(1):151, 2004.

[40] Francis Halzen and Dan Hooper. High energy neutrinos from the TeV Blazar 1ES

1959+650. Astroparticle Physics, 23(6):537 – 542, 2005.

[41] A. Reimer et al. Neutrino Emission in the Hadronic Synchrotron Mirror Model: The
“Orphan” TeV Flare from 1ES 1959+650. The Astrophysical Journal, 630(1):186, 2005.

[42] Pat Scott and et al. Direct Constraints on Minimal Supersymmetry from Fermi-LAT

Observations of the Dwarf Galaxy Segue 1. JCAP, 1001:031, 2010.

[43] E. Aliu et al. VERITAS deep observations of the dwarf spheroidal galaxy Segue 1. Phys.

Rev. D, 85:062001, Mar 2012.

128

[44] J. Aleksi´c et al. Segue 1: the best dark matter candidate dwarf galaxy surveyed by
In Proceedings, 32nd International Cosmic Ray Conference (ICRC 2011):

MAGIC.
Beijing, China, August 11-18, 2011, volume 5, pages 149–152, 2011.

[45] M. et al. Ackermann. Fermi-LAT Observations of the Gamma-Ray Burst GRB 130427A.

Science, 343(6166):42–47, 2014.

[46] R. Atkins et al. The High-Energy Gamma-Ray Fluence and Energy Spectrum of GRB
970417a from Observations with Milagrito. The Astrophysical Journal, 583(2):824, 2003.

[47] A. A. Abdo et al. Milagro Constraints on Very High Energy Emission from Short-

Duration Gamma-Ray Bursts. The Astrophysical Journal, 666(1):361, 2007.

[48] R. C. Gilmore, P. Madau, J. R. Primack, R. S. Somerville, and F. Haardt. GeV gamma-
ray attenuation and the high-redshift UV background. Monthly Notices of the Royal
Astronomical Society, 399:1694–1708, November 2009.

[49] M. A. S´anchez-Conde, D. Paneque, E. Bloom, F. Prada, and A. Dom´ınguez. Hints of
the existence of axionlike particles from the gamma-ray spectra of cosmological sources.
Physical Review D, 79(12):123511, June 2009.

[50] R. Brun and F. Rademakers. ROOT: An object oriented data analysis framework. Nucl.

Instrum. Meth., A389:81–86, 1997.

[51] Hao Zhou. Search for TeV Gamma-Ray Sources in the Galactic Plane with the HAWC

Observatory. PhD thesis, Michigan Technological University, 2015.

[52] D. Heck, J. Knapp, J. N. Capdevielle, G. Schatz, and T. Thouw. CORSIKA: a Monte

Carlo code to simulate extensive air showers. February 1998.

[53] S. et al. Agostinelli. GEANT4—a simulation toolkit. Nuclear Instruments and Methods

in Physics Research A, 506:250–303, July 2003.

[54] Andreas Hoecker, Peter Speckmayer, Joerg Stelzer, Jan Therhaag, Eckhard von Toerne,
and Helge Voss. TMVA: Toolkit for Multivariate Data Analysis. PoS, ACAT:040, 2007.

[55] K. M. G´orski, E. Hivon, A. J. Banday, B. D. Wandelt, F. K. Hansen, M. Reinecke, and
M. Bartelmann. HEALPix: A Framework for High-Resolution Discretization and Fast
Analysis of Data Distributed on the Sphere. The Astrophysical Journal, 622:759–771,
April 2005.

[56] D. C. Wells, E. W. Greisen, and R. H. Harten. FITS - a Flexible Image Transport

System. Astronomy and Astrophysics, 44:363, June 1981.

129

[57] R. Atkins, W. Benbow, D. Berley, E. Blaufuss, J. Bussons, D.G. Coyne, R.S. Delay,
T. DeYoung, B.L. Dingus, D.E. Dorfan, R.W. Ellsworth, and A. Falcone. Observation
of TeV Gamma Rays from the Crab Nebula with Milagro Using a New Background
Rejection Technique. The Astrophysical Journal, 595(2):803–811, 2003.

[58] D. E. Alexandreas, D. Berley, S. Biller, G. M. Dion, J. A. Goodman, T. J. Haines,
C. M. Hoffman, E. Horch, X.-Q. Lu, C. Sinnis, G. B. Yodh, and W. Zhang. Point source
search techniques in ultra high energy gamma ray astronomy. Nuclear Instruments and
Methods in Physics Research A, 328:570–577, May 1993.

[59] S. S. Wilks. The Large-Sample Distribution of the Likelihood Ratio for Testing Com-

posite Hypotheses. Ann. Math. Statist., 9(1):60–62, 03 1938.

[60] O. C. de Jager, A. K. Harding, P. F. Michelson, H. I. Nel, P. L. Nolan, P. Sreekumar,
and D. J. Thompson. Gamma-Ray Observations of the Crab Nebula: A Study of the
Synchro-Compton Spectrum. The Astrophysical Journal, 457:253, January 1996.

[61] D. Horns and F. A. Aharonian. The Crab Nebula: Linking MeV Synchrotron and 50
TeV Inverse Compton Photons. In V. Schoenfelder, G. Lichti, and C. Winkler, editors,
5th INTEGRAL Workshop on the INTEGRAL Universe, volume 552 of ESA Special
Publication, page 439, October 2004.

[62] A. U. et al. Abeysekara. The 2HWC HAWC Observatory Gamma-Ray Catalog. The

Astrophysical Journal, 843:40, July 2017.

[63] G. Cowan. Statistical Data Analysis. Oxford science publications. Clarendon Press,

1998.

[64] D. Gotz, S. Covino, A. Fernandez-Soto, P. Laurent, and Z . Bosnjak. The polarized

Gamma-Ray Burst GRB 061122. Mon. Not. Roy. Astron. Soc., 431:3550, 2013.

130