THE UNIVERSITY OF TOKYOGRADUATE SCHOOL OF SCIENCES

Thesis

TIME-INTEGRATED AND TIME-DEPENDENTDIRECTIONAL SEARCH FOR HIGH ENERGYASTROPHYSICAL NEUTRINO POINT SOURCES IN
SUPER-KAMIOKANDE(スーパーカミオカンデにおける高エネルギー天体ニュートリノ点
源の時間および方向の探索)

by

XUBIN WANG

Submitted in partial fulfillment of therequirements for the degree ofMaster of Science2022© Copyright by
XUBIN WANG
2022

ACKNOWLEDGMENTS
I sincerely appreciate the support provided by Prof. Kimihiro Okumura in both study
and life. His guidance which enlightens me through the whole research process is theindispensable key in completing this project. I am also truly grateful that he offeredme the opportunity for such a research so that I did not get lost in this miserabletime.
I also would like to express my appreciation to Dr. Mathieu Lamoureux, for the
analysis libraries he has established for astrophysical neutrino field, and his patienceanswering to my questions. With his work and suggestions, I was able to carry on
the research with less diﬀiculties.
Prof. Katsuaki Asano has also gave me a lot of help in the past two years. As my
mentor, he not only helped me in the life on campus, but also inspired me with his
insightful understanding in the field of astrophysics.
Also I want to say thanks to Dr. Linyan Wan, as she answered the question about
the UPMU sample which confused me for a long time.
Members in ICRR and Kamioka have provided me great help in getting used tothe study of neutrino physics and life in Japan. I would especially thank Junjie Xia,
Seungho Han from RCCN and Jiahui Feng from Kyoto University.
And my sincere thank give to all members in SK group, for none of the analysis
work could be done without the the experiment setup and software development.
Last but not least, thank you my parents, my family, and my friends. Without
your support, I would never be able to make this far.

iii

TIME-INTEGRATED AND TIME-DEPENDENTDIRECTIONAL SEARCH FOR HIGH ENERGYASTROPHYSICAL NEUTRINO POINT SOURCES IN
SUPER-KAMIOKANDE(スーパーカミオカンデにおける高エネルギー天体ニュートリノ点
源の時間および方向の探索)
XUBIN WANG
ABSTRACT
Astrophysical high-energy neutrino point sources has been considered of extreme interest in the field of both astrophysics and particle physics. The mysterious neutral
lepton with supper tiny mass, effected by neither the magnetic field nor the cosmic
microwave background, is believed to be able to carry the information of the ultrahigh-energy celestial objects and phenomena, and even the cosmic accelerator thathas confused scientists for a long time. Very few is known to us about the astrophysical neutrinos or their sources up until now, but the effort has never been stopped
and remarkable progress has been achieved.
Using the data collected by the water Cherenkov detector Super-Kamiokande
during 2008 and 2018, we performed both a time-integrated and a time-dependent
search for the potential candidates of the high-energy astrophysical neutrino sources
with upward-going muon events via unbinned maximum likelihood method. No excesswith respect to the background is found in either search.
The most probable location to have a point source in the full sky is found at(198.3◦ , 15.1◦ ) in J2000 equatorial coordinate system with pre-trial p-value of 1.18 ×
10−4 and post-trial p-value of 79.2% in the time-integrated search. The upper limitsiv

of the neutrino fluence and flux are set to Φ = 2.06 × 102 [GeVcm−2 ] and φupper =4.53 × 10−7 cm−2 s−1 , respectively.
In the time-dependent search, the most significant cluster is at (118.0◦ , −62.7◦ ) on
2010.10.27 13:09:42 JST and half emission duration σ̂t = 9.25 × 103 s. The pre-trial
p-value is 2.24 × 10−4 and post-trial p-value is 33.4%. Neutrino fluence upper limit isΦ = 7.86 × 101 GeVcm−2 , with the emission time interval estimated by 1.85 × 104 s(2σ̂t ).
The algorithm and the library of the numerical analysis established in this study
can be applied on future relative researches on topics regarding the astrophysicalneutrinos.

v

CONTENTS

Acknowledgementsiii

Abstractiv

List of Tablesix

List of Figuresx

1 Introduction1

1.1

Neutrino . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1

1.2

High-Energy Astrophysical Neutrinos . . . . . . . . . . . . . . . . . .

2

1.2.1

Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2

1.2.2

Cosmic Accelerator . . . . . . . . . . . . . . . . . . . . . . . .

4

Potential Neutrino Source Candidates . . . . . . . . . . . . . . . . . .

5

1.3.1

Active Galactic Nucleus . . . . . . . . . . . . . . . . . . . . .

5

1.3.2

Gamma-Ray Burst . . . . . . . . . . . . . . . . . . . . . . . .

7

1.3.3

Supernova Remnant . . . . . . . . . . . . . . . . . . . . . . .

8

1.3.4

Magnetar . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

8Multi-Messenger Astronomy . . . . . . . . . . . . . . . . . . . . . . .

9

1.4.1

Cosmic Rays . . . . . . . . . . . . . . . . . . . . . . . . . . .

9

1.4.2

Gamma rays . . . . . . . . . . . . . . . . . . . . . . . . . . .

11

1.4.3

Gravitational Waves . . . . . . . . . . . . . . . . . . . . . . .

12Previous Related Search . . . . . . . . . . . . . . . . . . . . . . . . .

13

1.3

1.4

1.5

vi

2 Super-Kamiokande

14

2.1

Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

14

2.2

Cherenkov Radiation . . . . . . . . . . . . . . . . . . . . . . . . . . .

16

2.3

Photomultiplier . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

17

2.4

Electronics and Data Acquisition System . . . . . . . . . . . . . . . .

19

2.5

Water and Air System . . . . . . . . . . . . . . . . . . . . . . . . . .

22

3 Upward-Going Muon

23

3.1

SK Data Sample . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

23

3.2

Data Reduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

24

3.3

Event Reconstruction . . . . . . . . . . . . . . . . . . . . . . . . . . .

26

4 Event Simulation

28

4.1

Atmospheric Neutrino Flux . . . . . . . . . . . . . . . . . . . . . . .

28

4.2

Neutrino Interaction . . . . . . . . . . . . . . . . . . . . . . . . . . .

30

4.3

Detector Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . .

32

5 Analysis

34

5.1

Data Sample . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

34

5.2

Unbinned Maximum Likelihood Method . . . . . . . . . . . . . . . .

36

5.2.1

Likelihood . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

36

5.2.2

Variables . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

38

5.2.3

Probability Density Function . . . . . . . . . . . . . . . . . .

39

5.2.4

Expected Number of Background Events . . . . . . . . . . . .

46

5.2.5

Test Statistics . . . . . . . . . . . . . . . . . . . . . . . . . . .

47

All-Sky Search . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

50

5.3.1

50

5.3

Parameters Fitting and Likelihood Maximization . . . . . . .
vii

5.4

5.5

5.3.2

First Guess of Time Parameters . . . . . . . . . . . . . . . . .

51

5.3.3

Search Directions . . . . . . . . . . . . . . . . . . . . . . . . .

53

Pseudo Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . .

54

5.4.1

Toy Monte Carlo . . . . . . . . . . . . . . . . . . . . . . . . .

55

5.4.2

Search Algorithm Confirmation . . . . . . . . . . . . . . . . .

56

5.4.3

Sensitivity . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

60

Neutrino Fluence . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

61

6 Results and Discussion
6.1

68

Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

68

6.1.1

Time-Integrated Search . . . . . . . . . . . . . . . . . . . . .

68

6.1.2

Time-Dependent Search . . . . . . . . . . . . . . . . . . . . .

71

6.2

Potential Candidate Check . . . . . . . . . . . . . . . . . . . . . . . .

73

6.3

Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

75

7 Summary and Prospects

83

Bibliography

84

viii

LIST OF TABLES
2.1

Cherenkov energy thresholds in water for different particles. . . . . .

18

6.1

Table of the best fitted parameters for time-integrated search. . . . .

69

6.2

Table of the best fitted parameters of the top three largest T S fortime-integrated search. . . . . . . . . . . . . . . . . . . . . . . . . . .

6.3

72

Energy, distance from best fitted position, time interval from the best
fitted t0 , and event type of the events within 15◦ to the position of the
top three largest T S with in 7σ̂t around t̂0 . . . . . . . . . . . . . . . .

74

6.4

List of the IceCube alert. . . . . . . . . . . . . . . . . . . . . . . . . .

80

6.5

List of potential astrophysical neutrino source candidates[21]. The column from left is name, declination, right ascension, distance measured
in z (redshift) or kpc, observed gamma ray flux, energy threshold forflux integration, gamma ray emission energy spectrum index[60], andneutrino flux 90 CL upper limit at that direction for time-integrated
search (averaged for total live time of SK-IV). . . . . . . . . . . . . .

ix

81

LIST OF FIGURES
1.1

Top: Results of time-dependent analysis. Black dashed line represents
IceCube-170922A. A significant peak over 3σ was discovered during
IC86b, at around the beginning of 2015. Bottom: Skymap of p-value.
Black cross at the middle is TXS 0506+056. Figures are taken from [30].

61.2Schematic of radio-loud AGN[59] . . . . . . . . . . . . . . . . . . . .

7

1.3

Crab Nebula[51] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

8

1.4

Multimessenger observation of TXS 0506+056 via neutrino and gamma
rays by IceCube, MAGIC and Fermi. Figure taken from [31] . . . . .

10

1.5

Energy spectrum of CR [38] . . . . . . . . . . . . . . . . . . . . . . .

11

2.1

A schematic of SK.[44] . . . . . . . . . . . . . . . . . . . . . . . . . .

14

2.2

Frame structure which supports PMTs.[35] . . . . . . . . . . . . . . .

15

2.3

Cherenkov radiation light cone.[9] . . . . . . . . . . . . . . . . . . . .

16

2.4

Spectrum of Cherenkov radiation in water.[48] . . . . . . . . . . . . .

18

2.5

20-inch PMT.[35] . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

19

2.6

20-inch PMT quantum eﬀiciency distribution as a function of wavelength.[32] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

20

2.7

FRP case and acrylic cover for ID PMT.[35] . . . . . . . . . . . . . .

20

2.8

Timing chart of QTC in QBEE.[52] . . . . . . . . . . . . . . . . . . .

21

2.9

Block diagram of the QTC.[52] . . . . . . . . . . . . . . . . . . . . .

21

3.1

Schematic of SK samples.[62] . . . . . . . . . . . . . . . . . . . . . .

23

x

3.2

Expected atmospheric neutrino energy distribution for different samples.[8] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

24

4.1

Schematic the creation of atmospheric neutrinos.[47] . . . . . . . . .

29

4.2

Left:Calculated atmospheric neutrino fluxes averaged over all directions for different models. Right: Neutrino and anti-neutrino flavor
ratios for different models.[42] . . . . . . . . . . . . . . . . . . . . . .

4.3

Total charged current cross sections for neutrinos (top 4 lines) andantineutrinos (bottom 4 lines).[39] . . . . . . . . . . . . . . . . . . . .

5.1

36

AS (∆ϕ|Eµ ) for signal events at different reconstructed energy. Left:
P (∆ϕ|Eµ ). Right: P (φ) . . . . . . . . . . . . . . . . . . . . . . . . .

5.5

35

Angular resolution corresponding to 1σ for FC, PC and UPMU fromatmospheric neutrino simulation. . . . . . . . . . . . . . . . . . . . .

5.4

34

Time distribution of different types of events. Time is presented in themodified Julian day (mjd). . . . . . . . . . . . . . . . . . . . . . . . .

5.3

32

A sky map of different types of UPMU events. All the sky map in thisand next chapter are shown in J2000 equatorial coordinate. . . . . . .

5.2

30

40

AB (dec, ra) for background events used in time-dependent search. Left:
P (alt). Right: P (azi). . . . . . . . . . . . . . . . . . . . . . . . . . .

41

5.6

Local sidereal time (LST) distribution. . . . . . . . . . . . . . . . . .

41

5.7

AB (dec, ra) for background events used in time-integrated search. Left:
P (dec). Right: P (ra). . . . . . . . . . . . . . . . . . . . . . . . . . .

5.8

42

Probability P (Eν , Eµmin ) for a neutrino of energy Eν creating a muonof energy larger than Eµmin . . . . . . . . . . . . . . . . . . . . . . . .

xi

44

5.9

Earth shadow effect S(z, Eν ) for νµ (left) and ν̄µ (right) as a functionof log10 Eν and z = cos θ. . . . . . . . . . . . . . . . . . . . . . . . . .

45

5.10 Top: The probability that a neutrino of energy Eν creates a muon ofreconstructed energy Eµ (f (Eµ |Eν )).Bottom: The probability that amuon of reconstructed energy Eµ is created by a neutrino of energy Eν 45
5.11 Expected reconstructed energy distribution for background events andsignal events with different spectrum indices. . . . . . . . . . . . . . .

46

5.12 T S (t|t0 , σt ) for signal (left) and T B for background (right) . . . . . .

47

5.13 Top left: Sky map of expected number of background events within
the search cone. Top right: Sky map of number of events within the
search cone for the data. Bottom: Expected number of background
events as a function of declination. . . . . . . . . . . . . . . . . . . .

48

5.14 T S distribution for time-integrated search in a pseudo experiment withsolely background toy MC. . . . . . . . . . . . . . . . . . . . . . . . .

51

5.15 T S as a function of t0 with fixed σt in a pseudo experiment using toy
MC. Vertical dashed blue line represents the assumed emission center
time. Simulated with nS = 10 on top of the background toy MC,
σt = 105 s, at dec = 30◦ . . . . . . . . . . . . . . . . . . . . . . . . . . .

52

5.16 A schematic of how HEALPix divide the full sky into pixels.[3] . . . .

53

5.17 A sky map of T S for time-dependent search in pseudo experiment withsolely background toy MC. . . . . . . . . . . . . . . . . . . . . . . . .

54

