Search for p → ν̄K + in the Super-Kamiokande Experiment
A Dissertation presented
byGabriel Chicca Santucci
to
The Graduate Schoolin Partial Fulfillment of theRequirements
for the Degree ofDoctor of Philosophy
inPhysics

Stony Brook University

December 2018

Stony Brook University
The Graduate School

Gabriel Chicca Santucci

We, the dissertation committee for the above candidate for theDoctor of Philosophy degree, hereby recommend
acceptance of this dissertation

Michael Wilking - Dissertation Advisor
Department of Physics and Astronomy, Stony Brook University

Chang Kee Jung - Dissertation Advisor
Department of Physics and Astronomy, Stony Brook University

John Hobbs - Chairperson of Defense
Department of Physics and Astronomy, Stony Brook University

Dmitri Kharzeev
Department of Physics and Astronomy, Stony Brook University

Rafael Coelho Lopes de Sa
Physics Department, UMass Amherst

Verena Martinez Outschoorn
Physics Department, UMass Amherst

William Marciano
Physics Department, Brookhaven National Laboratory
This dissertation is accepted by the Graduate School
Charles Taber
Dean of the Graduate Schoolii

Abstract of the Dissertation
Search for p → ν̄K + in the Super-Kamiokande Experiment
byGabriel Chicca Santucci
Doctor of Philosophy
inPhysics
Stony Brook University
2018

Super-Kamiokande is a large water Cherenkov detector located deep undergroundin Kamioka, Japan. Recently, a maximum likelihood algorithm, fiTQun, was developed to reconstruct the events in Super-Kamiokande. The present work consists ofa new search for the proton decay mode p → ν̄K + in the prompt-γ channel usingfiTQun as the event reconstruction algorithm.
A new event hypothesis was developed, dedicated to the prompt-γ channel. An
increase of 56% in signal efficiency was achieved using the new method comparedto the previous result. The expected number of background events was reevaluatedand the systematic uncertainties for this search were also updated to include an
uncertainty for kaon production in neutrino interactions.
The analysis was performed using 2867.19 days of data from the SK-IV era (Oct.
2008 - Apr. 2017). With an expected number of background events of 0.30 ± 0.13,
one candidate event was found. The number of observed events was consistent withthe expectation from background and a partial lifetime limit of τ (p → ν̄K + ) >1.7 · 1033 years was obtained at the 90% confidence level.

iii

A meus paisiv

Contents
1 Introduction
1.1 Unification and the Standard Model . . . . . . . . . . . . . . . . . . .
1.2 Grand Unified Theories . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3 Supersymmetry and Grand Unified Theories . . . . . . . . . . . . . .
1.4 Proton Decay . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.5 Proton Decay Searches . . . . . . . . . . . . . . . . . . . . . . . . . .
1.5.1 Proton Decay in Super-Kamiokande . . . . . . . . . . . . . . .

1
1
2
3
4
6
8

2 The Super-Kamiokande Detector
13
2.1 The Cherenkov Effect . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.2 Detector Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.3 SK Phases . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.4 Photomultiplier Tubes . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.5 Water and Air Purification Systems . . . . . . . . . . . . . . . . . . . 20
2.6 The Data Acquisition System . . . . . . . . . . . . . . . . . . . . . . 20
2.7 Calibration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.7.1 Water Properties . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.7.2 Relative PMT Gain Calibration . . . . . . . . . . . . . . . . . 22
2.7.3 Absolute PMT Gain Calibration . . . . . . . . . . . . . . . . . 24
3 Event Simulation25
3.1 Simulation of Proton Decay Events . . . . . . . . . . . . . . . . . . . 25
3.2 Simulation of Atmospheric Neutrino Events . . . . . . . . . . . . . . 27
3.2.1 Atmospheric Neutrino Flux . . . . . . . . . . . . . . . . . . . 27
3.2.2 Neutrino Interactions . . . . . . . . . . . . . . . . . . . . . . . 28
3.2.3 Elastic and Quasi-Elastic Scattering . . . . . . . . . . . . . . . 29
3.2.4 Single Meson Production . . . . . . . . . . . . . . . . . . . . . 30
3.2.5 Deep Inelastic Scattering . . . . . . . . . . . . . . . . . . . . . 32
3.3 Simulation of Neutrino Events with GENIE . . . . . . . . . . . . . . 32
3.4 Detector Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
4 Event Reconstruction
35
4.1 The Likelihood Function . . . . . . . . . . . . . . . . . . . . . . . . . 35
4.1.1 Predicted Charge . . . . . . . . . . . . . . . . . . . . . . . . . 36
4.1.2 Unhit Probability and Charge Likelihood . . . . . . . . . . . . 37
4.1.3 Time Likelihood . . . . . . . . . . . . . . . . . . . . . . . . . . 38
4.2 The Sub-Event Algorithm . . . . . . . . . . . . . . . . . . . . . . . . 39
v

4.3
4.4
4.5

4.2.1 The Peak Finder Algorithm . . . . . . . . . . . . . . . . . . .
The Single-Ring Fitter . . . . . . . . . . . . . . . . . . . . . . . . . .
The Multi-Ring Fitter . . . . . . . . . . . . . . . . . . . . . . . . . .
The Prompt-γ Fitter . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.5.1 The Seeding Algorithm . . . . . . . . . . . . . . . . . . . . . .
4.5.2 Performance of the µγ-fitter . . . . . . . . . . . . . . . . . . .
4.5.3 Vertex Position . . . . . . . . . . . . . . . . . . . . . . . . . .
4.5.4 Momentum . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.5.5 Time Difference . . . . . . . . . . . . . . . . . . . . . . . . . .
4.5.6 Vertex Direction . . . . . . . . . . . . . . . . . . . . . . . . .

40
40
42
42
43
44
45
49
51
53

5 Proton Decay Search
57
5.1 Signal and Background . . . . . . . . . . . . . . . . . . . . . . . . . . 57
5.2 Efficiency and Expected Number of Background . . . . . . . . . . . . 57
5.3 Samples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
5.4 Selection Criteria . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
5.4.1 Ring Counting . . . . . . . . . . . . . . . . . . . . . . . . . . 61
5.4.2 Particle Identification . . . . . . . . . . . . . . . . . . . . . . . 62
5.4.3 Number of Sub-Events . . . . . . . . . . . . . . . . . . . . . . 64
5.4.4 Muon Momentum . . . . . . . . . . . . . . . . . . . . . . . . . 66
5.4.5 Distance Between the Muon Decay and the Michel Electron
Vertex Positions . . . . . . . . . . . . . . . . . . . . . . . . . . 69
5.4.6 Time Difference Between the Muon and the Gamma Vertices . 74
5.4.7 Gamma Momentum . . . . . . . . . . . . . . . . . . . . . . . 76
5.4.8 Reconstruction of Events with no De-excitation Gamma . . . 78
5.5 Summary of Selection Efficiencies and Expected Number of Background Events . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
6 Systematic Uncertainties
83
6.1 Emission Probabilities of De-excitation Gamma Rays . . . . . . . . . 84
6.2 Atmospheric Neutrino Flux and Cross-Section Models . . . . . . . . . 84
6.3 Kaon Production in NC Neutrino Interactions . . . . . . . . . . . . . 85
6.4 Selection Criteria . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
6.4.1 Fiducial Volume . . . . . . . . . . . . . . . . . . . . . . . . . . 87
6.4.2 Number of Rings . . . . . . . . . . . . . . . . . . . . . . . . . 88
6.4.3 Particle ID . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
6.4.4 Number of Michel Electrons . . . . . . . . . . . . . . . . . . . 89vi

6.4.5

6.5
6.6

Distance Between the Muon Decay and the Michel Electron
Vertex Positions . . . . . . . . . . . . . . . . . . . . . . . . . .
6.4.6 Energy Scale . . . . . . . . . . . . . . . . . . . . . . . . . . .
Hybrid Sample . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Systematic Uncertainties Summary . . . . . . . . . . . . . . . . . . .

90
91
91
93

7 Data Results
95
7.1 Sidebands . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95
7.2 Data Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99
8 Results and Discussions
105
8.1 Lifetime Limit Sensitivity . . . . . . . . . . . . . . . . . . . . . . . . 105
8.2 Lifetime Limit Result . . . . . . . . . . . . . . . . . . . . . . . . . . . 106
8.3 Candidate Event . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107
Appendix A Lifetime Limit Calculations

vii

110

List of Figures
1.1
1.2
1.3
1.4
1.5

Unification of coupling constants. . . . . . . . . . . . . . . . . . . . .
4
Proton decay to e+ π 0 . . . . . . . . . . . . . . . . . . . . . . . . . . .
5+
Proton decay to ν̄K . . . . . . . . . . . . . . . . . . . . . . . . . . .
6+ 0
Event display for a p → e π simulated event in SK. . . . . . . . . .
9Momentum vs Mass distribution for proton decay and atmospheric ν
simulated events. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
9
1.6 Schematic representation of proton decay to ν̄K + . . . . . . . . . . . . 11
1.7 Hit timing distribution for proton decay with a prompt-γ. . . . . . . 11
2.1 Cherenkov wave fronts and the Cherenkov opening angle. . . . . . . . 14
2.2 Event display for showering and non-showering rings in SK. . . . . . . 16
2.3 Schematic representation of the SK tank . . . . . . . . . . . . . . . . 17
2.4 Schematic drawing of the 20-inch ID PMT in SK. . . . . . . . . . . . 19
2.5 Quantum efficiency of the ID PMTs in SK. . . . . . . . . . . . . . . . 19
2.6 Schematic representation of the laser injection system. . . . . . . . . 23
3.1 Atmospheric neutrino flux in SK. . . . . . . . . . . . . . . . . . . . . 28
3.2 Cross section of CC neutrino events. . . . . . . . . . . . . . . . . . . 29
3.3 Incoming amospheric neutrino energy spectrum. . . . . . . . . . . . . 30
4.1 The unhit probability correction of fiTQun. . . . . . . . . . . . . . . 37
4.2 FiTQun’s normalized charge pdf for different values of predicted charge. 38
4.3 FiTQun’s electron and muon PID as a function of electron momentum. 41
4.4 Event displays for a 235.5 MeV/c muon and a 6.3 MeV/c photon in SK. 44
4.5 Reconstructed and true muon vertex position using the µγ-fitter. . . . 45
4.6 Reconstructed and true muon vertex position in cylindrical coordinates using the µγ-fitter. . . . . . . . . . . . . . . . . . . . . . . . . . 46
4.7 Reconstructed and true muon vertex 2D position using the µγ-fitter.
46
4.8 Resolution of muon vertex position using the µγ-fitter. . . . . . . . . 47
4.9 Resolution of muon vertex position for events with no gamma. . . . . 48
4.10 Kaon track length. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
4.11 True µ and γ momentum distributions. . . . . . . . . . . . . . . . . . 49
4.12 Kaon lifetime for kaon decay in flight. . . . . . . . . . . . . . . . . . . 50
4.13 Reconstructed µ and γ momentum distributions. . . . . . . . . . . . . 50
4.14 Reconstructed minus true µ and γ momentum distributions. . . . . . 51
4.15 Resolution of µ and γ momentum distributions. . . . . . . . . . . . . 51
4.16 True and reconstructed ∆t distributions. . . . . . . . . . . . . . . . . 52
4.17 Exponential fit for the true and reconstructed ∆t distributions. . . . . 52
4.18 Reconstructed minus true ∆t distribution. . . . . . . . . . . . . . . . 53
viii

4.19 Angle between true and reconstructed µ and γ directions. . . . . . . . 54
4.20 Reconstructed and true opening angle between the µ and the γ. . . . 54
4.21 Opening angle vs true ∆t. . . . . . . . . . . . . . . . . . . . . . . . . 55
4.22 Reconstructed and true opening angle between the µ and the γ forevents with true ∆t > 10 ns. . . . . . . . . . . . . . . . . . . . . . . . 56
5.1 Number of rings distribution for signal and background. . . . . . . . . 62
5.2 Number of rings distribution for signal and background satisfying entire selection criteria. . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
5.3 PID distribution for signal and background. . . . . . . . . . . . . . . 63
5.4 PID distribution for signal and background satisfying entire selection
criteria. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
5.5 Number of sub-event distribution for signal and background. . . . . . 65
5.6 Number of sub-event distribution for signal and background satisfying
entire selection. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
5.7 Reconstructed µ momentum for signal and background for the 2-nse
sample. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
5.8 Entire momentum spectrum of subGeV background muons. . . . . . . 67
5.9 Reconstructed µ momentum for signal and background satisfying entire selection criteria. . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
5.10 Reconstructed µ momentum for signal and background for the 3-nse
sample. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
5.11 Reconstructed ∆X distribution for signal and background. . . . . . . 69
5.12 Reconstructed ∆X distribution for signal and background satisfying
the entire selection criteria. . . . . . . . . . . . . . . . . . . . . . . . 70
5.13 Reconstructed ∆X distribution for NC protons. . . . . . . . . . . . . 70
5.14 Reconstructed background momentum and direction for NC protons.
72
5.15 Vertex displacement that happens when a proton ring is mis-reconstructed
as a muon. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
5.16 Mis-reconstruction of the vertex position for NC protons. . . . . . . . 73
5.17 ∆X distribution for signal and background in the 3-nse sample. . . . 73
5.18 ∆t distribution for signal and background. . . . . . . . . . . . . . . . 74
5.19 ∆t distribution for CCQE events. . . . . . . . . . . . . . . . . . . . . 75
5.20 ∆t distribution for signal and background satisfying the entire selection criteria. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
5.21 ∆t distribution for signal and background in the 3-nse sample. . . . . 76
5.22 Reconstructed γ momentum for signal and background for the 2-nse
sample. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77

ix

5.23 Reconstructed γ momentum for signal and background satisfying theentire selection criteria. . . . . . . . . . . . . . . . . . . . . . . . . . . 77
5.24 Reconstructed γ momentum for signal and background for the 3-nse
sample. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
5.25 Reconstructed opening angle between the µ and the γ vs gamma momentum. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
5.26 Reconstructed opening angle between the µ and the γ vs gamma momentum (no γ). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
6.1 Kaon production cross section measured by MINERvA. . . . . . . . . 86
6.2 Kaon kinectic energy in background events. . . . . . . . . . . . . . . 87
6.3 Distribution of dwall for atmospheric neutrino MC. . . . . . . . . . . 88
6.4 Ring counting likelihood uncertainty estimation. . . . . . . . . . . . . 89
6.5 Particle ID likelihood uncertainty estimation. . . . . . . . . . . . . . 90
6.6 Example of merged files for hybrid sample. . . . . . . . . . . . . . . . 93
7.1 Number of rings distribution for data and atmospheric ν MC (sideband). 96
7.2 PID distribution for data and atmospheric ν MC (sideband). . . . . . 96
7.3 Number of sub-event distribution for data and atmospheric ν MC
(sideband). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
7.4 µ momentum distribution for data and atmospheric ν MC for the
2-nse sample (sideband). . . . . . . . . . . . . . . . . . . . . . . . . . 97
7.5 µ momentum distribution for data and atmospheric ν MC for the
3-nse sample (sideband). . . . . . . . . . . . . . . . . . . . . . . . . . 98
7.6 ∆X distribution for data and atmospheric ν MC for the 2-nse sample(sideband). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98
7.7 ∆t distribution for data and atmospheric ν MC. (sideband) . . . . . . 99
7.8 γ momentum distribution for data and atmospheric ν MC (sideband). 99
7.9 Number of rings distribution for data and atmospheric ν MC. . . . . 100
7.10 PID distribution for data and atmospheric ν MC. . . . . . . . . . . . 100
7.11 Number of sub-event distribution for data and atmospheric ν MC. . . 101
7.12 µ momentum distribution for data and atmospheric ν MC for the
2-nse sample. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101
7.13 µ momentum distribution for data and atmospheric ν MC for the
3-nse sample. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102
7.14 ∆X distribution for data and atmospheric ν MC for the 2-nse sample. 102
7.15 ∆X distribution for data and atmospheric ν MC for the 3-nse sample. 103
7.16 ∆t distribution for data and atmospheric ν MC. . . . . . . . . . . . . 104
7.17 γ momentum distribution for data and atmospheric ν MC. . . . . . . 104
8.1 Event display of the data candidate event. . . . . . . . . . . . . . . . 109x

8.2

Hit time distribution for the data candidate event. . . . . . . . . . . . 109

xi

List of Tables
2.1
3.1
5.1
5.2
5.3
6.1
6.2
8.1

Definition of the event triggers in SK. . . . . . . . . . . . . . . . . . . 21
Emission probabilities of the de-excitation gamma. . . . . . . . . . . 27
Summary of signal efficiencies for the 2 sub-event sample. . . . . . . . 81
Summary of signal efficiencies for the 3 sub-event sample. . . . . . . . 81
Summary of neutrino interaction modes for background. . . . . . . . 82
Summary of signal efficiency systematic uncertainties. . . . . . . . . . 94
Summary of background systematic uncertainties. . . . . . . . . . . . 94
Lifetime limit results. . . . . . . . . . . . . . . . . . . . . . . . . . . . 107

xii

Acknowledgements

Primeiramente eu quero agradecer a meus pais pelo apoio que recebi durante
toda minha vida, e em particular durante o tempo do meu doutoramento em que eu
estive fora de casa. Muito obrigado por tudo!
I would also like to thank my advisors, professors Michael Wilking and Chang
Kee Jung. For giving me the chance to work in such an exciting topic in a wonderfulexperiment. For all the time dedicated to discussions that made me (I hope) a
better physicist and person. For providing the conditions for me to travel to Japanand advance in my research. And overall, for providing such an excellent work
environment at Stony Brook University.
Thanks to my committee for participating in the entire process, including prof.
Shrock that could not be in my defense but provided me great knowledge during
all these years. Special thanks to professors Rafael and Verena that helped meimmensely during the final stretch and for their friendship. I met Rafael in the very
first day I came to Stony Brook and he has been a great friend and mentor eversince, muito obrigado!
Thanks to all the people I met and friends I made along these years. Gustavo,
Jose and Luana that I met from the beginning and thanks to them I survived theearly years. I am certain they will be lifelong friends no matter the distance. The
entire Spanish crew, in particular Inigo, Sara and Douglas, whose friendship helpedme all these years. I can not express how thankful I am to Sara and Douglas for thehelp and support I got these last few months before my defense.
Thanks to all NN group members for the friendship and mentoring. Thank you,professors Clark McGrew and Chiaki Yanagisawa for comments and suggestions,and thanks Chiaki for making my experience in Osawano even better. Thanks to
my lab mates that became friends along the years, Jay, Zoya, Xiaoyue, Yue and
Kevin. Thanks to all the postdocs that worked with me, Jose, Cristovao and Guang.
Specially, to Jose and Ruth and their family for the friendship and all the delicious
food! Obrigado, Jose e Cris pelas muitas horas de discussoes, conselhos e conversas!
Thank you to my collaboration colleagues and friends. ATMPD conveners, professors Ed Kearns and Shiozawa-san, for the guidance and support. Thank you
Miura-san and Mine-san, for all the help, guidance and mentoring I received every
week in our PDK meetings. Thanks to all the friends I met in Kamioka, Lluis, Lauren, Akutsu-san, Kai, Roger and many others with whom I could share many fun
moments, either inside the mine, in the field playing soccer or singing (very badly)
in karaoke rooms!
xiii

1

Introduction

1.1

Unification and the Standard Model

The standard model (SM) of particle physics (for a comprehensive description, see
Reference [1]) is based on gauge theories of the strong, weak and electromagnetic interactions among all known elementary particles. It has had enormous success in theagreement between its predictions and most experimental results. It is worth notingthat the SM is an empirical model, constructed using many different experimental
observations. Two of these empirical observations are of interest here.
The first observation is that the symmetry group of the SM is given by theexternal product of each gauge group separately as in Equation 1.1:
SU (3)C × SU (2)L × U (1)Y ,(1.1)

with each interaction having its own coupling constant. The weak and electromagnetic forces are related by the spontaneous symmetry breaking sector of the model,
but they still have different couplings. Thus, there are in total, three parameters
related to the coupling strength of each interaction that need to be experimentallymeasured.
The second observation is related to the strong interactions. The fermion content
of the model can be separated into quarks and leptons, depending on their interactions with gluons, the carriers of the strong force. Quarks participate in all the forces
including the strong force and are grouped into triplets of color following the SU (3)Csymmetry, while leptons are singlets under SU (3)C and do not interact strongly.
The matter content of the SM is split into different SU (2)L doublets, one forquarks and another for leptons as in Equations 1.2 and 1.3, respectively.
     
u
c
t
,
,
,
(1.2)
d
s
b
and

     
e
µ
τ
,
,
,
νe
νµ
ντ

(1.3)

