学位論文

Search for proton decay
into a muon and a neutral kaonin the Super-Kamiokande water Cherenkov detector（スーパーカミオカンデ水チェレンコフ検出器
によるミューオンと中性 K 中間子への陽子崩壊探索）2023 年 3 月

松本遼

Abstract
A grand uniﬁed theory (GUT) is an extension of the standard model (SM) of particle physics that
uniﬁes leptons, baryons and gauge bosons. The GUT is considered a strong candidate for new
physics beyond the SM as it provides an explanation to some fundamental questions in our universe
with a natural extension of the SM although it has not been conﬁrmed by any experiments. For
this reason, veriﬁcation of the GUT is one of the most important subjects in particle physics today.
Whereas protons are considered stable particles in the SM, GUT models predict proton decay
through interactions in which the baryon number is not conserved. Therefore, the measurement ofproton decay is direct evidence of the GUT.
In this thesis, a search for proton decay via p → µ+ +K 0 is described with 0.37 Mton·years of datacollected between 1996 and 2018 from the Super-Kamiokande detector. The Super-Kamiokande
detector is a large water Cherenkov detector located 1000 m underground in the Kamioka mine in
Hida City, Gifu Prefecture. It contains 50 kton of ultrapure water in a cylindrical tank with 39.3 m
in diameter and 41.4 m in height. The Cherenkov light emitted from charged particles in the water isdetected by approximately 11000 optical sensors (photomultiplier tubes) installed on the wall. The
analysis method was updated for the data taken between 2008 and 2018 (referred to as SK-IV) with
upgraded electronics that allows for the detection of neutron capture on hydrogen. In addition toneutron tagging, the analysis was performed with an improved event reconstruction algorithm, andevent selection criteria optimized to select proton decay signals against the atmospheric neutrino
background events. The eﬃciency of the proton decay signals and the expected background eventswere estimated based on the simulation. As K 0 exists as superpositions of KS0 and KL0 , the selection
criteria were deﬁned separately for KS0 and KL0 channels.
Observations of atmospheric neutrinos are themselves important physics targets of the SuperKamiokande detector to measure the transition of neutrino ﬂavor (neutrino oscillation) and reveal
their properties including as yet unknown neutrino mass ordering and CP phase. This thesis also
describes a new neutrino event classiﬁcation method using machine learning. The performance of
the identiﬁcation of neutrino ﬂavor was improved with the machine learning technique by accounting
for the correlations of the observable kinematic parameters.
As a result of the proton decay search, no signiﬁcant event excess has been observed over theiexpected background rates. According to the observation, the lower limit of the lifetime of proton
decay to a muon and a neutral kaon was set as 4.5 × 1033 years from the SK-IV data. Furthermore,this result was statistically combined with the result of the previous analysis [1] using the data
during 1996 and 2008, and the lower limit was set as 3.6 × 1033 years. This is the most stringent
limit in the world, updating the lifetime limit from the previous analysis more than factor two.
This study achieved giving additional constraints to the GUT models with the lower limit on the
lifetime.

Contents
Abstracti

1 Physics background1

1.1

Standard model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1

1.2

Uniﬁcation of gauge interactions . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3

1.3

Grand uniﬁed theories . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3

1.3.1

Minimal SU(5) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5

1.3.2

Other GUT models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6

Proton decay via p → µ+ + K 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7

1.4

2 Motivation

9

3 Super-Kamiokande detector

11

3.1

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

11

3.2

Detector structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

13

3.3

Cherenkov light . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

14

3.4

Photomultiplier tube . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

16

3.5

Water and air puriﬁcation system . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

20

3.6

Data acquisition system . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

21

3.7

Scattering and absorption in water . . . . . . . . . . . . . . . . . . . . . . . . . . . .

22

3.8

Energy scale calibration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

24

4 Simulation

25

4.1

Proton decay . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

25

4.2

Final state interaction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

26

4.3

Detector simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

27

5 Atmospheric neutrino

30

5.1

Neutrino ﬂux . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

30

5.2

Neutrino oscillation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

31

iii

5.3

Neutrino interaction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

33

5.4

Event classiﬁcation of atmospheric neutrino . . . . . . . . . . . . . . . . . . . . . . .

37

6 Data reduction

42

6.1

First reduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

44

6.2

Second reduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

44

6.3

Third reduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

45

6.4

Fourth reduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

48

6.5

Fifth reduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

49

6.6

Final fully contained dataset

51. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7 Event reconstruction
7.1

53

Reconstruction algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

53

7.1.1

Likelihood function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

54

7.1.2

Vertex pre-ﬁtter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

55

7.1.3

Hit clustering . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

57

7.1.4

Single-ring ﬁtter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

57

7.1.5

Multi-ring ﬁtter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

57

7.2

Performance of ﬁTQun . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

60

7.3

Displaced vertex ﬁter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

63

7.4

Performance of displaced vertex ﬁtter . . . . . . . . . . . . . . . . . . . . . . . . . .

65

7.5

Neutron tagging . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

66

8 Event selection

67

8.1

Overview of selection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

67

8.2

Selection for p → µ+ + KS0 , KS0 → 2π 0 . . . . . . . . . . . . . . . . . . . . . . . . . .

69

8.3

Selection for p → µ. . . . . . . . . . . . . . . . . . . . . . .

72

8.4

Selections for p → µ+ + KL0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

74

8.5

Selection summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

79++ KS0 , KS0 → π + + π −9 Lifetime limit

82

10 Systematic uncertainties

84

10.1 Correlated decay . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

84

10.2 Fermi momentum . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

86

10.3 Pion interaction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

86

10.4 Kaon interaction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

87

10.5 Neutrino ﬂux and cross section . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

90

10.6 Fiducial volume . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

90

10.7 Energy scale . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

91

10.8 Detector non-uniformity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

92

10.9 Ring counting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

92

10.10Particle type identiﬁcation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

93

10.11Michel electron tagging. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

94

10.12Vertex separation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

95

10.13Neutron tagging eﬃciency . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

96

11 Conclusion

98

Acknowledgments

99

A p → µ+ + K 0 search in SK-I to SK-IIIi

BibliographyiiChapter 1

Physics background
1.1

Standard model

The standard model (SM) is the most fundamental framework of particle physics (Table 1.1). In theSM, matter consists of Dirac particles and they are classiﬁed into two categories, that are quarks
and leptons. These quarks and leptons interact with each other via gauge particles. There are three
types of gauge particles: gluon, photon and weak boson which mediate interactions corresponding
to the fundamental forces in nature, i.e. strong, electromagnetic and weak forces, respectively.
These particles have a pair of antiparticles with opposite electric charges and equal masses. In
addition to these particles, there is a Higgs boson, whose vacuum expectation spontaneously breaks
the symmetry and gives mass to particles.

Quarks
Leptons
Gauge ﬁelds
Scalar ﬁelds

Table 1.1: Elementary particles in the SM.
1st generation
2nd generation
3rd generationup (u), down (d)charm (c), strange (s)top (t), bottom (b)electron (e), e-neutrino (νe ) muon (µ), µ-neutrino (νµ ) tauon (τ ), τ -neutrino (ντ )photon (γ), weak-boson (Z 0 , W ± ), gluon (g)Higgs (H 0 )Quarks and leptons have pairs of diﬀerent charges and form doublets for weak interactions. The
examples are up quark u (charge +2/3 in the unit of elementary charge) and down quark d (charge−1/3) for quarks, and electron e (charge −1) and electron neutrino νe (charge 0) for leptons. They
can change into each other forming doublets via weak bosons. Figure 1.1 shows examples of thecharge-exchange process via weak interactions. The states of quarks and leptons can be divided
by the helicity, which is deﬁned as the projection of the spin onto the direction of momentum,
into those in which the spin direction is the same as the direction of momentum (right-handed) or
opposite (left-handed) as shown in Figure 1.2. The weak interaction works only on the left-handed
component of particles.
Quarks have three degrees of freedom called colors (red, blue and green) for strong interactions.

1

u

νe

d

W+

e

W+

Figure 1.1: Diagram of the charge-exchange process via weak interactions.

Left-handed

Right-handed

Momentum

Momentum

Spin

Spin

Figure 1.2: Schematic diagram of helicity. When spin, the angular momentum intrinsic to theparticle, and momentum are in the same (opposite) direction, it is right-handed (left-handed).
A color becomes stable with its anti-color (e.g. red and antired), and an antiparticle has an anticolor.
Also, a set of two colors becomes an anticolor of a remaining color, so the combinations of redblue-green and red-antired are stable combinations for example. Therefore, unlike leptons, quarks
cannot be observed as a single state and appear in nature as composite particles, called hadrons,
in which quarks are bound by strong interactions. Among the hadrons, those composed of threequarks of antiquarks are called baryons, and those composed of one quark and one antiquark arecalled mesons (Figure 1.3). For example, a proton and neutron consist of three quarks as uud andudd, respectively. A pion and kaon are examples of mesons which are composed of u, d, s quarks and¯ π 0 : uū or dd,¯ π − : dū, K + : us̄, K 0 : ds̄ or K 0 : sd,¯ and K − : sū. Eachtheir antiquarks as π + : ud,
quark (antiquark) has a conserved additive quantum number 1/3 (−1/3) called baryon number B.
Proton is the lightest baryon with B = 1 and therefore stable in the SM due to the baryon number
conservation.
Furthermore, it is known that there is a hierarchy (called generation) of quarks and leptons:
three types of particles with equal electric charge but diﬀerent masses. In total, six types (called
ﬂavors) of quarks and leptons exist.
All particles in the SM were conﬁrmed by the experiments, among those Higgs boson was
discovered in 2012 using the Large Hadron Collider (LHC) accelerator of CERN. The SM has beentested up to the TeV energy scale by various experiments and the observed data are precisely
reproduced by the prediction from the SM. On the other hand, various fundamental problems
remain in the theory such as the reason why positrons (antiparticles of electrons) and protons haveexactly equal charges and the origin of the matter in the universe.

2BaryonMeson

Figure 1.3: Schematic diagram of color combination for baryons (left) and mesons (right). Baryons
consist of red, green and blue, and mesons consist of red and antired (cyan). Other combinations,i.e. blue and antiblue (yellow), and green and antigreen (magenta), are also possible for mesons.
Colorless combinations are stable.

1.2

Uniﬁcation of gauge interactions

The SM contains a gauge symmetry of SU(3)×SU(2)×U(1), and among the gauge interactions theelectromagnetic and weak interactions are uniﬁed as electroweak interactions SU(2)×U(1) at highenergy regions. Since the symmetry of the left-handed and right-handed is broken, the left-handed
fermions and right-handed fermion form SU(2) doublet and U(1) singlet, respectively.


νeL
, eR
eL

(1.1)

The electroweak interaction has SU(2)×U(1) gauge symmetry, which is spontaneously broken toU(1)EM symmetry of the electromagnetic interaction due to the Higgs mechanism with a vacuum
expectation value of the Higgs ﬁeld. On the other hand, the strong interaction is described by
SU(3) group and forms triplets for quarks with diﬀerent colors.
The coupling constants between the gauge ﬁelds and fermions show diﬀerent dependencies on
the energy scale. Electromagnetic interactions tend to be more strongly coupled at higher energies,
whereas strong interactions tend to have weaker coupling constants at higher energies with a small
distance and large momentum transfer. Figure 1.4 shows the energy dependence of the U(1),SU(2) and SU(3) coupling constants. These coupling constants approach the same value around1014−16 GeV. It indicates a uniﬁed model of these gauge interactions beyond the SM at high energies.

1.3

Grand uniﬁed theories

A grand uniﬁed theory (GUT) is an extension of the SM that uniﬁes leptons, baryons and gaugebosons [3]. In the GUT, strong, electromagnetic and weak forces are uniﬁed at high energy, typicallyabove 1014 GeV, which is much higher than the energy scale achieved by the current particle
3

SM
60

α1
α2
α3

50

-1

αi (Q)

40
30
20
10
SOFTSUSY 3.6.2

0
0

5
10
log10(Q/GeV)

15

Figure 1.4: Running couplings of U(1) (red, α1 ), SU(2) (green, α2 ) and SU(3) (blue, α3 ) for theenergy scale Q in the SM using two-loop RG evolution [2].
accelerators. The three forces are described by a single coupling constant in the GUT. Various
GUT models based on diﬀerent gauge groups have been theoretically proposed, while there is noevidence for the GUT given by the experiments so far. and their veriﬁcation is one of the mostimportant subjects in particle physics today.
In the SM, quarks and leptons are not uniﬁed and there are no interactions that convert into eachother, unlike electrons and neutrinos which are uniﬁed and transformed into each other throughthe exchange of weak bosons. Baryon number is conserved as a consequence of the absence of suchinteractions which convert quarks into leptons. Therefore, protons do not decay in the SM as they
are the lightest baryon. In contrast, many of the GUTs predict proton decay into particles including
leptons through interactions in which the baryon number is not conserved. The dominant proton
decay channel depends on the assumptions of the framework and parameters of the theoreticalmodels.
The existence of baryon number violation is one of the Sakharov conditions for the developmentof the universe after the Big Bang to break the symmetry between matter and antimatter and toestablish a matter-dominated universe [4]:1. Baryon number violation
2. C and CP symmetry breaking
3. Interactions out of thermal equilibrium

4

The GUTs accommodate interactions to violate the baryon number in the theory of particle physics
which satisﬁes one of the conditions for explaining the origin of matter.

1.3.1

Minimal SU(5)The minimal SU(5) model [5] is the simplest model of the GUT, which contains SU(3), SU(2) andU(1) gauge groups to unify the strong, weak and electromagnetic interactions.
The particles correspond to the representation of the SU(5) group. The gauge ﬁelds are assigned
to the adjoint representation (24 components), and the upper left 3×3 of the 5×5 matrix corresponds
to the SU(3) gauge and the lower right 2 × 2 corresponds to the SU(2) gauge.


g12
g13
X̄ 1
Ȳ 1
g11 − √2B
30


g21
g22 − √2B
g23
X̄ 2
Ȳ 2


30


3
3


g31
g32
g33 − √2B
X̄
Ȳ
30


3
W
3B


+
1
2
3
√
√
X
X
X
+
W


2
30
3
W
3B
1
2
3
−
√
√
Y
Y
Y
W
− 2 + 30

(1.2)

where g are gluons, and B and W are gauge bosons in the electroweak theory. The X and Y bosons
are not included in the SM and have the charge of 4/3 and 1/3, respectively. It can be explainedthat X and Y bosons have not been observed yet due to the extremely heavy mass compared tothe other gauge bosons. To avoid introducing new particles, quarks and leptons are assigned tobasic and antisymmetric representation. The sum of the number of components (15 = 5 (basic) +10 (antisymmetric)) is consistent with the number of fermions in the SM (15 = 2 (u, d) × 3 (color)−× 2 (left, right) + 2(e−L , eR ) + 1(νL ), where L and R represent the left-handed and right-handed

fermions). The 5∗ representation (complex conjugate representation of 5 representation) is
 c 
drR
 dcgR  dc 
 bR 
 −eL νL

(1.3)

to be the total charge is zero since SU(5) is traceless. Here r, g and b represent the colors, andc represents the charge conjugation of particles (it means antiparticles). Since the charge of thebasic representation is [−1/3, −1/3, −1/3, 1, 0], the charge for antisymmetric representation should
be assigned as

−2/3
 −2/3 −2/3

 2/3
−1/3

−2/3
−2/3
−2/3
2/3
−1/3

−2/3 2/3−2/3 2/3−2/3 2/3
2/3
2
−1/3
1

5


−1/3
−1/3 

−1/3 
.
1 
0

(1.4)

The remaining fermions are assigned to these components with attention to the colors:
0 ucbR −ucgR urL drL0
ucrR ugL dgL 0ubL dbL 

,

0
ecR 
0

(1.5)

where the opposite components are omitted since it is an antisymmetric matrix. Although themass of X and Y bosons are much higher than the mass of protons, the exchange of these bosons
occurs in protons with an extremely small probability. This model predicts proton decay with thedominant mode as p → e+ + π 0 via the exchange of X boson (Figure 1.5) with an extremely long
lifetime.
u

ē
X

d

ū

u

u

Figure 1.5: Feynman diagram of proton decay p → e+ + π 0 via the exchange of X boson. uud
quarks in the left-hand side form a proton and ūu in the right-hand side form a neutral pion (π 0 ).
The expected lifetime of protons is approximately 1030 years in this model assuming the massof X and Y are an order of 1014 GeV/c2 from the uniﬁcation of the coupling constants (Section
1.2). Although this lifetime is longer than the time scale of the universe, which started about 10 8
years ago from the Big Bang, it is still possible to search for the proton decay and test the GUT
by observing a large number of protons in the detector. For example, if the proton lifetime is 10 30
years, 10 protons out of 1031 protons decay on average in one year. This minimal SU(5) model
has already been rejected by the Kamiokande, IMB-3 and Super-Kamiokande experiments in the
searches for p → e+ + π 0 [6][7][8].

1.3.2Other GUT models

Supersymmetric SU(5) is an extension of the simplest SU(5) model. Supersymmetry (SUSY) isa framework beyond the SM based on symmetry between a boson and a fermion, which requires
the existence of partner fermions (bosons) in the bosons (fermions) appearing in the theory. The
theory is motivated by resolving ﬁne-tuning of parameters to explain the Higgs mass that areexperimentally determined in the SM. The contributions from the SUSY particles make the energy
scale of the GUT higher and extend the lifetime of the simple decay via X bosons. As a consequence,
the SU(5) gauge model still survives depending on the parameters in SUSY, such as the couplingof the SUSY particles with the SM particles and the ﬂavor mixing of the SUSY bosons. Minimal
6

SUSY SU(5) model [9] predicts proton decay with the dominant mode as p → ν̄ + K + and thelifetime ≤ 1.1 × 1035 years when the SUSY scalar masses are required to be ≤ 30 TeV, so that they
are within reach of collider experiments.
In addition to SU(5), SO(10) is also a group including SU(3)×SU(2)×U(1). Before breaking intoSU(3)×SU(2)×U(1) of the SM, SO(10) breaks into intermediate symmetries such as SU(5)×U(1)and SU(4)×SU(2)×SU(2). In the former case, the right-handed lepton and quark assignments
are ﬂipped compared to the standard SU(5) and it is called the ﬂipped SU(5) model [10]. As an
example of the latter case, minimal non-supersymmetric SO(10) model [11] predicts the proton
lifetime as O(1035 ) years and the dominant decay channels are p → e+ + π 0 (branching ratio
BR=47%), p → µ+ + π 0 (branching ratio BR=1.0%), p → ν̄ + π + (branching ratio BR=48%) andp → µ+ + K 0 (branching ratio BR=3.6%).
Supersymmetric GUT models often predict p → ν̄ + K + and p → µ+ + K 0 as the dominantchannel due to the contribution of the color triplet Higgs which are predicted in typical SUSY
models as an extension of the Higgs mechanism to introduce symmetry breaking at diﬀerent energy
scales [12]. Some models [12] predict the partial lifetime of the p → µ+ + K 0 channel can becomparable to p → ν̄ + K + or shorter and in a range that slightly exceeds the current experimental
limit of 1.6 × 1033 years [1]. These models have a prediction of the proton lifetime as O(1034 ) years
and the fraction is accounted for by p → µ+ + K 0 channel. p → e+ + π 0 , p → µ+ + π 0 , p → ν̄ + K +and p → e+ + K 0 appear as the other dominant decay channels.
Also, there is a R-symmetric ﬂipped SU(5) model [13] which predicts protons do not decay
via p → ν̄ + K + because the color triplet Higgs does not couple to the left-handed neutrino, andp → µ+ + K 0 becomes the dominant decay channel. This model predicts the proton lifetime asτ (p → e+ + π 0 ) ∼ 1035 years, τ (p → µ+ + π 0 ) ∼ 1035 years, τ (p → e+ + K 0 ) ∼ 1037 years andτ (p → µ+ + K 0 ) ∼ 1034 years if the color triplet Higgs mass is ∼ 1012 GeV and tan β = 3, where β
is a Higgs vacuum expectation ratio.
These models generally prefer shorter lifetimes which are accessible by the search with theSuper-Kamiokande detector, or the next-generation projects such as Hyper-Kamiokande. Several
decay channels of proton decay have been searched in Super-Kamiokande as shown in Table 1.2.

1.4

Proton decay via p → µ+ + K 0

The study described in this thesis focuses on a search for the p → µ+ + K 0 decay mode in the
Super-Kamiokande (SK) detector. This two-body decay generates a muon and a neutral kaon withmonochromatic momentum of 326.5 MeV/c. The K 0 is a composite state of KS0 and KL0 . KS0 decays
into π + π − (69.2%) and 2π 0 (30.7%) with a lifetime of 90 ps while KL0 decays into π ± e∓ ν (40.6%),
π ± µ∓ ν (27.0%), 3π 0 (19.5%) and π + π − π 0 (12.5%) with the lifetime of 51 ns. Proton decay searches

7Table 1.2: List of proton decay channels which have been searched in Super-Kamiokande.
ChannelLifetime limit (years) Referencep → e+ + π 0
2.4 × 1034
[8]+0p→µ +π1.6 × 1034[8]
p → ν̄ + K +8.2 × 1033
[14]
p → e+ + e+ + e−
3.4 × 1034
[15]
+
+
−
p→µ +e +e
2.3 × 1034
[15]
p → µ− + e+ + e+
1.9 × 1034
[15]
p → e+ + µ+ + µ−
9.2 × 1033
[15]
1.1 × 1034
[15]
p → e− + µ+ + µ+
p → µ+ + µ+ + µ−
1.0 × 1034[15]
p → ν̄ + π +3.9 × 1032
[16]
p → µ+ + K 0
1.6 × 1033
[1]
p → e+ + η
1.0 × 1034
[17]
p → µ+ + η
4.7 × 1033
[17]+0p→e +ρ
7.2 × 1032
[17]
p → µ + + ρ0
5.7 × 1032
[17]
p → e+ + ω
1.6 × 1033
[17]
p → µ+ + ω
2.8 × 1033
[17]
p → e+ + X
7.9 × 1032
[18]
p → µ+ + X
4.1 × 1032
[18]
p → e+ + ν + ν
1.7 × 1032
[19]
p → µ+ + ν + ν
2.2 × 1032
[19]
for the p → µ+ + K 0 mode have been performed in SK using the data collected in the SK-I, SK-II,and SK-III phases (0.17 Mton·years of exposure 1 from 1996 to 2008) and no signiﬁcant signal ofp → µ+ + K 0 was observed [20][1]. A lower limit on the partial proton lifetime p → µ+ + K 0 of
1.6 × 1033 years was set at the 90% conﬁdence level (C.L.). In this thesis, the search was extended
with the data collected in the SK-IV phase from 2008 to 2018 which corresponds to 0.20 Mton·yearsof exposure with the improved analysis method.