5.18 Sky map of T S of time-dependent search for solely background 5.18aand nS = 7 at dec = 30◦ 5.18b, 0◦ 5.18c, −30◦ 5.18d, −60◦ 5.18e, −90◦ 5.18f.
Red cross represents the location of assumed source, and yellow plus
represents the direction of T Smax . . . . . . . . . . . . . . . . . . . . .

xii

57

5.19 T S as a function of t0 and σt for nS = 2 σt = 0s 5.19a, nS = 3 σt = 10s
5.19b, nS = 5 σt = 100s 5.19c, nS = 7 σt = 1, 000s 5.19d, nS = 10 σt =104 s 5.19e, nS = 13 σt = 105 s 5.19f, nS = 17 σt = 106 s 5.19g. Verticaland horizontal dashed green line represents the assumed log10 σt and
t0 , respectively. Yellow cross represents the best fitted log10 σt and t0 .

58

5.20 Sky map of T S of time-integrated search for solely background 5.20aand nS = 20 dec = 30◦ 5.20b, nS = 20 dec = 0◦ 5.20c, nS = 25 dec =−30◦ 5.20d, nS = 30 dec = −60◦ 5.20e, nS = 40 dec = −90◦ 5.20f.
Red cross represents the location of assumed source, and yellow plus
represents the direction of T Smax . . . . . . . . . . . . . . . . . . . . .

59

5.21 T Smax distribution for solely background events. Vertical red solid line
indicates the 90% CL threshold. . . . . . . . . . . . . . . . . . . . . .

60

5.22 T Smax distribution with different nS . Vertical red solid line indicatesthe background 90% CL threshold. Vertical blue and green line represent the 50% CL for nS = 2 and nS = 3. . . . . . . . . . . . . . . . .

61

5.23 Sensitivity for time-integrated search in terms of nS (top) and Φ (bottom). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

62

5.24 Sensitivity for time-dependent search in terms of nS (top) and Φ (bottom). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

63

5.25 ψ(z) for νµ (left) and ν̄µ (right) . . . . . . . . . . . . . . . . . . . . .

65

5.26 Schematic of calculating effective area.[55] . . . . . . . . . . . . . . .

65

5.27 Effective area as a function of z. . . . . . . . . . . . . . . . . . . . . .

66

5.28 Product of ψ(dec) and Aef f (dec) as a function of declination. . . . . .

67

xiii

6.1

Sky map of T S for time-integrated search. Location of T Smax ismarked with red triangle. Positions with n̂S = 0 has been excluded
from the sky map, which is represented by the white color region. . .

6.2

68

Blue filled histogram: T S distribution for time-integrated search fromall searching directions from data. Red dashed line histogram: T S distribution from all searching directions obtained from 50 trials (roughly
2 million times searches) with solely background toy MC. The peakaround 0 represents the situation where n̂S = 0. The result from thedata basically agrees with the simulation for background. . . . . . . .

6.3

The histogram of post-trial T Smax distribution for time-integrated
search obtained from 1,000 trials with solely background toy MC. . .

6.4

70

Histogram of pre-trial T S distribution for time-integrated search bysimulating 106 times at the position of T Smax . . . . . . . . . . . . . .

6.5

69

70

Stopping (red plus), through-going (green dot) and showering (bluecross) UPMU events detected around the position of T Smax . Greyregion encircled by the black circle is the 15◦ search cone. . . . . . . .

6.6L(nS |γ̂) as a function of nS with γ fixed to the best fitted value γ̂.
The vertical solid blue line represent the upper limit for 90% CL. . .

6.7

71

72

Sky map of T S for time-dependent search. The locations of the top
three maximum T Ss are marked with red triangle, plus and cross,respectively. Positions with n̂S = 0 has been excluded from the skymap, which is represented by the white color region. . . . . . . . . . .

xiv

73

6.8

Blue filled histogram: T S distribution for time-integrated search fromall searching directions from data. Red dashed line histogram: T S
distribution from all searching directions obtained from 50 trials withsolely background toy MC. The peak around 0 represents the situations
where n̂S = 0 and σ̂t = (Tend − Tstart )/2, the maximum value set forσt . The result from the data basically agrees with the simulation forbackground. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6.9

74

The histogram of post-trial T Smax distribution for time-dependent
search obtained from 1,000 trials with solely background toy MC. . .

75

6.10 Histograms of pre-trial T S distribution for time-dependent search of
1st (top), 2nd (middle) and 3rd (bottom). Each one is simulated 106
times at the corresponding direction. . . . . . . . . . . . . . . . . . .

76

6.11 L(nS |γ̂, t̂0 , σ̂t ) as a function of nS with γ, t0 and σt fixed to the best
fitted value γ̂, t̂0 and σ̂t , for 1st (top), 2nd (middle) and 3rd (bottom).
The vertical line represents the upper limit for 90% CL. . . . . . . . .

77

6.12 Space orientation (left column) and time distribution (right column) of
the events within 15◦ to the position of the 1st (top), 2nd (middle) and3rd (bottom) largest T S with in 7σ̂t around t̂0 . Black dashed circle
in the left figures is the boundary of the search cone. The red plus
indicates the best fitted source position. Blue cross is the directionof the UPMU event, with the circle around them being the angularresolution. The black dashed line in the right figures is the Gaussian
distribution with the best fitted t0 and σt as the parameter. . . . . .

xv

78

6.13 T S sky map for time-integrated (top) and time-dependent (bottom)
search, with the the green squares represent the events from IceCube
alert catalog (Table 6.4) and yellow diamonds represent the candidates
from list in Table 6.5, and the red dot is the location of TXS 0506+056. 79

xvi

1
CHAPTER 1
Introduction
1.1

NEUTRINO

In 1930, Pauli first postulated the existence of neutrinos to explain the momentum
conservation of beta decay. The electron anti neutrino was first observed by Reines
and Cowan in 1956 at a nuclear reactor. Using the neutrino beam produced by the
decay of charged pions, the existence of the muon-neutrino was confirmed by Danby in
1962. The third kind of neutrino, tau-neutrino, was discovered by Direct Observationof the Nu Tau (DONuT) Experiment in 2001.
In the standard model, neutrinos were first thought to be massless, neutral leptons.
Although they have been proved to have mass now, it is negligible when considering
the effect of gravity. As a result, the only interaction channel for them is weak interaction. The neutrinos detected on Earth are mainly generated via nuclear reactionsin the core of the sun and high energy cosmic rays interacting with nucleus in theatmosphere. Except for natural sources, nuclear reactors and accelerators can also
generate neutrinos and are widely used in the studies of neutrino properties.
The discovery of neutrino oscillation can date back to the Homestake Experimentin 1960s, which observed a huge difference between the measured electron-neutrino
flux and the one predicted by the standard solar model. This solar neutrino problem had remained unsolved for about thirty years before physicists realized that theproperties of neutrino that had been assumed might not be accurate. If the threekinds of neutrinos are massive and have mass differences, their flavor eigenstates will
be different from their mass eigenstates, and thus the flavor can oscillate to one thatdiffers from what the neutrino is generated with. Such a phenomenon is called neu-

2trino oscillation. In 1998, Super-Kamiokande Experiment (SK) found the evidenceof neutrino oscillation via the asymmetry between the upward-going and downwardgoing atmospheric neutrino. Then with the measurement of solar neutrinos at SK
and Sudbury Neutrino Observatory, solar neutrino problem was well explained.
The first observation of astrophysical neutrinos outside solar system was on February 23rd, 1987, from a supernova explosion in Large Magellanic Cloud 51.4kpc awayfrom the Earth. This supernova was numbered as SN1987A, and SK, Irvine-MichiganBrookhaven Experiment (IMB) and Baksan Underground Scintillation Telescope (BUST)
detected 11, 8 and 5 neutrinos, respectively. In the supernova model, 99% of theenergy of the core collapse supernova will be released via neutrino emission. Thisobservation not only examined the model, but also remarked the beginning of using
neutrino as a tool to study astronomy.

1.2

HIGH-ENERGY ASTROPHYSICAL NEUTRINOS

1.2.1

Overview

High-energy astrophysical neutrinos refer to the neutrinos from galactic or extragalactic distant objects with energy above 1 GeV. Unlike solar neutrinos or supernova
neutrinos, high-energy neutrinos are not produced directly by the celestial bodies.
Instead, they are produced as a result of the cosmic accelerator, which is a kind ofelectromagnetic accelerator. More details will be discussed in 1.2.2.
Protons accelerated by the cosmic accelerator can scatter off another accelerated
proton or a photon and produce mesons through the following channel:

p + γ → π+ + n(1.1)

p + γ → π0 + p

(1.2)

3
p + p → π±, π0(1.3)

The charged and neutral pions, together with protons and neutrons, can undergo
the following decays:

π + → µ+ νµ → e+ + νe + ν̄µ + νµ(1.4)

π 0 → 2γ

(1.5)

Among the produced particles, photons, neutrinos and neutrons are neutral andwill escape from accelerator. Additionally, the neutron can decay into a proton:n → p + e− + ν¯e(1.6)

Energy distribution of neutrinos generated in this way will follow a power spectrum. This will be introduced in details in Section 1.4.1.
Although cosmic rays (mainly consists of protons), gamma rays and neutrinos
are all produced in this process, neutrino has advantages compared with the othertwo. First, it is electromagnetic neutral, meaning that it will not be deflected by
the magnetic field existing in the galaxy and so that the direction information of thesource is not lost. Second, its cross section is extremely small so that the probability
that it scatters off another particle is tiny. Gamma rays can be blocked by the opaqueobjects, and cosmic rays will interact with the cosmic microwave background (CMB)
at high-energy region.

4
1.2.2

Cosmic Accelerator

The cosmic rays we have observed on Earth have energies that is much higher than the
expected thermal energies of cosmic rays emitted by astrophysical objects. The mostcommonly accepted explanation of such a phenomenon is that the cosmic rays areaccelerated by the electromagnetic fields known as the cosmic accelerator. Althoughmany models have been established, the mechanism of acceleration still remains unknown to us.
The mechanisms of acceleration can be different for different astrophysical objects, and the process may consist of more than one mechanism. Among the mechanisms, shock acceleration is widely included in the scenarios of many astrophysicalobjects[34]. In this model, the accelerating region contains two parts: upstream anddownstream. The upstream is a region full of charged particles with high velocity(mainly protons) generated by some certain kind of driven force such as the supernova explosion. In the downstream, on the other hand, the velocity of the charged
particle is smaller. The boundary between the upstream and downstream is called
the shock front. When a charged particle moves across the shock front, it can scatteron a magnetic irregularity and change the momentum direction. Then it will moveacross the shock front and scatter again. In the process of bouncing back and forth,the energy of the particle will increase or decrease, depending on the scattering angles, but in average the particle gains energy. After the energy of the particle is large
enough, it has the probability to escape from the shock.
Such kind of acceleration can take place both galactically and extra-galactically.
There is an energy limit which the particle trapped within the magnetic field in thegalaxy can reach. As a result, the cosmic rays with the highest energies are probably
generated from extra-galactical sources.

5
1.3

POTENTIAL NEUTRINO SOURCE CANDIDATES

Except for solar neutrinos and supernova neutrinos, no astrophysical neutrino event
had been observed before 2017. However, on 22 September 2017, IceCube detected
an ultra-high-energy neutrino event of approximately 290 TeV, labeled as IceCube170922A. Following this observation, IceCube performed a search for neutrino point
sources at the position of TXS 0506+056, but independent of IceCube-170922A[30].
In that search, they used the time-dependent search method which is similar to what
is used in this study, and discovered an excess of neutrino emission of a confidential
level above 3σ between 2014 and 2015 (Figure 1.1). This provided a strong evidence
that blazars could be the sources of high-energy astrophysical neutrinos.
There are, however, many potential neutrino source candidates that are considered
to generate neutrinos theoretically. In this part I will briefly introduce some of them.
1.3.1

Active Galactic Nucleus

An Active Galactic Nucleus (AGN) is super massive black hole (BH) with an accretion
disk located at the center of a distant galaxy.
The BH and accretion disk is surrounded by clouds of gas moving rapidly in thepotential of the BH and a torus or warped disk of gas and dust well outside theaccretion disk. Relativistic collimated jets are emitted perpendicular to the disk.
Figure 1.2 shows a schematic of AGN. Those AGN with the jet direction pointing tothe observer are called the blazar[59].
When matter in the accretion disk fall into the BH, it creates an accretion shock,leading to the first order Fermi acceleration of charged particles including protons.
These ultrarelativistic protons can collide with intensive photon fields and produce
high-energy neutrinos[57].

6

(a)

(b)

Figure 1.1: Top: Results of time-dependent analysis. Black dashedline represents IceCube-170922A. A significant peak over 3σ was discovered during IC86b, at around the beginning of 2015. Bottom: Skymapof p-value. Black cross at the middle is TXS 0506+056. Figures aretaken from [30].

7Figure 1.2: Schematic of radio-loud AGN[59]1.3.2

Gamma-Ray Burst

Gamma-ray burst (GRB) are extremely luminous electromagnetic radiation events
observed in distant galaxies.
The GRB is considered to be produced by the dissipation of an expanding fireball. The ultrarelativistic particles in the fireball are generated in the collapse of a
massive rotating star core or the merger of a neutron star with another neutron staror blackhole, where large amount of gravitational energy is released in a short timein a small region[49].
In this fireball model, protons can undergo Fermi acceleration to energies up to
102 0 eV[61]. Then high-energy neutrinos can be generated through the photomeson
production.

8
1.3.3

Supernova Remnant

Supernova Remnants (SNR) are expanding materials ejected from supernovae explosions bounded by shock waves. Figure 1.3 is one of the most famous supernova
remnant, Crab Nebula, which is also a plerion (pulsar wind nebula).

Figure 1.3: Crab Nebula[51]The possible acceleration mechanism SNR have is the diffusive shock acceleration[33]. When the ejecta meets the surrounding interstellar medium, it is possible
to produce neutrinos via pion production.
1.3.4

Magnetar