with each column representing a family (or generation) of the known elementary
particles.
In conclusion, the SM does not attempt to unify all three forces into a single
interaction and there is no interaction vertex capable of mixing a quark and a lepton.
Both of these features are in agreement with current experimental results. Despite
this agreement, Section 1.2 presents hints that point to an unification of the coupling1

constants. Physicists have built extended models where the gauge interaction comesfrom a single simple Lie group (instead of the product of groups like in the SM). In
these models, the three SM interactions are simply different aspects of a single force
spontaneously broken at a very high energy scale. Due to this unified character,quarks and leptons are no longer distinguishable and lepton-quark mixing is allowed.

1.2

Grand Unified Theories

Despite the great success of the SM to explain many experimental results, it does
not address the unification of forces. The goal of Grand Unified Theories (GUTs) [2]is to describe the SM interactions in a unified way. The first GUT was developed in1974 by Pati and Salam [3]. The model is based in the SU (4)C × SU (2)L × SU (2)Rsymmetry, treating leptons as a fourth color, thus unifying quarks and leptons. As
noted in the previous section, total unification is realized when the gauge symmetry
is a simple Lie group. This is not the case for the Pati-Salam model, thus strictly
speaking it is not a unified theory of the SM forces.
Still in 1974, Georgi and Glashow wrote the first GUT based on a simple gauge
group, SU (5) [4]. To be able to embed the SM group into a larger one, a groupwith rank 4 or higher1 is needed and the minimal simple group that satisfies thisrequirement is SU (5). The theory has one universal coupling constant, αG as opposedto three as in the SM.
In GUTs, quark and lepton fields are not in separate representations of the gauge
groups any longer. Instead there is mixing in the irreducible representations. In
SU (5) for example, the two representations are 5̄ = [dR , LL ] and 10 = [QR , uR , eR ].
Equation 1.4 shows both representations with all the fields written explicitly.


 c 
0
uc3 −uc2 u1 d1
d1
−uc3 0
 dc2 u 2 d2 
uc1

 c
 c 
c


u 3 d3 (1.4)
5̄ = 
, d−3  , and 10 =  u2 −u1 0
−u1 −u1 −u3 0 ec e −d1 −d2 −d3 −ec 0−νewhere i = 1, 2, 3 represents the color indices and c is charge conjugation.
The gauge boson sector needs new bosons, since SU (5) requires 24 bosons andthere are only 12 in the SM2 . The adjoint representation of SU (5) (24) can beconstructed by placing the 3 × 3 block of SU (3) on the top left corner and the 2 × 2
block of SU (2) × U (1) on the bottom right. This way, the usual SM interactions
1
2

The SM group is rank 4 = 2+1+1 from SU (3), SU (2) and U (1), respectively.
Each SU (N ) group has N 2 − 1 gauge bosons and U (1) has one.

2

can be described without any mixing between quarks and leptons. The off-diagonalterms on the other hand, are filled with the new 12 gauge boson terms. These terms
will contribute to lepto-quark interactions leading to proton decay, explained in moredetail in Section 1.4.
Besides unification of all gauge interactions, SU (5) explains why electrical charge
is quantized by imposing the traceless condition on its generators. The charge operator Q is a linear combination of SU (5) generators, thus it is traceless3 . When Q
acts upon the 5̄ , one gets:
1
3Q(dc ) + Q(e− ) + Q(ν) = 0 ⇒ Q(dc ) = − Q(e− ),
3

(1.5)

where it becomes explicit that the charge of down quarks are a fraction of the chargeof the electron. While SU (5) has many virtues, it also predicts quantities that areinconsistent with experimental observation. Notably, the Weinberg mixing angle [5]and the lifetime of the proton [6, 7, 8, 9]. Other GUTs have been conjectured along
the past four decades, most of them sharing the same virtues that make GUTs so
attractive. Among the shared features, the prediction of nucleon decay stands out
as the only experimental evidence that can be pursued by current and near future
experiments.

1.3

Supersymmetry and Grand Unified Theories

Supersymmetry (SUSY) remains as one of the most appealing extensions of the SM.
It is beyond the scope of this brief introduction to describe most predictions of SUSY,
but two are directly connected to the present discussion (for a proper introduction
to supersymmetric theories, see Reference [10]).
One of these predictions is related to the unification of the gauge coupling constants. Before the constants were precisely measured by experiments such as theLarge Electron-Positron Collider (LEP) [11], it appeared as they would nicely converge to the same point at a very large energy scale, named the unification (orGUT) scale. Once more precise measurements were made, it was observed that thisunification was not achieved as can be seen in Figure 1.1, on the left.
The Minimal Supersymmetric Standard Model (MSSM) is an extension of the
SM that incorporates SUSY. Since new supersymmetric particles are part of themodel, it is necessary to include their contribution in the renormalization groupevolution equations. The new contributions were such that a new unification scale
3This is true for any generator of SU (N ).

3was predicted where the constants unify consistently with experimental observations
as in Figure 1.1, on the right.

Figure 1.1: Evolution of the inverse of the three gauge coupling constants in the SM(left) and in the (MSSM) (right) [12]. The dark spot at the unification on the right
represents model dependent corrections to the evolution equations.
The other prediction of supersymmetric theories happened when SUSY was incorporated in GUTs. Motivated by the unification of coupling constants, grandunification theories were extended to include supersymmetric particles. This newclass of GUTs is generally known as SUSY-GUTs and among other differences, they
typically predict different branching ratios for proton decay modes. The exact values
for lifetime and branching fractions depend on the particular model, but two mainmodes characterize the searches for proton decay.
When SUSY is absent in the theory, most GUTs predict the highest branchingratio of proton decay to be to p → e+ π 0 . The presence of SUSY in the theory
typically leads to the dominant decay mode p → ν̄K + . Section 1.4 explains in depth
the theoretical prediction of both modes and Section 1.5 presents an introduction of
how they can be experimentally observed.

1.4

Proton Decay

As discussed in the previous section, GUTs can be classified by the presence orabsence of SUSY in the theory. Both of these classes of theories predict proton
decay, although the exact lifetime and branching ratios depend on the details of each4

model. In general, non-SUSY GUTs favor the gauge mediated proton decay diagram
shown in Figure 1.2. As introduced in Section 1.2, new gauge bosons introduce
mixing between quarks and leptons. To preserve charge conservation, these bosons
have fractional charges of ± 43 and ± 31 .
This gauge mediated decay can be seen as an effective operator of dimension six,
analogous to the four-fermion theory before the discovery of the electroweak bosons.
In this mode, the partial lifetime of the proton is given byτp→e+ π0 ∼

MG4
,
m5p

(1.6)

where mp is the mass of the proton and MG is the mass of the X and Y gaugebosons. In most theories MG is on the order of the unification scale discussed in theprevious section. For minimal SU (5), this was about 5 · 1014 GeV, which led to a
predicted partial lifetime around 1031±1 years. This prediction has been excluded byexperimental results as discussed in Section 1.2.

Figure 1.2: Gauge mediated proton decay to positron and neutral pion. This mode
is typically favored in non-SUSY GUTs [13].
The presence of supersymmetric particles in the theory allows for another modeto have the largest branching ratio. As seen in Figure 1.1, the unification scale in
SUSY theories is higher, causing the mass of the gauge bosons to be around 1016 GeV
and τ (p → e+ π 0 ) → 1036 years. Figure 1.3 shows another diagram that leads to the
decay of the proton through the propagation of supersymmetric particles in the loop.
Despite being a loop diagram as opposed to the tree-level diagram from Figure 1.2,
this contribution can be much larger depending on the mass of the supersymmetricpartners.
In the scenario where the supersymmetric partners of the quarks and leptons
have a mass on the order of the electroweak scale, the only particle in the loop that
does not necessarily have an electroweak mass is the Higgsino (H̃), superpartner ofthe Higgs. The mass of the Higgsino is on the order of the GUT scale, causing the5diagram to be an effective dimension five operator instead of dimension six as before.
In this case, the proton partial lifetime is given byτp→ν̄K + ∼2
MG2 MSU
SY
,
5
mp

(1.7)

where again, mp is the proton mass and MG is on the order of the GUT scale. Thistime, the mass scales factorize and the light supersymmetric particles contribute with
2
MSU
SY , which is many orders of magnitude smaller than MG . Minimal SUSY SU (5)predicted a lifetime on the order of 1033 years for this mode.

Figure 1.3: Proton decays to anti-neutrino and charged kaon though a supersymmetric mediated loop. In most SUSY GUTs, this mode has the largest branchingratio.
As discussed in Section 1.2, many other theories have been proposed along thefive decades that follow the initial ideas of Pati and Salam followed by Georgi andGlashow. The next section describes how these predictions can be experimentallytested. In particular, the main topic of this thesis is discussed.

1.5

Proton Decay Searches

The previous section described some of the predictions made by GUTs for the lifetimeof the proton. As it was discussed, these predictions are usually on the order of1032−38 years, depending on the model. For an amusing comparison, the age of the
Universe is estimated to be about 1010 years, which gives an idea of how large thesepredictions are.
Particle decay processes are characterized by the following differential equation:tdN (t)= −αN ⇒ N (t) = N (t0 )e−αt ≡ N0 e− τ ,
dt

6

(1.8)

where N (t) represents the number of particles in a sample at a given time t and α issome proportionality constant. This equation shows that the number of particles inthe sample is decreasing with time at a rate proportional to the number of particles
itself. The solution is given in the first step, where N (t0 ) is the initial number ofparticles at time t0 . In the last step, N (t0 ) was defined as N0 and α as 1/τ .
For protons with τ larger than 1030 years, the observation time, t, is always muchtsmaller than τ . In that case e− τ can be very well approximated by 1− τt and Equation
1.8 can be rewritten ast

N (t) = N0 e− τ ≈ N0 −

N0 t
N0 t
⇒ ∆N (t) ≡ N0 − N (t) =
,
τ
τ

(1.9)

where ∆N (t) is the number of expected decays after time t.
This shows the importance of the quantity known as exposure defined as the initial number of particles in the sample times the observation time. For an observation
time, T , the exposure defined as λ is therefore:
λ ≡ N0 T.

(1.10)

If the number of expected decays is to be of the order of one event, it is necessary
to have an exposure similar to the lifetime. One possibility is to observe a proton
for more than 1032 years while waiting for it to decay. Another option is to observe
many protons at once for a more reasonable amount of time4 , in the hope that one
decay event can be seen.
The second approach was adopted by different scientists over the decades that
followed the initial ideas of unification5 . Many experiments were build to search for
proton decay, an observation that would have certainly shifted the SM paradigmof particle physics. Despite the lack of evidence for such new phenomenon, proton
decay still remains as one of the leading motivations for new experiments [14, 15].
The history of proton decay experiments is tightly connected to neutrino physics
and ultimately culminated in the discovery of neutrino oscillations. Given the vastand rich history of such experiments, a complete review will not be presented here
(for a comprehensive review, see Reference [16]). Instead, a more detailed overview
about a single experiment, that currently leads the proton decay search limits, will
be given.
4
5

Preferably, less than the age of the Universe.
The author is not aware of any attempt to observe a handful of protons for a long time.

71.5.1

Proton Decay in Super-Kamiokande

Super-Kamiokande (SK) is a large water Cherenkov detector, the full description of
the experiment is given in Chapter 2. This section discusses how a search for proton
decay events can be done in SK. As described in the previous section, the adopted
strategy to observe the decay of the proton is to build a large detector with many
protons and constantly monitor them. It is important to have a large exposure,defined previously as the product of the number of protons and the time duration of
the observation. The first is fixed by the volume of the SK tank, leaving the second
as the only degree of freedom to change the exposure. It is important to reduce deadtime in the detector to a minimum, so that the largest possible amount of data canbe collected. This observation time for which data can be recorded and analyzed isdefined as the livetime of the experiment.
The non-SUSY favored mode, p → e+ π 0 is the so called “Golden Mode” of water
Cherenkov experiments since all final state particles can be reconstructed with fairly
high efficiency. The positron being electrically charged creates a Cherenkov ringthat can be detected by photo-multiplier tubes and the neutral pion decays to two
photons, that produce one extra ring each. The photons are electrically neutral, but
as they propagate in the water, interactions such as e+ e− pair creation occur andthese particles produce Cherenkov light. Figure 1.4 shows a simulated e+ π 0 event.
The positron ring is located on the left and the two overlapping rings on the rightare from the photons.
Since all particles in the final state can be seen, it is possible to reconstruct the
total invariant mass and momentum of the event. For proton decay events, the total
invariant mass should be around the proton mass, while the total momentum should
be close to zero. Background events coming from atmospheric neutrino interactions
can mimic this type of event signature in the detector. One difference is that thesetype of events typically have a different phase space distribution of reconstructedinvariant mass and momentum. Using this information, it is possible to distinguishproton decay from atmospheric neutrino events.
Figure 1.5 shows the total momentum versus invariant mass distribution for proton decay and atmospheric neutrino simulated events. It can be seen that proton
decay events are in the region described above, centered around the proton mass withlow momentum6 . The distribution on the right shows the dominance of atmospheric
neutrinos in the low momentum-mass region, which allows for the separation of these
two types of events.
6

Effects such as Fermi momentum and correlated decay cause signal events to be reconstructed
outside the signal box.

8

Figure 1.4: Simulation of a p → e+ π 0 event in SK. The colored circles represent thetime of the hits observed by the photo-multiplier tubes that detect light coming fromthe Cherenkov rings of final state particles. The large ring on the left corresponds toto the positron and the two overlapping rings on the right are from the two photons
from the pion decay.

Figure 1.5: Total momentum vs invariant mass distribution for proton decay (left)
and atmospheric ν (right) simulated events. Free (light blue) and bound (dark blue)
protons are shown separately.

9For the SUSY-favored mode, p → ν̄K + , the situation is very different. Both final
state particles travel undetected in water Cherenkov experiments. The neutrino
travels through the detector without leaving any trace. The kaon, despite being
electrically charged, does not have enough kinetic energy to produce Cherenkov light.
This is because the initial energy of the process is given by the proton mass of about
938 MeV/c2 . The final state energy is the sum of the neutrino kinetic energy (itsmass is essentially zero) and the kaon total energy. With a mass of 494 MeV/c2 , thekaon kinetic energy is about 105 MeV, which is below the Cherenkov threshold of
257 MeV in water.
Not all hope is lost, since the kaon lifetime is short enough for it to decay insidethe detector. The largest branching fraction is to an anti-muon and a neutrino. The
neutrino is once again traceless, but this time, the muon does have enough energy to
produce Cherenkov light and leave a signal in the detector. Since the kaon undergoes
a two-body decay at rest (Chapter 5), the momentum of this muon is about 235.5
MeV/c. A single ring can be searched in SK, compatible with a muon with thismomentum.
Atmospheric neutrinos are once again the villains for this search, as it is certainly possible for a neutrino interaction to produce a muon in the final state with amomentum similar to the expected monochromatic signal. In fact, atmospheric neutrino events that produce a muon in the final state are abundant and it becomes very
challenging to separate proton decay signal from neutrino background on a statisticsbase.
The hero for this analysis comes from the realm of nuclear physics. If the decayof the proton happens inside an oxygen atom in the water, the remaining nucleus
can be left in an excited state. Once this nucleus transitions to its ground state, ade-excitation photon might be emitted. This whole process happens in the scale offemtoseconds, almost immediately after the proton decayed. The muon signal on theother hand, only becomes visible once the kaon decays a few nanoseconds later (thekaon lifetime is 12 ns). Figure 1.6 is a schematic representation of this process.
Performing a delayed time coincidence measurement between the monochromatic
muon of 235.5 MeV/c and the de-excitation photon allows for a strong discrimination
between signal and background events, since neutrino interactions do not present
this time difference between the particles. Figure 1.7 shows one example of the
distribution of the number of hits in the detector as a function of time. The firstpeak corresponds to the nuclear de-excitation photon, while the intense peak that
follows are the hits from the muon. This analysis is known as the prompt-γ channelof the ν̄K + mode. Chapter 5 describes this search in more details.
This thesis consists of the search for the SUSY-favored mode, p → ν̄K + in the10prompt-γ channel using an exposure of 177 kton·years of SK data.

Figure 1.6: Schematic representation of proton decay to ν̄K + . Immediately after the
proton decay, the de-excitation photon is emitted by the nucleus. The kaon travels
for a few nanoseconds before coming to rest and decay to an anti-muon.

Figure 1.7: Distribution of the number of hits in the SK detector for a proton decay
event with a prompt-γ present. The first peak corresponds to the nuclear photon,
while the intense peak that follows are the hits from the muon.
In the interest of completeness, it is necessary to precisely define the expression
11

“partial lifetime” used above. If a particle has multiple decay modes, then itsPtotal
decay width is given by the sum of each individual decay rate. In general, Γ = i Γi ,over all decay modes i. The predicted lifetime of the particle is given by the inverse of
the total decay width, Γ and partial lifetime is defined as the inverse of one particulardecay rate as in Equation 1.11.
τ≡1
1and τi ≡ ,ΓΓi(1.11)

where i represents the i-th decay mode. The goal of the present work is to infer alower limit on the partial lifetime τ (p → ν̄K + ) ≡ 1/Γ (p → ν̄K + ), in the absence ofa signal excess in the data results. The expressions “partial lifetime” and “lifetime”
are used as synonyms when no confusion should arise.

12

2

The Super-Kamiokande Detector

Super-Kamiokande (SK, Super-K) is a large water Cherenkov detector located deep
underground in Kamioka town - Gifu Prefecture, Japan. It consists of more thaneleven thousand photomultiplier tubes (PMTs) installed around the walls of a cylindrical tank holding 50 ktons of ultra-pure water. The experiment was first proposed
to search for nucleon decay and detect neutrinos from different sources. Initial neutrino analyses focused on solar and atmospheric neutrinos as well as astrophysicalneutrinos coming from supernovae. Later, SK also studied neutrinos coming fromaccelerators serving as the far detector for the K2K (KEK to Kamioka) experimentand currently, the T2K (Tokai to Kamioka) experiment.

2.1

The Cherenkov Effect

The Cherenkov effect is an electromagnetic phenomenon that describes the radiation
emitted as the coherent response of a dielectric medium to the passage of a fast
charged particle. When a charged particle travels in a dielectric medium such as
water, it polarizes the medium creating a local electromagnetic field. If the speed ofthis particle is faster that the phase velocity of light in that medium, the coherent
response of the medium emits the so called Cherenkov radiation. The phase velocity
of light in a medium of refractive index n is given by c/n, therefore the condition for
Cherenkov radiation is given by:
v = βc > c/n,

(2.1)

where v is the velocity of the particle. Light will then be emitted in a cone that
forms an angle θC with respect to the direction of travel, as can be seen in Figure
2.1. In a time interval t, the particle will travel a distance βct and light will travel
ct/n. Thus, the value of θC is given by:cos(θC ) = 1/nβ.

(2.2)

For water, nwater ≈ 1.33, so the opening angle of the radiation cone (sometimes
called the Cherenkov angle) is about 42◦ , for ultra-relativistic particles with β ≈ 1.
Using the condition given by Equation 2.1, it is possible to write the minimum
momentum or energy a particle must have to emit Cherenkov light as in Equations
2.4 and 2.5.




2
1
βmin
2
2
2
m =
m2 ,
(2.3)
pmin = (γmvmin ) =22
1 − βminn −1
13Figure 2.1: Cherenkov wave fronts and the Cherenkov angle. In a time interval t, theparticle travels a distance βct and light travels ct/n. The wave front forms a rightangle with the direction in which light is emitted.
orm≈ 1.14m,(2.4)
n2 − 1for n = nwater ≈ 1.33. Similarly we have,


m2
nm
2
≈ 1.52m,(2.5)
Emin =⇒ Emin = √2
1 − βminn2 − 1for n = nwater ≈ 1.33. Therefore, in water, a particle’s energy must be approximately
50% more of its mass for the Cherenkov effect to happen. For light particles suchas electrons this is always achieved in high energy events in SK. It is also satisfied
for Michel electrons7 . Muons and pions are heavier and not all of them can be seenin SK. Typical atmospheric neutrino charged current interactions will have enoughenergy to produce a muon above Cherenkov threshold, but muons from pion decays
are generally invisible. Protons are even more massive and require a lot more energy
to be seen. Nonetheless, as will be discussed later, the flux of atmospheric neutrinos
spans many decades of energy and in some interactions the proton leaves the nucleus
with enough energy to produce a Cherenkov ring.
In SK, the photons produced by Cherenkov radiation will travel the tank andsome will be detected by the PMTs. Using the pattern of hits, it is possible to infer
what kind of particle created that ring. As noted before, ultra-relativistic particles
create cones with θC ≈ 42◦ . Light particles like electrons will always satisfy thislimit, but heavier particles such as protons will have collapsed cones. As will bepmin = √7