1 In proton decay search, megaton·years, a product of the detector mass and observed time, is used as the unit of
the observed quantity because the sensitivity depends on the number of observed protons and observed time.

8Chapter 2

Motivation
As described in Chapter 1, the grand uniﬁed theory (GUT) is a framework beyond the standard
model (SM) of particle physics. The GUT provides explanations for some of the fundamental
questions as yet unresolved in the SM and therefore the veriﬁcation of the GUT is one of the mostimportant subjects in particle physics. The energy scale of the GUT is predicted as the region
around 1014−16 GeV, while this energy scale is too high to test by the current accelerators. On theother hand, the GUT predicts proton decay which is forbidden in the SM due to the baryon number
conservation. Therefore, observation of proton decay provides direct evidence for the GUT although
it has not yet been observed. The Super-Kamiokande (SK) experiment has been exploring various
proton decay channels. The SK experiment limits the lifetime of the proton decay via p → e+ + π 0
as 2.4 × 1034 years [8] which is predicted as a major decay channel in the nominal SU(5) model. The
proton decay via p → ν̄ + K + is also predicted as a major decay channel in SUSY GUT models,
and the lifetime of this channel is limited by SK as 8.2 × 1033 years [14]. In this thesis, the proton
decay via p → µ+ + K 0 has been searched for the data collected from the SK detector. While some
SUSY GUT models predict p → µ+ + K 0 as another major channel, the search of this channel has
not been updated in SK since 2012 although approximately twice of the data has been accumulatedwith respect to the previous publication [1]. Furthermore, the analysis method can be improved
with a new reconstruction algorithm, neutron tagging (Chapter 7) and optimized selection criteria(Chapter 8). These improvements make it possible to expand the lifetime limit by a factor morethan two from the previous result of 1.6 × 1033 years on the lifetime of p → µ+ + K 0 [1]. For thesereasons, there is a possibility to make a discovery of proton decay in this channel and this is the
motivation of this thesis. Since this decay channel has K 0 , the simulation of the kaon interaction
is important to search the events. The simulation was updated for this study by referring to theexperimental data of K-N scattering (Chapter 4 and 10). Combining the results in the previous
search [1], proton decay via p → µ+ + K 0 have been searched with the largest exposure in SK(Chapter 9). These studies are reported in [21].
In addition to the proton decay search, this thesis describes a new classiﬁcation method for
9atmospheric neutrino interactions using machine learning (Section 5.4). Atmospheric neutrinos are
the dominant background in the proton decay search. In addition, the study of neutrino oscillations
is also being conducted using the data from atmospheric neutrino observations. The new method
improves the classiﬁcation and could be used for future studies of the neutrino oscillation parameterssuch as the CP phase and mass ordering. These studies are reported in [22].

10Chapter 3

Super-Kamiokande detector
3.1

Overview

Super-Kamiokande (SK) is a large water Cherenkov detector that contains 50 kton of ultrapure
water in a cylindrical tank (Figure 3.1). It is designed to search for proton decay and reveal theproperties of neutrinos. The Super-Kamiokande detector is located 1000 m underground (2700
m.w.e.) in Gifu Prefecture, Japan. The detector is located underground to reduce the cosmic-ray
muons to 1/100000 of the ground surface by the rock overburden.

Figure 3.1: Schematic view of the SK detector [23].
The detector has two regions, the inner detector and the outer detector surrounding it, andboth regions have photosensors. The inner region works as a target of proton decay search andneutrino interactions while the outer detector is necessary to reject the backgrounds. When protons
in the water decay or neutrinos enter the SK detector and interact with electrons or nuclei inthe water tank, charged particles are produced. As explained in Section 3.3, Cherenkov light is11

emitted when charged particles travel in water with a velocity larger than the propagation of light.
These Cherenkov photons are detected by photomultiplier tubes (PMTs), which are photo sensors
installed on the wall of the water tank, and the particles are observed if the Cherenkov ring is
recognized from the PMT signals. Figure 3.2 shows examples of Cherenkov rings by an electron
and muon. Electrons and gammas generate showering rings called e-like rings because they produce
electromagnetic showers and generate diﬀuse Cherenkov rings. In contrast, muons and charged pions
do not produce showers and generate Cherenkov rings with clear edges called µ-like rings. Using
the event reconstruction algorithm, we estimate basic event properties such as vertex (generated
position), direction, particle type and momentum based on the PMT hit information.

Figure 3.2: Examples of the e-like ring by an electron (left) and µ-like ring by a muon (right) inthe simulation. Points correspond to the PMTs on the inside wall, and the colors correspond tothe observed electric charge of PMTs.
The data used in this analysis are divided into four periods as shown in Table 3.1 by the
conﬁguration of the detector. The ﬁrst period has started as SK-I in 1996 with 11146 PMTs.
During maintenance after SK-I in 2001, one of the PMTs ruptured while being reﬁlled with water,which triggered a chain reaction of shock waves and destroyed more than half of the PMTs. The
SK-II period has started in 2002 with half of the PMTs for SK-I. The locations of the remaining
PMTs are arranged to be uniform in the detector. Since SK-II, all PMTs are protected by cases
to prevent the chain reactions of PMT destruction. After the production of PMTs, SK-III has
started in 2006 with almost the same number of PMTs as SK-I. After the upgrade of the electronics
from ATM to QBEE [24], SK-IV period has started. By the upgrade of the electronics, neutron
tagging was realized even with pure water, as described in Section 7.5. The details of QBEE will
be explained in Section 3.6. After 2018 (the end of SK-IV), the SK detector was further upgradedwith gadolinium loading in water to enhance the detector. In this study, data collected in SK-IV(3244.39 live days) were newly analyzed and the result was combined with the previous result using
SK-I to SK-III (2805.9 live days in total) [1]. The details of the SK detector are described below.

12PhaseSK-I
SK-II
SK-III
SK-IV

3.2

Table 3.1: Detector conﬁguration for each period.
PeriodInner PMTs Outer PMTs PMT cover
1996–2001
11146
1885
no
1885
2002–2005
5182
no2006–2008
11129
1885
yes
2008–2018
11129
1885
yes

Electronics
ATM
ATM
ATM
QBEE

Detector structure

The detector is a cylindrical water tank with a 39.4 m diameter and 41.4 m height ﬁlled with 50 ktonof ultra-pure water. The water tank is divided into two concentric volumes that are opticallyseparated (Figure 3.3). The inner detector (ID) is a cylindrical volume with a 33.8 m diameter

Figure 3.3: A photo from the bottom during the tank renovation in 2018 [25]. The upper part of
the tank is the inner region and the lower part is the outer region, bordered by the center of thephoto.
and 36.2 m height. It contains 32 kton of water. There are 11129 inward-facing PMTs with 50 cm
diameter on the wall to detect the Cherenkov light from charged particles. The ID is surrounded
by a 2 m thickness layer called the outer detector (OD). The 2 m thick layer of water attenuates
γ and neutron radiation from the surrounding rock and suppresses the background contamination.
The OD is composed of 1885 outward-facing PMTs with 20 cm diameter. To increase the lightcollection eﬃciency, the walls of the OD are covered with a reﬂective sheet made of a material
called Tyvek and wavelength shifting plates which convert ultraviolet light to the wavelength in
the eﬀective sensitivity region of the PMTs (Figure 3.4). The PMTs in the OD are mainly used to
identify incoming cosmic ray muons and particles exiting the ID. For example, if charged particles
are generated in the ID due to proton decay or neutrino reactions and lose all kinetic energy inside

13

Figure 3.4: A photo from the bottom of the outer region. The PMTs with 20 cm diameter are
mounted with wavelength shifting plates on the wall on the side of the ID, and the OD is surrounded
by a white Tyvek sheet [25].
it, Cherenkov light is emitted only in the ID and there is no response in the OD. In contrast,
when cosmic ray muons enter the SK detector, Cherenkov light is detected in both the ID and OD.
Therefore, cosmic ray muons can be identiﬁed by observing the reaction in the OD. As explainedin Section 3.3, energy deposit is measured from the amount of the observed Cherenkov light in
the ID. Since the main purpose of the OD is not to measure the energy deposit but to identify
particles penetrating the OD, only the minimum necessary PMTs are installed. A dead zone exists
with 55 cm thickness between the ID and OD, where the frame supporting the PMT, signal cables,
and high-voltage cables are stored. The wall on the ID side is covered with a black sheet madeof polyethylene terephthalate to suppress reﬂection to clarify the Cherenkov ring and prevent lightfrom passing between the ID and OD (Figure 3.5).

3.3

Cherenkov light

Cherenkov light is emitted by charged particles moving faster than light in a medium.

The

Cherenkov light spreads in a conical shape along the trajectory of the charged particle as shown in
Figure 3.6. The angle θ between the charged particle and the cone is expressed using the ratio β ofthe charged particle velocity v to the light velocity c in a vacuum ascos θ =

1
,
nβ

(3.1)

where n is the refractive index of the medium. Since the refractive index in pure water is about
1.34, a particle moving with a speed close to light (β ≃ 1) in water emits Cherenkov light in the
direction of about 42 degrees. Since the velocity of a charged particle must be greater than the
speed of light in water to emit Cherenkov light, there is a threshold for the momentum of a particle

14

Figure 3.5: The frame supporting the PMT [23].

𝜃Figure 3.6: Diagram of Cherenkov light emission (blue arrow) from a charged particle (red arrow).

15

emitting Cherenkov light, as expressed inp≥ p

1
1 − 1/n2

mc/n.

(3.2)

The threshold of Cherenkov light emission depends on the mass m of a charged particle (Table 3.2).

Table 3.2: Momentum thresholds to emit Cherenkov light.
Particle Mass (MeV/c2 ) Momentum threshold (MeV/c)
Electron
0.511
0.57
Muon
105.6
118
Pion
140
156
Kaon
494
554
Proton
938
1051

The number of emitted photons N per wavelength λ per unit travel distance x of a charged
particle is described as follows.
2πα
d2 N
= 2
dxdλ
λ


1−

1
n2 β 2,(3.3)

where α is the ﬁne structure constant. This equation shows that Cherenkov light in pure water
typically peaks in 300–600 nm and emits about 340 photons per 1 cm particle trajectory. Theenergy of the particle can be estimated by the number of observed photons once the particle type isidentiﬁed. If the energy is relatively low and β is suﬃciently smaller than 1, the energy can also be
obtained from θ given by the Cherenkov radiation angle as Equation 3.1. Therefore, the accuracyof the measurement of the energy is improved by the calibration of the observed number of photons
and the opening angle of the Cherenkov ring.

3.4

Photomultiplier tube

A photomultiplier tube (PMT) is a photo sensor that detects signals by amplifying electrons (photoelectrons) emitted by the photoelectric eﬀect. Incident photons are converted to signals by thefollowing steps (Figure 3.7):
1. A photoelectron is emitted from the photocathode by the photoelectric eﬀect.
2. A photoelectron is accelerated by the electric ﬁeld and collides with the ﬁrst dynode.
3. Secondary electrons produced by the collision are accelerated by the electric ﬁeld and collidewith the second dynode, producing further electrons.
4. Electrons are multiplied by repeated collisions at the dynodes and extracted from the anodeas an electric current.

16

Focusing electrode
DynodeIncident
photonAnodeelectron

PhotocathodeFigure 3.7: Diagram of γ detection by PMT.
Since the amount of electric charge obtained by integrating the output current is proportional to thenumber of photoelectrons, the incident light intensity can be obtained by conversion from the PMT
signal with a relation determined by the calibration. Since the amount of charge generated froma single photoelectron is not always constant but statistically varies, the number of photoelectrons
estimated from the amount of charge is not an integer but a real number. Figure 3.8 shows a pulse
height distribution for single photoelectron signals.

Figure 3.8: Pulse height distribution by single photoelectron electric signals [28].
Super-Kamiokande detects Cherenkov light by PMTs on the wall (Figure 3.9). Large PMTs
with a photocathode diameter of 50 cm (20 inches) are installed in the ID (Figure 3.10), which wasdeveloped by Hamamatsu Photonics by a cooperation with the Kamiokande collaborators [27][28].
The characteristics of this PMT are summarized in Table 3.3. The photocathode area is covered
by Bialkali (Sc-K-Cs) which has high spectral sensitivity for Cherenkov light and low thermionic17

Figure 3.9: Schematic view of Cherenkov light detection [26].

Figure 3.10: Schematic view of the 20-inch PMT [28].

Table 3.3: Characteristics of the 20-inch PMT
Photocathode area
50 cm in diameter
Photocathode material
Bialkali (Sb-K-Cs)
Quantum eﬃciency
22% at λ =390 nm
DynodeVenetian blind type, 11 stages
Gain
107 at ∼ 2000 V
Dark current
200 nA at 107 gain
Dart noise rate
3 kHz at 107 gain
Transit time
90 ns at 107 gain
Transit time spread
2.2 ns (1σ)
Weight13 kg
Pressure tolerance
6 kg/cm2

18

Figure 3.11: Quantum eﬃciency distribution of 20-inch PMT and relative Cherenkov light distribution as it passes through pure water of 15 m [27].
emission. The quantum eﬃciency peaks in the region between 300–600 nm which corresponds tothe peak of the spectrum of Cherenkov light in water (Figure 3.11). At the wavelength of 390 nm,
the quantum eﬃciency is 22%. Magnetic ﬁelds over 100 mG aﬀect trajectories of photoelectrons in
PMT and lead to ineﬃciency and degradation of the time resolution. The geomagnetism of about
450 mG at the location of Super-Kamiokande is reduced to 50 mG by 26 sets of Helmholtz coils
which cover the horizontal direction and vertical direction of the tank [29].
After the PMT chain destruction in 2001, PMTs are covered with cases made of ﬁber-reinforced
plastics (FRPs) and acrylic (Figure 3.12). There are holes on the slide of the cases to allow water toenter and exit the case. The transparency to photons vertically incident in the water of acrylic coveris 96% at a wavelength of 350 nm and does not signiﬁcantly aﬀect the detection of the Cherenkov
light.

19

Figure 3.12: Photo of the case of the 20-inch PMT [25].

Figure 3.13: Schematic diagram of the SK-IV water system [30].

3.5

Water and air puriﬁcation system

Spring water from the Kamioka mine is used for Super-Kamiokande. It is circulated in a water
puriﬁcation system with 60 ton/hour ﬂow rate to get a long attenuation length of water (Figure3.14). Dust and radioactive materials such as radon are removed from water in this system. Radonis the main source of radioactive materials in water and the background of low-energy events suchas neutron capture and solar neutrino. The heat exchangers keep the water temperature at 13 ◦ C
to prevent bacterial growth and keep water transparency, reduce the dark current in PMT, andsuppress convection in the tank. Surviving bacteria are killed in the UV sterilizer.
To prevent radon in the air from dissolving into the water in the tank, radon-free air produced
by a puriﬁcation system is supplied to the top of the tank. Figure 3.14 shows the schematic diagramof the air puriﬁcation system.

20

Figure 3.14: Schematic diagram of the air puriﬁcation system [23].

3.6

Data acquisition system

Analog signals from PMTs are sent to electronics huts located on the top of the tank by coaxial
cables. Lengths of all cables are set to 70 m to make the signal smearing eﬀect and the propagation
time of the signals uniform for all PMTs. In the electronics huts, signals are digitized by front-end
modules, and the integrated charges (corresponding to the number of detected photons) and hit
timings are sent to data acquisition computers. After they are sorted by time, events are selectedby software trigger and recorded (Figure 3.15).

Figure 3.15: Schematic diagram of the data acquisition system for SK-IV [31].
From the SK-IV running period, a new front-end electronics module, QBEE (QTC-Based Electronics with Ethernet), has been used in the data acquisition system (Figure 3.16). In QBEE, QTC
converts analog signals from PMTs to rectangular pulses whose leading edge represents the signal
timing and width represents the integrated charge. The rectangular pulses are then digitized by
TDC and read by FPGA for each 17 µs and recorded after the data compression. About 560 QBEE
boards each with 24 input channels are operating synchronously at 60 MHz in SK.
If the number of hit PMTs in a sliding 200 ns time window exceeds the trigger threshold, asoftware trigger is issued and hits around the trigger timing are recorded. There are four trigger
21

Figure 3.16: Schematic diagram of QBEE [31].
types for the ID in SK-IV as shown in Table 3.4. If the SHE trigger is issued without the OD
trigger (which corresponds to ≥ 7 MeV energy deposit in the ID without response in the OD), a
special trigger called AFT trigger is issued, and all PMT hits are recorded up to 535 µs after the
SHE trigger timing. The AFT trigger information is used for a search to 2.2 MeV γ emitted fromneutron capture on hydrogen as described in Section 7.5.
Table 3.4: Thresholds and time widths of the triggers. The trigger thresholds for SHE and SLE
were changed during SK-IV, and only the initial and last values are shown.
Name Threshold
Time width
SHE
70 → 58[−5 µs, +35 µs]
HE50[−5 µs, +35 µs]
LE47[−5 µs, +35 µs]
SLE
34 → 31[−0.5 µs, +1.0 µs]
OD
22[−5 µs, +35 µs]

3.7Scattering and absorption in water

Light is scattered and absorbed in water, aﬀecting the hit pattern of the Cherenkov ring and thenumber of observed photoelectrons. In SK, lasers of four diﬀerent wavelengths are injected from the
top to the bottom of the detector to calibrate the coeﬃcients of photon scattering and absorption.

22

The coeﬃcients of photon scattering and absorption are parametrized aslI(λ, l) = I0 (λ) exp −,L(λ)1
L(λ) =,
aabs (λ) + asym (λ) + aasy (λ)
P1aabs (λ) = P0 4 + C,
λ

P4
P5
asym (λ) =
1.0
+
,
λ4
λ2


P7
aasy (λ) = P6 1.0 + 4 (λ − P8 )2 ,
λ

(3.4)
(3.5)
(3.6)
(3.7)
(3.8)
(3.9)

where l is the travel length of the light, I0 (λ) is the initial light intensity, I(λ, l) is the light intensity
at l, aabs (λ) is the absorption coeﬃcient, asym (λ) is the symmetric scattering coeﬃcient, aasy (λ)
is the asymmetric scattering coeﬃcient, and P0,...,8 is the constant coeﬃcient. Rayleigh scatteringand the symmetric component of Mie scattering contribute to the symmetric scattering, and theasymmetric component of Mie scattering corresponds to the asymmetric scattering. Figure 3.17
shows the observed water parameters [29].

Figure 3.17: Water parameters in the measurements (points) and ﬁtted functions (lines) [29]. The
total of coeﬃcients (black) consists of the absorption (red), symmetric (blue) and asymmetric(magenta) coeﬃcients.

23

3.8

Energy scale calibration

As explained in Section 3.3, the energy of charged particles can be estimated from the number ofobserved photoelectrons. Calibration of the energy scale for the event reconstruction algorithm,
ﬁTQun (details are described in Chapter 7), is performed using independent control samples withknown energies.
1. Comparison of the number of photoelectrons and track length for cosmic ray muons (1–10 GeV/c)
2. Comparison of the number of photoelectrons and Cherenkov angle for cosmic ray muons(200–500 MeV/c)
3. Invariant mass of π 0 produced by neutrino interactions ( ∼130 MeV/c for the total energy of
two γ’s)
4. Momentum distribution of electrons from muon decay (called as Michel electron) ( <50 MeV/c)
Figure 3.18 shows the results of the absolute energy scale measurements [32]. The total absoluteenergy scale error is 2.09%.

Figure 3.18: Deviations of the absolute energy scale between the data and simulation for eventswith the distance from the ID wall larger than 2 m [32]. The statistical uncertainty is denoted by
the vertical error bars. The momentum range of each sample is shown by the horizontal error bar.

24

Chapter 4

Simulation
Although the proton decay events have not been observed yet, their detection eﬃciency is estimatedby the Monte Carlo (MC) simulations. The proton decay simulation with 10 5 events was used toestimate the signal eﬃciency. Details of the simulations are described in this chapter.
The simulations consist of two steps, the primary interaction in the nucleus and the propagationof the particles in the detector. The interaction simulation determines the number of particles andthe particle types, positions and momenta. Then, their trajectories in the detector, emission of
Cherenkov light, propagation of photons in water, and the detection of photons are calculated by
the detector simulation.
After the simulation, data reduction (Chapter 6), event reconstruction (Chapter 7) and selectionof proton decay signals (Chapter 8) are applied to the output of the detector simulation as well as
the observed data.

4.1

Proton decay

In the proton decay simulation, particles were generated from nuclei of hydrogen and oxygen inwater within the inner detector (ID) of Super-Kamiokande (SK). Protons in a hydrogen atom arecalled free protons and the momenta are negligible to the momenta of outgoing particles from theproton decay. In this case, the momenta of µ+ and K 0 from the proton decay of the p → µ+ + K 0
channel are the same (326.5 MeV/c) and their directions are back-to-back due to the two-bodydecay. In contrast, protons in oxygen interact with other nucleons and are called bound protons.
These bound protons experience Fermi motion. Due to Fermi motion, momentum distributions ofµ+ and K 0 from the decays of bound protons are diﬀerent than those from free protons. Figure
4.1 shows momentum distribution for each proton state. The Fermi momentum distribution in thesimulation is implemented using the measurement of electron-12 C scattering [33]. The eﬀective massof the bound proton is smaller than that of free protons due to interactions with the surrounding
nucleons. The eﬀective proton mass m′p is calculated by subtracting the binding energy Eb from the25

proton rest mass mp , where Eb is simulated for each nuclear state as a Gaussian random variablewith a mean and a standard deviation of 39.0 MeV and 10.2 MeV for the s-state and 15.5 MeV and3.82 MeV for the p-state, respectively. The ratio of protons in the s-state and p-state is taken to be

Number of events

1:3 based on the nuclear shell model [34].

104
103
102
10
10

100 200 300 400 500 600 700 800
Proton momentum (MeV/c)

Figure 4.1: Proton momentum distributions of free (black), s-state (red), p-state (green) andcorrelated decay protons (blue).
Protons in an oxygen nucleus often decay as a pair with a nearby nucleon with an estimatedprobability of 10% for each proton decay, referred to as correlated decay [35].

4.2

Final state interaction