Magnetars are rotating neutron stars with strong magnetic fields (∼ 1015 G).
There are two main energy source to power the magnetars: the magnetic power
and the spin-down power[28]. For magnetars whose spin and magnetic momentum
are in opposite directions, spin-induced electric fields can accelerate the protons fromsurfaces. Then they will scatter off the photons generated in the decay of strong

9magnetic fields and produce pions, which will later decay into neutrinos.
The particles emitted by magnetars are not isotropic but will concentrate in acertain direction. The expected neutrino flux will be magnificent in that direction,
yet much lower in others. It is not very likely the direction happens to point at Earth.
However, if it does, it will be helpful for the study on astrophysical neutrinos andmagnetar structure.

1.4

MULTI-MESSENGER ASTRONOMY

After the detection of high-energy neutrino event IceCube-170922A in 2017, IceCube
sent out an alert to all gamma ray observatories. It was later found out that theblazar TXS 0506+56, whose position is coincident with the direction of the neutrino
event, was in a flaring state. From that time, a series of gamma-ray observation ina various of wavelength range followed, and high-energy gamma rays were confirmed
to be detected from that blazar (Figure 1.4)[31].
This kind of research, usually called multi-messenger research, has been proved
very useful, as various messengers can carry information of different process. Other
than neutrinos, cosmic rays, electromagnetic radiation and gravitational wave are also
the major messengers used in astronomy studies.
1.4.1

Cosmic Rays

Cosmic rays (CR) are high energy astrophysical particles mainly consist of protons(90%), nucleus (9%) and electrons (1%)[53]. I have briefly introduced their origins
and the acceleration mechanism in previous chapters so I will mainly introduce theirenergy spectrum in this chapter.
Energy distribution of cosmic rays follows a power spectrum, with a spectral index

10

Figure 1.4: Multimessenger observation of TXS 0506+056 via neutrino and gamma rays by IceCube, MAGIC and Fermi. Figure takenfrom [31]11
gamma, shown in Figure 1.5. There are two important structures: a knee, and anankle. The knee is at ∼ 101 5eV . Below the knee, gamma

2.7. The ankle is at∼ 101 9eV . Above that gamma 2.6. Between the knee and the ankle, gamma 3.1.

Figure 1.5: Energy spectrum of CR [38]The spectrum does not continue to infinity. At around 5 × 1019 , there is a rapid
steepening in the spectrum called the Greisen-Zatsepin-Kuzmin (GZK) cutoff[37, 63].
Protons with energies higher than this threshold undergo inelastic interactions withCMB:
p + γ → π0 + p
1.4.2

Gamma rays

Gamma rays can be generated both leptonically and hadronically.

(1.7)

12
Leptonic process refers to the synchrotron radiation of accelerated electrons propagating in strong magnetic fields. Since the synchrotron radiation will be suppressedfor particles with larger mass, heavier particles like protons does not have much contribution to this process.
The hadronic process is similar to neutrino creation, with only charged pion decay
being changed to neutral pion decay. Since neutrinos and gamma rays are produced
at the same time, and neither of them are deflected by the magnetic fields in thegalaxy, they can be used as the triggers of each other.
1.4.3

Gravitational Waves

Gravitational wave (GW) is a different kind of messenger in comparison with CR,
gamma rays and neutrinos. It is not particle, but the propagation of the change
in curvature of spacetime. On August 17th, 2017, LIGO and Virgo collaboration
observed the first evidence of gravitational wave. This not only proved the generalrelativistic, but also provided a new way and aspect to study the astrophysics andcosmology.
GWs are usually generated by the merger of two massive objects, such as neutronstar and neutron star, black hole and black hole, and neutron star and black hole.
After the merging, a magnificent amount of energy is released in the form of GW.
GW travels in the universe with the speed of light. Its direction does not change
and hardly lose any energy through its propagation, making it an extremely usefulmessenger.
Current observatories of GW including LIGO, Virgo and KAGRA. Each of thecan detect the GW signal from all directions. After a signal is received, an alert issent to all the astronomical observatories all over the world and following observation

13of other messengers can be carried out immediately.

1.5

PREVIOUS RELATED SEARCH• SKI-III:
In 2009, a time-integrated search on data from SKI to SKIII was conducted[58].
Only up-going muon sample was used. Sensitivity in this research is estimatedto 3 × 10−7 cm−2 s−1 in the region of declination of (−30◦ , 0◦ ). The lowest
neutrino flux limit averaged in right ascension is set to 6 × 10−8 cm−2 s−1 ataround dec= −54◦ .
• IceCube:
IceCube searched for the time and space clustering of high-energy neutrinos
time-dependently in 2020 using data from 2012 to 2017, and in 2015 using datafrom 2008 to 2012[5, 4].
• ANTARES:ANTARES searched for point-like neutrino sources time-integrally with 9-year
data and 11-year data in 2017 and 2020, respectively.[21, 43]

14
CHAPTER 2
Super-Kamiokande
2.1

OVERVIEW

Super-Kamiokande (SK) is a water Cherenkov detector. The detector lies 1,000 meters under Mt. Ikenoyama, Hida City, Gifu Prefecture, Japan. The main part of thedetector is a stainless steel cylindrical water tank with a diameter of 39 meters andheight of 42 meters. This water tank contains 50,000 tons of ultra pure water. On
the wall of the tank and outside of the tank, there are in total 13,000 photomultipliers (PMT). The rock above and around the detector create a good shelter whichshields most of the cosmic rays of energies below 1.3 TeV, and only about 1 × 10−5of the cosmic muons arrive at the earth surface can reach the detector. Figure 2.1 isa schematic of SK.

Figure 2.1: A schematic of SK.[44]The SK detector contains two parts: the inner detector (ID) and outer detector

15(OD). ID uses about 11,000 20-inch PMTs located 2.6 meters inward from the wallof the tank. OD is placed 2.2 meters from the wall and 2.1 meters from the top
and bottom of the tank, consisting of 1,885 8-inch PMTs. The OD wall is covered
with white Tyvek sheet to maximize light detection eﬀiciency. Additionally 26 setsof vertical and horizontal Helmholtz coils are lined up along the wall, reducing thegeomagnetic field inside the tank from 450 mG to 50 mG. The frame structure isshown in Figure 2.2.

Figure 2.2: Frame structure which supports PMTs.[35]SK started to operate since 1 April 1996, and has 6 completed phases (from SK-I

16
to SK-VI). During SK-II in 2001, more than half of the PMTs were destroyed due toan implosion of one bottom PMT, causing the ID photon coverage rate dropping to19%. In 2008, front-end electronics were upgraded. From SK-VI, gadolinium (Gd)was dissolved in the ultra-pure water to make it easier to separate neutrinos fromantineutrinos. The latest phase SK-VII started from 1 June 2022. In this phase theGd concentration will be increased from 0.01% to 0.03%.

2.2

CHERENKOV RADIATIONSK detect charged particles via the emitted Cherenkov lights that hit the PMTs. The
speed of light in dielectrics is smaller than that in the vacuum, making it possible
for ultra-relativistic charged particle to move with a velocity larger than that. Whensuch a situation occurs, an electromagnetic shock called the Cherenkov radiation iscreated. Using the Huygens-Fresnel principle, it is predicted that the emitted lights
form a cone centered at the particle traveling direction with an opening angle related
to the velocity of the particle:

Figure 2.3: Cherenkov radiation light cone.[9]17cos θc =

1
nβ

(2.1)

Here θc is the opening angle, n is the refractive index of the dielectric, and β = v/c.
The number of photons emitted per unit traveling length per unit wave length as afunction of wave length and β can be given by
d2 N
2πα
= 2
dxdλ
λ



1
1− 2 2
β n


(2.2)

where α is the fine structure constant. The β threshold of Cherenkov radiation thencan be obtained asβthr =

1n(2.3)

and using this the momentum and energy threshold is given bymc
n2 − 1mc2 nEthr = √n2 − 1pthr = √

(2.4)
(2.5)

In Table 2.1 thresholds for different particles in water are listed. PMTs in SK are
sensitive to light with wave lengths from 340 nm to 600 nm, and the spectrum in thisregion of water is shown in Figure 2.4.

2.3

PHOTOMULTIPLIER

Two types of PMT are installed in SK: 8-inch PMT for OD and 20-inch PMT for ID.
The ID PMT is Hamamatsu R3600 PMT developed by Hamamatsu Photonics
K.K., and the structure is shown in Figure 2.5. The photocathodes are made of18

Table 2.1: Cherenkov energy thresholds in water for different particles.
Particle Rest Mass [MeV] Threshold [MeV]
e±
µ±
π±
K±
p

0.511
105.7
139.6
493.7
938.3

0.775
160.3
211.7
748.8
1408.9Figure 2.4: Spectrum of Cherenkov radiation in water.[48]19

Figure 2.5: 20-inch PMT.[35]Bialkali (Sb-K-Cs). The quantum eﬀiciency peaks at around 360 nm to 400 nm to avalue of 22% as shown in Figure 2.6.
The OD PMT has been used in SK since the accident in 2001 is Hamamatsu
R5912 PMT. Each PMT is attached to a 60cm × 60cm × 1.3cm wavelength-shifting(WS) plate made of acrylic-doped panels to absorb ultraviolet and emit blue green
light which is closer to the sensitivity peak of the PMT.
Every PMT has been covered by a 12 mm ultraviolet transparent acrylic dome
and a fiberglass reinforced plastic (FRP) case to prevent chain explosion in SK-II(Figure 2.7).

2.4

ELECTRONICS AND DATA ACQUISITION SYSTEM

Since SK-IV, the front-end electronics used for data acquisition (DAQ) system ischarge-to-time converter (QTC) based electronics with ethernet (QBEE)[52].
Each QBEE board consists of 8 QTCs, and each QTC connects to 3 PMTs.
Charge signal from PMT is first integrated by QTC, and then the output signal is20

Figure 2.6: 20-inch PMT quantum eﬀiciency distribution as a functionof wavelength.[32]Figure 2.7: FRP case and acrylic cover for ID PMT.[35]21

Figure 2.8: Timing chart of QTC in QBEE.[52]sent to time-to-digital converter (TDC). The digital signal is then recorded by a fieldprogrammable gate array (FPGA) at 60 kHz for data-sort mapping (DSM) (Figure2.8). Four DSM chips are installed in each QBEE, and their outputs are sent to thesystem-interface controler (SIC) which sorts the data and stores them into first-in
first-out (FIFO) memory (Figure 2.9).

Figure 2.9: Block diagram of the QTC.[52]One QBEE has the dynamic charge range from 0.2 to 2,500 pC, and there are three
gains, 1, 1/7, 1/49, corresponding to small, medium, large, in each QTC channel.

22
Trigger used before SK-IV is the hardware trigger. Together with the upgradein 2008, the hardware trigger was replaced with software trigger, allowing the event
width to change according to the trigger type. The triggers include super low energy(SLE), low energy (LE), high energy (HE), special high energy (SHE), outer detector(OD) triggers, and after trigger (AFT). AFT is used for the neutron tagging system.

2.5

WATER AND AIR SYSTEM

The water used in SK is the underground water from the mine. The dust in thewater can scatter or absorb the light and cause the detector sensitivity to decrease.
The water may also be contaminated by the radioactive materials from the rocksuch as the radon and radium. These radioactive materials increase noise hits and is
the main background in low-energy neutrino measurements. The water purification
system is applied in SK for those reasons mentioned above. About 50 tons of water
are circulated and purified in this system per hour. Additionally a heat exchangesystem is installed to keep the water temperature at 13.6 degrees.
The air is also purified by an air purification system which pumps radon-free air
into the tank to prevent the radioactive materials being dissolved in the water.

23
CHAPTER 3
Upward-Going Muon
3.1

SK DATA SAMPLE

The events SK collected are classified into three types: fully contained (FC), partiallycontained (PC) and upward-going muon (UPMU). UPMU are additionally classified
as stopping, through-going and showering. Figure 3.1 shows the schematic of different
types of samples. For FC and PC, neutrinos interact with the water and create
secondary particles within the ID volume. The difference is that all charged particles
stop within the ID volume for FC, and at least one particle exits the ID volume
and deposit energy in OD for PC. For UPMU, the muon neutrinos interact with therock around the detector and create energetic muons that enters the detector fromdirections below horizon. If the muon stops and decays inside the detector, it isclassified as stopping muon. Both through-going and showering muons pass throughthe detector, but showering muons pass the shower selection.

Figure 3.1: Schematic of SK samples.[62]FC, PC and UPMU have mean energies of 1 GeV, 10 GeV, and 10 to 100 GeV
respectively. The simulated neutrino energy distribution is shown in Figure 3.2.

24

Figure 3.2: Expected atmospheric neutrino energy distribution fordifferent samples.[8]3.2DATA REDUCTION

SK receives a huge amount of events at the trigger level, yet only a tiny percent ofthem are the true signal events. The major parts of the background here is the CR
muons and PMT noises. In order to reject the background events and classify thesignal events, data reduction is of vital importance. Here I will briefly introduce themajor steps of UPMU reduction. More details for UPMU and FC, PC can be foundin [25].
In the first reduction, a charge cut which requires

6, 000 p.e. ≤ PE300 < 1, 750, 000 p.e.

(3.1)

is applied. Here PE300 is the total number of observed photoelectrons (p.e.) within 300
ns. This eliminate the low energy and extremely high energy events. The requirement
for the track length of UPMU is longer than 7 m, and a muon with a momentum of 1
GeV and a track length of 3.5 m produces 6,000 p.e.. If the ID charge is larger than
1,750,000 p.e., the muon fitters will fail because the electronics become saturated.
The second reduction downward-going muons are rejected, for they are mostly