Electrons coming from the decay of a muon.

14

described later, SK is able to distinguish showering rings, such as those created byelectrons, from non-showering, such as those coming from muons. Due to the factthat electrons are much lighter, they bounce more while traveling through the tankand therefore create the so called fuzzy rings when compared to muon rings whoseedges are much sharper. Figure 2.2 shows this difference between e-like and µ-likerings.

2.2

Detector Overview

The SK detector is located 1000 meters of rock below Mt. Ikenoyama or 2700 meterswater-equivalent mean overburden. It consists of a 50 kton tank of water opticallyseparated into a 32 kton inner detector (ID) and an outer detector (OD). Most of
the information used in event reconstruction comes from the ID while the OD is usedfor veto purposes only. It detects light activity coming primarily from cosmic raymuons, but also from γ-rays and neutrons produced in the rock. Figure 2.3 shows a
representation of the SK tank inside the mine.
The structure of the tank is a vertical cylinder of 41.4 m in height and 39.3 m
in diameter. Within the tank, a stainless steel framework of 55 cm separates the ID
and OD. It is located approximately 2 m from the walls of the OD and it holds the
structures where the photomultiplier tubes (PMTs) are mounted. Tyvek sheets [17]are placed in the space between PMTs to optically separate the ID and OD regions.
In the region facing the OD, white reflective Tyvek is placed to increase the chanceof light being detected in the veto region. In the region facing the ID on the other
hand, black plastic sheets are used to reduce reflected light to minimize the effect ofindirect light at the reconstruction stage.
Events that happen in the ID are detected by 11,129 20-inch photomultipliertubes (PMTs) pointing inwards, whereas the activity in the OD is detected by 1,885
outward-facing 8-inch PMTs. In this thesis we only consider events that are fully
contained inside the ID, i.e., with no activity in the OD.

2.3

SK Phases

SK started data taking in April 1996 and stopped for maintenance in July 2001,
with 1489.2 livetime days of data. This is known as the SK-I phase (or period orera). During maintenance, an accident happened when a PMT imploded during therefilling of water. Due to the vacuum inside the glass, water very quickly rushed insideand bounced back out creating a shock wave that destroyed 60% of the ID PMTs ina chain reaction. The remaining PMTs as well as spare PMTs were redistributed to15

Figure 2.2: Shower (top) and non-shower (bottom) rings in the fully contained data
set of SK. Each dot in the figure represents a PMT hit and the charge of the hit iscolor-coded. The small event display labelled OD shows the hits in the outer part ofthe detector. The small number of hits indicate that the event originated inside thetank, probably coming from a neutrino interaction instead of a cosmic muon thatentered from outside.

16

Figure 2.3: The Super-Kamiokande tank with its Inner and Outer detectors and theelectronic huts on top. The control room is located directly to left of the tunnel
entrance and the water purification system to the right.

17

have an effective photocathode coverage of 20% and data taking resumed in October
2002. This second period, SK-II, took data for 798.6 livetime days ending in October
2005.
The original configuration was restored in June 2006, starting the SK-III phasewith 518.1 livetime days until August 2008. Significant upgrades to the electronics
system were made until data taking resumed in September 2008. Notable improvements are the increase in Michel-electron tagging efficiency and the replacement of
the hardware trigger by a software trigger made possible by the improved electronics. This period is known as SK-IV and it represents more than 50% of all the datacollected by SK. It stoped in May 2018 for refurbishment with 3118.45 livetime days.
This thesis uses a partial dataset of SK-IV, until April 2017 with 2867.2 livetimedays or an exposure of 176.6 kton·years.

2.4

Photomultiplier Tubes

The ID uses 20-inch diameter Hamamatsu R3600 PMTs (Figure 2.4) developed incollaboration with Kamiokande scientists. Their photocathodes are coated with bialklai that has high sensitivity to Cherenkov light. The quantum efficiency peaks at
22% around the 360-400 nm region, as shown in Figure 2.5.
The 11,129 PMTs are evenly placed inside the ID and provide a 40% coverageof the tank. To prevent implosions such as the 2001 accident, each PMT is enclosed
in an acrylic cover. The transparency of the cover for photons at normal incidence
is higher than 96% for photons above 350 nm. The average transit time of a single
photoelectron is 2.2 ns. This time is achieved after a set of 26 Helmholtz coils wereplaced around the tank to reduce the effect of the geomagnetic field in the drift ofthe electrons. This reduces the geo-field of about 450 mG to about 50 mG.
In the OD, 1,885 8-inch R1408 Hamamatsu PMTs are evenly spaced facing outwards. Wavelength shifting plates are installed on the photocathodes of these PMTs,
which increases their efficiency by about 50% [18]. These 60 cm2 plates absorb ultraviolet Cherenkov light and reemit in the blue-green visible spectrum which the PMTs
are sensitive. However, the re-emission process degrades the timing resolution of theOD PMTs from 11 to 15 ns. In any case this is still desirable, because of the higher
detection efficiency and the fact that the main goal of the OD is to veto entering
particles and tag particles that left the ID. The white Tyvek sheets previously described help in this veto process because their high reflectivity increases the lightdetection efficiency even further. The top and bottom of the OD are optically separated from the barrel with white Tyvek to better identify particles entering/exiting
the OD region.
18

Figure 2.4: Schematic drawing of the 20-inch Hamamatsu PMT used in the ID [17].

Figure 2.5: Quantum efficiency of the SK ID PMT as a function of wavelength [17].

19

2.5

Water and Air Purification Systems

Given the dimensions of the SK tank, typical Cherenkov light travels large distances
before being detected by the PMTs. Therefore, a high water transparency is crucial
to allow detection of this light and to better reconstruct events. SK uses spring waterfrom the mine after multiple processes of filtering, sterilizing and degassing to remove
particles, bacteria and radioactive contaminants. The water is continuously purified
at a rate of 30 tons/hour (which corresponds to the entire tank volume in a period of
about 70 days). The purification process starts with 1 µm filters to remove dust andthen a heat exchanger to kill and suppress the growth of bacteria, keeping the water at
a stable temperature of 13◦ C. The last step uses reverse osmosis to remove dissolved
gases, particularly radon. This is especially important for low energy neutrino events
such as the ones produced in the Sun or supernovae. Radioactive processes initiatedby radon decay are the main source of background for these events.
The air inside the mine has naturally high radon levels. To decrease the contamination of the water with the radon air, filtered Rn-free air is continuously pumped
inside the mine from outside, from the so called Radon Hut. Also, the rock surrounding the SK tank is coated in a polyurethane material to contain the radon from therock and prevent it from being released in the air.

2.6

The Data Acquisition System

An event in Super-K starts with a charged particle producing Cherenkov light. This
light will then travel and reach the photocathode of one of the PMTs, producing aphotoelectron (p.e.) signal that is amplified by the dynode of the PMT resulting
in a current. The final product of this whole process is a set of times and charges
recorded by each PMT in the tank, which will then be used to infer what kind ofphysics event happened inside the detector.
In SK-IV, each PMT signal is first fed into a charge-to-time converter (QTC)
which was especially designed for SK [19]. When this signal is strong enough (above
a certain threshold), the QTC will integrate all the charge from the signal for the
next 400 ns and output a rectangular pulse. In SK, this is called a “hit”. This signal
is then digitized by a time-to-digital converter (TDC). The information in the signal
contains the time of the first photoelectron and the total number of photoelectrons
detected, which will finally be used in reconstruction for physics analyses.
The timing and charge resolutions of the QTC are 0.3 ns for 2 pC and about
0.2 ns up to 50 pC, respectively. The charge dynamic range of the QTC is 0.2 2500 pC, and the charge non-linearity is better than 1% for the overall dynamic
range. A software trigger is then applied to decide if a particular event should be20

recorded or not. The key quantity used in the trigger is the so called N200 variable
that measures the number of PMT hits in a sliding window of 200 ns. There are fivedifferent software triggers dedicated to different physics events, each with a different
threshold for N200 , they are summarized in Table 2.1. If an event is interestingenough to be recorded, the time window containing all the hit information relative
to that event depends on the physics case and it is also summarized in the table.
This work concentrates on events recorded by the High Energy (HE) trigger,
which is always satisfied by a 235 MeV/c muon. The relation between hits andenergy is non-linear, but for low energy events one can use a reference number of 5
hits per MeV. Thus, the HE trigger threshold is conservative given that most HE
physics analyses require at least 30 MeV of reconstructed visible energy.
Table 2.1: Definition of the event triggers: Super Low Energy (SLE), Low Energy(LE), High Energy (HE) and Special High Energy (SHE) triggers. The OD trigger
is based on OD hits alone.
Trigger Type
SLE
LE
HE
SHE
HE

2.7

N200 Threshold
34 → 3147
50
70 → 58
22 in OD

Event Time Window (µs)[-0.5, +1.0][-0.5, +35][-0.5, +35][-0.5, +35][-0.5, +35]

Calibration

The SK detector is calibrated using different sources, the most important are described in this section. The calibration analysis is used as input for the detector
simulation described in Section 3.4 and for data analysis. A complete description of
the detector calibration can be found in Reference [20].
2.7.1

Water Properties

It is necessary to measure some properties of the water in the SK tank in order tocorrectly simulate the propagation of Cherenkov light through the detector. The
main property is called water transparency related to the attenuation of light that
propagates in the tank. In particular, light can be absorbed or scatter when travelingthrough the water, the intensity of light as it travels a distance L can be modeled by21

Equation 2.6.
I(L) = I0 e−L/Latten ,(2.6)

where I is the intensity of light after traveling the distance L and I0 is the initialintensity when the light was created. Latten is the attenuation length of light in thewater. This parameter depends on the wavelength of light and it has to accountfor scattering and absorption. Two scattering processes are of most relevance in
SK, Rayleigh and Mie scattering. The former is relevant for the case when the
scattering object is much smaller than the wavelength of light and the latter isrelevant for scattering processes where the length scales are comparable, such as
some large impurity in the water with size of order 100 nm. Considering these two
scattering effects and absorption, the attenuation length parameter can be written
as in Equation 2.7
Latten (λ) =1αRay (λ) + αMie (λ) + αabs (λ)

,

(2.7)

where each parameter α corresponds to the three effects described above. In orderto determine each parameter, calibration data from a laser injector system was used.
Figure 2.6 is a schematic representation of the experimental setup. A laser beam
is injected vertically from the top of the SK tank and the the wavelength of thelaser can be changed to study the behavior of the α-parameters as a function of
wavelength. The apparatus is permanently mounted inside the detector so that thewater quality can be constantly monitored. The MC simulation is tuned to theobtained experimental data until good agreement is achieved. Measurements show
an attenuation length of approximately 120 m for light of 400 nm wavelength.
2.7.2

Relative PMT Gain Calibration

Each PMT in the SK detector has an individual high voltage (HV) supply such thatthe charge response of the PMTs for a given light intensity is approximately uniform
for all PMTs. In order to determine the individual PMT gains, pulsed laser light isinjected in the tank at a fixed position near the centre of the ID. The hits and thecharge of each PMT in the ID are measured twice by flashing the laser isotropicallywith two different intensities.
The first set of measurements is done for a high intensity laser flash, IH , suchthat all PMTs detect several photons at a time. The charge Qi at the i-th PMT is
proportional to the intensity IH and to the individual gain Gi . It also depends on theacceptance rate ai and the quantum efficiency i . The acceptance is a geometricalproperty of the PMT and depends on the position of the i-th PMT in the tank as
22

Figure 2.6: Schematic representation of the laser injector system in the SK tank formeasuring the water quality and reflectivity of the PMTs [20].
well as the incidence angle of the incoming photon. The entire dependence can be
written as in Equation 2.8.
Qi = ci × IH × ai × i × Gi ,(2.8)

where ci is some proportionality constant.
The second set of measurement happens for a flash with very low intensity IL , with
only a few PMTs getting a hit at a time. In this low intensity mode, a PMT detects
at most one photon per flash and the number of hits at each PMT is proportional
to the intensity as before but almost independent of the gain because of a low hit
discriminator threshold. In this case, the number of hits, Ni can be expressed as inEquation 2.9.
Ni = c0i × IL × ai × i ,(2.9)
where c0i is another proportionality constant. With everything else held the same between the two measurements, the relative individual variation of the i-th PMT gain,
Gi , can be obtained with the ratio between the two measurements as in Equation
2.10.
Qi
.
(2.10)Gi ∝Ni
The data obtained from these measurements show a standard deviation of 5.9%
across all PMT gains. The correction factor obtained from Equation 2.10 is used tomake the charge response of the PMTs uniform in data analysis.
23

2.7.3Absolute PMT Gain Calibration

For the absolute gain calibration it is necessary to obtain the single p.e. distribution
to determine the relation between measured charge and p.e.’s at each PMT. In orderto achieve single p.e. calibration, a low energy source is used, namely a Californium
252
Cf source, surrounded by a spherical mixture of nickel oxide and polyethylene.
The source emits neutrons that are thermalized by the polyethylene and absorbedby the nickel. Upon absorption by the nickel nuclei, a gamma ray of about 9 MeV
is produced. The emission yield is about 100 p.e. isotropically, which is low enoughfor most of the hit PMTs to produce single p.e. signals.
After correcting for the relative gain calibration described previously, the gain,or the ratio of the final PMT charge output to the p.e. input, can be measured.
The conversion factor between observed charge and number of p.e.’s obtained fromnickel calibration is 2.658 pC/p.e., which is the value used for detector simulation
described in Chapter 3.

24

3Event Simulation

The simulation of events in the SK detector is described here. The simulation process can be divided in two parts, namely, event generation and detector simulation.
This chapter describes these two parts for signal and background simulated events.
Event reconstruction is then applied to simulated or data events and it is described
separately in Chapter 4.

3.1

Simulation of Proton Decay Events

This analysis assumes equal probability for the decay of protons inside the SK detector. Protons inside the hydrogen nucleus (free) are assumed to have equal chance of
decay compared to protons inside the oxygen nucleus (bound). On the other hand,
the kinematics of the final state particles can be changed depending on the state of
the initial proton. Free protons are considered to be at rest and no interaction withother nucleons is simulated. Analyses that have a dedicated search for free protons
benefit from smaller systematic uncertainties and less expected background events,for example the e+ π 0 analysis of free protons shown in Figure 1.5. Simulation of the
decay of bound protons consider the effects of Fermi momentum, correlation withother nucleons, nuclear binding energy and meson-nucleon interactions for the final
state particles.
In this thesis, free and bound protons are automatically separated by the search ofa de-excitation gamma. The decay of a free proton is never followed by the emissionof a de-excitation gamma. Only the decay of protons inside the oxygen nucleus can
produce the signal signature of this analysis. As described in Section 1.5.1, if the
decay happens inside an oxygen nucleus, the remaining nucleus can be left in an
excited state and a de-excitation gamma might be emitted in the process. Table 3.1
summarizes the gamma ray emission probabilities and energies. The prompt gamma
emission process is based on Reference [21]. The position of the proton inside theoxygen nucleus is calculated following the Woods-Saxon nuclear density model [22].
Fermi momentum and nuclear binding energies are simulated according to measurements of electron scattering in carbon nuclei [23]. The nuclear binding energy
effect is simulated by modifying the proton mass. Ten percent of proton decay events
are assumed to have entangled wavefunctions with other nucleons [24]. These correlated decays cause the invariant mass of the final state particle to be smaller than the
proton mass, because a portion of the proton momentum is carried by the correlated
nucleon. Since this analysis does not reconstruct the proton mass, these effects haveno influence in the kinematics of the final state in this search. Correlated decay and25Fermi momentum certainly affect the momentum of the kaon, but as described in
Section 1.5.1, the kaon is below Cherenkov threshold and it is not observed. Figure
4.11 shows the true momentum distributions of muons after kaon decay. More than
89% of the kaons stop in the water and decay at rest.
Meson-nucleon interactions are also included in the proton decay simulation by
calculating the mean free path of the meson inside the nucleus [25]. For kaons, elastic
and inelastic scatterings are considered. Elastic scatterings change the momentum
and direction of the kaon, but as just discussed, this does not affect the final kinematics of the muon. Inelastic scattering is considered via charge exchange. The
estimated probability of charge exchange for kaons with the simulated momenta inproton decay events is 0.14%.
An interesting comparison between the two main proton decay modes can be
made. For e+ π 0 , bound protons are affected by all the nuclear effects described inthis section. In fact, the main source of signal efficiency loss results from pion-nucleus
interactions. Pions have a high probability of interacting before leaving the oxygen
nucleus, either by absorption or charge exchange. In both cases, the signal signature
is completely lost. In elastic scatterings, the charged pion leaves the nucleus witha different momentum, which affects the reconstruction of the invariant mass of the
proton. Only for the case when the pion leaves the nucleus without interacting, thesignal signature is unaffected. This causes signal efficiency loss and large systematic
uncertainties because of all the model dependent simulations. On the other hand,
free protons do not suffer any nuclear effect and produce a very clear signal in thedetector.
The SUSY favored mode analysis does not have the contribution of free protons
for a clear signature in the detector, since these are not accompanied by the emissionof a gamma. On the other hand, none of the nuclear effects described in this section
substantially affect the observed quantities in this analysis. There is a very small lossof signal efficiency caused by kaon charge exchange, but the other effects produce
no change on the muon kinematics. The only model dependent result used in thesignal simulation is the rate of emission of de-excitation gammas. In fact, this is thelargest source of systematic uncertainty in the signal efficiency of this analysis andit is described in Chapter 6.

26

Table 3.1: Summary of probabilities and energies of de-excitation gammas emitted
by the remaining nucleus after proton decay.
State Energy (MeV)
p3/2
6.3
9.9
p3/2
s1/2
7.03
s1/2
7.01
others3.5

3.2

Probability (%)
41
3
2
2
16

Simulation of Atmospheric Neutrino Events

A large number of neutrinos travel through the SK tank every second. Most of themwere originated in nuclear reactions in the Sun or in cosmic ray interactions with the
Earth’s atmosphere. Solar neutrinos are abundant, but their energy is much smaller
than the energy produced by a proton decay event and cannot produce the same
signature in the detector.
Atmospheric neutrinos are produced when a cosmic ray strikes the Earth’s atmosphere. Pions and kaons are created in these interactions and subsequently decay
to neutrinos and charged leptons. Neutrinos produced in these processes have energies that span many orders of magnitude, from tens of MeV to hundreds of TeV.
Despite the large number of neutrinos arriving in SK every day, the rate of observed
atmospheric neutrino events is about 8 per day.
The next sections describe the simulated atmospheric neutrino interactions withspecial emphasis to the ones capable of producing a signature similar to proton decay
events.
3.2.1

Atmospheric Neutrino Flux

The simulated atmospheric neutrino flux at SK follows the the Honda model [26,27]. The model calculates the propagation of cosmic rays through the atmosphereconsidering the effects of geomagnetic field and solar wind. Solar activity can change
the flux at 1 GeV by a factor of two during its variations between minimum andmaximum activities. The geomagnetic field introduces effects such as the up-down
and east-west asymmetries of the neutrino flux at SK. Other models are used forthe estimation of theoretical uncertainties, namely the FLUKA [28] and BARTOL[29] fluxes. Figure 3.1 shows the absolute flux used in this analysis in the region
where atmospheric neutrino interactions produce similar signatures to proton decay
27

events. The figure also shows the ratio between the adopted model and other flux
calculations used in the estimation of systematic uncertainties discussed in moredetail in Chapter 6.
The flux calculation is done assuming no neutrino oscillations. The effect ofoscillations is considered using an event-by-event reweighing method based on theneutrino type, direction and energy. The exact calculation of each weight is doneusing a two-flavor oscillation paradigm, in particular, the µ − τ oscillation according
to Equation 3.1.


1.27∆m2 L
2
2
,
(3.1)
P (νµ → ντ ) = sin (2θ23 ) sin
E
where θ23 is the mixing angle between the second and third neutrino mass eigenstates,∆m2 is the mass splitting between these eigenstates, L is the length (in km) of
the distance travelled by the neutrino from its creation in the atmosphere to theinteraction point inside SK and E is the neutrino energy (in GeV). The values for L
and E are obtained from the MC simulation and maximal mixing is assumed between
these eigenstates with 2θ23 = 90◦ and ∆m2 = 2.5 · 10−3 eV2 .

Figure 3.1: Atmospheric neutrino flux at SK as a function of energy (averaged overdirection) [27]. Absolute flux (left) for the four neutrino types and flux ratio (right)
between the Honda model (this work) and other flux predictions.

3.2.2

Neutrino Interactions

Neutrino interactions in the SK detector are simulated using the neutrino event
generator NEUT [25]. NEUT simulates how neutrinos interact with oxygen and hydrogen nuclei in the water. Interactions with electrons are neglected in atmospheric28

