Study of neutrino signal from astronomical object
in Super-Kamiokande

2019, March

Chenyuan Xu

Graduate School of
Natural Science and Technology(Doctor’s Course)
OKAYAMA UNIVERSITY

1

Abstract
This thesis focus on the study of neutrino signal search in SuperKamiokande(SK), the largest water Cherenkov detector in the world. This
thesis includes two topics: one is the coincident neutrino search in SK
with the first detected gravitational wave produced by a binary neutronstar merger, GW170817; another is the study of the de-excitation gamma
signal from Charged Current Quasi-Elastic(CCQE) interaction of atmospheric neutrino in SK.
For the first topic, we searched for coincident neutrino events in therange from 4 MeV to 100 PeV, in a time window of ±500 seconds aroundthe gravitational wave detection time, as well as during a 14 day period
after the detection. No significant neutrino signal was observed for either
time window. The upper limits on the neutrino fluence for GW170817
were calculated at 90% confidence level. From the upward-going-muon
events in the energy region from 1.6 GeV to 100 PeV, the neutrino flu+1.1
2◦
ence limit is 16.0+0.7−0.6 (21.3−0.8 )cm for νµ (ν̄µ ), with an error range of ±5
around the zenith angle of the source of GW170817, and the energy spectrum assumed as an index of -2. The fluence limit for neutrino energies
less than 100 MeV, for which the emission mechanism would be different
than for higher-energy neutrinos, is also calculated. The best limit is foranti-electron neutrinos under the assumption of a Fermi-Dirac spectrumwith average energy of 20 MeV.
The second topic focus on an analysis method to directly measure thebranching ratio of de-excitation gamma in CCQE interaction of atmospheric muon neutrino on Oxygen nucleus in water Cherenkov detector.
The branching ratio of de-excitation gamma has not been directly measured yet, which causes not only the unreducable background component
in Supernova Relic Neutrino search but also the systematic error in long
baseline neutrino experiment. In this thesis, details about the method
will be introduced, and the result of the direct search for de-excitation
gamma in SK-IV real data will also be shown.

2

Acknowledgement
First of all, I would like to express my sincere gratitude to my supervisor
Prof. Y. Koshio for giving me the opportunity to studying neutrino physics at
Super-Kamiokande. This thesis would never exist without his kind support andexcellent guidance.
I would like to appreciate Prof. M. Nakahata for giving me the opportunity
to participate in the Super-Kamiokande experiment.
I would like to extend my gratitude to Prof. Y. Takeuchi, Dr. M. Smy,
Prof. H. Sekiya, Prof. M. Vagins, Prof. L. Labarga, Prof. K. Okumura, Dr. L.
Marti, Dr. M. Ikeda, Dr. T. Yano, Dr. E. O’Sullivan, Dr. Y. Nakajima, Dr.
G. Pronost, Dr. P. Weatherly, Dr. S. Cao. They gave me great help and usefuladvice when I was working in LOWE and EGADS group.
I would like to appreciate Prof. M. Sakuda and Prof. H. Ishino, for theirkind guidance in the physics leature in Okayama University.
I would like to appreciate Prof. J. Beacom, who support me a lot in bothlife and research when I was in Ohio State University.
I would like to thank to laboratory members or friends in Okayama University, Dr. S. Ito, D. Fukuda, K. Hagiwara, M. Harada, A. Hourai, H. Nagata, T.
Hirashige, T. Kayano, Y. Wang, T. Mori, who gave me great help and support
in my 5 years’ life of PhD course.
Finally, I would like to express my greatest appreciation and gratitude to
my family, relatives, and all my other friends.

3

Contents
1 Physics Background
1.1 Neutrino signal in Super-Kamiokande . . . . . . . . . . . . . . .
1.2 Binary Neutron Star merger . . . . . . . . . . . . . . . . . . . . .
1.2.1 Neutrino Emission . . . . . . . . . . . . . . . . . . . . . .
1.2.2 Expected Neutrino Number in Detector . . . . . . . . . .
1.3 Solar Neutrino . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.4 Supernova Neutrino and Supernova Relic Neutrino . . . . . . . .
1.4.1 Supernova Explosion . . . . . . . . . . . . . . . . . . . . .
1.4.2 Supernova Relic Neutrino . . . . . . . . . . . . . . . . . .
1.5 Quasi-Elastic neutrino reaction and de-excitation gamma . . . .

6
6
6
6
7
9
10
10
12
13

2 Super-Kamiokande detector
15
2.1 Detector Structure . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.1.1 History and current status . . . . . . . . . . . . . . . . . . 16
2.2 Dection method . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.2.1 Cherenkov radiation . . . . . . . . . . . . . . . . . . . . . 16
2.2.2 Photomultiplier tube . . . . . . . . . . . . . . . . . . . . . 18
2.2.3 Neutrino interaction in water Cerenkov detector . . . . . 22
2.3 Data acquisition system . . . . . . . . . . . . . . . . . . . . . . . 23
2.3.1 Electronics . . . . . . . . . . . . . . . . . . . . . . . . . . 23
2.3.2 Trigger . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
3 Event Reconstruction25
3.1 low energy event . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
3.1.1 Vertex reconstruction . . . . . . . . . . . . . . . . . . . . 25
3.1.2 Direction reconstruction . . . . . . . . . . . . . . . . . . . 26
3.2 Mu event . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
4 Simulation30
4.1 SKdetSIM overview . . . . . . . . . . . . . . . . . . . . . . . . . 30
4.2 Particle generation . . . . . . . . . . . . . . . . . . . . . . . . . . 30
4.3 Particle tracking and Cherenkov photon emission . . . . . . . . . 30
4.4 Progration of Cherenkov photon . . . . . . . . . . . . . . . . . . 31
4.5 PMT response and electronics simulation . . . . . . . . . . . . . 32
4.6 Summary of SKdetsim . . . . . . . . . . . . . . . . . . . . . . . . 35
4.6.1 Tunable Parameters . . . . . . . . . . . . . . . . . . . . . 35
4.6.2 Some other explanation . . . . . . . . . . . . . . . . . . . 36
5 Calibration
37
5.1 Xe and Ni calibration . . . . . . . . . . . . . . . . . . . . . . . . 37
5.2 LINAC Calibration . . . . . . . . . . . . . . . . . . . . . . . . . . 38
5.3 DT Calibration . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
5.4 Water Transperency . . . . . . . . . . . . . . . . . . . . . . . . . 41
5.4.1 By Laser . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
5.4.2 By decay electron . . . . . . . . . . . . . . . . . . . . . . 44

4

6 Data Set
46
6.1 Solar neutrino sample and reduction . . . . . . . . . . . . . . . . 46
6.1.1 Bad run selection and bad channel cut . . . . . . . . . . . 46
6.1.2 Pre reduction . . . . . . . . . . . . . . . . . . . . . . . . . 46
6.1.3 Spallation cut . . . . . . . . . . . . . . . . . . . . . . . . . 47
6.1.4 16 N cut and other cut . . . . . . . . . . . . . . . . . . . . 52
6.2 Relic neutrino sample and reduction . . . . . . . . . . . . . . . . 54
6.2.1 Spallation cut . . . . . . . . . . . . . . . . . . . . . . . . . 54
6.2.2 2peak cut . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
6.2.3 Multiple ring cut . . . . . . . . . . . . . . . . . . . . . . . 59
6.2.4 Cherenkov angle cut . . . . . . . . . . . . . . . . . . . . . 60
6.2.5 Pion cut . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
6.3 Atmopheric neutrino sample . . . . . . . . . . . . . . . . . . . . . 62
7 Neutrino search for GW170817
63
7.1 GW170817 overview . . . . . . . . . . . . . . . . . . . . . . . . . 63
7.2 Search method and result . . . . . . . . . . . . . . . . . . . . . . 64
7.3 Fluence limit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
8 De-excitation gamma search from Charge Current Quasi-Elastic
reaction
74
8.1 Motivation and Plan . . . . . . . . . . . . . . . . . . . . . . . . . 74
8.2 Search Method for Branching Ratio of de-excitation gamma . . . 75
8.2.1 Simulation setup . . . . . . . . . . . . . . . . . . . . . . . 75
8.2.2 Light profile of mu and gamma . . . . . . . . . . . . . . . 76
8.2.3 Building Likelihood . . . . . . . . . . . . . . . . . . . . . 79
8.3 Direct measurement of de-excitation gamma in Super-K . . . . . 81
9 Conclusions and Future88
9.1 Neutrino search associated with GW170817 . . . . . . . . . . . . 88
9.2 De-excitation gamma from atmospheric neutrino CCQE interaction 89

5

1

Physics Background

1.1

Neutrino signal in Super-Kamiokande

Super-Kamiokande(SK) is a water Cherenkov detector. The purpose ofSuper-Kamiokande is various neutrino physics and the measurement of nucleon
life time. It observed solar neutrinos [1], atmospheric neutrinos [49] as well as
accelerator neutrinos [50], and several evidences for neutrino oscillations havebeen found. SK has been collecting data for more than twenty years and it
is still the water Cherenkov detector with best sensitivity for MeV and GeV
neutrino in the world.
This thesis is focused on two research topics. One is the coincident neutrino
search in SK with the first detected gravitational wave produced by a binary neutron star merger, GW170817; another is the study of the de-excitation gamma
signal from Charged Current Quasi-Elastic(CCQE) interaction of atmospheric
neutrino in SK.
For the first topic, the mechanism of neutrino emission from binary starmerger will be presented in this section. Especially for this research, two kindsof event reduction method are used for neutrino signal search below 100MeV:
supernova relic neutrino search method and solar neutrino search method. Inorder to explain the purpose of the two, the related physics background will also
be briefly introduced in this section, and the details about reduction steps will
be in Section 6.
For the second topic, the essential purpose for this research is to estimate akind of unreducable background in supernova relic neutrino search. This kindof background results from the atmospheric neutrino, and it not only be themost difficulty in current search but also remains unreducable in the next phaseof SK. The motivation will be explained in this section, while the research plan,
analysis status and current result will be shown in Section 8.

1.2

Binary Neutron Star merger

1.2.1

Neutrino Emission

On August 17th 2017 at 12:41:04 UTC, the Advanced LIGO and Advanced
Virgo experiment identified the first evident signal of a gravitational wave fromthe binary neutron star merger, named GW170817 [2]. The interpretation is amerger of two compact objects consistent with neutron stars having total system
mass of 2.74 solar masses and a luminosity distance of 40 Mpc.
The discovery of GW170817 marked not only a starting point of gravitional
wave astronomy with binary neutron-star mergers, but also the possibility ofmultimessenger observation which can lead us to study the equation of state ofsupranuclear density matter [2]. Different to the gravitational wave events frombinary black hole merger(such as the first detected gravitational wave event
GW150914 [40]), the one from binary neutron star merger can enable the signal
detection in electromagnetic channels [41]. In fact, associated observation ofshort gamma ray busrt GRB170817 as well as the brightening X-Ray emission
are already been suggested [1, 42].
Since the high-density material reaches several tens of MeV via the binary
neutron star merger, thermal neutrino emission could be expected from theremnant [7], unless the possibility that the remnant finally result into a black6

hole with prompt collapse, but fortunatelly this is not likely to happen whenthe maximum mass of spherical neutron stars exceeding ∼ 2M [43].
Direct detections of thermal neutrinos could be an important step in multimessenger astronomy, as the theoretical models of supernova explosions havebeen qualitatively confirmed by the detection of neutrinos from SN1987A [44].
However, little has been known about the realistic spectrum of neutrinos frombinary neutron-star mergers so far. Monte-Carlo neutrino-transport simulations suggest that the spectrum can be approximated by pinched Fermi-Dirac
distribution for as the case of supernova explosions [45].
The mechanism of neutrino emission from binary neutron star can be considered to follow these steps:
1. As the matter temperature getting high after the merger, pairs of e− ande+ are produced from thermal photons.
2. e− and e+ are captured on nucleons, from which νe and ν e are emittedwith a rise time of ∼ 10ms Fig 1.
e− + p → νe + n(1)

e+ + n → ν e + p(2)

In the case of binary neutron star merger, obviously the number of neutron
is much higher than proton in neutron star, thus ν e is emitted much morethan νe .
3. All neutrino type are made from e− e+ annihillation as:
e+ + e− → νx + ν x(3)

The peak luminosity of electron antineutrinos reaches 1 ∼ 3 × 1053 ergs−1
with the typical energy of 10 ∼ 30MeV depending on binary parameters andunknown equations of state for supranuclear-density matter.
1.2.2

Expected Neutrino Number in Detector

In Super-Kamiokande, the water Cherenkov detector with best sensitivity inthe world for neutrino signals from MeV to GeV, the most dominating reactionin the expected energy range(10 ∼ 30MeV for neutrino from binary neutron starmerger) is inverse beta decay(IBD) of electron antineutrino ν e (p+ν e → n+e+ ).
The expected number of observed neutrinos for a single merger can be simply
estimated as:
Z tf Z EmaxNν = NTφ(E, t)σ(E)dEdt
ti

(4)

Emin

Here NT is number of protons, which is the target for IBD in water Cherenkov
detector. φ(E, t) is the flux of ν e in unit energy. σ(E) is the cross section of
IBD for a ν e on a proton. ti and tf mean the starting and ending for time
integral, while Emin and Emax represent the integral energy range.

7