25
CR muons. To do so, a logic tree of 7 different muon fitters specialized to fit stoppingmuons, through-going muons and muons with bremstrahlung are involved. For one
fitter, those events classified as upward with a goodness-of-fit larger than the threshold
are accepted, and those classified as downward with a goodness-of-fit larger than thethreshold are rejected. If the event is classified as horizontal with a goodness-offit above the threshold or all of the goodness-of-fits are below the thresholds, theevent will be passed to the next fitter and repeat the steps above. If none of thefitters classify the event as upward or downward with a goodness-of-fit larger than
the threshold, but at least one of them classified it as horizontal with a goodness-of-fit
above the threshold, it will be accepted. Otherwise, it is rejected.
The third reduction is eye-scanning. And eventually the final reduction classifies
the events. To classify whether an event is a stopping or through-going (including
showering) muon, a criteria is used to determine the entry or exit. If the number
of OD hits within 800 ns and with 8 m from the entry or exit projected by the
track is greater than 10, then this entry or exit is considered to be an entry or exit
signal. If both entry and exit signal exist, and the distance between the entry andexit is greater than 7 m, the event is classified as a through-going muon. If only
entry signal exists and the fitted momentum using Cherenkov light is not less than
1.6 GeV, corresponding to 7 m track length, it is classified as a stopping muon. For
the showering muons, an extra chi-squared test based on the observed charge and the
expected charge of non-showering muon is applied. Events pass all selections will bereconstructed.

26
3.3

EVENT RECONSTRUCTION

Information including energy, time, direction and neutrino type can be reconstructedusing the time, charge and position of PMTs. The current fitting algorithm usedby SK is called APFit. I will only cover the reconstruction procedure of UPMU
reconstruction, which is same to part of the single-ring FC and PC reconstruction.
[25] introduces the reconstruction for all situations.
First part is the vertex fitting. This is done using the time information of hit
PMTs.
The first step is the point fit. In this step, it is assumed that all photons from theCherenkov radiation are emitted at the same time from one point. The goodness ofthe fit is dined as:
1 X(ti − t0 )2
G(t0 ) =exp −N i2(1.5σ)2(3.2)

where ti is the time residual equals to the photon arrive time minus time of fly of thephoton for the ith PMT, t0 is a free parameter to be fitted to maximize G(t0 ) andσ = 2.5 ns is the PMT timing resolution.
The second step is to estimate the opening angle and direction. A goodness-of-fit
is used to test the combination of the opening angle and direction:
R θedgePE(θ)dθQ(θedge ) = 0
×sin (θedge )dPE(θ)dθ
θ=θedge!2(θedge − θexp )2
× exp −
2σθ2


(3.3)

Here PE(θ) is the angular distribution of the observed charge as a function of opening
angle θ (not the opening angle to be estimated, which is denoted as θedge ) at a certain
direction, θexp is the expected Cherenkov opening angle given the charge within θedge ,and σθ is the angular resolution.
In the third step for vertex fitting, the vertex position will be fitted again using

27
the TDC-fit. This time photons are not assumed to be emitted from one point, but
along the track length. Photon scattering is also taken into consideration. Finally,
MS-fit is used to make a preciser fit[46]. In this fit, the information of ring pattern
is used to search for the vertex position parallel to the direction. The upward-going
muons are generated outside the ID, so their vertices are constrained on the surfaceof ID.
The momentum can be estimated using the number of observed Cherenkov photons inside a cone of 70◦ centered at the ring direction. However, the momentum
estimated this way only includes the energy deposit in the detector. For UPMU, thisis only the lower bound of the muon momentum, for the energy dissipated before
reaching the ID, or even after leaving the ID for through-going and showering muons,is unknown. The estimated momentum is used as the energy for the UPMU events.

28
CHAPTER 4
Event Simulation
In this study, the background events dominate the data set where we search forthe potential astrophysical neutrino events. The evidence for an existing neutrino
point source is usually considered to be a cluster of multiple high-energy neutrinos
in space and time. The spatial and time distribution of atmospheric neutrinos aresupposed to be quite uniform in general, but could have coincident fluctuations which
produce a false appearance looks like a signal. We cannot estimate the probability
that such a situation occurs either by theoretical calculation or using the data we havecollected, but we can use a powerful tool called Monte Carlo (MC) simulation to doit. MC simulates with three steps: atmospheric neutrino flux, neutrino interactions,and detector simulations. We will not directly use the simulation introduced in thischapter for this study. Instead, we use the already existed MC to generate a subset
of toy MC events. Details of this will be talked about later in Section 5.4.1.

4.1

ATMOSPHERIC NEUTRINO FLUX

Atmospheric neutrinos are produced by CR. As mentioned in Section 1.4.1, the major
components of CR are protons and nuclei such as α-particles. These high-energy protons and α-particles strike the atmosphere and interact with the nuclei isotropically,
producing pions and kaons. Charged pions decay into muons and muon neutrinos:
π − → µ− + ν̄µπ + → µ+ + ν µ29Figure 4.1: Schematic the creation of atmospheric neutrinos.[47]and muons further decay into electrons, electron neutrinos and muon neutrinos:µ+ → e+ + νe + ν̄µµ− → e− + ν̄e + νµ
The neutrino flux model SK uses is the Honda-Kajita-Kasahara-Midorikawa (HKKM)
flux model[41]. HKKM model uses flux primary taken from experiments CR as input.
Solar activity and geomagnetic field, which affect the flux of primary CR, are also
taken into consideration. The theoretic models used for hadronic interaction between
CRs and air nuclei are JAM and DPMJET-III. The particle interactions and propagation are calculated in three dimensions to enhance the near-horizon neutrino flux.
A summary of the results from HKKM in comparison with Bartol flux[26] and Fluka

30
flux[27] are shown in Figure4.2.

Figure 4.2: Left:Calculated atmospheric neutrino fluxes averaged overall directions for different models. Right: Neutrino and anti-neutrino
flavor ratios for different models.[42]4.2

NEUTRINO INTERACTION

Interaction between neutrinos and hydrogen and oxygen atom in water molecule androck (primarily made of SiO2 ) are simulated by NEUT[39]. According to the exchanged weak boson (W ± or Z), the interactions are classified into charge current(CC) or neutral current (NC) two categories. The simulated interactions include:

CC/NC Quasi-Elastic and Elastic Scatteringνl + N → l/νl + N ′ /N

CC/NC Single Meson and Gamma Production νl + N → l/νl + N ′ /N + (π, K, η, γ)
CC/NC Deep Inelastic Scatteringνl + N → l/νl + N ′ /N + hadrons

CC/NC Coherent Pion Production

νl + 16 O → l/νl + 16 O + π31where N (N ′ ) is the nucleon (proton or neutron), l is lepton (not include neutrinos)and νl is the neutrino for l.
For NC elastic scattering, no new particle is created, and the neutrino just scattersoff a nucleon. The charged lepton that can be detected by SK is generated throughCC quasi-elastic (CCQE) scattering. Llewellyn-Smith model[45] is used for CCQE
cross-section of free protons. For nucleons bounded in an oxygen nucleus, Fermi
motion and Pauli blocking of the nucleon are taken into account in model by Smith
and Moniz[56].
A single meson can be produced via baryon resonance in the CC/NC neutrinonucleon interaction. Such a simulation is based on the Rein and Sehgal model[54]. 18
resonances together with the pion-less delta-decay have been taken into consideration.
When neutrino energy exceed 10 GeV, the interactions between neutrino andquark in the nucleon which generate multiple hadrons via deep inelastic scattering(DIS) dominates the interactions. This is the dominant interaction for the upward2

d σgoing muons. The cross section is by integrating dxdy
, which is given by[39]:d2 σ
G2F MN Eν
1=(1 − y + y 2 + C1 )F2 (x, q 2 )
dxdy
π
2

1
2 2
±y(1 − y + C2 )[xF3 (x , q )]
2

C1 =

xyMN
M2
yMl2
Ml2
−
− l2 −
4MN Eν x
2Eν
4Eν
2MN Eν x
Ml2
C2 = −
4MN Eν x

(4.1)

(4.2)
(4.3)

32
Here, x = −q 2 /(2M (Eν − El )), y = (Eν − El )/Eν , Mn and Ml are the mass of thenucleon and lepton, respectively, Eν and El are energies of the incoming neutrino andoutgoing lepton in laboratory frame. F2 and xF3 are the nucleon structure functions
taken from GRV98[36]. Figure 4.3 shows the calculated total charged current crosssection in comparison with the experiment data.

Figure 4.3: Total charged current cross sections for neutrinos (top 4lines) and antineutrinos (bottom 4 lines).[39]Beside the baryon resonance, a pion can also be produced by the interaction
between neutrino and oxygen nucleus called coherent pion production. Such an interaction is also simulated in Rein and Sehgal model.

4.3

DETECTOR SIMULATION

The detector simulations include particle propagation, Cherenkov photon generationand propagation, and PMT response. These are simulated by a GEANT3 based
package, SKDETSIM.

33
For hadronic interactions except for low energy pions, GCALOR[64] is used. Low
energy pions are simulated using a custom program[50]. Charged particle propagation
are simulated based on what we discussed in Section 2.2.
In Cherenkov photon propagation simulations, Rayleigh scattering, Mie scatteringand absorption by water molecule are considered. Reflection on the wall of the tankis also calculated.
PMT charge and time response simulation use the measurements from calibration.

34
CHAPTER 5
Analysis
5.1

DATA SAMPLE

In this study, we look for the high-energy astrophysical neutrino point sources in SKIV using only UPMU events. We used 4547 UPMU events from SK-IV with a total
live time of 3288.8 days for the analysis, and also the MC equivalent to 500.0 years
simulation are used for the calculation of the PDF in Section 5.2.2 and the generationof the toy MC in Section 5.4.1. Figure 5.1 shows a sky map of all UPMU events andFigure 5.2 is the time distribution of all UPMU events.

Figure 5.1: A sky map of different types of UPMU events. All the sky
map in this and next chapter are shown in J2000 equatorial coordinate.
Since we are searching the full sky without relying on any other triggers, it is very
important to have good angular resolutions for all the events we use. As shown in
5.3, UPMU events have the best angular resolution among the all three samples.

35

Figure 5.2: Time distribution of different types of events. Time ispresented in the modified Julian day (mjd).

36

Figure 5.3: Angular resolution corresponding to 1σ for FC, PC and
UPMU from atmospheric neutrino simulation.
5.2

UNBINNED MAXIMUM LIKELIHOOD METHOD

5.2.1

Likelihood

In this study, events are classified into two categories: either a signal event or abackground event. For the data set we use, there are two hypothesis: a null hypothesis
which means no signal events is in the data set, and another hypothesis that at least
one signal event is in the data set. The direction, energy and time PDF for the signal
and background events are different.
For the binned method, events within a certain angular region and time intervalaround the assumed source direction and burst time are selected. This region andinterval can be chosen based on the detector angular resolution and burst model. By
comparing the number of selected events and the expected background events, we
can get the probability of having such a number of events without any signal events.
Then we can evaluate the expected number of signal events if the probability is lower

37
than a confidential level threshold.
For the unbinned method that is used in this study, each event does not have tobe taken as either a signal or a background. Instead, it is a combination of a signal
component and a background component. The probability density for ith event in aset of N events is:
nS
nB
Si +
Bi
nS + nB
nS + nB(5.1)

Then the likelihood for N events is:
enS +nB (nS + nB )N Y
L=
N!
i=1
N



nS
nB
Si +
Bi
nS + nB
nS + nB(5.2)

Here the nS (nB ) is the expected number of signal (background) events, Si (Bi ) is
the signal (background) probability density function (PDF). nS is the only unknown
variable here, so we can fit nS to maximize L(nS ) and obtain the expected numberof signal events.
Both binned and unbinned methods are top-down searches and are model independent. This is a much more eﬀicient choice compared with the bottom-up search,as there are too many possible models to test, but little progress has been obtained.
In this study we will conduct both a time-integrated search and a time-dependent
search for high-energy astrophysical neutrino point sources. These two methods haveadvantages in different neutrino emission duration. Time-integrated search have abetter performance for sources with emission duration longer than 100 days, while
time-dependent search will have a better discovery potential for shorter time duration[4].
When calculating the likelihood, if all events are included, it will take an extremely long time. This is neither eﬀicient nor necessary. The angular resolution for

38
UMPU have been discussed in Section 5.3 and it can be seen in Figure5.3 that eventsreconstructed to energy above 1 GeV have angular resolution smaller than 15◦ . As aresult, it is reasonable for us to calculate likelihood with only events within a searchcone centered at the source direction with an opening angle of 15◦ .
5.2.2

Variables

Before introducing PDF, it is important to take a look at the variables that will be
used in the searching algorithm. In time-integrated search, only angular and energy
PDF will be included in the likelihood Equation 5.3 and we calculate L(nS , γ) as afunction of nS and energy spectrum index γ. For time-dependent search, Equation
5.2 will also consist of a time PDF, and L(nS , γ, t0 , σt ) is calculated as a function of
two extra parameter, emission center time t0 and emission duration σt . Details on
this will be discussed in the following section.
The Si and Bi in Equation 5.3 is given as
Si = ASi (∆ϕi |Ei )EiS (Ei |γ)

(5.3)

B
B i = ABi (deci |dec0 )Ei (Ei )(5.4)

for time-integrated search and

Si = ASi (∆ϕi |Ei )EiS (Ei |γ)TiS (ti |t0 , σt )

(5.5)

B
B
B i = ABi (alti |dec0 )Ei (Ei )Ti

(5.6)

S(B)

for time-dependent search. Here Ai

S(B)

, Ei

and TiPDF and time PDF for signal (background).
The definition of each variable is listed below:S(B)are the angular PDF, energy

39• ∆ϕi : the angular separation between reconstructed direction of ith muon event
and assumed point source in equatorial coordinate.
• dec0 : declination of assumed point source in equatorial coordinate.
• Ei : reconstructed energy of ith muon event.
• γ: spectrum index.
• deci : declination of ith muon event in equatorial coordinate.
• alti : altitude of ith muon event in detector coordinate.
• ti : detection time of ith muon event.
• t0 : center of neutrino emission time.
• σt : duration of neutrino emission.
5.2.3

Probability Density Function