After the proton decay or neutrino interaction (described in Chapter 5), generated hadrons interact
with surrounding nucleons. A custom semi-classical intranuclear cascade model is used for thesimulation of hadron interactions in the nucleus. The scattering probability is calculated from themean free path for every 0.2 fm step in the nucleus assuming the nucleon density of oxygen [36].
Pion scattering cross section in the low momentum region is based on the model by Salcedo et al.
[37]. Fits to pion-nucleon scattering data are used at high momentum (more than 500 MeV/c) toestimate the mean free path. The kaon interaction model was updated since the previous analysis[1] as explained below.
In the nucleus, neutral kaon scattering is simulated assuming an eigenstate of K 0 . Since K +and K 0 form an isospin doublet, the eigenstates of isospin can be represented by the superpositionsof the scattering amplitudes of kaon and nucleon:|1, 1⟩=

K +, p ,(4.1)|1, 0⟩=√( K + , n + K 0 , p )/ 2,(4.2)|1, −1⟩=

K 0, n ,(4.3)|0, 0⟩=√( K + , n − K 0 , p )/ 2.

(4.4)

26

1

Fraction

0.8
0.6

No interaction

0.4↓ Charge exchange

0.2

Scattering
00.5
1Kaon initial momentum (GeV/c)

1.5

Figure 4.2: Fraction of K 0 interactions in the 16 O nucleus as a function of K 0 momentum. Elastic
scattering (red) and charge exchange (blue) are simulated in the nucleus.
Therefore, kaon-nucleon scattering amplitudes can be represented as follows:
K +, p=|1, 1⟩ ,K +, n=√(|1, 0⟩ + |0, 0⟩)/ 2,

K 0, n=|1, −1⟩ ,K 0, p=√(|1, 0⟩ − |0, 0⟩)/ 2.

(4.5)
(4.6)
(4.7)
(4.8)

The scattering amplitudes of K 0 and a nucleon in the simulation were evaluated using (4.7) and(4.8) in which the eigenstates of isospin were obtained by the measurements of K + scattering [38]corresponding to (4.5) and (4.6). Only elastic scattering K 0 + N → K 0 + N , where N is a nucleon(proton p or neutron n), and charge exchange K 0 + p → K + + n are simulated in the nucleus
as they are relevant for K 0 N interactions at the proton decay energy scale. If the momentum of
the scattered nucleon is below the Fermi surface momentum, scattering is suppressed by the Pauli
exclusion principle. The cross section of K 0 and 16 O scattering is calculated by scaling the crosssection of the K 0 and nucleon scattering assuming they are proportional to the number of nucleons,referring to the experimental results [39]. Figure 4.2 shows the fraction of K 0 interactions in the
nucleus as a function of K 0 momentum.

4.3

Detector simulation

After the hadrons leave the nucleus, the detector simulation determines how the generated charged
particles are observed by the detector. The simulation is performed by a program called SKDETSIM. It is composed of a Fortran-based toolkit called Geant3 [40].
Since it is known that the Geant model for π meson scattering does not reproduce the data
correctly, pion scattering is calculated by NEUT program [41]. The kaon interactions in water are27

Cross section (mb)600
400
Inelastic scattering
200
Elastic scattering
0

0.5

1
Momentum (GeV/c)

1.5

Figure 4.3: Cross section of KL0 in water. In the simulation of KL0 interactions, inelastic scattering(gray) and elastic scattering (red) are implemented.
simulated as the eigenstates of KL0 and KS0 . The K 0 leaving the nucleus is randomly assigned to
KL0 or KS0 . Figure 4.3 shows the cross section of KL0 in the simulation. KL0 and KS0 cross sections
are calculated from the K 0 and K̄ 0 scattering amplitudes.
|KL , p⟩=

|KS , p⟩=√
( K 0 , p + K¯0 , p )/ 2,√( K¯0 , p − K 0 , p )/ 2.

(4.9)(4.10)The K 0 N scattering amplitude is derived from the K + N cross section [38] as written above while
the K̄ 0 N scattering amplitude is calculated from the K − N cross section [42]. Figure 4.4 shows
the distance between the positions of proton decay and the end point of KL0 track, either by the
decay or inelastic scattering. The mean free path calculated from the cross section is about 1 m in
the water while the average travel length of KL0 with 326.5 MeV/c momentum is 13 m if it is in a
vacuum. This means KL0 from proton decay scatters in water before it decays in most cases. The
ﬂight length becomes shorter than that in a vacuum if hadronic inelastic scattering occurs before
KL0 decays.
In contrast, KS0 immediately decays in water with a lifetime of 90 ps. Coherent regenerationfrom KL0 to KS0 is also implemented both in the oxygen nucleus and in water. The regeneration
probability in water is evaluated based on the results of a kaon scattering experiment using a carbon
target [43]. In the oxygen nucleus, the regeneration probability is assumed to be proportional to its
density relative to water. The fraction with regeneration is about 0.1% of the total p → µ+ + KL0
events in the simulation.
The Cherenkov light generation and propagation in water and the PMT response are calculated
by a custom program developed by the SK collaboration. When a charged particle moves in waterwith a velocity larger than the propagation of light, photons are produced by Cherenkov radiation(Section 3.3) along the particle’s trajectory. These photons are propagated in water in which
28

Arbitary unit (area normalize)

1
10−1
10−2
10−3
10−4
−5

10

0

2

4
6
Distance (m)810Figure 4.4: Distance between vertices of primary proton decay and the end point of KL0 track.
Total p → µ+ KL0 events (black) and the breakdown of KL0 decay events (red) and hadronic inelastic
scattering events (blue) are shown.
Rayleigh scattering, Mie scattering, and absorption with water molecules are accounted for based on
the parameters measured by the calibration of the detector (Section 3.7). The detection of photons
by PMTs is simulated using quantum eﬃciency and collection eﬃciency based on measurements.
Once the photons are detected by the PMTs, each of them is converted to a single photoelectron
signal and accumulated as the output signal. The ﬂuctuation of the single photoelectron signal andthe timing resolution are also simulated based on the measurements.

29

Chapter 5

Atmospheric neutrino
As explained in Chapter 6, atmospheric neutrino events are the only possible background to thesearch for proton decay. The expected background rate in the proton decay search is also estimatedusing the atmospheric neutrino simulations as well as the proton decay simulation in Chapter 4.
Although the live time of SK used for this analysis is about 10 years, atmospheric neutrino MC
samples equivalent to 500 years of operation are produced for this analysis to reduce statisticaluncertainty in the background estimation. The number of the remaining background events in theMC simulation is then normalized by the live time to estimate the expected number of background
events in the observed data.
Furthermore, atmospheric neutrino events are not only background events for proton decay but
also important physics targets. Section 5.4 describes a new atmospheric neutrino event classiﬁcation method using machine learning. The sensitivity of the neutrino oscillation measurements is
improved by this method.

5.1

Neutrino ﬂux

Cosmic rays are continuously raining down on the earth from space. Among various particles,
the primary cosmic rays above several hundred MeV energy mainly consist of protons with a small
fraction (about 9%) of helium nuclei. Most of the primary cosmic ray particles collide with nuclei inthe atmosphere and produce secondary hadrons, mainly pions and kaons. The neutrinos produced
from their decay are called atmospheric neutrinos (Figure 5.1). Atmospheric neutrinos are thedominant source of background in the proton decay search.
Atmospheric neutrino ﬂux is calculated by independent groups such as Honda [46][47][48], Fluka[49] and Bartol [50]. Super-Kamiokande mainly employs Honda ﬂux for analyses, and it is comparedwith the other models to evaluate the systematic uncertainty caused by the theoretical calculation.
As is clear from the production process, atmospheric neutrino ﬂux depends on the amount ofprimary cosmic rays. The atmospheric neutrino ﬂux calculation uses the data observed above the30Figure 5.1: Schematic diagram of collisions between cosmic rays and nuclei in the atmosphere [44].
atmosphere to determine the primary cosmic ray ﬂux as the input. The primary cosmic ray ﬂux
decreases rapidly with increasing energy as shown in Figure 5.2, and the atmospheric neutrino ﬂux
also decreases. Figure 5.3 shows the simulation of the atmospheric neutrino ﬂux and the ﬂux ratios
at the Super-Kamiokande site. Collisions between primary cosmic rays and nuclei in the atmospheregenerate mainly π mesons and K mesons. Charged pions decay in ﬂight as follows:
π + → µ+ + νµ → e+ + νe + ν̄µ + νµ(5.1)

π − → µ− + ν̄µ → e− + ν̄e + νµ + ν̄µ(5.2)

The right-hand plot in Figure 5.3 shows the ratios of atmospheric neutrino ﬂux (νµ + ν̄µ )/(νe + ν̄e )in the simulation. The ratio is close to two below the neutrino energy of 1 GeV as expected from
Equation 5.1 and 5.2. At high energies, the ratio is greater than two and increases with the energy
because the eﬀective lifetime of the muon is longer for higher energy due to relativistic eﬀects andsuch high-energy muons reach the ground before they decay.

5.2

Neutrino oscillation

Flavor eigenstates of neutrinos (νe , νµ , ντ ) are superpositions of mass eigenstates (ν1 , ν2 , ν3 ). Neutrinos are produced in the ﬂavor eigenstates in the atmosphere, and the ﬂavor compositions change
before the observation in SK (known as neutrino oscillations). For example, the oscillation probability of two-generation oscillation is

P (να → νβ ) = sin 2θ sin
2

2

31


1.27∆m2 [eV2 ]L[km],
E[GeV](5.3)

Figure 5.2: Primary cosmic ray ﬂux of Honda model (solid line) and observations by several experiments (points) [46]. Calculations in a diﬀerent way discussed in the paper (dashed line) and theprevious model they used (dotted line) [45] are also shown.

Figure 5.3: The atmospheric neutrino ﬂux at Kamioka averaged over all directions (left) and theﬂux ratio (right) compared with other ﬂux models, HKKM06 (previous calculation of this model),model by Bartol group [50] and model by FLUCA group [49]. Taken from [48].

32

where να and νβ are ﬂavor eigenstates, θ is the mixing angle of mass eigenstates, ∆m2 is the
mass-squared diﬀerence of the neutrino mass eigenvalues (∆m2 := m22 − m21 , where m1 and m2 areneutrino masses of ν1 and ν2 ). L and E are ﬂight length and the energy of the neutrino. P (να → νβ )can be nonzero if neutrinos have ﬁnite mass. Although neutrinos have three generations νe , νµ andντ , if L/E is suﬃciently smaller than the mass-squared diﬀerence, oscillations between speciﬁc
generations become negligible. In this case, the oscillation probability is approximated by twogeneration oscillation as expressed by Equation 5.3. This approximation is valid for atmospheric
neutrino oscillation. The contribution from νµ → ντ oscillation is dominant for atmospheric neutrino observation in SK, and therefore the eﬀect to the background estimation of the proton decay
search is accounted for assuming two-ﬂavor neutrino oscillation with parameters of sin 2 2θ = 1.0
and ∆m2 = 2.5 × 10−3 eV2 . Figure 5.4 shows the survival probability P (νµ → νµ ) of neutrino
predicted from νµ and ντ two-generation neutrino oscillation with parameters of sin2 2θ = 1.0 and∆m2 = 2.5 × 10−3 eV2 . Down-going muon neutrinos coming from the atmosphere above the detector rarely oscillate around 1 GeV energy due to the short ﬂight length of 10–100 km. In contrast,
the up-going muon neutrinos coming through the earth appears as a mixture with tau neutrinos
and only half of them are observed as muon neutrinos on average around 1 GeV energy. As theenergy threshold of the tau neutrino interactions is 3.5 GeV, which is higher than the typical energy
of atmospheric neutrinos, they are rarely observed by the detector.

1
Probability

0.8
0.6
0.4
0.2
0
1

10

102
L/E (km/GeV)

103

104

Figure 5.4:
The survival probability P (νµ → νµ ) of neutrino predicted from two-generation
neutrino oscillation with parameters of sin2 2θ = 1.0 and ∆m2 = 2.5 × 10−3 eV2 .

5.3

Neutrino interaction

Neutrino interactions are simulated by the NEUT program [41]. The neutrino interactions in thedetector can be divided into two categories: charged current (CC) and neutral current (NC). In thecase of the (quasi) elastic scattering, these interactions are expressed as:

33

• Charged current quasi-elastic (CCQE) scattering: νl + n → l + p• Neutral current elastic (NCE) scattering: νl + N → νl + N
νl

νl

l

Z

W

n

νl

p

N

N

Figure 5.5: Feynman diagrams of νl + n → l + p (left) and νl + N → νl + N (right).

Here, n is a neutron, p is a proton, and l is a charged lepton (e, µ, τ ). Figure 5.5 shows the Feynman
diagrams for these interactions. In the CC interaction, a neutrino is converted to a charged lepton
via the exchange of a W boson. In contrast, the charge is invariant in the NC interaction via theexchange of a Z boson. Although neutrinos are not directly observed by the detector, it is possible toidentify the ﬂavor of the incident neutrino if the ﬂavor (e, µ, τ ) of the outgoing lepton is identiﬁed in
the case of the CC interaction. As described in Section 3.1, Cherenkov rings made by electrons andmuons can be distinguished in SK by their shapes, which are called e-like and µ-like, respectively.
In contrast, particles with diﬀerent charges, such as e− and e+ , cannot be distinguished from thering shape. Since the NC reaction is ﬂavor-independent and the outgoing neutrinos cannot be
observed, the original neutrino ﬂavor cannot be identiﬁed. The cross section of CCQE scattering
was formalized by the Llewellyn Smith model [52] and implemented in the simulation. As well as
proton decay simulation, the scattering with hydrogen nucleus can be treated as a scattering with afree particle. In contrast, the eﬀects of Fermi momentum and Pauli blocking are considered in the
scattering with oxygen nuclei. A nucleus model is based on the local Fermi-gas model by Nieves etal. [53] and includes updates by Bourguille et al. [54].
In addition to the (quasi) elastic scattering, multi-nucleon scattering is implemented in thesimulation:
ν + n + N → l + p + N.

(5.4)

This scattering is accompanied by several nucleons. Accelerator-based neutrino oscillation experiments such as K2K, MINOS and MiniBooNE found that the number of the observed CCQE-likeevents was larger than predicted by the models [55][56][57]. It indicates the existence of multinucleon interaction, and it is implemented as the Valencia model by Nieves et al. [53] in NEUT.
Elastic and quasi-elastic scatterings are dominant below 1 GeV (Figure 5.6). Above 1 GeV,

34

Figure 5.6: CC cross section of neutrino (a) and antineutrino (b) [51]. Total cross section (solid),
quasi-elastic scattering (dotted), single meson production (dashed), DIS (chain) and observations(points) are shown.

35

hadrons are generated in addition to the outgoing leptons. Single meson production
ν + N → l + N ′ + meson (π, K, ρ, . . . ),(5.5)

ν + N → ν + N ′ + meson (π, K, ρ, . . . )(5.6)

ν + n → l + p + γ,

(5.7)

ν+N →ν+N +γ(5.8)

and single gamma production

are implemented in NEUT, and the pion production is the main reaction around 1 GeV. N and
N ′ are the initial and ﬁnal states of the nucleon. Single pion production is based on the Berger
and Sehgal model [58]. The model consists of two steps. First, an intermediate baryon resonancestate (N ∗ ) is produced from hadronic current (ν + N → l + N ∗ ). Second, it decays into a pion anda baryon (N ∗ → π + N ′ ). Since the reaction is divided into two stages, cross sections for single
kaon, eta and gamma production also can be calculated by changing parameters such as the decay
probability.
Pion production is the major interaction among meson productions. Coherent pion production
ν + 16 O → l + 16 O + π + ,(5.9)

ν + 16 O → ν + 16 O + π 0(5.10)

ν + p → l + p + π+ ,(5.11)

ν + p → ν + p + π0(5.12)and diﬀractive pion production

are implemented in addition to the single pion production from baryon resonance decay. Coherent
pion production is a pion production without breaking the target nucleus. The model is basedon Berger and Sehgal [59]. Diﬀractive pion production is also a pion production with a proton
implemented by Rein [60].
In the high energy region of more than several GeV, deep inelastic scattering (DIS) is a dominantinteraction:
ν + N → l + N ′ + hadrons,(5.13)

ν + N → ν + N ′ + hadrons.

(5.14)The interaction can be described as the scattering between a neutrino and a quark. When theinvariant mass of the hadronic system (W ) is below 2 GeV/c2 , only multiple meson productions are
considered to avoid double counting of single meson production events with the cross section as a36

function of the pion multiplicity which depends on W . The mean multiplicity of charged hadrons
is estimated from the results of bubble chamber experiments [61][62]. Interactions with W larger
than 2 GeV/c2 are simulated by PYTHIA [63].

5.4

Event classiﬁcation of atmospheric neutrino

The discovery of neutrino oscillations revealed the existence of neutrino mass. However, it is
not determined yet whether the order of mass (called mass ordering or mass hierarchy) is normal(m1 < m2 ≪ m3 ) or inverted (m3 ≪ m1 < m2 ) although the normal ordering is indicated at 91.9%–
94.5% conﬁdence level by SK-T2K combined analysis [64]. Observation of atmospheric neutrinos
has a sensitivity to the mass ordering as the oscillation probability depends on it. Upward-going
atmospheric neutrinos come from the opposite side of the earth and are observed in the detector.
Figure 5.7 shows the atmospheric neutrino oscillation probabilities for each mass ordering. The
oscillation probabilities depend on the neutrino energy and zenith angle θ (Figure 5.8). If it isa normal ordering, the neutrino oscillation probability P (νµ → νe ) is enhanced for upward-going
events in a few GeV energy due to the neutrino interaction with the earth (matter eﬀects). On theother hand, P (ν̄µ → ν¯e ) is enhanced if it is an inverted mass ordering. Therefore, observation of νeand ν̄e events in a few GeV leads to the understanding of the neutrino mass ordering.
In a few GeV energy, multi-ring events are observed as DIS is a dominant interaction. The
identiﬁcation of the neutrino ﬂavor for multi-ring events is more diﬃcult than that for single-ring
events mainly from CCQE scattering below 1 GeV. The neutrino ﬂavor can be estimated from theparticle type of the lepton in CCQE scattering (single-ring event), while Cherenkov rings fromhadrons are observed in addition to the ring from the charged lepton. Currently, atmospheric
neutrino events are classiﬁed by the log-likelihood (LL) method based on the variables reconstructedfrom the PMT hits. To improve the event classiﬁcation eﬃciencies, an alternative method basedon the neural network (NN) has been developed for multi-GeV multi-ring event samples (multiGeV events have Evis above 1.33 GeV, where Evis is electron equivalent total energy deposit in thedetector) in the SK-IV period.
Events are classiﬁed by the NN from the input variables based on the knowledge accumulated
from the training samples. We employed input variables listed in Table 5.1, following the current
analysis based on the LL method [64]. Figure 5.9 shows the distributions of these variables. The
atmospheric neutrino simulation is compared with the data of downward-going neutrino events(cos θ > 0.4, θ is a zenith angle) to avoid bias in the resonance region.
Multi-GeV multi-ring events are classiﬁed into 4 types.
• νe -like: νe charged current (CC) interaction as a correct class• ν̄e -like: ν̄e CC interaction as a correct class37

P (νµ → νe ), normal ordering

P (ν̄µ → ν̄e ), normal ordering

P (νµ → νe ), inverted ordering

P (ν̄µ → ν̄e ) inverted ordering

Figure 5.7: Atmospheric neutrino oscillation probability P (νµ → νe ) and P (ν̄µ → ν̄e ) assuming
normal (top) and inverted (bottom) mass ordering. The horizontal axis shows the neutrino energy,
the vertical axis shows the zenith angle to the entering neutrino in the detector, and the color scaleshows the oscillation probability assuming |∆m232 | = 2.5 × 10−3 eV, sin2 θ23 = 0.5, sin2 θ13 = 0.021,
δCP = 0.

𝜃

Figure 5.8: Schematic diagram of the deﬁnition of the zenith angle θ for the atmospheric neutrino.
The large circle is the earth, the cylinder is the SK detector, and the arrow is the atmospheric
neutrino.

38

Figure 5.9: Distributions of the number of rings, Evis , Tmom , Fmom , PID, the number of Michelelectrons and Ldecay e in the SK-IV period. The data (points) with statistical error bar are comparedwith the atmospheric neutrino simulation of νe CC (red), ν̄e CC (green), νµ and ν̄µ CC (blue) andNC (gray) events. Downward-going events (cos θ > 0.4, θ is a zenith angle) after FCFV selection(described in Section 6.6) are shown. The atmospheric neutrino simulation is normalized by thelive time.
39

Name
Nring
Evis
Tmom
Fmom
Ndecay e
Ldecay e
PID

Table 5.1: Input variables used in NN
Deﬁnition
Number of rings
Electron equivalent total energy deposit
Fraction of transverse momentum to the total visible energy
Fraction of the momentum of the most energic ring to the total visible energy
Number of Michel electron
Distance between Michel e vertex and primary vertex
Likelihood of particle type identiﬁcation of the most energic ring

• µ-like: νµ and ν̄µ CC interaction as a correct class• Others: Neutral current (NC) interaction and ντ and ν̄τ CC interaction as a correct class
The neural network consists of three layers. The activation function is ReLU, and the batchnormalization [65] is adopted before each activation. The neural network is trained by 300 yearsof atmospheric neutrino simulation. The performance is evaluated with respect to the eﬃciency
and contamination by using 100 years of atmospheric neutrino simulation which is independent of
the training samples. The eﬃciency is deﬁned as the fraction of correctly identiﬁed events in eachinteraction (νe CC, ν̄e CC, νµ or ν̄µ CC, NC or tau neutrino events). The contamination is deﬁned
as the fraction of misidentiﬁed events in each category (νe -like, ν̄e -like, µ-like, others).
Results of the classiﬁcation are shown in Table 5.2. Compared to the LL method, νe eﬃciency
is increased by a factor of 1.3 by the NN method with the comparable contamination.
Table 5.2: Results of the selection by NN
Eﬃciency(%) Contamination(%)
LL
NN
LL
NN
νe -like 34.4
46.7
54.0
52.8
ν̄e -like 60.4
63.1
76.1
76.1
µ-like 77.0
85.9
7.6
8.9
others 55.6
47.4
56.7
39.6

Figure 5.10 shows νe event distributions of transverse momentum Tmom and fraction of mostenergic ring momentum Fmom . These show that the NN selected additional νe events in high Tmomand low Fmom region with respect to the LL method. This indicates that the NN is optimized to
identify events with energetic hadrons and improves νe selection eﬃciency.

40

Figure 5.10: Distributions of the fraction of the transverse momentum (left) and the most energic
ring momentum (right). From total νe CC events (black), the NN method (red) and LL method(blue) classiﬁed νe -like events.

41

Chapter 6