neutrino simulations since the cross section for these events is three orders of magnitude smaller then interactions with nuclei. Interactions are classified based on theoutgoing lepton or equivalently, the charge of the exchanged weak boson.
Neutrinos only interact through the weak force, which means only Z and W
bosons participate in the interactions. If an interaction is mediated by the Z 0 boson,the interaction is called a Neutral Current (NC) interaction and the outgoing lepton
is a neutrino. If the mediator is a W ± , the interaction is called Charged Current(CC) and the outgoing lepton is an electron or muon, generically represented by l± .
The lack of distinction between electrons and positrons is intentional, since there is
no external magnetic field in SK. Electrons and positrons create the same patternof Cherenkov rings and cannot be differentiated. The same is true for muons andanti-muons. The simulation of interactions between the final state particles and thenucleus is also made by NEUT. The next sections describe each simulated neutrino
interaction in this analysis. Figure 3.2 shows the cross section of CC interactions ofeach of the processes described in this section. Figure 3.3 shows the neutrino energy
spectrum for atmospheric neutrino events in the SK detector.

Figure 3.2: Cross section divided by energy as a function of energy for CC interactionsof neutrinos with a nucleon. The NEUT prediction for neutrinos (left) and antineutrinos (right) is compared to experimental data [30].

3.2.3

Elastic and Quasi-Elastic Scattering

Elastic interactions are defined as events in which a neutrino interacts with a nucleon
with some momentum transfer without producing new particles in the process. The
most general NC elastic interaction can be written as:
ν + N → ν + N,
29

(3.2)

Figure 3.3: Energy spectrum of the incoming neutrino in atmospheric neutrino eventsin the SK detector. The distribution peaks around 500 MeV and the spectrum
continues beyond the TeV scale, but it was truncated here for visualization clarity.
with N representing a proton or neutron and ν a neutrino or anti-neutrino. The same
particles are present in the initial and final states with no particle creation. For CC,
these interactions are called quasi-elastic (CCQE) because new particles must be
created. The outgoing lepton must be charged and the nucleon also changes, but foratmospheric neutrinos with energies of hundreds of MeV or more, the momentum
transfer is bigger than the lepton mass and the approximation is justified. Thus,
CCQE interactions have one outgoing charged lepton and can be represented by:
ν + N → l + N 0,

(3.3)

where N, N 0 represent different nucleons and l is a charged lepton.
More than half of all atmospheric neutrino events in the SK detector are CCQE
interactions. Most of these interactions are classified as single ring events and since
the outgoing lepton can be a muon, CCQE scattering is a potential background forthe proton decay search in this work. The discrimination between CCQE interactions
and proton decay events is done primarily by requiring a time difference between theµ-like ring and the de-excitation gamma as discussed in Chapter 5.
3.2.4

Single Meson Production

The simulation of single meson production in NEUT is done primarily through resonance production, where a baryon resonance produces a single meson in the final
30

state. Similarly to other interactions, the charge of the mediator determines the
charge of the the outgoing lepton for CC or NC interactions. Interactions with single
meson production can be summarized as:
ν + N → l + N ∗ → l + N 0 + m,

(3.4)

where N ∗ is the baryonic resonance that decays to a baryon N 0 and a meson m suchas a pion, kaon or eta. These interactions are simulated for W < 2 GeV/c2 , where
W is the invariant mass of the hadronic final state.
A similar interaction is modeled for coherent pion production, where the incident
neutrino interacts with the entire oxygen nucleus producing a pion in the process:ν + 16 O → l + 16 O + π,

(3.5)

where the outgoing lepton and meson are likely to be produced in the forward region
because of the low momentum transfer in this process.
Charged current interactions with single meson production are responsible forabout 30% of atmospheric neutrino interactions, while neutral current processes occur about 10% of the time. The most important background for the proton decay
analysis in this thesis is single meson production, in particular when a charged kaonis produced. Charged kaons can be produced in both CC and NC interactions but
CC processes can be easily distinguished from proton decay signal since there is an
extra lepton in the event.
Single kaon production in neutral current interactions constitute more than 80%of the expected number of background events in the search for proton decay to ν̄K + .
The production of kaons is associated with a lambda baryon in the process:ν + p → ν + Λ0 + K + ,(3.6)

where ν can be νe or νµ and the corresponding anti-particles. There is a similar
process where the initial nucleon is a neutron and the final state kaon is neutral.
In that case, the kaon might charge exchange to a charged kaon before leaving thenucleus. Each of these processes is very rare and occur only about 0.032% of thetime in all atmospheric neutrino interactions. And for the neutral kaon case, it still
has to charge exchange which makes the process even more unlikely.
Despite their extremely low probability, kaon production is the main source of
background for this analysis because it is completely irreducible in some cases depending on the lambda baryon decay mode. About two thirds of the time, it decays
to a proton and a negative pion and the rest of the time to a neutron and a neutralpion. In the second case, the π 0 can always be seen since it decays to two photons
31and these will produce extra rings in the event, independently of the decay modeof the kaon. On the other hand, most of the time the lambda decays to potentiallyinvisible particles. The proton is always below Cherenkov threshold since the lambda
and proton masses are similar and the proton is produced with low momentum. In
case the π − is also below Cherenkov threshold, the only visible particles will be thedecay products of the kaon producing the same signature as a proton decay event.
3.2.5

Deep Inelastic Scattering

In deep inelastic scatterings (DIS), the initial neutrino interacts with a nucleon’s
constituent quark instead of the entire nucleon. Such interactions happen when theneutrino has enough energy to probe the inner structure of the nucleon and asymptotic freedom is a good approximation for the constituent quarks. These processes
are simulated in NEUT for W > 1.3 GeV/c2 and becomes increasingly dominant
for multi-GeV interactions as seen in Figure 3.2. For W < 2 GeV/c2 , only pions
are considered as outgoing mesons while kaons and etas are also considered in theW > 2 GeV/c2 region. Production of hadronic final states in the high energy regionis simulated using PYTHIA/JETSET [31].
Deep inelastic events in CC interactions account for 5% of the total atmospheric
neutrino interactions, while NC events happen about 2.5% of the time. DIS is not
an important background for the proton decay analysis since the typical energy ofthese events is very high. Most events have multiple mesons in the final state anddiscrimination between proton decay and DIS events can be done by number of ringsand visible energy in the detector.

3.3

Simulation of Neutrino Events with GENIEAs discussed in Section 3.2.4, kaon production in NC interactions constitute the main
source of background events for the proton decay analysis in this thesis. In particular the largest contribution for this type of background event was the interaction
described in Equation 3.6, simulated in NEUT for events with W < 2 GeV/c2 ofmass of the final hadronic state. The interaction in Equation 3.6 is the most probable way of producing a charged kaon in atmospheric neutrino interactions with low
W , but it is not the only way. There is another process using another baryon isospin
multiplet as shown in Equation 3.7
ν + n → ν + Σ− + K + .

(3.7)

Despite the interaction in Equation 3.7 having a lower cross section (of orderhalf) than Equation 3.6, the sigma baryon always decays to nπ − . Similarly to the32lambda baryon case, the neutron is invisible and the π − can also be produced below
Cherenkov threshold producing an event with the same signature as in Equation
3.6 and proton decay events. This interaction is not simulated in NEUT, but it isincluded in another neutrino event generator called GENIE [32].
Understanding kaon production in neutrino interactions is a crucial part of thiswork. Thus, GENIE was chosen as the event generator for atmospheric neutrino
events that produce kaons in the final state, since it has a more complete model
than NEUT in the simulation of kaons. The simulation process for the events withkaons in GENIE follows the same procedure as the events simulated by NEUT.
The same atmospheric flux, the Honda model described in Section 3.2.1, is used.
Events with kaons in the final state were selected using the simulated informationof each interaction to reject events without kaons. Detector simulation was done in
the same way as events generated by NEUT, as described in Section 3.4. Similarly,
event reconstruction using the fiTQun algorithm was performed on events with kaons
simulated with GENIE and events with no kaons generated by NEUT as describedin Chapter 4.
In summary, the expected number of background events for the proton decay
search in this thesis was simulated using GENIE for all atmospheric neutrino interactions that produce kaons in the final state, while NEUT was used for all otheratmospheric neutrino interactions. Chapter 6 describes the estimation of systematic
uncertainty associated with kaon production in neutrino interactions, as well as allother systematic uncertainties associated flux and cross-section models described in
the previous sections.

3.4

Detector Simulation

Atmospheric neutrino events in the SK tank were generated using NEUT and GENIE
simulations as described in the previous sections. The output of these simulations
was a set of final state particles and their kinematics after leaving the nucleus. It
is then necessary to simulate how these particles propagate in the water producing
signals in the detector. Detector simulation was done using a custom software called
SKDETSIM based on the GEANT3 package [33].
Similar to the hadronic interactions in the nucleus, the simulation of hadronicinteractions in the water is divided in two momentum regions. For pions below
500 MeV/c, NEUT is used to simulate their propagation in the water, while all otherhadronic simulations are done using GCALOR [34]. Cherenkov photons produced bycharged particles in the tank are then simulated using the custom code developed for
the SK detector. The model uses scattering and absorption of light in water based
33

on the calibration measurements described in Chapter 2. Similarly, reflectivity of
the PMTs and the black Tyvek sheets are also considered in the model.
The electronics of the detector is also simulated to take into account the charge
and time responses of the PMTs. The final output of detector simulation is then aset of charge and time pairs for each PMT with the same data structure as observed
data, such that both can be reconstructed and analyzed in the same manner.

34

4

Event Reconstruction

As discussed in Chapter 2, each event in SK is a collection of charges and times
recorded for every single PMT inside the tank. In order to measure physics quantities
such as momentum, direction and position of particles inside the SK detector, it isnecessary to first reconstruct events. In other words, one has to infer the type ofparticles and its kinematic properties involved in an event, from the set of charge-timepairs recorded by each PMT8 .
The event reconstruction algorithm used in this analysis is a maximum likelihood
fitter called fiTQun. The algorithm is based on methods developed for the MiniBooNE experiment [35], but it has been developed for the first time for Super-K
with additional functionalities such as the multi-ring reconstruction.
Sections 4.1 to 4.4 describe the general methods used in fiTQun and Section
4.5 describes the dedicated fiTQun fitter developed exclusively for this analysis of
p → ν̄K + in the prompt-γ channel. The description of fiTQun in this chapter
is intended to be brief, but self-contained. For a more detailed description of thealgorithm see References [30] and [36].

4.1

The Likelihood Function

This section describes the fiTQun algorithm and the calculation of the likelihood
function. The discussion is for the case of a single ring, i.e, a muon or an electron
traveling in the tank and producing a Cherenkov ring, but the generalization for themulti-ring case is straightforward.
A single particle in the detector can be characterized by seven kinematic parameters and its particle ID (muon, electron, proton, etc) (PID). Given a PID, a ring isdefined by its vertex position and time (x, y, z, t) and its momentum, which can beseparated into magnitude and direction (p, θ, φ), where θ is the zenith angle withrespect to the vertical z-axis and φ is the polar angle in the xy-plane9 . Let x be thevector containing this set of 7 parameters describing the particle track. The goal offiTQun is to find the value of x that maximizes the likelihood function described
below, in other words, determine the x that best describes the data.
For a given event in SK, fiTQun calculates the likelihood function based on the
charge and time information of all the ID PMTs:
8

The PMTs measure charge and send the information to the DAQ (Section 2.6) that records thetime of the hit.
9
The z-axis is vertical and points up (main axis of the cylinder). The xy-plane is horizontal,
parallel to the top and bottom circumference walls of the tank. The origin is at the center, i.e., at
the middle of the detector height and at r = 0.

35

L(x) =unhit
Y

Pj (unhit|x)

j

hit
Y

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

(4.1)

i

where the first index j runs over all the PMTs that did not record a hit and thesecond, i, runs over the ones that did. For each PMT in the tank there are only two
possible outcomes in each event: either it does not register a hit or it does, in thatcase with a charge and time associated with it. Pj (unhit|x) is the probability thatthe j-th PMT does not register a hit and {1 − Pi (unhit|x)} is the chance that thei-th one does. If there is a hit on the i-th one, then fq (qi |x) is the charge probability
function (pdf) for that hit, and ft (ti |x) the time pdf. Thus, fq (qi |x)ft (ti |x) is theprobability density function for observing charge qi at time ti for the i-th PMT given
parameters x.
The product of all these terms is the likelihood function that fiTQun will use as
the optimization metric.
4.1.1

Predicted Charge

The calculation of the likelihood function can be simplified by introducing a parameter called the predicted charge, which is the mean number of photoelectrons at a
PMT given the set of parameters, x. Using this variable, the likelihood function canbe rewritten as:
L(x) =unhit
Y

Pj (unhit|µj )

j

hit
Y

{1 − Pi (unhit|µi )} fq (qi |µi )ft (ti |x),

(4.2)

i

where fq (qi |µi ) now only explicitly depends on the i-th PMT and electronics response,while Pi (unhit|µi ) depends indirectly on x through µi .
Thus, the likelihood calculation is made in two separate steps: (i) the calculationof predicted charge for each PMT for a given set x, and (ii) the likelihood calculation
based on the predicted charges.
The calculation of predicted charge has to take into account direct light as wellas indirect light, such as scattered or reflected light. In particular, µ is given by thesum of direct and indirect light and contribution from dark noise. In the case of amulti-particle hypothesis, µ is also summed over each particle’s predicted charge asin Equation 4.3.
X

sct
dark
µi =
µdir
,
(4.3)
i,n + µi,n + µin36Figure 4.1: The unhit probability P (unhit|µ) with (red) and without (blue) the
correction of the PMT threshold effect. The data points show the values obtained
from detector simulation [30].
where µi is the predicted charge in the i−th PMT, µdir and µsct represent the direct
and indirect light components respectively, and µdark is the contribution from darknoise (5.72 kHz in SK-IV). For a complete description of the calculation of predictedcharge see Reference [30].
4.1.2

Unhit Probability and Charge Likelihood

Given the total predicted charge calculated as in Equation 4.3, it is possible tocalculate the unhit probabilities and charge likelihoods. Since µ is the predicted
number of photoelectrons in a PMT, the actual number obeys a Poisson distribution
with mean µ and the probability for seeing zero photoelectrons is given by e−µ .
However, in order to take into account the effect of the PMT threshold, a third ordercorrection is applied:P (unhit|µ) = 1 + a1 µ + a2 µ2 + a3 µ3 e−µ ,(4.4)

where the coefficients were obtained from detector simulation. Figure 4.1 shows
the unhit probability with and without the correction applied compared to the true
values from detector simulation.
Similarly, the charge likelihood is obtained by generating photoelectrons following
Poisson statistics with mean µ directly at the PMTs in the detector simulation, and37

Figure 4.2: The normalized charge pdf fq (q|µ) for different values of µ. The data
points are taken from detector simulation and the lines are the polynomial fits usedin the evaluation of the likelihood in the fitter [30].
using the distribution of the measured charge for the hit PMTs. Figure 4.2 shows theobtained normalized charge likelihood distribution for different values of µ. Whenevaluating the likelihood for a given event, the distributions are interpolated in termsof the observed charge q and the resulting polynomial is evaluated at µ.
4.1.3

Time Likelihood

The time likelihood ft (ti |x) depends on x and the PMT position in a complicatedway. To simplify the calculation of this pdf, the residual time is used instead of thehit time ti . The residual time is given by:
smid
−
tres
i = ti − t −~ P M T − ~x − smid p̂
R
i

,
(4.5)
c
cn
where ~x, t are the vertex position and time, p̂ is the particle direction, smid is the~ iP M T is the position
midway point in the particle track (starting at the track vertex), R

38

of the i-th PMT and cn is the speed of optical photons in water10 . Thus, tres is theresidual hit time after subtracting the expected direct photon arrival time from theraw hit time, assuming all the photons are emitted when the particle is at the middleof its track.
Using this approximation, the time pdf is rewritten as:
dir res
sct resft (tresi ) = wft (ti ) + (1 − w) ft (ti ) ,

(4.6)

µdir