It is important to notice that the PDF time-integrated search and time-dependent
search do not only differ in the time part, but also have a difference in the angularpart for background. The PDF are calculated using the MC.
The angular PDF of signal events depends on their reconstructed energy. Theenergy interval used for this study is divided into bins with same number of eventsin each bin. The distribution of the scattering angles between the reconstructed
direction of muons and the true direction of the neutrinos are fitted with a Gaussianlike function in each bin. Then the PDF is normalized between 0◦ and 15◦ . Azimuth(φ) distribution is taken to be uniform. Plots of AS (∆ϕ|Eµ ) = P (∆ϕ|Eµ )P (φ) isshown is in Figure 5.4.

40

(a)

(b)

Figure 5.4: AS (∆ϕ|Eµ ) for signal events at different reconstructedenergy. Left: P (∆ϕ|Eµ ). Right: P (φ)
It is more complicated for the background angular PDF. Using the MC, we canobtain the angular PDF in the detector coordinate. But when we are searching for
astrophysical point sources, we look for them in the equatorial coordinate. There isan one-to-one correspondence between the altitude and azimuth, and declination andright ascension for a given local sidereal time. Thus we can calculate the declination
PDF in the following way:Z24
P (dec, ra) =P (dec, ra, t)dt
0
Z24

=

P (alt, azi, t)dt

(5.7)

0
Z24

=

P (alt)P (azi)P (t)dt0Here P (alt) and P (azi) are the altitude and azimuth PDF in the detector coordinate(Figure 5.5) and P (t) is the local sidereal time (lst) distribution (Figure 5.6). The
integration of P (t) and Jacobian matrix equal to 1, hence are omitted here.

41

(a)

(b)

Figure 5.5: AB (dec, ra) for background events used in time-dependent
search. Left: P (alt). Right: P (azi).

Figure 5.6: Local sidereal time (LST) distribution.

42
Taking the right ascension distribution to be uniform, we have AB (dec) = P (dec)P (ra)(Figure 5.7). The probability of being inside the search cone is:

(a)

(b)

Figure 5.7: AB (dec, ra) for background events used in time-integrated
search. Left: P (dec). Right: P (ra).

◦Z

P ((dec, ra) ∈ {∆ϕ(dec, ra) ≤ 15 }) =P (dec, ra) cos (dec)d(dec)d(ra)∆ϕ(dec,ra)≤15◦(5.8)
And eventually we have the angular PDF with the precondition of being in the searchcone for time-integrated search:
AB (dec|dec0 ) = P (dec, ra|∆ϕ(dec, ra) ≤ 15◦ })

(5.9)
◦

= P (dec)P (ra)/P ((dec, ra) ∈ {∆ϕ(dec, ra) ≤ 15 })Things are, on the other hand, a bit different for the time-dependent search, sincevariables in P (dec, ra, t) are correlated. We have to calculate the PDF in the following

43
way:
P (dec, ra, t|∆ϕ(dec, ra) ≤ 15◦ }) = P (dec, ra|(dec, ra) ∈ {∆ϕ(dec, ra) ≤ 15◦ }, t)P (t)P (dec, ra, (dec, ra) ∈ {∆ϕ(dec, ra) ≤ 15◦ }|t)P (t)P ((dec, ra) ∈ {∆ϕ(dec, ra) ≤ 15◦ })P (dec, ra|t)P (t)=P ((dec, ra) ∈ ∆ϕ(dec, ra) ≤ 15◦ })P (dec, ra, t)=
P ((dec, ra) ∈ {∆ϕ(dec, ra) ≤ 15◦ })P (alt, azi, t)=
P ((dec, ra) ∈ {∆ϕ(dec, ra) ≤ 15◦ })P (alt)P (azi)P (t)=P ((dec, ra) ∈ {∆ϕ(dec, ra) ≤ 15◦ })(5.10)=

The angular PDF with the precondition of being in the search cone for time-dependent
search is given by:

AB (alt|dec0 ) = P (alt)P (azi)/P ((dec, ra) ∈ {∆ϕ(dec, ra) ≤ 15◦ })(5.11)The energy distribution of high-energy astrophysical neutrinos has not been confirmed experimentally, so we assume they follow a spectrum distribution, but leave
the spectrum index free to change in a range. With this, we can calculate the energy
PDF of signal event:Z Z
E(Eµ , γ) =Eν−γ P (Eν , Eµmin )S(z, Eν )f (Eµ |Eν )dEν dz/N

(5.12)

N is the normalization factor.
The first part we need is the probability that a neutrino of energy Eν is scatteredoff into a muon with an energy larger than the minimum energy Eµmin to be detected
as an UPMU event P (Eν , Eµmin ) (Figure 5.8). We also need to take into account the44earth shadow effect S(z, Eν ) (Figure 5.9), where z ≡ cos(θ) and θ is the zenith angle
in the detector coordinate.

Figure 5.8: Probability P (Eν , Eµmin ) for a neutrino of energy Eν creating a muon of energy larger than Eµmin .
Then we calculate the energy response function f (Eµ |Eν ) (Figure: 5.10), whichis the probability of a neutrino of Eν being reconstructed to an UPMU event of Eµ .
Just like what we did when calculating the angular PDF for signal events, we divide
the energy interval of MC into small bins, yet the reconstructed energy is changed toneutrino energy. In each bin, we calculate the PDF of reconstructed energy.
Background energy PDF can simply be obtained from the MC. Figure 5.11 shows
a comparison of the E S (Eµ , γ) and E B (Eµ ).
Shape of signal time PDF is also unknown to us, so here we use the Gaussian
distribution
1 ti −t0 2
1
T (ti |t0 , σt ) = √ e− 2 ( σ )σ 2π(5.13)and leave t0 and σt as free parameters to be fitted. For background time PDF, we

45

(a)

(b)

Figure 5.9: Earth shadow effect S(z, Eν ) for νµ (left) and ν̄µ (right)
as a function of log10 Eν and z = cos θ.

Figure 5.10: Top: The probability that a neutrino of energy Eν creates a muon of reconstructed energy Eµ (f (Eµ |Eν )).Bottom: The probability that a muon of reconstructed energy Eµ is created by a neutrinoof energy Eν46

Figure 5.11: Expected reconstructed energy distribution for background events and signal events with different spectrum indices.
1
assume it to be uniformly distributed Tend −T, where Tend − Tstart is the operation
starttime of SK-IV (not live time). (Figure 5.12)
5.2.4

Expected Number of Background Events

SK receives atmospheric neutrino events at a quite stable event rate in each direction.
This again is in the detector coordinate. In order to calculate the expected number ofbackground events in the search cone, we need to transform to equatorial coordinate.
By assigning local sidereal time to each event in MC according to the distribution, we
can calculate their corresponding declination and right ascension. Then we count the
number of events within the search cone at all searching directions. Finally, all nB
are weighted by TSK−IV, data and MC live time ratio. In Figure 5.13, a comparison
TM C
between the calculated expected nB and the true event distribution of the SK-IV
data is shown, together with the nB as a function of declination averaged over right

47

(a)

(b)

Figure 5.12: T S (t|t0 , σt ) for signal (left) and T B for background(right)ascensions.
5.2.5

Test Statistics

We have defined the likelihood in a general sense in Section 5.2.1. In order to calculate
the probability that such a situation is generated by the fluctuation of background
events, we need to make a comparison with the null hypothesis (nS = 0). To do so,
we can calculate the ratio of two likelihood and define it as the test statistics (T S).
The T S for time-integrated search is defined as:L(n̂S , γ̂)T S = 2 log
L(nS = 0)NXn S Silog 1 +
− 2nS
=2
n
B
B
i
i=1

(5.14)

where n̂S and γ are the best fitted number of signal events and the energy spectrum
index. The probability of observing a given T S with null hypothesis will approximately accord to the χ2 distribution.
In the time-dependent search, we must pay attention to a not so obvious fact which

48

(a)

(b)

(c)

Figure 5.13: Top left: Sky map of expected number of background
events within the search cone. Top right: Sky map of number of events
within the search cone for the data. Bottom: Expected number ofbackground events as a function of declination.

49
is easy to be neglected. For background events with a uniform time distribution, theprobability for nB events to be within 2σt is smaller for shorter σt . On the other
hand, if σt becomes shorter, T (ti |t0 , σt ) will have a larger peak value. As a result,the algorithm will prefer a shorter σt , which is reasonable. Yet we must also realizeanother aspect that for the same searching period, the number of independent waysof choosing t0 , or in another way of saying, the number of non-overlapping region
with a time width of 2σt is larger for shorter σt , meaning that the probability thatnB background events appear within at least one of the time region becomes larger.
In order to include this effective trail factor in the algorithm, we can integrate thelikelihood over all possible t0 for a given σt :1
L(nS , γ, σt ) =Tend − Tstart

T
Zend

L(nS , γ, σt , t0 )dt0

(5.15)

Tstart

The integral can be approximated by a penalty term times the maximized original
likelihood[29]:√L(n̂S , γ̂, σ̂t ) ∼2πσ̂tL(n̂S , γ̂, σ̂t , t̂0 )Tend − Tstart(5.16)

Where Tend and Tstart are the ending and starting time of SK-IV. Ignoring the constant√
2π, the modified T S then becomes:L(n̂S , γ̂, σ̂t )σ̂t×T S = 2 log
Tend − Tstart
L(nS = 0)




N
X
n S Siσt=2
log 1 +
− 2nS + 2 log
nB B i
Tend − Tstarti=1
(5.17)

50
5.3

ALL-SKY SEARCH

So far we have discussed the construction of T S and its elements in details. Our next
step is to obtain the value of the free parameters that maximize the likelihood andthe corresponding T S at each searching direction. In general, we will carry out the
all-sky search in the following steps:
1. Obtain the best fitted n̂S , γ̂, t̂0 and σ̂t which maximize L and the corresponding
T S at each direction.
2. Find the direction with the overall maximum T Smax , together with its corresponding n̂S , γ̂, t̂0 and σ̂t .
We will discuss the detail procedures for each step in the following sections.
5.3.1

Parameters Fitting and Likelihood Maximization

We will maximize the likelihood numerically because the energy PDF of signal andits dependence on γ is obtained numerically, and there are to many terms containedin the likelihood. But the PDF of signal part can be extremely large and exceed the
maximum limit of number storage for a signal like event, so instead we will minimize
the minus logarithm likelihood, which has opposite monotonicity and same extreme
point as the original likelihood.
As has been mentioned in previous sections, the events we are looking at arebackground dominated, which means in most of the cases, nS will be fitted to be 0.
The number of events within the search cone at a given direction N will fluctuate so
it is possible for that number to be smaller than the expected number of background
events nB , and the likelihood will maximize at negative nS to make up the ”missing”

51
events in the search cone. We are not interested in such a situation, so the lower limit
of nS is set to be zero (see Figure 5.14).

Figure 5.14: T S distribution for time-integrated search in a pseudo
experiment with solely background toy MC.

5.3.2

First Guess of Time Parameters

The minimizer we use here is iminuit. This minimizer is powerful to find the local
maximum (it is supposed to be the minimum not maximum, but since the minimum of minus logarithm likelihood corresponds to the maximum of likelihood, I will
use maximum instead for convenience) if given an accurate initial value. For thetime-integrated search, it is expected that there is only one local maximum, so thismaximum is the global maximum. On the other hand, a large amount of local minima may exist in the time-dependent search, especially when σt is small, as shownin Figure 5.15. If we set the initial value of t0 to be Tstart or Tend , this minimizer
will return the first or last local maximum in the search time region, which does not52necessarily to be the global maximum. As a result, we must make an accurate enoughfirst guess of t0 and σt , and use them as the initial value for the fit.

Figure 5.15: T S as a function of t0 with fixed σt in a pseudo experiment using toy MC. Vertical dashed blue line represents the assumed
emission center time. Simulated with nS = 10 on top of the backgroundtoy MC, σt = 105 s, at dec = 30◦ .
A method to obtain the first guess values is intro[29] and is adjusted for this study.
To do so, only a part of the events in the search cone are needed because those faraway from the source does not have much effect on fitting the time parameters. Afterbalancing between the accuracy and calculation speed, we decide to use events within
7.5◦ from the searching direction. Then we arrange these events in time order, anddivide them into groups. The ith group contains n consecutive events from the ith
to the (i + n − 1)th. Here n is an integer in range [2, 8], as the results hardly change
for n > 8. For each group of each n, we calculate the likelihood with n as nS , t0 andσt as the mean value and standard deviation of ti respectively and fixing γ to 2. And
finally we use the t0 and σt corresponding to the overall largest likelihood as the first

53
guess value.
5.3.3

Search Directions

In order to do the all-sky search for the point sources, we need to divide the fullsky into small bins. A usual way to do this is to use the equally sliced grids inequatorial coordinate, for example 1◦ × 1◦ grids. This, however, has a disadvantagethat the corresponding solid angle of each grid will decrease as the declination getslarger, making the angular uncertainties of the search nonuniform. An improved waythat has been widely used by astronomy research is HEALPix method[3], which will
divide the full sky into pixels with the same solid angle. Fig5.16 shows an exampleof how it works. The searching directions are set to be the positions of pixel centers.
In principle, the trial factor will become larger if we more searching directions, as

Figure 5.16: A schematic of how HEALPix divide the full sky into
pixels.[3]the probability for at least one T S in all directions is larger than a certain value
increases. The number of pixels that full sky is divided into is eventually decided to

54
be 49152 with an separation between two pixels being 0.5◦ . This is smaller than thebest angular resolution for the highest-energy events of SK, and reduce the possibilityof the clustering of point sources in one pixel as well.
After calculating T S for all searching directions, we get a sky map of T S, as shownin Figure 5.17.

Figure 5.17: A sky map of T S for time-dependent search in pseudo
experiment with solely background toy MC.

5.4

PSEUDO EXPERIMENTS