Data reduction
As described in Section 3.2, if a proton in water decays and the generated charged particles deposit
all energy in the ID, Cherenkov light is observed only in the inner detector (ID) and no response
should be observed in the outer detector (OD). These events are called fully contained events (FCevents), and they are selected for the proton decay analysis. Similarly, if a neutrino interacts witha nucleus or electron in water and generated charged particles that deposit all of their energy in
the ID, they are also selected as the FC events and used for the study of neutrino oscillation andother physics. The FC atmospheric neutrino event can be the background of the proton decay
search. The number of FC atmospheric neutrino events observed in SK is about ten events perday. In addition to neutrinos, cosmic ray muons are observed by the SK detector. Two cosmic raymuons per 10 cm2 per second are poured on the Earth, but since SK is located 1000 m underground(which is equivalent to 3000 m of water), the muon rate is suppressed by a factor 10 5 and two muons
per second are observed in the SK detector with a diameter of about 40 m. The emission of light
due to radioactive isotopes in the detector and discharges in the circuit of PMTs (called ﬂasher
events) is also observed as background events. These events must be removed to search for proton
decay events. The FC reduction in the SK-IV period consists of several steps as described in thischapter. To eﬃciently remove backgrounds, simple conditions are applied ﬁrst, and in later steps,
the remaining backgrounds are removed by the conditions optimized for the characteristics of eachtype of background. In the SK-II running period, some cut thresholds were diﬀerent because thenumber of photomultiplier tubes (PMTs) was about half of the other periods. The FC reduction
is applied not only for the data but also for the simulation which is used for the estimation of thesignal eﬃciency and expected backgrounds in Chapter 8.
Figure 6.1 shows examples of FC atmospheric neutrino, cosmic ray muon, ﬂasher and low-energy
events.

42

Figure 6.1: Typical event displays of FC atmospheric neutrino (left top), cosmic ray muon (righttop), ﬂasher (left bottom) and low-energy (right bottom) events. The points show the position of
PMTs in the ID, and the size of the circles corresponds to the observed charge. The small plot in
the top left side for each ﬁgure is the OD hit PMT distribution.

43

6.1

First reduction

About 106 events are recorded in a day and most of them are due to cosmic ray muons or low-energy
radioactive backgrounds. The following conditions are applied as the ﬁrst step in the FC reduction
to remove such backgrounds.
1. PE300 >200 photoelectrons (p.e.’s)
The total number of p.e.’s observed in the ID within 300 ns (PE300 ) is more than 200 p.e.’s.
2. NHITA800 ≤ 55
The number of hit PMTs in the OD in ±400 ns around the event trigger time (NHITA800 ) is
below 55.
The threshold of the ﬁrst condition corresponds to a momentum of 22 MeV/c for an electron. Mostof the low-energy events such as radioactive backgrounds (several MeV energy scale) are rejectedby this cut. The second condition is intended to cut cosmic ray muons passing through the OD.
The event rate in SK is reduced from 106 events/day to about 3000 events/day by these cuts.

6.2

Second reduction

In the second step of the FC reduction, events satisfy the following conditions are selected.
1. PEmax /PE300 < 0.5
The ratio of the highest number of p.e.’s in all hits of the ID PMTs (PE max ) and the total
number of photoelectrons (PE300 ) should be below 0.5.
Figure 6.2 shows the PEmax /PE300 distribution.

Figure 6.2: PEmax /PE300 distribution of data (left) and atmospheric neutrino simulation (right)
after the ﬁrst reduction in SK-IV [67]. The atmospheric neutrino events are normalized by SK-IV
live time.

44

2. NHITA800 ≤ 30 or PEtot ≥ 100000 p.e.′ s
The number of hits in the OD (NHITA800 ) should be below 25, or the total number of p.e.’s
in the ID (PEtot ) should be more than 100000. Figure 6.3 shows the NHITA800 distribution.

Figure 6.3: NHITA800 distribution of data (left) and atmospheric neutrino simulation (right) after
the ﬁrst reduction in SK-IV [67]. The atmospheric neutrino events are normalized by SK-IV livetime.

The ﬁrst condition rejects events that have one PMT with signiﬁcantly larger signals than other
PMTs such as due to the electric noise of PMTs. The second condition rejects cosmic ray events
left from the ﬁrst FC reduction with a tighter cut while very high energy FC events with lightleakage into the OD side are accepted. The event rate is reduced to about 600 events/day by these
cuts.

6.3

Third reduction

The third reduction was designed to reduce speciﬁc backgrounds.
• Hard muon cut
Cosmic ray muons which have energy more than 1 TeV, hard muon, are rejected by NHITA500 ≥40. Here NHITA500 is the number of p.e.’s within the OD sliding 500 ns time window set at
a time to maximize it.
• Through-going muon cut
Through-going muon ﬁtter is applied to events that have more than 1000 ID PMT hits with
PEmax > 230 p.e.’s. The ﬁtter ﬁnds the entering and exiting points of the through-going
muon in the ID and returns goodness which represents the consistency of the observed PMT
signals and the prediction from muons based on the time information. Using this ﬁtter, therejection criteria are deﬁned as follows:45

– Goodness > 0.75
– NHITAin800 ≥ 10 or NHITAout800 ≥ 10
NHITAin800 and NHITAout800 are the number of hit PMTs located within 8 m from enteringand exit points, respectively. These criteria reject events with responses in the OD at theentering and exiting points of the muon.
• Stopping muon cut
As well as through going muon cut, the following cuts are applied using the ﬁtter to ﬁnd theentering point.
– Goodness ≥ 0– NHITAin800 in ≥ 10
It rejects cosmic ray muons that enter from outside the detector and stop in the ID. In thiscase, the OD shows a response only near the entering point. Figure 6.4 shows examples ofthrough-going muon and stopping muon events.

Figure 6.4: Typical event display of through-going muon (left) and stopping muon (right) events.

• Cable hole muon cut
The high voltage and signal cables from the PMTs are bundled into twelve and pulled out
through four holes in the top of the tank (Figure 6.5). If a cosmic ray muon passes through thiscable bundle, it will enter the ID without leaving PMT hits in the OD. In that case, it cannot
be removed by the ﬁrst and second reductions. To reject these events, plastic scintillatorswith 2 m×2.5 m are installed on the cable holes as VETO counters. Events are rejected if thefollowing criteria are satisﬁed.
– There is a signal from one of the plastic scintillators
46

Figure 6.5: Schematic diagram of cable hole muon events and a VETO counter [66].
– A distance from a reconstructed vertex (position where the charged particle is generated)in the ID to the cable hole is within 4 m.
• Coincidencce muon cut
If a cosmic ray muon event happens to coincide with a low energy event, the timing of the OD
signals is shifted and the muon is not removed by the former cuts. These events are rejectedby the following rejection criteria.
– NHITAoﬀ ≥ 20
The number of OD hits within the time window of 400–900 ns from the event trigger(called oﬀ-time window), NHITAoﬀ , is more than 20.
– PEoﬀ > 5000 p.e.’s
The number of p.e.’s in the ID within the oﬀ-time window, PEoﬀ , is more than 5000.
Figure 6.6 shows an example of coincidence muon events.
• Flasher event cut
The spontaneous discharge from PMTs causes events called ﬂasher events. Flasher events
often have a broad timing distribution compared to events produced by charged particles.
Events are rejected if they satisfy the following condition.
– NMIN100 ≥ 20
The minimum number of ID hits in a sliding 100 ns time window from +300 ns to +800 ns
after the trigger, NMIN100 , is more than 20.
47

Figure 6.6: Typical event display of coincidence muon events.
• Low energy event cut
The remaining low-energy backgrounds are rejected if the events satisfy the following criteria.
– The number of ID hits is less than 500– N50 ≤ 50 (25 for SK-II),N50 is the maximum number of ID PMT hits in a sliding 50 ns window after subtracting
the photon time-of-ﬂight. The threshold of these cuts corresponds to roughly 10 MeV
for an electron, which is suﬃciently lower than the energy expected from proton decay.
The event rate in SK is reduced to about 50 events/day by these cuts.

6.4

Fourth reduction

The fourth reduction is performed to remove ﬂasher events that could not be removed by the third
reduction. Since ﬂasher events often have similar hit patterns (Figure 6.7), the remaining eventsare compared with the other events. The ID is divided into 1450 regions of 2 m×2 m squares, andthe correlation r between events is evaluated using the following equation.
1 X (QA,i − ⟨QA ⟩)(QB,i − ⟨QB ⟩)
r=
,
N i=0
σA σB
N

(6.1)

where QA,i (QB,i ) is the number of p.e.’s of event A (B), ⟨QA ⟩ (⟨QB ⟩) is the mean of them, andσA (σB ) is the standard deviation. If the correlation between two events exceeds a threshold, the48

Figure 6.7: Typical event displays of ﬂasher events. The same ID PMT causes the ﬂasher events.
hit patterns are determined to be similar. The Kolmogorov-Smirnov test statistic dks is also used
for the reduction. It is performed for the accumulated charge distribution normalized by the total
observed charge. The cut parameter of the fourth reduction is introduced by the likelihood method.
The likelihood value is calculated from probability density functions for these two parameters using
the atmospheric neutrino simulation.
Prob =
Probr=

Probdks=

1(Probr + Probdks ),
2

!
top
10
X10
X
1
i
−2
ln P r(ri , Qi ) + 2
ln+ Cr ,
10
Nα
i
i=1

!
top
10
X10
X
i
1
−2
ln P r(dksi , Qi ) + 2
ln+ Cdks ,
10
Nβ
i
i=1

(6.2)
(6.3)

(6.4)

where P r is a probability function, Qi is the average total charge of two compared events, N is thenumber of compared events, and α, β, Cr and Cdks are tuning parameters with the values of 0.909,
0.674, −0.31 and −3.39, respectively. Only 10 combinations with the largest correlations are used
for the likelihood calculation. Events with a likelihood value of more than 3.0 are rejected.

6.5

Fifth reduction

The remaining backgrounds are rejected by the ﬁfth reduction.
• Stopping muon cut
Events satisfying either of the following conditions are rejected as remaining stopping muon
events.
– Stopping muon ﬁt goodness ≥ 0.5 and NHITAin800 ≥ 5

49

– NHITAin200 ≥ 5
NHITAin200 is the maximum number of OD hits located within 8 m from the entering
point and in a sliding 200 ns time window.
– The vertex reconstruction goodness is below 0.77, the number of ID hits is more than7000 with more than 70000 p.e.s deposited therein, and NHITAin800 ≥ 6
The second criterion rejects events with the OD PMT hits localized in terms of time andlocation. The third criterion rejects high-energy cosmic-ray muons which are incorrectly
reconstructed in the ﬁducial mass region by the event reconstruction algorithm.
• Cable hole muon cut
Events satisfying the following conditions are rejected as the events near the cable hole without
signals from VETO counters.
– Stopping muon ﬁt goodness ≥ 0.4– PEtot > 4000 p.e.’s– cos θz > 0.6
cos θz is the reconstructed muon direction, and cos θz = 1 is vertically downward-going.
– dhole < 250 cm
where dhole is the distance between the reconstructed entering point and the closest cablehole position.
• Coincidence muon cut
Events satisfying the following criteria are rejected as accidental coincidence of muon andother events:– PE500 < 300 p.e.’s
PE500 is the total number of p.e.s in a ﬁxed 500 ns time window from −200 ns to +300 ns
after the trigger.
– NHITALATE200 ≥ 16
NHITALATE200 is the maximum number of the OD hits in a 200 ns sliding time window
from +300 ns to +1500 ns after the trigger.
• Invisible muon cut
Invisible muons are cosmic ray muons that have momentum below the Cherenkov thresholdin the ID. They reach the inner water tank and decay into electrons after 2.2 µs on average.
Cherenkov light emitted from Michel electron is observed in the ID, while the OD PMT hits50

are produced earlier than the ID trigger timing. Therefore, the ID hits from Michel electron
are misidentiﬁed as an FC event. These events are rejected by requiring either of the following
conditions to identify the visible muons:– PEtot ≤1000 p.e.s (500 p.e.s for SK-II) and NHITAearly200 ≥ 5 and NHITAsum ≥ 10.
NHITAearly200 is the maximum number of OD hits in a sliding 200 ns time window from
−9000 ns to −200 ns after the trigger. NHITAsum is deﬁned as

NHITAearly200 + NHITA500 (dOD < 500 cm)
NHITAsum =NHITAearly200(otherwise)(6.5)

where NHITA500 is the number of OD hits in a ﬁxed 500 ns time window from −200 ns to+300 ns after the trigger, and dOD is the distance between the position of the OD PMT
hit cluster in NHITAearly200 and the position of the OD PMT hit cluster in NHITA500 .
– NHITAearly200 ≥ 55
• Long tail ﬂasher cut.
The vertex position of the charged particle is reconstructed by the hit pattern of PMTs
assuming the light is emitted from a certain point in the ID (Point-ﬁt). While the goodnessof the Point-ﬁt tends to be small for ﬂasher events due to broad timing distribution by thedischarge of the PMT. The remaining ﬂasher events are rejected using this characteristic.
Events satisfying one of the following conditions are rejected.
– Point-ﬁt goodness < 0.3– Point-ﬁt goodness < 0.4 and NHITmin100 ≥ 6
The event rate in SK is reduced to about 30 events/day by these cuts.

6.6

Final fully contained dataset

After the data reduction, quality cuts are applied to select events for the proton decay search: theevents should be fully contained (FC) with the vertex position in the ﬁducial volume (FV) deﬁned
as a region more than 2 m away from the ID wall as shown in Figure 6.8 (FCFV selection). The
fully contained events are selected by the condition that the number of hit PMTs in the largest OD
hit cluster should be less than 16 hits. In addition, we require that the total visible energy (electron
equivalent total energy deposit in the detector) should be greater than 30 MeV (Evis selection).
Events away from the ID wall are selected as the analysis sample because events close to the ID
wall have fewer hit PMTs, resulting in the degradation of the event reconstruction performance.
Also, these events often contain ﬂasher events and radioactive backgrounds from the wall that could
not be removed by the reduction. Figure 6.9 shows the distribution of the vertex position of the data

51

Total volume: 50 kton
Outer detector:
18 kton
Fiducial volume:
22.5 kton

Inner detector:
32 ktonFigure 6.8: Schematic diagram of the ID, OD, and ﬁducial volume in SK.

Figure 6.9: Distributions of the radial (left) and central axial (right) positions (R and Z) of thedata after the ﬁfth reduction in SK-IV. The FC and Evis selection are applied for the data (black)and atmospheric neutrino simulation (red). For the R2 (Z) distribution, events are rejected if
the distance from the ID wall along the Z (R) direction is within 2 m. The atmospheric neutrino
simulation is normalized by the live time.
after the FCFV selection in SK-IV. The energy and positions are reconstructed from the observed
PMT signals by the ﬁTQun algorithm as explained in Chapter 7. In the FV region, the positions
are uniformly distributed as expected from the atmospheric neutrino events, and the backgroundcontamination from the events near the ID wall is rejected. After the all FC reduction, the event
rate is about 10 event/day. Signal eﬃciency after FCFV selection is estimated to be 97% from thesimulation. In the recent p → e+ + π 0 and p → µ+ + π 0 searches [8], an enlarged ﬁducial volume
cut (vertex more than 1 m away from the ID wall) was used for the proton decay search. It was not
employed in this analysis because the KL0 decays a few meters away from the point of the proton
decay as discussed in Section 4.1, and the KL0 decay point could be close to the ID wall or outside
the ID if the proton decays within 2 m of the ID wall.

52

Chapter 7

Event reconstruction

Figure 7.1: Event display of a typical neutrino event candidate. The solid line shows the reconstructed Cherenkov ring.

7.1

Reconstruction algorithm

Figure 7.1 shows an event display of a typical neutrino event candidate selected by the FC reduction. As explained in Section 3.1, hit patterns of PMTs are diﬀerent for particle properties. The
event reconstruction algorithm estimates basic event properties such as vertex (generated position),
direction, particle type and momentum based on the PMT hit information. As an example, theevent shown in Figure 7.1 was reconstructed as a muon with 484.7 MeV/c momentum. In the previously published SK-I to SK-III search [1], these properties are reconstructed by several steps (APFit
algorithm). For example, the vertex position is reconstructed ﬁrst by the timing information. Then,the number of rings, the particle type and momentum are reconstructed from the observed photoelectrons (p.e.’s) and the pattern of PMT hits in separate steps. In this study, the improved
event reconstruction algorithm ﬁTQun [68] is used. It is based on the maximum likelihood method
for both hit PMTs and unhit PMTs (PMTs with no signal). The vertex, momentum and particle
53

type of the ring are simultaneously determined from the time and hit information based on several
hit patterns created from the prior simulations. The performance of reconstruction is signiﬁcantly
improved with ﬁTQun [69].

7.1.1

Likelihood function

The event reconstruction maximizes the following likelihood function to determine the parameters.
L(x) =unhit
Y

Pj (unhit|x)

j

hit
Y

(1 − Pi (unhit|x))fq (qi |x)ft (ti |x),

(7.1)

i

where x is a set of reconstruction parameters (three parameters for vertex position, one parameterfor time, two parameters for zenith angle and azimuth angle, and one parameter for the absolutevalue of the momentum). Pj (unhit|x) is an unhit probability for each PMTs. fq (qi |x) and ft (qi |x)
are the probability density function of the observed number of p.e.’s and time, respectively.
Pj (unhit|x) is evaluated based on the predicted number of p.e.’s µ at parameter x. Since theobserved number of p.e.’s should follow a Poisson distribution with the predicted number of p.e.’sµ, the unhit probability Pj (unhit|x) follows exp(−µ). In addition, the probability is corrected up
to the third order to take into account the ineﬃciency due to the threshold of each PMT:Pj (unhit|x) = (1 + a1 µ + a2 µ2 + a3 µ3 ) exp(−µ).

(7.2)

The coeﬃcients (a1 , a2 , a3 ) are determined from the detector simulation.
The predicted number of p.e.’s for each PMT µi is estimated from the sum of expectations
sct
from the direct light µdir
(which arrives PMT photocathode after scattering ori , indirect light µi

reﬂection) and dark rate µdark . The dark rate is due to the PMT noise and the average rate foreach PMT was evaluated as 5.7195 kHz for SK-IV based on the measurements and implemented in
the simulation.
The prediction by the direct light from Cherenkov light emission is given as an integral along
the particle track length s:Zµdiri = Φ(p)dsg(p, s, cos θ)Ω(R)T (R)ϵ(η),(7.3)

where Φ(p) is a normalization factor proportional to the average total number of photons as afunction of initial momentum p. g(p, s, cos θ) is the distribution function of the number of emitted
photons which depends on the momentum p, track length s and angle to PMT θ. Ω(R), T (R) andϵ(η) are correction factors which depends on the distance to PMT R and the angle from PMT η.
Each factor will be explained below. Figure 7.2 shows the geometrical relations between the particle
and PMT.
g(p, s, cos θ) is called Cherenkov emission proﬁle and normalized asZg(p, s, cos θ)dsdΩ = 1,
54

(7.4)

Figure 7.2: Schematic diagram describing the variables relevant to calculate the predicted numberof p.e.’s. The white dot is the initial position of the particle. Taken from [71].
where Ω represents a solid angle. Figure 7.3 shows examples of the Cherenkov emission proﬁle forelectrons. Ω(R) is the correction factor of the PMT solid angle at a distance of R. The correction
to a circle of radius a = 25.4 cm (the radius of the PMT) from a location R away is assumed:Ω(R) =

πa2
.
R2

(7.5)

T (R) is the attenuation factor of direct light due to absorption and scattering in water:
T (R) = exp(−R/L),(7.6)

where L is the wavelength averaged attenuation length 7496.46 cm obtained from the detector
simulation with the absorption parameters determined by the calibration (Section 3.7). ϵ(η) is the
PMT angular acceptance obtained from the detector simulation (Figure 7.4).
The indirect light includes the light scattered in water and the reﬂected light coming fromdetector components such as the black sheet and surface of PMTs. The prediction of the indirect
light is similar to the one of the direct light:Zρ(p, s)
µsct
=
Φ(p)
ds
Ω(R)T (R)ϵ(η)A(s),i4πwhere ρ(p, s) is the integral of the Cherenkov emission proﬁle to the solid angle:Zρ(p, s) = g(p, s, cos θ)dΩ.

(7.7)

(7.8)

A(s) is a correction table for the contribution of scattered light obtained from the detector simulation
with the scattering and absorption parameters determined by the calibration (Section 3.7).

7.1.2

Vertex pre-ﬁtter

As the ﬁrst step of the event reconstruction algorithm, the generated position of the charged particles
is roughly determined by the vertex pre-ﬁtter. The vertex pre-ﬁtter is performed to start optimizing
55

Figure 7.3: Cherenkov emission proﬁle for electrons at diﬀerent initial momentum. Horizontal axes
represent the cosine of the angle from the particle direction, and vertical axes are the distancetraveled from the initial position of the particle. Taken from [71].

Figure 7.4: Angular acceptance of the PMT plotted as a function of cos η which is obtained fromdetector simulation with the reﬂection indices of water and PMT surface. A ﬁtted curve is overlaid
with a solid line. Taken from [71]. The vertical axis is shown in an arbitrary unit.

56

parameters from nearby the global minimum and avoid converging at a local minimum. Assuming
a point light source, the following likelihood function is used for the estimation.
" 
2 #
hitX
PMTs
i
Tres
G(x, t) =exp − √
,
2σ
i
where
i
Tres
= ti − t −

i
|RPMT
− x|
.
c/n

(7.9)

(7.10)

Here, x and t are the position and time of the vertex to ﬁt, σ is a hyperparameter of the ﬁt, ti is aihit timing of each PMT, and c/n is the light velocity in water. Tres
indicates the residual hit timeiof the i-th PMT. Tres
becomes smaller when x and t approach the true value, and G(x, t) decreasesiwith Tres. The pre-ﬁtter is applied to the hits around the trigger (−100, +400 ns). The grid size andσ in the G(x, t) are gradually shrunk, and the grid search in space and time is iteratively performed.

7.1.3

Hit clustering

If there are multiple hit clusters by the secondary reactions such as Michel electrons from muon
decay in the same time window, they are counted as subevents and separated by the hit clustering
algorithm.
First, the goodness, a return value of Equation 7.9, is evaluated for each 8 ns time window from−200 ns to 15,000 ns around the trigger with σ = 6.3 ns and ﬁxed vertex given by the pre-ﬁtter.
After ﬁnding the subevent candidate, the time window is set around the peak [−180 ns, +800 ns].
Finally, the vertex pre-ﬁtter and the peak search are conducted again for a precise peak ﬁnding.
The ﬁnal time windows are ﬁxed for each peak, and the maximum likelihood is calculated for eachof them. Figure 7.5 shows the goodness distribution as a function of time as an example of hit
clustering.

7.1.4

Single-ring ﬁtter