where w = 1−e1−e. The direct and indirect light residual time pdfs are repreµdir e−µsct
dirsct
sented by ft and ft , respectively. Both pdfs depend on the predicted charge µ
and the particle momentum.
To calculate ftdir , particle gun simulations are used for various values of particle
momenta and predicted charges and each distribution is fitted by a gaussian. The
final evaluation is done using interpolation of these fits. For the indirect pdf, asimilar procedure is employed, but the simulated distributions present a long tail on
the right (larger residuals) compared to the direct light distributions. This is becauselight reflects on the walls of the tank and in the PMT surfaces and these produce
hits later in time. Therefore, a different model is used to fit these distributions:( τ21e− 2σ2 , (τ < 0)
p
×
(4.7)
ftsct (tres
)
=
i
π
− γτ
τ
σ
+
2γ
e
,
(τ
>
0)
2
γ+1
where τ = tres − 25 ns, σ = 8 ns and γ = 25 ns were obtained as parameters of the
fit to the simulated distributions.

4.2

The Sub-Event Algorithm

FiTQun reconstructs any event in SK using the algorithm described in the previous
sections. For some events it is necessary to define different time windows, each with
their own PMT hits. This is useful for example, in the case of a muon that decays to
a Michel electron. Both particles are in the same event, but the detector will measure
the hits from the muon and the Michel with a time separation, allowing fiTQun to fit
each of these collections of hits independently. Each time window found by fiTQun
is called a sub-event. This section describes how fiTQun finds and treats different
sub-events.
10

We use n = 1.38 to calculate the group velocity of light in water. Note that this is different
from the typical value nwater ≈ 1.33 which is related to the phase velocity of light.

39

4.2.1The Peak Finder Algorithm

The sub-event algorithm starts by performing a quick fit searching for the position
and time of the vertex that maximizes the following goodness-of-fit metric:
hits
T2
X
− res,iG(~x, t) =e 2σ2 ,

(4.8)

iwith σ = 4 ns and Tres,i is the residual time of the i-th PMT given by:Tres,i = ti − t −~ iPMT − ~x
R
cn

,

(4.9)

~ PMT is the position of the i-th PMT, ~x is the vertex position of thewhere again, Ritrack and cn is the group velocity of optical photons in water. It is similar to theresidual time calculated in the time pdf, but here it is assumed that the light iscoming from a point source in the vertex position. When the values of ~x and t get
close to their true values, T res is distributed around zero resulting in large values forthe goodness G.
For each peak found by this method, fiTQun tries to assign a time window corresponding to that cluster of hits. The condition to define a time window is thatany PMT must have at most one hit per window. Once all the time windows arefound, fiTQun can perform the likelihood optimization as described in Section 4.1,
the vertex position and time found in the peak finder algorithm serve as the initial
seed for the final fit.

4.3

The Single-Ring Fitter

After finding all the sub-events in an event using the method described in Section
4.2, the single-ring fitter is run on each one of them.
The fitter assigns initial values for the track parameters as seeds for the minimization. These values are then used by MINUIT [37] to minimize11 the negative
logarithm of the likelihood (− ln L), described in Section 4.1.
For the electron hypothesis, the values of ~x and t found by the peak finder algorithm are used as the initial seeds. It is still necessary to provide initial values
11This is typical of numerical optimization problems. For stability purposes, it is better tominimize −f (x) than to maximize f (x). Also, since the values of the likelihood function are
typically very small, taking the logarithm of L improves numerical stability. But algebraically,
maximizing L and minimizing − ln L are equivalent problems.

40Figure 4.3: Likelihood separation of electron (left) and µ (right) particle gun events.
The vertical axis shows the likelihood ratio ln LLµe and the x-axis shows the reconstructed momentum in the electron-hypothesys.
 The solid line is one example ofmomentum-dependent separation criteria, ln LLµe > 0.2 · prece . Figure taken from
Reference [30].
for the momentum and direction of the track. The first step is to roughly estimate
the momentum using the total charge measured in the event. By keeping this valueof momentum fixed, a likelihood scan is performed over 400 points equally spaced
along the unit sphere. The direction (θ and φ) chosen is the one that maximizes thelikelihood. With this fixed direction, a new likelihood scan is done, this time varying
the momentum.
Once this set of seven initial parameters is found, MINUIT can perform theoptimization. When this is done, it returns the best set of parameters, x, thatminimizes − ln L. FiTQun’s final output consists of the maximum likelihood foundand the parameters of the fit for each hypotheses, for example, electron and muon.
It is the analyzer’s job to decide which hypothesis  best describes the event.
As an example, take the case where ln LLµe > 0. This indicates that the event
is more e-like than µ-like, in other words, it is more likely to be an electron thana muon. But the exact criteria for PID separation can be optimized differently for
different analyses. As it will be shown later, in this work the separation value chosen
is ln LLµe > 30, but generally, this ratio varies with momentum as it can be seen in
Figure 4.3.
For the muon hypothesis, it uses the electron’s best fit value of position, time anddirection as its initial seed. To determine the seed for momentum, a likelihood scan
41

with varying momentum is done, similarly to the electron case. After which the fullminimization of − ln L is done by MINUIT exactly the same way as before.

4.4

The Multi-Ring Fitter

FiTQun can also reconstruct events with more than one ring. This is essential foratmospheric neutrino analyses where multi-GeV events often have more than one
particle in the final state. The multi-ring algorithm is applied only to the first subevent and it can identify up to six rings.
The way fiTQun looks for extra rings is similar to the single-ring muon fit discussed earlier. It starts by using the same vertex and time of the ring found in the
single-ring hypothesis. The best direction of the second ring is found by arbitrarily
fixing the momentum at 50 MeV/c and a scan in the unit sphere is performed. Finally, fiTQun varies the momentum of the second ring keeping the other parameters
fixed at their best values. Once all of these initial seven parameters are found, MINUIT performs the optimization. Note that now there are fourteen parameters in
the fit, seven for each ring. But the vertex and time of both rings are assumed to be
the same, thus the number of free parameters in the fit is reduced to ten.
Once the two-ring fit is done, fiTQun has to decide if the extra ring is indeed atrue ring or if it is fake. It does that by comparing the likelihood before and after
adding the new ring. If the likelihood improvement is large enough, then the eventis classified as a two-ring event, otherwise it considers the new ring to be fake andreturns the single-ring result instead.
This iterative process is done up to six rings, always adding one ring at a timeand comparing the likelihood before and after to judge the quality of the fit. The
thresholds for classifying the new ring as true or not vary from sample to sample,
but the procedure is always the same. For more details on the multi-ring fitter, see
References [30] and [36]4.5

The Prompt-γ Fitter

As described in Chapter 1, the signal for p → ν̄K + in the prompt-γ channel iscomposed of a low-energy de-excitation gamma followed in time by a 235.5 MeV/c
muon. This signal presents two difficulties to be reconstructed with fiTQun, namely,
the energy of the gamma is too low to be found by the multi-ring fitter and the timedifference between the particles can not be taken into account using the standard
multi-ring algorithm from Section 4.4.
The likelihood of the single-ring µ-hypothesis and the two-ring hypothesis are very
42

similar, since the gamma produces very few extra hits in the detector. Therefore, themulti-ring algorithm can classify the second ring as fake. But even if fiTQun wereable to find the extra hits from the gamma and tag the event as two-rings, there
is still the issue that the multi-ring algorithm assumes that both rings are comingfrom the same vertex position and time. This is an issue since the biggest power toreject background in this search is to apply a time difference selection between thetwo rings.
Because of these difficulties, a new fitter was developed, dedicated to this search.
In this section we present the seeding algorithm for the new fitter (called the µγ or
pmg-fitter12 ) and its performance in signal Monte Carlo simulation (MC). From hereon, the definition of signal MC is the following:• There is a true de-excitation γ present in the event.
• The K + undergoes a two-body decay to a µ+ and a νµ .
If both conditions are met, the event is called a true µγ signal event.
4.5.1

The Seeding Algorithm

There are two ways that the µγ-fitter can be used, namely, one can force the position of the muon and the gamma vertices to be the same or to have them as freeparameters in the fit. In the first case, there are eleven free parameters (seven foreach track minus three for the position of the gamma). While the second case is thefull fourteen parameter fit. Since the γ-ring is very weak, it is hard to reconstruct its
position with the same resolution as the muon vertex position, so in this work only
the first option was used. All the results shown here are for the eleven-parametercase assuming the same vertex position for both the muon and the gamma.
The seeding algorithm for the µγ-fitter starts by using the result of the single-ringµ-hypothesis fit. For the seven parameters of the muon it uses the entire result of
the single-ring fit as the seed values. The gamma vertex position is also seeded at
the same point as the muon vertex position. Its momentum is seeded at 6 MeV/c
and a likelihood scan is performed to determine the best direction and vertex time
to be used as initial values.
Similarly to before, a likelihood scan in the unit sphere is performed using 400
equally spaced steps. This scan is repeated eleven times, each time choosing a different time offset between the muon and the gamma vertices.
12

For proton to mu gamma, which is clearly a jargon, since the de-excitation gamma does notcome directly from the decay of the proton, but no confusion should arise from this convention.

43

Figure 4.4: Event displays for a 235.5 MeV muon (left) generated at (5, 0, 0) m
pointing at (1, 0, 0) and a 6.3 MeV gamma (right) generated at (10, 0, 0) m pointing
at (1, 0, 0).
In summary, fourteen parameters are needed for the initial seed (seven for themuon and seven for the gamma tracks). The seeding algorithm is the following:• Set the single-ring µ-hypothesis result as the seed for the seven parameters ofthe µ-track.
• Set the γ vertex position at the same point as the µ vertex.
• Set pγ = 6 MeV/c for the γ-momentum seed.
• Perform a likelihood scan in the unit sphere for each value of ∆t in (2.5, 5.0,
7.5, 10.0, 12.5, 15.0, 20.0, 25.0, 30.0, 40.0, 50.0) ns.
• The (θ, φ, t) combination that minimizes − ln L is chosen to be the seed values.
Figure 4.4 illustrates the difference between the µ and the γ rings. In the muoncase one can clearly see the sharp-edge ring characteristic of muons. Whereas in the
gamma example, only a few scattered hits are found in the detector and a clear ringis not formed.
4.5.2

Performance of the µγ-fitterIn this section we evaluate the performance of the µγ-fitter by comparing reconstructed (Rec, Reco) quantities and the true MC information from simulation. Biases and resolutions of the algorithm are quantified here, both for the muon and the
gamma parameters.
44

4.5.3Vertex Position

As explained in Section 3.1, proton decay events are generated one meter away fromthe detector walls to take into account the effect of migration in and out of thefiducial volume, defined to be two meters away from the walls of the tank. Besides
that constraint, proton decay events are expected to be scattered uniformly acrossthe tank. Figure 4.5 shows the reconstructed and true values of the muon vertex
position. Since SK has a cylindrical geometry, it is easier to check uniformity incylindrical coordinates as in Figure 4.6 and Figure 4.7 that shows the vertex positionof simulated events as seen from above.

Figure 4.5: Reconstructed and true vertex position: x (top left), y (top right), z(bottom left).

45

Figure 4.6: Reconstructed and true vertex position: r2 (top left), φ (top right) andz (bottom left).

Figure 4.7: True (left) and reconstructed (left) radial vertex position: y versus x.
Figure 4.8 shows the distance between the true and reconstructed vertices of the46

muon and the gamma, in other words, the resolution of vertex position. It is useful
to compare the position resolution of the µγ-fitter with the single ring µ-fitter for thesame events and also for events with a mono-chromatic muon but no de-excitation
gamma present. This is shown in Figure 4.9.
Despite the energy of the de-excitation gamma being very low compared to themuon energy, its presence is enough to worsen the resolution of the single ring fitter.
The µγ-fitter on the other hand, takes into account the presence of the nuclear
gamma and its resolution is comparable to the single ring fitter in the pure muon
case.

Figure 4.8: Resolution of the muon vertex position (left), defined as the distance
between the reconstructed and true muon vertices. Similarly, for the gamma vertex
position (right).

47

Figure 4.9: Resolution of the muon vertex position using the single ring µ-fitter intrue µγ events (left) and true monochromatic muon events with no de-excitation
gamma (right).
As discussed in Section 4.5.1, the gamma vertex position is constrained to be
the same as the position of the muon. But despite this constraint, Figure 4.8 showsthat the resolution for the gamma vertex is comparable to the muon vertex. Thiscan be understood by looking at the typical distance travelled by kaons before they
decay or come to a stop in the water as in Figure 4.10. Most kaons travel about
13 cm before decaying, so the resolution of the gamma vertex position is dominatedby the resolution of the muon vertex position convoluted with the kaon track length.
A quick calculation shows an agreement√ between the observed resolution of about
30 cm and the expected resolution of 132 + 272 ≈ 30 cm.

Figure 4.10: Kaon track length - distance travelled by the kaon before decaying or
coming to rest (log scale on the right).

48

4.5.4

Momentum

Figure 4.11 shows the true distributions of muon and gamma momentum. The peakat 235.5 MeV/c in the muon momentum distribution corresponds to kaons thatdecayed at rest and the tails correspond to kaons that decayed in flight. Figure 4.12
shows the time kaons typically travel when decaying in flight, the majority of kaons
travel for at most about 2 ns before coming to a stop in the water and decaying at
rest. The gamma distribution shows the emission lines in N15 discussed in Chapter
3, the most important line corresponds to the 6.3 MeV emission and the second mostimportant lies at 3.5 MeV.
Figure 4.13 shows the reconstructed distributions of muon and gamma momentum. The muon momentum curve is fitted by a gaussian that was used in the selection
criteria of the analysis in Section 5.4.4. A small bias is observed in the central valueof the fit compared to the expected monochromatic value of 235.5 MeV/c. Figure
4.14 shows the difference between reconstructed and true momentum for the muonand the gamma. Section 5.4.4 also shows the behavior of the background momentum(defined in Chapter 5). Since the background distribution is essentially flat in thisregion, this small bias has no effect in signal efficiency or background rejection.
Figure 4.15 shows the momentum resolution defined as (Reco-True)/True momentum for the muon and the gamma, respectively.

Figure 4.11: True momentum distribution for the muon (left) and the gamma (right)
tracks.

49

Figure 4.12: Lifetime distribution for kaons that decayed in flight. These events are
the reason for the tails in the true muon momentum distribution of Figure 4.11

Figure 4.13: Reconstructed momentum distribution for the muon (left) and the
gamma (right) tracks.

50

Figure 4.14: Reconstructed minus true momentum distribution for the muon (left)
and the gamma (right) tracks.

Figure 4.15: Reconstructed minus true over true momentum distribution for themuon (left) and the gamma (right) tracks.
4.5.5Time Difference

As seen in Section 1.5.1, the time difference (∆t ≡ tµ − tγ ) between the muon andthe gamma vertices follows the decay curve characterized by the kaon lifetime of
12.380 ± 0.020 ns according to PDG. Figure 4.16 shows the true and reconstructed
time difference between the muon and the gamma vertices. The µγ-fitter does notconstrain ∆t to be positive and because of resolution negative values can be seenin the figure. Section 5.4.6 presents further discussions about the reconstructed ∆tdistribution.

51

Figure 4.17 shows an exponential fit to the true and reconstructed curves respectively and their best fit value for the kaon lifetime. Figure 4.18 shows the difference
between reconstructed and true ∆t and the resolution.

Figure 4.16: True and reconstructed distributions of the time difference between themuon and the gamma vertices (∆t) (log scale on the right).

Figure 4.17: True (left) and reconstructed (right) distributions of the time difference
between the muon and the gamma vertices (∆t).

52

Figure 4.18: Reconstructed minus true (left) and resolution (right) of the time difference between the muon and the gamma vertices (∆t).
4.5.6

Vertex Direction

Since there is no preferred direction for the emission of the de-excitation gamma, theopening angle between the muon and gamma is expected to be isotropic. As can beseen in Figure 4.19, fiTQun has no problem in finding the muon direction correctly.
99.7% of the time, the cosine of the angle between the true and reconstructed muon
directions is larger than 0.9. On the other hand, finding the gamma direction correctly is much more difficult. Using the same metric as before, only 47.2% of theevents are reconstructed in the correct direction. As Figure 4.20 indicates, there isa preference for locating the gamma in the same direction as the muon.

53

Figure 4.19: Angle between true and reconstructed directions for the muon (left) andthe gamma (right) tracks. The muon direction is well reconstructed in more than
99% of the events (shown in log scale), while the pointing resolution for the gamma
is much worse.

Figure 4.20: Cosine of true (black) and reconstructed (red) opening angle between
the muon and the gamma directions. The curves are normalized by area.
Instead of being isotropic as the true curve, the reconstructed distribution showsa concentration of events in the low opening angle region. This can be understoodin terms of the true time difference between the particles as shown in Figure 4.21. In
Section 4.5.1 it was described that the first initial value for the time offset between
the two particles is 2.5 ns, which is the timing resolution of ID PMTs in SK. For
54

Figure 4.21: Cosine of reconstructed opening angle versus the true time difference
between the muon and the gamma tracks.
very low values of true ∆t, it is hard for fiTQun to find the hits corresponding to
the gamma since these are very close in time with the hits coming from the muon.
As Figure 4.4 exemplified, the hits from a low energy gamma are very scattered andif they are close in time with a bright muon ring, the likelihood function in fiTQun
peaks when both particles point in the same direction, even though the true opening
angle might be large.
For the events with true ∆t > 10 ns, when the two particles are well separated
in time, fiTQun has a better performance as can be seen in Figure 4.22.
In Chapter 5 it will be shown that the most abundant atmospheric neutrino
interactions are CCQE. This type of background favors the low ∆t region, making
it impossible for this search to have any significance at very low values of ∆t. Thus,
the fact that it is hard for fiTQun to find the gamma when the kaon decays very fast
is not a problem, given that the rejection of these events is unavoidable if a large
signal to background ratio is to be achieved.

55

Figure 4.22: Cosine of true (black) and reconstructed (red) opening angle between
the muon and the gamma directions. The blue curve is the fraction of reconstructed
events that satisfy true ∆t > 10 ns. The curves are normalized by area.

56

5

Proton Decay Search

As it was discussed in Sections 1.5.1 and 4.5.2, the experimental signature of the
SUSY-favored proton decay mode, p → ν̄K + in the prompt-γ channel is composedof a low-energy de-excitation gamma (typically 6.3 MeV) and a monochromatic muonof 235.5 MeV/c, product of the two-body decay of the kaon. Besides the momentaof both particles, it is also known that there is a time difference between the hitscorresponding to the gamma and the the hits corresponding to the muon. The
process of nuclear de-excitation is essentially instantaneous, while the muon only
appears when the kaon decays. The exact time difference varies from event to event,but it is on the order of the kaon lifetime of about 12 ns. The order of the particles
was defined such that ∆t ≡ tµ − tγ is strictly positive for true signal events.
In this section, we discuss the selection criteria (cuts) applied in the analysis toseparate true signal events from background.

5.1

Signal and Background

Section 4.5.2 showed the performance of the µγ-fitter in signal simulation (MC), in
particular in the sample defined in Section 4.5 as true signal. These were defined as
proton decay events that contain a de-excitation gamma and the leptonic decay of thekaon to νµ µ+ . The choice of the selection criteria shown in the present chapter wasdone in order to select these true signal events. Despite this effort, other proton decay
events that do not meet both conditions can also satisfy the entire selection criteria.
These events are typically proton decay events with a monochromatic muon, but node-excitation gamma present. Because of that, signal is defined as any p → ν̄K +simulated event that satisfy all selection criteria.
Atmospheric neutrino interactions in the water can mimic this kind of signature as
explained in Chapter 3. Neutrinos do not produce any Cherenkov light when entering
the tank, therefore they are not vetoed by the OD. Depending on the interaction
between the neutrino and the nucleus, the final state particles can be the same as insignal events. Any interaction induced by an atmospheric neutrino is considered to
be background and the selection criteria were chosen to reject such events.

5.2

Efficiency and Expected Number of Background

In Chapter 1 it was discussed how the expected number of background events plays arole in a search for rare events. Ideally, a background free search could be developed
and the sensitivity for discovery would grow with the exposure of the detector. This57

is not possible for the prompt-γ analysis since there is an irreducible backgroundcoming from NC interactions that produce a charged kaon as defined in Chapter 3.
However, even though a background-free search is not possible, it is still desirable
to reduce the amount of background as much as possible while keeping the signal
efficiency high. Here we define the different metrics used to compare this work withprevious results of the SK collaboration and with other experiments.
In the previous search for p → ν̄K + in SK [38], the choice of the selection criteria was done to optimize the one event discovery sensitivity13 . In other words, theselection requirements minimize the number of expected background events while
maintaining the signal efficiency as high as possible, in order to maximize the discovery potential with a single event. In the present work, the optimization strategy
adopted was to obtain a similar expected number of background events as the previous analysis while increasing the signal efficiency, if possible.
Once the selection criteria was chosen based on the method described above, 90%of the remaining background events were related to kaon production in atmospheric
neutrino interactions, so called irreducible background for the present proton decay
search. The remaining 10% of background events were related to charged pions
produced in the primary neutrino interaction or secondary interactions of nucleons
propagating in the water. The particle signature of these events (PID and kinematics)
are different than signal events, but in some cases, such background events cansatisfy all the selection requirements. Given that the fraction of reducible background
events was negligible compared to events with kaons in the final state, a dedicatedsample to study kaon production in atmospheric neutrino interactions was generated(discussion in the next section). After this study was done, the final backgroundexpectation was larger (almost double) compared to the previous analysis. This canbe attributed to the analyzed background sample size used in the previous search,that was not large enough to account for statistical fluctuations of rare events andunderestimation of the production of kaons in neutrino interactions. Section 5.5
summarizes the background expectation as well as the composition of the background
divided into the different simulated neutrino interactions for reducible and irreducible
background.
The first key quantity used to compared both analyses is the signal efficiency,
defined in Equation 5.1.
=

number of signal events passing the selection,
total number of signal events inside true FV(5.1)

where true FV is defined as the events whose true vertex position are inside the13

The optimization did not include the effect of systematic uncertainties in the sensitivity.

58

fiducial volume. Another useful metric for signal is the efficiency using true µγ
events only, defined as µγ .
For the background, it is useful to have the expected number of events in two
different exposures. One is the exposure used in the data analysis and the otheris an exposure of 1 Mton·year. The first case is simply the amount of background
events that are expected in the analyzed data. The second number can be directlycompared to other experiments with different volumes and/or livetimes. In this work,
the analyzed data was taken during the SK-IV period, with a livetime of 2867.19 days.
Since the fiducial volume of SK has 22.5 ktons of water, the corresponding exposure
was 0.177 Mton·years.

5.3

Samples

Chapter 3 described the dedicated proton decay simulation of signal events. For the
background sample, NEUT was used as the event generator of atmospheric neutrino
interactions. Traditionally, analyses in SK use a background sample corresponding
to an exposure of 500 years for each SK period. Such simulated exposure is morethan sixty times larger than the SK-IV exposure used in this analysis. Despite thelarge number of simulated events, the actual background sample that can be misidentified as signal is very small. Using the traditional 500 year exposure sample,
only 10 background events were classified as signal, 9 out of those being NC kaons.
A larger sample of simulated background was necessary to reliably estimate theexpected number of background events in this search. For this reason, a sample
corresponding to an exposure of 5000 years was created.
In order to estimate the number of background events coming from kaons, adifferent sample was used, as described in Chapter 3. This sample used GENIE as
the neutrino event generator and corresponds to an exposure of 12,550 years. The
detector simulation was done with SKDETSIM and the events were reconstructedwith fiTQun, following the same simulation procedure as the 5000 year sample.
Reconstruction of proton decay events in SK requires two basic requirements usedin all distributions throughout this analysis. The first requires the event to be fully
contained (FC) and its vertex must be reconstructed inside the fiducial volume (FV).
Events that satisfy these conditions are part of the FCFV sample (see Chapter 2).
The other basic selection requires the events to have a visible energy14 of at least 30
MeV and at most 1330 MeV, these events are known as the subGeV sample.
14This is defined as the energy above Cherenkov threshold. It is calculated assuming the electron
hypothesis since the Cherenkov threshold for electrons is very low and essentially all their energy
is visible.

59

A sample of 578042 proton decay events were generated within 100 cm of thedetector walls, 142218 of those are true µγ-events inside the true FCFV region. From
here on the signal and background samples used in this analyses will be referred toas:• Proton decay signal: PDK• True µγ-events: µγ• 5000 year background: NEUT-5000• Kaon background: GENIE• SK-IV data (livetime: 2867.19 days): DATA
The estimation of signal efficiency was derived from the PDK MC sample andthe efficiency for observing the true µγ-events used the µγ-sample. The estimationof all background events that contain a kaon in the final state was obtained fromthe GENIE sample, while the expected number of background events for all otherneutrino interactions was derived from the 5000 year NEUT sample.

5.4

Selection Criteria

The method to select signal events while removing background is presented in thissection. Each subsection defines the variable used in the discrimination of signal
from background and explains the selection criterion by comparing signal and background MC. Signal and background samples are PDK and NEUT-5000 respectively,
as defined in the previous section.
In each subsection, two types of distributions are presented. The first uses a cumulative set of cuts, which means that each distribution is selected using the criteria
defined in the precedented subsections. The other is the so called n-1 distributions,where the entire selection criteria is applied with the exception of the selection thatis being presented. The first section on the number of rings already uses the basic
criteria defined previously as the FCFV, subGeV sample.
As it will be explained in Subsection 5.4.3 the analysis is split between the socalled two and three sub-event samples.
The selection for the two sub-event sample is:
1. Fully contained inside the fiducial volume (FCFV) sub-GeV event
2. Number of rings = 160

3. Number of sub-events = 2
4. PID is µ-like
5. 215 < pµ < 255 MeV/c
6. ∆X < 100 cm
7. 7 < ∆t < 80 ns
8. 2 < pγ < 10 MeV/c
The selection for the three sub-event sample is:
1. Fully contained inside the fiducial volume (FCFV) sub-GeV event
2. Number of rings = 1
3. Number of sub-events = 3
4. PID of the first sub-event is γ-like
5. PID of the second sub-event is µ-like
6. 215 < pµ < 255 MeV/c
7. ∆X < 100 cm
8. 7 < ∆t < 100 ns
9. pγ < 20 MeV/c
5.4.1

Ring Counting

Figure 5.1 shows the distribution of the number of rings identified in the events. It
can be seen that the PDK sample contains events with more than one ring, but thesetypically correspond to kaons that decay hadronically (K + → π + π 0 ). True µγ events
are reconstructed as one ring events the majority of the time, since the γ signal is
too weak to create another ring (Section 4.5).
Figure 5.2 presents the distribution of the number of rings after applying all
selection criteria except the one-ring selection. One can see that the remaining PDK
events are essentially all coming from the true µγ sample. The 2-ring bin has a large
contamination of background, therefore only the single ring sample is chosen.
61Figure 5.1: Reconstructed number of rings distribution (log scale on the right). The
PDK sample (black) contains events with more than one ring due to hadronic decayof the kaon. True µγ events (blue) are predominantly single ring. Background events(red) are more abundant than signal in the region with more than one ring. All curves
are area normalized.

Figure 5.2: Reconstructed number of rings distribution after applying the entire
selection except the number of rings (log scale on the right). The remaining PDK(black) sample is almost entirely true µγ events (blue) with a single ring identified.
Background is shown in red. All curves are area normalized.
5.4.2

Particle Identification

Selecting µ-like while rejecting e-like events is done using the likelihood of each hypotheses as explained in Section 4.1. Figure 5.3 shows the logarithm of the likelihood
ratio between µ and e hypotheses, known as the PID cut. At this stage, the back62

ground is dominated by CCQE interactions which result in a charged lepton in the
final state. For electron neutrinos the final state lepton is an electron and for muon
neutrinos it is a muon.
To
all the background coming from electron rings the applied selection removeLµ
is ln Le > 30. As in the number of rings case, the PDK signal has e-like rings
originated by the photons of hadronic decay of the kaon. But a very small fractionof true µγ events is mis-identified as electrons (1.3%).

Figure 5.3: Logarithm of the likelihood ratio between the µ and e hypotheses (logscale on the right). Proton decay signal (black) contains e-like events due to thehadronic decay of the kaon. True µγ events (blue) are predominantly µ-like. All
curves are area normalized.

63

Figure 5.4: Logarithm of the likelihood ratio between the µ and e hypotheses (logscale on the right) after applying the entire selection but PID. PDK (black) events
are almost entirely true µγ (blue) with few events mis-identified as e-like. Remaining
background (red) has a large fraction in the showering region, which shows thenecessity of the PID cut. All curves are area normalized.
5.4.3

Number of Sub-Events

In Section 4.2, the algorithm for finding sub-events in fiTQun was defined. In mostanalyses, the number of sub-events (nse) is a proxy for the number of Michel electrons
in the main event. For example, in a CCQE event with an electron in the final state,
there is no Michel electron and the number of sub-events is typically 1 (the main
event itself). For a muon that decays inside the tank, one Michel electron will bepresent in the event and the number of sub-events is two. Thus, in most cases thefollowing relation is valid: NMichel−e = Nsub-events − 1. Analyses that require a singleµ-like ring would normally require nse= 2.
In this work, an interesting counter-example for this proxy occurs. In a small
number of cases, fiTQun is able to identify the de-excitation gamma as a sub-event(Figure 5.5). This happens when the hits from the muon and the gamma are sufficiently separated in time (usually more than 20 ns). In those cases, true signal
events have three sub-events, the first being the prompt gamma, the second themonochromatic muon and finally, the Michel electron.
When fiTQun is able to find these three clusters of hits separately, the dedicatedµγ-fitter is not necessary and the standard single ring fitter is used. Although this
3-nse sample is small, it is extremely pure and the only neutrino interaction capableof producing a similar signature is the usual irreducible background of NC kaon
production.

64

Therefore, a selection is made for the 2-nse and 3-nse samples and only the
1-nse sample is removed. Figure 5.6 shows the number of sub-events distribution
after applying the entire selection except the requirement on the number of subevents. It can be seen that explicitly removing the 1-nse sample is necessary, assome background events still remain in that bin.
From this point, the analysis is separated in two samples, for the 2-nse and 3nse cases. The 2-nse sample uses the dedicated µγ-fitter to reconstruct the events.
While, the 3-nse sample uses the standard single ring fitter from fiTQun, since thesub-event algorithm was capable of identifying the gamma.

Figure 5.5: Number of sub-events (nse) distribution for PDK (black), true µγ (blue)and background (red) samples (log scale on the right). Some true µγ events have
3-nse when fiTQun is able to identify the γ as a sub-event before the µ. All curves
are area normalized.

65

Figure 5.6: Number of sub-events (nse) distribution for PDK (black), true µγ (blue)and background (red) samples (log scale on the right) after applying the entire selection but the number of sub-events requirement. The 1-nse bin shows the need to
select only the 2 and 3-nse samples. All curves are area normalized.
5.4.4

Muon Momentum

Section 4.5.4 showed the resolution for the reconstructed muon momentum for truesignal events. As previously discussed this muon is expected to be monochromaticwith pµ = 235.5 MeV/c. Figure 5.7 shows the distribution of reconstructed muon
momentum for signal and background. The momentum spectrum for CCQE events
is essentially flat in the region close to signal as can be seen in Figure 5.8.
The selection criteria for pµ is based on the muon momentum resolution studiedin Section 4.5.4. The center value and width of the gaussian fit was pµ = 236.6 ±5.8 MeV/c. A 3-σ interval was chosen as the selection criteria for muon momentum,
such that events satisfying 215 < pµ < 255 MeV/c are selected.

66

Figure 5.7: Reconstructed muon momentum distribution for PDK (black), true µγ(blue) and background (red) in the 2-nse sample (log scale on the right). PDK
events are reconstructed with a 5.8 MeV/c resolution centered around 236.6 MeV/c
as shown in Section 4.5.4. All curves are area normalized.

Figure 5.8: Same as Figure 5.7 with the entire momentum spectrum of sub-GeV
muons created in atmospheric neutrino interactions.

67

Figure 5.9: Reconstructed muon momentum distribution for PDK (black), true µγ(blue) and background (red) samples after applying the entire selection except themuon momentum requirement (log scale on the right). All curves are area normalized.
The same selection criteria is applied for the 3-nse sample, but this time thevariable used as the reconstructed muon momentum comes from the second subevent of the single ring µ-hypothesis of fiTQun. Figure 5.10 shows the distributionof pµ in the 3-nse sample. A small shift in the center value can be observed comparedto the 2-nse case, but the chosen interval still covers more than 99% of the signal
region, thus it was kept the same.

Figure 5.10: Reconstructed muon momentum distribution for PDK (black), true µγ(blue) and background (red) in the 3-nse sample. All curves are area normalized.

685.4.5Distance Between the Muon Decay and the Michel Electron Vertex
Positions

Using the reconstructed muon position, direction and momentum, it is possible tocalculate the expected position where it decays. The distance between the expectedmuon decay position and the Michel electron vertex position is defined as ∆X. The
true value of this distance is zero, since the decay position of the muon is the same
position where the Michel electron is created. Resolution effects cause this distanceto be non-zero. At first, this variable seems to be a weak discriminator between
signal and background as shown in Figure 5.11.

Figure 5.11: Reconstructed distance between the muon decay and the Michel electron
vertices (∆X) for PDK (black), true µγ (blue) and background (red) samples (logscale on the right). All curves are area normalized.

69

Figure 5.12: Reconstructed distance between the muon decay and the Michel electron
vertices (∆X) for PDK (black), true µγ (blue) and background (red) samples after
applying the entire selection except the ∆X requirement (log scale on the right) .
All curves are area normalized.
On the other hand, Figure 5.12 shows the ∆X distribution after applying theentire selection. It can be seen that more than 80% of the background is in theregion of ∆X greater than 100 cm. Figure 5.13 starts to explain why this is the case.
The green curve shows the distribution of Neutral Current atmospheric neutrino
elastic interactions that produced protons above the Cherenkov threshold.

Figure 5.13: Reconstructed distance between the muon decay and the Michel electron
vertices (∆X) for PDK (black), true µγ (blue), background (red) and NC elastic(green) samples after applying the entire selection but the ∆X requirement (logscale on the right) . All curves are area normalized.
In these cases, a proton is emitted above Cherenkov threshold by the nucleus
70after the NC interaction. Since protons are very heavy (compared to electrons)
they also produce a sharp non-showering ring and they can be mis-identified as amuon. Figure 5.14 shows the typical proton momentum for NC elastic events thatare mis-identified as single-ring µ-like events with reconstructed momentum in thesignal region. It also shows the angle between the true proton direction and thereconstructed ring direction, it can be seen that fiTQun correctly finds the direction
in most cases.
There are two main differences between events with a proton and a muon ring.
The first is the fact that, just like pions, protons interact hadronically, thus some
proton tracks might have a kink where the proton interacted, shifting the directionof the ring. For protons with momentum approximately 1500 MeV/c the chance ofa hadronic interaction is about 40% [39]. The second main difference is the opening
angle of the Cherenkov cone.
1, since protons
As seen in Chapter 1, the Cherenkov angle is given by cos θC = βnare much heavier, they are typically less relativistic then muons. A muon with
momentum of 235.5 MeV/c has βµ = 0.91, while a 1500 MeV/c proton has βp = 0.85.
The typical Cherenkov angle for signal muons is θC,µ = 34.5◦ , while protons thatmimic this signal have Cherenkov angles that range from 25◦ to 35◦ .
If a proton ring is reconstructed using the µ-hypothesis, its vertex has to be
displaced such that the opening angle is larger. Figure 5.15 is a cartoon showing
the vertex displacement that has to happen so that the opening angle in the µhypothesis becomes larger, i.e, there must be a displacement in the positive directionof propagation of the proton. The displacement is given by Equation 5.2:tan θp,(5.2)
d = dp · 1 −tan θµwhere d is the displacement, dp is the distance of the true proton vertex to thenearest wall and the angles are the Cherenkov angles. This example only works for
tracks that point directly to the wall, but it can be generalized to any case. For
dp = 800 cm (proton between the middle of the tank and the wall), the displacement
is about 200 cm. Figure 5.16 shows the vertex displacement for protons that have
∆X greater than 100 cm. It can be seen that essentially all the mis-reconstruction
happens in the longitudinal direction.
Besides NC elastic interactions, proton rings can occur in CCQE interactions
when the neutrino is energetic enough to put the proton above Cherenkov threshold.
In such cases, the muon may be below Cherenkov threshold and the event has thesame signature as the NC case. All NC proton rings in the 5000 year sample can beremoved with the selection ∆X < 100 cm.
71

Figure 5.14: True proton momentum (left) and cosine of angle (right) between thetrue proton direction and reconstructed ring direction for NC elastic interactions
that mimic µ-like events.

Figure 5.15: Vertex displacement that happens when a proton ring is misreconstructed as a muon.

72

Figure 5.16: Displacement: distance between the true proton vertex and the reconstructed vertex in the µ-hypothesis (left). Displacement projected in the longitudinal
direction of the proton (right). This effect happens because of wrong PID hypothesis
as explained in the text and represented in Figure 5.15.
Similarly to the momentum requirement, the 3-nse sample behaves in the same
way as the 2-nse sample. The position and momentum used to calculate ∆X arefrom the single ring µ-fitter instead of the dedicated µγ-fitter and the Michel electron
position used the third sub-event vertex. Figure 5.17 shows the ∆X distribution forthe 3-nse sample, the selection is the same as before, ∆X < 100 cm.

Figure 5.17: Reconstructed distance between the µ decay and the Michel electron
vertices (∆X) for PDK (black), true µγ (blue) and background (red) samples (logscale on the right). All curves are area normalized.

73

5.4.6Time Difference Between the Muon and the Gamma Vertices

As explained in Chapter 1, there is a time difference between the emission of thede-excitation gamma by the remaining nucleus and the µ-ring coming after the kaon
decays leptonically. The expected distribution of the time difference, ∆t defined astµ − tγ , is given by an exponential distribution with a time constant equal to thelifetime of the Kaon at rest of about 12 ns. This is the key feature of the signal thatallows for the major rejection of background events coming from CCQE interactions.
Figure 5.18 shows the ∆t distribution for signal and background.
Two features are worth noting. The first is that some small fraction of true signal
events is reconstructed with negative ∆t. Even though this is unphysical for true µγ
events, the algorithm does not constrain ∆t to be positive to not bias the likelihood
function. These events are only removed at analysis stage. Another feature is thatbackground events are heavily concentrated around ∆t = 0.

Figure 5.18: Time difference between the muon and the gamma vertices for PDK(black), true µγ (blue) and background (red) in the 2-nse sample (log scale on theright). All curves are area normalized.
This can be understood by looking at the fraction of background events corresponding to CCQE interactions as in Figure 5.19. Even if a neutrino interaction produces a muon with the right signal momentum and the nucleus emits a de-excitation
gamma, this type of background event can still be removed since the emission of bothparticles happens at the same time. Unfortunately, due to timing resolution of the
PMTs and reconstruction failures, sometimes a CCQE event with no time difference
can still be reconstructed with a finite ∆t between the two particles. Thus, only
requiring a positive ∆t is not enough to remove all the CCQE background.

74

Figure 5.19: Time difference between the muon and the gamma vertices for PDK(black), true µγ (blue), background (red) and CCQE-only (green) samples in the2-nse case (log scale on the right). All curves are area normalized.

Figure 5.20: Time difference between the muon and the gamma vertices after applying all other selection criteria for PDK (black), true µγ (blue), background (red) and
CCQE-only (green) samples in the 2-nse case (log scale on the right). All curves arearea normalized.
As seen in Figure 4.18, the reconstructed ∆t resolution is very close to the PMT
resolution of 2.2 ns. By requiring ∆t to be greater than 7 ns, the events are alreadymore than 3σ away from the ∆t resolution and signal can be separated from CCQE
background. Some CCQE events still remain after this requirement as it is shownin the log-scale plot of Figure 5.19. These events are removed after the gamma
momentum selection is applied (Figure 5.20) as will be discussed in the next section.

75

The final selection is given by 7 < ∆t < 80 ns, where the upper bound was arbitrarily
chosen in the tail of the signal distribution15 .

Figure 5.21: Time difference between the muon and the gamma vertices for PDK(black), true µγ (blue) and background (red) in the 3-nse sample (log scale on theright). All curves are area normalized.
Figure 5.21 shows the ∆t distribution for the 3-nse sample. Signal events start at
late times (about 20 ns) as expected from the discussion of Section 5.4.3. FiTQun
is only able to find the gamma sub-event if the two particles are well separated intime. The same lower bound can be applied for both 2 and 3-nse samples, but sincethe 3-nse sample is enriched for large values of ∆t, the upper bound is increased to
100 ns. For this sample, selection requirement is given by 7 < ∆t < 100 ns.
5.4.7

Gamma Momentum

Figure 5.22 shows the reconstructed gamma momentum for the remaining background events that passed all the selection criteria so far. More than 90% of thisremaining background is concentrated in the pγ < 2 MeV/c region. Figure 5.23 includes the fraction of events that do not have a gamma present in the interaction(about 97%).
FiTQun is a maximum likelihood fitter, which means that given a particularhypothesis it will optimize the parameters for the given particle content. In the caseof this analysis, the hypothesis is that there are two particles in the event: a photon
and a muon. In case fiTQun is not able to find the gamma (when there is none, forexample), it returns a very low value for its momentum. This is best understood by
15

No significant change happens when this selection is varied by ±10 ns.

76

looking at the correlation between pγ and opening angle. For clarity, let us postpone
this discussion to the next section and finalize the selection criteria here.
A lower bound of 2 MeV/c removes more than 90% of the remaining background
while keeping most of the true signal. The few PDK events in the region with lessthan 2 MeV/c are monochromatic muons with no nuclear de-excitation gamma. A
symmetric interval around 6 MeV/c was chosen as the signal region for the 2-nse
sample: 2 < pγ < 10 MeV/c.

Figure 5.22: Reconstructed gamma momentum for PDK (black), true µγ (blue) andbackground (red) samples in the 2-nse case (log scale on the right). All curves arearea normalized.

Figure 5.23: Reconstructed gamma momentum for PDK (black), true µγ (blue),
total background (red) and background events with no de-excitation gamma present(green) in the 2-nse sample (log scale on the right). All curves are area normalized.

77

Figure 5.24: Reconstructed gamma momentum for PDK (black), true µγ (blue) andbackground (red) in the 3-nse sample (log scale on the right). All curves are area
normalized.
Figure 5.24 shows the gamma momentum distribution for the 3-nse sample. The
remaining background events in this sample that are close to the signal region aretrue kaon events produced in NC interactions. Thus, an upper bound on gamma
momentum removes all reducible background for this search. The remaining NC
kaon events are completely irreducible as they can not be distinguished from kaons
that were originated in proton decay events (discussion in Chapter 3). The final
selection cut for this sample is pγ < 20 MeV/c.
Section 5.5 summarizes the selection criteria and quantifies the signal efficiency
loss and background reduction for each cut. The final efficiency and expected numberof background events can then be used to calculate the sensitivity for the lifetime
limit.
5.4.8

Reconstruction of Events with no De-excitation Gamma

The difference between reconstruction of events with and without a de-excitation
gamma present is presented in this section. Simulated proton decay events where
the kaon decayed leptonically to a mono-energetic muon are used in this study.
Let α be the reconstructed opening angle between the muon and the gamma
directions. Figures 4.21 and 4.22 show the behavior of α as a function of the true
time difference between the muon and the gamma. It was discussed in Section 4.5.6,
that for events with a very low ∆t it was less likely for fiTQun to correctly reconstruct
the gamma direction.
Figure 5.25 shows the cosine of α as a function of reconstructed gamma momentum for events where a de-excitation gamma was present. Figure 5.26 shows the78same distribution for the case where no de-excitation gamma is present in the event.
It can be seen that when no gamma is present in the event, the reconstructed gamma
hypothesis tends to be in the same direction as the muon and with very low momentum. The requirements of lower bounds for the time separation and the gamma
momentum remove these events where either the gamma reconstruction failed or
there was no de-excitation gamma present in the event.

Figure 5.25: Cosine of the reconstructed opening angle between the muon and the
gamma directions as a function of the reconstructed gamma momentum for eventswith a de-excitation gamma present (log scale).

79

Figure 5.26: Cosine of the reconstructed opening angle between the muon and the
gamma directions as a function of the reconstructed gamma momentum for events
without a de-excitation gamma present (log scale).

5.5

Summary of Selection Efficiencies and Expected Number of Background Events

Tables 5.1 and 5.2 summarize the signal efficiency and expected number of background events for each selection criterion in the two and three sub-event samples,
respectively. For each selection, the signal efficiency and the efficiency for true µγ
events are shown. The expected number of background events is shown for the SK-IV
livetime used in data analysis of 2867.19 days and for an exposure corresponding to
1 Mton·year.
The combined signal efficiency was estimated to be 14.2% and the expected number of background events was 0.30. In both cases the statistical uncertainty is muchsmaller than the systematic uncertainty and the final uncertainties of signal efficiency
and expected number of background events are discussed in the next chapter.

80

Table 5.1: Signal efficiency, true µγ efficiency, expected number of background events
in SK-IV and Mton·year rate for the 2 sub-event sample.
Selection
FCFV+subGeV
Single-Ringµ-like2-nse
215 < Pµ < 255 MeV/c∆X < 100 cm
7 < ∆t < 80 ns
2 < Pγ < 10 MeV/c

 (%)
92.3
63.7
59.5
57.9
54.3
52.5
14.8
13.5

µγ (%)
94.5
93.2
91.9
88.0
82.4
79.7
45.2
43.3

#Bckg (SK-IV) Mton·year Rate
17560.899425.1
12585.2
71254.4
5896.233383.1
5029.028472.7
432.4
2447.9125
400.0
2263.6
3.2
18.1
0.29
1.62

Table 5.2: Signal efficiency, true µγ efficiency, expected number of background events
in SK-IV and Mton·year rate for the 3 sub-event sample.
Selection
FCFV+subGeV
Single-Ring
3-nse
1st se γ-like2nd se µ-like
215 < Pµ < 255 MeV/c∆X < 100 cm
7 < ∆t < 100 ns
Pγ < 20 MeV/c

 (%)
92.3
63.7
1.4
1.3
0.80
0.80
0.78
0.77
0.76

µγ (%)
94.5
93.2
2.3
2.66
2.63
2.61
2.54
2.51
2.50

#Bckg (SK-IV) Mton·year Rate
17560.899425.1
12585.2
71254.4
308.5
1746.4
32.8
185.9
0.38
2.16
0.085
0.49
0.028
0.15
0.020
0.11
0.015
0.085

Table 5.3 summarizes all neutrino interaction modes that produced background
events satisfying all the selection criteria and did not have a kaon in the final state(reducible background). These events constitute 15.2% of the final number of expected background events, the remaining 84.8% are events with a kaon in the final
state (irreducible background).

81

Table 5.3: Composition of background events separated by neutrino interaction
modes.
Mode
CCQE
CC 1-πNC 1-πNC elastic
Total

82

(%)
5.4
0.27
5.3
4.2
15.2

6

Systematic Uncertainties

Chapter 5 described the selection criteria used in this analysis to discriminate between signal and background events. The obtained values for signal efficiency andexpected number of background events were derived using simulated MC events. A
number of assumptions were made during the entire process of event simulation and
it is necessary to estimate the uncertainty associated with the obtained parameters.
The uncertainties can be classified in two, both will be referred to as systematic
uncertainties.
The first type of systematic uncertainties is related to the models used in thesimulation, such as the gamma emission probabilities in de-excitation processes inTable 3.1 or the cross-section model used in neutrino interactions in NEUT. Thesemodels are approximate descriptions of natural phenomena and it is necessary toestimate their validity. These uncertainties are referred to as theoretical systematic
uncertainties.
Similarly, the model used to simulate the propagation of particles in the SK
detector is also an approximate description of the real propagation of particles in a
medium such as water. Despite the efforts to tune these models to best agree with experimental data, there are still limitations and discrepancies and an estimation of thedetector model validity is also needed. Even though the uncertainties associated withthe detector simulation come from a model just as the theoretical uncertainties, they
were treated differently. The reason is that the effect of theoretical uncertainties canbe easily propagated to the final numbers by applying weights in an event-by-event
basis, without the need to run a new MC simulation. For the detector parameters
this can not be done since changing even a single parameter of the detector, for example, the water attenuation length, would require a completely new MC simulation,
which is impractical. In order to address this difficulty, the estimation was done at
the level of analysis quantities such as the reconstructed momentum, position, etc.
The idea is that an uncertainty in the detector model propagates to an uncertainty inreconstructed quantities and these are easier to treat. The uncertainties associatedwith selection criteria are referred to as selection systematic uncertainties.
The method used to estimate the systematic uncertainties associated with thiswork are presented in this chapter. Uncertainties associated with theoretical models
are presented first for signal and background parameters, respectively. Selection
uncertainties are either common between the two samples or estimated using a similar
method and are treated together in the following sections.

83

6.1

Emission Probabilities of De-excitation Gamma Rays

The process of simulation of proton decay events was described in Section 3.1. It was
discussed how nuclear effects such as Fermi momentum or binding energies do notaffect the kinematics of the observed signal in the SK detector. The reason was that
nuclear phenomena affect directly the final state particles of proton decay events,
but for ν̄K + , the visible particles come from the decay of the kaon, unaffected by theremaining nucleus. The only modeled effect that has a direct impact in this work is
the emission probabilities for de-excitation gammas.
The emission probabilities were first derived in Reference [21] and are summarizedin Table 3.1. In the last published proton decay analysis by SK [38], the estimationof the uncertainties related to this calculation was done using experimental datafrom Reference [40]. The experiment measured de-excitation gamma rays from thes-hole state in 15 N produced via the 16 O (p, 2p)15 N reaction. The difference between
experimental results and the theoretical prediction was used as the systematic uncertainty for the emission probabilities. The 6 MeV line was estimated to have 15%
uncertainty in its emission probability while all other lines were estimated to have30%.
The final systematic uncertainty in signal efficiency due to the uncertainty in theemission probabilities of the gamma lines was estimated to be 23.2% and it is thedominating systematic uncertainty for signal efficiency.

6.2

Atmospheric Neutrino Flux and Cross-Section Models

The systematic uncertainties associated with the atmospheric neutrino flux and crosssection models are evaluated using an event-by-event weighting method for each
uncertainty based on the atmospheric neutrino analyses by SK [41].
For the neutrino flux model, the following sources of uncertainty are consideredin this analysis:• Energy dependent normalization• Neutrino flavor ratio• ν̄/ν ratio• Up/down asymmetry• Horizontal/vertical ratio• K/π production ratio84• Neutrino flight length
The largest contribution for the final systematic uncertainty of the neutrino flux
comes from the energy normalization, which is estimated to be 10.1%.
The following uncertainties related to the neutrino cross-section model were considered:• Quasi-elastic scattering for bound nucleons (total cross section)• Quasi-elastic scattering for bound nucleons (ν̄/ν)• Quasi-elastic scattering for bound nucleons (flavor ratio)• Single (multi)-meson production (total cross section)• Single (multi)-meson production (model dependence)• Coherent pion production• NC/CC ratio
The largest uncertainties related to the neutrino cross-section model are the total
cross section for CCQE events and single pion production with 9.3% and 6.3%,
respectively. As discussed in Section 3.2.4, the largest background for the proton
decay search comes from charged kaons produced in NC interactions. The uncertainty
in the cross section of this interaction was considered separately and it is describedin Section 6.3.
The total systematic uncertainty in the expected number of background events
related to the neutrino flux and cross-section models is estimated to be 21.9%. Eachcontribution was taken to be independent from the others and the uncertainties wereadded in quadrature.

6.3

Kaon Production in NC Neutrino Interactions

As reported in Section 3.2.4, kaons produced in neutral current atmospheric neutrino interactions constitute the largest background for the proton decay search inthis work. The cross section for this interaction was recently measured by the MINERvA collaboration [42]. Figure 6.1 shows the differential cross section of charged
kaon production as a function of the kaon kinetic energy. The figure also showsthe expected cross section according to NEUT and GENIE event generators. As
described in Section 3.2.4, the cross-section model for kaon production in NEUT85

is known to be incomplete as it does not include contributions from the associatedproduction with sigma baryons. GENIE was chosen as the event generator for atmospheric neutrino interactions that produce kaons in the final state. Figure 6.2 shows
the kinetic energy of the kaons produced in NC neutrino interactions that satisfied all
selection criteria. It can be seen that most kaons have kinetic energies less than 300
MeV and the distribution is essentially flat in this region. The Cherenkov thresholdfor the kinetic energy of kaons in water is about 250 MeV. The events with morethan 300 MeV kaon kinetic energy produce extra light in the event and do not satisfy
the selection criteria.
The expected number of background events is then the addition of the number of
events with kaons in the final state predicted by GENIE and the events with no kaonspredicted by NEUT. The relative uncertainty in the GENIE prediction was estimatedto be 32.9%. The uncertainty in the measurement itself used the total error bar in thefirst bin of Figure 6.1, which corresponds to 32.7%. This bin was chosen since it has
the largest uncertainty out of the relevant bins for this work. The difference between
the central value of the observed data and the GENIE prediction was also considered
and it corresponds to a relative uncertainty of 3.6%. Both uncertainties were addedin quadrature to estimate the final uncertainty used for the kaon production crosssection.

Figure 6.1: Differential cross section of charged kaon production in neutral current
neutrino interactions as a function of kaon kinetic energy [42].

86Figure 6.2: Kinetic energy of final state kaons produced in NC atmospheric neutrino
interactions that satisfy all selection criteria of the proton decay analysis.

6.4

Selection Criteria

As discussed in the beginning of this chapter, the systematic uncertainties associated with detector simulation were estimated through their effect on reconstructed
quantities. The uncertainties related to the muon reconstruction are described here.
The selection of the gamma is discussed separately in Section 6.5, as a new method
was implemented.
6.4.1

Fiducial Volume

The fiducial volume systematic was estimated using the reconstructed vertex positionof atmospheric neutrino data and simulated events in the region outside the fiducial
volume. A sample of cosmic muons was used to estimate the contamination of cosmic
muons in the atmospheric neutrino sample. It was found that for dwall greater than
100 cm, the contamination from cosmic muons is less than 0.1%, where dwall isdefined as the distance between the vertex position to the nearest wall of the tank.
The atmospheric neutrino sample was then used to estimate the fiducial volume
uncertainty by quantifying the difference between the number of events in the data
and the simulated samples.
The standard fiducial volume in the SK tank is defined to be two meters awayfrom the detector wall. The expanded region defined in this study was chosen to
be from one to two meters away from the wall to avoid contamination from cosmic
87

muons as discussed above. The number of simulated atmospheric neutrino eventswas normalized to the number of data events inside the standard FV region and thesame weight was used in the expanded region. Figure 6.3 shows the reconstructedvertex position for data and simulated events inside the FV and in the expanded
regions. The difference between data and simulation in the expanded region was
estimated to be 1.2%. The same uncertainty was used for the signal efficiency andexpected number of background events.

Figure 6.3: Distribution of dwall for atmospheric neutrino MC simulated events,
where dwall is the distance of the reconstructed vertex to the nearest tank wall.

6.4.2

Number of Rings

The number of rings in fiTQun reconstruction is determined by a complicated algorithm. Essentially there is a relation between the number of rings and the logarithmof the likelihood ratio between the single-ring and multi-ring hypotheses. Figure 6.4
shows the logarithm of the ratio of likelihood hypotheses for atmospheric neutrino
data and simulation. A fit was performed in the simulated distribution to best agreewith the data and the parameters of the fit were used to estimate the systematic
uncertainty.
The observed number of rings for events in the phase space of this analysis is very
well determined by fiTQun and the uncertainty is very low. Single muon events with
momentum less than 500 MeV/c are well reconstruct by the single ring hypothesis
and the uncertainty was estimated to be 0.1% for signal efficiency and expected
number of background events.

88Figure 6.4: Logarithm of the likelihood ratio between the single-ring and multiring hypotheses for atmospheric neutrino events that satisfy µ-like and one Michelelectron criteria.
6.4.3

Particle ID

A similar procedure for the estimation of the PID uncertainty was adopted following
the number of rings case. Atmospheric neutrino simulation was fit to data and theparameters that optimize the χ2 distribution were chosen to estimate the uncertainty.
Figure 6.5 show the data and simulated distributions before and after performing thefit for the logarithm of the ratio between the e-like and µ-like hypotheses for single
ring events with a Michel electron.
Similarly to before, the events that satisfy the entire selection criteria are single
ring muons with a likelihood value well above the threshold of the analysis. No
candidate event for signal or background samples shift after applying the parametersoptimized in the fit and the systematic uncertainty for PID was estimated to be 0%
for signal efficiency and expected number of background events.
6.4.4

Number of Michel Electrons

The same sample with cosmic muons used in the number of rings uncertainty estimation was used to estimate the uncertainty on the number of Michel electrons. Cosmic
muons that satisfied the subGeV, single-ring µ-like criteria were selected. The requirement of one Michel electron was applied to the data and simulated samples, andan efficiency for this selection was calculated for both samples. The ratio between
89

Figure 6.5: Logarithm of the likelihood ratio between e-like and µ-like hypotheses
for single ring events with a Michel electron.
both efficiencies was used to calculate the systematic uncertainty, as in Equation 6.1
syst = 1 −

MC
DATA

,

(6.1)

where  is the selection efficiency of the Michel electron requirement. The systematic
uncertainty was estimated to be 1.5% for signal efficiency and expected number ofbackground events.
6.4.5

Distance Between the Muon Decay and the Michel Electron Vertex
Positions

Section 5.4.5 defined ∆X as the distance between the muon decay vertex and theMichel electron vertex positions. The estimation for the systematic uncertainty associated with this selection was done in the same way as the estimation for the numberof Michel electrons. The cosmic muon sample was used and muons satisfying thesubGeV, single µ-like ring with a Michel electron were selected. The same method
as in Equation 6.1 was applied and the systematic uncertainty associated with ∆X
was estimated to be 0.8% for signal efficiency and expected number of background
events.

90

6.4.6

Energy Scale

The systematic uncertainty on the energy scale is estimated by comparing dataand simulation in three control samples. The first control sample used is the same
cosmic muon sample as before. Michel electrons and neutral pion events are usedas the other control samples. The Michel electron momentum distribution can becompared to the well known Michel spectrum and the π 0 sample is used to reconstruct
the π 0 mass peak. The largest difference between data and simulation is used as theenergy scale systematic uncertainty and it was determined to be 2.2%. Propagating
this uncertainty into the reconstructed momentum of the muon and the gamma, asystematic uncertainty of 2.5% was estimated in the signal efficiency and 10.3% for
the expected number of background events.
The difference between the two uncertainties can be understood by analyzing
Figure 5.22. The proton decay signal events that passed the selection criteria up
to the gamma momentum selection are almost entirely true µγ signal events, sincemost of the fake gamma events were removed by the ∆t selection. The reconstructedgamma momentum distribution of signal events in Figure 5.22 does not change significantly when a 2.2% shift is applied to the gamma energy. On the other hand,the distribution for background events is different since there is a big contribution
coming from events with no gamma. These events can enter the signal region whenthe energy shift is applied causing a larger change in systematic uncertainty.

6.5

Hybrid Sample

In order to estimate the systematic uncertainty in the time difference selection, ahybrid sample was constructed since there were no control samples available to estimate this value. A hybrid sample for this analysis was defined as a muon selectedfrom the atmospheric neutrino data or MC simulation, whose hits were merged withhits from a simulated gamma. The idea of this study was to estimate the effectof detector simulation on the reconstruction of the time difference between the twoparticles, ∆t. The sample just described is called the primary hybrid sample, where
the brighter ring of the event is selected from data or MC simulation, while the faintring is always simulated. A secondary sample can also be defined by selecting the
gamma ring from data or MC simulation and merging the hits of these events withhits from a simulated muon. The difficulty with the latter is the lack of a control
sample where single gammas can be selected.
Since no isolated gamma can be found in control samples of the SK detector, nosecondary sample was used in this study. Previous analyses used this technique of
constructing such hybrid samples but this was the first time where the hits from each91

ring needed to be merged with a time difference. Figure 6.6 is a proof of concept thatan artificial ∆t between the hits could be introduced. In this example a non-physicaltime difference of 500 ns was applied between the hits of the muon and the gamma,
both from simulated samples.
Following the same procedure as before, the estimation of the systematic uncertainty associated with the reconstructed time difference between the gamma andthe muon was done by calculating the difference in the ∆t selection efficiency fordata and MC simulation samples. Using the same criteria as in Equation 6.1, theestimated uncertainty was 10.5% for signal efficiency and expected number of background events.

92

Figure 6.6: Example of merging hits for the construction of the hybrid sample. A
muon (top left) with momentum of 236 MeV/c was generated at the center of thetank pointing at (-1, 0, 0). A gamma (top right) with momentum 6 MeV/c wasgenerated at the same position pointing backwards at (1, 0, 0). The hits of both
events were merged together to build a hybrid event with a muon and a gammapresent (bottom left). An artificial time difference of 500 ns was applied between the
gamma and muon hits. The time distribution on the bottom shows a small clusterof hits in the 500 ns region corresponding the the gamma hits.

6.6

Systematic Uncertainties Summary

The final systematic uncertainty for signal efficiency was estimated to be 25.7%.
Table 6.1 summarizes the theoretical and selection systematic uncertainties for signal
efficiency. The systematic uncertainty in the expected number of background events93

was estimated to be 42.2% and Table 6.2 summarizes each component individually.
For comparison, the systematic uncertainties in the last published analysis [38] were
22% for signal efficiency and 25% for the expected number of background events.
The signal efficiency uncertainties are similar as expected since no significant change
was introduced between the two analyses besides the change in event reconstruction.
The number of background events was treated differently, not only the estimation
was done with a different event generator, but also the kaon production cross-section
uncertainty included in the present analysis is new. The previous analysis did nothave an estimation for this uncertainty since the MINERvA measurement was notavailable at the time. Excluding the kaon production uncertainty, the systematic
uncertainty for the background expectation is 26.5%, comparable to the previous
number.
Table 6.1: Systematic uncertainty summary for signal efficiency.
Systematic Uncertainty (%)
Sourceγ emission probability
23.2
11.0
Selection efficiency
Total
25.7

Table 6.2: Systematic uncertainty summary for the expected number of background
events.
Systematic Uncertainty (%)
SourceAtmospheric ν flux and cross section
21.9
32.9
Kaon production cross section
14.9
Selection efficiency
Total
42.2

94

7

Data Results

Chapter 5 described the selection criteria used to discriminate signal from background
events (summary in Section 5.4). The selection method presented was applied forthe data and simulated background samples. Data was collected during 2867.19 daysof SK-IV as described in Section 5.3 and the expected number of background eventswas derived using the simulated atmospheric neutrino sample with an exposure of5000 years.
Two different sets of requirements are presented in this chapter. Initially, dataand simulated events are compared in the so called sideband regions, defined in thenext section. Finally, Section 7.2 presents the comparison between events from thedata and simulated samples using the selection criteria from Chapter 5.

7.1

Sidebands

This section presents the comparison between events from the data and the simulated
samples in two sideband regions, defined as a function of muon momentum. As
described in Chapter 5, the region of muon momentum in the signal sample is definedas 215 < pµ < 255 MeV. The regions pµ < 215 MeV and 255 < pµ < 355 MeV werechosen as sidebands for the analysis to perform comparisons between data and MC
simulation before studying the selected signal region. The lower momentum region
does not have a lower bound since it reaches the muon Cherenkov threshold before
reaching 115 MeV.
Figure 7.1 shows the distribution for the number of identified rings in the FCFV,subGeV sample for events satisfying pµ < 215 MeV or 255 < pµ < 355 MeV. Only
the one ring sample is selected and Figure 7.2 shows the PID distribution between
the muon and electron hypotheses for the selected ring.

95

Figure 7.1: Number of identified rings distribution for data (black) and simulated
atmospheric neutrino background (red) samples (log scale on the right) for events inthe sideband regions. Events with one identified ring are selected.

Figure 7.2: PID distribution for data (black) and simulated atmospheric neutrino
background (red) samples (log scale on the right) for events in the sideband regions.
Events with the logarithm of the likelihood ration greater than 30 (blue) are selected
as µ-like events.
The distribution of the number of sub-events is shown in Figure 7.3. The twoand three sub-event samples are selected and Figures 7.4 and 7.5 show the muon
momentum distributions for the relevant sub-event. Only one event is expected fromthe simulated three sub-event sample and none was observed in the sideband regions.

96

Figure 7.3: Number of sub-event distribution for data (black) and simulated atmospheric neutrino background (red) samples (log scale on the right) for events in the
sideband regions. Events with two or three sub-events are selected.

Figure 7.4: Muon momentum distribution for data (black) and simulated atmospheric neutrino background (red) events in the 2-nse sample (log scale on the right)
for events in the sideband regions.

97

Figure 7.5: Muon momentum distribution for data (black) and simulated atmospheric neutrino background (red) events in the 3-nse sample (log scale on the right)
for events in the sideband regions. There are no data entries in these figures since
there are no data events that satisfy the selection criteria.
The ∆X distribution is shown in Figure 7.6 for the two sub-event sample. Eventswith ∆X > 100 cm are selected. The distribution of the time difference between themuon and the gamma is shown in Figure 7.7. Events with 7 < ∆t < 80 ns are
selected and Figure 7.8 presents the gamma momentum distribution. Both figures
are for the two sub-event only since there were no data events in the three sub-event
sample.

Figure 7.6: Distribution of ∆X (defined in Section 5.4.5 ) for data (black) andsimulated atmospheric neutrino background (red) events in the 2-nse sample (logscale on the right) for events in the sideband regions.

98

Figure 7.7: Distribution of ∆t (defined in Section 5.4.6 ) for data (black) and simulated atmospheric neutrino background (red) events in the 2-nse sample (log scaleon the right) for events in the sideband regions.

Figure 7.8: Distribution of gamma momentum for data (black) and simulated atmospheric neutrino background (red) events in the 2-nse sample (log scale on the right)
for events in the sideband regions.
Following the discussion in Section 5.4.8, Figures 7.7 and 7.8 show that the remaining data events in the sideband regions are likely to be events with no deexcitation gamma present.

7.2

Data Results

The final results for the proton decay search are shown in this section following theselection criteria defined in Chapter 5.
99

Figure 7.9 shows the distribution for the number of identified rings in the FCFV,subGeV sample. Events with exactly one identified ring were selected.
The PID distribution is shown in Figure 7.10 for the one-ring sample. The verticalblue line shows the threshold for selecting µ-like events, based on the logarithm of
the likelihood ratio between the muon and electron hypotheses.

Figure 7.9: Number of identified rings distribution for data (black) and simulated
atmospheric neutrino background (red) samples (log scale on the right). Events withexactly one identified ring are selected.

Figure 7.10: PID distribution for data (black) and simulated atmospheric neutrino
background (red) samples (log scale on the right). Events with the logarithm of thelikelihood ratio greater than 30 (blue) are selected as µ-like events.
Figure 7.11 shows the number of sub-events distribution, as described in Section
5.4.3. The analysis is performed differently for events with two and three sub-events.
100

The 2-nse sample is reconstructed with the dedicated µγ-fitter. While, the 3-nse
sample is reconstructed with fiTQun’s standard single-ring µ-fitter.
The distribution of reconstructed muon momentum is shown in Figure 7.12 forthe 2-nse sample and in Figure 7.13 for the two data events in the 3-nse sample.

Figure 7.11: Number of sub-event distribution for data (black) and simulated atmospheric neutrino background (red) samples (log scale on the right). Events with twoor three sub-events are selected.

Figure 7.12: Muon momentum distribution for data (black) and simulated atmospheric neutrino background (red) events in the 2-nse sample (log scale on the right).
Signal region (blue) is defined in the region 215 < pµ < 255 MeV/c.

101

Figure 7.13: Muon momentum distribution for data (black) and simulated atmospheric neutrino background (red) events in the 3-nse sample (log scale on the right).
Signal region (blue) is defined in the region 215 < pµ < 255 MeV/c.
The ∆X distribution (defined in Section 5.4.5) is shown in Figures 7.14 and 7.15
for the two and three sub-event samples, respectively. Only one data event remained
after applying the selection criteria on the reconstructed muon momentum. The
reconstructed value of ∆X for this events was 466 cm, which is outside the signal
region (less than 100 cm).
No data candidate was observed in the three sub-event sample.

Figure 7.14: Distribution of ∆X (defined in Section 5.4.5 ) for data (black) andsimulated atmospheric neutrino background (red) events in the 2-nse sample (logscale on the right). Signal region (blue) is defined in the region ∆X < 100 cm.

102

Figure 7.15: Distribution of ∆X (defined in Section 5.4.5 ) for data (black) andsimulated atmospheric neutrino background (red) events in the 3-nse sample (logscale on the right). Signal region (blue) is defined in the region ∆X < 100 cm.
This concludes the distributions that select single-ring µ-like events with thereconstructed momentum consistent with the expected monochromatic signal fromproton decay. As explained in Chapter 5, CCQE events are abundant in this sample
and additional requirements are necessary to reject the remaining background. In
particular, the selection on the time difference between the muon and the gamma
vertices and the reconstructed gamma momentum are strong discriminators to reject
CCQE background.
Figure 7.16 shows the time difference distribution for the two sub-event sample.
Four events in the data sample satisfy the requirement of ∆t to be greater thanseven nanoseconds. The gamma momentum distribution for these four events isshown in Figure 7.17. The final selection criterion requires the reconstructed gamma
momentum to be in the region between 2 and 10 MeV/c. One event satisfied all theselection criteria for the two sub-event sample and it is the only signal candidate inthis search.

103

Figure 7.16: Distribution of ∆t (defined in Section 5.4.6 ) for data (black) andsimulated atmospheric neutrino background (red) events in the 2-nse sample (logscale on the right). Signal region (blue) is defined in the region 7 < ∆t < 80 ns.

Figure 7.17: Distribution of gamma momentum for data (black) and simulated atmospheric neutrino background (red) events in the 2-nse sample (log scale on theright). Signal region (blue) is defined in the region 2 < pγ < 10 MeV/c.

104

8

Results and Discussions

A search for proton decay in the SUSY favored mode, p → ν̄K + in the prompt-γ
channel, has been performed in the Super-Kamiokande detector. A total of 176.6
kton·years ( 2867.19 days) of fully contained data was used for the search. Data
was collected during the SK-IV period (from October, 2008 to April, 2017) . One
signal candidate event was found with an expectation of 0.30 ± 0.13 events fromatmospheric neutrino background.
With no statically significant excess found, a lower limit on the proton decay
partial lifetime was calculated and the results are presented in this chapter. At
the same time, a sensitivity result is also presented, which allows for comparisonswith previous searches. The detailed description of the calculation can be found in
Appendix A, which follows the Bayesian method of Reference [43].

8.1

Lifetime Limit Sensitivity

A sensitivity study for the lifetime limit is presented in this section. In order tocalculate the lifetime limit sensitivity, an assumption of zero observed candidate
events is made. Given the analysis parameters estimated using signal and backgroundsimulated events, the assumption of zero candidate events characterizes the strongest
limit that can be infered by the analysis.
Using the parameters summarized in Table 8.1 and zero observed events, thesensitivity for the lifetime limit is 2.9 · 1033 years for an exposure of 176.6 kton · years.
The last published result by the SK collaboration for the same proton decay modein the prompt-γ channel used an exposure of 260 kton · years [38]. That analysis
was performed using data of all SK periods and a different reconstruction algorithm.
The published lifetime limit result was 2.5 · 1033 years. The combined result using
other search methods16 was 5.9 · 1033 years.
When comparing the present work with the previous SK search, it is notable thata higher sensitivity was achieved by the result presented here, despite using dataobtained during a smaller exposure. This is a reflection of the significant increasein signal efficiency. Previously, the signal efficiency was 9.1% for SK-IV with an
expectation of 0.13 background events. For searches with a non-negligible number of
expected background events, the sensitivity of the experiment grows with the squareroot of the exposure as opposed to background free searches that grow linearly with
exposure. The lifetime limit obtained in the previous analysis can be scaled to16This includes the other kaon decay mode (π + π 0 ) and a search for an excess of monochromaticmuons in a region around 235.5 MeV/c.

105

the livetime used in this work correcting also for the difference in efficiencies. The
sensitivity scales linearly with the signal efficiency and a rough approximation canbe made to the previous limit such as in Equation 8.1.
r
λnew
new
Expected lifetime = previous lifetime ×
×
,
(8.1)
old
λold
where λ is the exposure and  the signal efficiency. Using 1.56 as the ratio of signal efficiencies and 0.82 as the square root of the ratio of exposures, the expected
lifetime for this search is 3.2 · 1033 years. This quick estimation neglects the difference in the expected number of background events which is one of the causes for thediscrepancy between this estimation and the actual calculated lifetime sensitivity of
2.9 · 1033 years.
This exercise shows the importance of improving data analysis for current andfuture experiments. An increase in exposure by a factor of ten results in an improvement of about 3 in sensitivity, but such an increase can be very expensive. The
SK detector requires relatively low maintenance compared to other experiments, but
running it for 10 years can still be expensive. Building an entire new detector is evenmore expensive and possibly not worth it unless it presents a significant increasein active volume, which makes it even more expensive and reconstruction can be apotential challenge. The same factor of 3 in increase can be obtained by increasing
the exposure by a factor of 4 and the signal efficiency by 50%, saving six years oftime and money.

8.2

Lifetime Limit Result

The estimated signal efficiency using simulated signal events was 14.2% for the combined two and three sub-event samples. The systematic uncertainty for signal efficiency was estimated to be 25.7%. The expected number of background events inthe observed exposure was 0.30 with a systematic uncertainty of 42.2%. One data
candidate event was observed.
This leads to a partial lifetime limit of 1.7 · 1033 years at the 90% confidence level.
For comparison, the world’s most stringent lifetime limit for the same search was the
published result by Super-Kamiokande collaboration of 2.5 · 1033 years using a total
exposure of 260 kton·years [38]. This result is summarized in Table 8.1.

106

Table 8.1: Summary of the proton decay search with selection efficiency and expectedbackgrounds with the correspondent systematic uncertainties.
Exposure
176.6 kton·years
Signal Efficiency
14.2%
25.7%
Signal Syst. Uncertainty
BKG in SK-IV
0.30
BKG Syst. Uncertainty
42.2%
BKG/Mton·year
1.7
Lifetime Limit (90% CL) 1.7 · 1033 years

8.3

Candidate Event

Details of the candidate data event are presented in this section. Figure 8.1 showsthe event display for the candidate. The white ring represents the reconstructedring using fiTQun’s single µ-like ring hypothesis. The single ring fitter reconstructed
the event with a momentum of 235.9 MeV/c. The number os sub-events algorithm
detected 2 sub-events, thus the dedicated µγ-fitter was used to reconstruct the eventas discussed in Chapter 5. The following list summarizes the reconstructed variables
used in event selection, the numbers in parentheses show the selection criteria.
• Distance to the wall: 334.0 cm (dwall > 200 cm).
• Visible energy: 58.3 MeV (30 < evis < 1330 MeV).
• Number of rings: 1 (1 ring).
  • PID likelihood: 309.8 ln LLµe > 30 .
• Number of sub-events: 2 (nse = 2 or nse = 3).
• Muon momentum: 233.7 MeV/c (215 < pµ < 255 MeV/c).
• Distance between muon decay vertex and Michel electron vertex positions:34.3 cm (∆X < 100 cm).
• Time difference between the muon and the gamma: 14.3 ns (7 < ∆t < 80 ns).
• Gamma momentum: 3.9 MeV/c (2 < pγ < 10 MeV/c).

107

In the interest of description, more information (not used for data analysis) about
this particular event is listed below. Figure 8.2 shows the distribution of hit times
for the data candidate event, time of flight corrected using the position of the reconstructed muon vertex in the µγ-hypothesis.
• Vertex position: (x, y, z) = (−361.6, 1306.9, 19.9) cm.
• Muon direction : (êx , êy , êz ) = (0.97, 0.15, −0.17).
• Muon time (with respect to trigger): 935.8 ns.
• Gamma direction : (êx , êy , êz ) = (−0.75, −0.62, 0.19).
• Gamma time (with respect to trigger): 921.4 ns.
• Opening angle (µγ): 150.3◦• Michel electron vertex position: (x, y, z) = (−331.3, −1310.4, 13.1) cm.
• Michel electron direction : (êx , êy , êz ) = (0.44, 0.89, −0.14).
• Michel electron time (with respect to trigger): 3692.2 ns.
• Michel electron time (with respect to the muon time): 2756.4 ns.
• Opening angle (µe): 53.7◦

108Figure 8.1: Event displays of the data candidate event in time (left) and charge(right) modes. The white ring in both figures represent the reconstructed ring in
fiTQun’s single µ-like ring hypothesis. In time mode, the colored scale represents
the different times of each hit and charge mode shows the scale for the integrated
charge (in p.e.) in each PMT.

Figure 8.2: Distribution of hit times (time of flight corrected using the reconstructedmuon vertex position) for the data candidate event. The intense peak around 935 ns
corresponds to the hits coming from the muon, while the small cluster of hits around
920 ns are identified as the gamma hits.
109

A

Lifetime Limit Calculations

The goal of the present analysis is to measure the lifetime, or equivalently, the decayrate of protons. If the number of observed candidate events presented a significant
excess over the expected number of background events, the rate Γ would be calculated
as in Equation A.1.
nsignal
,
(A.1)
Γ=
NT 
where the numerator is the number of observed events, N is the total number ofmonitored particles, T is the running time of the experiment and  the selection
efficiency. The uncertainty of the rate can be found by the standard method of
propagation of uncertainties.
In the absence of an excess of candidate events, an upper limit on the rate is
calculated within some confidence level. The frequentist approach to this problem
lacks the ability to incorporate any uncertainties on the parameters used in thecalculation, namely, the signal efficiency and the expected number of background
events. Since these uncertainties are related to the concept of confidence level, it isnatural to perform the calculation within a Bayesian framework. This work follows
the procedure of Reference [43].
Equation A.2 shows the probability of observing n events given a process withtrue rate Γ, a true exposure λ = N T , a true selection efficiency  and a true numberof expected background events b.
P (n|Γλb) =e−(Γλ+b) (Γλ + b)n
.
n!

(A.2)

The desired parameter, Γ can be obtained by applying Bayes’ theorem as inEquation A.3.
1
(A.3)
P (Γλb|n) = P (n|Γλb) P (Γλb) ,A
where A is a normalization constant to be determined insisting that P (Γλb|n) is
normalized. It must be noted that an implicit assumption was made regarding theconditions of the experiment. Let all the experimental conditions used in the experiment be represented by I. Formally, P (Γλb|n) should be P (Γλb|nI) to indicatethat the calculation of the rate is being performed under such conditions. The conditions I are omitted for clarity, since no ambiguity should arrise from this omission.
Given a particular set of selection criteria, the parameters Γ, λ,  and b are
independent and P (Γλb) can be separated, and Equation A.3 can be rewritten as
P (Γλb|n) =1
P (n|Γλb) P (Γ) P (λ) P () P (b) ,
A
110

(A.4)

where again the experimental conditions are implicit. The quantities P (Γ), P (λ),
P () and P (b) are known as the prior probability distributions (priors) for Γ, λ,  andb, respectively. The priors codify the knowledge of each parameter before the outcomeof the experiment was known. Finally, the joint probability given by Equation A.4
can be marginalized with respect to the so called nuisance parameters, λ,  and b.
The marginalization (integration) is performed since the nuisance variables are notparticularly interesting in this analysis. Instead, the goal is to determine what canbe inferred about the decay rate, Γ, given the experimental observation of n events.
The marginalization procedure is shown in Equation A.5.
Z Z Z −(Γλ+b)e(Γλ + b)n1
P (Γ) P (λ) P () P (b) dλddb.
(A.5)
P (Γ|n) =An!
R∞The normalization constant A is determined by requiring 0 P (Γ|n) dΓ = 1. A
confidence level is obtained by integrating P (Γ|n) from zero to an upper limit Γlim
as in Equation A.6.
Z
Γlim

P (Γ|n) dΓ,

CL =

(A.6)

0

where CL is the confidence level and the lower bound for the lifetime τlim is given by
the inverse of Γlim .
The Bayesian method described above, automatically allows for systematic uncertainties to be incorporated in the calculation. In this work, the prior distributionsfor  and b are assumed to be Gaussian distributions centered at their expected values
from the MC simulated sample and with widths given by their systematic uncertainties. The distribution for the exposure, λ is assumed to have negligible uncertainty
and a Dirac delta distribution is used. Finally, the prior distribution for the rate is
assumed to be uniform for Γ > 0 and zero otherwise. Technically, the distributionsfor  and b are truncated Gaussian distributions since negative values are not allowed.
Explicitly, the distribution for the signal efficiency is given by Equation A.7.
(22e−(−0 ) /2σ , 0 ≤  ≤ 1
P () ∝,
(A.7)
0
, otherwise
where 0 is the estimated efficiency in the MC simulated sample and σ is the estimated systematic uncertainty of this efficiency. The proportionality constant is not
explicitly shown since it can be incorporated in the total normalization constant, A.
For the expected background rate, the situation is similar, but for very low background searches as the one in the present analysis, it is also necessary to account for111

the statistical uncertainty on the expected background rate. Let b0 be the true number of background events evaluated using the un-normalized MC simulated sample,
i.e., not normalized to the data exposure. The statistical uncertainty of b0 is taken
into account by convoluting a Poisson distribution to the truncated Gaussian prioras in Equation A.8.
Z −b0 nbe b0 −(bC−b0 )2 /2σb2
e
db0 ,
(A.8)
P (b) =nb !
for b > 0 and zero otherwise, where nb is the number of background events evaluatedusing the un-normalized MC simulated sample and C is the oversampling factor, i.e.,
the ratio of the exposures for the simulated and the data samples.
Using the priors presented in Equations A.7 and A.8, it is possible to calculate
the confidence limit given by Equations A.5 and A.6. The sensitivity study and thepartial lifetime limit presented in Chapter 8 used a confidence level of 90%.

112References[1] P. Langacker, The Standard Model and Beyond, CRC Press (2009).
[2] W. de Boer, Prog. Part. Nucl. Phys. 33:201-302 (1994).
[3] J.C. Pati and A. Salam, Phys. Rev. D8, 1240 (1973).
[4] H. Georgi and S.L. Glashow, Phys. Rev. Lett. 32, 438 (1974).
[5] P. Langacker and N. Polonsky, Phys. Rev. D47, 4028 (1993).
[6] T. Haines et. al. (IMB Collaboration), Phys. Rev. Lett. 57, 1986 (1986).
[7] W. Gajewski et. al. (IMB Collaboration), Phys. Rev. D42, 2974 (1990).
[8] K.S. Hirata et. al. (Kamiokande Collaboration), Phys. Lett. B220, 308 (1989).
[9] M. Shiozawa et. al. (Super-Kamiokande Collaboration), Phys. Rev. Lett. D81,
3319 (1998).
[10] S.P. Martin, Adv. Ser. Direct. High Energy Phys. 18, 1 (1998).
[11] Myers, Stephen Geneva : CERN, 33 (1991).
[12] S. Raby, Rept. Prog. Phys. 74, 036901 (2011).
[13] Takhistov, V. (2016). Searches for Exotic Baryon Number-Violating Processes at Super-Kamiokande, UC Irvine.
[14] Hyper-Kamiokande Proto-Collaboration Hyper-Kamiokande Design Report arXiv:1805.04163[15] DUNE Collaboration Long-Baseline Neutrino Facility (LBNF) and Deep
Underground Neutrino Experiment (DUNE) Conceptual Design Report Volume 2: The Physics Program for DUNE at LBNF
[16] D. H. Perkins Proton Decay Experiments Annual Review of Nuclear andParticle Science Vol. 34:1-50 (1984).
[17] Super-Kamiokande collaboration, Nucl. Instrum. Meth., A 501 418-462 (2003).
[18] R. Becker-Szendy, et al., Nucl. Instrum. Meth., A324 363 (1993).

113

[19] H. Nishino et al., Nucl. Instrum. Meth., A610:710-717 (2009).
[20] K. Abe, et al., Nuclear Instruments and Methods in Physics Research A737 253(2014).
[21] H. Ejiri, Phys. Rev. C vol. 48, 3 (1993).
[22] R. D. Woods and D. S. Saxon, Phys. Rev. 95, 577 (1954).
[23] T. Yamazaki and Y. Akaishi, Phys. Lett. B 453, 1 (1999).
[24] K. Nakamura et al., Nucl. Phys. A 268, 381 (1976).
[25] Y. Hayato, Nucl. Phys. Proc. Suppl. 112, 171 (2002); G. Mitsuka, AIP Conf.
Proc. 967, 208 (2007); G. Mitsuka, AIP Conf. Proc. 981, 262 (2008).
[26] M. Honda et al., Phys. Rev., D75:043006 (2007).
[27] M. Honda et al., Phys. Rev., D83:123001 (2011).
[28] G. Battistoni et al., Astropart. Phys., 19:269-290 (2003). [Erratum: Astropart.
Phys. 19, 291 (2003)].
[29] G. D. Barr et al., Phys. Rev., D70:023006 (2004).
[30] Shimpei Tobayama An Analysis of the Oscillation of Atmospheric Neutrinos. PhD thesis, THE UNIVERSITY OF BRITISH COLUMBIA (2016).
[31] T. Sjostrand. CERN Program Library CERN-TH-7112-93 (1994).
[32] C. Andreopoulos et al., Nucl. Inst. Meth. A614:87 (2010).
[33] R. Brun, F. Carminati, S. Giani. CERN Program Library CERN-W5013 (1994).
[34] C. Zeitnitz, T. A. Gabriel. Nucl. Instrum. Meth., A349:106 (1994).
[35] R. B. Patterson, E. M. Laird, Y. Liu, P. D. Meyers, I. Stancu, and H. A. Tanaka.
Nucl. Instrum. Meth., A608:206-224, (2009).
[36] S. Berkman et al., T2K-TN-146.
[37] F. James and M. Winkler MINUIT User’s Guide, (2004).
[38] Super-Kamiokande collaboration, Phys. Rev. D 90:072005 (2014).
114

[39] Super-Kamiokande collaboration, Phys. Rev. D 79:112010 (2009).
[40] K. Kobayashi et al, arXiv:nucl-ex/0604006 (2006).
[41] Y. Ashie et al., Phys. Rev., D 71:112005 (2005).
[42] C.M. Marshall et al., Phys. Rev. Lett. 119, 011802 (2017).
[43] R. M. Barnett, et al., Physical Review D54 1 (1996).

115