After performing the full sky time-integrated and time-dependent search, we needto estimate the significance of the result with respect to the background. To do so,
we calculate two different p-values using the results from simulation. The pre-trial
p-value is calculated by comparing T Smax with the distribution of T S from the same
position obtained by simulating using toy MC 106 from trials. In order to take into
account this trial factor, we also need to calculate the post-trial p-value. To do so, we
first run many pseudo experiments called the trials and get a distribution of T Smax

55
in each trial. Then we compare the T Smax we obtain from the data with this andcalculate the post-trial p-value. In this section, I will introduce how to run the trials.
5.4.1

Toy Monte Carlo

A set of simulated events is needed to be thrown for each trial. We can generate theso called toy MC, which only contains the required information for the PDF. Eachvariable in the toy MC will follow the distribution calculated from the MC in Chapter
4.
We generate one background toy MC in the following procedure:
1. Assign energy according to the distribution shown in Figure 5.11.
2. Assign altitude and azimuth according to the distribution in Figure 5.7 or 5.5.
3. Assign Julian date with the following steps:(a) Choose which run the event is in. Each run is weighted by its duration
times its live time rate.
(b) Assign a time in the duration of that run assuming a uniform distribution.
4. Calculate the equatorial coordinates (right ascension and declination).
The number of background toy MC events follows Poisson distribution, with the meanvalue being the event rate times total live time for SK-IV.
Besides the threshold of T S, we also what to evaluate the sensitivity of this searching algorithm for SK. This not only requires events sets containing solely backgrounds,
but also several signals on top of backgrounds. Before generating signal toy MC, weneed to determine the direction and time information for the assumed point source.
Procedure for the generation of one signal toy MC for time-integrated search and56time-dependent search are mostly the same, with only a difference in Julian dataassignment:
1. Assign Julian date:• For time-integrated search: the same as the Julian date assignment forbackground toy MC.
• For time-dependent search: assign a random value according to the Gaussian distribution of t0 and σt . If this lies in the dead time region (the time
between two runs or the time for test runs), we repeat this again.
2. Assign true neutrino energy according to Eν−γ

R

P (Eν , Eµmin )S(z, Eν )dz.

3. Assign reconstructed muon energy using f (Eµ |Eν ).
4. Assign scatter angle ∆ϕ based on A(∆ϕ|Eµ ).
5. Calculate the equatorial coordinates.
6. Calculate altitude and azimuth in the detector coordinate.
7. If the event is above horizon, we repeat the steps again.
5.4.2

Search Algorithm Confirmation

By running the trials with toy MC, we can check if the algorithm is working properly to find the position and time information of the assumed source. Figure 5.18
and Figure 5.20 are sky maps of T S for time-dependent and time-integrated search,
respectively. It is can be seen that the T S distribution for trial with background
only toy MC is quite uniform, while those with signals have a T Smax located close to
the neutrino source position we have set at different declination. In Figure 5.19, we
calculated the T S as a function of t0 and σt for different emission duration.

57

(a)

(b)

(c)

(d)

(e)

(f)

Figure 5.18: Sky map of T S of time-dependent search for solely
background 5.18a and nS = 7 at dec = 30◦ 5.18b, 0◦ 5.18c, −30◦ 5.18d,−60◦ 5.18e, −90◦ 5.18f. Red cross represents the location of assumed
source, and yellow plus represents the direction of T Smax .

58

(a)

(b)

(c)

(d)

(e)

(f)

(g)

Figure 5.19: T S as a function of t0 and σt for nS = 2 σt = 0s 5.19a,
nS = 3 σt = 10s 5.19b, nS = 5 σt = 100s 5.19c, nS = 7 σt = 1, 000s
5.19d, nS = 10 σt = 104 s 5.19e, nS = 13 σt = 105 s 5.19f, nS = 17σt = 106 s 5.19g. Vertical and horizontal dashed green line represents
the assumed log10 σt and t0 , respectively. Yellow cross represents thebest fitted log10 σt and t0 .

59

(a)

(b)

(c)

(d)

(e)

(f)

Figure 5.20: Sky map of T S of time-integrated search for solely background 5.20a and nS = 20 dec = 30◦ 5.20b, nS = 20 dec = 0◦ 5.20c,
nS = 25 dec = −30◦ 5.20d, nS = 30 dec = −60◦ 5.20e, nS = 40
dec = −90◦ 5.20f. Red cross represents the location of assumed source,
and yellow plus represents the direction of T Smax .

60
5.4.3

Sensitivity

By running many trials for a given nS , we can find the T Smax in each trial and geta distribution of T Smax . Figure5.21 is the T Smax distribution for pure backgrounds.
The detection threshold is defined the T Smax value of 90% CL due to lack of previous
detection. We then define the sensitivity as the nS needed so that 50% of the T Smax
is larger than the detection threshold. If none of the nS matches the threshold, weinterpolate linearly between two successive values of nS , as shown in Figure5.22.

Figure 5.21: T Smax distribution for solely background events. Vertical red solid line indicates the 90% CL threshold.
Sensitivity for time-integrated search is calculated as a function of the source
declination and converted to neutrino fluence using Equation 5.24, as shown in Figure
6.6. And sensitivity for time-dependent search is calculated at different declinationand σt , shown in Figure 5.24. Using Equation 5.24, we can convert the sensitivity
obtained from Section 5.4.3 into fluence (shown in the right panel of Figure 6.6 and5.24).

61

Figure 5.22: T Smax distribution with different nS . Vertical red solidline indicates the background 90% CL threshold. Vertical blue andgreen line represent the 50% CL for nS = 2 and nS = 3.
The calculated sensitivity strongly relies on the number of trials we run, especially for the cases where the distribution for two successive nS values are similar.
Such a situation occurs when calculating sensitivity for the time-integrated search at
declination above 85◦ . Running the trials is supper time-consuming, for we only get
one T Smax value in one trial, but need to perform nearly 50,000 times search in one
trial. As a result, we can only provide a less accurate result here.

5.5

NEUTRINO FLUENCE

The final results we are interested in is not nS but neutrino fluence Φ(GeVcm−2 ).
The relationship between the neutrino fluence and muon fluence is given by:
EZ
ν,max

dEν

Φµ =
Eν,maxX
ν=νµ ,ν̄µP ν (Eν , Eµmin )S ν (z, Eν )

dΦν
dEν

(5.18)

62

(a)

(b)

Figure 5.23: Sensitivity for time-integrated search in terms of nS(top) and Φ (bottom).

63

(a)

(b)

Figure 5.24: Sensitivity for time-dependent search in terms of nS(top) and Φ (bottom).

64
By assuming spectrum distribution, we have:dΦν= Φ0 Eν−γdEν

(5.19)and muon fluence can be expressed by:

Φµ =

nS
Aef f

(5.20)

Substitute Equation 5.19 and Equation 5.20 into Equation 5.18, we have:Φ0 (Eµ > Eµmin ) =
Here

EZ
ν,max

dEν

ψ(z) =
Eν,maxXnS
Aef f (z)ψ(z)P ν (Eν , Eµmin )S ν (z, Eν )Eν−γ(5.21)(5.22)

ν=νµ ,ν̄µis the detection eﬀiciency and Aef f (z) is the effective area for SK, γ is assumed to be
2 and, Eν,min = 1.6 GeV, Eν,max = 105 GeV[58].
We have introduced P ν (Eν , Eµmin ) and S ν (z, Eν ) in 5.2.2. With these two parts,dnν
we average them with respect to Eν together with dE(Eν ) ∝ Eν−γ and get ψ as aνfunction of z (shown in Figure 5.25).
The effective area of SK can be calculated using a numerical technique from [55]:• Take a large plane next to a cylinder with the same dimensions of SK.
• Construct a grid of 10cm × 10cm squares on the plane.
• Draw a ray perpendicular to the plane at each grid vertex through the cylinder.
• If the path length inside the cylinder is larger than 7m, the grid is counted as

65

(a)

(b)

Figure 5.25: ψ(z) for νµ (left) and ν̄µ (right)
a part of the effective area contributing 100cm2 .
Figure 5.26 is a schematic of how effective area is calculated. This process is repeatedFigure 5.26: Schematic of calculating effective area.[55]for every zenith angle. Then we get the effective area as a function of z, as shown in
Figure 5.27
There is one last step we need to do, that is convert these from functions of zenith
angle to functions of declination using the relation between zenith angle, declination

66

Figure 5.27: Effective area as a function of z.
and local sidereal time:Rψ(dec)Aef f (dec) =dtψ(z)Aef f (z)P (t){t|z(t|dec)<0}R

dtP (t)(5.23){t|z(t|dec)<0}

Figure 5.28 shows the result of Equation 5.20 This is because the zenith angle isalways changing for a certain declination as the Earth is rotating. Then we have:Φ0 (Eµ > Eµmin ) =nS
Aef f (dec)ψ(dec)(5.24)(5.25)

67

Figure 5.28: Product of ψ(dec) and Aef f (dec) as a function of declination.
We can also calculate the flux by:
EZ
ν,max

φ = Φ0dEν Eν−γ /TEν,max−γ+1−γ= Φ0 (Eν,min− Eν,max)/T

For time-integrated search, T = 2.84 × 108 s, the total live time of SK-IV.

(5.26)

68
CHAPTER 6
Results and Discussion
6.1

RESULTS

6.1.1

Time-Integrated Search

T S sky map and histogram of time-integrated search are shown in Figure 6.1 andFigure 6.2. T Smax = 13.32 occurs at (198.3◦ , 15.1◦ ). This value is much smaller than
the 90%C.L. T S 90%C.L. = 18.76, and the post-trial p-value is 79.1% (see Figure.6.3).
The pre-trial p-value is 1.18 × 10−4 (Figure 6.4).

Figure 6.1: Sky map of T S for time-integrated search. Location of
T Smax is marked with red triangle. Positions with n̂S = 0 has been
excluded from the sky map, which is represented by the white color
region.
Details of the fitted parameters are listed in Table6.1. Figure 6.5 shows a zoomedin map of the events around T Smax .
We calculated L(nS |γ̂) as a function of nS with γ fixed at best fitted value γ̂, and69

Figure 6.2: Blue filled histogram: T S distribution for time-integrated
search from all searching directions from data. Red dashed line histogram: T S distribution from all searching directions obtained from50 trials (roughly 2 million times searches) with solely background toy
MC. The peak around 0 represents the situation where n̂S = 0. The result from the data basically agrees with the simulation for background.

Table 6.1: Table of the best fitted parameters for time-integrated
search.
nSγra [◦ ]

dec [◦ ]pre-trial p

post-trial p13.32 7.913

2.58

198.3

15.1

1.18 × 10−4

79.1%

TS

70

Figure 6.3: The histogram of post-trial T Smax distribution for timeintegrated search obtained from 1,000 trials with solely background toy
MC.

Figure 6.4: Histogram of pre-trial T S distribution for time-integrated
search by simulating 106 times at the position of T Smax .

71

Figure 6.5: Stopping (red plus), through-going (green dot) and showering (blue cross) UPMU events detected around the position of T Smax .
Grey region encircled by the black circle is the 15◦ search cone.
calculated upper limit by
R nupperS
90% = R0 +∞−∞dnS L(nS , ...)

dnS L(nS , ...)

(6.1)

The result nupper= 13.86, with corresponding fluence upper limit Φupper= 2.06 ×
0S
102 GeVcm−2 and flux upper limit φupper = 4.53 × 10−7 cm−2 s−1 :6.1.2

Time-Dependent Search

Figure6.7 and Figure 6.8 show the sky map and histogram of T S. The maximum T S
occurs at (118.0◦ , −62.7◦ ), with T Smax = 14.53. This value is smaller than the thethreshold T S 90%C.L. = 17.8. The post-trial p-value is 33.4%, with the pre-trial p-value
found to be 2.24 × 10−4 , as shown in Figure 6.9 and Figure6.10. The best fitted timeof the emission is 13:09:42 JST, October 27th , 2010, with a emission duration (2σ̂t )72Figure 6.6: L(nS |γ̂) as a function of nS with γ fixed to the best fitted
value γ̂. The vertical solid blue line represent the upper limit for 90%
CL.
of around 5.1 hours. Rest of the best fitted parameters of the top three T S are listed
in Table 6.2.
Table 6.2: Table of the best fitted parameters of the top three largest
T S for time-integrated search.
TS

nS

γ

t0 [mjd]σt [s]

14.53 1.996 1.56 2010.10.27 04:04:48 9.25 × 103
13.96 1.999 2.50 2009.09.27 03:07:12 3.37 × 103
12.07 3.829 3.00 2016.03.03 21:24:30 1.36 × 105

ra [◦ ]

dec [◦ ]pre ppost p

118.0
3.4
277.9

-62.7
12.0
-60.2

2.24 × 10−4
5.08 × 10−5
6.09 × 10−4

33.4%
39.5%
69.2%

Since no excess was found with respect to the background, we calculate the upper
limit of the neutrino flux. In Figure 6.11, we calculated the L(nS |γ̂, t̂0 , σ̂t ) as afunction of nS at the position where T Smax occurs. The results are nupper= 5.29,S
with corresponding Φupper= 78.6 GeVcm−2 for the 1st largest T S, nupper= 5.27, with
0S
= 78.1 GeVcm−2 for the 2nd largest T S and nupper= 7.84, withcorresponding Φupper
0S
corresponding Φupper= 117.3 GeVcm−2 for the 3rd largest T S. The emission time0

73

Figure 6.7: Sky map of T S for time-dependent search. The locationsof the top three maximum T Ss are marked with red triangle, plus andcross, respectively. Positions with n̂S = 0 has been excluded from thesky map, which is represented by the white color region.
interval is estimated by 2σ̂t (1.85 × 104 s, 6.74 × 104 s and 2.72 × 105 s for 1st , 2nd
and 3rd , respectively).
It is clear that most of the events are more background-like, but we do findsomething interesting at the largest three T S. At least two events with offset lessthan 5◦ from the fitted source direction in a time region of ±σ̂t are found at eachof these three positions. The information of are listed in Table 6.3. Figure 6.12
shows the orientation of the events with their angular uncertainties and their time
information.

6.2