After the hit clustering, the single-ring ﬁtter is applied to each subevent. The negative log-likelihood
function (Equation 7.1) is minimized to obtain the optimal parameters such as vertex position andmomentum for the assumed particles e, µ, and π + which reproduce the observed PMT signals inthe detector. MINUIT [70] package is used to minimize the negative log-likelihood in the multiple
parameter space.

7.1.5

Multi-ring ﬁtter

More than one charged particle can be generated from proton decay and each makes a Cherenkov
ring. Figure 7.6 shows an example of such an event. Multi-ring ﬁtter is an algorithm to ﬁndsuch multiple Cherenkov rings. It is applied to only the ﬁrst subevent because multiple Cherenkov
rings are expected for the ﬁrst subevent, while the secondary particles are generated one by one at57

Figure 7.5: Distribution of the goodness as a function of time for an example event with a parent
muon and a Michel electron. The vertical dashed lines indicate the true generation time of theparticles in the simulation, the black points are the scanned goodness points, and the blue andgreen curves represent the threshold curves used for the subevent peak search. The red verticallines show the time of the subevent peaks found by the algorithm. Taken from [71].

Figure 7.6: Typical event display of multi-ring events.

58

diﬀerent timing and the subsequent subevents appear as single-ring events, such as Michel electron
events. The algorithm searches the rings by sequentially increasing the number of rings as shown in
Figure 7.7. The expected PMT hits of the single rings are added to the prediction and the negative
log-likelihood is evaluated in comparison with the observed PMT signals.

Figure 7.7: A tree diagram showing how the algorithm evolves the tree of the multi-ring hypotheses
for the increased number of rings. The diagram is for the case assuming the ﬁrst ring is an electron,and the equivalent procedure is developed for the case assuming π + (which is equivalent to the µhypothesis) as the ﬁrst ring. Taken from [71].
Additional rings are identiﬁed by the following procedure. First, the algorithm searches the
direction of the additional ring assuming an electron with 50 MeV/c momentum from the samevertex as the ﬁrst ring. Second, the momentum of the additional ring is adjusted to reproduce
the observed PMT signals while all the other parameters are ﬁxed. Third, the direction and themomentum of the additional ring are simultaneously optimized, while the vertex is ﬁxed to be
the same as the ﬁrst ring, and the direction and the momentum of the original rings are ﬁxed.
Fourth, the momenta of all rings are ﬁtted simultaneously while the vertex and their directions are
ﬁxed. Finally, the directions and momenta of all rings are simultaneously ﬁtted to get the ﬁnal
best-ﬁt result. This procedure is repeated until the sixth ring. To save the computation time, thesimultaneous ﬁt of all parameters (ﬁnal step) is performed only for the two-ring case. The best ﬁt
hypothesis for the number of rings and the combination of particle types is determined by choosing
the smallest negative log-likelihood.
After that, the best ﬁt hypothesis with the optimized parameters is corrected because it often
contains misidentiﬁed rings (fake rings). The angles between the directions of the most energetic
ring and the other rings are calculated, and the lower energy ring is merged with the highest energy
ring by adding the energy if the angle is below 20◦ . After that, all ring parameters are ﬁtted
again assuming three particle types: e, µ and π + . In this step, all vertices of the rings are ﬁtted
independently.

5950000

Number of events

Number of events

70000
40000
30000
20000
10000

60000
50000
40000
30000
20000
10000

0
0

500

1000

1500

2000

2500

0
0

3000

Momentum (MeV/c)

500

1000

1500

2000

2500

3000

Momentum (MeV/c)

Figure 7.8: Reconstructed momentum distributions in APFit (black) and ﬁTQun (red) for FCFV
single-ring events. The left (right) plot shows the e-like (µ-like) events for the νe , ν̄e (νµ , ν̄µ ) CCQE
sample.

7.2

Performance of ﬁTQun

Figure 7.8 shows the reconstructed momentum distributions in APFit (the algorithm used in theprevious search in SK) and ﬁTQun for FCFV (fully contained ﬁducial volume, as deﬁned in Section
6.6) single-ring events of the atmospheric neutrino simulation. Only charged current quasi-elastic
scattering (CCQE) sample is used for the distributions in which one charged lepton (e of µ) is
generated and most of them are expected as single-ring events. For some true lepton momentum
regions (200–400 MeV/c, 400–600 MeV/c and 600–1000 MeV/c), the momentum resolution (relative
error between the true and reconstructed momentum) is evaluated as shown in Figure 7.9. Twosided 68% intervals are calculated from them and summarized in Figure 7.10 as the momentum
resolution for each true lepton momentum. For e-like rings, the momentum resolution is clearly
improved with ﬁTQun in all momentum regions. For µ-like rings, the momentum resolutions ofboth algorithms are comparable but slightly improved with ﬁTQun in the 600–1000 MeV/c region.
Similarly, the vertex resolution (distance between true and reconstructed vertex) is evaluatedas shown in Figure 7.9 and 7.10. For both e-like and µ-like rings, the vertex resolution is clearly
improved with ﬁTQun in all momentum regions. As well as the momentum resolution, the improvement of the performance is clearer in the e-like rings compared to the µ-like rings.
π 0 meson decays into two γ’s with short lifetime (8.5 × 10−2 fs). Therefore π 0 production events
can be used to evaluate the performance of reconstruction for multi-ring events. The invariant
mass of π 0 was also reconstructed in APFit and ﬁTQun as shown in Figure 7.13. The distributions
consist of FCFV 2-ring e-like events for the neutral current single π 0 production (NC1π 0 ) sample
in the atmospheric neutrino simulation. The means and width of the distributions are evaluated by
Gaussian ﬁts as (137.2 ± 0.1) MeV/c2 and (13.9 ± 0.2) MeV/c2 for APFit and (135.9 ± 0.1) MeV/c2
and (13.5 ± 0.1) MeV/c2 for ﬁTQun, respectively. The width of the invariant mass distribution isslightly improved in ﬁTQun.

60

800

Number of events

Number of events

1000
800
600
400
200

700
600
500
400
300
200
100

0
− 0.5 −0.4 − 0.3 −0.2 − 0.1

0

0.1

0.2

0.3

0.4

0
−0.2 − 0.15 − 0.1 − 0.05

0.5

Momentum resolution

0

0.05

0.1

0.15

0.2

0.15

0.2

0.15

0.2

Momentum resolution
400

Number of events

Number of events

450
500
400
300
200
100

350
300
250
200
150
100
50

0
− 0.5 −0.4 − 0.3 −0.2 − 0.1

0

0.1

0.2

0.3

0.4

0
−0.2 − 0.15 − 0.1 − 0.05

0.5

Momentum resolution

0

0.05

0.1

Momentum resolution
400

500

Number of events

Number of events

450

400
300
200
100

350
300
250
200
150
100
50

0
− 0.5 −0.4 − 0.3 −0.2 − 0.1

0

0.1

0.2

0.3

0.4

0
−0.2 − 0.15 − 0.1 − 0.05

0.5

Momentum resolution

0

0.05

0.1

Momentum resolution

0.1
0.09
0.08
0.07
0.06
0.05
0.04
0.03
0.02
0.01
0

Momentum resolution

Momentum resolution

Figure 7.9: Relative error ((reconstructed − true)/true) of momentum in APFit (black) and ﬁTQun(red) for FCFV single-ring events. The left (right) plots show the e-like (µ-like) events for the νe , ν̄e(νµ , ν̄µ ) CCQE sample. The top, center and bottom plots are in the true e± or µ± momentum
regions of 200–400 MeV/c, 400–600 MeV/c and 600–1000 MeV/c, respectively.

200

300

400

500

600

700

800

900

1000

Momentum (MeV/c)

0.1
0.09
0.08
0.07
0.06
0.05
0.04
0.03
0.02
0.01
0

200

300

400

500

600

700

800

900

1000

Momentum (MeV/c)

Figure 7.10: Summary of the momentum resolution in APFit (black) and ﬁTQun (red) for FCFV
single-ring events. The left (right) plot shows the e-like (µ-like) events for the νe , ν̄e (νµ , ν̄µ ) CCQE
sample.

61

500

Number of events

Number of events

600
500
400
300
200
100
0
0

10

20

30

40

50

60

70

80

90

400
300
200
100
0
0

100

10

20

Vertex resolution (cm)

30

40

50

60

70

80

90

100

80

90

100

80

90

100

Vertex resolution (cm)300

Number of events

Number of events

400
250
200
150
100
50

350
300
250
200
150
100
50

0
0

10

20

30

40

50

60

70

80

90

0
0

100

10

20

30

40

50

60

70

Vertex resolution (cm)
450

240
220
200
180
160
140
120
100
80
60
40
20
0
0

400

Number of events

Number of events

Vertex resolution (cm)

350
300
250
200
150
100
50

10

20

30

40

50

60

70

80

90

0
0

100

10

20

Vertex resolution (cm)

30

40

50

60

70

Vertex resolution (cm)

50
45
40
35
30
25
20
15
10
5
0

Vertex resolution (cm)Vertex resolution (cm)Figure 7.11: Distance between the true and reconstructed vertex in APFit (black) and ﬁTQun(red) for FCFV single-ring events. The left (right) plots show the e-like (µ-like) events for the νe , ν̄e(νµ , ν̄µ ) CCQE sample. The top, center and bottom plots are in the true e± or µ± momentum
regions of 200–400 MeV/c, 400–600 MeV/c and 600–1000 MeV/c, respectively.

200

300

400

500

600

700

800

900

1000

Momentum (MeV/c)

50
45
40
35
30
25
20
15
10
5
0

200

300

400

500

600

700

800

900

1000

Momentum (MeV/c)

Figure 7.12: Summary of the vertex resolution in APFit (black) and ﬁTQun (red) for FCFV singlering events. The left (right) plot shows the e-like (µ-like) events for the νe , ν̄e (νµ , ν̄µ ) CCQE sample.

62

Arbitrary unit (area normalize)

0.08
0.07
0.06
0.05
0.04
0.03
0.02
0.01
0
0

50

100

150

200

250

300

2

Invariant mass (MeV/c )Figure 7.13: Reconstructed π 0 mass distributions of APFit (black) and ﬁTQun (red). Events are
FCFV 2-ring e-like of the NC1π 0 sample in the atmospheric neutrino simulation.

𝜋±
𝑝
𝜇

!

𝑙∓

𝐾"#
𝜈

Figure 7.14: Schematic diagram of p → µ+ + KL0 , (KL0 → π ± + l∓ + ν) channel.

7.3

Displaced vertex ﬁter

As explained in Section 4.1, a muon and a neutral kaon are produced at the same vertex from proton
decay. Among these, the muon immediately emits Cherenkov light. Neutral kaon is a superpositionof KS0 and KL0 . KS0 ’s decays at almost the same position as the proton decay. On the other hand,
KL0 ’s undergo hadronic interactions while propagating in the water, and decay about 1 m away witha lifetime of 51 ns. Figure 7.14 shows the schematic diagram of the process with KL0 . In ﬁTQun,
multiple rings are searched for assuming a single vertex in the initial step of the multi-ring ﬁt.
This leads to the degradation of the reconstruction performance for charged particles generated at
the displaced vertex. For example, Figure 7.15 shows the comparison of the number of rings forp → µ+ + KL0 events in the simulation with the nominal lifetime of 51 ns for KL0 decay (black) andthe intentionally shortened lifetime of 90 ps for KL0 decay (red). Considering the fact that morethan one particle is generated from p → µ+ + KL0 and following KL0 decay, a fraction of single-ring
events tells the misreconstruction of KL0 decay (as a muon creates clear ring which is reconstructed
as the ﬁrst ring for most events). A large fraction of single-ring events is reconstructed as multi-ringif the lifetime is shorter.

63

Arbitrary unit (area normalize)

0.4
0.3
0.2
0.1
0

0

2

4

6

Number of rings
Figure 7.15: Number of rings reconstructed by ﬁTQun. p → µ+ + KL0 simulation with nominal KL
lifetime (black) and KL lifetime of 90 ps (red) are shown.

BeforeAfter

1. Reconstruct the first ring

1. Reconstruct the first ring
2. Search second vertex
3. Search the second ring assuming
the second vertex.
4. Search the third ring assuming
the second vertex.

2. Search the second ring assuming
the primary vertex.
3. Search the third ring assuming
the primary vertex.

Figure 7.16: Schematic diagram of multi-ring reconstruction before and after the modiﬁcation.
For this reason, the multi-ring reconstruction algorithm in ﬁTQun was modiﬁed for p → µ+ +KL0
event selection to reconstruct a muon from proton decay and particles from KL0 decay. Hereafter,
the modiﬁed algorithm is called ”displaced vertex ﬁtter”.
Figure 7.16 shows the procedure of multi-ring reconstruction before and after the modiﬁcation.
After the ﬁrst ring reconstruction, the hit PMTs which belong to the ﬁrst ring are masked and thevertex pre-ﬁt is performed again to search the secondary vertex using the remaining hits. Masked
PMTs are not used for the calculation of the likelihood. Once the secondary vertex candidate
is determined, PMTs are unmasked and the reconstruction is performed to ﬁnd the third rings.
The direction of the additional ring is searched in grid space assuming an electron with 50 MeV/cmomentum generated at the second vertex determined in the previous step. After that, as well as
the multi-ring ﬁt, the parameters are ﬁtted in steps with the independent second vertex. Whenﬁtting the subsequent rings, the search is performed assuming the same vertices as the second ring.

64

7.4

Performance of displaced vertex ﬁtter

To check the performance of the displaced vertex ﬁtter, it was applied to the simulation of p →µ+ + KL0 (→ π ± + e∓ + ν) events with π ± momentum below the Cherenkov threshold. These eventsshould be identiﬁed as two rings for µ+ and e∓ from the decay of KL0 because the π ± does not emit
the Cherenkov light. The number of ring distribution is shown in Figure 7.17. This plot showshigher eﬃciency to ﬁnd the second ring by the displaced vertex ﬁtter. Then, the performance ofthe vertex reconstruction is evaluated by the 2-ring sample as shown in Figure 7.18. It shows the
distance between the true primary (secondary) vertex and the reconstructed primary (secondary)vertex. It was conﬁrmed that the most of primary and secondary vertices are reconstructed within

Number of events

50 cm of the true vertex with the displaced vertex ﬁtter.
800
600
400
200
0

0

2

4

6

Number of rings

150

150

Number of events

Number of events

Figure 7.17: Number of rings reconstructed by original multi-ring ﬁtter (black) and displaced vertex
ﬁtter (red) for p → µ+ KL (→ π ± e∓ ν) simulation with the π ± momentum below the Cherenkov
threshold.

100

50

0
0

20

40

60

80

100

50

0
0

100

Distance (cm)

20

40

60

80

100

Distance (cm)Figure 7.18: Left (right) shows the distance between the true primary (secondary) vertex and thereconstructed primary (secondary) vertex.
Decay electrons from muons are often reconstructed as the second vertex although these are not
the target of the search for the kaon decay. To remove the contamination of decay electrons, rings
belonging to the second vertex are discarded if the diﬀerence in the reconstructed times between
65the primary and secondary candidates is more than 100 ns.

7.5

Neutron tagging

Neutrons generated in water are captured by the surrounding hydrogen (n + p → d + γ) after about
200 µs on average. Neutrons can be identiﬁed by detecting the 2.2 MeV γ emitted by this capture.
Until SK-III running period, only the hits within 1.3 µs around the trigger were recorded, so theneutrons could not be identiﬁed. Since SK-IV running period, the electronics modules were updated
to record all PMT hits from 35 µs to 535 µs after the primary trigger as described in Section 3.6.
This enables the detection of neutron capture on hydrogen in that time window. The detection
eﬃciency of neutrons is about 20% in the SK-IV period because the 2.2 MeV γ are easily buried by
the dark hits of PMT and environmental radioactivity which mainly produces events below 4 MeV[30]. The signals from neutron captures by hydrogens are identiﬁed based on the algorithm of theneural network [72]. Despite the relatively low eﬃciency, the information on the number of detected(tagged) neutrons is useful to improve the sensitivity of proton decay searches.
Neutrons are often generated by atmospheric neutrino interactions and by the secondary interactions of the hadrons in water, especially at relatively high energy above the GeV scale, which are
the main component of the background in the proton decay search. On the other hand, neutrons
are generated with about 10% probability after proton decay in oxygen due to de-excitation of15

N [73]. As a result, neutrino interactions tend to have tagged neutrons while it is not the casefor proton decay. Therefore, the background events are suppressed by requiring no neutron in the
proton decay search as written in the next section.

66

Chapter 8

Event selection
As proton decay is an extremely rare process, even if it happens, the signal events are overwhelmedby the atmospheric neutrino backgrounds. To further reduce atmospheric neutrino events after
the data reduction as described in Chapter 6, selection criteria were determined based on thespeciﬁcations. Selection criteria were ﬁxed before they are applied to the data to avoid biases inthe deﬁnition of the cuts and the threshold values. To suppress the statistical ﬂuctuation of theremaining backgrounds, background events must be reduced to O(1) events or below with the tightcut. The detection eﬃciencies of proton decay events and the number of expected background eventswere estimated by applying the selection criteria to the proton decay simulation and atmospheric
neutrino simulation.

8.1

Overview of selection

There are ﬁve selections to extract each K 0 decay mode in p → µ+ + K 0 search.
A: KS0 → 2π 0
B: KS0 → π + + π − 0KL → π ± + l ∓ + ν l





KL0 → 3π 0
C:




 0
KL → π + + π − + π 0

where l is an electron or muon

Here, A, B and C are the labels of the selections. The three selections for KL0 are grouped as C andbranch oﬀ in the middle of the selection criteria. The ﬁnal samples of selection A (A and B) areexcluded in selection B (C) to avoid double counting of events. Figure 8.1 shows the ﬂow diagramof the event selection.
Figure 8.2 shows the event display of the proton decay simulation for p → µ+ + K 0 , KS0 → 2π 0(left) and p → µ+ + K 0 , KL0 → 3π 0 (right). For both KS0 and KL0 decay modes, the existenceof monoenergetic primary muon from proton decay is required in the selection. In the search for67

Events after data reduction
Criteria A
Criteria B
Criteria C
Criteria C
Criteria C

Final sample for 𝐾!" → 2𝜋 "Final sample for 𝐾!" → 𝜋 # 𝜋 $Final sample for 𝐾%" → 𝜋 ± 𝑙 ∓ 𝜈
Final sample for 𝐾%" → 3𝜋 "Final sample for 𝐾%" → 𝜋 # 𝜋 $ 𝜋 "Figure 8.1: Flow diagram of the proton decay selection.

Figure 8.2: Event displays of the proton decay simulation for p → µ+ + K 0 , KS0 → 2π 0 (left) andp → µ+ + K 0 , KL0 → 3π 0 (right). The MC true rings of the muon (red) the γ’s from π 0 decay(blue) and are also shown.

68

KS0 decay mode, background events are largely reduced by requiring the reconstructed invariant
masses are consistent with kaon mass and proton mass. On the other hand, invariant mass is notreconstructed in the search of KL0 decay as some of the outgoing particles are invisible in the water
Cherenkov detector. Therefore, instead of the invariant mass cut, the distance between the primaryand secondary vertices is reconstructed by the displaced vertex ﬁtter and used in the selection to
suppress the background. Details of each selection are described in the following sections.
There are two major updates to the selection criteria from the previous p → µ+ + K 0 search [1].
In the previous study, events with two and three Cherenkov rings were selected as candidates for
the KS0 → π + + π − decay mode, in which three charged particles (µ+ , π + and π − ) are generated.
While only three-ring events are selected in this study, as the background rate in two-ring events ishigher than the other selections and the sensitivity from the selection is not signiﬁcant comparedto the others. As the other update, three selection criteria are applied for each KL0 decay mode
instead of a single selection for all KL0 decay used in the previous analysis.

8.2

Selection for p → µ+ + KS0 , KS0 → 2π 0

The events which satisfy the following criteria (A-1 – A-9) are selected as candidates for KS0 → 2π 0 .
Each π 0 decay into two γ’s, so the primary muon and four γ’s appear in this decay mode.
A-1: Events should pass FCFV and Evis selection.
As explained in Section 6.6, we require that events are fully contained (FC) deﬁned as there
is no cluster of hit PMTs in the outer detector (OD), the vertex is in the ﬁducial volume (FV)deﬁned as a region more than 2 m away from the inner detector (ID) wall, and the electron
equivalent total energy deposit in the detector, Evis , should be greater than 30 MeV. These
cuts are applied to avoid the degradation of the event reconstruction performance and toreject remaining ﬂasher events and radioactive backgrounds from the wall.
A-2: The number of rings should be three, four or ﬁve.
The top left plot in Figure 8.3 shows the number of rings after selection A-1 is applied. One
muon from the primary proton decay and four γ’s from the decay of two π 0 ’s (π 0 → 2π 0 )appear in this decay channel. Two γ’s from π 0 decay make overlapped showering rings, andthey are often reconstructed as one ring with the energy close to the sum of two γ’s. Or,
if one of two γ’s was small energy, it cannot be recognized as a ring but the eﬀect on thereconstruction of the invariant mass is limited. For these reasons, events with three, four andﬁve rings are kept as candidates for this decay channel.
A-3: There must be one non-showering ring.

69

300

2

200

100

0

200

400

600

Muon momentum (MeV/c)

Number of events

Number of events

Number of events

80

60

40

20

0
0

200

400

600

800 1000

K0 invariant mass (MeV/c2)

1.5

1

0.5

0

0

2

4

6

Number of neutrons

Figure 8.3: From top left to bottom right, reconstructed number of rings after selection A-1, numberof non-showering rings after selection A-2, number of Michel electrons after selection A-3, muon
momentum after selection A-4, kaon invariant mass after selection A-5, and number of tagged
neutrons after selection A-8. The data (black dots) are compared with the atmospheric neutrino
MC events (red) normalized by the live time, and the signal MC events (blue) normalized to theatmospheric neutrino MC events.

70