Figure 1: The expected neutrino energy after binary neutron star merger [8].
Considering the Fermi-Dirac distribution with temperature T, the averageenergy of neutrinos is hEi ≈ 3.15kB T , while kB T is Boltzmann constant. Ignoring time evolution, φ(E, t) can be given by:
φ(E) =cE2
2π 2 (~c)3 exp[E/kB T ] + 1(5)Where ~ is the reduced Planck constant. By taking the leading-order crosssection from [46]:E)2
10M eV
Equation 4 can be approximately written into:σLO (E) = 9.5 × 10−42 cm2 (

(6)

MT
E∆T
hEi
D
)×(
)(
)(
)−2
1M t 3 × 1052 erg
10M eV 100M pc(7)where MT is detector volume and D is the distance from the source tothe earth. In Equation 7, time integral
R is removed and the total energy of ν eemitted during ∆tobs ≈ 1s as E∆t = Lν dt (Lν is the Luminosity of ν e ), whichis ∼ 3×1052 erg [47]. One reason for this is that the time evolution of luminosity
has not been well understood, especially when ∆tobs > 1s.
The other three parameters in Equation 7 is fE , fSE and fOSC . They
respectively represent for the factor of the fraction of energy range, detection
efficiency and the neutrino oscillation effect. Considering a integral range of
10 ∼ 50MeV, fE is 0.77 for typical energy hEi = 10M eV [8]. And in orderto make an order estimation for neutrino events in Super-Kamiokande, here
we ignore fSE and fOSC (the detection efficiency will be discussed in later
section). Considering the fiducial volume 22.5kton of Super-Kamiokande andNν ≈ 1.0 × 10−3 × fE fSE fOSC (8the luminosity distance of 40 Mpc, the expected order of neutrino events is∼ 10−4 , while siginificant signal could be detectable when the source is in a
distance of ∼ 100kpc. The result of neutrino signal search carried out in SuperKamiokande will be introduced in later section.

1.3

Solar Neutrino

The Standard Solar Model (SSM) is the well-established theoretical model
which has been constructed to explain the stellar evolution of the Sun [51]. In
the core of the sun, the origin of the energy is from the the fusion of four protons
into a 4 He as:
4p → α + 2e− + 2νe + 26.73M eV(8)Actually, reaction 8 is realized by two process: one is pp-chain; the otherone is Carbon-Nitrogen-Oxygen (CNO) cycle. According to SSM, 98.4% of thetotal solar luminosity is produced by pp-chain process while the rest 1.6% is
produced by CNO cycle.
The pp-chain includes reactions of:p + p →2 H + e+ + νe (≤ 0.420M eV, pp)(9)p + e− + p →2 H + νe (1.442M eV, pep)
7

−

7

(10)
7

Be + e → Li + νe (0.861M eV (90%), 0.383M eV (10%), Be)8
3

B →8 Be∗ + e+ + νe (≤ 14.06M eV,8 B)+

He + p → α + e + νe (≤ 18.77M eV, hep)

(11)
(12)
(13)

Figure 2: The pp-chain reactions
While CNO cycle include the reactions of:13

N →13 C + e+ + νe (≤ 1.27M eV )

(14)

15

O →15 N + e+ + νe (≤ 1.73M eV )

(15)

17

F →

17

+

O + e + νe (≤ 1.74M eV )

9

(16)

Figure 3: The CNO cycle reactions
The energy spectrum of different kinds of solar neutrinos predicted by SSM
BP05(OP) [52](a revised version of SSM) is shown in Fig 4. The electron neutrino from the three reactions in CNO cycle, and pp, pep, 7 Be in pp-chain arewith energy lower than the detection threshold of SK. Thus the detectable ones
are 8 B and hep. Since hep flux is too low, the observation of 8 B spectrum
becomes the most reliable way to study SSM.

Figure 4: The solar neutrino energy spectrum predicted by the BP05 (OP).

1.4Supernova Neutrino and Supernova Relic Neutrino

1.4.1

Supernova Explosion

On February 23rd 1987, supernova burst of Large Magellanic Cloud wasobserved. It is a supernova from about 51.4kpc distance from the earth and it is10numbered with SN1987A. It is the first time that human detected the supernova
by neutrinos and that is the sign of the birth of “Neutrino Astronomy”. In SN
1987A, the neutrino detector KamiokaNDE observed 11 neutrinos [44], while
American detector IBM observed 8 neutrinos [53]. The detection of supernova
neutrino is not only a new page in the physics history but also proved that thefundamental model of core collapse supernova is correct. Study about supernova
was based on optical observation in the past and it is limited by a lot of reasons.
Neutrino has offered a new way to approach to the the mechanism of supernova
explosion and cosmic evolution. Now experiments of supernova neutrino searchare being carried out in the world wide.
To explain the mechanism of supernova explosion, the process can be divided
as follows.
Beginning of core collapse
When the mass of the core exceed the Chandrasekhar Limit, the gravitational
collapse will start and the density of the core will get larger and larger. νe will
be released by electron capture(2) by the nuclei and the photodisintegration ofiron nuclei(3) also occurs in this step.
e− + A(N, Z) → νe + A(N + 1, Z − 1)(17)

γ +56 F e → 13α + 4n − 124.4M eV(18)

Both of the two interactions will reduce the degenerate pressure and become
the trigger of a rapid core collapse.
Neutrino trapping
When the density of the core exceeds 3 × 109 g/cm3 , the neutrinos will follows
the coherent scattering(4) with the nuclei and be unable to go outside. This iscalled “neutrino trapping” and the boundary is called “neutrino sphere”.
ν + A(N, Z) → ν + A(N, Z)(19)Core bounce
When core density continue to increase and finally reaches the level of nuclear
density(∼ 1014 g/cm3 ), the shock wave will be repulsed and driven towards
outside.
Neutrino cooling
After the shock wave passed through, the matter outside the core will drop to
the center and ∼ 1053 erg of gravitational energy will be converted into thermalenergy. All flavors of neutrinos will be produced by pair annihilation 20 andbremsstrahlung 21.
e− + e+ → ν + ν̄(20)
N + N → N + N + ν + ν̄(21)ν̄e and νe will also be produce by electron capture 22 and positron capture 23.
e− + p → νe + n
11

(22)

e+ + n → ν̄e + p(23)

Most of the energy are released by neutrinos and the star cools down in about
10s.
Supernova burst
The shock wave may takes several hours until it arrives the surface of the star.
The penumbra part will be blasted off and the remnant of the star will become
neutron star. Or if the mass of the core is heavy enough, it will finally becomea black hole.
Summarizing this session, only νe can be emitted in the very beginning of supernova explosion, and later after that neutrinos of all flavors will be produced.
1.4.2

Supernova Relic Neutrino

As 99% of the energy is released by neutrino flux in supernova process, it
is an important way to study supernova models by neutrinos. In Milky Way,
core collapse supernova is estimated to be 2 or 3 times per century [54], and SK
would observe ∼ 104 events at a typical distance of 10kpc. From a supernova
in the nearest galaxy Andromeda, SK would detect ∼1 event [29]. Although
there is no supernova neutrino detected since SN1987A, all of the supernova
exploded through the history of the universe have released neutrinos, and they
are refered as Supernova Relic Neutrino(SRN) or Diffuse Supernova Neutrino
Background(DSNB). The SRN spectrum is a convolution of a single supernova
spectrum, by considering the supernova rate as a function of z, and neutrinoenergy redshifted. The first models for SRN were founded even before SN1987A[55], while the studies about SRN largely increased after SN1987A. Figure 2
shows the ν̄e fluxes predicted by the following theoretical models:Constant SN rate model [56] (Totani et al. 1995)
The supernova rate was assumed as a constant in the whole universe in thismodel.
Population synthesis [57] (Totani et al. 1996)
It is founded by the same author by using a population synthesis of the stars in
the galaxy evolution.
Cosmic gas infall model [58] (Malaney et al. 1997)
This model uses the supernova rate which is obtained from the density distribution of interstellar gas as a function of red shift parameter.
Chemical evolution model [59] (Hartmann et al. 1997)
This model considered the history of the cosmic star formation which is obtained
from the chemical evolution.
Heavy metal model [60] (Kaplinghat et al. 2000)
In this model, the theoretical upper limit of SRN flux is predicted by using theabundance of heavy elements.
12

LMA neutrino oscillation [61] (Ando et al. 2002)
The parameter for neutrino oscillasion is considered in this model.

Figure 5: Colour lines are the ν̄e flux predicted by each model. Black lines
represent other background neutrinos.
The reason why only ν̄e is considered is that the cross section of ν̄e reaction
in water Cherenkov detector is two orders larger than the others in the mainenergy range of SRN. Even though ν̄e reaction has large cross section, SRN flux
is still hidden in various background and has not been detected yet. As canbeen seen from the Figure 2, Solar and reactor neutrinos are dominant below
10MeV, and atmospheric neutrinos become dominant above 30MeV.

1.5

Quasi-Elastic neutrino reaction and de-excitation gamma

In water Cherenkov Detector, there are two kind of processes which can
produce gamma. One is inelastic scattering ν +16 O → ν +16 O∗ , by which
the excited Oxygen nucleus emit a nucleon or gamma rays between 1∼10 MeV.
Another is quasielastic scattering, where nucleon is knock out from the Oxygen nucleus. Neutral Current Quasi-Elastic Scattering(NCQE) 24 and Charged
Current Quasi-Elastic Scattering(CCQE) 25 both have de-excitation process,but through CCQE a muon will be produced. Calculation result shows that QE
process overwhelms the inelastic process at neutrino energy of Ev >200MeV [15].
ν +16 O → ν + p +15 N ∗(24)
ν +16 O → ν + n +15 O∗

13


16
−
15
∗

 νe +16 O → e+ + p +15 O ∗νe + O → e + n + N
νµ +16 O → µ− + p +15 O∗ν µ +16 O → µ+ + n +15 N ∗(25)

The CCQE signal from atomospheric neutrino is a kind of unreducable background for SRN search. Refering to the recent SRN analysis in SK [62], the sample has been applied with the coincident cut with a previous gamma or muon toremove atomospheric ν̄µ CCQE events, but for those CCQE events with invisible
muon and no de-excitation gamma, they still remains as the largest backgroundin the final spectrum. Since the muon product can be invisible in SK, we hope
to know the branching ratio of de-excitation gamma to estimate the remaining
background component.
Oxygen has 8 protons and neutrons, 2 in 1p1/2 state, 4 in 1p3/2 state, and2 in 1s1/2 state. Gamma can be released by de-excitation when 1p3/2 or 1s1/2
nucleon is knock out, for 1p1/2 case there is no gamma because it is already theground state. Calculated result of the spectroscopic factors for each states issummarize in Table 1. Branching ratio for each state is from some indirect measurement, which uses electron or proton beam to make proton hole in Oxygen
nucleus. These measurement has a constrain of proton observation, but since
the energy level of proton and neutron is similar and 16 O is an isoscalr nucleus,
we can roughly consider neutron case as the same.
NCQE cross section for on Oxygen was also measured by Tokai-to-Kamioka (T2K)
long baseline oscillation experiment with a median neutrino energy of 630MeV,
while in the result, de-excitation gamma was observed in SK and the NCQE
and CC component was estimated from simulation [16].
Unlike NCQE, CCQE interaction gives a muon signal and a coincident searchfor decay electron can ensure the correctness of CCQE event selection. However,the difficulty is muon signal and de-excitation gamma appears from the same
position and they occur in a time difference less than 1ns. The analysis method
will be introduced in Section 8.

Spectroscopic Factor
Br(γ > 6M eV ) from p-holeBr(γ > 6M eV ) from n-hole

1p1/2
0.632
0%
0%

1p3/2
0.703
91.8%
86.9%

1s1/2
0.422
14.7%
14.7%

Table 1: Spectroscopic factors and branching ratios for 16 O de-excitation γ ray
above 6MeV. Spectroscopic factors are from calculated result [15], branchingratios are from measurement which make proton or neutron holes in 16 O byelectron or proton beam [17–19].

14

2Super-Kamiokande detector

2.1

Detector Structure

Super-Kamiokande (SK) is a water Cherenkov detector built for neutrino
observation [20]. The last three letters can be explained as two meanings: ”Neutrino Detection Experiment” and ”Nucleon Decay Experiment”. SK is locatedat latitude 36◦ 250 N and longitude 137◦ 180 E, Mt.Ikenoyama of Gifu Prefecture,
Japan. It is a water tank of cylinder shape, with 39.3m diameter and 41.4m
height, filled with 5kton ultra water. SK is built in 1000m underground, where
was once one of the world’s largest Zinc mines but now used for neutrino detection experiment. Since the underground environment, SK avoid suffering fromthe high rate background of down-going cosmic muons. SK observe an averagerate of cosmic muon by ∼2Hz, which is 10−5 comparing with the mountainsurface.
Additionally, SK has a structure of two layers to reject the muon backgrounds: An inner detector(ID) with 11139 (in SK-IV) inward-facing 20-inch
Photomultiplier Tubes (PMT), and an outer detector(OD) with 1885 outwardfacing 8-inch PMTs. The inner volume is 32.5kton, the PMT coverage is ∼ 40%(except SKII) and the rest ∼ 40% is covered by polyethylene black sheet to
reduce the unneccessary photon reflection. The outer detector part is 2m thickness, also filled with water, and can be used as veto signal for cosmic rays as wellas the gamma rays from the surrounding rock, because those signals will give
Cherenkov light in the outer detector before they enter the inner part. Also, thewall of the rock in the mine is covered by “Mineguard”, a polyurethane material which can prohibits the radon from the rock. The dome above the tank isused for 5 electronics huts, a linear accelerator used for detector calibration, andstorage for various equipment. The data acquisition system will be explained in
Section2.3, while the calibration methods will be in Section 5.

Figure 6: The detector appearance. The bottom right shows the location within
the mountain. (cutaway view)15

2.1.1

History and current status

SK has five stages(including the latest one), called SKI, SKII, SKIII, SKIV
and SK-Gd. SK start running from April 1996 with 40% ID cathode coverageand 11146 ID PMTs, and stopped in July 2001 for routine maintanance andupgrades. The period from April 1996 to July 2001 is referred as SK-I period. During this period, in 1998, the SK collaboration released the first highly
significant experimental evidence for neutrino oscillations as well as the solar
deficit.
On November 12, 2001, after the maintenance was completed and the tank
was being refilled with water, a catastrophic accident occurred. One of the
PMTs imploded when the tank was ∼75% full, causing a chain reaction resulting in the implosion of almost all the underwater PMTs. The chain reaction
was violent enough that it registered as a seismological event to a seismologicalrecorder 8.8 km distant. 6665 out of 11146 ID tubes were found to be destroyed.
The water was drained again, and the broken glass cleaned. Unfortunately remaking so many PMTs is a time consuming process, and it was decided to run
the detector in the meantime. The surviving tubes were redistributed throughout the detector, after being fitted with new acrylic coverings (tube front) and‘FRPs’ (fiberglass reinforced plastic, tube back), to prevent the shock wave ofany future tube implosions from causing another chain reaction. The detector
was turned back on with 19% cathode coverage (47% that of SK-I), and ran
until the new PMTs were ready in July 2005. This period of data is referred toas SK-II.
The refitting of the PMTs took almost a year, and the detector was operational again in June 2006. In August 2008, the SK electronics system wascompletely overhauled and a new upgraded system was put in place. The data
period from June 2006 to August 2008 is called SK-III. Since September 2008,
after the electronics upgrade was completed, data taking has continued, andthis period is referred as SK-IV.
During SK-IV, a project proposed to disolve Gadolinum into SK tank has
been planned [29]. Gadolinum can enable the detector to distingush neutrino
and anti-neutrino, by making the neutron signal produced in anti-neutrino interaction to be visible. Various works have been carried out, such as Gadolinum
quality check, impurity search, Gadolinum circulation system design and so on.
In order to test Gadolinum performance, a 200ton water Cherenkov detector
called Evaluating Gadolinium’s Action on Detector System (EGADS) was built
inside the same mine with SK. Since the success of EGADS experiment and theprovement of feasibility of Gadolinum project, SK was shut down in May 2018
for Gadolinum disolving work. With Gadolinum-load detector, SK started a
new stage from October 2018, called SK-Gd.
Overall, SK-IV is the most stable and long-lasting period of SK, and in thisthesis, most of the analysis will be based on the data set of SK-IV.

2.2

Dection method

2.2.1

Cherenkov radiation

SK is a water Cherenkov detector, it detects the neutrino by observing
Cherenkov light. The speed of light in medium is c/n, when the index of refraction is n. If a charged particle passes through the medium with velocity
16

PhaseStart
End
ID PMT
OD PMT
PMT CoverageElectronics
Energy Thre.

SK-I
Apr 1996
Jul 2001
11146
1885
40%
ATM
4.5MeVSK-II
Oct 2002
Oct 2005
5182
1885
19%
ATM
6.5MeV

SK-IIIJul 2006
Aug 2008
11129
1885
40%
ATM
4.0MeV

SK-IV
Sep 2008
Jun 2018
11129
1885
40%
QBEE
3.5MeV

SK-Gd
Nov 2018
Running
11129
1885
40%
QBEE
-

Table 2: The summary of characteristics of each SK phase. Energy thresholdis the recoil electron kinetic energy, adopted for the analysis. The threshold of
SK-Gd is not shown because currenly when this thesis is being written, SK-Gd
is still in test.

larger than c/n, the particle will emit a forward cone of light. This is named
as Cherenkov radiation, while P.A.Cherenkov is the person who discovered thiseffect in 1934. To theoretically explain it, when a charged particle travels in themedium, it will emit spherical wave by electric polarization. If the velocity of theparticle exceeds c/n, the spherical waves can interfere with each other and form
radiation. Otherwise, they will never meet each other and stay independent.

Figure 7: when velocity of thecharged particle is smaller thanc/nFigure 8: when velocity of thecharged particle is larger thanc/nThus, the angle θ between Cherenkov light and the direction of the charged
particle can be calculated with the ratio of light velocity in the medium c/n andthe velocity of incident particle as follows:cos θ =

1
nβ

(26)

Here β = v/c. Considering the critical condition of Cherenkov radiation is
β > 1/n ,the energy threshold for incident particle to emit Cherenkov light is:nm
Ethr = √(27)
n2 − 1
17

Particlee±µ±π±Threshold(MeV)
0.767
157.4
207.9

Table 3: Cherenkov threshold for charged particles in water.
For water Cherenkov detector, the index of refraction is 1.33, so the thresholdof total energy for the main particles is shown below.
In most cases, the signals observed in Super-K are from electrons or muons
with β ≈ 1, so the Cherenkov angle is close to 42◦ . The number of photon (N)produced by Cherenkov radiation in unit length dL can be calculated as:1
d2 N
2παZ 2(1 − 2 2 )
=
dλdL
λ2
n β

(28)

Here λ is the wavelength of Cerenkov radiation, Z is the charge of incident
particle, and α is fine-structure constant. By integrating it, the number ofphoton in specific wavelength range (λ1 , λ2 ) is:dN
2παZ 2 1
11=( − )(1 − 2 2 )
dL
n
λ1
λ2
n β

(29)

For example, when an electron travels 1cm in Super-K tank with β ≈ 1,
about 340 photons would be produced in Super-K PMT sensitive range of
280nm∼660nm. It can hardly be observed and that is why we need photonmultiplier tubes (PMT).
To consider muons, they can be more complicated because of the continuous
emission of light as the muon is usually energitic to travel through the tank.
Some muons are created by neutrino interactions, but most of the muons arecosmic ray muons (∼2 Hz) which is able to penetrate the rock overburden andenter the detector. These events can be easily distinguished by use of the OD,
and the entire muon track can usually be reconstructed. Other charged particles,
such as pions, are also seen.
SK is capable to observe neutrinos from the sun, from atmospheric cosmic
ray interactions, from supernovae, from neutrino beams, and most likely from
nuclear reactors. Due to the large volume of SK, the detector is also useful for
proton decay and dark matter searches.
2.2.2

Photomultiplier tube

The PMTs being used in Super-K are produced by Hamamatsu Photonics
c.o, and the development work was done under the cooperation of Super-K
group. The inner detector PMTs are 20-inch Hamamatsu R3600 PMTs, andthe outer detector PMTs are 8-inch Hamamatsu R1408 PMTs [21]. The designof inner detector PMTs are based on the PMTs being used in KamiokaNDE
experiment, with several improvements.
Figure 10 shows the structure of a PMT. The photocathode material isbialkali(Sb-K-Cs) which has a sensitive region of 300-600nm and the quantum
efficiency peaks at 360-400nm as shown in Figure 12. Electron photons could
be multiplied with a gain of ∼10 7 by 11 chain dynodes. A good time resolution
18

Photoelectronic Surface
Material
DynodeQuantum Efficiency
Dark Rate
Time Resolution

50cm
Sb-K-Cs
Venetian Blind type
22% at 390nm
3.5kHz3ns RMS at 1pe

Table 4: Properties of Super-K inner PMT
is critical to have better vertex resolution for low energy events. For single
photoelectron the transit time spread is about 2.2 ns. The desired 1 p.e. peakcan be clearly seen in Figure 13, where the peak close to zero ADC count is dueto the dark current. Super-K PMTs not only have better time resolution than
Kamiokannde PMTs, but also become able to identify single photon signals. It
is an improvement of great significance because most of the signals in low energy
events like supernova neutrinos are due to single photon.

Figure 9: Photo of a bare PMT

19

Figure 10: Inner structure of PMT

Figure 11: The acrylic cover and FRP case used from 2001 to protect the
PMT from shock wave [20].
The threshold is set to 1/4 pe and the rate of dark noise above this thresholdwas ∼3.5 kHz when in Super-K first phase. The magnetic field of the earth canlargely affect the photoelectrons so Helmholtz coils are used for compensation.
The residual geo-magnetic field is kept less than 100 mG in every position of
the tank.
In the accident happened in 2001, about half of the PMTs were broken as aresult of shock wave. Since that, PMTs in Super-K are all installed with acrylic
covers and FRP cases as shown in Figure 11. The transmittance of the acrylic
cover is about 96% and it doesn’t affect the event reconstruction.

20

Figure 12: The quantum efficiency of 50cm PMT as a function of wavelenth.

Figure 13: Single photoelectron pulse height distribution.

21

Figure 14: Neutrino interaction cross sections in a large water Cherenkov detector, including resolution and threshold effects [27].
2.2.3

Neutrino interaction in water Cerenkov detector

The neutrinos react with electrons or Oxygen in water and finally charged
particles are observed by Cherenkov radiation. The main reactions of neutrinos
in water Cerenkov detector is shown below.
Inverse Beta Decay:ν e + p → n + e+
(30)
This is the most dominant interaction and it has the largest cross section in theenergy region of Supernova neutrinos. It occurs between the anti-neutrinos ofelectron type and free protons in water. Most of supernova neutrinos will be
detected via this interaction. However, the direction information will be lost in
this interaction.
Neutral Current Interaction:
νx + 16 O → νx + γ + X(31)

ν + e− → ν + e−

(32)

Elastic Scattering:This is the elastic scattering between neutrinos and electrons in water. It
happens in all flavours, but the electron-neutrino has about 6 times larger crosssection than others. This is main reaction for solar neutrino detection, as thesolar neutrino only has νe . For supernova neutrino, the cross section of elastic
scattering is two order smaller than inverse beta decay but it can conserve thedirection information of the incident neutrinos, which is very important for22the early announcement to the other measurement such as optical observation,
gamma-ray and so on.

2.3

Data acquisition system

2.3.1

Electronics

There are 5 huts on the top of the SK tank, named with a ”Central hut”
and four ”Electronics huts”. In the electronics huts, HV power supplies and theelectronics system are placed; the former supplies high voltage to ID and OD
PMTs and the latter digitizes the analog signal from the PMTs. In the central
hut, the trigger system and the control electronics are placed and the signals
digitalized in the four Electronics huts are merged. Electronics of SK changed
from SKIV. The difference is front-end part. From SKI to SKIII, the frontend electronics is called the Analog Timing Module (ATM) [23, 24], which isbased on TKO (Tristan KEK Online) standards [25]. From SKIV, QTC-BasedElectronics with Ethernet (QBEE) [25] is used. The most difference between
ATM and QBEE is, in ATM period data is recorded with hardware trigger,
while QBEE collect all the data and apply software trigger later.

Figure 15: The schematic view of the DAQ system used in SK-IV [28].

2.3.2

Trigger

The detector has triggers for ID and OD. ID triggers includes 4 types: Super
Low energy (SLE) trigger, Low energy (LE) trigger, High Energy (HE), Super
High Energy (SHE). The four types of ID triggers are seperated by different
HITSUM threshold.
The HE trigger is mostly used for atmospheric neutrino studies. The LE
trigger is good for SN relic events and solar neutrino events. The SLE trigger was
introduced to look for low energy solar neutrino events. The lower threshold of
the SLE trigger causes a much higher trigger rate (roughly an order of magnitude
more triggers per MeV the threshold is lowered). Much of the increase comesfrom background, such as gammas from the rock walls and radioactivity from23

Figure 16: Left: The picture of the new front end electronics QBEE. Right:
The block diagram of the signal processing in the QBEE.
the PMT glass itself. In order to effectively lower the low energy threshold andnot be swamped by background, a prefilter was introduced for SLE events called
the intelligent trigger (IT).
The IT does a quick vertex fit to SLE events and only saves those within
the fiducial volume (inner 22.5 ktons), eliminating much of the background. As
computing power increased throughout SK-I, the IT system was upgraded tohandle a higher event rate, and the SLE threshold was lowered multiple times
as SK-I progressed. During SK-II, all thresholds were lowered to reflect thedecreased cathode coverage. Trigger thresholds are summarized in Table 2.1.
The OD trigger has remained constant at 19 OD hits.

Threshold (hit)Gate width (µs)

SLE
31-0.5∼1.0

LE
34
-5∼40

HE
45
-5∼40

SLE
51
-5∼40

Table 5: The threshold for each trigger and its event time width.

24

3Event Reconstruction

3.1

low energy event

3.1.1

Vertex reconstruction

In SK analysis, usually the events below 100MeV are considered as low energy events, and those above 100MeV are considered as high energy events.
For low energy events, the expected signals are the Cherenkov light of the electrons and positrons around 10MeV. Those electrons or positrons are most fromneutrino-electron elastic scattering or inverse beta decay with proton. Since
the energy loss of electron in water is ∼2MeV/cm, the travel distance before
stopping in the tank is usually very short (∼10cm). As a result, for low energy
event, the Cherenov light source can be considered as a point. Another fact is,in the low energy as ∼10MeV, most of the hit PMTs can only get one photon,
thus the charge response of the PMTs become insignificant.

Figure 17: The timing distribution of hit-PMTs in one event. Horizontal axis
is the relative timing of the each hit-PMT obtained from data.
There are several vertex reconstruction fitter in SK, Kai fit, Clus fit, andcurrently the one with best vertex resolution is called BONSAI (Branch Optimization Navigating Successive Annealing Interactions) fitter [35].
BONSAI search for the vertex by 2 steps :1. Many possible starting points are constructed out of combinations of four
hit PMTs. Then substract the time of flight for the photon from the test
point ~x to each hit PMTs as t − ttof − t0 , while t is the original time
response from PMT, ttof is the time of flight, and t0 is the time offset.
Figure 17 shows the time response of all the hit PMT in one typical event.
2. Calculate the likelihood for the test point ~x. The likelihood is defined as:
L(~x, t0 ) =

N
hit
X

log(P (t − ttof − t0 ))

(33)

i

P (t − ttof − t0 ) is the probability density function (PDF) of the timing
residual for a single photoelectron signal as shown in Figure 18. The PDF
is made from LINAC calibration data (will be explained later).
25

3. Find the best point by maximizing the likelihood after checking against
12 nearby points on a grid of successively smaller sizes, which is down to 1
cm. The dark noise background is considered to be flat in time axis. The
fitter also returns a fit goodness value which can be used to evaluate thequality of the fit and identify misfits.
(i)−t0 2t(i)−t0 2
Pallhit −( tt−tof√√) −( t−tof
)
2ω
2δ
e
e
gv =
tt−tof (i)−t0
Pallhit −(
√
)2
2ω
e

(34)

Figure 18: Timing PDF used for vertex reconstruction in BONSAI fitter. The
PDF is made from LINAC calibration data. The peaks around 40nsec and110nsec are caused by the after pulses.

3.1.2

Direction reconstruction

Sine the Cherenkov light makes a ring-like pattern of photons, it is possible
to reconstruct the direction of the event. The direction reconstruction uses a
maximum likelihood function which is defined as:
L(~a) =

N20
X

log(f (cosθi , E)) ×

i

cosΘi
a(Θi )(35)

1. N20 : The number of hit PMTs inside 20nsec window near t − ttof − t0 = 0
2. f (cosθi , E): θi is the angle between the particle direction and the vector
from the vertex to i-th hit PMT. f (cosθi , E) is the expected PDF of θi ,which depends on the energy because the particle with different energy will
suffer from multiple Coulomb scattering at different degrees. The energy
dependence is evaluated by SK simulation of mono-energetic electrons.

26

i
3. cosΘa(Θi ) : Θi here means the angle between the verticle direction from i-th
PMT surface and the vector from reconstructed vertex to i-th PMT, whichmeans Θi is the expected incident angle. a(Θi ) is the photon acceptanceof i-th PMT, depending on the incident angle.

The directions are scanned by grid search, each direction is estimated andfinally the best one will be selected as reconstructed direction. To explain more
about the estimation method: firstly, the vertex and tentative direction are
already known, a hit PMT(either is fine) is choosed as the starting point andthe vector from vertex to starting point is used as a reference; secondly, for theother hits, the azimuth angle to the reference vector is caculated and then theazimuth angle distribution is made(example event shown in Fig 19); thirdly,
the uniformity of the azimuth angle distribution is defined by the arrow length
shown in Fig 19. If the event is normal, the distribution is uniform and thearrow length is small. On the other hand if the hits intend to be a cluster, thenthe arrow length is large.
After the recontructed direction is fixed, the recontruction precision is calculated by Equation 36, which means a Kolmogorov–Smirnov test between thereal event and a ideal event with the same number of hits uniformly distributedon the Cherenkov cone.
max[θunif orm (i) − θdata (i)] − min[θunif orm (i) − θdata (i)](36)
2πFrom the defination we could know that when gd = 0.0 it means the best fit,
and when gd = 1.0 it means the most poor fit. gd will be recorded as well as thereconstructed direction, to be used for background removing in later process.
gd =

3.2

Mu event

Muons in the detector were fitted by software tool called Muboy. Muboy isdesigned for a passing-through muon that it is searching for the entry positionand exit position of the muon. However, it not only fits simple single throughgoing muons with a single entry and exit position, but actually looks for different
types of muon events as it fits, and categorizes the event accordingly. Muboy
categorize each event as being one of the following:
1. Single through-going(∼ 82%). This is the most common case. The muon
simply entry the tank, and then goes out.
2. Stopping muons(∼ 7%). The muon stops inside the tank, which means it
has an entry position but no exit position. Usually the muon will decay
to electron after it stops.
3. Multiple tracks(∼ 7%). More than one track are found.
4. Corner clipper(∼ 7%). The muon with very short track length(L < 7m).
Usually it refers to the muons entry from the top and then exit from thebarrel.
Muboy is not the only muon fitter in SK, but it is capable to categorize
muon into different types and fit the multiple tracks. A fitting goodness is also27

Figure 19: Left figure is shown as normal event and the right one is the eventwith small cluster. Lower distribution fill the azimuth angle of each hit-PMT
from the start of one point [31].

28

given, while in case of multiple track, each track will have a fitting goodness.
SK has other fitters which are also able to fit muon. One is APfit, currently theofficial reconstruction tool for high energy analysis. Another is fitQun, a newly
developped recontruction tool with better precision, but have not been used
as official yet. However, since this thesis is mainly about low energy analysis,
the details about APfit or fitQun, as well as the fitting algorithm, will not be
discussed.

29

4

Simulation

4.1

SKdetSIM overview

SKdetsim is the Monte Carlo(MC) simulation software of Super-K detector.
It’s based on Fortran language and GEANT3.21 simulation framework [39],which is developed by CERN, widely used in high energy experiments, andqualified by many physicists. Most of the physics process in SKdetsim follow
GEANT3, except that the parts of Cherenkov photon generation and photon
propagation in water are customized by SK group. This is because when SKdetsim was been developed, the GEANT3 version at that time didn’t have thepackage for Cherenkov process.
The simulation process can be seperated into four steps to help to understand:
1. Particle generation.
2. Particle tracking and Cherenkov photon emission
3. Progration of Cherenkov photon
4. PMT response and electronics simulation
The four steps will be explained in the following subsections.

4.2

Particle generation

SKdetsim is used for different purposes: solar neutrino, supernova relic neutrino, high energy analysis, and calibration analysis. It should be mentionedthat SKdetsim does not treat the neutrino reaction itself. For solar neutrino,
the events are made by 8 B and hep flux at earth, then the cross section of neutrino reaction is taken into calculation and finally the recoil electrons are given
into SKdetsim to simulate the signal in the tank. Simulation for relic neutrino
does the similar thing, but the Ando LMA model and HBD 6 MeV model areused for neutrino source, and Strumia and Vissani’s calculation for cross sectionof Inverse Beta Decay is used. For the simulation of LINAC calibration, the pipe
is defined in SKdetsim program and the electron is generated from the pipe atfixed energy. For the simulation of Ni calibration, the Ni ball is defined but the
gamma is generated by a pre-loaded spectrum, not from the 252 Cf source andNi-neutron reaction as the real case.
However, this thesis does not focus on solar neutrino spectrum or relic neutrino spectrum, though the data set of those two are used later. So the detailof neutrino spectrum generation and cross section of neutrino reaction will not
be discussed here.

4.3

Particle tracking and Cherenkov photon emission

When the particle is generated into tank water, SKdetsim track the particle
by each small step. The dominant processes for low energy region(a few MeV
to a few tens MeV) are multiple scattering, ionization loss, δ-ray production,bremsstrahlung and annihilation of positrons for electrons, pair creation, Compton scattering and photoelectric effect for gammas. In each step, SKdetsim calculate the flying distance, energy loss, direction change, and the physics process30happened in this step. If there is any secondary particle(including Cherenkov
photon) produced in this step, the information of the secondary particles will
be saved into the buffer memory and be tracked later. Usually when a charged
particle goes below the Cherenkov threshold, SKdetsim will stop the tracking.
In particle tracking, the other physics processes follow GEANT3 and detail
information can be found in GENAT3 manual. Special attention should be
paied to Cherenkov photon production, as this part is the exception that it isdeveloped by SK group. For each tracking step, the Cherenkov photon number
N is give by equation 37.
d2 N =

12πα(1 − 2 2 )dxdλ
2
n(λ)λ
n β

(37)

Here the n is the refractive index of water(n ≈ 1.33), α is the fine structure
constant, and β is the velocity of the electron in unit of the light velocity invacuum. In per track length dx, dN is calculated on per wava length dλ. It
should be mentioned that only the range of 300-700nm is intergrated because
PMT is only sensitive to this region. The number of Cherenkov photon produced
is approximately proportional to the energy of electron, this property is used
to determine the energy scale as explained in previous section. The direction of
Cherenkov photon is given by equation 38.
cosθ =

1nβ(38)

θ here is the open angle.

4.4

Progration of Cherenkov photon

The velocity of Cherenkov photons depends on its wavelength. The groupvelocity can be evaluated as:vg =1n(λ − λ dn(λ)
dn )

(39)

where c is light velocity in vacuum, λ is light wavelength, was used for thelight velocity in the water. The denominator of this equation can be treated as”effective index”.
Though the refractive index of water is ∼1.33, actually it depends on the
wave length as:s
a1n(λ) =+ a2 + a3λ2 + a4λ3 + a5λ6(40)
λ2 − λ2a
The parameters are summarized below, and they are obtained from realmeasurements.
When the photon is reflected or scattered, the direction vector changed andthe tracking will continue. While the photon is absorbed in water, it will bekilled and the tracking will be stopped. The intensity of the light traveling inwater exponentially decreases as:
I(x) = I0 (λ)exp(

31

−x
)
L(λ)

(41)

λ2a
a1
a2
a3
a4
a5

0.018085µm
5.7473534 × 10−3 µm
1.769238µm−2.797222µm
8.715348µm−1.413942 × 10−3 µm
Table 6: 111

Here x is the traveling length, L(λ) is the total attenuation length, I0 (λ)
is the initial light intensity. The total traveling consist of absorption, Rayleighscattering and Mie scattering. The three coefficients in total attenuation can beunderstood as:
Lattn. =1
αabs (λ) + αRayleigh (λ) + αM ie (λ)(42)

From the measurement result, it is difficult to know the Mie scattering partfrom the total scattering. However one feature we know is that Rayleigh scattering is symmetric in forward and backward, while Mie scattering is asymmetricand it scatter more light to the forward direction. So the equation above canbe modified to:Lattn. =1
αabs (λ) + αsym (λ) + αasym (λ)(43)αsym and αasym are the sysmmetric and asysmmetric component in scattering respectively. αsym consist of Rayleigh and symmetric Mie scattering, while
αasym consist of forward Mie scattering.
In SKdetsim, they are calculated as empirical function as:αabs (λ) = P0 × Pλ41 + C(44)
αsym (λ) = Pλ44 × (1.0 + Pλ25 )
αasym (λ) = P6 × (1.0 + Pλ47 × (λ − P8 )2 )P0 ∼ P8 are fitting parameters tuned by calibration data, the parameter C
in αabs is the amplitude based on the experimental data.

4.5

PMT response and electronics simulation

When the photon arrive the PMT surface, it can be reflected or absorbed.
When it is absorbed, then it will follow the possibility of quantum efficiency to
be randomly transfered into a signal or not.
The probability of the PMTs to produce the one photoelectron when thegenerated photon reaches its surface is defined as:
P rob(λ, i) = QE(λ) ×A(λ, θi )× COREP M T × qetable(i)A(λ, 0)The parameters in this equation mean:

32

(45)

1. QE(λ): the quantum efficiency of 20-inch PMT depending on wavelengthof incident photon.
2. A(λ, θi ): the acceptance of photon on PMT surface, which depends on the
wavelength and incident angle. In other words, it represent the possibility
for a photon to be absorbed by PMT cathode and no to be reflected ortransmitted.
3. COREPMT: a global parameter used for all PMT to correct the average collection efficiency. This is a tunnable parameter and adjusted with
LINAC calibration data.
4. qetable(i): the relative quantum efficiency for i-th PMT. This is measured
by Ni calibration.
After a photon is successfully transfered into a photoelectron, the output
charge follows 1p.e. charge distribution which is made from Ni calibration
data. If the output charge exceeds the electronic threshold(∼ 0.2p.e.), it will berecorded into the data. Finally the time and charge information of all the hit
PMTs will be applied with software trigger, which is the same one used for real
data.

Figure 20: The solid line is refractive index depends on wavelength. The dashedline is the ”effective index”.

33Figure 21: The possibility of absorption, reflection and transmission dependingon the incident angle by the case of Unpolarized wave of 345nm.

Figure 22: An incident photon hit on PMT surface, it can be reflected or transmitted, or absorbed to be changed into a photonelectron.

34

Figure 23: The component of Mie scattering, Rayleigh scattering and absorption
in attenuation length in SKdetsim.

4.6

Summary of SKdetsim

4.6.1

Tunable Parameters

Overall, the tunnable parameter in SKdetsim can be summarized below:
Water parameters
The parameters explained in previous section are measured by laser calibration
which inject the laser light from tank top, barrel and bottom. There is usuallyno need to modify these parameters unless the measurement has been updated.
Top-Bottom asymmery
This is not mentioned in previous section. From Ni calibration data it canbe found that the hit rate for PMTs in top and bottom has small difference.
To reproduce this effect in simulation, SKdetsim has a parameter to correct
the top-bottom asymmetric(TBA table). The TBA table is calculated from Nidata, and double-checked by auto-laser calibration measurement.
PMT time resolution and after pulse
Both of the two are tuned in SKdetsim by LINAC calibration data to make
PMT time response to be consistent with data. After pulse is the later peakwhich appear after the main peak in PMT timing distribution. The reason for
after pulse is considered as a photoelectron which is back scattered at the firstdynode, loses its velocity against the electric field and then comes back to thefirst dynode to produce a decay hit. The timing distribution tuned by LINAC
can be double-checked by Ni calibration data.

35

COREPMT(global parameter of collection efficiency)
As explained in previous section, COREPMT is a global parameter represent
for quantum efficiency of all PMTs. The small difference between PMTs iscorrected by the parameter of relative quantum effciency. COREPMT is tuned
by LINAC, and relative QE table is made by counting PMT hit rate with Ni
calibration data
4.6.2

Some other explanation

For the reader who is interested to make further modification in SKdetsim,
some special treatment in SKdetsim program is explained here:
1. SKdetsim is based on Fortran Language, and zebra stream(an old data format with ”.zbs” as the suffix of output data file) is used as buffer memory.
When a new particle is generated, the information including the particle
ID, position, direction, momentum, will be stored in zbs buffer. In tracking step, SKdetsim will read the particle information from zbs buffer andmake the tracking one by one. If any secondary particle is produced intracking, the information will also be stored into zbs buffer, and the tracking of secondary particle will start after the tracking of parrent particle
is finished. The only exception is Cherenkov photon. When a Cherenkov
photon is made, the tracking of parrent particle will be suspended andthe new Cherenkov photon will be tracked immediately. This is because
Cherenkov photon can be largely produced and if they are tracked later,
the memory is not enough.
2. To save running time, the number of Cherenkov photon is half reduced,
and as compensation, the absolute QE, which affected by COREPMT, isdoubled from the real number. Attention should be paid to this if someone
try to define a High QE PMT in SKdetsim, because the absolute QE isdoubled and can exceed 100% to cause bug.
3. When the photon is reflected on PMT surface or black sheet, the photon
is actually killed and a new photon with the reflected direction is store
into the zbs buffer. Differs from GEANT3, GEANT4 does not kill thephoton when reflected.

36

5Calibration

5.1

Xe and Ni calibration

Xe and Ni calibration are used to calibrate PMT parameters [30]. In PMT
calibration, Quantum Efficiency(QE) and gain are two important factors. Xecalibration is to determine the High Voltage(HV), and Ni caliration is to obtain
the absolute gain of each PMT.

Figure 24: Setup of Xe scintillation light. The optical fibers are connectted to
APD and a 20-inch PMT, both for monitoring the light indensity [30].
QE also includes the collection efficiency of photoelectrons onto the firstdynode of the PMT. Gain is the conversion factor from the number of collected
photoelectrons to output charge (in units of pC). Low energy events like relicneutrino or solar neutrino usually consist of single photoelectron (pe) hits sotheir analysis heavily depends on QE calibration, while high energy events depend more on gain calibration. The output of the PMT depends on the voltage
applied on it. For further explaination, the charge response of one PMT can be
defined as:
Qobs (i) ∝ Nphoton (i) × QE(i) × A(i)(46)

Qobs is the charge response finally observed in i-th PMT. Nphoton (i) is thenumber of photons hit on the photocathode of i-th PMT. QE(i) and A(i) are
the QE and gain of i-th PMT respectively.
So first of all, it is necessary adjust the HV of the PMTs to make sure the
PMTs with same geometry will give the same Qobs . For this calibration, Xe flash
lamp is used as a light source. Light from the lamp passes through a UV filter
and then be injected into the tank through optical fiber, finally arrive at a scintillator ball placed near the center of the tank((x, y, z) = (353.5, 70.7, 0.0)cm).
The other two fibers go to an avalanche photodiode(APD) module and a darkbox installed with a 2 inch PMT, which are both used to monitor the lightintensity of Xe lamp. The scintillator ball is an acrylic ball containing 15ppm
POPOP as a wavelength shifter and 2000ppm MgO as a diffuser to make thelight emission from the ball as uniform as possible. The light emitted from thescintillator ball has the wavelength near to Cherenkov light(peak at 440nm).
The PMTs with the same geometry are divided into one group. In eachgroup, there are standard PMTs which have been pre-calibrated to be used as37

Figure 25: left shows the locations of ”standard” PMTs which are indicated by
the red points. Right shows the PMT grouping. PMTs with the same geometry
are divided into the same group [30].
references for the other ones. After HV calibration, finally the relative difference
for the PMT charge response is ∼ 2%.
After HV calibration, Ni-Cf source is used to measure the absolute gain. This
source is made by 58 Ni ball and a 252 Cf stick inserted inside the ball. Neutrons
are provided by 252 Cf, and 9MeV γ rays are emitted isotropically by neutron
capture on 58 Ni. The source is made like a ball as shown in the Figure 26.
Since more than 99% of observed signals are due to single photoelectron,
the result of Ni calibration can be used to determine the factor to convert the
charge(pC) to the number of photonelectron. At the begining of SKIV, thefactor is determined as 2.645pC/p.e..

5.2

LINAC Calibration

The full name of LINAC device is Mitsubishi ML15MIII electron linear accelerator (LINAC) [37], which is used for calibration of the absolute energy. It
was originally used to be for medical purpose and was acquired from a hospitalthen converted into a calibration tool for SK in 1996. It is permanently housedin the dome area above the tank, consisting of a special electron gun, steering
magnets and collimators.
LINAC is purposed to produce single electron event, so the steering magnets
and collimators are used to control the mono-energetic electron beam. Beamintensity is tunable, and the beam can be collimated such that electrons are
injected into the tank at a low enough rate(∼0.1 per bunch). Before injecting thebeam into the tank, a germanium detector with an hight energy resolution((1.92
keV at 1.33 MeV electron)) is set to measure the beam energy. The Ge detector
itself was calibrated by 0.662 MeV monochromatic gamma rays from 137 Cs and9.0MeV gamma from N i(n, γ)N ∗ . Multiple calibration port holes exist on the
top of the detector, and the beampipe can be lowered to a variable depth in thedetector. In this fashion, many different areas in the detector are accessible andcan be calibrated.
38

Figure 26: The Ni ball for absolute gain measurement. The 252 Cf source isinserted as a stick inside the ball [30].

Figure 27: Distribution of 1p.e. charge response in SKIII. The right one is inlog scale and left one is in linear scale. Due to the hardware threshold, it is
impossible to make measurement to 0pC, so the red line in left figure is a linear
fitting.

39

LINAC calibrations were performed in 2009, 2010, 2012, 2017 and 2018
during SK-IV period. Since the calibration hole positions are fixed and thelength of the beam pipe is also fixed, the calibration data were taken at the
fixed 6 points as shown in Fig 28.
After the LINAC calibration data is taken, the SK simulation is tuned tomatch the detector response. Comparing the peak positions in the Neff distributions, COREPMT is determined. After this value is fixed, the difference of
Neff between the calibration data and the MC simulation was evaluated andtreated as systematic uncertainty.

Figure 28: The fixed positions for LINAC data taking.

5.3

DT Calibration

A deuterium-tritium (DT) generator is used to check the absolute energy
scale for low energy events [36]. It operates in conjunction with the LINAC,
providing a cross check on the energy calibration. It is less time consuming
to use then the LINAC and better for monitoring long term stability with frequent checks. Furthermore, the device output is isotropic, and isn’t limited toonly downward going events like the LINAC, thus allowing study of direction
dependent effects.
40

The generator emits neutrons with the following reaction:3

He +2 H →4 He + n(47)

The energy of the generated neutron is 14.2 MeV and this energy is large
enough to create 16 N in water as
16

O + n →16 H + p

(48)

16

The N decays via several channels with its half-life of 7.13 seconds. The
main decay channels produce a 6.1 MeV γ ray and a 4.3MeV β(66%) and a
10.41MeV β(28%).
16
N →16 O + e− + νe
(49)
Although the beta spectrum is spread and it is not good for absolute energy
calibration, there are several advantages of the DT calibration over the LINAC
calibration. This provides the directional and the positional dependence ofenergy scale. The directional and positional dependences are important for thesolar neutrino analysis, and this is also treated as the systematic uncertainty insolar analysis.
DT data is taken every few months, and the energy scale as determined by
the DT agrees with the results of the LINAC calibration to better than 1%.
Various positions in the tank are monitored to assure the energy scale at all
regions are consistent.

Figure 29: The DT source is like a rocket, and it is put into the tank by alifter [36].

5.4

Water Transperency

5.4.1

By Laser

Lasers are used to directly measure light scattering and absorption parameters in a wavelengths of 337nm,375nm, 405nm, 445nm and 473nm. The three41

parameters αabs , αsym , αasym in SK simulation are determined by this measurement. The light is injected every one minute and the hit rate in the different
region of the SK ID are monitored. The measurement methods are as following:
1. The detector is divided into 6 regions, top and barrel (B1 to B5 as shownin Fig 30).
2. The scattered hit rate Fig 31 and the observed total charge Qtotal aremeasured, then the ratio of hits and Qtotal (hits/Qtotal ) distributions are
obtained for each divided region. The observed charges in the bottom
PMT are used for the reference to monitor the intensity of the laser light.
3. Several sets of MC simulations are generated with varying the coefficientsand the same (hits/Qtotal ) distributions are prepared.
4. The distributions of the calibration data and the MC simulation are compared using the χ2 defined asχ2 =X (Data − M C)2
2
2
δdata
+ δmc
region(50)

Data (MC) is the peak position of (hits/Qtotal ) distribution, and δdata(δmc ) is the standard variation of that.
5. Scanning the water transparency coefficients to minimize the value of χ2 .

Figure 30: The schematic view of the laser injection system for water parametermeasurement [30].

42Figure 31: The typical hit rate distribution for each layer. The horizontal axis
is the timing after subtracting TOF. The vertical axis is the hit rate. Red line
is simulation with best tune, black dot is the real measurement [31].

435.4.2

By decay electron

Cosmic ray muons enter the detector at a rate of about 2 Hz, and about
one in 20 will stop in the detector. These stopping muons decay into electrons
with the well known Michel spectrum. These decay electrons can be identified
by associating the event with a muon that closely near it in time, and whose
reconstructed muon stopping point is near to the location of the decay electron
vertex. With 10 years of data in SK-IV, a very large sample of the events has
been accumulated, which can be used as an independent means of monitoring
detector stability, as well as the water transarency change. The selection criteriaof decay electron events are following:
1. The time difference between the parent muon event and the decay electron
candidate event ∆T , selected as 3.0µsec < ∆T < 8.0µsec.
2. The reconstructed vertex of the decay electron candidate event is within
the 22.5kton fiducial volume.
3. The distance between the stopping point of cosmic ray muon and decay
electron candidate event is within 250cm.

Figure 32: A typical histgram of ln(Qobs ) vs distance r. The dashed vertical
black line corresponds to 1200 cm and a linear fitting can be made above that.
The blue line is fitted with a fixed intercept of 1.486 in y axis, red line is fittedwith intercept [32].
After selecting the decay electron events, the observed charge Qobs (i) of thei-th PMT and the distance(r) between the i-th PMT and the decay electron
candidate are calculated. Then, mean of the observed charge is evaluated as
P
Qobs (i)Qobs =(51)
Ntotal
44

The typical histograms of ln(Qobs ) and the mean distance r are shown in Fig 32.
To obtain the water transparency, the histograms are fitted by a linear fitting
function from 1200cm to 3500cm. The inverse of the slope of the linear fitting
function is used as water transparency. Though the result is affected by PMT
gain, because if the PMT gain increases, Qobs (i) will also increase, but the
PMT gain can be measured by other calibration so this is still a useful method
to double-check the water transparency change.

45

6Data Set

6.1

Solar neutrino sample and reduction

6.1.1

Bad run selection and bad channel cut

In SK system, it is required to change the run by every 24 hours. And atypical run, which is 24 hours, is usually divided into ∼1100 subruns, whose time
length are ∼70 seconds (depending on the trigger rate). Before the reduction
process for solar neutrino data sample, the data runs are selected by removing
the bad runs to guarantee the data quality. The defination of a bad run is:
1. The run less than 5 minutes will be discarded because the pedestal data
is not enough.
2. Sometimes hardware or software problems occured and the DAQ system
was stopped. Then the corresponding subrun will be discarded.
3. When the calibration or maintenance is carried out, the DAQ system will
be changed to calibration run or test run, these runs will be removed fromthe solar data sample. Also in normal run, sometimes fake Supernova
burst generated by LED is used for training, the correponding subruns
are also removed.
4. Sometimes it is necessary to open the tank (for example, LINAC or Ni
calibration), thus the high voltage of the electronics will be turned off.
When the high voltage is turned on, the dark noise will be higher thanusual. Those subruns will be discarded until the dark rate become stable.
Additionally, SK-IV has been running for 10 years, some of the PMTs havebad performance. Currently there are 140(160) bad PMTs in ID(OD) detector.
For every month, an average number of bad PMT channels Nbad will be calculated, the subruns with bad channel Nbad < 10 or Nbad > Nbad + 1.5 × 24 areremoved. This is because most of the short runs do not have enough data toestimate bad channels so the calculated Nbad is small, and for some case, theQBEE board get out of order and the 24 PMTs connected to the QBEE board
are treated as bad channels.
6.1.2

Pre reduction

The events observed in SK is ∼ 2 × 105 per day. Most of them are backgrounds, with 295 elastic scattering events from solar neutrinos are expected in
one day. The size of raw data is ∼ 100GB per day, so it is necessary to remove
most of the obvious background events in the first reduction step before go into
the final sample. Though the cut criteria applied in this step is very loose, data
size and event number are largely reduced. The details of the cut criteria is:fiducial volume cut The fidualcial volume cut refers to the 2m distance cut
from ID wall. This is intended to remove the events whose reconstructedvertex position is close to the wall, because a good reconstruction quality
can not be expected near the wall, as the photon reflection on the PMT
and black sheet surface can cause the problem. The inner tank has aradius of 16.9m and a height of 18.1m, which means this cut removes the46

events with R > 14.9m or |z| > 16.1m, and keep the fiducial volume as
22.5kton, equal to 70% of the inner volume.
energy cut and OD cut
To keep the low energy range(4 ∼ 20M eV ) for solar neutrino, events withmore than 2000 hit PMTs will be removed, as those events are too energitic
and they are usually cosmic muons. For the same purpose, events with
OD trigger flag are also removed, because an expected elastic scattering
event of solar neutrino is impossible to make an OD trigger. The eventswith reconstructed energy lower than 4MeV are also removed, since low
energy background is dominate below that.
loose external gamma cut
The major background in low energy range is the gamma ray from therock around SK detector, or the radioactivity from the detector structure
itself. These background events have two features: The vertex are close
to the wall; The direction goes inward. To use these two features, firstly
the vertex and direction are reconstructed. Then the position Pwall is
calculated by going from the vertex to the wall through the opposite reconstructed direction. The distance def f , which is the distance from thevertex to Pwall , is used for loose external gamma cut as the events withdef f < 400cm are removed. Here we use the word ”loose” because more
tighter and optimaized cut will be applied in later process.
loose reconstruction quality cut
As mentioned before, the parameter gV is used to estimate the qualityof vertex reconstruction and the possibility of misfit. gd is for that of
direction reconstruction. In first reduction, events with gv2 − gd2 > 0.1 isremoved.
Time difference cut
This cut removes the event which is too close to the previous LE trigger by∆T < 50µs. Those events can be considered as three kinds of background:
1. The ringing effect from very high energy muon. This usually happens
in 1µs after the muon event.
2. Decay electron from cosmic muon. The half-life of muon is 2.2µs.
3. After pulse caused by PMT electronics. In time distribution, after
pulse usually makes a peak right after the real event, but it can also
make the peak as late as 15µs.

6.1.3

Spallation cut

The cosmic muons are energitic enough to break the oxygen nucleus andleave radioactive elements in tank water. Such radioactive impurity decays
in the water, and they refer to spallation backgrounds. The reactions occur
between cosmic muon and oxygen nucleus are complicated and can be writtenas:µ +16 O → µ + X + · · ·

47

(52)

Figure 33: Explanation of the defination of def f , pwall , fwall and θP M T .

48

Figure 34: The distribution of time difference to previous LE triggered event.
The red dashed line represents for the cut point of 50µsec cut. The peak around
1µsec is caused by ringing effect after high energy cosmic muons, the peakaround 15µsec is considered to be caused by after pulse [32].
Here X is the radioactive necleus produced by breaked oxygen. In the casethat hadrons(π ± ,n,p,etc) are knocked out, they can hit on other oxygen nucleus
and cause secondary or thirdary interactions. The radioactive products decay torelease β or γ rays, and especially some of the decay products are neary to thesolar neutrino energy range, which cause the difficulty to remove the spallation
background. The products already known are summarized in the table below.
The most long-life one can be up to 13.9sec.
The basic method for spallation cut is:
For each low energy event, find all the muons in previous 50µs window.
For each pair of muon and spallation candidate, calculate the time difference∆T and the distance ∆L from the vertex of spallation candidate to thereconstructed track of the corresponding muon.
For the corresponding muon, calculate the residual charge Qref by substracting
the expected charge along the reconstructed muon track. In Equation 53,
L is the length of reconstructed muon track, while in most cases the muon
pass through the tank and then L is equal to the distance of enter pointand exit point. Qunit is the expected charge in unit length of muon track.
Qunit is from the average of muon sample, by projecting all the hit PMTs
to the track. As the PMT gain increases by time, the parameter Qunit
also need to be updated to the lastest one.
Qres = Qtotal − Qunit × L

49

(53)

When spallation or oxygen nucleus breaking occurs, usually some additional photons can be expected, and thus the residual charge Qres will belarger than a clean muon.
Make the likelihood from ∆T , ∆L and Qres with Equation 54.
Lspa = P (∆t) × P (∆L) × P (Qres )(54)

The PDF for ∆L and Qres are made from a spallation sample whick is
selected by ∆L < 0.1 and E > 8M eV . The PDF for ∆T is generated
by considering all the known radioactive products. To set the cut pointof likelihood, the spallation sample is compared with a random sample in
Figure 37. The random sample consist of the pairs of a low energy eventwith E < 5M eV and a random muon. The reason for selecting random
sample by E < 5M eV is that below 5MeV the events are mostly causedby radioactive sources from the wall or the detector structure itself, but
not by spallation. As shown in Figure 37, in the lastest solar sample, thecut point for spallation likelihood is set at Lspa = 4.52 and 88.8%(20.0%)of the spallation-like(non-spallation) events are removed by this cut.

Figure 35: The list of Spallation production [38].

50

Figure 36: The increasing of Qunit by year.

Figure 37: The distriution of final likelihood result for random sample andspallation sample. The horizontal axis is the log of likelihood. By setting thecut point at spaloglike > 4.52, 88.8% of spallation sample are removed [34].

516.1.416

N cut and other cut

When a low energy µ− enter the tank, it can be captured by oxygen nucleus
and produce 16 N . 16 N has a half-life of 7.13s and goes beta decay in water, themain decay channel gives a 6.1MeV γ and a 4.3MeV β(66%) and a 10.41MeV
β(28%). Both of them are made from muon, but different from spallation, 16 N
backgrounds have the following features: the muon should be a stoppping one
in the tank, which means it has an enter point but no exit point; the vertex andtiming of the following low energy event should be close to the stopping muon.
16

N →16 O + e− + νe

(55)

Since a stopping muon is much more rare than a passing-througn one, asimple cut is more efficient than building a likelihood. So the cut criteria is set
as follows:
1. For each low energy event, pick up the previous muon by Qtotal > 1000p.e.,
as long as the constrain that the muon has no exit point.
2. Calculate the muon stopping position by the enter point and reconstructed
direction/momentum. Cut the low energy event if the vertex is in 250cm
range from muon stoppint position when their time difference is also lessthan 30sec.
Since this is the last step before the solar final samle, some other tight cuts
are also applied in this step, which includes:
1. PMT hit number cut
This cut removes the events whose total number of hit PMTs(Ntotal ) islarger than 400. This is because Ntotal =400 corresponds to ∼60MeV incase of a recoil electron, when Ntotal exceeds 400, it becomes too energitic
for solar neutrino.
2. Tight fiducial volume cut
Though the 2m distance cut from the ID wall has already been applied infirst reduction step, in 3.5-5.0MeV energy range, there still remains background near the bottom and the barrel. As an non-uniform distributed
background will cause a large uncertainty in solar analysis, additional vertex cut is applied in 3.5-5.0MeV range:z > −7.5m (4.5 − 5.0M eV )(x2 + y 2 ) + (150.0
4
× |z − 4.25| ) < 150.0 (3.5 − 4.5M eV )11.754(56)(57)

After the tight cut, fiducial volume is 8.85kton and 16.45kton for 3.54.5MeV and 4.5-5.0Mev respectively. There is no additional cut needed
above 5MeV.
3. Tight external cut
As discussed in the part of tight fiducial volume cut, the radioactive background is non-uniformly distributed and most of them are near the wall.

52

Those low energy background can be considered from the decay of Rn
daughters.
To reduce the external signals, the same parameter def f is used for cut
above 5MeV:

def f > 650cm (5M eV < E < 7.5M eV )def f > 400cm (E > 7.5M eV )(58)

Which can be known from Figure 38, the radioactive backgrounds in 3.55.0MeV concentrate near the bottom and barrel. Thus the additional cut
in this range is defined by pwall . pwall means the position on the wallwhen tracking from the reconstructed vertex along the opposite direction.
pwall is categorized into top, bottom or barrel, and the cut is applied as: def f > 1000cm (pwall = top)
def f > 1200cm (pwall = barrel)
def f > 1300cm (pwall = bottom)(59)

Figure 38: The vertex distribution for 3.5 ∼ 4.0 MeV(left), 4.0 ∼ 4.5
MeV(middle) and 4.5 ∼ 5.0MeV (right) in kinetic energy [32].
4. Cluster hit cut
When a radioactive event occurs on the PMT surface or FRP cover, it cancause hits on neighbour PMTs. Usually the number of hit PMTs is small
but sometimes the event accidentally coincides with dark noise hits andit can be recorded in data. Such event has the features of: Some of thehits clusters near one PMT; About half of the hits are from dark noise, so
the sharpness of timing distribution is worse than a real neutrino signal
event.
To cut these events, two parameter are defined:(a) R02 : The minimum radius containing more than 20% of hit PMTs
within 20 nsec time window.
(b) N20rawT : The maximum number of hits in a 20ns timing window(without
subtracting Ttof ).
For convienence, another parameter Clik is defined as:
Clik = R20 × N20rawT /Nef f
53

(60)

By comparing the background with solar neutrino MC simulation, the cut
criteria is set as:

Clik < 75.0, r2 > 155m2 , z < −7.5m (4.5M eV < E < 5.0M eV )Clik < 75.0, r2 > 120m2 , z < −3.0m or z > 13.0m (3.5M eV < E < 4.5M eV )(61)

5. Tight event quality cut
Tight event quality cut uses the same parameter(gv and gd ) with first
reduction step, and the cut criteria is optimized with MC simulation as a
Signal.
signal to ensure the maximum significance Signif icance = √Background 2
gv − gd2 > 0.29 (3.5M eV < E < 5.0M eV )g 2 − gd2 > 0.25 (5.0M eV < E < 7.0M eV ) v
gv2 − gd2 > 0.20 (E > 7.0M eV )(62)

6. Hit pattern cut
Solar analysis is searching for the recoil electron from the neutrino, so it
is necessary to find a way to distingush γ ray or signals with multiple
rings. The γ ray can induce Compton scattering by many times and it
will final give a very dirty ring, which actually consist of many rings fromthe scattered electron. While for a recoil electron event, it usually has aclean ring and the vector from vertex to hit PMT has ∼ 42◦ angle to thetrack. To use this feature, likelihood for hit patter is defined as:N

Lpattern (E, ~v ) =50
1 X
log(Pi (E, cosθpmt , fwall ))
N50 i

(63)

Here N50 is the maximum number of hit PMTs in 50ns time window, E isreconstructed energy, θpmt is the angle between the track and the vector
from vertex to hit PMT. fwall is the distance from vertex to nearest wall.
Pi is the PDF made from MC simulation of single electron events. By
maximum the significance, the cut point is set as: Lpattern > −1.88 (6.0M eV < E < 7.5M eV )Lpattern > −1.86 (7.5M eV < E < 11.5M eV )
Lpattern > −1.95 (E > 11.5M eV )6.2

Relic neutrino sample and reduction

6.2.1

Spallation cut(64)

The spallation background has already been introduced in solar section, but
the spallation cut method in relic analysis is different from that in solar analysis.
The reason is solar analysis concentrate on the energy range of 4∼20MeV, while
relic analysis is looking for neutrino between 16∼30MeV. Since the energy range
is different, so the expected radioactive isotopes are also differenet. This means
the likelihood built for solar analysis can not been used for relic sample directly.
Thus, spallation cut method for relic is developped independently. The basic
54

Figure 39: Histogram of r02 vs N20rawT . The left is the background sample andthe right is the solar simulation result [32].

Figure 40: Distribution of hit pattern likelihood. Black is for real data, and redis the simulation of solar neutrino, for 5.5 ∼ 7.5MeV in left, 7.5 ∼ 11.5MeV in
middle, and 11.5 ∼ 19.5MeV in right.

55

thinking is similar, but the likelihood and the variables used for likelihood aredefined in another way.
The highest energy spallation product believed to be produced in the detector is 11 Li and 14 B, each with beta decay end point energies of up to 20.6 MeV.
Because of energy resolution effects, events with reconstructed energy up to 24
MeV are reasonable to be suspect to be spallation background.
Luckily, 11 Li and 14 B have short half-life of around 0.01 seconds, so they
are rather easy to associate with the muon which created them. But when
we go lower in energy, the number of possible isotopes increases rapidly, as
do their half-lives. The basic idea to remove the spallation events is similar
to solar analysis, to check the correlations between the electron candidate andthe proceeding muon tracks. There are four variables used in searching for the
spallation activity:
1. ∆T . The time difference between the electron candidate and the proceeding muon. This is very powerful in identifying short-lived spallation
products. The smaller the ∆T , the greater the likelihood that an event
being spallation.
2. Ltran . The transverse distance from the electron candidate to the muontrack. The spallation products are created along the muon track and ingeneral they can not travel far before decaying. A short Ltran indicatesmore possibility of spallation.
3. Llong . Firstly, the point where the electron candidate lies on the muontrack will be found. Then the peak position of dE/dx plot on muon track
will be found. Llong is the distance between the two. This is because thespallation occurs on muon track intends to make a peak on dE/dx plot.
Llong is not as powerful as Ltran , but still useful. Since the physics behind spallation has not been fully understood, especially that the hardron
physics occur in spallation process can be complicated. So the spallation
longitudinally may deserve further investigation.
4. Qpeak . The amount of light seen in the peak of dE/dx plot in a bin of
4.5m width. The strong peak Very strong peaks, which show that lotsof light was deposited in the detector originating from that region of themuon track, usually indicate spallation.
Since the relic candidates is not as much as solar candidates, Brute Force
Fitter(BFF) is also used for muon reconstruction. BFF is also a muon fitter but
the reconstructed vertex will be tried in all of the tank by steps which is smallenough. BFF cost much more time than Muboy, so it is not unrealistic to be
applied to the large number of solar candidates, but can be tried on the limited
number of relic candidates.
The result of Muboy will be first tried, when the recontruction precision of
Muboy is poor(Goodness < 0.4), the result of BFF will be used instead. If
both of Muboy and BFF(Goodness < 0.3) are not trustable, all the events in
2s window behind the muon will be banned as dead time. Additionally, whenthe muon is too energitic(totallight > 400, 000p.e.), all data after the muon will
also be rejected for a few seconds. If both of Muboy and BFF are good, thenboth of their result will be input to calculate the likehood and either of the twolooks like spallation, the event will be rejected.
56

Figure 41: Explanation of the defination of LT RAN and LLON G .

57

In order to build the likehood, a spallation sample and a non-spallation
sample are needed. Since the amount of spallation events are much more larger
than the real relic events we can expect, all the proceeding muons occuringwith 30 seconds before the candidate events are used as the spallation sample.
And the muons occuring within 30 seconds after the candidate events are usedas the random sample as they are impossible to have spallation correlations
with the candidate events. Distributions of ∆T , Ltran , Llong , Qpeak are madeindependently and be normalized to make PDFs for the four variables.

Figure 42: PDF of the four varibles. Red is random sample, and black isspallation sample [63].
Muboy categorize the muon as four kinds: single through-going, stopping,
multiple, and corner clipper. The likelihood is built seperately in different energy
range, as well as for the four categories as:
Lspal = log[

spal
PQspalPLspalPLspal
P∆T
P EAK
T RAN
LON G
×
×
×
]
rand
rand
rand
P∆T
PLrandP
P
LLON GQP EAK
T RAN(65)For muon with multiple track, the likelihood is also seperated for each track
and the one with largest Lspal will be used for cut. The likelihood is the primary
way that spallation events are eliminated, and some cuts on specific quantities
also exist.
6.2.2

2peak cut

Although most decay electrons are eliminated by the 50µsec pre cut, occasionally it occurs that the muon decays quickly enough that both the muon andthe decay electron information are captured in the same 1.3µsec timing windowthat comprises a single event. In these cases, the timing information for thatevent will have a double peak structure. Other, even rarer events can cause58multiple peaks as well. The cut purposed to remove such kind of events is theso-called 2peak cut or pre/post activity cut.
In order to find and reject such events, the timing information of the event
was searched for any indication of a peak outside the main timing peak. The
algorithm uses the BONSAI fit time to define the main peak. Any hits withtime of flight subtracted times earlier than 12ns before main peak is considered
as possible pre-activity; any time after 20ns after the BONSAI event time isconsidered post activity.
For each of the 3 regions (pre-activity, main, and post-activity), the peakis searched by scanning the maximum hit inside a 15ns time window. The
maximum hit number in each region is recorded as N15 . In the final stafe, if
N15 in pre-activity region is more than 12, or N15 in post-activity region is morethan 15, then the event is rejected.

Figure 43: An example event withe 2 peaks.

6.2.3

Multiple ring cut

Atmospheric neutrino events can produce multiple charged particles by theinteraction with water. As the multiple charged particles may share the sametiming peak, such events can not be removed by the pre/post activity cut. These
multiple charged events can have multiple Cherenkov rings.
To remove these specifically, we use the fact that these multi-particle events
have multiple Cherenkov rings, A ring counting method has already been developed for SK atmospheric neutrino analyses utilizing Hough transforms, We
used the same method to look for multi-ring events. The software determines
an angle between the rings we can use. In rare cases overly fuzzy electron rings
can be mistaken for two different rings by the software; for this reason events
determined to be multi-particle, but separated by less then 60 degrees, were59

kept, while events with multiple rings separated by more then 60 degrees were
rejected.

Figure 44: An example event with multiple rings.

6.2.4

Cherenkov angle cut

For a typical single electron event, or other Charged particles traveling close
to the speed of light in pure water emit Cherenkov light in a cone shape with an
opening angle of about 42 degrees. While for low momentum heavier particles,
such as muons, pions, and nucleons created by atmospheric neutrinos, they canalso be reconstructed in the relevant energy region by tagged as an electron.
However, in that case those particles produce Cherenkov light with an opening
angle that is smaller then 42 degrees. For γ ray, since multiple scattering canconfuse the ring, so the Cherenkov angle distribution is not as sharp as we could
hope. As a result, the reconstructed Cherenkov angle for γ event is usually muchlarger than 42 degrees.
Fig 45 shows the Cherenkov angle distribution of SN relic MC (LMA model).
The main peak is determined to be between 38 and 50 degrees, and this is the cut
applied to the data to eliminate remaining heavy particles at lower angles, andmore isotropic events at higher angles. As the signal region is defined between
38 and 50 degrees, the Cherenkov angle cut is also useful in further removing
spallation events that are due to a β decay plus a γ emission.
6.2.5

Pion cut

When it goes to relatively high energy range(E > 30M eV ), most of theevents are actually pions created by atmospheric neutrino interactions. Since
the pions are much heavier and they are usually quickly captured by Oxygenand disappear, pions intend to give the Cherenkov ring to be cleaner and sharper
60Figure 45: Distributions of the Cherenkov angles for the data(black) and SRN
simulation(red).

61

than electron, because when an electron is travelling in water it will suffer frommultiple scattering and the direction will be slightly changed, which results ina relatively poor and dirty ring, but much better than a γ signal.
The ’fuzziness’ or ’sharpness’ of the Cherenkov ring can be considered for discriminating pions and electrons. To quantitatively estimate that, the distribution of the cone angles from three-hit combinations which is made in Cherenkov
angle fitting can be used. This distribution will be narrower for pions thenfor electrons, though the two peak at the same angle. Using this logic, a pion
likelihood variable was constructed, where pion likelihood is defined as follows:
◦

◦

±3
±10
LP ION = Nentries/Nentries
◦

(66)
◦

±3
±10
Nentries
is the number of entries in ±3◦ range from the peak. Nentries
is the◦same meaning for ±10 . By looking at SN relic MC compared to pion MC, thecut criterion was determined to be 0.58, which makes the cut approximately 1%
inefficient.

6.3

Atmopheric neutrino sample

The high-energy data sample has three different categories: fully-contained(FC), partially-contained (PC), and upward-going muon (UPMU) [66]. FC
neutrinos have reconstructed interaction vertices inside the fiducial volume of
the inner detector, combined with low light levels in the outer detector. PC
neutrinos also have interaction vertices inside the fiducial volume, but havesignificant light in the outer detector volume indicating exiting particles. UPMU
neutrinos are the highest-energy SK sample; they result from muon-neutrino
interactions in the rock surrounding the detector, which produce penetrating
muons. These muons either stop in the inner-detector volume as stopping events,
or go through the inner detector as through-going events. The energy range forneutrino parents in FC and PC sample is 100 MeV–10 GeV, and for UPMU it
is 1.6 GeV–100 PeV.

62

7

Neutrino search for GW170817

7.1

GW170817 overview

On August 17th 2017 at 12:41:04 UTC, the Advanced LIGO and Advanced
Virgo experiment identified the first evident signal of a gravitational wave fromthe binary neutron star merger, named GW170817 [2]. The interpretation isa merger of two compact objects consistent with neutron stars having total
system mass of 2.74 solar masses and a luminosity distance of 40 Mpc. Associated with this gravitational wave signal, the Fermi Gamma-ray Burst Monitor and International the Gamma-Ray Astrophysics Laboratory also detected ashort gamma-ray burst, GRB170817A, which has a consistent location with themerger and a 1.7-s delay to the merger time [1]. Subsequent extensive electromagnetic follow-up observations in ultra-violet, optical and infrared wavelengths
were performed. These observations led to the conclusion that the merger happened in galaxy NGC4993 and was followed by a short gamma-ray burst anda kilonova/macronova [3, 4]. High-energy neutrino signals associated with themerger were also searched for by the ANTARES, IceCube, and Pierre Auger
Observatories. It was concluded that no significant neutrino signal was observed [5].
A search for neutrinos in Super-Kamiokande (SK) associated with this gravitational wave signal produced by the binary neutron star merger in NGC4993
will be reported in this section. The analysis method is similar to that for theprevious neutrino search in SK for GW150914 and GW151226 [6].
As mentioned in previous sections, in SK analysis neutrino events with reconstructed energies above 100 MeV are categorized as the ‘high-energy data
sample’, which are typically used to study atmospheric neutrinos and to searchfor proton decay. Neutrino events with reconstructed energies below 100 MeV
and down to 3.5 MeV are categorized as the ‘low-energy data sample’. They are
typically used to study solar neutrinos and to search for core-collapse supernova
neutrinos. The directional determination accuracy varies according to sample
and direction, but can be as accurate as ∼1 degree for upward-going muons.
Some theoretical predictions of neutrino emission mechanism via binary neutron star mergers have been proposed; for example, some fraction of the kinetic
energy in relativistic ejecta from gamma-ray bursts could convert to high-energy(∼ 1014 eV) neutrinos, or a similar mechanism as for core-collapse supernovae
could produce few-tens-of-MeV neutrinos [7, 8], the expected fluence is roughly
estimated to be 104 cm−2 for 10 MeV neutrinos within 1 second after merger.
Neutrino observations associated with a binary neutron star merger using theunique characteristics in SK would validate such proposed mechanisms. Coincident search was made in the full data sample using the same time window
as ANTARES-IceCube-Pierre Auger, i.e., ±500 s around the merger time andin a 14-day time window relevant for longer-lived emission processes [5]. The
primary background events for this search in the high-energy data sample arealmost entirely atmospheric neutrinos, while radioactive impurities, spallation
products from cosmic ray muons, atmospheric and solar neutrinos are the main
backgrounds in the low-energy data sample. We note that SK carried out a
LINAC calibration from August 3-22, 2017. Fortunately, physics data-taking
operated at the time when the neutron star merger occurred; however, there
were unavoidable radioactive impurities adhered on the surface of the LINAC

63beam pipe present in the low-energy data sample.

7.2

Search method and result

As mentioned before, the energy range for neutrino parents in FC and PC
sample is 100 MeV–10 GeV, and for UPMU it is 1.6 GeV–100 PeV. All thethree event topologies are considered for this search. Further information about
the selection and reconstruction methods for the three categories can be foundin reference.

FC
PC
UPMU
FC
PC
UPMU

observed num. of eventin ±500 s000
following 14 days for all sky
76 ± 8.72
8 ± 2.83
following 14 days for 5◦ solid angle0

expected num. of event(9.36 ± 0.06) × 10−2
(7.52 ± 0.23) × 10−3
(1.64 ± 0.02) × 10−2
91.44 ± 0.57
7.35 ± 0.23
(6.11 ± 0.04) × 10−2

Table 7: The numbers of expected and observed events in FC, PC and UPMU
data sets, respectively, for a ±500-s time window around GW170817 and for
14 days after GW170817. The errors √on the observed number of events in thefollowing 14 days for the entire sky are N . For UPMU event search in following
14 days, the event number in a solid angle of ±5◦ around NGC4993 was shown
instead. The energy range for FC and PC is 100 MeV–10 GeV, and the energy
range for UPMU is 1.6 GeV–100 PeV. The total livetime for the following 14
days is 11.30 days.

A ±500-s window search around the LIGO detection time of GW170817 has
been conducted. Since the observation of optical light and gamma ray burst
lasted for 14 days, a 14-day window search following the GW detection was also
carried out. These two kinds of time windows are consistent with those neutrino
experiments in reference.
±500-s window search In the ±500-s window around GW170817, no neutrino event was found in the FC, PC, or UPMU data sets. This null result
is used in the calculation of the upper limit on neutrino fluence in the next
sections.
For lowe energy signals, there are 7 events found after solar reduction(3.5 ∼16M eV ) and no event was found after SRN reduction(16 ∼ 100M eV ). However, the LINAC calibration was carried out in that period. At the timing of
GW170817, beam was not running but the pipe was inside the tank.
In Super-K regular calibration, the events reconstructed inside 2m spacefrom the calibration source are removed for solar or relic neutrino search. Afterchecking the vertex of the 7 events in solar sample Fig 46, all the 7 events foundin solar sample are considered to be from the impurities of the pipe surface andthus no signal was for GW170817.
64Figure 46: The left shows the timing(horizontal axis) and reconstructed energy(vertical axis) of the 7 events found in ±500-s window(shadow box) of solar
sample. The rightpis the reconstructed vertex of the 7 events, while the horizontal axis is r = x2 + y 2 and the vertical axis is z. The yellow shows the 2m
space from the center of LINAC pipe. All the 7 events found are considered tobe from the impurities of the pipe surface.
14 days window search Since the LINAC calibration was carried out in that
period, there are too many noise events in low energy range. These noise events
are not only from the beam, but also can be from the impurities of the pipe or
the electronics noises, because the tank was opened and closed very frequently,
as well as the pipe was moved in and out, and the hardware electronics wereturned on and off. Thus, it is difficult to make physics study on the low energy
sample for the 14 days window and the result will not be used for this search.
For high energy sample, the event numbers found in 14 days window issummarized in Table 7. The expected number of events in Table 7 based on
2976.01 days of SK data. The livetime for the high-energy analysis is 11.30 days,
after removing LINAC beam runs in the following 14 days after GW170817.
Unlike the FC and PC samples, the UPMU sample only contains upgoing
muons, so it is sensitive to only half of the sky. In 60.4% of the following 14days, NGC4993 is within the sensitive half. Since the direction of NGC4993 iswell known [3,4], for UPMU data, for which the angular resolution is better thanfor the other two samples, we concentrated on a ±5◦ cone around NGC4993 forthe event search in the following 14 days. This method was previously used in
SK to search for neutrino signals associated with astrophysical objects [9]. The
5◦ constraint was not used for the ±500-s search in Table 7 because no event
was observed in all sky during this window, and unlike the 14-day-window case,the zenith angle change of NGC4993 in ±500-s can be ignored. All the results
listed in Table 7 are consistent with our expected event rates and no significant
signal was found for GW170817.

7.3

Fluence limit

Though there is no event observed within a ±500-s window, either in lowenergy data nor in high-energy data, the null number can be converted to an
upper limit on neutrino fluence. The calculation for neutrino fluence limit isdone separately for the low-energy, FC+PC, and UPMU data sets. The calculation uses the same procedure laid out in [10], which follows from [11].
FC+PC (100M eV ∼ 10GeV ) For the FC and PC data set, the neutrino flu65

ence can be calculated using equation (67),
ΦF C,P C =

NT

R

N90
,
dEν σ(Eν )(Eν )λ(Eν−2 )(67)

N90 is the 90% C.L. limit calculated from a Poisson distribution, for the observed neutrino events in a ±500-s window with the expected backgroundnumber. Since there is no neutrino event found in a ±500-s window, N90
can be fixed as N90 = − ln(0.1) = 2.3. NT is the number of target nuclei
relevant to the neutrino interactions, while for Charge Current QuasiElastic reaction and Neutral Current Quasi-Elastic reaction, the target is
Oxygen. σ is the combined cross section for all interactions as plotted in
Fig 47.  is detection efficiency as shown in Fig 48, and λ is the densityof Eν assuming an energy spectrum with index of −2 as shown in Fig 49.
This spectral index is commonly assumed for astrophysical neutrinos accelerated by shocks [12].

Figure 47: The cross section for all interactions combined in the range of100M eV ∼ 10GeV . Colors represents the four differenet neutrino types included in FC and PC data.
Fluence limits are calculated separately for each neutrino type because
the cross section and detection efficiency depend on neutrino type.
Cross sections in equation (67) are from NEUT 5.3.5 [13]. NEUT 5.3.5 isalso used to produce mono-energetic neutrino interactions in the SK Monte
Carlo detector simulation in order to determine the detection efficiency.
low energy (4M eV ∼ 100M eV ) The fluence calculation for low-energy neutrinos uses an expression similar to (67) but with different energy spectra.
Here we assume two kind of spectra, one with an index of 0, i.e., a flat
spectrum, and another being a Fermi-Dirac distribution with average energy of 20 MeV. The two kind of spectrum are plotted in Fig 50.
66

Figure 48: Detection efficiency for the four neutrino types in FC and PC data.
Left is for FC, and fight is for PC. It need to be noticed that the Cherenkov
threshold for muon is 157.4MeV so there is no sensitivity for νµ and ν µ near
100MeV in FC. To produce PC events, usually more energetic neutrino is needed
so the sensitivity is near 0 below 1GeV.

Figure 49: The energy spectrum used in FC+PC fluence limit calculation. It is
assumed with a index of −2 and normalized to make total are to be 1.

67Φlowe =

NT

R

N90
,
dEν λ(Eν )σ(Eν )R(Ee , Evis )(Evis )(68)

R is the response function to convert electron or positron energy (Ee ) to
kinetic energy in SK (Evis ). R is estimated by SK detector Monte Carlosimulation explained in previous section. (Evis ) is the detection efficiency.
(Evis ) depends on the detected energy, and was estimated for 4 ∼ 20M eV
by solar reduction method and 20 ∼ 100M eV by SRN reduction method,respectively. The combination of R(Ee , Evis ) and (Evis ) is shown in
Fig 51. Again, no event was observed in ±500 s, so N90 is 2.3.

Figure 50: The energy spectrum used in low energy fluence limit calculation.
A flat spectrum and a Fermi-Dirac distribution with average energy of 20 MeV
are used. Both the two are normalied to make total area to be 1.
We also express the fluence limit which is calculated for monochromatic
neutrino energy Eν . The fluence limits at various energies are shown in
Figure 52.
UPMU (1.6GeV ∼ 100P eV ) For the UPMU data set, the neutrino fluence iscalculated using equation (69),ΦU P M U =N90
R.
Aef f (z) dEν P (Eν )S(z, Eν )λ(Eν−2 )(69)

The fluence of UPMU events depends on zenith angle. Aef f (z) is thezenith-dependent effective area (shown in Fig 54), where z is the the zenith
angle of the incoming neutrino. P (Eν ) is the probability for a neutrino
to create a muon which is higher than Eνmin and thus can be detectablein SK. P (Eν ) for νµ and ν µ are shown in Fig 55 respectively. S(z, Eν ) is
the shadowing of the neutrinos due to interactions in the Earth (shown68

90% fluence limit(cm-2)

Figure 51: Detection efficiency for lowe energy range. This refers to the combination of R(Ee , Evis ) and (Evis ) in Equation 68.

νx
νx
νe
νe

1012
1011
1010
109
108
107
0

20

40

60

80

100
total energy(MeV)

Figure 52: The 90% C.L. limits for GW170817 events on fluence obtained formono-energetic neutrinos at 4 MeV, 7 MeV, 10 MeV, 14 MeV, 20 MeV, 30 MeV,
50 MeV, 80 MeV and 100 MeV.

69

in Fig 56). S(z, Eν ) also depends on zenith angle because when cos(z) isgetting close to 1, the travel distance in the earth is getting longer andless neutrino can reaches SK detector. Similar to FC and PC limits, λhere is the number density of Eν in a spectrum with index of −2.

Figure 53: The E −2 spectrum used in the calculation for fluence limit of UPMUdata.
The results of fluence limits for FC+PC, UPMU, and low-energy data aresummarized in Table 8. The UPMU upper limit fluence values range from (14–37) cm−2 for neutrinos and from (18–50) cm−2 for antineutrinos, depending
on zenith angle from 90◦ to 0◦ . The zenith-dependent upper limit of neutrino
fluence from UPMU events are shown as a sky map in Figure 57. As mentioned
before, UPMU data set collect the upgoing muon events and thus it only has
sensitivity to half of the sky. Fortunately at the detected timing of GW170817,
the source NGC4993 is in the sensitive half.
−2
To focus on the direction of NGC4993, UPMU limit is 16.0+0.7
and−0.6 cm+1.1−2
21.3−0.8 cm for neutrinos and antineutrinos, while the error is calculated bya range of ±5◦ around the zenith angle of NGC4993.
We note that the present study is sensitive to neutrinos between 1.6 GeV
and 100 GeV, which is not covered in other searches [5]. Our UPMU data may
be compared or combined directly with that of other neutrino telescopes. We
provide an UPMU fluency limit in Figure 58. It need to be stressed that SK has
the best limit for neutrinos below 100MeV, which is the expected energy rangeof thermal neutrinos from binary neutron star merger.
Considering dGW as the distance from the detector to NGC4993, our upper
limit on fluence of UPMU data can be converted into an upper limit on total
radiated energy in neutrinos, by weighting by 4πd2GW in equation (69). The
resulting upper limit on total energy is Eνtot ∼(1–6)×1053 ergs for GW170817
assuming the luminosity distance of 40 Mpc.

70

Figure 54: Effect area which depends on zenith angle.

Figure 55: The possibility for a neutrino to produce a detectable muon in SK,
which refers to P (Eν ) in Equation 69. Both the detection efficiency and crosssection are included.

71

Figure 56: S in Equation 69, which can be considered as the possibility for theneutrino to interact in the earth. The left is for νµ and the right is for ν µ . S
depends on neutrino energy and zenith angle, as different zenith angle means
different distance for neutrino to travel from the back of the earth.

Figure 57: The 90% C.L. limit on fluence for neutrinos(left) and antineutrinos(right) in UPMU data set, overlaid with the 90% C.L. contour for the location of GW170817 according to LIGO and Virgo released data (solid red line).

72

1012
11

E2F [GeV2cm-2]

10

νµ
νµ

1010
109
108
107
106
1.6

10

102

103

104

105

106

107

108

E [GeV]

Figure 58: The 90% C.L. limits for the UPMU data set during a ±500-s window,
in the direction of NGC4993. Limits are calculated separately for each energy
range, by assuming a spectrum with index of −2.

νµ
ν̄µ
νe
ν̄e

ν̄e
νe
ν̄x
νx

GW170817 Φν (cm−2 )from FC+PC only
5.6 × 104
1.3 × 105
4.8 × 104
1.2 × 105
from low-energy only
flat spectrum
1.2 ×107
1.0 ×109
7.5 ×109
6.3 ×109

from UPMU only
16.0+0.7−0.6
21.3+1.1−0.8
Fermi-Dirac with Eave =20 MeV
6.6 ×107
3.4 ×109
2.6 ×1010
2.1 ×1010

Table 8: Limits at 90% C.L. on the fluence of neutrinos from GW170817 given aspectral index of −2 and a range of 100 MeV–10 GeV for FC+PC and 1.6 GeV–
100 PeV for UPMU data. The error of UPMU limit is made with ±5◦ range
around zenith angle of NGC4993. Low-energy limits assume a flat spectrum
as well as a Fermi-Dirac spectrum with Eaverage =20 MeV from 3.5 MeV to
100 MeV. νx (ν̄x ) represents νe (ν̄e ) and νµ (ν̄µ ).

73

8De-excitation gamma search from Charge Current Quasi-Elastic reaction

8.1

Motivation and Plan

This research focus on the search for de-excitation γ from CCQE interaction of atmospheric neutrino and the direct measurement of the Branchingratio(Br(γ)). De-excitation γ has been briefly introduced in Section 1. To explain the motiavtion more clearly, when an atmospheric neutrino induce CCQE
interaction with Oxygen inside SK, a muon and a possible ∼ 6M eV γ will be
produced. As the product muon can be above or below Cherenkov threshold
and the de-exicitation γ can be given or not given, the final signal observed in
the detector can be considered to have four categories:
1. The muon energy is below the Cherenkov threshold and thus invisible,and no de-excitation γ is emitted. In this case, no prompt signal can beobserved, and the only signal is the decay electron from the invisible muon,but since there is no visible signal to tag the electron, it is impossible
to know whether it is a decay electron from muon or an electron fromneutrino, and thus this will be a unreducable background component for
low energy neutrino analysis.
2. The muon is invisible, but de-excitation γ is emitted. In this case, firstly
the de-excitation γ is observed, after the invisible muon decays with ahalf-live of 2.2µsec, there will be a decay electron. By searching for thepair of a prompt γ and a delay electron, a sample of the de-excitation γ
can be obtained
3. The muon is energetic enough to give Cherenkov light in SK water, butthere is no de-excitation γ emitted. Since the muon is produced inside
the tank, such kind of events are saved in FC set of ATMPD data sample,
and the decay electron is saved in the −5 ∼ +35µsec window.
4. The muon is visible, and de-excitation γ is emitted. similar to 3 case, thiskind of events is saved in FC data set and the decay electron is saved inthe −5 ∼ +35µsec window.
For case 3 and 4, a coincident search for a delayed electron can ensure thatthe muon is from CCQE interaction. Additionally, by search the γ signal inside
the prompt muon event, case 3 and case 4 could be distingushed and thus thebranching ratio of de-excitation γ could be obtained. A possible method which
can be used for γ search inside muon will be introduced in following section.
For case 1, only decay electron can be observed, thus it is impossible todistingush it from the expected electron signal from either the IBD interaction or
the elastic scattering. The energy of decay electron follows the Michel Spectrum
and can be up to ∼ 60M eV , which overlaps the search window of Supernova
Relic Neutrino(16 ∼ 30M eV ). As a result, this is an unreducable backgroundcomponent in SRN search. Besides, this component does not affect the solar
neutrino analysis, because solar neutrino search focus on 4 ∼ 20M eV window
which is rather little in Michel Spetrum, and further more they could be removed
by solar angle cut(cos(θsun ) > 0.8).
74

The background component from case 1 is the motivation of this study.
This is one of the two most difficulties in SRN search(another is spallation
background), and can also make an uncertainty term in long base long neutrino
experiment. Though this kind of background is not countable, if we can directly measure the branching ratio of de-excitation γ, then it will be possible toestimate this background component.
For case 2, a pair of a prompt γ and a delayed electron is expected. By thecoincident search of the prompt and delayed signal, this case is countable and asample of de-excitation γ can be made. Supposing the Br(γ) is already known,
by counting case 2, case 1 can be esimated.
To specify a detailed plan for this study, the following steps are considered:
1. Use simulation to establish the method of searching for a γ inside a muon
event.
2. Search for the coincident signal of a prompt γ and a delayed electron to
make a sample of de-excitation γ events. This refers to the case 2 above.
3. Considering the flux of atmospheric neutrino, the theoretical expectationof CCQE cross section, and the efficiency of SK detector, including all of
those into the caculation to give an expectation of CCQE event number
observed in SK according to the live time. Compare the expected number
with the de-excitation γ sample made from case 2.
4. By useing the method of searching a γ inside the muon, and the deexcitation γ sample, apply them onto the muon events saved in FC data
set and finally give the Br(γ) of de-excitation γ in CCQE interaction andestimate the background component refers to the case 1 above.
Details of the first step and second step will be explained in the following
sections. It need to be mentioned that our target is the Br(γ) of invisible
muon, however Br(γ) obtained by this way is the one of visible muon. This canbe solved by slicing the muon energy into every 50MeV or 100MeV bins. To
calculate the Br(γ) for each bin, the Br(γ) for invisible bin could be estimatedby considering proper model.

8.2

Search Method for Branching Ratio of de-excitation
gamma

8.2.1

Simulation setup

Simulation is made by SKdetsim, the SK official Monte Carlo, which is introduced in early sections. In SKdetsim neutrino interactions are modeled by
NEUT [13] and hadronic interactions are calculated by GCALOR. As explained
before, detector dependent parameters as timing resolution, photon acceptanceon PMT, water transparency, reflection on PMT and black sheet surface, are
input by custom SK code, based on real measurement. The performance of
SKDETSIM have been well tested and verified, and the energy scale are determined by electron linear accelerator LINAC from 5MeV to 18MeV [14]. Since
our purpose is atmospheric CCQE de-excitation γ, we simply generate a muon

75

from the random position inside SK tank fiducial volume, with a 6MeV γ generated from the same position by the same timing, in order to avoid the dependence of neutrino flux and cross section. We hope to use weak muon near
Cherenkov threshold, because our final motivation is to get the Br(γ) of inivisible muon. However when the muon is getting close to Cherenkov threshold, thering angle become small and thus the reconstruction precision become low. So
here we take 300MeV/c as a typical value of muon product from CCQE interaction of sub-GeV neutrino, and the γ is generated uniformly to all directions
at 6MeV.
8.2.2

Light profile of mu and gamma

The number of Cherenkov photons from a sub-GeV muon can be from a
few houndreds to a few tens thousands, while for a 6MeV γ, the number is a
few tens. Typically for a 300MeV/c muon, the muon will induce 500 ∼ 600 hit
PMTs, and a 6MeV γ will induce only ∼ 30 hits.
In order to search for the limited γ-induced hits from muon-induced hits,
the light profile need to be studied. Considering the photons finally reach PMT
surfaces, the origin of the photons can be the following process : Cherenkov
photon prodcution, delta rays, scattering in water, reflection on PMT or black
sheet surface. A typical simulation event of the pair of 300MeV/c muon and
6MeV γ is shown in Figure below. The hit PMTs are tagged with the source ofthe photons, while for multi-photon hits, the PMTs are tagged with the earliest
photon detected.

Figure 59: A typical simulation event of the pair of 300MeV/c muon and 6MeV
γ. The hit PMTs are tagged with the source of the photons. In this event, theangle between muon and γ is 154.7◦ .
From the display, the muon-induced Cherenkov photons intend to be a ring76

and the γ-induced photons intends to be a cluster. As explained in early section,
this is because a γ will surfer from Compton scattering and electron/position
pair production for several times, the finally detected one is the Cherenkov lightof the electrons so the signal consist of several electron rings and intends to becluster rather than a clean ring.

Figure 60: The upper two are the distance from reconstructed vertex to realvertex. The down two are the angle from recontructed direction to real direction.
Left two are the simulation of 400MeV/c muon with 6MeV γ generated fromthe same position. Right two are the simulation of a single 400MeV/c muon.
From the comparison of left and right, it can be found that a 6MeV γ from thesame position will not affect but will improve the reconstruction precision by alittle bit.
Besides, an additional γ will not affect the reconstruction of muon event,both for vertex and direction. On the contrary, more hits from the same position can help to fit the starting point and timing, thus finally improve thereconstruction a little bit.
Supposing the muon direction is already known by reconstruction, a cosθ −(T − T OF ) map as shown in Fig 62, can be used to deeply study the light
profile. Here the points on the map means each hit PMT, and cosθ is the angle
between the muon direction and the vector from vertex to corresponding PMT.
The vertical axis is the hit timing of each PMT with TOF substracted. Similar
to the event display, the hit PMTs are tagged with the the source of the earliest
photon.
From the event displayed and the cosθ − (T − T OF ) map, we could see
definitely the easiest one to be distingushed is muon-induced Cherenkov photons, without suffering from scattering, reflection or absorption. Hit PMTs fromthese photons have strong timing and position constraint because the photons
generated by muon go directly onto the PMTs in the angle of 0 ∼ 50◦ . Energy
loss of muon in water is ∼2MeV/cm, so the Cherenkov angle can decrease to 0
when the sub-GeV muon loss the energy and finally stop in the tank. While for
300MeV/c muon, the Cherenkov angle is still 42◦ and the corresponding hits

77cluster at cosθ = 0.74.

Figure 61: The angle between the direction of δ ray and parrent muon. The
horizontal axis is muon momentum. Red and blue dots are the mean value andmaximum value for δ − muon angle respectively.
The delta rays are intrinsically the electrons which obtain energy from muon
when it is lossing energy by ”ionization” in the water. Those electrons can bea few eV to a few tens MeV, with the angle to the muon track by 40 ∼ 65◦ ,depending on the muon momentum. Most of those electrons are invisible, but
a small part of them can be above Cherenkov threshold and give Cherenkov
photons. The delta ray photons have a wide angle range and they are themost difficulty in gamma searching. For a sub-GeV muon, the delta ray has
a maximum angle of ∼ 60◦ to the muon direction. Since the delta ray gives
Cherenkov light, the photons induced by delta ray can go ∼ 102◦ from themuon direction, while a small part of them can even go to the back half by
scattering.
The dark hits are not real photons, but be due to the dark noise of PMTelectronics. They are randomly distributed in timing and have a frequency of∼3.5kHz, which means for 11129 PMTs in 1.3us gate we could expect a typicalnumber of 51 dark hits.
For the muon-induced photons which are scattered or reflected before detected, they travel for longer distance and thus the hit PMTs will have later
responce. From the cosθ − (T − T OF ) map, it can be found that the scatteredphotons, reflected photons, and dark hits can be removed by a −10ns + T0 <T − T oF < +10 + T0 cut, while T0 is the time offset given by software trigger
and usually represents the timing of maximum hits in 1.3µsec window.
A histogram of cosθ is made below. After −10ns + T0 < T − T oF <+10 + T0 cut, three main component remain: photon induced by gamma, muon
Cherenkov effect, and delta ray. De-excitation gamma goes to random direction
so it is uniform in cosθ. Cherenkov photon of muon clusters at cosθ u 0.74.
The most difficult one to remove is delta ray. In order to distingush the muons78

T-TOF(ns)γ [0.18,0.31,-0.93]6.0MeV/c µ[-0.45,0.04,0.89]300.0MeV/c [-431.1cm,609.4cm,-768.1cm]

1100

1050

1000

950
µ Cherenkovµ δ rayµ reflectionµ scattering
γ hit
dark hit

900

850
-1

-0.8

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

1
cosθ

Figure 62: A typical simulation event of the pair of 300MeV/c muon and 6MeV
γ. The hit PMTs are tagged with the source of the photons. In this event, theangle between muon and γ is 154.7◦ .
with 6MeV gamma, a likelihood built on the hit numbers after −10ns + T0 <T − T oF < +10 + T0 cut will be explained in next section.
In Fig 63, the muon is 300MeV/c. To study further more on sub-GeV muon,
the photon component of the pair of a 200-500MeV/c muon and a 6MeV/c
gamma is shown in Figure 64.
8.2.3

Building Likelihood

In real data, the muon vertex and direction are supposed to be known, as
they can be reconstrued with precise fitter. As explained in previous section,
scattering and reflection photons have later timing responce, they can be removed by a simple ±10ns limit on T-ToF. Dark hits are randomly disturbed intime, so they can also be removed by the −10ns + T0 < T − T oF < +10 + T0
cut. To deal with the remaining three component, gamma hits, delta ray hits
and muon Cherenkov hits, we divide the cosθ into three regions.
• −1 < cosθ ≤ −0.34. Delta ray photons largely decrease in this range but
a few percent of them still remains. This is the best region to search fora extra gamma.
• −0.34 < cosθ ≤ 0.6. Most of the hits in this range are from delta rays.
Muon Cherenkov photons can be ignored and 47% gamma photons are inthis region.
• 0.6 < cosθ ≤ 1. In this range, most of the PMT hits are muon Cherenkov
photons. Part of the delta ray photons also goes into this region. Since
gamma is uniform, 20% gamma photons are in this region.
79

300MeV/c µ and 6MeV γ
102

300MeV/c µ and 6MeV γ
102µ Cherenkov
δ rayreflection
scatteringdark hitγµ Cherenkov
δ rayreflection
scatteringdark hit
γ

10
10

1

1
10-1
-1

-0.8

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

1
cosθ

-1

-0.8

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

1
cosθ

Figure 63: The component of hit PMT numbers for µ events of 300MeV/c.
6MeV γ is generated from the same position to uniform directions. X-axis iscosθ to µ direction. Left/right figure is before/after ±10ns T-Tof constrain.
Cherenkov photons without scattering or reflection will arrive PMT surface
without traveling more distance, this is the reason for the dip in cosθ = 0.8
because the PMT hit is tagged by the earliest photon.

1000

hit component
hit PMT numbers

hit PMT numbers

hit componentµ Cherenkov
dark hit
γ hit
δ ray
scatteringreflection

800

700
600

µ Cherenkov
dark hit
γ hit
δ ray
scatteringreflection

500
400

600

300
400
200
200

0
200

100
250

300

350

400

0
200

450
500
µ momentum(MeV/c)

250

300

350

400

450
500
µ momentum(MeV/c)

Figure 64: The component of hit PMT numbers for µ events in the momtemum
range of 200-500MeV/c. 6MeV γ is generated from the same position to uniform
directions. Left/right figure is before/after ±10ns T-Tof hit timing constraint.
For the case of multi photon hit, the hit is tagged by the earliest photon.

80Figure 65: Comparison of PMT hit numbers of random 300MeV/c µ events
w/o 6MeV γ. Nhits1 , Nhits2 , Nhits3 are the hit numbers in −1 < cosθ ≤ −0.34,−0.34 < cosθ ≤ 0.6, 0.6 < cosθ ≤ 1, divided by the angle between µ and γ
respectively.
Figure 65 shows the comparison of the PMT hit numbers in the three regions in case of a 300MeV/c muon with an extra 6MeV gamma(signal) or no
gamma(background). To build a likelihood, we can simply multiply the hit
numbers in region 1 and 2, and Figure 66 shows the likelihood√response. By
cutting the likelihood at the maximum significance Sig = Ns / Ns + Nb , weobtain the cut efficiency of gamma case(Ef fs ) and no gamma case(Ef fb ) ofEf fs = 13% and Ef fb = 89%.
Defining Ntotal as the total number of stopping muon(including both gamma
0case and no gamma case) before likelihood cut, and Ntotal
as the total number
0of muons after likelihood cut, Ntotal and Ntotal are supposed to be known inanalysis of real data. By knowing Ef fs and Ef fb from simulation, the Br(γ) =Ns /(Ns + Nb ) can be calculated by obtaining Ns and Nb from Equation 70:

Ns + Nb = Ntotal(70)0
Ns × Ef fs + Nb × Ef fb = Ntotal
When applying this simple method onto real data, an expected problem is,the de-excitation gamma can be a mix of 6.18MeV, 6.32MeV, 9.93MeV gamma,
and even secondary gamma that the knock-out nucleon hit on another Oxygen
nucleus and knock out a nucleon again. This means the likelihood must be
modeled for each case. Fortunately, a pure de-excitation gamma sample canbe made by coincident search of a gamma and a delayed electron in a typicaltime difference of 2.2us, which refers to the case 2 in section 8.1 that the muon
product from CCQE is invisible.

8.3

Direct measurement of de-excitation gamma in SuperKAs mentioned before, in the real case, the de-excitation γ can be quite complicated then a ideal 6MeV γ that it should be a mix of 6.18MeV, 6.32MeV,
9.93MeV gamma, or even secondary gamma. Therefore, if a de-excitation γ
sample can be made from case 2 in section 8.1 by searching for the pair of aprompt γ and a decay electron, then the sample can be useful in the future
study of finding the de-excitation γ inside a muon event.

81

Signal
Background25

20
U/O-flow (S,B): (0.0, 0.0)% / (0.0, 0.0)%(1/N) dN / dx

TMVA response for classifier: Likelihood

15

10

5

0

0

0.2

0.4

0.6

0.8

1

Likelihood response

Figure 66: Likelihood built by multiplying Nhits1 and Nhits2 in Figure 65. Cut
Efficiency for signal(with γ) and background(without γ) is 13% and 89% respectively
when cutting at L = 0.33 with maximum significance of Sig =√Ns / Ns + Nb .
Considering a typical ∼6MeV γ ray and the energy resolution, it should bereconstructed in the range of 5 ∼ 7M eV . As mentioned before, though a bitof ∼6MeV γ may be tagged as SLE only and thus −0.5 ∼ +1.0µsec gate issaved instead of −5 ∼ +35µsec gate, but due to the large amount of low energy
background tagged with SLE only(they are mainly the radioactive impurity
from the tank structure or surrounding rock in 3 ∼ 5M eV ), it is not efficient tosearch the de-excitation γ by trigger of SLE only.
To start from LE trigger, since the expected energy is in the same range
with solar neutrino(4 ∼ 20M eV ), the search can be started from the spallation cut step of solar reduction process. The pre reduction and spallation cut
are common with solar analysis, which means not only spallation cut but also
22.5kton fiducial volume cut, loose external gamma cut(def f < 400cm), loose
reconstruction quality cut(gv2 − gd2 > 0.1) have been applied. 16 N cut is notincluded because 16 N goes beta decay and does not give delay signals.
By requiring a delay signal in 0.5 ∼ 35µsec window after a LE trigger up
to 50MeV, the energy and time difference distribution is shown in Fig 68 and
Fig 69. As can be seen from energy distribution of delay signal, it shows a good
agreement with Michel specturm above 15MeV. However, there are still large
amount of low energy background exsiting in delay signal. The most of them
are from radioactive impurities which accidently detected in the 35µsec gateof the prompt LE trigger. Those background form the flat component in timedifference distribution. The accidental pair can be efficiently removed by set acontrain on the energy of both prompt and delay signal by Eprompt < 10M eV
and Edelay > 15M eV . Though some nuclear reaction in CCQE may give γ
ray higher than 10MeV, the dominating one should be ∼ 6M eV and hard to82

Figure 67: The reconstructed energy of the prompt event.

Figure 68: The reconstructed energy of the delay signals. Above 15MeV, it
shows a good agreement with Michel spetrum.

83

Figure 69: The decay time fit of the pre-selected pairs of a prompt LE trigger
and a delay signal. The flat component is due to the accidental backgrounds.
extend above 10MeV. So here both Eprompt < 10M eV and Edelay > 15M eV
are required, and can be seen from the time difference distribution after the cut,the flat componenet has been removed and the decay time fit is 1.7 ± 1.8µsec.
The decay time is less than 2.2µsec, but has a good agreement with the decay
time fitted from visible muon sample.
Though the remaining pairs of prompt and decay signal are muons anddecay electorns, most of the prompt events are not yet the de-excitation γ weare looking for. This is because when the weak muon whose energy is near
Cherenkov threshold, is possible to be reconstructed as an electron or a γ,
the weaker it is, more likely this happens. So up to here, though the energy
distribution of delay events has good agreement with Michel spectrum and decay
time fit is also consistent with visible muon sample, but it is still a problem thatthe prompt signal may be an invisible muon and a de-excitation γ as we want,
but there are also just weak muons near to Cherenkov threshold tagged with
LE trigger and mis-reconstructed.
Usually even the muon is close to Cherenkov threshold, they are hard to berecontructed below 10MeV(this is also the reason for Eprompt < 10M eV selection). There is another useful to distingush γ ray and weak muon: Cherenkov
angle. This is because when the weak muon give Cherenkov light and lose energy
in water, it will soon go below the Cherenkov threshold and become invisible.
Thus the Cherenkov angle fit of a weak muon will be rather smaller than 42◦
electron. As mentioned in relic section, for γ ray, the Cherenkov light is fromscattered electrons, so it can be considered as multi-rings. While the scatteredelectrons will suffer from multi-scattering which causes the ring to be more dirty
and fuzzy. As a result, the Cherenkov angle of a ∼ 6M eV γ ray will be larger,
usually above 50◦ .
The Cherenkov angle fit method is the same used in relic section, and thefitted angle for remaining prompt signal is shown below. In relic analysis, the84

Figure 70: The decay time fit of the remaining events after Eprompt < 10M eV
and Edelay > 15M eV constrain.

Figure 71: The distribution of opening angles obtained from all 3-hit PMT
combination of a 170M eV /cµ event. The peak is around 25◦ .

85

Figure 72: The distribution of opening angles obtained from all 3-hit PMT
combination of a typical 6MeV γ event. The peak is above 50◦ .

Figure 73: The Cherenkov angle fit result of the remaining events after
Eprompt < 10M eV and Edelay > 15M eV constrain.

86electron region is defined as 38◦ ∼ 50◦ and > 50◦ part is considered to be γ ray.
Here we make the same defination that the < 38◦ part should be weak muons
which are close to Cherenkov threshold, 38 ∼ 50 is a mix region of muon and γ,
and > 50◦ is the γ region we are looking for. Finally after Cherenkovangle >50◦ cut, for 10 years SKIV data, there are 58 events found for the search ofde-excitation γ.

879Conclusions and Future

9.1

Neutrino search associated with GW170817

This thesis made a coincidence search for neutrino signals with the gravitational wave, GW170817, produced by a binary neutron star merger in NGC4993,
in the Super-Kamiokande detector in an energy range from 3.5 MeV to ∼100 PeV.
The analysis was performed within a time window of ±500 s of GW170817 and
14 days after the neutron star merger.
In the high-energy data sample, three neutrino interaction categories are
considered: FC, PC and UPMU. No neutrino candidate was found in the ±500s window. The numbers of candidates in a 14-day time window in the entire
sky, as well as in a limited spatial region around NGC4993, are consistent withthe expectation.
Low-energy neutrino events were also examined using the SRN and the solar
neutrino data samples in the same window. No neutrino candidate was found in
the SRN and solar neutrino data samples in the ±500-s window. Two candidates
were found in the SRN data sample in the 14-day search window, which isconsistent with the estimated background rate.
Considering the observation of no significant neutrino signal associated with
the GW170817 in SK, we calculated the neutrino fluence limits. The obtained
results give the most stringent limits for neutrino emission in the energy region
below 100 GeV, especially that the thermal neutrino from binary neutron starmerger is expected to be 10 ∼ 30MeV.
The binary neutron star mergers are expected to occur more than once a
year [64] at the distance of ≥ 100M pc, while a single merger in such distanceis difficult to make detectable signal in SK. Comparing to the expected rate ofsupernova, the rate of binary neutron star merger is two order smaller, and thus
”the diffuse neutron star merger neutrino background” will be hidden under the
diffuse supernova neutrino background. However, with the precise detection for
the timing of gravitational wave, the timing constrain of neutrino searching in≈ 1s from merger time could substantially reduce the contamination from othersources of neutrinos. An idea of how to stack the results of neutrino search frommultiple mergers has been suggested in [8].
On the other hand, a next generation underground water Cherenkov detector
with 260kton volume has been planned [65], which will bring larger target andhigher sensitivity for the neutrino signal from binary neutron star merger.

88

9.2

De-excitation gamma from atmospheric neutrino CCQE
interaction

SK-Gd project can efficiently reduce background and increse ν̄e sensitivity
for SRN siganl search. However, atomospheric νµ CCQE interaction which gives
invisible muon and decay electron, still remains as an unreducable component
in DSNB spectrum. In order to estimate CCQE background component in observed SRN spectrum, direct measurement for the branching ratio of de-excation
gamma is needed. This thesis introduced the method to directly measure thebranching ratio of de-excitation gamma based on the CCQE event simulation of
SK detector. The direct measurement can not only benefits SRN analysis, butalso reduces the uncertainty in long baseline neutrino ossillation experiment.
As the first step, this thesis start from simulation and an analysis methodof searching for 6MeV gamma inside a muon event has been introduced.
As the second step, this thesis also reported the current analysis result of
CCQE de-excitation gamma with SK-IV data. The result indicates that the deexcitation gamma sample includes the muon events which are close to Cherenkov
threshold and tagged as low energy events. However the contamination can be
efficiently removed by Cherenkov angle cut and thus a pure sample of CCQE
de-excitation gamma could be made. The analysis reported by this thesis is up
to here.
For the next step in the future, by using the flux of atmospheric neutrino
and the theoritical expectation of CCQE cross section, the expected number of
CCQE de-excitation gamma events could be estimated and compared with thecurrent sample.
For the final step of this study, the de-excitation gamma sample could be
used for likelihood building and direct searching for de-excitation gamma inside
the Fully-Contained muon events of SK real data.

89

References[1]

B.P.Abbott et al. The Astrophysical Journal Letters, 848, L13 (2017)[2]

B.P.Abbott et al. Physical Review Letters, 119, 161101 (2017)[3]

B.P.Abbott et al. The Astrophysical Journal Letters, 848, L12 (2017)[4]

D.Coulter et al. Science, 358, 1556 (2017)[5]

A.Albert et al. The Astrophysical Journal Letters, 850, L35 (2017)[6]

K.Abe et al. The Astrophysical Journal Letters, 830, L11 (2016)[7]

Y.Sekiguchi et al. Physical Review Letters, 107, 051102 (2011)[8]

K.Kyutoku and K.Kashiyama, Physical Review D, 97, 103001 (2018)[9]

K.Abe et al. The Astrophysical Journal, 850, 166 (2017)[10] E.Thrane et al. The Astrophysical Journal, 704, 503 (2009)[11] M.Swanson et al. The Astrophysical Journal, 652, 206 (2006)[12] T.Gaisser et al. Physics Reports, 258, 173 (1995)[13] Y.Hayato et al. Acta Physica Polonica B, B40, 2477 (2009)[14] M.Nakahata et al. Nuclear Instruments and Methods in Physics Research
Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, 421, 113 (1999)[15] A.Ankowski et al. Physical Review Letters, 108, 052505 (2012)[16] K.Abe et al. Physical Review D, 90, 072012 (2014)[17] M.Leuschner et al. Physical Review C, 49, 955 (1994)[18] F.Ajzenberg-Selove et al. Nuclear Physics A, 523, 1 (1991)[19] K.Kobayashi et al. ArXiv Nuclear Experiment e-prints, nucl-ex/0604006(2006)[20] S.Fukuda et al. Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment,
501, 418 (2003)[21] http://sales.hamamatsu.com/assets/pdf/parts R/R3600-06.pdf[22] K.Kobayashi et al. ArXiv Nuclear Experiment e-prints, nucl-ex/0604006(2006)[23] H.Ikeda et al. Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment,
320, 310 (1992)[24] T.Tanimori et al. IEEE Transactions on Nuclear Science, 36, 497 (1989)90

[25] T.K.Ohsuka et al. KEK Report 85-10 (1985)[26] H.Nishino et al. Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment,
610, 710 (2009)[27] Y.ZHANG, PhD Thesis, Tsinghua University (2015)
http://www-sk.icrr.u-tokyo.ac.jp/sk/ pdf/articles/2015/thesis yangzhang 20150615.pdf[28] S.Yamada et al. IEEE Transactions on Nuclear Science, 57, 428 (2010)[29] J.Beacom and M.Vagins, Physical Review Letters, 93, 171101 (2004)[30] K.Abe et al. Nuclear Instruments and Methods in Physics Research Section
A: Accelerators, Spectrometers, Detectors and Associated Equipment, 737,
253 (2004)[31] Y.Koshio, PhD Thesis, University of Tokyo (1998)
http://www-sk.icrr.u-tokyo.ac.jp/sk/ pdf/articles/koshio-dthesis.pdf[32] Y.Nakano, PhD Thesis, University of Tokyo (2016)
http://www-sk.icrr.u-tokyo.ac.jp/sk/ pdf/articles/2016/doc thesis naknao.pdf[33] M.Ikeda, PhD Thesis, University of Tokyo (2010)
http://www-sk.icrr.u-tokyo.ac.jp/sk/ pdf/articles/ikeda-dron.pdf[34] T.Yokozawa, PhD Thesis, University of Tokyo (2013)
http://www-sk.icrr.u-tokyo.ac.jp/sk/ pdf/articles/yokozawa thesis.pdf[35] M.Smy. Proceedings of the 30th 30th International Cosmic Ray Conference,
5, 1279 (2007)[36] E.Blaufuss et al., Nuclear Instruments and Methods in Physics Research
Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, 458, 638 (2001)[37] M.Nakahata et al., Nuclear Instruments and Methods in Physics Research
Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, 421, 113 (1999)[38] J.Hosaka et al., Physical Review D, 73, 112001 (2006)[39] ”GEANT - Detector Description and Simulation Tool”, CERN Program
Library (1993)[40] B.P.Abbott et al., Physical Review Letters, 116, 061102 (2016)[41] B.D.Metzger and E. Berger, The Astrophysical Journal, 746, 48 (2012)[42] J.J.Ruan et al., The Astrophysical Journal Letters, 853, L4 (2018)[43] K.Hotokezaka et al., Physical Review D, 83, 124008 (2011)[44] K.Hirata et al., Physical Review Letters, 58, 1490 (1987)91

[45] S.Richers et al., The Astrophysical Journal, 813, 38 (2015)[46] A.Strumia and F.Vissani, Physics Letters B, 564, 42 (2003)[47] S.Fujibayashi et al., The Astrophysical Journal, 846, 114 (2017)[48] S.Fukuda et al., Physics Letters B, 539, 179 (2002)[49] Y.Fukuda et al., Physical Review Letters, 81, 1562 (1998)[50] K.Abe et al., Physical Review Letters, 107, 041801 (2011)[51] J.N.Bahcall et al., Reviews of Modern Physics, 54, 767 (1982)[52] J.N.Bahcall et al., The Astrophysical Journal, 621, 85 (2005)[53] R.M.Bionta et al., Physical Review Letters, 58, 1494 (1987)[54] E.Cappellaro et al., Astronomy & Astrophysics, 351, 459 (1999)[55] G.S.Bisnovatyi-kogan and Z.F.Seidov, Annals of the New York Academyof Sciences, 422, 319 (1984)[56] T.Totani and K.Sato, Astroparticle Physics, 3 ,367 (1995)[57] T.Totani, K.Sato and Y.Yoshii, The Astrophysical Journal, 460, 303 (1996)[58] R.A.Malaney, Astroparticle Physics, 7, 125 (1997)[59] D.H.Hartmann and S.E.Woosley, Astroparticle Physics, 7, 137 (1997)[60] M.Kaplinghat, G.Steigman, and T.P.Walker, Physical Review D, 62,043001 (2000)[61] S.Anto, K.Sato, and T.Totani, Astroparticle Physics, 18, 307 (2003)[62] K.Bays et al., Physical Review D, 85, 052007 (2012)[63] H.Zhang, PhD Thesis, University of Tsinghua (2012)
http://www-sk.icrr.u-tokyo.ac.jp/sk/ pdf/articles/zhanghb thesis.pdf[64] M.Dominik, The Astrophysical Journal, 806, 263 (2015)[65] K.Abe et al., arXiv:1109.3262[66] Y.Ashie, Physical Review D, 71, 112005 (2005)

92