POTENTIAL CANDIDATE CHECK

After obtaining all the results, we tried to find potential sources which is coincident
with the largest three T S. The first candidate is the TXS 0506+56 discussed in
Chapter 1, but none of these three matches either with direction (77.36◦ , 5.69◦ ) or

74

Figure 6.8: Blue filled histogram: T S distribution for time-integrated
search from all searching directions from data. Red dashed line histogram: T S distribution from all searching directions obtained from 50
trials with solely background toy MC. The peak around 0 represents
the situations where n̂S = 0 and σ̂t = (Tend − Tstart )/2, the maximum
value set for σt . The result from the data basically agrees with thesimulation for background.
Table 6.3: Energy, distance from best fitted position, time intervalfrom the best fitted t0 , and event type of the events within 15◦ to theposition of the top three largest T S with in 7σ̂t around t̂0 .
Ei [GeV] ∆ϕ [◦ ]∆t [s]

Event Type

1st

10.26
37.74

1.18
1.90

−8.64 × 103
9.50 × 103

Through
Showering

2nd

7.04
9.29

1.26
0.57

−3.46 × 103
2.59 × 103

ThroughThrough

3rd

0.78
8.99
2.67
9.66
4.45

5.90
11.01
2.99
1.18
5.26

−2.37 × 105−1.37 × 105
2.25 × 104
2.94 × 104
1.52 × 105

Through
Showering
Stop
Showering
Stop

75

Figure 6.9: The histogram of post-trial T Smax distribution for timedependent search obtained from 1,000 trials with solely background toy
MC.
time (2017.09.22). Then we compared the alerts of high-energy neutrino in IceCube
catalogue[2] and 27 potential candidates taken from [21], as shown in Figure 6.13,
Table 6.4 and Table 6.5, and still found no match. Finally, we checked the catalogs
and astronomical archives of HEASARC[1].

6.3

DISCUSSION

In comparison with the previous time-integrated search for astrophysical neutrino
point sources from SK-I to SK-III in [58], this search additionally included the energy distribution in the likelihood to consider different power-law emission spectrawith varying spectral index. The sensitivity 3.63 × 10−7 [cm−2 s−1 averaged in the
declination region (−30◦ , 0◦ ) is comparable to 3 × 10−7 cm−2 s−1 in previous search.
There is no events or celestial objects that has an obvious coincidence with the position or time of the largest T S in either time-integrated search and time-dependent

76

(a)

(b)

(c)

Figure 6.10: Histograms of pre-trial T S distribution for timedependent search of 1st (top), 2nd (middle) and 3rd (bottom). Eachone is simulated 106 times at the corresponding direction.

77

(a)

(b)

(c)

Figure 6.11: L(nS |γ̂, t̂0 , σ̂t ) as a function of nS with γ, t0 and σt fixed
to the best fitted value γ̂, t̂0 and σ̂t , for 1st (top), 2nd (middle) and 3rd(bottom). The vertical line represents the upper limit for 90% CL.

78

(a)

(b)

(c)

(d)

(e)

(f)

Figure 6.12: Space orientation (left column) and time distribution(right column) of the events within 15◦ to the position of the 1st (top),
2nd (middle) and 3rd (bottom) largest T S with in 7σ̂t around t̂0 . Blackdashed circle in the left figures is the boundary of the search cone. The
red plus indicates the best fitted source position. Blue cross is the
direction of the UPMU event, with the circle around them being theangular resolution. The black dashed line in the right figures is theGaussian distribution with the best fitted t0 and σt as the parameter.

79

(a)

(b)

Figure 6.13: T S sky map for time-integrated (top) and timedependent (bottom) search, with the the green squares represent theevents from IceCube alert catalog (Table 6.4) and yellow diamonds
represent the candidates from list in Table 6.5, and the red dot is thelocation of TXS 0506+056.

80

Table 6.4: List of the IceCube alert.

UTCra [◦ ]

dec [◦ ]UTCra [◦ ]

dec [◦ ]2010-09-12 01:41:53
2011-02-16 19:40:59
2011-12-28 12:45:49
2013-04-08 04:31:49
2013-12-02 13:39:09
2014-04-20 01:35:28
2014-12-09 03:26:05
2015-09-11 13:34:02
2016-04-27 05:52:32
2016-11-03 09:07:31
2017-05-06 12:36:56
2010-10-28 07:16:47
2011-03-04 22:54:57
2012-03-01 19:21:57
2012-09-22 13:11:01
2013-09-07 19:02:03
2013-12-04 11:16:53
2014-01-09 12:04:17
2014-06-10 16:54:23
2015-08-12 18:13:05
2015-09-23 00:38:35
2015-11-22 12:45:31
2016-05-10 15:56:05
2016-08-06 12:21:33
2017-03-21 07:32:21

127.9
182.4
345.6
167.17206.63238.98160.05240.20240.57
40.83
221.8
88.68
116.37238.0170.75
129.81
289.16292.85110.30328.19103.27262.18352.34122.81
98.30

-31.2
-21.2
-24.8
20.67
-22.02-37.73
6.57
-0.45
9.34
12.56
-26.0
0.46
-10.72
18.60
19.79-10.36-14.25
33.06
11.57
6.21
3.88
-2.38
2.09
-0.81
-15.02

2010-11-12 13:14:21
2011-07-14 02:42:41
2012-05-01 13:41:24
2013-06-27 02:38:57
2014-01-22 04:54:26
2014-11-09 04:58:36
2015-04-28 11:20:46
2015-12-07 10:36:57
2016-08-14 21:45:54
2017-03-12 13:49:39
2010-10-09 09:08:02
2011-01-28 01:30:24
2011-09-30 10:41:00
2012-05-23 13:46:58
2012-10-11 18:30:09
2013-10-23 13:24:14
2014-01-08 07:23:22
2014-02-03 18:50:34
2014-09-23 17:18:22
2015-08-31 05:13:40
2015-09-26 21:37:43
2016-01-28 04:24:11
2016-07-31 01:55:03
2016-12-10 20:06:41
2017-09-22 20:54:30

110.56
67.86
164.8
93.43
219.64
55.63
80.77
0.4
200.3
305.15331.09307.53266.48171.03205.22301.82344.53349.54169.7254.85
194.50263.40
214.5
46.58
77.43

-0.37
40.32
-71.5
14.02
-86.16-16.50-20.757.5-32.4-26.61
11.10
1.19
-4.41
26.36
-2.39
11.49
1.57
-13.71
-1.34
33.96
-4.34
-14.79
-0.33
14.98
5.72

81

Table 6.5: List of potential astrophysical neutrino source candidates[21]. The column from left is name, declination, right ascension,
distance measured in z (redshift) or kpc, observed gamma ray flux,
energy threshold for flux integration, gamma ray emission energy spectrum index[60], and neutrino flux 90 CL upper limit at that direction
for time-integrated search (averaged for total live time of SK-IV).
Namedec [◦ ]ra [◦ ]

distanceΦγ [×10−12 cm−2 s−1 ]Eth [TeV]γ

Φ90CL[×10−7 cm−2 s−1 ]ν

1ES0347-121
RBS0723H2356-309
PKS2005-489
PKS0548-322
RGBJ0152+017
RGBJ2243+203
VERJ0521+211
PKS1440-389
1ES1215+303
PKS1424+240
PKS2155-304
Markarian421
1ES1218+304
PSRB1259-63
Terzan5
Crab
RCW86
MSH15-52
SNRG327.1-01.1
RXJ0852.0-4622
RXJ1713.7-3946
W28
SNRG015.4+00.1
HESSJ1912+101
W51C
IC443

-11.9911.56-30.63-48.82-32.27
1.79
20.35
21.21
-39.1430.10
23.79-30.2238.19
30.19-63.83-24.9022.01-62.48-59.16-55.08-46.37-39.75-23.34-15.47
10.15
14.19
22.50

57.35
131.80359.78302.3787.67
28.17
340.9880.44
220.99184.45216.75329.72166.08185.36
195.70266.95
83.63
220.68228.53238.65133.00258.25270.43274.52288.21290.7594.21

z=0.19z=0.20z=0.17z=0.07z=0.07z=0.08z=0.14z=0.13z=0.12z=0.03z=0.182.7 kpc
5.9 kpc
2.0 kpc
2.5 kpc
5.2 kpc
9.0 kpc
0.2 kpc
1.0 kpc
2.0 kpc
4.8 kpc
5.4 kpc
1.5 kpc

3.9[18]2.6[14]4.1[17]2.6[15]2.7[16]2.7[19]78[10]
19[23]
19[7]8.0[22]
21[24]
43[11]
310[20]12[13]
0.4[6]
0.3[6]
21[40]
2.4[6]
5.7[6]
0.29[6]23[6]17[6]0.45[6]0.23[6]2.5[6]
0.69[6]4.6[12]

0.25
0.2
0.2
0.4
0.25
0.3
0.112
0.2
0.147
0.2
0.12
0.2
0.4
0.2
1
1
1
1
1
1
1
1
1
1
1
1
0.3

3.1
3.6
3.1
3.2
2.9
3.0
4.6
3.4
3.7
3.6
3.8
3.5
2.2
3.1
2.6
2.4
2.5
1.6
2.3
2.2
1.8
2.1
2.7
2.1
2.7
2.6
3.0

1.36
1.39
1.07
1.09
2.15
1.43
1.42
2.19
1.12
3.03
1.40
1.96
3.21
3.11
2.48
1.72
3.28
1.58
1.04
1.13
2.24
1.34
1.54
0.95
1.02
1.29
1.58

82
search. For the direction of the second largest in time-dependent search, the closest
event both in space and in time is a GRB numbered as GRB090428552 detected
157.9′ away (at (0.8◦ , 11.5◦ )) and about 5 month earlier (on 2009.04.28). The GRB
typically last from several seconds to several thousand seconds, and no evidence hassuggested that the neutrino emission from the GRB has a time delay compared withgamma ray emission up to several month, so we just keep this as a record for future
study. As for the third one, a known pulsar PSR J1833-6023 is located within 18.5’
away from it at (278.3◦ , −60.4◦ ). It is possible that these are related, but further
study regarding this is needed to confirm it.

83
CHAPTER 7
Summary and Prospects
In this study, we performed both time-integrated search and time-dependent searchin the full sky for high-energy (above 1.6 GeV) astrophysical neutrino point sourcesin SK-IV with UPMU. Unfortunately, no excess with respect to the background isfound in either search.
The most significant cluster in the full sky for the time-integrated search is at(198.3◦ , 15.1◦ ) in J2000 equatorial coordinate system with pre-trial p-value of 1.18 ×
10−4 and post-trial p-value of 79.2% . The upper limit of the neutrino fluence isΦ0 = 2.06 × 102 [GeVcm−2 ] and corresponding flux is φupper = 4.53 × 10−7 cm−2 s−1 .
For the time-dependent search, the most significant cluster is at (118.0◦ , −62.7◦ )on 2010.10.27 13:09:42 JST and σ̂t = 9.25 × 103 s. The pre-trial p-value is 2.24 × 10−4
and post-trial p-value is 33.4%. Neutrino fluence upper limit is set to Φ0 = 7.86 ×
101 [GeVcm−2 ], with the emission time interval estimated by 1.85 × 104 s (2σ̂t ).
There are still a lot of work that could be done which is related to this study.
Although not many, part of the FC and PC data still have a good-enough angularresolution which could be used for the search. And till now the SK-V data is completely ready and a joint search from SK-I to SK-V can be performed. Except for thefull sky search, trigger based search could also help to reduce the huge trial factor
raised when looking for too many directions at a time.
In addition to the present experiment facility, the next-generation neutrino detector Hyper-Kamiokande (HK) has already been constructed. With larger water tankand more PMTs, we can expect a better detection performance in the future.

84

BIBLIOGRAPHY[1] Nasa’s heasarc: Archive. https://heasarc.gsfc.nasa.gov/docs/archive.html.
[2] Icecube
catalogofalert
eventsupthrough
icecube170922a.
https : / / icecube . wisc . edu / data-releases / 2018 / 07 /icecube-catalog-of-alert-events-up-through-icecube-170922a/, 2018.
[3] Healpix. https://healpix.sourceforge.io/, 2019.
[4] M. Aartsen, M. Ackermann, J. Adams, J. Aguilar, M. Ahlers, M. Ahrens, D. Altmann, T. Anderson, M. Archinger, C. Arguelles, et al. Searches for timedependent neutrino sources with icecube data from 2008 to 2012. The Astrophysical Journal, 807(1):46, 2015.
[5] R. Abbasi, M. Ackermann, J. Adams, J. Aguilar, M. Ahlers, M. Ahrens, C. Alispach, A. Alves, N. Amin, K. Andeen, et al. A search for time-dependent
astrophysical neutrino emission with icecube data from 2012 to 2017. The Astrophysical Journal, 911(1):67, 2021.
[6] H. Abdalla, A. Abramowski, F. Aharonian, F. A. Benkhali, E. Angüner,M. Arakawa, M. Arrieta, P. Aubert, M. Backes, A. Balzer, et al. The hessgalactic plane survey. Astronomy & Astrophysics, 612:A1, 2018.
[7] H. Abdalla, R. Adam, F. Aharonian, F. Ait Benkhali, E. Angüner, M. Arakawa,
C. Arcaro, C. Armand, T. Armstrong, H. Ashkar, et al. Very high energy γ-ray
emission from two blazars of unknown redshift and upper limits on their distance.
Monthly Notices of the Royal Astronomical Society, 494(4):5590–5602, 2020.
[8] K. Abe, C. Bronner, Y. Haga, Y. Hayato, M. Ikeda, K. Iyogi, J. Kameda,
Y. Kato, Y. Kishimoto, L. Marti, et al. Atmospheric neutrino oscillation analysis with external constraints in super-kamiokande i-iv. Physical Review D,
97(7):072001, 2018.
[9] K. Abe, C. Bronner, G. Pronost, Y. Hayato, M. Ikeda, K. Iyogi, J. Kameda,
Y. Kato, Y. Kishimoto, L. Marti, et al. Search for nucleon decay into charged antilepton plus meson in 0.316 megaton· years exposure of the super-kamiokande
water cherenkov detector. Physical Review D, 96(1):012003, 2017.
[10] A. Abeysekara, S. Archambault, A. Archer, W. Benbow, R. Bird, R. Brose,
M. Buchovecky, J. Buckley, V. Bugaev, M. Cerruti, et al. Discovery of veryhigh-energy emission from rgb j2243+ 203 and derivation of its redshift upper
limit. The Astrophysical Journal Supplement Series, 233(1):7, 2017.

85[11] A. Abramowski, F. Acero, F. Aharonian, A. Akhperjanian, G. Anton, U. B.
De Almeida, A. Bazer-Bachi, Y. Becherini, B. Behera, W. Benbow, et al. Vhe
γ-ray emission of pks 2155–304: spectral and temporal variability. Astronomy &Astrophysics, 520:A83, 2010.
[12] V. Acciari, E. Aliu, T. Arlen, T. Aune, M. Bautista, M. Beilicke, W. Benbow,
S. Bradbury, J. Buckley, V. Bugaev, et al. Observation of extended very high energy emission from the supernova remnant ic 443 with veritas. The AstrophysicalJournal, 698(2):L133, 2009.
[13] V. Acciari, E. Aliu, T. Arlen, M. Beilicke, W. Benbow, S. Bradbury, J. Buckley,
V. Bugaev, Y. Butt, K. Byrum, et al. Veritas observations of the bl lac object
1es 1218+ 304. The Astrophysical Journal, 695(2):1370, 2009.
[14] V. A. Acciari, S. Ansoldi, L. Antonelli, A. A. Engels, K. Asano, D. Baack,A. Babić, B. Banerjee, U. B. de Almeida, J. Barrio, et al. New hard-tev extreme blazars detected with the magic telescopes. The Astrophysical Journal
Supplement Series, 247(1):16, 2020.
[15] F. Acero, F. Aharonian, A. Akhperjanian, G. Anton, U. B. De Almeida, A. BazerBachi, Y. Becherini, B. Behera, W. Benbow, K. Bernlöhr, et al. Pks 2005-489
at vhe: four years of monitoring with hess and simultaneous multi-wavelength
observations. Astronomy & Astrophysics, 511:A52, 2010.
[16] F. Aharonian, A. Akhperjanian, G. Anton, U. B. De Almeida, A. Bazer-Bachi,
Y. Becherini, B. Behera, W. Benbow, K. Bernlöhr, A. Bochow, et al. Discovery ofvhe γ-rays from the bl lacertae object pks 0548–322. Astronomy & Astrophysics,
521:A69, 2010.
[17] F. Aharonian, A. Akhperjanian, A. Bazer-Bachi, M. Beilicke, W. Benbow,
D. Berge, K. Bernlöhr, C. Boisson, O. Bolz, V. Borrel, et al. Discovery ofvery high energy γ-ray emission from the bl lacertae object h 2356-309 with thehess cherenkov telescopes. Astronomy & Astrophysics, 455(2):461–466, 2006.
[18] F. Aharonian, A. Akhperjanian, U. B. De Almeida, A. Bazer-Bachi, B. Behera,
M. Beilicke, W. Benbow, K. Bernlöhr, C. Boisson, O. Bolz, et al. Discovery ofvhe γ-rays from the distant bl lacertae 1es 0347-121. Astronomy & Astrophysics,
473(3):L25–L28, 2007.
[19] F. Aharonian, A. Akhperjanian, U. B. De Almeida, A. Bazer-Bachi, B. Behera,
M. Beilicke, W. Benbow, K. Bernlöhr, C. Boisson, V. Borrel, et al. Discoveryof vhe γ-rays from the high-frequency-peaked bl lacertae object rgb j0152+ 017.
Astronomy & Astrophysics, 481(3):L103–L107, 2008.

86[20] M. Ahnen, S. Ansoldi, L. Antonelli, P. Antoranz, A. Babic, B. Banerjee, P. Bangale, U. B. De Almeida, J. Barrio, J. B. González, et al. Long-term multiwavelength variability and correlation study of markarian 421 from 2007 to 2009.
Astronomy & Astrophysics, 593:A91, 2016.
[21] A. Albert, M. André, M. Anghinolfi, G. Anton, M. Ardid, J.-J. Aubert, T. Avgitas, B. Baret, J. Barrios-Martí, S. Basa, et al. First all-flavor neutrino pointlike source search with the antares neutrino telescope. Physical Review D,
96(8):082001, 2017.
[22] E. Aliu, S. Archambault, T. Arlen, T. Aune, M. Beilicke, W. Benbow, R. Bird,
A. Bouvier, J. Buckley, V. Bugaev, et al. Long term observations of b2 1215+
30 with veritas. The Astrophysical Journal, 779(2):92, 2013.
[23] S. Archambault, T. Arlen, T. Aune, B. Behera, M. Beilicke, W. Benbow, R. Bird,
A. Bouvier, J. Buckley, V. Bugaev, et al. Discovery of a new tev gamma-ray
source: Ver j0521+ 211. The Astrophysical Journal, 776(2):69, 2013.
[24] S. Archambault, T. Aune, B. Behera, M. Beilicke, W. Benbow, K. Berger,
R. Bird, J. Biteau, V. Bugaev, K. Byrum, et al. Deep broadband observations of the distant gamma-ray blazar pks 1424+ 240. The Astrophysical Journal
Letters, 785(1):L16, 2014.
[25] Y. Ashie, J. Hosaka, K. Ishihara, Y. Itow, J. Kameda, Y. Koshio, A. Minamino,
C. Mitsuda, M. Miura, S. Moriyama, et al. Measurement of atmospheric neutrino
oscillation parameters by super-kamiokande i. Physical Review D, 71(11):112005,2005.
[26] G. Barr, T. Gaisser, P. Lipari, S. Robbins, and T. Stanev. Three-dimensional
calculation of atmospheric neutrinos. Physical Review D, 70(2):023006, 2004.
[27] G. Battistoni, A. Ferrari, T. Montaruli, and P. Sala. The fluka atmospheric
neutrino flux calculation. Astroparticle Physics, 19(2):269–290, 2003.
[28] P. Blasi, R. I. Epstein, and A. V. Olinto. Ultra-high-energy cosmic rays fromyoung neutron star winds. The Astrophysical Journal, 533(2):L123, 2000.
[29] J. Braun, M. Baker, J. Dumm, C. Finley, A. Karle, and T. Montaruli. Timedependent point source search methods in high energy neutrino astronomy. Astroparticle Physics, 33(3):175–181, 2010.
[30] I. Collaboration, M. Aartsen, M. Ackermann, J. Adams, J. A. Aguilar, M. Ahlers,
M. Ahrens, I. Al Samarai, D. Altmann, K. Andeen, et al. Neutrino emission fromthe direction of the blazar txs 0506+ 056 prior to the icecube-170922a alert.
Science, 361(6398):147–151, 2018.

87[31] I. Collaboration, Fermi-LAT, MAGIC, AGILE, ASAS-SN, HAWC, HESS, INTEGRAL, Kanata, Kiso, et al. Multimessenger observations of a flaring blazar coincident with high-energy neutrino icecube-170922a. Science, 361(6398):eaat1378,
2018.
[32] S. Desai. High energy neutrino astrophysics with Super-Kamiokande. Boston
University, 2004.
[33] P. Diamond and M. Malkov. Dynamics of mesoscale magnetic field in diffusive
shock acceleration. The Astrophysical Journal, 654(1):252, 2007.
[34] E. Fermi. On the origin of the cosmic radiation. Physical review, 75(8):1169,
1949.
[35] S. Fukuda, Y. Fukuda, T. Hayakawa, E. Ichihara, M. Ishitsuka, Y. Itow, T. Kajita, J. Kameda, K. Kaneyuki, S. Kasuga, et al. The super-kamiokande detector.
Nuclear Instruments and Methods in Physics Research Section A: Accelerators,
Spectrometers, Detectors and Associated Equipment, 501(2-3):418–462, 2003.
[36] M. Glück, E. Reya, and A. Vogt. Dynamical parton distributions revisited. The
European Physical Journal C-Particles and Fields, 5(3):461–470, 1998.
[37] K. Greisen. End to the cosmic-ray spectrum?
16(17):748, 1966.

Physical Review Letters,[38] P. D. Group et al. Review of particle physics. Chinese Physics C, 40(10):100001,
2016.
[39] Y. Hayato. A neutrino interaction simulation program library neut. Acta Physica
Polonica B, 40(9), 2009.
[40] A. Hillas, C. Akerlof, S. Biller, J. Buckley, D. Carter-Lewis, M. Catanese,
M. Cawley, D. Fegan, J. Finley, J. Gaidos, et al. The spectrum of tev gamma
rays from the crab nebula. The Astrophysical Journal, 503(2):744, 1998.
[41] M. Honda, T. Kajita, K. Kasahara, and S. Midorikawa. Calculation of the fluxof atmospheric neutrinos. Physical Review D, 52(9):4985, 1995.
[42] M. Honda, T. Kajita, K. Kasahara, and S. Midorikawa. Improvement of low
energy atmospheric neutrino flux calculation using the jam nuclear interaction
model. Physical Review D, 83(12):123001, 2011.
[43] G. Illuminati, J. Aublin, and S. Navas. Searches for point-like sources of cosmic
neutrinos with 11 years of antares data. In 36th International Cosmic Ray
Conference (ICRC2019), volume 36, page 920, 2019.

88[44] Y. Itow, T. Kajita, K. Kaneyuki, M. Shiozawa, Y. Totsuka, Y. Hayato, T. Ishida,
T. Ishii, T. Kobayashi, T. Maruyama, et al. The jhf-kamioka neutrino project.
arXiv preprint hep-ex/0106019, 2001.
[45] R. Jaffe and C. L. Smith. Deep-inelastic scattering, the subtraction of divergent
sum rules, and chiral-symmetry breaking in the gluon model. Physical Review
D, 7(8):2506, 1973.
[46] J. Kameda. Detailed studies of neutrino oscillations with atmospheric neutrinos
of wide energy range from 100 MeV to 1000 GeV in Super-Kamiokande. PhD
thesis, 2002.
[47] B. Louis, V. Sandberg, G. Garvey, H. White, G. Mills, and R. Tayloe. The
evidence for oscillations. Los Alamos Sci, 25:16, 1997.
[48] P. F. Menéndez. Neutrino Physics in Present and Future Kamioka WaterČerenkov Detectors with Neutron Tagging. Springer, 2018.
[49] P. Meszaros. Gamma-ray bursts. Reports on Progress in Physics, 69(8):2259,
2006.
[50] M. Nakahata, K. Arisaka, T. Kajita, M. Koshiba, Y. Oyama, A. Suzuki,
M. Takita, Y. Totsuka, T. Kifune, T. Suda, N. Sato, K. Takahashi, andK. Miyano. Atmospheric neutrino background and pion nuclear effect for
kamioka nucleon decay experiment. Journal of the Physical Society of Japan,
55(11):3786–3805, 1986.
[51] J. H. NASA, ESA and A. L. A. S. University). Messier 1 (the crab nebula).
https://www.nasa.gov/feature/goddard/2017/messier-1-the-crab-nebula, 2021.
[52] H. Nishino, K. Awai, Y. Hayato, S. Nakayama, K. Okumura, M. Shiozawa,
A. Takeda, K. Ishikawa, A. Minegishi, and Y. Arai. High-speed charge-to-time
converter asic for the super-kamiokande detector. Nuclear instruments and methods in physics research section a: accelerators, spectrometers, detectors and associated equipment, 610(3):710–717, 2009.
[53] A. Panov. Electrons and positrons in cosmic rays. In Journal of Physics: Conference Series, volume 409, page 012004. IOP Publishing, 2013.
[54] D. Rein and L. M. Sehgal. Neutrino-excitation of baryon resonances and single
pion production. Annals of Physics, 133(1):79–153, 1981.
[55] C. Saji. Study of upward-going muons in Super-Kamiokande. PhD thesis, Niigata
University, 2002.