The non-showering ring is assumed to be due to the primary muon. The top center plot in
Figure 8.3 shows the number of non-showering rings (µ± or π ± like) after selection A-2.
A-4: There must be one Michel electron.
The top right plot in Figure 8.3 shows the number of Michel electrons after selection A-3.
This cut requires the existence of a Michel electron from the decay of the primary muon.
A-5: The reconstructed momentum of the non-showering ring should be 150 < Pµ < 400 MeV/c.
The bottom left plot in Figure 8.3 shows the momentum of the non-showering ring (muon)
after selection A-4. Muon momentum from p → µ+ + K 0 decay is expected to peak at
326.5 MeV/c since it is two-body decay. It should be monochromatic momentum if proton
decay occurs in hydrogen, while it slightly ﬂuctuates if the proton decay occurs in oxygen
nuclei due to the Fermi motion of protons.
A-6: The reconstructed invariant mass of the showering rings should be 400 < MK < 600 MeV/c2 .
The bottom center plot in Figure 8.3 shows the reconstructed invariant mass of the showeringrings after selection A-5. γ’s from π 0 decay make showering rings, and the invariant mass ofthese rings is expected to be 497.6 MeV/c2 , equivalent to K 0 mass. This cut selects the events
around the peak of K 0 mass. A peak around 150 MeV/c2 in the invariant mass distribution
is due to π 0 either from atmospheric neutrino interactions (red) or from the decay of kaons
from proton decay (blue).
A-7: The reconstructed total momentum should be Ptot < 300 MeV/c.
A-8: The reconstructed total invariant mass should be 800 < Mtot < 1050 MeV/c2 .
Figure 8.4 shows scatter plots of the reconstructed total invariant mass and total momentum
after applying all the cuts except those on the plotted variables (the neutron tagging cut,
A-9, is also applied in this plot). The black box shows the region selected by both cuts A-7
and A-8. The total momentum and total invariant mass are reconstructed from all rings andtherefore assumed to be the proton momentum and mass. In the proton decay simulation,
the momentum of the proton in hydrogen is close to zero, and the invariant mass has a peakaround 938 MeV/c2 , equivalent to proton mass. In contrast, protons have Fermi momentumin oxygen nuclei, and the total momentum and invariant mass of the proton decay simulation
have broad distribution. Free proton events in the signal MC simulation which appear aroundthe top left corner of the box are due to secondary KS0 decay generated by KL0 scattering. As
shown in the right-hand plot in Figure 8.4, no candidates remained in this sample.
A-9: There should be no tagged neutrons.

71

Total momentum (MeV/c)

1000
800
600
400
200
0

400

600

800

1000

400

600

800

1000

400

600

800

1000 1200

Total invariant mass (MeV/c2)

Figure 8.4: Scatter plot of the reconstructed total invariant mass and total momentum after
applying all cuts in KS0 → 2π 0 selection except those on the plotted variables. From left to right,
signal MC events, atmospheric neutrino MC events (500 years-equivalent) and data (3244.39 livedays) are shown. In the signal MC events, cyan shows free protons and blue shows bound protons.
The bottom right plot in Figure 8.3 shows the number of tagged neutrons after selection A-8.
For this decay mode, 93% of the signal events have no neutrons, while background events
tend to have tagged neutrons. A candidate event remaining after selection A-8 was rejectedby requiring no tagged neutrons.

8.3

Selection for p → µ+ + KS0 , KS0 → π + + π −Both π + and π − make non-showering Cherenkov rings. As the momentum of these secondary
muons from π ± decay are below the Cherenkov threshold, they are not observed in the detector.
In this selection, three non-showering rings are selected, and the most energetic ring is recognized
as the muon from primary proton decay. The other rings are taken to be charged pions from KS0
decay. The events which satisfy the following criteria (B-1 – B-7) are selected as candidates forKS0 → π + + π − .
B-1: Events should pass FCFV and Evis selection.
B-2: The number of rings should be three.
B-3: All rings should be non-showering.
The top left and top right plots in Figure 8.5 show the number of rings after selection B-1 andthe number of non-showering rings after selection B-2, respectively. These rings are assumed
to be due to a muon from proton decay and charged pions from KS0 decay.
B-4: The number of Michel electrons should be one or two.
The bottom left plot in Figure 8.5 shows the number of Michel electrons after selection B-3.
π + generates a muon from the decay as π + → µνµ . On the other hand, π − from KS0 decay
72

Number of events

40

20

0

400

600

800

K0 invariant mass (MeV/c2)

Figure 8.5: Reconstructed invariant mass distribution for kaons after selection B-4. Reconstructed
number of rings after selection B-1, number of non-showering rings after selection B-2, number ofMichel electrons after selection B-3, kaon invariant mass after selection B-4. The data (black dots)
are compared with the atmospheric neutrino MC events (red) normalized by the live time, and thesignal MC events (blue) normalized to the atmospheric neutrino MC events.

73

is often captured on a 16 O nucleus and does not generate a muon. Therefore, the number oftagged Michel electrons in the signal event is expected to be one or two, one from the primarymuon and the other from π + . Considering the ineﬃciency of Michel electrons, events with
only one Michel electron are kept.
B-5: The reconstructed invariant mass of the second and third energetic non-showering rings shouldbe 450 < MK < 550 MeV/c2 .
The bottom right plot in Figure 8.5 shows the reconstructed invariant mass of the second
and third energetic non-showering rings after selection B-4. A peak at around 500 MeV/c2
corresponds to the K 0 mass (497.6 MeV/c2 ). The lower edge of the kaon invariant mass
distribution is due to the Cherenkov threshold of charged pions.
B-6: The reconstructed total momentum should be Ptot < 300 MeV/c.
B-7: The reconstructed total invariant mass should be 800 < Mtot < 1050 MeV/c2 .
Figure 8.6 shows scatter plots of the reconstructed total invariant mass and total momentum
after applying all cuts except those on the plotted variables. No candidates remained in thissample.

Total momentum (MeV/c)

1000
800
600
400
200
0

400

600

800

1000

400

600

800

1000

400

600

800

1000 1200

Total invariant mass (MeV/c2)

Figure 8.6: Scatter plot of the reconstructed total invariant mass and total momentum after
applying all cuts in KS0 → π + π − selection except those on the plotted variables. From left to right,
signal MC events, atmospheric neutrino MC events (500 years-equivalent) and data (3244.39 livedays) are shown. In signal MC events, cyan shows free protons and blue shows bound protons.
Criterion with the number of tagged neutrons is not applied because about 45% of signal events
have tagged neutrons from pion capture.

8.4

Selections for p → µ+ + KL0

There are three selection criteria for KL0 decay, that is KL0 → π ± + l∓ + ν, KL0 → 3π 0 and
KL0 → π + + π − + π 0 . l represents charged leptons either electron or muon. In the KL0 decay
selections (C), the primary non-showering ring is chosen as the primary muon from proton decay.
74

Number of events

10000

5000

0

0

2

4

6

Number of rings

Figure 8.7: The number of p.e.’s after selection C-1 and the number of rings after selection C-2.
The data (black dots) are compared with the atmospheric neutrino MC events (red) normalized by
the live time, and the signal MC events (blue) normalized to the atmospheric neutrino MC events.
C-1: Events should pass FCFV and Evis selection.
C-2: Total observed photoelectrons (p.e.’s) should be 500 < Qtot < 8000 p.e.
Since KL0 decay is separated from the proton decay in time, and rings tend to overlap withsmaller momentum in three-body decay, not all of the rings can be always reconstructed.
Even if some rings are not reconstructed, the signal events can be roughly distinguished fromthe background events by the total observed p.e.’s. The left-hand plot in Figure 8.7 shows
the number of p.e.’s after selection C-1.
C-3:
C-3-1: The number of rings should be two or three (for KL0 → π ± l∓ ν).
C-3-2: The number of rings should be four, ﬁve or six (for KL0 → 3π 0 ).
C-3-3: The number of rings should be three or four (for KL0 → π + π − π 0 ).
The right-hand plot in Figure 8.7 shows the number of rings after selection C-2. As explained
above, not all charged particles are reconstructed with the Cherenkov ring. Therefore, loose
cuts are applied for each decay mode.
C-4:
C-4-1: The number of showering rings should be zero or one (for KL0 → π ± l∓ ν).
The left plot in Figure 8.8 shows the number of showering rings after selection C-3-1.
The rings are assumed to be the primary muon and charged pion and lepton from KL0
decay. If charged pions have momentum below the Cherenkov threshold, it cannot be

75

observed and the number of rings becomes two. One (zero) showering ring is observedfor KL0 → π ± e∓ (µ∓ ) ν channel.
C-4-2: The number of non-showering rings should be one (for KL0 → 3π 0 ).
The center plot in Figure 8.8 shows the number of non-showering rings after selection
C-3-2. The rings are assumed to be the primary muon and γ’s from π 0 decay.
C-4-3: The number of showering rings should be two (for KL0 → π + π − π 0 ).
The right plot in Figure 8.8 shows the number of showering rings after selection C-3-3. In
this decay mode, ﬁve particles are generated as candidates of Cherenkov rings (primarymuon, two charged pions and two γ’s from π 0 decay). Among those, at least one ofthe pions is invisible due to momentum below the Cherenkov threshold. Four rings are
assumed to be due to the primary muon, two γ’s from π 0 decay and one of the charged
pions from KL0 decay. If both charged pions are invisible, the number of rings becomes
three.

Figure 8.8: The number of showering rings or non-showering rings after C-3 in each selection labeledin the plots. The data (black dots) are compared with the atmospheric neutrino MC events (red)normalized by the live time, and the signal MC events (blue) normalized to the atmospheric neutrino
MC events.

C-5:
C-5-1: The number of Michel electrons should be two or three (for KL0 → π ± l∓ ν).
Three Michel electrons are generated as the secondary particles from the decay of the
primary muon and the pion (and muon for KL0 → π ± µ∓ ν) from KL0 decay.
C-5-2: The number of Michel electrons should be one (for KL0 → 3π 0 ).
Only one Michel electron is generated as the secondary particle from the decay of the
primary muon.
C-5-3: The number of Michel electrons should be two or three (for KL0 → π + π − π 0 ).

76

Three Michel electrons are generated as the secondary particles from the decay of the
primary muon and two charged pions from KL0 decay. If π − is captured, the number ofMichel electrons becomes two.
The three plots in Figure 8.9 show the number of Michel electrons after selection C-4-1 (left),
C-4-2 (center) and C-4-3 (right).

Figure 8.9: The number of Michel electrons after C-4 in each selection labeled in the plots. The
data (black dots) are compared with the atmospheric neutrino MC events (red) normalized by thelive time, and the signal MC events (blue) normalized to the atmospheric neutrino MC events.

C-6: The reconstructed muon momentum should be 260 < Pµ < 410 MeV/c
This cut selects the primary muon from p → µ+ +K 0 decay with the momentum of 326.5 MeV/c2 .
Figure 8.10 shows the reconstructed muon momentum after C-5 selection. The combined muon
momentum distribution for all KL0 decay selections is also shown in Figure 8.11.

Figure 8.10: Reconstructed muon momentum after C-5 in KL0 → π ± l∓ ν (left), KL0 → 3π 0 (center)
and KL0 → π + π − π 0 (right) selections. The data (black dots) are compared with the atmospheric
neutrino MC events (red) normalized by the live time, and the signal MC events (blue) normalized
to the atmospheric neutrino MC events.

C-7: The reconstructed vertex separation should be 1.5 m < vsep (vsep is deﬁned in the text below).

77

Number of events

80
60
40
20
0

200
400
Momentum (MeV/c)

600

Figure 8.11: Reconstructed muon momentum of all events after C-5-1, C-5-2 and C-5-3 selections.
The data (black dots) are compared with the atmospheric neutrino MC events (red) normalized by
the live time, and the signal MC events(blue) normalized to the atmospheric neutrino MC events.
Since these KL0 decays are three-body decays with relatively long lifetimes, it is diﬃcult toreconstruct all secondary particles. Therefore, strict cuts to the invariant mass are not appliedin these selections. Instead, the vertex separation vsep is deﬁned as the distance between the
primary (muon) and secondary (kaon) vertices along the opposite direction of the primarymuon as illustrated in Figure 8.12, and a cut is applied on it to distinguish the signal eventsfrom the backgrounds. Typical p → µ+ + KL0 events have positive vertex separation due toopposite directions of muon and kaon from proton decay, while the average vertex separationof atmospheric neutrino events is zero.

Figure 8.12: A schematic view of vertex separation.
Figure 8.13 shows vertex separation vsep distributions after applying all cuts including C-8
but except for C-7. The combined vertex separation distribution for all KL0 decay selections
is also shown in Figure 8.14. The right side region is the ﬁnal sample of KL0 decay selection.
The peak positions between the signal and atmospheric neutrino MC simulation diﬀer by
about 0.5 m. The distribution of the atmospheric neutrino MC simulation has a larger tailin the negative region. This is mainly due to Michel electrons from muons produced by νµ
CC interaction and the scattering of charged pions in water which causes multiple Cherenkov
rings. The value of the vertex separation cut was chosen to optimize for the sensitivity of the78

proton decay search. One candidate remained in KL0 → π ± l∓ ν selection.

Number of events

Figure 8.13: Reconstructed vertex separation in KL0 → π ± l∓ ν (left), KL0 → 3π 0 (center) andKL0 → π + π − π 0 (right) selections. All cuts except vertex separation are applied. The data (blackdots) are compared with the atmospheric neutrino MC events (red) normalized by the live time,
and the signal MC events (blue) normalized to the atmospheric neutrino MC events.

30
20
10
0

−5

0
5
Vertex separation (m)Figure 8.14: Reconstructed vertex separation distribution of all events in KL0 → π ± l∓ ν, KL0 → 3π 0
and KL0 → π + π − π 0 selections. All cuts except vertex separation are applied. The data (black dots)
are compared with the atmospheric neutrino MC events (red) normalized by the live time, and thesignal MC events (blue) normalized to the atmospheric neutrino MC events.
C-8: There should be no tagged neutrons.
The three plots in Figure 8.15 show the number of tagged neutrons after selection C-7. About
51%, 73% and 71% of the remaining backgrounds are rejected by requiring no tagged neutrons
while 38%, 26% and 42% of the signal events remain for KL0 → π ± l∓ ν, KL0 → 3π 0 and
KL0 → π + π − π 0 selections, respectively.

8.5

Selection summary

Table 8.1 shows a summary of the event selections. It shows the signal eﬃciency, the number ofbackground events and the number of candidate events after all cuts. The eﬃciency of each selection
is estimated from the proton decay simulation with a deﬁnition as a fraction of the number of signal
events after selection to the number of generated events (including all decay modes of K 0 ) within the79

Figure 8.15: The number of tagged neutrons after C-7 in each selection labeled in the plots. The
data (black dots) are compared with the atmospheric neutrino MC events (red) normalized by thelive time, and the signal MC events (blue) normalized to the atmospheric neutrino MC events.
ID more than 2 m away from the ID wall. The total eﬃciency for p → µ+ + K 0 decay is 17.0±1.2%
and the total number of expected backgrounds is 15.5 ± 2.9 events/(Mton·year) where the unit of
Mton·year means a statistics equivalent to one-year observation with one Mton mass. In the SKIV period, 3.1 ± 0.6 background events are expected in total for 3244.39 live days (0.2 Mton·years
exposure of the Super-Kamiokande detector with 22.5 kton ﬁducial mass). Eﬃciencies for KL0 decay
are lower than those for KS0 decay due to the cut on the vertex separation, which is necessary to
suppress the background. Systematic uncertainties are also shown for the signal eﬃciencies and thebackground rates. The details of the systematic uncertainties are explained in Chapter 10. As aresult of the selections, one candidate remains in the ﬁnal samples for KL0 → π ± l∓ ν selection, while
no candidates remained in the ﬁnal sample of the other selections. As the expected background
events for KL0 → π ± l∓ ν and total channels are 1.7 ± 0.5 and 3.1 ± 0.6, respectively, these numberof candidates are consistent with the background. As there are no signiﬁcant excesses beyond the
expected backgrounds, limits of the proton lifetime are set in Chapter 9.
Figure 8.16 shows the signal eﬃciency, the number of background events and the number ofcandidates for each step of the selections. For KS0 decay modes, about 95% of background events
are rejected by the total (proton) momentum cut while keeping about 90% of the signal events. For

Table 8.1: Summary of the p → µ+ + K 0 search in the SK-IV period. Uncertainties are quadraticsums of the MC statistical uncertainties and the systematic uncertainties.
Search mode (branching ratio)
KS0 → 2π 0(30.7% for KS0 )0+ −KS → π π(69.2% for KS0 )0± ∓KL → π l ν(67.6% for KL0 )00KL → 3π(19.5% for KL0 )0+ − 0
KL → π π π(12.5% for KL0 )Total

Eﬃciency (%)
9.9 ± 1.0
5.2 ± 0.6
1.4 ± 0.3
0.37 ± 0.05
0.18 ± 0.04
17.0 ± 1.2

80

Background (events)
0.3 ± 0.1
0.8 ± 0.2
1.7 ± 0.5
0.12 ± 0.06
0.16 ± 0.07
3.1 ± 0.6

Candidates
0
0
1
0
0
1

KL0 decay modes, more than 95% of the remaining background events are rejected by the vertex
separation cut whereas the signal eﬃciencies are about 10%.
K S →π+π-±

K S →2π0

K L→π±l νl

K L→3π0

K L→π+π-π0

Efficiency (%)

102

10

1

Number of events

10−1
104

102

1

10−2 A1 A2 A3 A4 A5 A6 A7 A8 A9 B1 B2 B3 B4 B5 B6 B7 C1 C2 C3 C4 C5 C6 C7 C8 C1 C2 C3 C4 C5 C6 C7 C8 C1 C2 C3 C4 C5 C6 C7 C8Selection criteria

Figure 8.16: Signal eﬃciencies (blue), number of background events (red) and number of candidates(black) for each step of the selections for SK-IV. Error bars show statistical uncertainty. The numberof atmospheric neutrino MC events is normalized by the live time (3244.39 live days).
Table 8.2 shows the breakdown of the remaining background events by the interaction modeof the atmospheric neutrinos (an explanation of neutrino interactions is written in Section 5.3).
The major background source is not kaon production but single pion and eta production, anddeep inelastic scattering DIS with multiple pions, due to much larger cross sections than the kaon
production.
Table 8.2: Breakdown of remaining atmospheric neutrino backgrounds by the interaction mode(%). CC, NC, QE and DIS stand for charged-current, neutral-current, quasi-elastic and deep
inelastic scattering, respectively. Uncertainties are statistical. ”-” means no events remained in thesimulation corresponding to the statistics of the 500 years.
Modes
CCQE
CC1π
CC1K
CC1η
CCDIS
NC

KS0 → 2π 0
1.9 ± 1.9
23.4 ± 9.3
37.3 ± 10.6
31.6 ± 11.5
5.8 ± 5.6

KS0 → π + π −70.2 ± 6.3
3.4 ± 2.5
11.3 ± 4.2
15.0 ± 5.5

KL0 → π ± l∓ ν
9.7 ± 2.8
82.9 ± 3.6
3.2 ± 1.6
4.2 ± 2.081

KL0 → 3π 08.2 ± 8.1
28.9 ± 15.4
48.6 ± 17.7
14.2 ± 13.0

KL0 → π + π − π 057.3 ± 16.0
42.7 ± 16.0-

Chapter 9

Lifetime limit
Since no statistically signiﬁcant excesses were observed as described in Section 8.5, a lower bound on
the proton lifetime was calculated by using a Bayesian method [78][2]. The probability distributionfunction of the decay width of proton decay is expressed by a Poisson distribution convolved withthe systematic uncertainties as follows:
ZZZ
P (Γ|ni ) =e−(Γλi ϵi +bi ) (Γλi ϵi + bi )ni
P (Γ)P (λi )P (ϵi )P (bi )dϵi dλi dbi ,
ni !

(9.1)

where i is an index for each selection. Γ is the decay rate which is a reciprocal of the lifetime. ni is
the number of observed events, λi is the exposure (proportional to the live time of the observation),
ϵi is the signal eﬃciency and bi is the number of expected background events. The probability
distribution function of decay rate P (Γ) is assumed to be uniform as there is no indication of its
value. P (λi ), P (ϵi ) and P (bi ) are the probability distribution functions for the exposure, eﬃciency
and backgrounds, respectively:
P (λi ) ∝
exp −(λi −µλi )2(λi > 0)(otherwise)

2
2σλ

i

0



exp −(ϵi −µϵi )222σϵi
P (ϵi ) ∝0(ϵi > 0)(otherwise)R ∞ e−B B nbi exp −(Ci bi −B)2 dB (b > 0)i
nb i !
0
2σb2
P (bi ) ∝i
0(otherwise)

(9.2)

(9.3)

(9.4)

where σλi , σϵi and σbi are the systematic uncertainties and µλi and µϵi are the expected exposureand eﬃciency, respectively. nbi is the number of expected backgrounds without live time normalization and Ci is a factor to normalize the MC events to the live time of the data. The proton
decay rate Γlimit at a 90% C.L. satisﬁes the following equation using this probability distributionfunction,

Z Γlimit
0.9 =

dΓ
0

Y
i

82

P (Γ|ni ).

(9.5)

The results of all selections are combined in this equation and the decay rate is evaluated. The
lower limit of the proton lifetime is expressed as the inverse of the decay rate:τlimit /Br = 1/Γlimit ,(9.6)

where Br is the branching ratio of the proton decay mode.
The lifetime limits by this method are summarized in Table 9.1 with the eﬃciencies, expected
number of backgrounds and number of candidates for each selection. The limit from SK-IV data is
4.5 × 1033 years at 90% C.L.
Table 9.1: Summary of the lifetime limits in p → µ+ K 0 search. The signal eﬃciency, number ofbackground events and number of candidate events for each selection are taken from Table 8.1.
Search modeEﬃciency (%) Background (events)
KS0 → 2π 09.9 ± 1.0
0.3 ± 0.1
KS0 → π + π −5.2 ± 0.6
0.8 ± 0.2
1.4 ± 0.3
1.7 ± 0.5
KL0 → π ± l∓ ν
0.37 ± 0.05
0.12 ± 0.06
KL0 → 3π 00.18 ± 0.04
0.16 ± 0.07
KL0 → π + π − π 0
SK-IV combined (199 kton·years)
SK-I+SK-II+SK-III+SK-IV combined (372 kton·years)

Candidates
0
0
1
0
0

Lower limit (10 33 years)
2.8
1.5
0.3
0.1
0.05
4.5
3.6

Then this result from SK-IV data (2008–2018) was combined with the previous result using thestatistically independent data in SK-I to SK-III period (1996–2008) [1] and the limit from SK-I to
SK-IV data was evaluated. For SK-I to SK-III, all channels in Table A.1 (Appendix) were used.
The signal and background estimations for each period and decay mode were added as independent
terms in Equation 9.1 to 9.4 and combined in Equation 9.5. Systematic uncertainties were assumed
to be fully correlated for the entire period from SK-I to SK-IV assuming common sources. As aresult, the limit of 3.6 × 1033 years at 90% C.L. was obtained from 0.37 Mton·years of data. This
limit is more than twice as long as the previous result, 1.6 × 1033 years, which uses data for SK-I
to SK-III. As explained in Chapter 10, the systematic uncertainties on the kaon scattering for SK-I
to SK-III periods were re-evaluated to cover the update of the kaon interaction model when the
combined lifetime limit was calculated. This aﬀects mainly the signal eﬃciency of KL0 selection.
The lifetime limit for SK-I to SK-III becomes 1.2 × 1033 years with this update. The lower limit
given by the combination of SK-I to SK-IV turns out to be smaller than that from only SK-IV
data. This can be explained due to two candidates in the ﬁnal KS0 → 2π 0 selection sample in SK-II
data compared to a background expectation of 0.20 events, which corresponds to a local p-value of
1.8%. On the other hand, no candidate events were found in the KS0 → 2π 0 selection in SK-IV.

83

Chapter 10

Systematic uncertainties
Tables 10.1 and 10.2 summarize the systematic uncertainties for the selection eﬃciencies and the
expected background rates, respectively. As the systematic uncertainties in the signal eﬃciencies
associated with the physics process, the correlated decay probability, Fermi momentum models,pion interaction and kaon interaction are considered. Since the ﬂight length of KL0 mostly depends
on the cross section as described in Section 4.1, the contribution of kaon interaction is larger thanthe other physics sources in KL0 selections. The contribution from pion interaction is negligible in
KS0 → 2π 0 selection as π 0 decay immediately.
Uncertainties in the neutrino ﬂux, neutrino cross section and pion interaction were accounted
for in the estimation of the atmospheric neutrino backgrounds. Since the major background events
are not kaon production as described in Section 8, the contribution from the kaon interaction is
negligible.
As the systematic uncertainties associated with the detector performance and reconstruction,
uncertainties in the number of events in the ﬁducial volume, detector non-uniformity, energy scale,ring counting, particle type identiﬁcation, decay electron tagging and neutron tagging are consideredin the signal eﬃciencies and background rates. In addition, uncertainties in the vertex separation
were evaluated for KL0 → π ± l∓ ν, KL0 → 3π 0 and KL0 → π + π − π 0 selections. Uncertainties inthe vertex separation are among the largest systematic uncertainties for the detector performanceand reconstruction in KL0 selections. Details of these systematic uncertainties are described in thischapter.

10.1

Correlated decay

As explained in Section 4.1, some theoretical models predict the proton decay process as a pair of
two nucleons in the oxygen nucleus (correlated decay). The probability of the correlated proton
decay within oxygen is predicted as 10% [35]. Considering the limitation of the reliability of themodel, the fraction of correlated decay events was changed to 0% and 20% by weighting events to84

Table 10.1: Systematic uncertainties on the signal eﬃciencies (%)
Sources
Physics
Correlated decay
Fermi momentum
Pion interaction
Kaon interaction
Reconstruction
Fiducial Volume
Energy scale
Non-uniformity
Ring counting
Particle identiﬁcation
Michel electron tagging
Vertex separation
Neutron tagging
Total

KS0 → 2π 0

KS0 → π + π −KL0 → π ± l∓ ν

KL0 → 3π 0

KL0 → π + π − π 0

6.1
1.0
3.1

6.3
1.8
4.5
3.1

5.6
1.3
2.7
11.6

3.9
0.1
1.8
5.3

2.7
1.0
6.3
10.3

4.0
1.9
0.2
0.4
0.4
5.4
0.6
9.9

3.4
1.6
0.2
0.4
0.2
5.4
10.8

0.5
0.9
2.8
0.6
5.4
10.5
2.6
18.2

0.5
4.7
2.8
5.1
5.4
4.8
2.2
12.6

0.5
2.3
2.8
0.8
5.4
9.2
3.0
17.0

Table 10.2: Systematic uncertainties on the number of the expected background events (%)
Sources
Physics
Neutrino ﬂux
Neutrino interaction
Pion interaction
Reconstruction
Fiducial Volume
Energy scale
Non-uniformity
Ring counting
Particle identiﬁcation
Michel electron tagging
Vertex separation
Neutron tagging
Total

KS0 → 2π 0

KS0 → π + π −KL0 → π ± l∓ ν

KL0 → 3π 0

KL0 → π + π − π 0

8.7
20.0
17.4

6.8
21.5
9.9

6.6
22.6
9.0

8.0
20.0
8.4

7.9
24.0
13.3

4.0
19.9
5.8
0.4
<0.1
5.4
2.7
35.5

3.4
2.7
1.1
0.4
<0.1
5.4
25.6

0.5
5.4
2.8
<0.1
5.4
13.9
2.7
28.9

0.5
14.6
2.8
28.9
5.4
10.7
2.9
41.4

0.5
0.1
2.8
<0.1
5.4
16.6
2.6
32.5

85

evaluate the systematic uncertainty.
ϵsys =nfree
NS + NP ∓ αNcorrnS + nP
ncorr+
×+ (1 ± α) ×
Ntot
NS + NP
Ntot
Ntot(10.1)

where n is the number of signal MC events passing the selection criteria, N is the number ofgenerated proton decay MC events in the ﬁducial volume,

freestands for free proton decays, S(P ) stands for proton decays in the S (P ) state in 16 O, andcorr

stands for correlated decays. A

factor α = 1, −1 corresponds to correlated decay probability of 20% and 0%, respectively. The mostdiscrepant ϵsys between the α = 1 and α = −1 from normal eﬃciency is taken for the uncertainty
estimation.

10.2

Fermi momentum

As explained in Section 4.1, protons have Fermi momentum up to a few hundred MeV in theoxygen nucleus. The Fermi momentum distribution is implemented in the proton decay basedon the experimental data [33]. On the other hand, the Fermi gas model [41] is employed in theatmospheric neutrino MC simulation (NEUT). This model uncertainty ϵsys is evaluated from thediﬀerence between the two models shown in Figure 10.1. The events in the ﬁnal sample are weighted
by the fraction of the Fermi momentum distributions with 11 bins (bin width is 25 MeV/c), andthe uncertainty was estimated by the deviation from the number of events before the weighting.
ϵsys =

11
nfree X ni
NNEUT,i+×,Ntot i=1 Ntot
NPD,i(10.2)

where i is the index of each momentum bin, NPD,i and NNEUT,i are the fraction of events in the
Fermi momentum distribution of the proton decay MC and atmospheric neutrino MC simulation,
respectively.

10.3

Pion interaction

The pion interaction model is implemented in the simulation with six internal parameters basedon the pion scattering experiments [74][75][76]: absorption probability fABS , quasi-elastic (QE)
scattering probability fQEL , charge exchange fraction of QE scattering fCXL , elastic scattering
probability fELH , charge exchange probability fCXH and inelastic scattering probability fHAD . The
applicable momentum region is below 500 MeV/c for the ﬁrst three parameters and above 500 MeV/c
for the other parameters. The nominal values of these parameters are(fQEL , fELH , fHAD , fABS , fCXL , fCXH ) = (1.069, 1.824, 1.002, 1.404, 0.697, 1.800).
To cover the possible parameter space from measurements, 24 patterns of the parameters areprepared (Table 10.3) in addition to the nominal one. For each case, events in the ﬁnal sample are86Arbitrary unit (area normalize)

0.2
0.18
0.16
0.14
0.12
0.1
0.08
0.06
0.04
0.02
0
0

50

100

150

200

250

300

Momentum (MeV/c)
Figure 10.1: Fermi momentum distributions of the proton decay MC (black) and atmospheric
neutrino MC simulation (red).
weighted by the fraction of the interaction probability to the nominal one, and the variations of
the signal eﬃciencies and background rates are calculated. The standard deviations of 24 biased
eﬃciencies and background rates over the nominal values are assigned as systematic uncertainties.

10.4

Kaon interaction

The simulation of the kaon interaction consists of two steps: interaction inside the oxygen nucleus inwhich proton decay or atmospheric neutrino interaction occurs and in water with the other nuclei.
For short distance (O(10−14 ) m) in the nucleus, the kaon interaction is simulated as an eigenstateof K 0 . Although there are no experimental data for the kaon-oxygen cross sections, the kaonoxygen cross section in NEUT was compared with the estimation by extrapolation from the kaondeuteron cross section measurement scaled by the atomic number to oxygen as shown in Figure 10.2.
This discrepancy is due to the interaction of kaons in the oxygen nucleus as described in Section
4.1. The deviation was estimated as ±25% to cover the discrepancy between the extrapolationfrom the measurements and the simulation around the 300 MeV/c momentum, and it was used
for the evaluation of the uncertainty. By the same method as the evaluation of the systematic
uncertainty for the Fermi momentum, events in the ﬁnal sample were weighted and the uncertainty
was estimated by the deviation from the number of events before the weighting. The weights were
determined by the fraction of the type of interaction (elastic or charge-exchange), kaon momentum
and direction between distributions with normal and biased cross sections. The maximum deviationof the eﬃciency from the nominal value was investigated with ±25% variation of the kaon interaction

87

Table 10.3: Parameter sets for the estimation of pion interaction uncertainty.
fQEL
0.6
0.6
0.7
0.7
1.4
1.3
1.5
1.6
0.6
0.6
0.7
0.7
1.4
1.3
1.5
1.6
0.6
0.6
0.7
0.7
1.4
1.3
1.5
1.6

fELH
1.8
1.8
1.8
1.8
1.8
1.8
1.8
1.8
1.1
1.1
1.1
1.1
1.1
1.1
1.1
1.1
2.3
2.3
2.3
2.3
2.3
2.3
2.3
2.3

fHAD
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.5
1.5
1.5
1.5
1.5
1.5
1.5
1.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5

fABS
0.7
0.7
1.6
1.6
0.6
0.7
1.5
1.6
0.7
0.7
1.6
1.6
0.6
0.7
1.5
1.6
0.7
0.7
1.6
1.6
0.6
0.7
1.5
1.6

88

fCXL
0.5
1.6
0.4
1.6
0.6
1.6
0.4
1.6
0.5
1.6
0.4
1.6
0.6
1.6
0.4
1.6
0.5
1.6
0.4
1.6
0.6
1.6
0.4
1.6

fCXH
1.8
1.8
1.8
1.8
1.8
1.8
1.8
1.8
2.3
2.3
2.3
2.3
2.3
2.3
2.3
2.3
1.3
1.3
1.3
1.3
1.3
1.3
1.3
1.3

Cross section (mb)

300
250
200
150
100
50
0

0.4

0.6

0.8
1
1.2
Momentum (GeV/c)

1.4

Cross section (mb)Figure 10.2: K 0 , 16 O cross section (elastic+charge exchange). Extrapolation from the experimental
data (black) and simulation in NEUT (Red) are shown. The shaded region is a ±25% area of thecross section in NEUT.

600
500
400
300
200
100
0

0.4

0.6

0.8
1
1.2
Momentum (GeV/c)

1.4

Figure 10.3: Total KL0 , 16 O cross section. Extrapolation from the experimental data (black) andsimulation in water (Red) are shown. The shaded region is a ±15% area of the cross section inwater.
cross section and taken as the uncertainty.
For longer distances in water, the kaon interaction is simulated as eigenstates of KS0 and KL0 . The
uncertainty of the KL0 cross section in water was evaluated from the comparison of the simulation
model and the independent measurement of KL0 -nucleon scattering [77] as well as K 0 cross section(Figure 10.3). The deviation was estimated to be ±15% from this comparison. The systematic
uncertainties were estimated by the weights determined by the fraction of the type of interaction(hadronic inelastic scattering, KS0 decay or KL0 decay), kaon momentum before interaction, andvertex separation between distributions from the variation with nominal and biased cross sections.
A quadratic sum of the uncertainties in nuclei and water is assigned as a systematic uncertaintyof the kaon interaction in Table 10.1. As the kaon interaction model is updated from the previous
analysis [1], systematic uncertainties for eﬃciencies of the SK-I to SK-III search in the previous

89analysis were re-evaluated to cover the diﬀerence between the old and new models. They aresummarized in Table A.1 (Appendix).

10.5

Neutrino ﬂux and cross section

The uncertainties of the neutrino ﬂux and cross section are based on the atmospheric neutrino
oscillation analysis in SK [64]. Events in the ﬁnal sample are weighted with the ratio of two models
or the estimated range of the variations, and the uncertainties are evaluated as the diﬀerence
between the number of events with and without the weighting. Tables 10.4 and 10.5 summarize
the systematic uncertainties in the background rates associated with the neutrino ﬂux and crosssection models, respectively.
Table 10.4: Systematic uncertainties on the background rates associated with the neutrino ﬂux
model. Eν is neutrino energy.
SourceKS0 → 2π 0 KS0 → π + π − KL0 → π ± l∓ ν KL0 → 3π 0 KL0 → π + π − π 0
Absolute normalization(Eν < 1 GeV)
<0.1%
3.2%
3.6%
1.2%
<0.1%
(1 < Eν GeV)
7.0%
4.3%
4.0%
6.4%
7.0%
νµ /νe ratio(Eν < 1 GeV)
<0.1%
0.4%
0.4%
<0.1%
<0.1%
(1 < Eν < 10 GeV)
0.6%
0.8%
0.7%
0.5%
1.2%
(10 < Eν GeV)
<0.1%
<0.1%
<0.1%
<0.1%
<0.1%
ν̄e /νe ratio(Eν < 1 GeV)
<0.1%
<0.1%
<0.1%
<0.1%
<0.1%
(1 < Eν < 10 GeV)
0.7%
<0.1%
<0.1%
0.7%
0.3%
(10 < Eν GeV)
<0.1%
<0.1%
<0.1%
<0.1%
<0.1%
ν̄µ /νµ ratio(Eν < 1 GeV)
<0.1%
0.4%
0.4%
<0.1%
<0.1%
(1 < Eν < 10 GeV)
1.8%
1.5%
1.5%
1.9%
2.7%
(10 < Eν GeV)
<0.1%
<0.1%
<0.1%
<0.1%
<0.1%
Asymmetry
Up/down
<0.1%
<0.1%
<0.1%
<0.1%
<0.1%
Horizontal/vertical
0.1%
<0.1%
0.2%
0.3%
<0.1%
K/π ratio
4.7%
3.7%
3.5%
4.2%
2.3%
Neutrino ﬂight length
<0.1%
<0.1%
<0.1%
<0.1%
<0.1%
Total
8.7%
6.8%
6.6%
8.0%
7.9%

10.6

Fiducial volume

A systematic uncertainty for the FV is estimated for the possible diﬀerence in the number of
events in the FV between data and the MC simulation. Sub-GeV (Evis < 1330 MeV) multi-ringsamples of the data and the atmospheric neutrino MC simulation were used for this estimation.
To determine a stable region where the vertex reconstruction accuracy does not aﬀect the numberof events, reconstructed Dwall (distance between the interaction point and the nearest ID wall)90

Table 10.5: Systematic uncertainties on the background rates associated with the neutrino crosssection model.
SourceKS0 → 2π 0 KS0 → π + π − KL0 → π ± l∓ ν KL0 → 3π 0 KL0 → π + π − π 0
Axial vector mass
10.7%
11.1%
13.7%
5.2%
11.2%
NCEL and CCQE ratio
0.2%
<0.1%
1.1%
1.2%
<0.1%
ν̄/ν ratio of CCQE
<0.1%
<0.1%
0.3%
0.2%
<0.1%
µ/e ratio of CCQE
<0.1%
<0.1%
<0.1%
<0.1%
<0.1%
Single meson production
12.1%
16.4%
17.2%
5.8%
11.5%
ν̄/ν ratio of single π
2.8%
4.2%
4.6%
8.0%
4.5%
DIS (model dependence)
11.2%
6.4%
0.8%
16.0%
17.2%
DIS (total cross section)
1.9%
0.9%
0.2%
3.1%
2.1%
Coherent pion production
<0.1%
<0.1%
<0.1%
<0.1%
<0.1%
NC/CC ratio
1.2%
3.0%
0.8%
2.8%
<0.1%
20.0%
22.6%
24.0%
Total
21.5%
20.0%

distribution was compared with the true values in the MC simulation, MC true values were smeared
by Gaussian with σ = 40 cm (green) and σ = 100 cm (blue) to account for the resolution of the vertex
reconstruction as shown in Figure 10.4. As a result, the ratio was stable in 600 < Dwall < 1000 cm
region. Therefore, the MC simulation was normalized to data in this region, and the number of
events in the FV (Dwall > 200 cm) was compared between the data and MC simulation (Figure 10.5).
Since the deﬁnition of the reconstructed primary vertex is diﬀerent for selections, uncertainties were
evaluated for three cases: KS0 → 2π 0 , KS0 → π + π − and KL0 decay selections.

Figure 10.4: Left: distance from the wall of the ID. Reconstructed distribution (black) is comparedwith MC true (red), MC true smeared by Gaussian with σ = 40 cm (green) and σ = 100 cm (blue).
Right: ratio of MC true distributions to the reconstructed one.

10.7

Energy scale

As described in Section 3.8, the energy calibration is performed using independent control samples
with known energies. The energy scale uncertainty was evaluated as 2.17% from the diﬀerence91Figure 10.5: Reconstructed distance between the wall of the ID and the position of the particle
generation in the sub-GeV multi-ring sample. The simulation (red) are normalized to the data
(black) in 600 < Dwall < 1000 cm region. From left to right, the distribution with the deﬁnition of
the reconstructed primary vertex for KS0 → 2π 0 , KS0 → π + π − and KL0 decay selections.
between the data and MC simulation[69]. It aﬀects the muon momentum, total momentum andinvariant mass of kaon and proton. The uncertainties were evaluated from the variation of thenumber of events in the ﬁnal samples with a bias in the energy scale to these variables.

10.8

Detector non-uniformity

Non-uniformity of the detector response, such as due to PMT and water transparency, is one of
the sources of uncertainty in the measurement. Such non-uniformity causes bias in the energy scaledepending on the direction of the particles (corresponds to the location of PMTs where Cherenkov
light is detected). The uncertainty was evaluated to be 0.58% in the energy scale [69]. Thisaﬀects the total momentum since the directions of µ+ and K 0 from proton decay are back-to-back.
The impact on the signal eﬃciency and the expected background rates were evaluated by the MC
simulation as shown in Table 10.1 and 10.2.

10.9

Ring counting

The number of rings is used in all selections in the p → µ+ + K 0 search. In the reconstruction,
likelihood values are compared for two hypotheses with n rings and n+1 rings and the additional ringis adopted if the likelihood with n + 1 rings exceeds the threshold (see Chapter 7). In this analysis,
systematic uncertainty in the number of rings was evaluated from the likelihood distributions.
Figure 10.6 shows the ring counting likelihood of data and simulation in the sub-GeV FCFV
sample. The scale and shift parameters for the simulation were determined to reproduce data. The
best scale and shift parameters in multi-ring (displaced) ﬁtter are 0.998 ± 0.006 and 0.17 ± 0.02
(1.026 ± 0.005 and 0.38 ± 0.02), respectively. The uncertainties are evaluated as the variation of
the multi-ring events fraction of the sub-GeV FCFV sample in this bias.

92

Figure 10.6: Ring counting likelihood of data (black), simulation (red) and simulation with thebest-ﬁt parameters for the correction (green) in sub-GeV FCFV sample. The likelihood by multiring ﬁtter (left) and displaced vertex ﬁtter (right) are shown.

10.10

Particle type identiﬁcation

Systematic uncertainty in the particle type identiﬁcation was also estimated by a similar approachas that for the number of rings (Section 10.9) from the likelihood distributions calculated by theevent reconstruction algorithm. Figure 10.7 shows the particle type likelihood for the ﬁrst ringof the data and simulation in sub-GeV (Evis < 1330 MeV) FCFV multi-ring sample. The scaleand shift parameters for the simulation were determined to reproduce the data for each ring aswell as the ﬁrst ring. Table 10.6 shows the best-ﬁt parameters for each ring. Fifth and sixth ring
distributions are merged to accumulate statistics. These biases were applied to the simulation anduncertainties are evaluated from the variation of the number of events in the ﬁnal sample of the
proton decay event selections.
Table 10.6: The best-ﬁt parameters of the particle type likelihood for each ring.
Multi-ring ﬁtter
Displaced vertex ﬁtter
Ring
Scale
Shift
Scale
Shift
First
0.96 ± 0.01
2.7 ± 2.4
1.01 ± 0.01
1.8 ± 3.0
Second
0.97 ± 0.01
0.6 ± 1.2
1.00 ± 0.01 −1.9 ± 1.3
Third0.97 ± 0.01 −2.2 ± 1.3 0.99 ± 0.02 −2.0 ± 1.6
Fourth
0.97 ± 0.02 −0.6 ± 1.8 1.06 ± 0.04 −3.9 ± 3.5
Fifth and Sixth 0.95 ± 0.03 −2.5 ± 2.3 1.14 ± 0.10 −23.2 ± 7.0

93Figure 10.7: Particle type likelihood for the ﬁrst ring of data (black), simulation (red) and simulation
with the best-ﬁt parameters for the correction (green) in sub-GeV FCFV multi-ring sample. The
likelihood by multi-ring ﬁtter (left) and displaced vertex ﬁtter (right) are shown.

10.11

Michel electron tagging

The number of Michel electrons after the primary event is related to the number of muons and also
charged pions as the parent particle of muons. Therefore, the information is used in the selectionsof the p → µ+ + K 0 search. Figure 10.8 shows the number of Michel electrons of the data andthe atmospheric neutrino simulation in the sub-GeV FCFV multi-ring sample. Both distributions
are normalized by their areas. Systematic uncertainty is evaluated to account for the relative
error between the data and MC simulation in the number of events with zero, one and two Michelelectrons.

Figure 10.8: The number of Michel electrons of data (black) and atmospheric neutrino simulation(red) in sub-GeV FCFV multi-ring sample.

9410.12

Vertex separation

A systematic uncertainty of vertex separation reconstruction is evaluated by comparison between
true and reconstructed vertex separation in the proton decay simulation. All cuts before the vertex
separation cut were applied. Figure 10.9 shows the diﬀerence between true and reconstructed vertex
separation. Table 10.7 shows the results of the Gaussian ﬁtting to these distributions. The meanvalues from the ﬁts are applied to the simulation as bias, and the variations of the number of events

Arbitrary unit (area normalize)in the ﬁnal sample were taken as the systematic uncertainties.

Reconstructed − True vertex separation (cm)Figure 10.9: Distribution of reconstructed − true vertex separation in the proton decay simulation
to which all cuts before the vertex separation cut were applied. The red line shows the Gaussianﬁt result.