89[56] R. Smith and E. J. Moniz. Neutrino reactions on nuclear targets. Nuclear Physics
B, 43:605–622, 1972.
[57] F. W. Stecker, C. Done, M. H. Salamon, and P. Sommers. Erratum:‘‘highenergy neutrinos from active galactic nuclei’’[phys. rev. lett. 66, 2697 (1991)].
Physical Review Letters, 69(18):2738, 1992.
[58] E. Thrane, K. Abe, Y. Hayato, T. Iida, M. Ikeda, J. Kameda, K. Kobayashi,
Y. Koshio, M. Miura, S. Moriyama, et al. Search for astrophysical neutrino point
sources at super-kamiokande. The Astrophysical Journal, 704(1):503, 2009.
[59] C. M. Urry and P. Padovani. Unified schemes for radio-loud active galactic nuclei.
Publications of the Astronomical Society of the Pacific, 107(715):803, 1995.
[60] S. P. Wakely and D. Horan. Tevcat: an online catalog for very high energy
gamma-ray astronomy. In International Cosmic Ray Conference, volume 3, pages
1341–1344, 2008.
[61] E. Waxman. Cosmological gamma-ray bursts and the highest energy cosmic rays.
Physical Review Letters, 75(3):386, 1995.
[62] R. J. Wilkes. Results on neutrino oscillations from super-kamiokande. Advancesin Space Research, 26(11):1813–1822, 2000.
[63] G. T. Zatsepin and V. A. Kuz’min. Upper limit of the spectrum of cosmic rays.
Soviet Journal of Experimental and Theoretical Physics Letters, 4:78, 1966.
[64] C. Zeitnitz and T. Gabriel. The geant-calor interface and benchmark calculationsof zeus test calorimeters. Nuclear Instruments and Methods in Physics Research
Section A: Accelerators, Spectrometers, Detectors and Associated Equipment,349(1):106–111, 1994.