Table 10.7: The results of the Gaussian ﬁttings to the reconstructed − true vertex separation
distributions. µ is a mean and σ is a standard deviation.
Selectionµ (cm)σ (cm)
KL0 → π ± l∓ ν 10.3 ± 0.4 40.0 ± 0.4
5.9 ± 1.0 42.4 ± 1.0
KL0 → 3π 0
6.8 ± 0.9 39.1 ± 0.8
KL0 → π ± l∓ ν
It was also checked if the deviation between the data and the simulation is consistent with these
uncertainties. The agreement of the data and simulation was evaluated using muons and their
Michel electrons. The displaced vertex ﬁtter was applied to atmospheric neutrino events with thefollowing criteria to select νµ CCQE events around the lifetime of KL0 (50 ns).
• Events should pass FCFV selection and are sub-GeV.
• The number of rings should be two with one non-showering (µ) and one showering (e) rings.
• The number of Michel electrons is one or zero.
• The momentum of the non-showering ring is 450 < pµ < 550 MeV/c.
• The momentum of the showering ring is 10 < pe < 80 MeV/c.
95• Time diﬀerence between two rings is less than 100 ns.
Figure 10.10 shows the vertex separations of the events after the selection. It was found fromthe comparison that the distribution of the simulation best reproduces the data with a shift of(7.0 ± 3.9) cm. This is a comparable size to the systematic uncertainties in Table 10.7.

Figure 10.10: The vertex separations of the data (black), the simulation (red) and the simulation
with the correction (green) after νµ CCQE event selection. The simulation is normalized by thelive time.

10.13

Neutron tagging eﬃciency

Uncertainty in the neutron tagging eﬃciency is evaluated as ±9% [72]. Systematic uncertainty
on the proton decay eﬃciencies and the expected backgrounds are evaluated from the variation of
the number of events in the ﬁnal sample with diﬀerent neutron tagging eﬃciency by ±9%. The
neutron tagging eﬃciencies in the ﬁnal sample are determined by applying binominal probabilities tothe true neutron distributions to reproduce the tagged neutron distributions in the sample before
the neutron cut. For example, Figure 10.11 shows the number of neutron distributions in theKL0 → π ± l∓ ν sample. Several binominal probabilities are applied to the true number of neutrons,
and the variation of the event fraction between before and after the neutron cut in ±9% around thebest tagging eﬃciency (about 28.7%) is evaluated as the systematic uncertainty (Figure 10.12).

96Figure 10.11: The number of tagged neutrons (black), true neutrons (red) and true neutrons with28.7% tagging eﬃciency (green) in the KL0 → π ± l∓ ν sample before the neutron cut.

Figure 10.12: The fraction of events with no tagged neutron in the KL0 → π ± l∓ ν sample before
the neutron cut as a function of the several tagging eﬃciencies.

97Chapter 11

Conclusion
Proton decay into a muon and a neutral kaon (p → µ+ + K 0 ) was searched for using the data fromthe Super-Kamiokande water Cherenkov detector. The sensitivity was improved from the previous
search [1] with a new event reconstruction algorithm, neutron tagging and optimized selection
criteria for each KS0 and KL0 decay mode. In addition, the simulation of the kaon interaction was
updated based on the external measurement of the kaon scattering.
K 0 is a superposition of KS0 with a lifetime of 90 ps and KL0 with a lifetime of 51 ns. Therefore,p → µ+ + K 0 was searched by diﬀerent approaches for KS0 and KL0 decay modes. KS0 decays atalmost the same position as proton decay and most of the particles from KS0 decay can be measured
in the Super-Kamiokande detector. Therefore, signal events from p → µ+ + K 0 are selected byrequiring an invariant mass of KS0 and proton. On the other hand, KL0 decays at a distance ofO(1) m away from the position of the proton decay after scattering in water and the secondary
particles from KL0 decay create Cherenkov rings at the secondary vertex. Since the primary muonof p → µ+ + K 0 creates a clear Cherenkov ring at the primary vertex of the proton decay, twovertices from the muon and the secondary particles from KL0 decay exist in p → µ+ + KL0 channel.
For this reason, the multi-ring reconstruction algorithm was modiﬁed for p → µ+ + KL0 event
selection to reconstruct a muon from proton decay and particles generated from KL0 decay at thedisplaced vertex. To eﬀectively distinguish the p → µ+ + KS0 and p → µ+ + KL0 events from thebackgrounds, ﬁve selection criteria are developed for each KS0 and KL0 decay channel.
As a result, no signiﬁcant signal has been observed in the ﬁnal sample of 0.2 Mton·years of datain the SK-IV period (2008–2018). From this result, a lower limit of 4.5 × 1033 years on the lifetimeof p → µ+ + K 0 was obtained at 90% C.L. By combining this with the previous results using thedata in SK-I to SK-III period (1996–2008), a lower limit of 3.6 × 1033 years was obtained from
0.37 Mton·years data. This limit is more than twice as long as the previous result of 1.6 × 1033 years[1] and is the most stringent for this channel. The result of this study provides knowledge inelementary particle physics by the constraints on the grand uniﬁed theories.

98

Acknowledgments
I am deeply grateful to my supervisor Masaki Ishitsuka. He spent a great deal of time and eﬀort
giving me many comments and advice which were always thoughtful and suggestive and helped me
grow. This thesis would never completed without his support.
I would like to thank all members of the SK group. In particular, Makoto Miura and Shunichi
Mine made many speciﬁc and important comments on the technical details in the analysis of proton
decay. I received many useful suggestions not only for analysis but also for writing the paper. I
also thank Edward Kearns and Jost Migenda. They gave me many very essential suggestions forthe analysis and writing the paper. I would like to thank Yoshinari Hayato and Roger Wendellfor giving me many suggestions, especially for the simulation of proton decay and atmospheric
neutrino. I express my gratitude to Akira Takenaka for advice on the analysis of proton decay andFC reduction.
This work was supported by Grant-in-Aid for JSPS Research Fellow Grant Number JP22J14531.
I thank the ﬁnancial support.
Finally, I am very grateful to my family for supporting my life as a student for a long period.

99

Appendix A

p → µ+ + K 0 search in SK-I to
SK-III
Table A.1 shows the summary of p → µ+ + K 0 search in SK-I to SK-III. These systematic uncertainties are re-evaluated due to the updated kaon interaction model. The selection criteria of SK-I,
SK-II and SK-III are diﬀerent from those of SK-IV.
Table A.1: Summary of the p → µ+ K 0 search in SK-I to SK-III [1].
Search mode
SK-I (91.7 kton·years)
KS0 → 2π 0
KS0 → π + π − (2 ring sample)
KS0 → π + π − (3 ring sample)
KL0
SK-II (49.2 kton·years)
KS0 → 2π 0
KS0 → π + π − (2 ring sample)
KS0 → π + π − (3 ring sample)
KL0
SK-III (31.9 kton·years)
KS0 → 2π 0
KS0 → π + π − (2 ring sample)
KS0 → π + π − (3 ring sample)
KL0

Eﬃciency (%)

Background (events)

Candidates

7.0 ± 0.7
10.6 ± 1.1
2.5 ± 0.3
3.8 ± 1.8

0.37 ± 0.05
3.0 ± 0.5
0.12 ± 0.08
3.5 ± 1.10602

6.2 ± 0.8
10.3 ± 1.3
2.4 ± 0.3
3.3 ± 1.6

0.20 ± 0.05
1.6 ± 0.4
0.23 ± 0.08
1.4 ± 0.5201
06.7 ± 0.9
10.3 ± 1.9
3.0 ± 0.3
3.8 ± 1.8

0.19 ± 0.04
1.2 ± 0.2
0.09 ± 0.02
1.3 ± 0.6

0
1
0
1

i

Bibliography[1] C. Regis et al., ”Search for Proton Decay via p → µ + K 0 in Super-Kamiokande I, II, and III”,
Phys. Rev. D 86, 012006 (2012).
[2] R. L. Workman et al. (Particle Data Group), Prog. Theor. Exp. Phys. 2022, 083C01 (2022).
[3] H. Georgi and S. L. Glashow, ”Unity of All Elementary Particle Forces”, Phys. Rev. Lett. 32,
438–441 (1974).
[4] A. D. Sakharov, ”Violation of CP-invariance, C-asymmetry, and baryon asymmetry of the Universe”, In The Intermissions… Collected Works on Research into the Essentials of Theoretical
Physics in Russian Federal Nuclear Center, Arzamas-16, 84–87 (1998).
[5] P. Langacker, Grand uniﬁed theories and proton decay, Phys. Rep. 72, 185–385 (1981).
[6] C. McGrew et al. (IMB Collaboration), Search for nucleon decay using the IMB-3 detector,
Phys. Rev. D 59, 052004 (1999).
[7] K. S. Hirata et al. (Kamiokande Collaboration), Experimental limits on nucleon lifetime for
lepton + meson decay modes, Phys. Lett. B 220, 308–316 (1989).
[8] A. Takenaka et al. (The Super-Kamiokande Collaboration), Search for proton decay via p →e + π 0 and p → µ + π 0 with an enlarged ﬁducial volume in Super-Kamiokande I-IV, Phys. Rev.
D 102, 112011 (2020).
[9] K. S. Babu, Ilia Gogoladze and Cem Salih Ün, ”Proton lifetime in minimal SUSY SU (5) inlight of LHC results”, J. High Energ. Phys. 2022, 164 (2022).
[10] S. M. Barr, ”A new symmetry breaking pattern for SO(10) and proton decay”, Phys. Lett. B
112, 219–222 (1982).
[11] K. S. Babu and S. Khan, ”Minimal nonsupersymmetric SO(10) model: Gauge coupling uniﬁcation, proton decay, and fermion masses”, Phys. Rev. D 92, 075018 (2015).
[12] K. S. Babu, J.C.Pati and F.Wilczek, Suggested new modes in supersymmetric proton decay,
Phys. Lett. B 423, 337–347 (1998); K. S. Babu, J. C. Pati, and F. Wilczek, Fermion masses,iineutrino oscillations, and proton decay in the light of SuperKamiokande, Nucl. Phys. B 566,
33–91 (2000); J. C. Pati, Confronting the conventional ideas of grand uniﬁcation with fermion
masses, neutrino oscillations and proton decay, arXiv preprint hep-ph/0204240 (2002).
[13] K. Hamaguchi, S. Hor and N. Nagata, ”R-Symmetric Flipped SU(5)”, J. High Energ. Phys.
2020, 140 (2020).
[14] K. Abe et al. (The Super-Kamiokande Collaboration), Search for proton decay via p → νK +using 260 kiloton·year data of Super-Kamiokande, Phys. Rev. D 90, 072005 (2014).
[15] M. Tanaka et al. (The Super-Kamiokande Collaboration), Search for proton decay into threecharged leptons in 0.37 megaton-years exposure of the Super-Kamiokande, Phys. Rev. D 101,
052011 (2020).
[16] K. Abe et al. (The Super-Kamiokande Collaboration), Search for Nucleon Decay via n → ν̄π 0
and p → ν̄π + in Super-Kamiokande, Phys. Rev. Lett. 113, 121802 (2014).
[17] K. Abe et al. (The Super-Kamiokande Collaboration), Search for nucleon decay into charged
antilepton plus meson in 0.316 megaton·years exposure of the Super-Kamiokande water
Cherenkov detector, Phys. Rev. D 96, 012003 (2017).
[18] V. Takhistov et al. (The Super-Kamiokande Collaboration), Search for Nucleon and Dinucleon
Decays with an Invisible Particle and a Charged Lepton in the Final State at the SuperKamiokande Experiment, Phys. Rev. Lett. 115, 121803 (2015).
[19] V. Takhistov et al. (The Super-Kamiokande Collaboration), Search for Trilepton Nucleon Decay via p → e+ νν and p → µ+ νν in the Super-Kamiokande Experiment, Phys. Rev. Lett. 113,
101801 (2014).
[20] K. Kobayashi et al., ”Search for nucleon decay via modes favored by supersymmetric granduniﬁcation models in Super-Kamiokande-I”, Phys. Rev. D 72, 052007 (2005).
[21] R. Matsumoto et al., ”Search for proton decay via p → µ+ K 0 in 0.37 megaton-years exposureof Super-Kamiokande”, Phys. Rev. D 106, 072003 (2022).
[22] R. Matsumoto and M. Ishitsuka for the Super-Kamiokande Collaboration, ”Classiﬁcation ofSuper-Kamiokande atmospheric neutrino events by using neural network”, J. Phys.: Conf. Ser.
1468, 012161 (2020).
[23] S. Fukuda, et al. (Super-Kamiokande Collaboration), ”The Super-Kamiokande Detector”,
Nucl. Instr. Meth. A 501 418–462 (2003).

iii

[24] S. Yamada et al., ”Commissioning of the New Electronics and Online System for the SuperKamiokande Experiment”, IEEE Trans. Nucl. Sci. 57, 428–432 (2010).
[25] Kamioka Observatory, ICRR, The University of Tokyo, ”Super-Kamiokande”, http://wwwsk.icrr.u-tokyo.ac.jp/sk/[26] Kamioka Observatory, ICRR, The University of Tokyo, ”Hyper-Kamiokande”, http://wwwsk.icrr.u-tokyo.ac.jp/hk[27] H. Kume et al., ”20 INCH DIAMETER PHOTOMULTIPLIER”, Nucl. Instrum. Meth. 205(1983) 443–449.
[28] A. Suzuki et al., ”Improvement of 20-inch diameter photomultiplier tubes”, Nucl. Instrum.
Meth. A329 (1993) 299–313.
[29] K. Abe et al., ”Calibration of the Super-Kamiokande detector”, Nucl. Instrum. Meth. A 737
253–272 (2014).
[30] Y. Nakano et al., Measurement of the radon concentration in puriﬁed water in the SuperKamiokande IV detector, Nucl. Instrum. Meth. A 977, 164297 (2020).
[31] T. Yokozawa, ”Precision solar neutrino measurements with Super-Kamiokande-IV”, Ph.D.
thesis, University of Tokyo, (2012).
[32] M. Jiang, ”Study of the neutrino mass hierarchy with the atmospheric neutrino data collected
in Super-Kamiokande IV”, Ph.D. thesis, Kyoto University, (2019).
[33] K. Nakamura et al., ”The reaction 12 C(e,e’p) at 700 MeV and DWIA analysis”, Nucl. Phys.
A 268, 381–407 (1976).
[34] M. G. Mayer and J. H. D. Jensen, ”Elementary Theory of Nuclear shell structure”, Wiley, NewYork (1955).
[35] T. Yamazaki and Y. Akaishi, ”Nuclear medium eﬀects on invariant mass spectra of hadrons
decaying in nuclei”, Phys. Lett. B 453, 1–6 (1999).
[36] C. W. de Jager, H. de Vries and C. de Vries, ”Nuclear charge- and magnetization-densitydistribution parameters from elastic electron scattering”, Atomic Data and Nuclear Data Tables 14, 479–508 (1974).
[37] L. L. Salcedo, E. Oset, M. J. Vicente-Vacas and C. Garcia-Recio, ”Computer simulation ofinclusive pion nuclear reactions” Nucl. Phys. A 484, 557–592 (1988).

iv

[38] J. S. Hyslop, et al., ”Partial-wave analysis of K + -nucleon scattering”, Phys. Rev. D 46, 961–969(1992).
[39] R. Weiss et al., ”Measurement of low energy K + total cross sections on N = Z nuclei”, Phys.
Rev. C 49, 2569–2577 (1994).
[40] R. Brun, F. Bruyant, F. Carminati, S. Giani, M. Maire, A. McPherson, G. Patrick, and
L. Urban, ”GEANT Detector Description and Simulation Tool”, CERN Report No. CERNW5013 (1994).
[41] Y. Hayato, ”NEUT”, Nucl. Phys. Proc. Suppl. 112 (2002) 171-176; Y. Hayato, L. Pickering,”The NEUT neutrino interaction simulation program library”, Eur. Phys. J. Spec. Top. 230,
4469–4481 (2021).
[42] B. R. Martin and M. K. Pidcock, ”K̄N interactions in the resonance region (I) Analysis ofdata”, Nuc. Phys. B 126, 266–284 (1977); B. R. Martin and M. K. Pidcock, ”K̄N interactionsin the resonance region (II) Amplitudes”, Nuc. Phys. B 126 285–297 (1977).
[43] P. H. Eberhard and F. Uchiyama, ”Neutral-kaon regeneration probabilities at asymmetricϕ-factory energies”, Nucl. Instrum. Meth. A 350, 144–149 (1994).
[44] T. Kajita, ”Atmospheric Neutrinos And Discovery Of Neutrino Oscillations”, Proc. Japan
Acad. B 86, 303–321 (2010).
[45] M. Honda, T. Kajita, K. Kasahara and S. Midorikawa, ”Calculation of ﬂux of atmospheric
neutrinos”, Phys. Rev. D 52, 4985–5005 (1995).
[46] M. Honda, T. Kajita, K. Kasahara and S. Midorikawa, ”New calculation of the atmospheric
neutrino ﬂux in a 3-dimensional scheme”, Phys. Rev. D 70, 043008 (2004).
[47] M. Honda, T. Kajita, K. Kasahara, S. Midorikawa, and T. Sanuki, ”Calculation of atmospheric
neutrino ﬂux using the interaction model calibrated with atmospheric muon data”, Phys. Rev.
D 75, 043006 (2007).
[48] M. Honda, T. Kajita, K. Kasahara, and S. Midorikawa, ”Improvement of low energy atmospheric neutrino ﬂux calculation using the JAM nuclear interaction model”, Phys. Rev. D 83
123001 (2011).
[49] G. Battistoni, A. Ferrari, T. Montaruli, and P. R. Sala., ”The FLUKA atmospheric neutrino
ﬂux calculation”, Astropart. Phys. 19, 269–290 (2003).
[50] G. D. Barr, T. K. Gaisser, P. Lipari, Simon Robbins, and T. Stanev. ”A Three-dimensional
calculation of atmospheric neutrinos”, Phys. Rev. D 70, 023006 (2004).
v

[51] Y. Ashie, J. Hosaka, K. Ishihara, Y. Itow, J. Kameda, Y. Koshio et al., ”Measurement of
atmospheric neutrino oscillation parameters by Super-Kamiokande I”, Phys. Rev. D 71, 112005(2005).
[52] C.H.Llewellyn Smith, ”Neutrino reactions at accelerator energies”, Phys. Rep. 3C, 261–379(1972).
[53] R. Gran, J. Nieves, F. Sanchez and M. J. Vincente Vacas, ”Neutrino-nucleus quasi-elastic and2p2h interactions up to 10 GeV”, Phys. Rev. D 88, 113007 (2013).
[54] B. Bourguille, and J. Nieves, and F. Sanchez, ”Inclusive and exclusive neutrino-nucleus crosssections and the reconstruction of the interaction kinematics”, J. High Energ. Phys. 2021, 4(2021).
[55] R. Gran et al. (K2K Collaboration), ”Measurement of the quasielastic axial vector mass inneutrino interactions on oxygen”, Phys. Rev. D 74, 052002 (2006).
[56] A. A. Aguilar-Arevalo et al. (MiniBooNE Collaboration), ”First measurement of the muon
antineutrino double-diﬀerential charged-current quasielastic cross section”, Phys. Rev. D 88,
032001 (2013).
[57] P. Adamson et al. (MINOS Collaboration), ”Study of quasielastic scattering using chargedcurrent νµ -iron interactions in the MINOS near detector”, Phys. Rev. D 91, 012005 (2015).
[58] C. Berger and L. M. Sehgal, ”Lepton mass eﬀects in single pion production by neutrinos”,
Phys. Rev. D 76, 113004 (2007).
[59] C. Berger and L. M. Sehgal, ”Partially conserved axial vector current and coherent pion production by low energy neutrinos”, Phys. Rev. D 79, 053003 (2009).
[60] D. Rein, ”Diﬀractive pion production in neutrino reactions”, Nucl. Phys. B 278, 61–77 (1986).
[61] M. Derrick et al., ”Properties of the hadronic system resulting from ν̄µ p interactions”, Phys.
Rev. D 17, 1–15 (1978).
[62] B. JongeJans et al., ”Multiplicity distributions of charged hadrons produced in (anti)neutrinodeuterium charged- and neutral-current interactions”, Nuovo Cim. A 101, 435–453 (1989).
[63] T. Sjostrand, ”High-energy-physics event generation with PYTHIA 5.7 and JETSET 7.4”,
Comput. Phys. Commun. 82, 74–89 (1994).
[64] K.Abe et al., ”Atmospheric neutrino oscillation analysis with external constraints in SuperKamiokande I-IV”, Phys. Rev. D 97 072001 (2018).

vi

[65] S. Ioﬀe and C. Szegedy, ”Batch Normalization: Accelerating Deep Network Training by Reducing Internal Covariate Shift”, Proc. ICML 37 448-456 (2015).
[66] C. Ishihara, ”Full three ﬂavor oscillation analysis of atmospheric neutrino data observed inSuper-Kamiokande”, Ph.D. thesis, University of Tokyo (2010).
[67] A. Takenaka, ”Search for Proton Decay via p → e + π 0 and p → µ + π 0 with an Enlarged
Fiducial Mass of the Super-Kamiokande Detector”, Ph.D. thesis, University of Tokyo (2020).
[68] R. B. Patterson, E. M. Laird, Y. Liu, P. D. Meyers, I. Stancu, and H. A. Tanaka, The extendedtrack reconstruction for MiniBooNE, Nucl. Instrum. Meth. A 608, 206–224 (2009).
[69] M. Jiang et al. (Super-Kamiokande Collaboration), Atmospheric neutrino oscillation analysis
with improved event reconstruction in Super-Kamiokande IV, Prog. Theor. Exp. Phys. 2019,
053F01 (2019).
[70] F. James,

MINUIT - Function Minimization and Error Analysis - Reference Manual Version

94.1 , CERN-D-506 (1994).
[71] S. Tobayama, ”An Analysis of the Oscillation of Atmospheric Neutrinos”, Ph.D. thesis, The
University of British Columbia, (2016).
[72] K. Abe et al., ”Neutron Tagging following Atmospheric Neutrino Events in a Water Cherenkov
Detector”, JINST 17 P10029 (2022).
[73] H. Ejiri, ”Nuclear deexcitations of nucleon holes associated with nucleon decays in nuclei”,
Phys. Rev. C 48, 1442–1444 (1993).
[74] L. L. Salcedo, E. OsetM, J. Vicente-Vacas and C. Garcia-Recio, ”Computer simulation ofinclusive pion nuclear reactions”, Nucl. Phys. A 484, 557–592 (1988).
[75] P. de Perio, ”NEUT pion FSI”, AIP Conf. Proc. 1405, 223–228 (2011).
[76] P. de Perio, ”Joint Three-Flavour Oscillation Analysis of νµ Disappearance and νe Appearance
in the T2K Neutrino Beam”, Ph. D. thesis, University of Toronto (2014).
[77] W. E. Cleland et al., Measurement of the KL0 -p and KL0 -d total cross sections, Phys. Rev. D
12, 1247–1259 (1975).
[78] B. P. Roe and M. B. Woodroofe, Setting conﬁdence belts, Phys. Rev. D 63, 013009 (2000).

vii