BOSTON UNIVERSITYGRADUATE SCHOOL OF ARTS AND SCIENCES

Dissertation

A SEARCH FOR DINUCLEON DECAY INTO KAONSUSING THE SUPER–KAMIOKANDEWATER CHERENKOV DETECTORbyMICHAEL D. LITOS
B.S., Michigan State University, 2003

Submitted in partial fulfillment of therequirements for the degree ofDoctor of Philosophy
2010Approved byFirst Reader
Edward T. Kearns, Ph.D.
Professor of Physics

Second Reader
Martin Schmaltz, Ph.D.
Associate Professor of Physics

Acknowledgements
It has been a long and winding road to my PhD, and I certainly couldn’t havemade it all the way to the end without the support of my family, friends, and mentor
figures.
First and foremost, I must thank the kind folks who made me in the first place:
Thank you, Mom and Dad! You provided me with the opportunities that made it
possible for me to even follow a path toward a PhD in the first place. Without yourlove and encouragement throughout my life, I would have never gotten this far.
Next, I must thank my siblings: Doug, Nick, Brian, Melody and Sarah. Especial
thanks to Doug for always being there to put things in perspective with a little doseof the “Doog truth.” It has been an indispensable tool for survival in graduate school
and in life.
I want to thank my step-parents, Marta and Bill, who have been nothing but
supportive of of me throughout the years.
I’d like to thank my B.U. Super-K colleagues and the other fine folks who haveshared workspace with me on the second floor of the PRB over the years: Jen Raaf,Dan Gastler, Fanny DuFour, Wei Wang, Scott Clark, Shantanu Desai, Scott Stackley,
Claire Thomas, Brian Henning, Johnny Penwell, Aaron Herfurth, Betsy Dusinberre,
Jason St. John, Hide Tomita, Luis Ibanez, Kevin Lynch, and Nigel Nation. You guys
were always a big part of why I loved my job.
iii

A huge thank you to all of the great friends I’ve made over the years here at B.U.,
who are simply too numerous to name, though they each and every one deserve far
more than a simple listing in my acknowledgements section, anyway.
I want to thank my Super-K friends from other institutions: Maxim Fechner, Naho
Tanimoto, Parker Cravens, Fumi Kato, Chris Regis, Ittoku Ikeda, Piotr Mijakowski,
Chizue Ishihara, Josh Albert, Gaku Mitsuka, Ryan Terri, Shoei Nakayama, Eric
Thrane, Yumiko Takenaga, Kevin Connolly, and Haruki Nishino.
Thank you to the top tier of the tireless custodial staff of the B.U. Physics Department: Maria and Bart, who always put a shine on our floors and a smile on our
faces.
I want to thank Bob Kingsland in the B.U. machine shop, for teaching me theways of the machine and for just being awesome.
Thanks to Eric Hazen, Chris Lawlor, and Paul Bohn of the B.U. Electronics
Design Facility for helping me in the design and construction of the amazing VORTEX
Monitor System.
Thanks to the fine B.U. Physics front office staff: First, Winna Somers for always
dealing with my messed up travel reimbursement forms. And of course I must thankour own Wonder Woman incarnate, Mirtha Cabello, upon whom the operation of theentire Physics Department depends.
I would like to thank the terrific Super-K faculty for their guidance and tutelage
throughout the years: Jim Stone, Takaaki Kajita, Larry Sulak, Masato Shiozawa,
Hank Sobel, Kenji Kaneyuki, Kate Scholberg, Makoto Miura, Yoshinari Hayato, Shunichi Mine, and Yusuke Koshio.
A huge thanks to my thesis committee: Claudio Chamon, John Butler, Rob Carey,
and especially Martin Schmaltz for helping me grapple with the theory behind myanalysis.
iv

I want to extend a special thanks to Chris Walter for mentoring me during myfirst few years as a graduate student and for all of the continued support throughout
my time on Super-K.
I offer a very large thank you to my advisor, Ed Kearns, who has spent so muchtime, effort, money, and probably hair on my education over the past seven years.
You are responsible for making me the scientist I am today, and for that I thank you.
Finally, I thank my eternal source of motivation and passion: Leah Angstman, myfiance and my best friend. I’d be aimless without you. And now that I’m a doctor, I
promise I’ll make an honest woman out of you!
Oh, and of course this entire thesis would simply be incomplete if I forgot to thankmy best pal Roger Wendell...’s mom. Thanks, Mrs. Wendell!

v

A SEARCH FOR DINUCLEON DECAY INTO KAONSUSING THE SUPER-KAMIOKANDEWATER CHERENKOV DETECTOR(Order No.

)

MICHAEL LITOS
Boston University Graduate School of Arts and Sciences, 2010
Major Professor: Edward T. Kearns, Professor of Physics
ABSTRACT
Dinucleon decay is the simultaneous decay of two nucleons bound within thesame nucleus. This dissertation describes a search for dinucleon decay into kaons
via 16 O(pp) → 14 C K + K + using the first 91.7 kiloton · years of data collected by
the Super-Kamiokande water Cherenkov detector. This is an interesting alternative
to single proton decay because the process violates only baryon number, not lepton
number; single proton decay violates both. In a supersymmetric framework, dinucleon
decay into kaons is predicted to be the most sensitive experimental search that canprobe the magnitude of the R-parity violating parameter λ′′uds .
The complex signature of the signal required the creation of several new event
reconstruction techniques, including tools for fitting Cherenkov rings produced bycharged kaons, a multiple vertex fitter for assigning a unique vertex to each ring,
and a ring classification algorithm that incorporated the overall event topology. Another technique new to the Super-Kamiokande collaboration was the use of a boosted
decision tree, which was utilized in the final stage of the analysis.
No signal event candidates were found in the data. Monte Carlo studies predicteda background of 0.28 ± 0.19 events induced by atmospheric neutrino interactions,and a signal detection efficiency of 12.6% ± 3.2%. Based on these results, the firstviexperimental lower limit on the dinucleon partial lifetime for the decay mode p p →K + K + has been placed at 1.7 × 1032 years at 90% confidence level. This is about twoorders of magnitude better than the only other dinucleon lifetime limits on record,made by the Frejus experiment, which include only decays with final state pions or
leptons.
The lifetime limit may be interpreted as an upper limit of 7.8 × 10−9 on λ′′uds .
This new experimental limit is more restrictive than the limit of ∼ 10−7 estimatedby J.L. Goity and M. Sher in their paper “Bounds on ∆B = 1 couplings in the
supersymmetric standard model” (Phys. Lett. B 343:1-2, 1995), which was based on
typical nucleon and dinucleon lifetime limits of ∼ 1030 years.

vii

Contents
Acknowledgementsiii

Abstractvi

Table of Contents

viii

List of Tables

xiv

List of Figures

xvi

List of Abbreviationsxx

1 Introduction

1

2 Motivation and Theory

6

3 Experimental Apparatus

14

3.1

Super-Kamiokande Detector Overview . . . . . . . . . . . . . . . . .

14

3.2

Structure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

15

3.2.1

Inner Detector (ID) . . . . . . . . . . . . . . . . . . . . . . . .

16

3.2.2

Outer Detector (OD) . . . . . . . . . . . . . . . . . . . . . . .

17

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

17

3.3

viii

3.4

Photomultiplier Tubes . . . . . . . . . . . . . . . . . . . . . . . . . .

21

3.4.1

Inner Detector Photomultiplier Tubes . . . . . . . . . . . . . .

21

3.4.2

Outer Detector Photomultiplier Tubes . . . . . . . . . . . . .

23

Data Acquisition (DAQ) System . . . . . . . . . . . . . . . . . . . . .

24

3.5.1

Inner Detector DAQ . . . . . . . . . . . . . . . . . . . . . . .

24

3.5.2

Outer Detector DAQ . . . . . . . . . . . . . . . . . . . . . . .

27

3.5.3

Event Trigger System . . . . . . . . . . . . . . . . . . . . . . .

28

3.6

Water Purification System . . . . . . . . . . . . . . . . . . . . . . . .

29

3.7

Calibration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

31

3.7.1

Relative Gain Calibration . . . . . . . . . . . . . . . . . . . .

31

3.7.2

Absolute Gain Calibration . . . . . . . . . . . . . . . . . . . .

32

3.7.3

Energy Scale Calibration . . . . . . . . . . . . . . . . . . . . .

33

3.7.4

Relative Timing Calibration . . . . . . . . . . . . . . . . . . .

36

3.7.5

Water Transparency Measurement . . . . . . . . . . . . . . . .

38

3.5

4 Atmospheric Neutrino Monte Carlo

44

4.1

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

45

4.2

Neutrino Interactions . . . . . . . . . . . . . . . . . . . . . . . . . . .

46

4.2.1

Elastic and Quasi-Elastic Scattering . . . . . . . . . . . . . . .

47

4.2.2

Single Meson Production . . . . . . . . . . . . . . . . . . . . .

48

4.2.3

Coherent Pion Production . . . . . . . . . . . . . . . . . . . .

50

4.2.4

Deep Inelastic Scattering . . . . . . . . . . . . . . . . . . . . .

50

4.3

Nuclear Effects . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

52

4.4

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

53

5 Dinucleon Decay Monte Carlo
5.1

55

Intra-Nuclear Simulation . . . . . . . . . . . . . . . . . . . . . . . . .
ix

55

5.2

Extra-Nuclear Simulation. . . . . . . . . . . . . . . . . . . . . . . .

6 Data Reduction Process

58
60

6.1

First Reduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

60

6.2

Second Reduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

61

6.3

Third Reduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

62

6.3.1

Through-Going Muon Cut . . . . . . . . . . . . . . . . . . . .

62

6.3.2

Stopping Muon Cut . . . . . . . . . . . . . . . . . . . . . . . .

63

6.3.3

Cable Hole Muon Cut . . . . . . . . . . . . . . . . . . . . . .

63

6.3.4

Flasher Cut . . . . . . . . . . . . . . . . . . . . . . . . . . . .

64

6.3.5

Accidental Muon Cut . . . . . . . . . . . . . . . . . . . . . . .

64

6.3.6

Low Energy Muon Cut . . . . . . . . . . . . . . . . . . . . . .

65

6.4

Fourth Reduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

65

6.5

Fifth Reduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

66

6.5.1

Final Stopping Muon Cut . . . . . . . . . . . . . . . . . . . .

66

6.5.2

Invisible Muon Cut . . . . . . . . . . . . . . . . . . . . . . . .

66

6.5.3

Final Accidental Muon Cut . . . . . . . . . . . . . . . . . . .

67

6.5.4

Long Tail Flasher Cut . . . . . . . . . . . . . . . . . . . . . .

67

Final Reduced Data Sample . . . . . . . . . . . . . . . . . . . . . . .

68

6.6

7 Physics Reconstruction
7.1

70Standard Reconstruction Algorithm . . . . . . . . . . . . . . . . . . .

70

7.1.1

Vertex Position and Ring Direction . . . . . . . . . . . . . . .

72

7.1.2

Ring Counting. . . . . . . . . . . . . . . . . . . . . . . . . .

75

7.1.3

Showering Likelihood . . . . . . . . . . . . . . . . . . . . . . .

79

7.1.4

Ring Momentum . . . . . . . . . . . . . . . . . . . . . . . . .

84

7.1.5

Ring Number Correction . . . . . . . . . . . . . . . . . . . . .

88x7.1.6

Decay Electron Finding . . . . . . . . . . . . . . . . . . . . .

88

7.2

Multiple Vertex Fitter . . . . . . . . . . . . . . . . . . . . . . . . . .

89

7.3

Reconstructed Ring Variables . . . . . . . . . . . . . . . . . . . . . .

91

8 Dinucleon Decay Signal

92

8.1

Choice of Dinucleon Decay Mode . . . . . . . . . . . . . . . . . . . .

92

8.2

Choice of Kaon Decay Modes . . . . . . . . . . . . . . . . . . . . . .

94

8.3

Cherenkov Ring Combinations . . . . . . . . . . . . . . . . . . . . . .

96

8.3.1

Ring Counting Inefficiencies . . . . . . . . . . . . . . . . . . .

96

8.3.2

True Particle Type Frequency . . . . . . . . . . . . . . . . . .

98

8.3.3

Final List of Ring Combinations . . . . . . . . . . . . . . . . . 100

8.4

Signal Characteristics . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
8.4.1

Multiple Vertices . . . . . . . . . . . . . . . . . . . . . . . . . 102

8.4.2

Kaon Rings . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103

8.4.3

Back-to-Back Kaons . . . . . . . . . . . . . . . . . . . . . . . 104

8.4.4

Features of Kaon Decay Modes . . . . . . . . . . . . . . . . . 104

9 Dinucleon Decay Search
9.1

108

Ring Classification . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109
9.1.1

Two Kaon Ring Classification . . . . . . . . . . . . . . . . . . 111

9.1.2

Kaon Decay Product Ring Classification . . . . . . . . . . . . 112

9.1.3

Pion Classification . . . . . . . . . . . . . . . . . . . . . . . . 115

9.1.4

Low Gamma Momentum Correction

9.1.5

Single Kaon Ring Classification . . . . . . . . . . . . . . . . . 117. . . . . . . . . . . . . . 116

9.2

Event Categorization . . . . . . . . . . . . . . . . . . . . . . . . . . . 119

9.3

Precuts

9.4

Multi-Variate Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . 122. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120

xi

9.5

9.4.1

Input Variables . . . . . . . . . . . . . . . . . . . . . . . . . . 125

9.4.2

Training and Testing of the Boosted Decision Tree . . . . . . . 139

9.4.3

Final Cut Placement for Boosted Decision Tree Output . . . . 142

9.4.4

Boosted Decision Tree Monte Carlo Results . . . . . . . . . . 144

9.4.5

Final Background Characteristics . . . . . . . . . . . . . . . . 145

9.4.6

Final Signal Characteristics . . . . . . . . . . . . . . . . . . . 147

Systematic Uncertainties . . . . . . . . . . . . . . . . . . . . . . . . . 148
9.5.1

Monte Carlo Simulation Uncertainties . . . . . . . . . . . . . . 149

9.5.2

Detector Related Uncertainties . . . . . . . . . . . . . . . . . 150

9.5.3

Reconstruction Uncertainties . . . . . . . . . . . . . . . . . . . 151

9.5.4

Search Method Uncertainties . . . . . . . . . . . . . . . . . . . 151

9.5.5

Total Systematic Uncertainty . . . . . . . . . . . . . . . . . . 152

10 Results and Conclusions

154

A MVFIT Resolutions

157

A.1 Charged Kaon Ring Resolutions . . . . . . . . . . . . . . . . . . . . . 159
A.2 Muon Ring Resolutions . . . . . . . . . . . . . . . . . . . . . . . . . . 161
A.3 Gamma Ring Resolutions . . . . . . . . . . . . . . . . . . . . . . . . 163
B Boosted Decision Tree165

B.1 Final Output . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 166
B.2 Separation Index . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 166
B.3 Adaptive Boost Method . . . . . . . . . . . . . . . . . . . . . . . . . 167
B.4 Pruning . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 168
B.5 Training Distributions . . . . . . . . . . . . . . . . . . . . . . . . . . 168
B.6 Variable Ranking . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 171
xii

C Limit Calculations

174

C.1 Lower Lifetime Limit on Decay Mode . . . . . . . . . . . . . . . . . . 174
C.2 Upper Limit on R-Parity Violating Parameter . . . . . . . . . . . . . 177
Bibliography

180

Curriculum Vitae

185

xiii

List of Tables
3.1

Specifications for the 20-inch ID PMT [1]. . . . . . . . . . . . . . . .

24

3.2

Specifications for the ATM module [1]. . . . . . . . . . . . . . . . . .

25

4.1

Simulated physics processes. . . . . . . . . . . . . . . . . . . . . . . .

54

6.1

Data reduction process event rates. . . . . . . . . . . . . . . . . . . .

69

8.1

Dinucleon decay final state branching ratios. . . . . . . . . . . . . . .

94

8.2

Final state branching ratios for 16 O(pp) → 14 C K + K + . . . . . . . . .

95

8.3

Hadronic inefficiencies for K + detection. . . . . . . . . . . . . . . . .

95

8.4

Frequency of found ring count by final decay state. . . . . . . . . . .

97

8.5

True particle type frequency. . . . . . . . . . . . . . . . . . . . . . . .

98

8.6

Ring combinations by true particle type. . . . . . . . . . . . . . . . .

99

9.1

Frequency of rings produced by kaons. . . . . . . . . . . . . . . . . . 110

9.2

Frequency of rings produced by muons. . . . . . . . . . . . . . . . . . 110

9.3

Frequency of rings produced by gammas. . . . . . . . . . . . . . . . . 110

9.4

Signal, background, and data after applying precuts. . . . . . . . . . 122

9.5

Background after precuts by neutrino mode. . . . . . . . . . . . . . . 123

9.6

Background after precuts by interaction type. . . . . . . . . . . . . . 123

9.7

Division of the Monte Carlo. . . . . . . . . . . . . . . . . . . . . . . . 125
xiv

9.8

Configuration of the boosted decision tree. . . . . . . . . . . . . . . . 141

9.9

Signal and background for different cut placements. . . . . . . . . . . 144

9.10 Signal and background after applying the boosted decision tree cut. . 146
9.11 Final background by neutrino mode. . . . . . . . . . . . . . . . . . . 146
9.12 Final background by interaction type. . . . . . . . . . . . . . . . . . . 147
9.13 True PID of all rings in final background sample. . . . . . . . . . . . 147
9.14 True PID of all rings in final signal sample. . . . . . . . . . . . . . . . 148
9.15 Event categorization purity for the final signal sample. . . . . . . . . 149
9.16 Systematic uncertainties in event reconstruction variables. . . . . . . 151
9.17 Systematic uncertainty of boosted decision tree performance. . . . . . 152
9.18 Systematic uncertainties in event reconstruction variables. . . . . . . 153
A.1 MVFIT reconstructed resolutions for true charged kaon rings. . . . . 159
A.2 MVFIT reconstructed resolutions for true muon rings. . . . . . . . . . 161
A.3 MVFIT reconstructed resolutions for true gamma and electron rings.

163

B.1 Input variables for the boosted decision tree. . . . . . . . . . . . . . . 173

xv

List of Figures
2.1

Feynman diagram for fast proton decay via p → e+ π 0 . . . . . . . . . .

10

2.2

Feynman diagram for dinucleon decay via pp → K + K + . . . . . . . . .

11

3.1

The Super-Kamiokande detector [1]. . . . . . . . . . . . . . . . . . . .

15

3.2

Support structure for mounting photomultiplier tubes [1]. . . . . . . .

16

3.3

An illustration of Cherenkov radiation. . . . . . . . . . . . . . . . . .

18

3.4

Cherenkov light spectrum. . . . . . . . . . . . . . . . . . . . . . . . .

19

3.5

SK event display. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

20

3.6

The 20-inch ID PMT [1]. . . . . . . . . . . . . . . . . . . . . . . . . .

22

3.7

Single p.e. distribution for the 20inch PMTs. . . . . . . . . . . . . . .

22

3.8

Quantum efficiency of the 20-inch ID PMT [1]. . . . . . . . . . . . . .

22

3.9

Timing resolution for the 20-inch PMTs. . . . . . . . . . . . . . . . .

23

3.10 The ID DAQ system [1]. . . . . . . . . . . . . . . . . . . . . . . . . .

26

3.11 The OD DAQ system [1]. . . . . . . . . . . . . . . . . . . . . . . . . .

27

3.12 Overview of the trigger system [1].

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

29

3.13 The water purification system [1]. . . . . . . . . . . . . . . . . . . . .

30

3.14 The relative gain calibration system [1]. . . . . . . . . . . . . . . . . .

32

3.15 The LINAC calibration system [1]. . . . . . . . . . . . . . . . . . . .

34

xvi

3.16 Decay electron spectrum from stopping muons used for calibration.
The points represent the data and the histogram represents the Monte
Carlo [1]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

35

3.17 Left: Calibration using lower energy stopping muons. The double ratioof Monte Carlo over data for the ratio Rp , of measured momentum
to inferred momentum from Cherenkov angle with the Monte Carloprediction [1]. Right: Calibration using high energy stopping muons.
The ratio of the momentum loss to the range. [1] . . . . . . . . . . .

36

3.18 π 0 mass calibration [1]. . . . . . . . . . . . . . . . . . . . . . . . . . .

37

3.19 Summary of energy scale calibration [1]. . . . . . . . . . . . . . . . .

37

3.20 TQ map for a 20-inch PMT. . . . . . . . . . . . . . . . . . . . . . . .

38

3.21 Water transparency laser system [1]. . . . . . . . . . . . . . . . . . .

39

3.22 Attenuation length as a function of time. . . . . . . . . . . . . . . . .

41

3.23 Light scattering laser system [1]. . . . . . . . . . . . . . . . . . . . . .

42

3.24 Attenuation coefficient as a function of wavelength. . . . . . . . . . .

43

4.1

Honda flux model for primary cosmic rays. . . . . . . . . . . . . . . .

45

4.2

Predicted atmospheric neutrino flux at Super-Kamiokande. . . . . . .

46

4.3

Cross sections for quasi-elastic scattering. . . . . . . . . . . . . . . . .

48

4.4

Single pion production cross sections. . . . . . . . . . . . . . . . . . .

49

4.5

Deep inelastic charged current cross sections. . . . . . . . . . . . . . .

52

5.1

Fermi momentum distribution. . . . . . . . . . . . . . . . . . . . . . .

56

5.2

Kinematics of decay protons. . . . . . . . . . . . . . . . . . . . . . . .

57

5.3

Kinematics of outgoing kaons. . . . . . . . . . . . . . . . . . . . . . .

57

5.4

Kaon - nucleon interaction cross-section. . . . . . . . . . . . . . . . .

59

xvii

7.1

TDC fit vertex resolution. . . . . . . . . . . . . . . . . . . . . . . . .

76

7.2

Illustration of the Hough Transformation. . . . . . . . . . . . . . . . .

77

7.3

Example of a Hough Transformed charge distribution. . . . . . . . . .

78

7.4

Example of a showering ring. . . . . . . . . . . . . . . . . . . . . . . .

80

7.5

Example of a non-showering ring. . . . . . . . . . . . . . . . . . . . .

80

7.6

RT OT - momentum look-up table. . . . . . . . . . . . . . . . . . . .

86

8.1

Illustration of an ideal p p → K + K + event. . . . . . . . . . . . . . . . 101

8.2

Kaon flight distances. . . . . . . . . . . . . . . . . . . . . . . . . . . . 102

8.3

True momentum distribution for dinucleon decay kaons. . . . . . . . . 103

8.4

Opening angle between kaon directions. . . . . . . . . . . . . . . . . . 105

8.5

True momentum distribution for dinucleon decay muons. . . . . . . . 106

8.6

True momentum distribution for dinucleon decay gammas. . . . . . . 107

9.1

Variables used in two-kaon classification process. . . . . . . . . . . . . 113

9.2

Illustration of the impact parameter classification variable. . . . . . . 114

9.3

Impact parameter using muon-like and gamma-like vertices. . . . . . 114

9.4

Difference between incorrect and correct impact parameter. . . . . . . 115

9.5

Variables used in the kaon decay product classification process. . . . . 116

9.6

Kaon-like Cherenkov angles for dinucleon decay rings. . . . . . . . . . 117

9.7

Illustration of the linearity variable. . . . . . . . . . . . . . . . . . . . 118

9.8

Difference in linearity between kaons and other rings. . . . . . . . . . 119

9.9Reconstructed variables used in the precuts. . . . . . . . . . . . . . . 121

9.10 Boosted decision tree input variables for final analysis. . . . . . . . . 131
9.11 A single decision tree. . . . . . . . . . . . . . . . . . . . . . . . . . . . 140
9.12 Training, testing, and analysis boosted decision tree output. . . . . . 142
9.13 Significance curve for the boosted decision tree output. . . . . . . . . 143
xviii

9.14 Boosted decision tree output for the Monte Carlo. . . . . . . . . . . . 145
9.15 Boosted decision tree output near the final cut. . . . . . . . . . . . . 145
10.1 Final output of boosted decision tree. . . . . . . . . . . . . . . . . . . 155
A.1 MVFIT resolutions for charged kaons. . . . . . . . . . . . . . . . . . . 160
A.2 MVFIT resolutions for muons. . . . . . . . . . . . . . . . . . . . . . . 162
A.3 MVFIT resolutions for gammas and electrons. . . . . . . . . . . . . . 164
B.1 Input variables used to train the boosted decision tree. . . . . . . . . 169

xix

List of Abbreviations
ADC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Analog to Digital Converter
ATM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Analog Timing ModuleATM ν . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Atmospheric Neutrino
BKG . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . BackgroundCC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Charged Current
CCD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Charge-Coupled Device
CCQE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Charged Current Quasi-Elastic
CL . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Confidence Level
DAQ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Data Acquisition
DIS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Deep Inelastic Scattering
DNDK . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Dinucleon Decay
EM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Electromagnetism
FRP . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Fiber-Reinforced Plastic
GONG . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .GO-NoGo
GUT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Grand Unified Theory
ID . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .Inner Detector
IMB . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Irvine Michigan Brookhaven
KAMLAND . . . . . . . . . . . . . . . . . . . Kamioka Liquid Scintillator Anti-Neutrino Detector
LHC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .Large Hadron Collider
LINAC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Linear Accelerator
MC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Monte Carlo
MSFIT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Muon-Shower Fitter
MVFIT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Multiple Vertex Fitter
NC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Neutral Currentxx

OD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Outer Detector
PDF . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Probability Density Function
PE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Photoelectrons
PID . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Particle identification
PMT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Photomultiplier tube
QAC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Charge-to-Analog Converter
QE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Quasi-Elastic
QTC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Charge-to-Time Converter
SIG . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Signal
SK . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Super-Kamiokande
SK-I . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Super-Kamiokande I
SK-II . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Super-Kamiokande II
SK-III . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Super-Kamiokande III
SM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Standard Model
SNO . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Sudbury Neutrino Observatory
Super-K . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .Super-Kamiokande
SUSY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Supersymmetry
TAC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Time-to-Analog Converter
TDC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .Time-to-Digital Converter
TKO . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Tristan KEK Online
VME . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Versa Module Europaxxi

Chapter 1

Introduction
On a long enough timeline, the survival rate for everyone drops to zero.
–Chuck Palahniuk

Though Mr. Palahniuk’s words certainly hold true for people, it remains a mysteryas to whether or not his maxim can be applied equally well to protons. Protons havenever been observed to decay. So far.
The age of the Universe and the very existence of life on Earth say that the
proton must be stable up to at least cosmological time scales in order for the complex
structure that is observed in the Universe to have arisen and continue to exist to thisday. This begs the question: If the proton lifetime is (at least) cosmologically long,
how can it possibly be measured experimentally? The answer lies in the fact that
protons are indistinguishable from one another, thus statistically identical. So themeasured lifetime will scale with the product of the number of target protons and theamount of time that they are observed. Thus one need not observe a single proton
for an inhumanly long time, one need only observe a large collection of protons foran experimentally reasonable amount of time. This principle is the crux of the basic
1

2

methodology employed by large-scale detectors such as Super-Kamiokande.
The Super-Kamiokande water Cherenkov detector and its direct predecessors,
Kamiokande [2] and IMB [3], have a long history of probing for physics beyond theStandard Model through searches for various modes of nucleon decay [4–7]. The
target volume of water, or fiducial volume, in the center of the Super-Kamiokandetank contains 22.5 kilotons of ultra-pure water, which corresponds to about 7.4 × 1033
protons. As suggested previously, it is important that this number be as large as
possible, because the sensitivity to the lifetime of the proton scales directly with thenumber of protons under observation. The detector is also used for studies of neutrino
oscillations, for which it is also important to have a large target volume in order toincrease the neutrino interaction rate.
One of the most impressive feats accomplished by large-scale water Cherenkov
detectors was to decisively rule out one of the earliest and simplest GUTs, minimal
SU(5), proposed by Georgi and Glashow [8] in 1974. IMB [3] and Kamiokande [2]accomplished this by searching for evidence of single proton decay via the modep → e+ π 0 . That reaction is predicted by minimal SU(5) to occur at a rate that would
yield a proton lifetime of 1031±1 years [9, 10]. IMB and Kamiokande were able toset a lower limit on the proton lifetime that was greater than the SU(5) prediction.
Super-Kamiokande has since set even longer limits on the lifetime of this decay. The
current limit set by Super-Kamiokande is in excess of 1034 years [11]. This limit is
two orders of magnitude larger than the upper bound predicted by SU(5).
This dissertation describes a search for evidence of dinucleon decay, which is thesimultaneous decay of two nucleons bound within the same nucleus. One could imagine such a reaction being mediated by the inter-nucleon exchange of a heavy, nonstandard particle. Only a few experiments have previously searched for dinucleon
decay: Frejus, DAMA, and the BOREXINO Counting Test Facility (CTF).

3

The Frejus detector is a 900 ton tracking iron calorimeter located beneath the alps
in France near the town of Frejus. The Frejus collaboration searched for dinucleon
decay modes involving final state pions or leptons (e.g., pp → π + π + , pp → e+ e+ ).
They did not, however, search for modes involving final state kaons. Their limits forpion and lepton dinucleon decay modes are the only such limits on record, and areon the order of 1030 years [12].
DAMA [13] and the BOREXINO CTF [14] are both liquid scintillator detectors
located in Gran Sasso, Italy. Both experiments have searched for dinucleon decay
into neutrinos or invisible modes by looking for the radioactive decay of an excitednucleus that would be left behind after the disappearance of two nucleons. Frejus [15]was also able to search for the neutrino modes by using the Earth as a target volume
and looking for neutrino interactions in their detector at the correct energy. All threeexperiments had comparable limits, and the current neutrino and invisible dinucleon
decay mode limits are on the order of 1025 years [12].
In this dissertation, the first 91.7 kiloton · years of data taken by the SuperKamiokande water Cherenkov detector are analyzed in the first experimental attempt
to find evidence for dinucleon decay into kaons.
If dinucleon decay into kaons were discovered, it would provide the first evidence
for baryon number violation. This would be very interesting because baryon number isa conserved quantity in the Standard Model. Thus dinucleon decay into kaons would
provide clear evidence for physics beyond the Standard Model. Another reason that
it would be interesting is that baryon number violation is one of the three Sakharov
conditions required to explain the matter/anti-matter asymmetry of the Universe [18]Dinucleon decay into kaon modes is also interesting in the context of supersymmetry, in that the corresponding diagram for the decay requires the presence of baryon
number violating interaction vertices, characterized by the coupling parameters λ′′ijk

4(where i, j, and k are generation indices) which are forbidden under R−Parity.
Thus, the positive detection of dinucleon decay into kaons could provide evidence
for R−Parity violation. If one of the λ′′ijk were shown to be non-zero, it could haveimportant implications in certain supersymmetric models, potentially reducing therequired amount of fine tuning by allowing a lighter bound on the Higgs mass [19].
It has been suggested that a search for dinucleon decay into kaons via p p → K + K +(the mode studied in this analysis) would be the most sensitive experimental probeof the magnitude of one of the λ′′ijk parameters, λ′′112 [20–22].
There are several key differences between this particular search and the moreconventional single nucleon decay searches, such as p → e+ π 0 . The first is that the
p p → K + K + signal included Cherenkov ring patterns created by charged kaons,
which were ejected from the parent nucleus at momenta above Cherenkov threshold.
The second difference is that the signal events all had multiple secondary vertices
due to the subsequent decay of the kaons, separated by distances resolvable by thedetector. These first two distinct features necessitated the development of a new
Cherenkov ring reconstruction method able to fit a unique vertex to each ring in an
event. This multiple vertex fitter is described in detail in Section 7.2.
The next key difference is that this search incorporated event topology in the particle type classification process, which is a technique that has never been incorporated
in a Super-Kamiokande analysis before. Typically, a likelihood analysis predicting
whether a given ring was produced by an electromagnetic showering particle (i.e., e±or γ) or non-showering particle (e.g., p, K ± , π ± , µ± ) is sufficient for particle identification. However, in the case of the p p → K + K + signal, this method alone would notallow one to distinguish charged kaon rings from muon rings produced by the decayof the kaons via the highest branching ratio decay mode K + → µ+ ν. Thus the unique
event topology of the signal was exploited and compared side-by-side with other ring-5

specific reconstructed variables, such as the reconstructed momentum, reconstructed
Cherenkov angle, etc., in order to classify the rings in an event.
The final key difference between this search and searches for single nucleon decay
is that this search scanned over a very large number of discriminating variables,
spanning several event categories by incorporating multivariate analysis techniques.
Unlike single nucleon decay, a reasonable signal efficiency to background ratio could
not be achieved by simply looking at a few variables alone, such as the invariant massand total momentum of the original (di)nucleon decay products. While these twovariables were taken into consideration, many other useful discriminatory variables
were also considered, such as variables describing the event geometry, the overall
symmetry of the event, the kaon decay products, and the goodness (or badness, as itwere) of fitting the event using a single vertex hypothesis. Though this was not thefirst study from Super-Kamiokande to use multivariate analysis techniques (artificial
neural networks have been used before [23]), this was the first study to utilize aboosted decision tree to analyze Super-Kamiokande data.

Chapter 2

Motivation and Theory
There are several motivations behind this search for dinucleon decay into kaons.
One motivation is simply that such reactions are allowed by the conservation of energy, momentum, angular momentum, and charge, and have never been searched for
before. A more interesting motivation is that the observation of dinucleon decay
would provide evidence for baryon number violation. Baryon number (along withlepton number) is a conserved quantity in the Standard Model 1 , and has never been
observed before.
To give this some context, other discrete symmetries once thought to be inviolate
have been found not to hold in certain sectors, such as the symmetry of parity.
Following the suggestions of Lee and Yang [16], Wu (amongst others) found evidence
for parity violation in an experiment that measured the direction of electrons emitted
by the beta decay of 60 Co relative to the direction of spin of the polarized cobalt
atoms [17]. The results of the experiment showed that parity is violated through the
weak interaction, and this discovery brought about a veritable paradigm shift in theworld of physics.
1

An exception to baryon and lepton number conservation arises in non-perturbative electroweak
effects at very high energies [24].

67Such historical precedents make it seem plausible that baryon number may also
not constitute an exact symmetry. This would be fascinating to discover in andof itself, but it could also provide insight into the matter-antimatter asymmetry of
the Universe. Sakharov proposed three conditions that, if satisfied, could explain
the matter-antimatter asymmetry, one of which is the non-conservation of baryon
number [18].
Though no evidence has ever been found for ∆B = 1 baryon number violation
through nucleon decay modes studied at Super-Kamiokande or elsewhere (e.g., p →e+ π 0 , or p → K + ν), this does not preclude that a ∆B = 2 interaction could nottake place. Hence, a search for dinucleon decay into kaons is in fact complimentary
to searches for single proton decay in terms of seeking evidence for baryon number
violation.
Another motivating factor is that searches for dinucleon decay into any purely
mesonic state, such as kaons, provides an interesting experimental alternative tosearches for single nucleon decay, in that the former does not violate lepton number,
yet the latter must in order to conserve angular momentum.
Protons are fermions, which means that they have odd spin. If they decay, theymust decay into at least two lighter particles in order to conserve energy and momentum. The only known particles that are lighter than the proton are the leptons (withthe exception of τ , mτ = 1.8 GeV/c2 ) and certain mesons, such as pions. Leptons are
fermions (odd spin), and mesons are bosons, which have even spin. This means thatany viable proton decay mode must include an odd number of leptons in order toconserve angular momentum. Typical experimental searches focus on decay modes
such as p → e+ π 0 . A two-nucleon system has even spin, however, and therefore does
not require any leptons in the final state (though an even number of leptons would beconceivable). From the perspective of angular momentum conservation, it is perfectly

8reasonable to imagine a two-nucleon system decaying purely mesonically.
Additionally, observation of dinucleon decay into kaons could have interesting
implications in a supersymmetric framework. Supersymmetry (SUSY) represents alarge family of Grand Unified Theories (GUTs) that attempt to unify the strong,weak, and electromagnetic forces. One distinguishing characteristic of SUSY is theintroduction of superpartners to the known particles of the Standard Model (SM). The
superpartners have the same quantum numbers as their corresponding SM partners(e.g., color, electric charge, and hypercharge), but have the opposite spin state.
Under supersymmetry, SM fermions have bosonic superpartners, and SM bosons
have fermionic superpartners. The nomenclature for the superpartners adds an ‘s’
to the beginning of the name of the particle’s corresponding SM partner. The superpartners are symbolically denoted with a tilde placed over the corresponding SM
partner symbol. For example, a Standard Model fermion such as the up quark, u, isgiven a scalar superpartner called the sup, denoted as ũ.
In the minimal supersymmetric Model (MSSM), the number of superpartners produced or annihilated in a given reaction must conserve a symmetry called R−Parity.
SM particles are assigned an R−Parity of +1, and their superpartners are assigned
an R−Parity of −1. The R−Parity value of a particle, PR , is calculated as:
PR = (−1)2S (−1)3B+L ,(2.1)

where S is its spin, B is its baryon number, and L is its lepton number.
R−Parity is a symmetry introduced to the MSSM in order to remove a selectgroup of operators from the superpotential that is allowed by all gauge symmetries.
The terms in the superpotential, W , prohibited under R−Parity are:

9W ⊃ µi Li H̄ + λijk Li Lj EkC + λ′ijk Li Qj DkC + λ′′ijk UiC DjC DkC ,(2.2)

where L, E C , H̄, Q, U C , and D C are lepton doublet, lepton singlet, up-type Higgs,
quark doublet, up-type quark singlet, and down-type quark singlet superfields. µ, λ,λ′ , and λ′′ are the coupling strengths. The ijk indices are generation indices, andgauge indices have been suppressed. The first three terms violate lepton number, andthe final term violates baryon number.
Without their removal through the introduction of R−Parity, the terms in Equation 2.2 would allow fast proton decay to occur via the diagram shown in Fig. 2.1, or
one very similar. The lifetime of the proton would scale as:
1 X 4
τ
∼ 5
m /|λ′ λ′′ | < 1 sec
B.R. p→e+ π+
mp i=2,3 d˜i 11i 11i(2.3)

where mp is the proton mass, λ′ and λ′′ are the R−Parity violating couplings, assumed
to be of order unity, and md˜i , are the squark masses, assumed to be of order 1 TeV [25].
Such a disastrously small prediction for the proton lifetime would ruin the relevanceof the model.
It should be noted here that the diagram shown in Fig. 2.1 does not provide themotivation behind the experimental searches for p → e+ π 0 at experiments like SuperKamiokande. Such searches are primarily motivated by different decay mechanisms
allowed by different models, such as SU(5).
Now, if either the baryon number violating term or the lepton number violating
terms barred by R−Parity were allowed to exist in the model, the proton lifetime
would still be protected. A decay channel like that of Fig. 2.1, for example, would
not be permitted. Any proton decay mode must violate baryon number, as protons
are the lightest baryons, but recall that an odd number of leptons are required in the

10

e+

u
s̃
p

d

λ′′

λ′

ū
π0

u

u

Figure 2.1: Feynman diagram for fast proton decay via p → e+ π 0 in a supersymmetricframework without R−Parity. There is one lepton number violating interaction vertex, λ′ , and one baryon number violating interaction vertex, λ′′ , in the
diagram.

final state in order to preserve angular momentum, thus lepton number must also
be violated. Therefore a model which does not provide both a mechanism to violatebaryon number and a mechanism to violate lepton number will not allow single proton
decay to occur.
It has been noted by Carpenter, Kaplan and Rhee that in certain supersymmetric models which allow either baryon number violation or lepton number violation
through the R−Parity violating couplings (but not both baryon number and lepton
number violation) a lighter bound on the Higgs mass could be allowed, which in turn
could reduce the amount of fine tuning required in the model [19]. This study isconcerned with the models in which the baryon number violating coupling, λ′′ijk , isallowed to be non-zero. If such were the case, one would expect dinucleon decay intokaon modes to arise through a channel such as that shown in Fig. 2.2. In this diagram,the λ′′112 = λ′′uds coupling appears twice. Analogous diagrams exist for p n → K + K 0
and n n → K 0 K 0 . And because dinucleon decay into kaons does not say anything
about lepton number violation, evidence for the signal could potentially be consistent
with—and perhaps even hint at the existence of—one of these R−Parity violating,11

reduced fine tuning SUSY models.

u
K+

u
p

s̄
λ′′uds

d

ũ
u
g̃
u
ũ
p

d
λ′′udsus̄K+u

Figure 2.2: Feynman diagram for pp → K + K + in a supersymmetric framework. There
are two λ′′uds interaction vertices in the diagram.

Yukawa terms in the Lagrangian, L, are derived from the superpotential, W , by
the relation:
1
L=−
2




δ2W
Ψ̄i Ψ̄j + H.C. ,δAi δAj(2.4)

where Ai,j are scalar fields, Ψ̄i,j are fermion fields, and H.C. represents the Hermitian
conjugate [26]. This gives a Lagrangian term for the ∆B = 1 R−Parity violating
term:
1
C
C
˜∗
˜∗
LUiC DjC DkC = − λ′′ijk (ũ∗iR d¯jR dC
kL + dkR ūiR djL + djR ūiR dkL ) + H.C.,
2

(2.5)

where ijk again represent generation indices, ū and d¯ are anti-up and anti-down-type
quark fields, ũ and d˜ are superpartner up and down-type fields, and dC is a down-type12

singlet field [27].
The coupling λ′′ijk UiC DjC DkC requires that the last two fields, DjC and DkC , havedifferent generations (i.e., j 6= k). This happens because the suppressed color indices
are combined anti-symmetrically, but the superfields themselves commute. Thus if
DjC and DkC were swapped and j were equal to k, then the coupling would have tobe equal to −1 times itself, which would be impossible unless it were 0. Thus theintroduction of a second generation field is necessary, explaining why the coupling
violates strangeness by one unit.
For this reason a search for final state kaon dinucleon decay modes, such as p p →K + K + are sensitive to the R−Parity violating SUSY parameters λ′′ijk , but final statepion modes, such as pp → π + π + , are not. A similar argument explains why minimal
SUSY SU(5) favors the single proton decay mode p → K + ν̄, but not p → π + ν̄ (though
in that case R−Parity violation is not required).
Many limits on the R−Parity violating λ′′ijk coupling constants come from experiments that can only measure a limit on a product of at least two different λ′′ijk , or a
product of one λ′′ijk with one of the other R−Parity violating coupling constants. A
summary of experimental constraints on these parameters is given in [27, 28]. Dinucleon decay into kaons, however, would be sensitive to a single coupling alone, λ′′uds ,and has been reasoned to be the most sensitive means of experimentally probing themagnitude of the coupling λ′′uds [20–22].
An upper limit on λ′′uds has been estimated based on typical experimental limits
on nuclear lifetimes from single nucleon decay experiments and the Frejus limits on
lepton pion dinucleon decay modes. Goity and Sher used an approximate figure of∼ 1030 years as the lower lifetime limit for dinucleon decay into kaons, and calculated
an upper limit on the magnitude of λ′′uds to be 10−7 . Appendix C.2 explains thecalculation of the upper limit on λ′′uds as obtained by a lower limit on the p p →13

K + K + partial lifetime. It is the goal of this dissertation to conduct a more sensitive,
direct search for dinucleon decay into kaons in order to provide evidence for dinucleon
decay, or to set a more authoritative experimental bound on the R−Parity violating
parameter λ′′uds .

Chapter 3

Experimental Apparatus
3.1

Super-Kamiokande Detector Overview

Super-Kamiokande is a large water Cherenkov detector located about 1 km beneath the peak of Mt. Ikenoyama in the Mozumi zinc mine near the town of Kamiokain the Gifu prefecture of Japan. Figure 3.1 shows a schematic drawing of the detector
and its surrounds.
The overburden of rock is roughly equivalent to 2700 m of water shielding, lowering
the flux of cosmic ray muons to about 6 × 108 cm −2 s −1 sr−1 . The resultant muon
event rate is five orders of magnitude lower than at that the surface of the earth.
Super-Kamiokande is filled with 50 kilotons of ultra-pure water, making it the largest
detector of its kind in the world.
Data taking began at Super-Kamiokande in April 1996 and continued for five years
until a scheduled shutdown for maintenance and upgrading in July 2001. This initial
five year period of data taking is referred to as Super-Kamiokande I, or simply SK-I.
Regrettably, during the shutdown an accident occurred which required a significant
alteration of the detector for the next data taking period, SK-II, reducing the number
14

15

Figure 3.1: The Super-Kamiokande detector [1].

of ID PMTs to about half that of original design. The SK-II data taking period began
in October 2002 and lasted for three years, ending in October 2005. The detector
has since been restored to full PMT capacity beginning with the SK-III data taking
period, which began in July 2006. This thesis is concerned only with data taken during
the SK-I period, corresponding to 1489.2 days of livetime and 91.5 kiloton · years of
exposure.

3.2

Structure

The water containment structure of Super-Kamiokande is a stainless steel, cylindrical water tank oriented in an upright position. It has a full diameter of 39.3 m and a16height of 41.4 m. The tank is optically segmented into two separate regions: a smaller
concentric cylinder called the Inner Detector and the region between the outer wallsof the tank itself and the outside of the support structure, called the Outer Detector.
These regions are physically separated from one another by a 55 cm thick stainlesssteel support structure on which the photomultiplier tubes are mounted. The tworegions have been optically isolated from one another by covering the support structure with a plastic black sheet. Figure 3.2 illustrates how the photomultiplier tubes
are arranged on the support structure.
Tyvek
8" PMT

Top

Black sheet

Barrel20" PMT

BottomFigure 3.2: Support structure for mounting photomultiplier tubes [1].

3.2.1

Inner Detector (ID)The Inner Detector (ID) is a cylinder 33.8 m in diameter and 36.2 m tall, concentric
to the full Super-Kamiokande tank. It contains 32 kilotons of ultra-pure water. For

17

SK-I, the ID is instrumented with 11,146 inward facing 20-inch photomultiplier tubes
arranged in a uniform grid-like pattern. Roughly 40% of the ID wall surface is thus
photosensitive.

3.2.2

Outer Detector (OD)The Outer Detector (OD) is a 2 m thick region surrounding the Inner Detector on
all sides. It is instrumented with 1,885 outward facing 8-inch photomultiplier tubes.
Each of the OD PMTs is embedded in the center of a 60 × 60 cm 2 wavelength shifting
plate to increase the effective photosensitive area. The OD is used to veto cosmic
ray muon events and tag partially contained physics events. It also acts as a shield
against radioactivity in the rock surrounding the detector. White Tyvek lines thewalls of the OD to maximize reflection, thus increasing efficiency in detecting andidentifying cosmic rays and partially contained events.

3.3

Cherenkov Radiation

The process by which Super-Kamiokande observes and reconstructs particle physics
events is through the detection and interpretation of Cherenkov radiation. Cherenkov
radiation is the light emitted by a relativistic charged particle as it travels through a
medium while its speed remains above that of light in the given medium.
In the case of Super-Kamiokande the medium is pure water, which has an index
of refraction n = 1.33. The requirement for a charged particle to produce Cherenkov
radiation in Super-Kamiokande is:β>1= 0.75,
n

(3.1)

18

c/n

θ
v > c/nFigure 3.3: An illustration of Cherenkov radiation.

where β = v/c, v being the magnitude of the particle’s velocity and c being the speedof light in vacuum. This can be re-expressed in terms of more practical quantities as:

p > 1.14 × m/c,

(3.2)

where p is the particle’s momentum, m is the particle’s mass, and c is the speed oflight in vacuum.
The Cherenkov light is emitted at an acute angle relative to the particle’s momentum vector. This produces a cone-like projection forward along the particle’s
direction of travel, as illustrated in Fig. 3.3. The cone’s half-opening angle, referred
to as the Cherenkov angle θC , is dependent upon the particle’s relativistic speed, β,as such:cos(θC ) =

1
.
nβ

(3.3)

If a particle is highly relativistic such that β ≃ 1, then its Cherenkov angle will
be maximal, θC ≃ 42◦ . This is almost always the case for extremely light particles(i.e., electrons and positrons) produced in particle physics events of the type studied

19

at Super-Kamiokande. More massive particles (e.g., charged kaons) tend to produce“collapsed” cones with smaller Cherenkov angles at typical Super-Kamiokande energy
scales.
The spectrum of light emitted by a Cherenkov radiating particle is given by thefollowing formula:
d2 N
2πα
= 2
dxdλ
λ


1−

1
n2 β 2



,

(3.4)

where λ is the wavelength of light, x is the distance traveled by the particle, and α isthe fine structure constant. In the sensitive region of the photomultiplier tubes usedin the Super-Kamiokande detector (300–600 nm), roughly 340 photons are emitted
per centimeter traversed by the particle. Figure 3.4 shows the relative spectrumof Cherenkov light in pure water is covered well by the sensitive range of SuperKamiokande’s photomultiplier tubes.

Figure 3.4: Relative spectrum of Cherenkov light in pure water and quantum efficiency ofthe 20-inch ID photomultiplier tubes [29].

20

The Cherenkov light that is produced by charged particles in an event at SuperKamiokande travels through the water with little attenuation until it finally reaches
the wall of the ID. There the photomultiplier tubes collect the light and convert it
into an electronic signal which is interpreted as a “hit”. The Cherenkov light cones
produce ring-shaped hit patterns when projected onto the wall, and from these ringpatterns one is able to count the number of charged particles produced in an event, as
well as discern some other crucial information about each of those particles. Figure 3.5
shows an example of a Super-Kamiokande event display.
NUM
4
RUN
9000
SUBRUN
396
EVENT
64753459
DATE **-Jul-15
TIME
13:41:42
TOT PE:
3858.0MAX PE:
24.0
NMHIT : 1456
ANT-PE:
25.4
ANT-MX:
4.2
NMHITA:
30

OD

ID

Figure 3.5: Example of an event at Super-Kamiokande using an “unrolled” event display.
The circular regions at the top and bottom represent the top and bottom
endcaps, respectively. The horizontal rectangle in the middle represents theunrolled side walls of the barrel. Each tiny circle represents a hit PMT, where
the size of the PMTs circle is proportional to the amount of charge observedby that PMT. The smaller display in the corner shows the activity in the OD.
Here, a single ring has been produced by a muon. The hits outside of the ringare from scattered and reflected light.

21

The momentum of the particle is roughly proportional to the total amount of lightin the ring it produced. One may distinguish between an electromagnetic showering
type particle (i.e., an electron, positron, or gamma) and a more massive particle(e.g., a charged kaon or a muon) by taking note of the relative sharpness of the outeredge of the ring. A showering type particle will produce a fuzzier ring pattern thana non-showering type. Particle vertex positions may also be reconstructed by taking
note of the timing of each of the PMT hits as well as their position relative to one
another and relative to the fixed geometry of the tank.

3.4

Photomultiplier Tubes

3.4.1

Inner Detector Photomultiplier Tubes

The 20 inch diameter photomultiplier tubes (PMTs) used in the ID during the SKI experimental period were originally developed specifically for Super-Kamiokande by
the Hamamatsu Corporation. Figure 3.6 shows a schematic drawing of the ID PMT.
These PMTs were designed to have a large photosensitive area, a good timing resolution, the ability to detect a single photo-electron, and long term stability. The“Venetian blind” dynode was adjusted to improve timing resolution and photon collection efficiency, ultimately achieving a timing resolution of about 2.5 ns. Single
photo-electron peaks can clearly be distinguished from dark noise in the PMT signal.
Figure 3.7 shows the single photo-electron ( p.e.) distribution.
The quantum efficiency of the ID PMTs is shown in Fig. 3.8. The PMTs havea significant quantum efficiency over the range of roughly 300–600 nm with a maximum efficiency of 20% at 390 nm wavelength. This range of wavelengths covers the
spectrum of Cherenkov light emitted in pure water quite well. The photocathode of

22~ 70000

cable length
φ 116

φ 254 10φ 82 2glass multi-seal

photosensitive area > φ 460< φ 520

cable

water proof structure( 610 20 )

~ 720

(mm)

Figure 3.6: The 20-inch ID PMT [1].

Quantum efficiency

Figure 3.7: Single p.e. distribution for the 20-inch ID PMT. A single p.e. is equivalent toabout 2 pC at 107 gain. The peak near 0 corresponds to the dark current [1].

0.2

0.1

0

300

400

500

600
700
Wave length (nm)Figure 3.8: Quantum efficiency of the 20-inch ID PMT [1].

23

Figure 3.9: Transit time spread for single p.e. equivalent signals in the 20-inch ID PMT.
The spread is interpreted as the timing resolution of the PMT [1].

the PMTs is Bialkali (Sb-K-Cs) for its high sensitivity to Cherenkov light and low
thermionic emissivity. The transit time spread for a 1 p.e. signal is about 2.2 ns,
as shown in Fig. 3.9. The average dark noise rate is around 3 kHz at the thresholdof 1/4 photo-electrons, which is the threshold value used during normal operational
conditions at Super-Kamiokande. A more complete list of specifications for the ID
PMTs is shown in Table 3.1.
The response of the ID PMTs is sensitive to magnetic fields, so it is critical thatthe 450 mG produced by the geomagnetic field at Super-Kamiokande be reduced.
This is achieved by a system of 26 Helmholtz coils arranged around the detector,
reducing the net magnetic field to approximately 50 mG.

3.4.2

Outer Detector Photomultiplier Tubes

The 8 inch diameter PMTs used in the OD during the SK-I experimental period
were model R1408 Hamamatsu PMTs procured from the decommissioned IMB experiment. The photocathodes of the OD PMTs are each fitted with a 60 cm × 60 cm
× 1.3 cm acrylic wavelength shifter plate, increasing the total OD light collection ef-

24

Product Name
Shape
Photocathode area
Window material
Photocathode material
Quantum efficiency
DynodesGain
Dark current
Dark pulse rate
Cathode non-uniformity
Anode non-uniformity
Transit time
Transit time spread
WeightPressure tolerance

R3600
Hemispherical
50 cm (20 in) diameter
Pyrex glass (4 ∼ 5 mm)
Bialkali (Sb-K-Cs)
22 % at λ = 390 nm
11 stage Venetian blind type107 at ∼ 2000 V
200 nA at 107 gain
3 kHz at 107 gain (1/4 p.e threshold)< 10 %< 40 %
90 nsec at 107 gain
2.2 nsec (1 σ) for 1 p.e. equivalent signals
13 kg
6 kg/cm2 water proof

Table 3.1: Specifications for the 20-inch ID PMT [1].

ficiency by about 50%. The presence of the wavelength shifter plates degrades thetiming resolution of the OD PMTs from 11 ns to 15 ns, but this is still sufficient toallow the OD to function well as a cosmic ray veto counter.

3.5

Data Acquisition (DAQ) System

3.5.1

Inner Detector DAQ

The front end electronics of the ID DAQ are comprised of custom built Analog
Timing Modules (ATMs), into which are fed the analog signals of the ID PMTs.
The ATMs are housed in Tristan KEK Online (TKO) crates located in the four
outer electronics huts sitting on top of the Super-Kamiokande tank. The ATMs
are responsible for recording and digitizing the arrival time and integrated charge
information of each PMT signal via a 12-bit Analog to Digital Converter (ADC).

25

The ATMs each have two separate channels, allowing them to store two successive
PMT pulses when required. This function makes it possible to detect muon-decay
electrons produced by muon decays, and also allows the front end electronics to handle
high rate events that could be generated by nearby supernovae or possibly other rare
phenomena.
Both channels in an ATM consist of a time to amplitude converter (TAC) forrecording the timing information of a PMT pulse, and a charge to analog converter(QAC) for recording the integrated charge of the pulse. The full time range is 1.6 µs
with a resolution of 0.3 ns. The charge saturation level is 600 pC, which corresponds
to approximately 300 p.e., and the charge resolution is 0.2 pC. Table 3.2 contains amore complete list of specifications for the ATM.
Number of channels
One hit processing time
Charge dynamic range
Timing dynamic range
Charge resolution (LSB)
Charge resolution (RMS)
Timing resolution (LSB)
Timing resolution (RMS)
Temperature dependence (QAC)
Temperature dependence (TAC)
Event number
Data size of one hit
FIFO

12 ch/board∼ 5.5 µsec
∼ 400 ∼ 600 pC (12-bit)∼ 1300 nsec (12-bit)
0.2 pC/LSB
0.2 pC (RMS)
0.3 ∼ 0.4 nsec/LSB
0.4 nsec (RMS)
3 Count/deg. ↔ 0.6 pC/deg.
2 Count/deg. ↔ 0.8 nsec/deg.
8 bit
6 Byte2 kByte (∼ 340 hits)

Table 3.2: Specifications for the ATM module [1].

Each ATM typically handles 12 PMTs, and there are (with rare exception) 20
ATMs slotted into each TKO crate. Each of the four outer electronics huts is home
to 12 TKO crates. Along with the ATMs, there is a go/no-go (GONG) trigger
control module and a super control header (SCH) in every TKO crate. The SCH

26

acts as an interfacing device between the TKO crate and its corresponding supermemory partner (SMP), sending the ATM information to the SMP when a trigger isreceived. The SMP sits in a VME crate located within the same hut as the TKO.
The duty of the SMP is to act as a data buffer to be read out by an online computer
workstation. There are two of these workstations in every electronics hut, each ofwhich is responsible for handling 6 SMP modules (half of the number of SMPs per
hut).
Inner DAQ
×6

20-inch PMT

×8• PMT × 11146

SCH
ATM

• ATM × 946

× 20

× 240

ATM
GONG

TKO

• SMP × 48

Ultra-Sparc
interface
SMP
×6• online computer
server × 9
host × 1

onlinecomputer(server)

Ultra-Sparc

trigger signal(to ATMs)HITSUM(from ATMs)

VME

trigger logic

onlinecomputer(host)

VME
interface
TRG

Ultra-Sparc

FDDI

onlinecomputer(server)

NIM

Figure 3.10: The ID DAQ system [1].

The eight total online computer workstations in the four outer electronics huts all
feed into a ninth online workstation which is located in a smaller central electronics
hut. This workstation assembles all of the accumulated information into events andthen sends those events away to an offline computer system for storage and analysis.
Figure 3.10 shows a schematic of the complete ID DAQ system.

27

3.5.2

Outer Detector DAQ

Unlike the ID PMTs, the signals from the OD PMTs are taken from the high voltage lines which provide them power. The high voltage to an OD PMT is distributedby a “paddle card” through a coaxial cable, and the signal is picked off by a capacitor
located on the paddle card.
The signals are fed from the paddle cards into custom built Charge to Time
Converter (QTC) modules, which consist of LeCroy MQT200 chips and comparatorsand comprise the heart of the OD DAQ’s front end electronics. There the signals areconverted into a rectangular pulses each having a width proportional to the input
signal’s charge. There are 5 QTCs distributed over 2 TKO crates within each of theouter electronics huts shared by the ID DAQ electronics. Each QTC reads in signals
from 48 PMTs served by 4 paddle cards (12 PMTs per paddle card).
Local Hitsum

QTC
Modulex40/quadx10/quad

TDC
x5
4 LSBsHigh Voltage

16 bit event #TDC

FASTBUSPMT

Paddle Card

StruckSupply~470/quadrant

FSCC

AUX

84 pin data ribbon cables

Ethernet FSCC Control

Figure 3.11: The OD DAQ system [1].

When the QTC receives a PMT signal, it simultaneously generates a rectangular
HITSUM signal (see Section 3.5.3) that is sent to the global trigger system, which has
a threshold set to 1/4 p.e. If a global trigger signal is received by the OD electronics,
the leading edge and width of the rectangular pulse are converted to the correspondingtime and charge information by a LeCroy 1877 multi-hit TDC module. The TDC
module can record up to 8 QTC output pulses at a resolution of 0.5 ns. The dynamic
range of the TDC module is set to 16 µs before the global trigger time. The digitized28

data stored in the TDC modules are read by a slave computer through a VME
memory module called the dual port memory (DPM) and are then sent to the online
host computer where they are added to the overall event information. Figure 3.11
shows a schematic of the OD DAQ system.

3.5.3

Event Trigger System

The global event trigger system looks for a large number of hit PMTs within a
200 ns window. It does this by checking the added HITSUM pulses from all of thedetector’s ATMs (or QTCs in the case of the OD). The HITSUM pulses are 200 ns
wide signals with a height of 15 mV × the number of hit PMTs belonging to a given
ATM (QTC). Thus, the pulse height of the global HITSUM signal is proportional to
the number of hit PMTs throughout the detector within a 200 ns window.
If the global HITSUM signal for the ID (OD) exceeds 29 (19) hit PMTs, then aglobal trigger signal is issued to the hardware trigger (TRG) module. The TRG module records the trigger type and trigger generation time using a 50 MHz clock, counts
an event number with a 16-bit counter, and generates a global trigger signal thatgets sent to the detector electronics, causing the data to be read out and processed.
Figure 3.12 shows an overview of the global trigger generation scheme.
The trigger system also includes an “Intelligent Trigger” designed specifically forlow energy events, but such events are not considered in this study and thus thistrigger’s description is irrelevant for the purposes of this thesis and will be omitted
here.

29

Analog signal
from PMT
ATM threshold

ADC gate
TDC startand stop
HITSUM

400 nsec width

start

stop

200 nsec width

Global triggerATM
Sum of HITSUM
Master threshold
Global trigger
CentralHut

Figure 3.12: Overview of the trigger system [1].

3.6

Water Purification System

The water in Super-Kamiokande is continuously cycled through a sophisticatedpurification system (see Fig. 3.13) at a rate of about 35 tons per hour, completely
recirculating all of the water in the tank over a period of about 70 days. The purposeof the purification system is to minimize light attenuation and to keep out radioactive
impurities, such as radon, which substantiate much of the background in low energy(<10 MeV) solar neutrino studies.
Maintaining a high level of water purity is critical to one’s ability to accurately
reconstruct physics events, as it increases the collection efficiency of the photomultiplier tubes, thus increasing the amount of data available in any given event. Temporalstability in the water purity level is also a crucial requirement for accurate modellingof the Cherenkov radiating process in the physics event reconstruction algorithms.
The purification cycle begins with a filter to remove dust larger than ∼1 µm,
followed by a heat exchanger which cools the water to 13◦ C. The water is then passedalong to a cartridge polisher to remove any heavy ions. From there, it is passed to30

Figure 3.13: The water purification system [1].

an ultraviolet sterilizer which kills bacteria, and then a reverse osmosis and vacuum
degassifier system for removing radon and other dissolved gasses. Finally the water ispassed through a second filtering stage, removing particles larger than ∼10 nm, andthen on to a membrane degassifier to reduce the radon levels even further, after which
it is returned again to the tank.
The purification system keeps the radon level in the tank to around 0.4 mBq/m3 .
The number density of particles in the water larger than ∼2 µm is approximately
6 /cm3 . The attenuation length of light is ∼100 m at a wavelength of 420 nm. The
resistivity of the water before and after passing through the filtration system is about
11 MΩ·cm and 18.20 MΩ·cm, respectively.

31

3.7

Calibration

3.7.1

Relative Gain Calibration

It is important that the gain of the PMTs be kept as uniform as possible inorder to minimize position and direction dependent systematic errors in reconstructed
momentum. Accordingly, the high voltage is adjusted individually on each PMT tohomogenize the gain throughout the detector.
To calibrate the relative gains of the individual PMTs, a Xe-powered scintillatorball is lowered into the tank at several different fixed positions. The scintillator ball ismade of an acrylic containing BBOT wavelength shifter and MgO powder diffuser. A
Xe lamp followed by an ultraviolet (UV) filter and a neutral density (ND) filter provide
the light to the ball via an optical fiber. The BBOT shifts the light to the visible
range, peaking at 440 nm, and the light is diffused through the ball isotropically into
the detector. The UV light from the Xe source is simultaneously sent to a monitoring
system outside of the tank to monitor the light intensity and to generate a calibration
trigger. Figure 3.14 shows a drawing of the relative gain calibration system.
The relative gain of the ith PMT is given by:

Gi =

ri
1 Qi 2
ri exp( ),f (θi ) Q0
L

(3.5)

where Qi is the charge of the ith PMT, f (θi ) is the angular acceptance of the ith
PMT, which depends on the angle θ between the axis of symmetry of the PMT and
the direction of the incident photon, ri is the distance from the Xe ball to the ith
PMT, L is the attenuation length of the water, and Q0 is a normalization factor.
In the relative gain calibration process, the high voltage of each PMT is adjusted inorder to minimize the spread of Gi over all of the ID PMTs, yielding an overall spread32

UV filter ND filter
Xe Flash Lamp

Photo
Diode

SK TANK

Photo
Diode

Optical fiber

ADC

Scintilator Ball
PMT
Scintilator

Monitor20inchPMT

Trigger

Figure 3.14: The relative gain calibration system [1].

in the relative gain of about 7%.

3.7.2

Absolute Gain Calibration

For measurement of the charge produced by a single photoelectron in a PMT, alow energy source must be used. To serve this purpose, a Cf 252 source surroundedby Ni wire is used. The Cf source and Ni wire are encased in a polyethylene vessel
and lowered into the tank, where the Cf emits neutrons through spontaneous fissionthat are captured by the Ni wire causing the emission of γ rays at various energies,
most commonly 9 MeV. Each γ ray then produces an electromagnetic shower which
in turn generates Cherenkov light that is observed by about 50–80 PMTs in total,

33

with no more than one p.e. per PMT. Single p.e. distributions are created for every
PMT to set the conversion factor from raw observed charge in pC to a number of p.e.
The mean value of the peaks, 2.055 pC, is used as the conversion factor from charge
to p.e.

3.7.3

Energy Scale Calibration

The energy scale calibration is critical for many of the physics studies performed
at Super-Kamiokande. Because of its importance and because of the large range inenergies that are seen in the detector, several different methods are employed in theenergy scale calibration process.
LINAC
For very low energies, a LINAC is used to inject electrons into the detector withenergies ranging from 5 to 16.3 MeV. The LINAC setup can be seen in Fig. 3.15. The
Monte Carlo is tuned to agree with the LINAC data to within 1%.
16

N Decay16

N is used as another low energy calibration source. By lowering a deuterium-

tritium generator into the detector 16 N can be produced in the water. When it decays,
an electron with an energy of up to 4.3 MeV, and a photon with an energy of 6.1 MeV
are produced. The energy scale from this method agrees with the LINAC-tuned
Monte Carlo to within 0.64%.

34

Figure 3.15: The LINAC calibration system [1].

Decay Electrons
All of the cosmic ray muons which stop and subsequently decay in the detector
provide a large sample of muon-decay electrons, which have a well known spectrum
in the range of tens of MeV. The data from decay electrons and the Monte Carloshow the energy scale in the relevant range agrees to within 2.0%.
Stopping Muons
Stopping muons (muons which stop somewhere in the detector) are used for energy
calibration in two different ranges. For muons with momenta of up to approximately
400 MeV/c, the Cherenkov angle can be used to accurately measure the momentum.
This momentum is compared to the momentum determined by the number of photoelectrons collected. The energy scale agreement between data and Monte Carlo for

35

Figure 3.16: Decay electron spectrum from stopping muons used for calibration. The
points represent the data and the histogram represents the Monte Carlo [1].

stopping muons in the range of 200–400 MeV is within 1.5%.
At higher momenta, the Cherenkov angle is indistinguishable from the maximal
value of 42◦ , thus is independent of the muon energy and cannot be used for momentum determination. Instead, the distance traveled by the muon inside the detector isused in the calibration. This distance is calculated by locating the entry point of themuon and the vertex position of the decay electron it produced. A momentum can beinferred from the travel distance, which is compared to the number of photoelectrons
collected. For the energy range 700 MeV to 3.5 GeV, the data and Monte Carlo agree
to within 2.6%.
Figure 3.17 shows the relevant calibration distributions for the two stopping muonenergy ranges.
π 0 Reconstruction
A single π 0 is often produced by neutral current neutrino interactions in SuperKamiokande. The π 0 always decays immediately into two γ s, each of which will

36

Figure 3.17: Left: Calibration using lower energy stopping muons. The double ratioof Monte Carlo over data for the ratio Rp , of measured momentum to inferred momentum from Cherenkov angle with the Monte Carlo prediction [1].
Right: Calibration using high energy stopping muons. The ratio of the momentum loss to the range. [1]produce a Cherenkov ring in the detector. The invariant mass of these γ s is reconstructed and compared between data and Monte Carlo, as shown in Fig. 3.18. These
distributions were made by calculating the invariant mass for all events which were
found to contain two showering (e-like) rings. A Gaussian fit over the range 100–
200 MeV/c2 was then applied to both distributions. This calibration method finds theenergy scale agreeing in the range of 150–600 MeV to within 1%.
Energy Scale Summary
All of the energy scale calibrations are summarized in Fig. 3.19. The different
methods all agree to within about 2%.

3.7.4

Relative Timing Calibration

The relative timing calibration of the PMTs is critical for vertex position reconstruction. The length of the signal cables and the amplitude of the signal itself areboth factors that have an effect on the timing response of the PMTs. The signal am-

37

Figure 3.18: π 0 mass calibration [1].

(MC - DATA) / DATA (%)

SK-I
8
6
4
2
0
-2
µ rangeµ angle0π massdecaye

-4
-6
-8
2

10

10

3

10

4

10

5

10

momentum range (MeV/c)
Figure 3.19: Summary of energy scale calibration [1].

plitude plays a role because of the time-walk effect, where smaller signals reach thediscriminator threshold at a later time relative to their peak than do larger signals.
The hit timing of the PMTs is calibrated using an N2 laser. The laser emits 3 ns
light pulses at 337 nm, which is shifted to 384 nm by a dye laser module. The lightis then carried to a diffuser ball inside the tank by an optical fiber. The light isemitted into a TiO2 diffuser tip inside the diffuser ball, which is made of LUDOX

38

silica gel, all of which allows the light to be very uniformly distributed throughout
the detector without significantly increasing the timing spread. Time-charge (TQ)
maps are created from the laser calibration process and are implemented in the MonteT [ns]

Carlo and applied to the data. A typical TQ map can be seen in Fig. 3.20.
990
980
970
960
950
940
930
920
910
900
890

Q [p.e.]Figure 3.20: TQ map for a typical Super-Kamiokande PMT. Larger time (T) values correspond to earlier hits [1].

3.7.5

Water Transparency Measurement

The water transparency measurement gives the attenuation length of light acrossthe Cherenkov light spectrum observed at Super-Kamiokande. The two effects that go
into the attenuation length calculation are the amount of scattering and the amountof absorption as a function of wavelength. Understanding the attenuation length iscritical to correctly model the light yield in the Monte Carlo. The water transparency
is measured using three different methods: a laser ball, cosmic ray muons, and a laser
fiber.

39

Laser Ball
To measure the water transparency, a Nd:YAG laser is used to pump a Ti:Sa laser,
the light from which is carried via optical fiber to a diffuser ball in the detector. A
CCD camera placed at the top of the tank measures the intensity of the light from thediffuser ball through the water. The light is simultaneously sent along another opticalfiber to an external PMT where the beam intensity is monitored. The intensity of thebeam can be adjusted with a variable attenuation filter, and a range of wavelengths
can be chosen. A drawing of the water transparency laser system can be seen in
Fig. 3.21.
CCD Camera
Laser Box

Nd:YAG

Ti:Sa

DiffuserBall
Beam Splitter
Lens
Integrating
SphereSK Tank

2-inch PMT
Optical Fiber

Figure 3.21: Water transparency laser system [1].

The ball is lowered to various depths and several measurements are taken at several
wavelengths for each depth. After taking data at each depth and each wavelength,
the attenuation length as a function of wavelength L(λ) is calculated by fitting thefunction:

40

log

ICCD
ld
=A−
,
IlaserL(λ)(3.6)

where A is a constant, ICCD is the intensity of the light seen at the CCD camera,
Ilaser is the intensity of the light at the laser monitor PMT, and ld is the depth of the
diffuser ball.
The laser ball measurement is repeated over time to account for variation in theattenuation length correlated with the water purification cycle.
Cosmic Ray Muons
Vertical, downward-going cosmic ray muons which traverse the entire length ofthe detector are also used to measure water transparency. Because these muons are
all highly energetic, they produce a constant number of Cherenkov photons per unit
length traveled in the detector. Thus, the water transparency can be measured bycorrelating the charge observed by each PMT and its distance from the muon track.
The muon track is reconstructed by creating a line between its entrance and exit
points in the detector.
The charge Q observed at each PMT is given by:


l
f (θ)
,
· exp −Q = Q0 ·
l
L

(3.7)

where Q0 is a constant, f (θ) is the PMT acceptance as a function of the incident angleof the photon, l is the path length of the photon, and L is the attenuation length.
One of the advantages of the cosmic ray muon measurement of the water transparency is that it can be done during normal data taking. The time variation of theattenuation length as measured by cosmic ray muons can be seen in Fig. 3.22.

atteunuation length (m)

41

140
120
100
80
60
40

1996

1997

1998

1999

2000

2001

2002

2003

2004

2005

20
0

0

1000

2000
3000
elapse days from 1996/Jan./1

Figure 3.22: Attenuation length as a function of time, attained from downward-goingcosmic ray muons [1].

Light Scattering
The attenuation length is given by:

L = (αabs + αscat )−1 ,(3.8)

where αabs is the coefficient corresponding to absorption and αscat is the coefficientcorresponding to scattering.
To separately measure the effects of scattering in the water, a N2 laser with awavelength selecting dye module is used. Monochromatic light at 337 nm, 371 nm,
400 nm, and 420 nm is produced and injected into the detector via an optical fiber
pointing straight down toward the bottom of the tank. The laser is fired every 6
seconds during normal data taking. Figure 3.23 shows a drawing of the light scattering
laser system.
As can been seen in Fig. 3.23, the barrel of the detector is divided into 5 regions
for this measurement. With this setup, the PMT hits that occur at the top or at thebarrel of the detector are clearly due to scattered or reflected photons. The absorption
and scattering coefficients are tuned such that the PMT hit timing distributions are in

42

LASER
Optical Fiber

Top

337nm
371nm
400nm
420nm

Scattering

Bottom

NUM
2
RUN
8399
EVENT
3919
DATE **-Mar- 1
TIME
5:34:38
TOT PE:
820.0
MAX PE:
19.6
NMHIT :
331
ANT-PE:
13.6
ANT-MX:
2.7
NMHITA:
17

RunMODE:NORMALTRG ID :00010111
T diff.: 2.08
us
:0.208E-02ms
FSCC:
2FF90
TDC0: 8899.2
Q thr. :0.0
BAD ch.: masked
SUB EV : 0/ 0

Figure 3.23: Light scattering laser system [1].

agreement with the light scattering calibration data taken with the laser fiber setup.
Figure 3.24 shows the final results of the attenuation coefficient calculation. The
markers in the plot represent the data from the laser system and the solid and dashedlines represent the model used in the Monte Carlo, which includes the effects ofRayleigh scattering, Mie scattering, and absorption. The Monte Carlo is tuned to fit
the laser data.

43

attenuation coefficient(1/m)

Light scattering measurement

-1

10

-2

10

absorption

Rayliegh scattering-3

10
Mie scattering

300

400

500

600

700

wavelength(nm)Figure 3.24: Attenuation coefficient as a function of wavelength. Markers indicate datataken by the laser system. Solid and dashed lines indicate the model used in
the Monte Carlo [1].

Chapter 4

Atmospheric Neutrino MonteCarlo
The only source of background in this search arises from atmospheric neutrino
interactions taking place within the detector. Atmospheric neutrinos are one of the
final by-products of cosmic rays which strike the atmosphere. Primary cosmic rays
interact hadronically with the nuclei of air molecules in the upper atmosphere tocreate a shower of mesons mostly comprised of pions and kaons. The pions and kaonsthen produce neutrinos as they progress down their respective decay chains. For themost part, these neutrinos pass through the earth without interacting at all. However,a small fraction of them do interact inside the inner detector of Super-Kamiokande at
a rate of about 8 per day. When summed over the total livetime of the experiment,
this tallies up to a substantial amount of background events to sort out, thus it isimportant to understand and model their production and interaction rates.

44

45

4.1

Atmospheric Neutrino Flux

The atmospheric neutrino flux model used by Super-Kamiokande is the Honda
flux model [30]. The flux of primary cosmic rays used in the Honda flux modelis parametrized and fitted to experimental data as shown in Fig. 4.1. The effectsof Earth’s geomagnetic field and solar wind are taken into account in the primarycosmic ray flux calculation, the latter of which can cause fluctuations as large as afactor of two or more for 1 GeV cosmic rays and about 10% at 10 GeV.

Figure 4.1: Flux of primary cosmic rays as predicted by the Honda flux model and observed experimentally. The data points are from the following experiments:
Webber [31] (crosses), LEAP [32] (upward triangles), MASS1 [33] (open circles), CAPRICE [34] (vertical diamonds), IMAX [35] (downward triangles),
BESS98 [36] (circles), AMS [37] (squares), Ryan [38] (horizontal diamonds),
JACEE [39] (downward open triangles), Ivanenko [40] (upward open triangles), Kawamura [41] (open squares), and Runjob [42] (open diamonds).

The simulation model for primary cosmic ray interactions with the nuclei of air
molecules was based on the NUCRIN [43] model for cosmic ray energies below 5 GeV
and DPMJET-III [44] for energies above 5 GeV. The interactions generated in the

46

simulation give a profile of the secondary particles, including pions and kaons. The
atmospheric neutrino flux is then obtained from the decays of these secondary particles.
The final atmospheric neutrino flux at Super-Kamiokande as calculated by Honda,
et al. is shown in Fig. 4.2. The Honda flux was compared to models by Barr [45](Fluka) and Battistoni [46] (Bartol), and all three models show agreement to within
about 10% for neutrino energies up to 10 GeV. The energy range of interest in thisstudy is around the mass of two nucleons, roughly 2 GeV. The overall systematic uncertainty of the absolute atmospheric neutrino flux is about 25% due to uncertainties

νµ + νµ102νe + νe2−2

Flux*Eν (m sec −1sr −1GeV)in the absolute flux and interactions of the primary cosmic rays.

Honda flux

10

Fluka flux
Bartol flux

1 −1
10

1

10

10
Eν (GeV)2

Figure 4.2: Atmospheric neutrino flux at Super-Kamiokande as predicted by the Honda
flux model (solid), Fluka flux model (dashed), and Bartol flux model (dotted).

4.2

Neutrino Interactions

The NEUT [47] model was used to simulate the atmospheric neutrino interactions in Super-Kamiokande. This model includes the following four basic types ofinteraction:

47(quasi-)elastic scattering :ν + N → l + N′

single meson production :ν + N → l + N ′ + meson

coherent pion production :ν + N → l + N ′ + hadrons

deep inelastic scattering :ν +16 O → l +16 O + π

where ν represents the incoming atmospheric neutrino or anti-neutrino of either flavor(µ or e), N and N ′ represent the original and outgoing nucleon (p or n), respectively,
and l represents the outgoing lepton. Note that l can be either a charged lepton in
the case of a charged current (CC) interaction, or a neutrino in the case of a neutral
current (NC) interaction. All of the four interaction types above can be of CC or NC
type.
Fermi momentum and Pauli blocking are taken into account for all of the neutrinonucleon interactions. Fermi momentum of nucleons in 16 O is calculated using the same
technique described in Section 5.1. Pauli blocking is simulated by requiring that the
momentum of the recoil nucleon exceed the Fermi surface momentum of the 16 O
nucleus, taken to be 225 MeV/c in NEUT.

4.2.1

Elastic and Quasi-Elastic Scattering

The V − A theory calculated by Llewellyn-Smith [48] is used to simulate theelastic and quasi-elastic interactions. These events typically produce only one ringfrom the outgoing lepton, or in some cases, a recoil proton with a momentum above
Cherenkov threshold. For this reason, they are a rare background in this search,though occasionally more rings can be produced by such means as multiple scatteringor subsequent hadronic interactions of the recoil nucleon.
The charged and neutral current quasi-elastic scattering cross sections from the
Llewellyn-Smith model are shown in Fig. 4.3. cross sections for both neutrinos and

48

2
1.8
1.6
1.4
1.2
1
0.8
0.6
0.4
0.2
0

-

(a) νµn→µ p

ANL
BNL
SKAT

Serpukhov
GGM (1977)
GGM (1979)

σ (10-38 cm2)

σ (10

-38

2

cm )

anti-neutrinos are shown in the figure.

0

2

4

6
8
10
Eν (GeV)

12

14

2
1.8
1.6
1.4
1.2
1
0.8
0.6
0.4
0.2
0

–

GGM (1977)
GGM (1979)
Serpukhov
SKAT+(b) νµp→µ n

0

2

4

6
8
10
Eν (GeV)12

14

Figure 4.3: Cross sections for quasi-elastic scattering of νµ (left) and ν̄µ (right) used by
NEUT (lines) overlaid on data (markers) from ANL [49], Gargamelle [50–52],
BNL [53], Serpukhov [54] and SKAT [55].

4.2.2Single Meson Production

Cross sections for single meson production are taken from Rein and Sehgal’s
model [56]. This type of interaction becomes significant for neutrino energies above∼1 GeV, where it is possible to form baryonic resonances which decay pionically:ν + N → l + N∗
N ∗ → N ′ + meson ,where N and N ′ are the initial and final state nucleons, respectively, N ∗ is an intermediate state baryonic resonance, and l is the outgoing lepton (charged lepton or
neutrino).
The ∆(1232) resonance is the dominant resonance below 1.4 GeV. For this resonance only, the angular distribution of pions is treated specially by the Rein and
Sehgal model. An isotropic angular distribution is assumed for all other resonances.
Though rare, the production of η and K mesons is also included in the model. ηand K single meson production account for only 4% and 0.6% of all the atmospheric

49

neutrino Monte Carlo events, respectively. Nonetheless, these are important interactions to consider in this analysis, which explicitly searches for rings from charged
kaons and their decay products.
The cross sections for charged current single meson production interactions usedin the Monte Carlo are shown in Fig. 4.4, overlaid on the relevant experimental data.
Single meson ∆ resonance interactions can create events with multiple rings, including rings from charged pions and recoil protons, both of which may be classified
as muon candidates in the analysis (see Section 9.4.5). For these reasons, this typeof atmospheric neutrino interaction was found to be the second largest source ofcm2)

νµ n → µ− n π+

−38

1
0.9

σ(ν µ )(10σ(ν µ )(10

−38

cm2)

background in the dinucleon decay search.

0.8
0.7

1.2

νµ p → µ− p π+

1
0.8

0.6
0.5

0.6

0.4
0.4

0.3
0.2

0.2

0.1

σ(ν µ )(10

−38

cm2)

0

1
0.9

0
1

10

Eν (GeV)νµ n → µ− p π0

1

10Eν (GeV)ANL Radecky, Phys.Rev.D 25, 1161 (1982)

0.8

ANL Campbell, Phys.Rev.Lett. 30, 225 (1973)

0.7

ANL Barish, Phys.Rev.D 19, 2521 (1979)

0.6

BEBC Allen Nucl.Phys.B 264, 221 (1986)

0.5
BEBC Allen Nucl.Phys.B 176, 269 (1980)

0.4
BEBC Allasia Nucl.Phys.B 343, 285 (1990)

0.3
BNL Kitagaki Phys.Rev.D 34, 2554 (1986)

0.2

FNAL Bell Phys.Rev.Lett. 41, 1008 (1978)

0.1
0

1

10

Eν (GeV)Figure 4.4: Single pion production cross sections from the Monte Carlo model (lines)
and experimental data (markers). Top left: charged pion production off ofa proton. Top right: charged pion production off of a neutron. Lower left:
neutral pion production. Lower right: Key for the experimental data points.

50

4.2.3

Coherent Pion Production

When a neutrino scatters off of the entire 16 O nucleus as a whole, rather thanoff of a single one of its constituent nucleons, a pion can be produced in a process
called coherent pion production. The pion will carry away the same charge as theincoming weak current. The nucleus absorbs very little of the momentum exchangedin the interaction, thus the angular distribution of the outgoing pion is peaked in the
direction of the incoming neutrino.
Because the pion never experiences the nuclear environment, only a single pion
is produced in the interaction. More than 1 or 2 rings can be produced by coherent
pion production if the outgoing pion interacts hadronically in the water, though such
events are rare. Because this analysis is concerned only with events that have 3 rings
or more, coherent pion production constitutes a small portion of the total background,
and will not be discussed in detail in this dissertation. For more information on thetreatment of coherent pion production in the atmospheric neutrino Monte Carlo,see [57].

4.2.4

Deep Inelastic Scattering

Two different models are used to calculate cross sections for deep inelastic scattering. A custom model [58] is used for interactions with an invariant mass W of theresulting hadrons from 1.3 to 2.0 GeV/c2 , where only pions are considered as outgoingmesons. Due to their ability to produce several non-showering Cherenkov rings in asingle event, multi-pion production from deep inelastic scattering of atmospheric neutrinos in this energy range was found to be the most prevalent source of backgroundin this search.
For energies in excess of 2.0 GeV/c2 , PYTHIA/JETSET [59] is used, which consid-

51

ers a wider variety of outgoing mesons. In all deep inelastic scattering calculations,
the GRV98 [60] parton distribution function is used, including corrections from Bodek
and Yang [61].
Studies have been conducted at Fermilab [62] using a 15-foot hydrogen bubble
chamber to find the mean number of pions, n̄π , produced by deep inelastic interactions, found to be:n̄π = 0.09 + 1.83 ln W 2 .

(4.1)

The KNO (Koba-Nielsen-Olsen) scaling is used to determine the number of pions ineach deep inelastic scattering neutrino event. To avoid conflicts with the single pion
production model in the overlapping range of W , at least two pions must be produced
to invoke the deep inelastic scattering calculation.
The forward-backward asymmetry of the pion production was studied at the
BEBC [63] experiment, and is given by:
0.35 + 0.4 ln W 2
n̄fπ
=
.
n̄bπ
0.50 + 0.09 ln W 2(4.2)

The charged current cross sections used in the Monte Carlo are shown in Fig. 4.5
along with the relevant experimental data.
Based on experimental results that are reviewed in [64] and [65], neutral current
cross sections were assumed to have the following relationship to the charged current
cross sections:52

σ/Eν (10

-38

2

-1

cm GeV )




0.26
( Eν < 3 GeV )σ(νN → νX)=
0.26 + 0.04 (Eν /3 − 1) ( 3 GeV ≤ Eν < 6 GeV ) ,σ(νN → µ− X) 
 0.30( Eν ≥ 6 GeV )



0.39
( Eν < 3 GeV )σ(ν̄N → ν̄X)=
0.39 − 0.02 (Eν /3 − 1) ( 3 GeV ≤ Eν < 6 GeV ) .
σ(ν̄N → µ+ X) 
 0.37( Eν ≥ 6 GeV )

1
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0

(4.3)

(4.4)

GRV94LOGRV94LO (w/ B-Y corr.)
GRV98LO
GRV98LO (w/ B-Y corr.)

CCFR 90
CDHSW 87
GGM-PS 79

0

10

CHARM 88
BNL 80
CRS 8020

30 50

BEBC-WBB 79
IHEP-JINR 96
IHEP-ITEP 79

CCFRR 84
SKAT

100

200

150

250

Eν (GeV)Figure 4.5: cross sections for deep inelastic charged current interactions. Calculations
are shown for models with and without Bodek-Yang corrections for neutrinos(top curves) and anti-neutrinos (bottom curves). Experimental data points(markers) are overlaid.

4.3

Nuclear Effects

Mesons and nucleons created or scattered inside the 16 O nucleus may interact
hadronically by inelastic scattering, charge exchange, and/or absorption as they prop-

53

agate through the nuclear environment. To simulate these effects, a mean free path
is first calculated using a point of origin determined by the Woods-Saxon distribution(Equation 5.1). A cascade model is then used to simulate inelastic scattering, charge
exchange, and absorption. The results from π − N scattering experiments [66] areused in the determination of the final momentum from inelastic scattering or charge
exchange. Pauli blocking is taken into account for all interactions.

4.4

Detector Simulation

Once the particles have left the parent nucleus, their propagation through thedetector, the emission of Cherenkov light, and the response of the PMTs is handled
by a GEANT3 [67] based Monte Carlo simulation of the Super-Kamiokande detector.
Data acquisition electronics and trigger systems are also simulated in the Monte Carlo.
A list of the simulated physics processes for each particle type is shown in Table 4.1.
Hadronic interactions above 500 MeV are simulated using the CALOR [68, 69] program, and a custom program developed originally for the Kamiokande experiment [58]based on experimental data from π-16 O and π-p scattering [70] is used for hadronicinteractions below 500 MeV. Simulation results are compared to the CALOR package
and the Fluka model to estimate systematic uncertainties for hadronic interactions.
A Poisson distribution with a mean given by Equation 3.4 is used to determine thenumber of Cherenkov photons generated at each wavelength, which are emitted at anangle given by Equation 3.3. Rayleigh scattering, Mie scattering, and absorption as
described in Section 3.7.5 are all simulated for the Cherenkov photons.
PMT response is based on the quantum efficiency curve shown in Fig. 3.8 and thesingle photoelectron distribution.

54

particle type
γ

e±

µ±

hadrons

simulated processes
e+ e− pair production
Compton scattering
Photoelectric effectMultiple scattering
Ionization and δ-ray production
Positron annihilation
Cherenkov radiation
Decay in flight and at rest
Multiple scattering
Ionization and δ-ray production
Direct e+ e− pair production
Nuclear interactions
Cherenkov radiation
Decay in flight and at rest
Multiple scattering
Ionization and δ-ray production
Hadronic interactions
Cherenkov radiation

Table 4.1: Physics processes simulated in the GEANT3 Super-Kamiokande Monte Carlo.

Chapter 5

Dinucleon Decay Monte Carlo
The dinucleon decay mode 16 O(pp) → 14 C K + K + was studied in this analysis.
The reasons for choosing this mode in lieu of other dinucleon decay modes with final
state kaons (i.e., 16 O(pn) → 14 N K + K 0 and 16 O(nn) → 14 O K 0 K 0 ) are described in
Chapter 9. A total of 75,000 events of 16 O(pp) → 14 C K + K + were simulated.
The 16 O(pp) → 14 C K + K + Monte Carlo (also referred to as simply the “p p →K + K + Monte Carlo”) was generated in two stages. In the first stage, the kinematicsof the kaons are determined, and any intra-nuclear interactions that take place before
the kaons leave the parent nucleus are simulated. In the second stage, the kaons
are propagated through the water using the official Geant3 [67] Super-Kamiokande
detector simulation software described in Section 4.4.

5.1

Intra-Nuclear Simulation

In the dinucleon decay Monte Carlo, the Fermi gas model is used to determine themagnitude of each parent nucleon’s momentum. Figure 5.1 shows the distributionof momenta for the s and p wave states that were used. Because the vector sumof the decaying nucleons will always be non-zero, the outgoing kaons will never be55

56

ejected with exactly back-to-back trajectories. The kinetic energy gained from theFermi momentum is relatively small compared to that gained from the rest mass ofthe parent nucleons, however, so the kaons will still have a very large opening angle
between them the great majority of the time.

Figure 5.1: Calculated Fermi momentum distribution used in the Monte Carlo (lines)
overlaid on experimental data points from electron scattering off 12 C [71](markers). Left: p-state. Right: s-state.

The momentum of the kaons can also be modified by an effect called “correlated
decay”. The Monte Carlo assumes that there is a 10% chance of a correlated decay [72] occurring in a given dinucleon decay event. Correlated decay is the namegiven to the decay process in which a third, spectator nucleon becomes involved in
the decay interaction, absorbing a portion of the available energy. The resulting momentum distribution of the kaons and the spectator nucleon is akin to a three-bodydecay. Correlated decay is speculated to happen on theoretical grounds, based on
the overlapping of wave functions. As this process is not well understood, a 100%
uncertainty was taken with this effect when calculating the systematic errors in theanalysis.

57

Figure 5.2 shows the simulated kinematics of the original protons before they
decay, and Fig. 5.3 shows the kinematics of the outgoing kaons from the decay. The
effect of correlated decay can be seen in the tails of the plots in the latter figure.
Proton 1

Proton 1

3000

5000
Proton 2

Proton 2
4000

protons

protons

2500
2000
1500

3000

2000
1000
1000

500
0
0

50

100

150

200

250

0
935

300

940

945

950

955

960

965

970

975

980

total energy (MeV/c2)

momentum (MeV/c)

Figure 5.2: Kinematics of decay protons in simulated 16 O(pp) → 14 C K + K + events. The
momentum (left) and total energy (right) of the first (solid lines) and second(dashed lines) protons.

6000
5000

kaon pairs

kaon pairs

104

3

10

4000
3000
2000

2

1000

10

1000

1100 1200

1300 1400 1500

1600 1700 1800

1900 2000

00invariant mass (MeV/c2)

200

400

600

800

1000

1200

total momentum (MeV/c)

Figure 5.3: Kinematics of outgoing kaons in simulated 16 O(pp) → 14 C K + K + events.
The invariant mass (left), and total momentum (right) are shown.

The final task of the intra-nuclear simulation is to propagate the kaons throughthe nuclear environment. It is possible, for example, that a kaon may undergo charge
exchange (K + n → K 0 p) before exiting the parent nucleus in which the dinucleon decay took place. The kaon is stepped through the nucleus by increments of 0.07 fm, andat each step the cross-sections for charge exchange, elastic, and inelastic interactions
are calculated using the model of Oset, et al. [73].

58

The distance that the kaon travels through the nuclear environment is determinedby the location of the dinucleon decay point within the nucleus. A Woods-Saxon
distribution is used to determine the dinucleon decay point:

ρ(r) =

ρ(0),
1 + exp r−a
b

(5.1)

where ρ(r) is the nuclear density as a function of r, the radial distance from thecenter of the nucleus, a = 1.07A1/3 = 2.69 fm is the maximum radius for 16 O, and2b = 0.82 fm is the thickness of the nuclear surface
Charge exchange of the K + is considered an intrinsic source of inefficiency inthis search. In the simulation, charge exchange occurred inside the parent nuclear
environment with a frequency of 0.05% per kaon.

5.2

Extra-Nuclear Simulation

In the next stage of the Monte Carlo simulation, the kinematic vectors of thekaons calculated from the intra-nuclear simulation are used as the inputs to the official
Geant3-based Super-Kamiokande Monte Carlo software. The kaons are propagatedthrough the detector, and hadronic interactions, ionization energy loss, Cherenkov
light production, subsequent particle decays, and detector response are all simulated.
See Section 4.4 for more details on the Super-Kamiokande detector simulation.
The hadronic interaction cross-section for K + calculated from the simulation isshown in Fig. 5.4. The types of interactions that were simulated are listed in Table 4.1.
The most important types of hadronic interactions that were simulated were elastic
and inelastic interactions, charge exchange, Σ creation, and Λ creation. The latter
three were sources of intrinsic inefficiency in this search, as they destroy the K + state

59

before it has a chance to decay. The simulated dinucleon decay kaons underwent
charge exchange in water at a rate of about 10% per kaon. Σ creation and Λ creation
occurred at a rate of about 1% and .5% per kaon, respectively.

30+

K p PDG+

Cross Section (mb)25

K n PDG++++

avg(K p & K n) pp → K K MC
20

15

10

5

0 -1
5×10

6×10-1

7×10-1

Momentum (GeV/c)

8×10-1Figure 5.4: Kaon-nucleon interaction cross-section calculated from the 16 O(pp) →14 C K + K + Monte Carlo overlaid on K + − p and K + − n cross sections takenfrom the PDG [12].

Chapter 6

Data Reduction Process
The data reduction process separates events caused by interesting physical phenomena from events induced by uninteresting incidental sources, such as radioactivity
in and around the walls or cosmic ray muons. Before reduction, there are on the orderof 106 events (≈ 10 gigabytes of data) recorded per day. Low energy triggers comprisethe majority of these events, occurring at a rate of about 11 Hz. Cosmic ray muons
dominate the high energy triggers, occurring at a rate of about 3 Hz.
The physically interesting events are categorized as either fully contained (FC),partially contained (PC) or upward going muons (UPMU). This search is concerned
only with the FC events, thus only the FC reduction process is described here. FC
reduction is performed by a five step automated process. The description of the PC
and UPMU reduction processes is beyond the scope of this thesis, but may be readabout in further detail in [57].

6.1

First Reduction

The first reduction step aims to remove events from cosmic ray muons and radioactive decays. The cuts in the step must be robust and relatively loose, as it is60

61

applied during the online data taking process.
The requirements for the first reduction are as follows:• P E300 ≥ 200: The number of photoelectrons in the ID that fall within a 300 ns
sliding window in time must be greater than 200.
• NHIT A800 ≤ 50: The number of hits in the OD that fall within a fixed window
in time running from 400 ns before the ID trigger to 400 ns after must be lessthan 50, or else the OD trigger must have been disabled at the time.
• T DIF F ≥ 100 µs: The time interval between the previous event an the current
event must be greater than 100 µs. This is to reject Michel electrons from thedecay of stopping muons.
The event rate is reduced to about 3000 events per day after the first reduction
stage.

6.2

Second Reduction

The aim of the second reduction stage is the same as that of the first with theadditional goal of removing flasher events, which are events caused by flashing PMTs.
Essentially the same cuts are applied, but in this stage they are slightly tighter. The
second reduction process and all later reduction processes are applied offline.
• P Emax /P E300 < 0.5: here PEmax is the maximum number of photoelectrons
collected by a single PMT in a given event, and PE300 is the same quantity as
defined in the first reduction process. This cut removes flasher events.

62• NHIT A800 ≤ 25 if P Etotal < 100000 p.e.: Where P Etotal is the total numberof photoelectrons collected in the ID. As in the first reduction stage, this cut
may be circumvented if the OD trigger was disabled at the time.

6.3

Third Reduction

The goal of the third reduction stage is similar to that of the first two, however thecuts and tools used in this stage are more elaborate than before. The third reduction
stage is broken up into six sub-stages:• through-going muon cut• stopping muon cut• cable hole muon cut• flasher cut• accidental muon cut• low energy muon cut

6.3.1

Through-Going Muon Cut

Through-going cosmic ray muons enter and exit the detector, leaving behind a
substantial amount of charge in the ID along with a cluster of charge at both theentrance and exit points in the OD. A special muon fitter is used to make the following
through-going muon cuts.
Events which satisfy the following criteria are categorized as through-going muon
events and are rejected:

63• P Emax > 250 p.e.: If a PMT in the ID collected more than 250 p.e., then thefollowing cuts are applied.
• mugood > 0.75: The goodness value from the special muon fitter must be
greater than 0.75.
• NHIT Ain ≥ 10 or NHIT Aout ≥ 10: The number of hit OD PMTs inside an
800 ns time window and located within 8 m of the entrance (in) or exit (out)point must be greater than 10.

6.3.2

Stopping Muon Cut

The special muon fitter is again used in the stopping muon cut.
Events which satisfy the following criteria are categorized as stopping muon eventsand are rejected:• {mugood > 0.5 and NHIT Ain ≥ 5} or NHIT Ain ≥ 10: Where the samedefinitions from the through-going muon cut are used.

6.3.3

Cable Hole Muon Cut

This cut removes cosmic ray muons which penetrate the detector through thecable holes used to run cable bundles down to the PMTS. When this happens, little
to no signal is left in the OD, thus making the standard through-going or stopping
muon cuts ineffective. The cable hole muon cuts are as follows:
Events which satisfy the following criteria are categorized as cable hole muon
events and are rejected:• NHIT A > 1: There must be at least one hit in the OD.

64• lveto < 4 m: The distance between the vertex and the cable hole (lveto ) must be
less than 4 m.

6.3.4

Flasher Cut

Flasher events occur when there is a discharge from the dynode of a PMT. They
are characterized by their unusually wide timing distribution.
Events which satisfy the following criteria are categorized as flasher events andare rejected:• NMIN100 ≥ 14 or NMIN100 ≥ 10 if NHITID < 800: Where NMIN100 is
the minimum number of hit ID PMTs in a 100 ns time window, and NHITID
is the total number of hit PMTs in the ID.

6.3.5

Accidental Muon Cut

Cosmic ray muons will occasionally enter the detector shortly after a low energy
event, but within the same trigger gate. Because the light deposited in the OD fromthe muon is relatively late compared to the initial trigger timing, these accidental
muon events require a special cut.
Events which satisfy the following criteria are categorized as accidental muon
events and are rejected:• NHIT Aof f > 20: The number of hit OD PMTs in the 500 ns time window
from 400 ns to 900 ns after the trigger must be greater than 20.
• P Eof f > 5000 p.e.: The total number of photoelectrons collected in the ID in
the same 500 ns time window as above must be greater than 5000.

65

6.3.6

Low Energy Muon Cut

Low energy events caused by radioactive decays and noise in the electronics are
removed in this step.
Events which satisfy the following criteria are categorized as accidental muon
events and are rejected:• NHIT50 < 50: The number of hit ID PMTs within a 50 ns sliding window intime must be less than 50. The residual hit time of each PMT is used, whichtakes into account the time of flight of the photons from a single vertex point.
The vertex is defined here as the point which most strongly peaks the residual
time distribution. 50 PMT hits roughly correspond to a 9 MeV event.
At the end of the third reduction stage, the event rate has been reduced to 45
events per day.

6.4

Fourth Reduction

The fourth reduction stage is dedicated to removing the remaining flasher events.
An important characteristic of flasher events is that the light pattern created in these
events tends to be repeated over long periods of time. This repetition is exploited in
the fourth reduction process, which is explained in further detail in [57].
The event rate has been reduced to about 18 events per day after applying thefourth reduction stage.

66

6.5

Fifth Reduction

The fifth and final stage in the FC reduction process is designed to remove various
remaining cosmic ray muon and flasher events. It consists of four steps:• final stopping muon cut• invisible muon cut• final accidental muon cut• long tail flasher cut

6.5.1

Final Stopping Muon Cut

This cut is similar to the stopping muon cut applied earlier during the third
reduction stage, but in this step the entrance point is now computed by extrapolating
backward along the fitted track of the muon, rather than using the earliest hit PMT.
Events which satisfy the following criteria are categorized as stopping muon eventsand are rejected:• NHIT Ain ≥ 5: Where NHIT Ain is defined as before.

6.5.2

Invisible Muon Cut

This cut targets muons which are below Cherenkov threshold (hence “invisible”),
but emit a decay electron in the detector which can be seen.
Events which satisfy the following criteria are categorized as invisible muon eventsand are rejected:• P Dtotal < 1000: The total number of photoelectrons in the ID is less than 1000.

67• NHIT Amax
early > 4: The maximum number of hit OD PMTs in a 200 ns time
window allowed to move from 100 ns before the trigger to 8900 ns after thetrigger is greater than 4.
• NHIAmax
early + NHIT A500 > 9 if lcluster < 500 cm or else NHIT Aearly > 9:
Where NHIT A500 is the number of hit OD PMTs in a 500 ns time window
from −100 ns to 400 ns, and lcluster is the distance between the two OD clusters
used during the calculation of NHIT Aearly and NHIT A500 .

6.5.3

Final Accidental Muon Cut

The final accidental muon cut further eliminates muons which accidentally sneakin after a low energy event.
Events which satisfy the following criteria are categorized as accidental muon
events and are rejected:• P E500 < 300 p.e.: The total number of photoelectrons in the ID within a 500 ns
time window from −100 ns to 400 ns must be less than 300.
• NHIT Amax
late > 20: The number of hit OD PMTs in a 200 ns sliding time window
allowed to move from 400 ns to 1600 ns after the trigger must be greater than
20.

6.5.4

Long Tail Flasher Cut

The long tail flasher cut is similar to the flasher cut applied during the third
reduction phase, but with tighter cut values.
Events which satisfy the following criteria are categorized as long tail flasher eventsand are rejected:

68• NMIN100 > 5 if GOODpoint < 0.4: The minimum number of hit ID PMTs ina 100 ns sliding time window allowed to move from 300 ns to 800 ns after thetrigger gate must be greater than 5 if the goodness of the point fit, GOODpoint,
is greater than 0.4.
After applying the fifth reduction process, the event rate has been reduced toabout 16 events per day.

6.6

Final Reduced Data Sample

Approximately 16 events per day remain in the data after applying all five of the
FC reduction processes. Three final cuts are applied to reduce the data even further
to create the final sample:• NHIT AC < 10: The number of OD PMT hits in the cluster with the highest
charge must be less than 10.
• EV IS > 30 MeV: The total amount of visible energy must be greater than
30 MeV, where the visible energy is defined as the sum of the energy of eachring when assuming all rings to have been produced by electrons.
• W ALL > 200 cm: The distance from the vertex point to the nearest wall must
be greater than 200 cm. This is the fiducial volume (FV) cut.
The above cuts reduce the final event rate in the data to about 8.2 events per day.
The contamination from undesired backgrounds in the data at this stage, such ascosmic ray muons and flashers, is estimated to be negligible (< 1%). The estimatedefficiency of the reduction process for both atmospheric neutrino events and dinucleon

69

Reduction Stage
Pre-Reduction
First Reduction
Second Reduction
Third Reduction
Fourth Reduction
Final Reduction (FV)

Events per Day
1,200,000
3,000
200
45
16
8.2

Table 6.1: Event rates in the data after each stage of the fully contained reduction process.

decay events is about 99.9%. Table 6.1 summarizes the event rates at each stage ofthe reduction process.

Chapter 7

Physics Reconstruction
Reconstruction on a ring-by-ring basis was done in a two-step, bootstrapping
process. First, the standard Super-Kamiokande event reconstruction algorithm wasapplied. From this, the reconstructed event vertex, ring directions, and Cherenkov
angles were used as inputs for the second step of ring reconstruction, which utilized
a custom built algorithm suitable for identifying multiple particle vertices in a single
event. The name of the second reconstruction algorithm is MVFIT, for multiple
vertex fitter.

7.1

Standard Reconstruction Algorithm

The standard Super-Kamiokande event reconstruction software is made of several
specialized algorithms that combined make up the full reconstruction process. The
algorithms are listed below in sequential order:(1) Vertex Position and Ring Direction
Calculates a single event vertex based on the point where the residual timing
distribution of hit PMTs has the sharpest peak. Also calculates the direction70

71

and Cherenkov angle of the dominant ring.
(2) Ring Counting
Additional rings are searched for using a likelihood method. The event vertex
position and direction of the dominant ring calculated in the previous step areused in the calculation. The total number of rings is determined in this step.
(3) Showering Likelihood
A showering likelihood variable is calculated for all rings in the event based on
the ring’s opening angle and a comparison between the observed light pattern
and the expected light patterns for showering and non-showering particle types.
(4) Ring Momentum
The momentum of each ring is calculated by comparing the amount of charge
accumulated inside a 70◦ cone with respect to the ring’s direction. A look-up
table based on MC simulations and detector calibration is used for the conversion
from charge to momentum. The observed charge in each PMT is also divvied up
into separate rings in this step.
(5) Ring Number Correction
Possible fitting mistakes are corrected by eliminating rings of very low momentum
that overlap with other more energetic rings.
(6) Decay Electron Finding
Rings from decay electrons following the primary event are identified.

72

7.1.1

Vertex Position and Ring Direction

There are three parts in the initial event vertex position and dominant ring fitting
procedure: (1) point fit, (2) direction fit, and (3) TDC fit.
Point Fit
The event vertex position estimation is performed with the assumption that allof the light observed in the detector was emitted simultaneously from a single point
source. To find the point that best matches this hypothesis, a three-dimensional gridof test points throughout the volume of the ID is considered, and the single best ofthese points is kept. Then the process is repeated with a finer granularity in a small
region around the best point of the previous granularity level.
The best test point in the grid is determined by maximizing a goodness calculation,
which is based on the hit times of the PMTs. The point fit goodness is defined asfollows:(ti − t̄)2
1 X
exp −
Gpoint =Nhit i2 (a × σt )2

!

,

(7.1)

where Nhit is the total number of hit PMTs, σt is the timing resolution of the PMTs
(2.5 ns), a = 1.5 is a factor to crudely account for light scattering, and t̄ represents
the average value of the residual times of the PMTs, ti , which are defined:

ti = t0i −di,v (di , qi )(7.2)

where t0i is the recorded absolute time of the PMT hit, di is the distance from the test
vertex point to the PMT, and v is the effective velocity of light in water as a functionof di and the recorded charge, qi , accounting for wavelength and acceptance effects.

73This vertex fitting algorithm is guided by the assumption that all light in theevent was originated simultaneously from the same point. In cases where this is
not true, however, the reconstructed event vertex found by the point fit algorithm
will tend to approximate an “average” event vertex (e.g., a point along the trackover which a single muon radiates Cherenkov light, or a point between two spatiallyseparated muons). The point fit vertex can still be safely used as a springboard formore accurate vertex fitting algorithms applied later in the reconstruction process,however.
Direction Fit
In the next step, the direction and Cherenkov angle of the dominant ring is estimated. The initial direction guess is found by calculating a charge-weighted vector
sum of all light in the detector using the event vertex found by the point fit in theprevious step. From there, a directional test grid in (θ, φ) space is generated, analogous to the test grid in the vertex point fit. A goodness calculation is again used, thistime to find the best direction in the grid, along with the best Cherenkov opening
angle.
The goodness used in the directional fit is defined as follows:

Gdir =R θC0
Q(θ)dθ
(θC − θmax )2,
exp −sin θC
σa2(7.3)

where θC is the test opening angle (allowed to vary), Q(θ) is the charge distribution
as a function of angle relative to the test direction, θmax = 42◦ is the maximum
Cherenkov angle assuming the particle’s velocity β = 1, and σa is the estimated rms
spread of PMT hits around θC . The test direction and opening angle which yield the
maximum goodness are taken as the best fit.

74

TDC Fit
The purpose of the TDC fitting step is to more precisely fit the vertex position bytaking into account the finite track length of a massive charged particle and the effectof indirect light caused by scattered Cherenkov photons. The vertex, ring direction,
and Cherenkov angle of the dominant ring found in the point fit and direction fit
steps are used as inputs to this step.
In this step, a modified residual time calculation for the PMT hits inside theCherenkov cone is used:
n ~
1 ~
~
~
ti = t0i − |X
i − O| − |Pi − Xi |,
c
c

(7.4)

~ is the test vertex position, X~ i is the position along the track where photons
where O~i iswould be emitted toward the i-th PMT, n is the index of refraction for water, P
the position of the i-th PMT, and t0i is the recorded absolute time of the hit of thei-th PMT.
The goodness for PMTs inside (GI ) and outside (GO ) the Cherenkov cone arecalculated slightly differently, as shown below:!
(ti − t̄)2
,
GI =
exp −
σi2
2 (a × σ̄)2
i
"
!
#

X 1
(ti − t̄)2ti − t̄max exp −,GO =, 0.8 exp −
σi2
20 ns
2 (a × σ̄)2
i
X 1

(7.5)
(7.6)

where σi is the timing resolution of the i-th PMT, σ̄ is the average PMT timing
resolution, a = 1.5 is a factor to account for indirect scattered light, ti is the residual
time of the i-th PMT, t̄ is the average residual time of the PMTs, and 20 ns is theaverage time difference between direct and scattered light. The GI goodness value

75

and the GO goodness value are calculated separately for PMT hits inside (θi < θC or
ti < t̄) and outside (θi > θC or ti > t̄) of the Cherenkov opening angle, corresponding
to the direct and scattered light, respectively.
The final overall goodness for the test vertex position is then calculated using both
GI and GO , as shown below:

Gtotal =GI + GO
P
.
i σi

(7.7)

The track length used in the calculation is estimated by summing the charge
within 70◦ of the direction of the ring and then calculating the corresponding muon
momentum required to produce that amount of light. Using this track length, thetotal goodness, Gtotal , is maximized as a function of vertex position and ring direction. Then the procedure is repeated, calculating a new track length and varying thevertex position and ring direction to maximize the total goodness once again. Thisis repeated in an iterative process until a final, stable fit is reached.
The resolution of the TDC fit vertex, determined by the distance between thereconstructed vertex position and the true vertex position of Monte Carlo events, isshown in Fig. 7.1 for various subsamples of single-ring atmospheric neutrino events.
The resolution falls between 50 cm and 90 cm, depending on the sample.

7.1.2

Ring Counting

After fitting an event vertex and identifying the dominant ring, the ring counting
algorithm is applied to search for other possible rings in the event. There are twobasic steps to the ring counting procedure: (1) search for ring candidates and (2) test
the best candidate found in step (1). If the ring candidate passes the test in step (2),
it is removed from the search map and the process is repeated from step (1) again.

76

1600

1400

1400
Number of events

Number of events

1200
1000
800
600
400

1000
800
600
400

200
0

1200

200
0

50

100

150

200

0

0

Vertex displacement (cm)

50

100

150

200

Vertex displacement (cm)

350

250

300

Number of events

Number of events

400
300

200
150
100

200
150
100

50
0

250

50
0

50

100

150

Vertex displacement (cm)

200

0

0

50

100

150

200

Vertex displacement (cm)Figure 7.1: Distance between the TDC fit reconstructed vertex and the event true vertex
for sub- GeV (<1.33 GeV) and multi- GeV (>1.33 GeV) single-ring atmospheric
neutrino Monte Carlo events.

This continues until all ring candidates fail the test in step (2).
Ring Candidate Search
To search for ring candidates, a Hough transform is applied to the PMT hits. This
involves first mapping the PMT hits in (θ, φ) space relative to the event vertex anddominant ring direction. Then, the charge in each hit is distributed uniformly along
a virtual circle corresponding to a projection of a cone with a 42◦ half-opening angle
centered on the position of the hit PMT itself. This results in the accumulation of
charge at the center of a true Cherenkov ring, where the Hough-transformed charge

77

distributions overlap the most. See Fig. 7.2 for an illustration of this concept.
42 deg. ringhit PMT(possible center)

center(most probable)

Cherenkov ring

Figure 7.2: The basic concept of the Hough Transform. The charge in a PMT is distributed along a circle corresponding to a 42◦ Cherenkov opening angle as seenfrom the event vertex. The circles of the hit PMTs from the same Cherenkov
ring will then overlap in the center of the actual ring.

In practice, virtual circles are not actually drawn out, but instead an expectedcharge distribution function f (θ) is mapped out for each hit PMT, weighted by theobserved charge distribution. When plotted in a 2-D histogram, the ring directions
become visible as peaks of charge as shown in Fig. 7.3.
Ring Candidate Test
After identifying the candidate rings, they are tested through the use of a loglikelihood method. If there are N rings already confirmed to be good candidates inthe event, then the (N + 1)-th candidate ring is tested for its validity.
The ring candidate test likelihood is defined as follows:

F =

5
X
i

{log [(Pi )N +1 ] − log [(Pi )N ]} ,(7.8)

where Pi is the probability density function for the i-th evaluation function. So,

78

600
400
200
0
150

Θ(

de 100
gr
ee 50
s)

300
200
100
00

ees)

gr
Φ (de

Figure 7.3: A two dimensional histogram of a Hough transformed charge distribution.
Two distinct charge peaks can be seen, corresponding to the directions of two
candidate Cherenkov rings in the event.

PN +1 and PN are the probability functions for the event to have (N + 1) and N rings,
respectively. Note that the algorithm is limited to a maximum of five rings in an
event.
The following quantities are used in determining the evaluation functions:• The difference of L(N + 1) − L(N)• The average value of the expected charge from the (N + 1)-th ring near theedges of the N confirmed Cherenkov rings• The average value of the expected charge outside the (N + 1)-th ring

79

• The residual charge left after subtracting the expected charge from the N confirmed rings• The difference between the charge in peak of the candidate ring and the averagecharge inside and outside the candidate ring
Here, the function L(N) is defined:

L(N) =

Xlog prob qiobs ,

i

N
X
n=1

exp
αn · qi,n

!!

,

(7.9)

where prob(qiobs , qiexp ) is the probability of detecting qiobs given the expected value of
qiexp . The probability density functions were determined using Monte Carlo studies.

7.1.3

Showering Likelihood

The showering likelihood algorithm is intended to separate rings into two different types: a showering type and a non-showering type. The former type describes
the diffused ring patterns created by electrons, positrons and gamma-rays due to theeffects of electromagnetic showering and multiple scattering. See Fig. 7.4 for an example of how a showering ring appears in the event display. The latter type describes
the sharper-edged ring patterns generated by more massive charged particles, suchas muons and charged mesons, including kaons. See Fig. 7.5 for an example of how anon-showering ring appears in the event display.
The Cherenkov opening angle provides another feature of distinction between
showering and non-showering rings. The Cherenkov angle is expected to always be
maximal (42◦ ) for the lighter, showering type particles, but it may be smaller for more
massive particles if they are not traveling at highly relativistic speeds or after they
have lost energy through ionization. The fuzziness or crispness of the light pattern

80

Figure 7.4: Example of a showering ring.

Figure 7.5: Example of a non-showering ring.

81

and the Cherenkov opening angle are both exploited simultaneously by the showering
likelihood algorithm.
The showering likelihood calculation depends on the ability to formulate expectedcharge distributions for the two ring type hypotheses. The details of those expectedcharge distributions will be explained first, followed by the description of the likelihood itself.
Expected Charge Distributions
The expected charge distributions describe the amount of charge expected to be
seen in in each PMT given that the ring was produced by an electron (q exp (e)) torepresent the showering type, or a muon (q exp (µ)) to represent the non-showering
type. They are defined as:

qiexp (e) = αe Qexp (pe , θi )qiexp (µ) =where:

 1.5
1Rf (Θi ) + qiscat ,riexp( rLi )sin2 θxi
 + qiknock
αµdθri sin θxi + ri · dx
|x=xi

!

1
f (Θi ) + qiscat ,exp( rLi )(7.10)(7.11)

82

αe , αµ : normalization factors
ri : distance from the vertex to the i-th PMT
θi : opening angle between the i-th PMT direction
and the ring direction
L : light attenuation length in waterf (Θi) : correction for the PMT acceptance as a functionof the photon incidence angle Θi
R : radius of the virtual sphere (16.9 m)
Qexp (pe , θi ) : expected p.e. distribution from an electron as afunction of the the electron momentum and theopening angle
x : position of the muon along its track
xi : position of the muon along its track where
Cherenkov photons are emitted toward the i-th
PMT
qiscatt (qiknock ) : expected p.e.s for the i-th PMT from scatteredphotons (knock-on electrons)θ (θxi ) : Cherenkov opening angle of the muon at track
position x (xi )The expected p.e. distribution distribution for an electron, Qexp (pe , θi ), was obtained through Monte Carlo studies. Note that scattered photons are accounted for
as well as direct photons by qiscat .
The expected p.e. distribution for a muon is calculated analytically. The sin2 θ dependence arises from the Cherenkov angle dependence of the intensity of the Cherenkov
photons. The term r(sin θ + r(dθ/dx)) takes into account the shrinking size of the83

Cherenkov angle as the particle loses momentum while it travels through the water.
Likelihood Calculation
The showering (e-like) and non-showering (µ-like) likelihood functions for the n-th
ring are defined as:

Ln (e, µ) =Y

exp
prob qiobs , qi,n(e, µ) +X

n′ 6=nθi <(1.5×θC )

exp
qi,n
′

!

,

(7.12)

where the product is over the PMTs inside the Cherenkov cone of the n-th ringexp
(1.5 × θC ). qiobs is the number of observed p.e.s in the i-th PMT, qi,n(e or µ) is the

expected number of p.e.s in the i-th PMT coming from the n-th ring when assuming
exp
the n-th ring was produced by either an electron or a muon. qi,n
′ is the same quantity,

only from the n′ -th ring. The function prob(qiobs , qiobs ) gives the probability of detecting
exp
exp
qiobs p.e.s in the i-th PMT given the expected amount, qiexp . The qi,n(e) and qi,n(µ)

expectation values are optimized by altering the direction and opening angle of then-th ring to yield the maximum likelihood value.
The likelihood is translated into a χ2 parameter to allow it to be combined withanother estimator that uses the Cherenkov opening angle. The χ2 value is shown
below:

χ2n (e, µ) = −2 log Ln (e, µ) + const.

(7.13)The probability for a ring to be of a particular type based on the light pattern isthen given by the following:2(χ2 (e, µ) − min[χ2n (e), χ2n (µ)])
Pnpattern (e, µ) = exp − n2σχ2 2n!

,(7.14)

84

where σχ2 2n =√2N is the approximate resolution of the χ2 distribution, N being thenumber of PMTs used in the calculation.
The probability for a ring to be of a particular type based on the Cherenkov
opening angle is given by:θobs − θnexp (e, µ)Pnangle (e, µ) = exp − n2(δθn )2

2 !

,(7.15)

where θnobs and δθn are the reconstructed opening angle of the n-th ring and the fitting
error, respectively, and θnexp (e or µ) is the expected opening angle of the n-th ringbased on the reconstructed electron-like or muon-like momentum.
The final probability calculation is simply the product of the light pattern probability calculation and the Cherenkov angle probability calculation:

P (e, µ) = P pattern (e, µ) × P angle (e, µ).

(7.16)

The final showering likelihood variable used to determine whether a ring is of ashowering or non-showering type is defined as:

Lshow ≡pp
− log P (µ) − − log P (e).

(7.17)

Positive values of Lshow indicate a preference for the non-showering type, and negativevalues indicate a preference for the showering type.

7.1.4

Ring Momentum

To calculate the momentum of the particle which produced the Cherenkov ring,the integrated charge within a 70◦ half-opening angle is used. This integrated charge
is called RT OT , and is described below:

85

RT OTn =



GMC 
α
Gdata X

θi,n <70◦−50 ns<ti <250 ns r  cos Θ
X

i
i
obs
Si 
−
qi,n exp
,L f (Θi )θ <70◦(7.18)

i,n

where:
α:normalization factor

Gdata , GMC :relative PMT gain parameter for data and MC

θi,n :

opening angle between n-th ring direction andi-th PMT directionti :residual (= time of flight subtracted) hit timeof i-th PMT

L:

light attenuation length of waterri :distance from vertex position to i-th PMT

f (Θi ) :correction function for PMT acceptance as afunction of photon incidence angle Θi

Si :expected amount of p.e.s from scattered photons for i-th PMT

The summation is restricted to a time window spanning −50 ns to +250 ns, where
0 ns represents the peak of the residual hit time distribution. The purpose of this window is to exclude light that may have originated from muon-decay electrons createdby muon decay.
After determining the RT OT value of a ring, the corresponding momenta forvarious particle hypotheses is calculated by linearly extrapolating between points on
the look-up table shown in Fig. 7.6. This typically means finding the showering-like

86

momentum (equivalent value for electron, positron, or gamma), and the muon-like
momentum. For the purposes of this study, a third particle type was added to thelist: the charged kaon. To generate the points in the look-up table, Monte Carlosimulations of single ring events for each of the particle types were generated atvarious true momentum values, and the mean RTOT value for the correspondingtrue momentum was calculated and entered into the table. The kaon points weregenerated specially for this study.

Figure 7.6: RT OT - momentum look-up table. The blue, upper set of crosses are forelectrons, the red, middle set of crosses are for muons, and the black, lower
set of crosses are for charged kaons.

The energy scale stability of the detector was tested by observing the mean reconstructed energy of stopping cosmic ray muons and the decay electrons which they
produced (see Section 3.7.3 for details). It varied within ±0.88% during the SK-I
runtime.
The absolute energy scale was adjusted separately by observing the number ofphotoelectrons generated by through-going cosmic ray muon events. A variety of

87calibration data were used to check the scale, such as stopping cosmic ray muons, thedecay electrons they produce, and the reconstructed invariant mass from π 0 particles
produced in atmospheric neutrino interactions. By comparing these data samples toMonte Carlo simulations, the absolute calibration error was estimated to be less than
±0.74% for the SK-I data taking period.
Ring Separation
In order to properly reconstruct the momenta of particles in a multi-ring event, thecharge collected by the PMTs must be correctly divided amongst the different rings.
This process is called ring separation. To perform the ring separation, a likelihood
function is used which describes the likelihood of the charge observed in a given PMT
to belong to a particular ring:L=X

θi′ ,nlog prob qi′ ,

<70◦

X
n′

αn′ · qiexp
′ ·n′

!!

,

(7.19)

where qi′ is the observed amount of charge in the i′ -th PMT, qiexp′ ,n′ is the expectedamount of charge in the i′ -th PMT coming from the n′ -th ring, prob(qiobs , qiexp ) is theprobability function for detecting qiobs charge in the i-th PMT given the total expectedamount qiexp , and αn′ is the scaling factor for each ring. The scaling factor is used as
an optimization parameter, increasing or decreasing the amount of expected charge
from each ring such that the total amount expected best matches the observed value.
The observed p.e.s in the i-th PMT belonging to the n-th ring are obtained as:exp
αn · qi,n
obs
obs
qi,n = qi P
exp ,n′ αn′ · qi,n′where the optimization parameters, αn , have been chosen to maximize L.

(7.20)

88

7.1.5

Ring Number Correction

Occasionally, a fake ring is identified and must be removed. To eliminate theserings, the ring number correction algorithm is applied. Ring number correction removes rings with a very small reconstructed momentum relative to any overlapping
rings.

7.1.6

Decay Electron Finding

The reconstruction algorithm for finding and identifying decay electrons categorizes them in three ways:• Sub-event type:The decay electrons are observed in a separate event (sub-event).
• Primary-event type:The decay electrons are observed in the primary event.
• Split type:The decay electrons occurred around the end of the event timing window, andis recorded in the primary event and sub-event.
Sub-event type decay electron events require the following criteria:(1) The time interval from the end of a primary event is less than 20 µs.
(2) The total number of hit PMTs is greater than 50.
(3) The goodness of the vertex fit is greater than 0.5.
(4) The number of hit PMTs in a 50 ns time window is greater than 30.

89(5) The total number of p.e.s is less than 2000.
(6) The number of hit PMTs in a 50 ns time window is greater than 60.
where the requirement of 60 hit PMTs in 50 ns corresponds to about 11 MeV ofelectron energy. The first criterion rejects gammas emitted from µ− captured on 16 O
nuclei.
Primary-event type decay electron events require another peak after the primary
peak with more than 20 hits above the background level in a 30 ns time window.
Additional requirements for decay electron events are:(7) the number of hit PMTs in a 30 ns time window is greater than 40 (primary-eventand split types only)(7) 0.1 µs < ∆t < 0.8 µs or 1.2 µs < ∆t < 20 µs (all types)
where ∆t is the time between the primary peak and the decay electron peak. The
final criterion rejects decays in the inefficient time interval around 1 µs. The overall
efficiency for detecting decay electrons is 80% for µ+ and 63% for µ− .

7.2

Multiple Vertex Fitter

The MVFIT (multiple vertex fitter) algorithm is simple in principle. It first starts
a loop over each ring found by the standard Super-Kamiokande event reconstruction
algorithm used in the first step of ring reconstruction described in Section 7.1. As it
considers each ring in turn, it masks all light outside of the ring’s Cherenkov angle+10◦ and subtracts any remaining light expected to have been contributed by otheroverlapping rings. Section 7.1.3 and Section 7.1.4 describe how the expected charge iscalculated for each ring. With the masking applied, selected portions of the standard

90Super-Kamiokande event reconstruction algorithm are then applied to the remaining
light pattern, treating the unmasked light pattern as if it were actually from a single
ring event. After one ring is considered, the masking is undone and the process isrepeated for the next ring until all found rings in the event have been exhausted.
MVFIT reconstructs each ring according to three particle type hypotheses andrecords the reconstructed particle information from each of the three hypotheses forevery ring in the event. These three particle hypotheses are (1) electromagnetic
showering-type particles (e± or γ), referred to in this study as gamma-like, (2) µ± ,referred to in this study as mu-like, and (3) K ± , referred to in this study as kaon-like.
A summary of the MVFIT procedure that is repeated for each ring is described
below:For each ring in an event:(1) Mask all light outside θC + 10◦(2) Subtract remaining expected light from overlapping rings(3) For each PID hypothesis (gamma-like, muon-like, and kaon-like):(a) Apply the vertex position and ring direction reconstruction algorithm described in Section 7.1.1(b) Apply the showering likelihood reconstruction algorithm described in Section 7.1.3(c) Apply the ring momentum reconstruction algorithm described in Section 7.1.4(4) Store reconstructed variables for all three PID hypotheses(5) Unmask all light in the event; continue on to next ring91

See Appendix A for the reconstructed variable resolutions obtained by MVFIT.

7.3

Reconstructed Ring Variables

By the end of the ring reconstruction process, each ring has associated with it alist of variables with values corresponding to all three of the potential PID hypotheses,
along with the PID independent showering likelihood.
The final list of reconstructed ring variables is shown below:
1. Vertex Position (kaon-like, muon-like, and gamma-like)
2. Momentum (kaon-like, muon-like, and gamma-like)
3. Cherenkov Angle (kaon-like, muon-like, and gamma-like)
4. Direction (kaon-like, muon-like, and gamma-like )5. Showering Likelihood (PID independent)Chapter 8

Dinucleon Decay Signal
Before performing the search for dinucleon decay into kaons, three main decisions
had to be made in order to limit the types of experimental signal that would bestudied to a reasonably small set. The three decisions that ultimately defined thesignal in this search were:• Choice of dinucleon decay mode
• Choice of kaon decay modes• Choice of Cherenkov ring combinations
Once the acceptable signal event types were narrowed down, a comprehensive list ofdescriptive characteristics was set down. This information then provided a platform
for the search methodology described in Chapter 9.

8.1

Choice of Dinucleon Decay ModeAt first glance, there appear to be three dinucleon decay modes with final statekaons from which to choose: p p → K + K + , p n → K + K 0 , and n n → K 0 K 0 . However, each K 0 will decay with equal probability as either a KS0 or a KL0 . This is92

93

important because the subsequent decay of the kaons is a critical part of the experimental signature.
The neutral kaons themselves emit no Cherenkov light. They can be detected
and reconstructed only through the signature of their respective decay products. It
turns out, however, that the majority of KL0 produced by dinucleon decay would often
travel outside of the detector before decaying (cτ = 15.3 m and β ∼ .8), thus being
rendered undetectable. The rate of regeneration was estimated to be too small toappreciably increase the possible detection rate of KL0 . Thus, all of the KL0 dinucleon
decay modes can be immediately eliminated from consideration.
Table 8.1 shows the branching ratios for the most probable final kaon decay states
for all combinations of K + and KS0 . The two most favorable decay modes of K + andKS0 were included in this analysis: K + → µ+ νµ (B.R. 64%); K + → π + π 0 (B.R. 21%);
KS0 → π + π − (B.R. 69%); KS0 → π 0 π 0 (B.R. 31%). A 50% penalty is taken forevery appearance of KS0 , since the listed value is intended to measure the branchingratio of the decay mode of the K 0 . A statistical factor of two has been folded into
the branching ratio of the final decay states where two identical nucleon decay via
different channels.
From Table 8.1 it can be seen that p p → K + K + has the greatest advantage
in final decay state branching ratios, with a maximum theoretical efficiency of 70%,
given the kaon decay modes that are considered. This is almost twice as large as
the total theoretically achievable efficiency of either of the other two dinucleon decay
modes, p n → K + K 0 and n n → K 0 K 0 , with maximum efficiencies of 42% and 36%,
respectively. Thus, p p → K + K + was chosen for study in this analysis.

94

Decay Modep p → K+ K+

Final Stateµ+ νµ µ+ νµµ+ νµ π + π 0
π+π0 π+π0

total
p n → K+ K0µ+ νµ π + π −µ+ νµ π 0 π 0π+π0 π+π−
π+π0 π0π0

total
n n → K0 K0
total

π+π− π+π−
π+π− π0π0
π0π0 π0π0

B.R.
40%
26%
4%
70%
22%
10%
7%
3%
42%
22%
12%
2%
36%

Table 8.1: Final state branching ratios for dinucleon decay into kaons. 100% represents
the sum of all possible final decay states for a given dinucleon decay mode.

8.2

Choice of Kaon Decay ModesAt the next level of detail, the decision of which kaon decay modes to consider had
to be made. The complexity of this choice was compounded by the fact that the final
decay states of the signal events are a combination of two kaon decays. Helpfully, thetheoretical branching ratios of the kaon decay modes are well known, and could be
used to inform this decision.
Only the two most favorable K + decay modes, K + → µ+ νµ and K + → π + π 0 , wereselected for this search. The next most favorable decay modes are three-body decays++(e.g., Ke3
and Kµ3), which are significantly more challenging from an experimental

perspective. The three-body decay modes would have a total Cherenkov ring countthat would be beyond even the large number expected from the two chosen decay
modes, significantly complicating and reducing the accuracy of the event reconstruction process. Further, the three-body decay modes lack the powerful discriminating

95

Decay Modep p → K+ K+

Final Stateµ+ νµ µ+ νµµ+ νµ π + π 0
π+π0 π+π0
sum of above

B.R. (Calc.)
40%
26%
4%
70%

B.R. (MC)
30%
19%
3%
52%

Table 8.2: Final state branching ratios for p p → K + K + . The second branching ratio(B.R.) column was obtained from the Monte Carlo, and thus includes inefficiencies due to hadronic interactions of the kaons.

Hadronic Interaction
Charge Exchange (K + n → K 0 p)
Σ Resonance (K + N → Σ π)
Λ Resonance (K + N → Λ π)

Probability per K +
9.5%
1.0%
0.5%

Table 8.3: Hadronic interactions leading to detection inefficiencies in the simulated dinucleon decay events.

power of the monochromatic momenta of the decay products seen in the chosen twobody decay modes. All other K + decay modes have insignificantly small branchingratios, and were thus not considered in this study.
The combined final state branching ratios for the two chosen kaon decay modes canbe seen in Table 8.2, where the first numerical column shows the calculated branching ratio, and the second numerical column shows the branching ratio as obtained
from the Monte Carlo. The latter takes into account inefficiencies due to hadronicinteractions of the kaons, which are displayed in Table 8.3. The leading source of
hadronic interaction inefficiencies was found to be charge exchange occurring as thekaon traveled through the water.

96

8.3

Cherenkov Ring Combinations

Having chosen which dinucleon decay mode and which final kaon decay states
to consider, one last choice must be made to narrow down the signal characteristics
to a set that can be reasonably studied. The basis of the final choice is strictly
experimental in nature; it is to determine which of the various possible Cherenkov
ring combinations that could arise in a signal event to consider in the analysis.

8.3.1

Ring Counting Inefficiencies

An important experimental limitation that had to be considered was that theSuper-Kamiokande ring counting software has a built in maximum limit of five foundrings. That means that an ideal event of the final state 16 O(pp) → 14 C K + K + →π + π 0 π + π 0 (two kaon rings, four gamma rings; see Section 8.4), could never have all
six rings reconstructed.
More confoundingly, the ring counting software also occasionally misses rings completely, as it was originally designed and optimized for atmospheric neutrino and single
proton decay (p → e+ π 0 ) events. Such events typically have at most two or threeeasily distinguishable rings, and have less than half the amount of visible light in thedetector as would be seen in a dinucleon decay event.
Additionally, there was about a one in four chance that a kaon will decay in flight(see Section 8.4.2), thus boosting its decay products along its direction of movement,
and therefore increasing the chance of creating significant ring overlap. Ring overlap
was another complication that taxes the abilities of the ring finding and reconstruction
software. Overlapping rings have a chance of being mistakenly identified as one single
ring with a reconstructed momentum larger than that of either of its true constituent
rings.

97

Nring
1
2
3
4
5

µ+ νµ µ+ νµ µ+ νµ π + π 0 π + π 0 π + π 0
5.6%
6.2%
2.8%
33.1%
23.5%
15.5%
44.9%
41.9%
38.9%
15.6%
24.5%
32.6%
0.8%
4.0%
10.2%

Table 8.4: Frequency of found ring count by final decay state. Each column sums to 100%.

Another factor to be considered was that about one fifth of the kaons are produced
below Cherenkov threshold, as shown in Fig. 8.3, due to the boost received as a resultof the Fermi momentum of the parent protons.
Table 8.4 shows the total ring count frequency for each of the chosen final decaystates. The combined effect of the sources of ring counting inefficiencies mentioned
above can be seen in the table. Ideally, one would expect to find 4 rings for K + K + →µ+ νµ µ+ νµ , and the maximum 5 rings for K + K + → µ+ νµ π + π 0 and K + K + →π + π 0 π + π 0 . Regardless, the ring counting software worked well enough with thedinucleon decay events to continue the search and ultimately achieve a respectable
final signal efficiency.
Though a large fraction of the chosen signal events had only two found rings,
this category was dropped from further consideration across all final decay states.
The reason for this choice depended largely upon the fact that many of these events
contained in-flight kaon decays (hence the low ring count), which in turn decrease theperformance of the reconstruction process described in Chapter 9. In addition, theamount of background to sort through would have been increased dramatically by
including two-ring events. Thus, only events with 3–5 rings were considered in thissearch.

98

Final Stateµ+ νµ µ+ νµµ+ νµ π + π 0
π+π0 π+π0

Nring
3
4
3
4
5
3
4

1 K+
60%
17%
61%
63%
35%
51%
62%

2 K+
36%
79%
9%
28%
54%
6%
15%

1 µ+
42%
14%
58%
71%
78%
–
–

2 µ+
54%
79%
–
–
–
–
–

1γ 2γ 3γ
–
–
–
–
–
–
45% 49%–22% 64%–
9% 54%–
9% 52% 39%
2% 22% 59%

4γ
–
–
–
–
–
–
17%

Table 8.5: Breakdown of true particle type frequency for rings generated in dinucleon
decay events. 100% represents all dinucleon decay Monte Carlo events of agiven final decay state.

8.3.2

True Particle Type Frequency

The next step in the process of narrowing down the accepted signal was to createa table such as Table 8.5, which documents the frequency of rings generated by agiven true particle type in an event of a given final decay mode.
The table immediately revealed a feature common to all final decay states, which
was that there was at least one true kaon ring in an large majority of the signal events.
This turned out to be quite useful, as the kaon rings became the crux of the event
reconstruction process described in Chapter 9.
Next, only the events with at least one true kaon ring were considered, and allof the possible ring combinations of these events were examined in Table 8.6. The
frequency of ring combinations are shown in the table as a percentage of all 16 O(pp) →14

C K + K + Monte Carlo events.
In rank order by frequency of appearance in the p p → K + K + Monte Carlo, thepossible ring combinations with at least one true kaon are:(1) K + K + µ+(3) K + γ γ(2) K + µ+ µ+(4) K + µ+ γ

99

3 Rings
K K + µ+ 7.4%
K + µ+ µ+ 6.8%
K+ K+ γ1.4%
K + µ+ γ
4.1%
K+ γ γ4.9%+

4 Rings
K K µ+ µ+
K + K + µ+ γ
K + µ+ γ γ
K+ K+ γ γ
K+ γ γ γ
+

+

3.5%
1.0%
3.2%
1.0%
2.2%

5 Rings
K K µ+ γ γ 0.4%
K+ K+ γ γ γ
0.2%+
K γγγγ0.6%++

Table 8.6: Table of found ring combinations broken down by true particle type. Numbers are taken as a percentage of the total number of simulated 16 O(pp) →14 C K + K + events. Note that the gammas originate from neutral pion decay(π 0 → γγ).

(5) K + K + µ+ µ+(10) K + K + γ γ(6) K + µ+ γ γ(11) K + γ γ γ γ(7) K + γ γ γ(12) K + K + µ+ γ γ(8) K + K + γ++(13) K + K + γ γ γ+(9) K K µ γ
To perform the final pruning of this list, the event reconstruction performance andthe characteristics of the atmospheric neutrino background had to be considered.
In order to be able to make a strong claim that even a single candidate event
found in the data is indeed a p p → K + K + event, strong evidence for two back-toback K + was required, though the evidence did not necessarily need to come in theform of two found K + rings. For lack of a second found K + ring in an event, a µ+ring could take its place, as the µ+ could be required to fall into a narrow momentum
window corresponding to K + decay. The µ+ ring would also be required to have areconstructed vertex point located where one would have expected the missing kaon
to have decayed, derived from the single found kaon ring’s vertex position and ring
direction.
A γ ring, however, was not considered a suitable substitute for a missed kaonring. Firstly, the momentum window of the true γ rings was too large to be a reliable100

indicator of K + decay. Secondly, the reconstructed vertex resolution for γ rings was
too large (see Section A.3) to be a reliable indicator of the expected event topology
corresponding to the back-to-back dinucleon decay kaons.

8.3.3

Final List of Ring Combinations

The final consideration weighed the gain in signal efficiency achieved by including
a particular ring combination against the estimated amount of effort required toaccurately reconstruct the event, along with the estimated background that would
accrue by including said ring combination. As a general rule for this search, the moregamma rings there were in the event, the more difficult it would be to accurately
reconstruct, and the more background one would have to contend with.
Thus, with all of the above considerations in mind, the ring combinations K + µ+γ, K + γ γ, K + γ γ γ, K + K + γ γ γ, and K + γ γ γ γ were eliminated. This left the
final list of signal event ring combinations, referred to as event categories in the event
reconstruction process, that were searched for in this analysis:• K + K + µ+

• K+ K+ γ• K + µ+ µ+• K + K + µ+ γ

• K + K + µ+ µ+• K+ K+ γ γ• K + µ+ γ γ• K + K + µ+ γ γ

8.4

Signal Characteristics

The characteristics of the p p → K + K + signal events were very distinct comparedthose of the atmospheric neutrino background events. In a dinucleon decay event,101

both of the kaons are ejected from the parent nucleus in a back-to-back configuration,
traveling a distance of up to 1.3m before coming to a stop and eventually decaying.
The kaons each would have a sufficient momentum to emit a detectable amount ofCherenkov light, and at least one of each kaon’s decay products would in turn emit
Cherenkov light of its own.
Figure 8.1 shows an illustration of an idealized p p → K + K + event, where one
kaon has decayed via K + → π + π 0 , and the other via K + → µ+ νµ . Note that theπ 0 is represented as two γ’s. This is because the π 0 decays at a rate fast enoughthat the distance it travels before it decays is unresolvable in a detector, therefore
the resultant γ’s can be considered to originate from the same vertex point as theirparent π 0 .

γ!())#*+,-.#!"#!())#*+,-.#
!"#

γ

$%&'#

$%&'#

Figure 8.1: A cartoon illustration of an ideal p p → K + K + event. Dashed lines represent
undetectable particles.

In total, there are three discernible particle vertices in a p p → K + K + event: one
shared by the kaons at the point of the dinucleon decay, and one additional vertex
at the point of each kaon decay. These vertices would be separated by distances thatare resolvable in Super-Kamiokande, and they all would lie approximately along theaxis of the back-to-back kaons.
All of this together presents an incredibly unique event geometry. When com-102

bined with kinematic information and other characteristics of the event’s constituent
rings, the result is an event signature that is quite difficult to mistake, when wellreconstructed.

8.4.1

Multiple Vertices

One of the most distinguishing aspects of the dinucleon decay signal is the factthat it contains three distinct vertices where Cherenkov light emitting particles areproduced, separated by distances resolvable by the detector (∼ 1 m). The first vertex
comes from the point where the dinucleon decay itself occurs. This vertex point isshared by both of the outgoing kaons. The second and third vertex points occur atthe decay points of the two kaons.
Figure 8.2 shows the distance between the kaon decay point and the original
dinucleon decay vertex. Also shown is the distance between the two kaon decay
points in each event. As expected, the distance between the kaon decay points hada peak near twice the distance of the first plot due to the back-to-back directionalityof the kaons, with a sharp fall off at larger distances.
3000
8000
2500

6000

2000

eventskaon decays

7000

5000
4000
3000

1500

1000

2000
500
1000
0
0

50

100

150

200

250

true distance from dincleon decay (cm)

300

0
0

100

200

300

400

500

true distance between kaon decays (cm)Figure 8.2: Left: The true distance between the kaon decay point and the dinucleon decay
point. Right: The true distance between the two dinucleon decay points within
the same event.

103

8.4.2

Kaon Rings

Another distinctive feature of the signal events was the presence of Cherenkov
rings produced by charged kaons. Figure 8.3 shows the true momentum distribution
for the K + produced in the dinucleon decay Monte Carlo. About 77% of the kaons
were produced above Cherenkov threshold. This feature is atypical not only whencompared to the atmospheric neutrino background, but also when compared to any
study previously performed at Super-Kamiokande. For this reason it required the
addition of new reconstruction tools, as described in Section 7.1.4.

8000
7000

kaons

6000
5000
4000
3000
2000
1000
0
0

100

200

300

400

500

600

700

800

900

1000

true kaon momentum (MeV/c)

Figure 8.3: The momentum distribution of the outgoing kaons produced by the dinucleon
decay reaction 16 O(pp) → 14 C K + K + in the Monte Carlo. The line represents
the Cherenkov threshold for K + , 563 MeV/c.

Unfortunately, it was difficult to distinguish kaon rings from muon rings on an
individual basis, as they are both non-showering particles. There was a further unfortunate similarity between the kaon and muon rings in the dinucleon decay signal
events that prevents their respective well-determined momenta from being used as an
identifying feature.
It happens that the kaon-like reconstructed momentum of the true µ+ rings in the

104

dinucleon decay signal fell right in the expected true K + momentum window, and viceversa for the muon-like momentum of the true K + rings. This can be understood byconsidering the RTOT-momentum look-up table used in the momentum reconstruction process, shown in Fig. 7.6. The peak kaon momentum is about 800 MeV/c, while
the peak muon momentum is 236 MeV/c. Both of these yield a very similar RTOT
value of roughly 500 p.e. Fortunately, the event topology was able to help informthe ring classification, as described in detail in Chapter 9, and the kaon rings in thesignal were yet able to be identified with very high fidelity.

8.4.3

Back-to-Back Kaons

Besides the mere presence of Cherenkov rings generated by charged kaons, theback-to-back directionality of the kaon rings is also a distinct feature of dinucleon
decay. Figure 8.4 shows the cosine of the opening angle between the true directions of
the kaons. It peaks sharply at −1, though there is a tail due to the Fermi momentumof the parent protons, as well as the effect of correlated decay (see Section 5.1).
Roughly 93% of the kaon pairs had an opening angle between them of greater than154◦ (i.e., cos(θ) < −0.9).

8.4.4

Features of Kaon Decay Modes

Each of the final decay modes of the kaons had its own set of unique features,
described in this section.
K + → µ+ νµ Features
In the muon decay mode, K + → µ+ νµ , the outgoing particles are emitted backto-back from the kaon decay point, each with a momentum of precisely 236 MeV/c.

105kaon pairs

104

3

10

102

-1

-0.9

-0.8

-0.7

-0.6

-0.5

cos θ

-0.4

-0.3

-0.2

-0.1

0

Figure 8.4: Cosine of the opening angle between the true directional vectors of the outgoing kaons.

This is above Cherenkov threshold for the muon (pthresh [µ± ] = 120 MeV/c), giving
the particle a Cherenkov angle of 34.6◦ . The neutrino leaves the tank without being
detected.
Figure 8.5 shows the distribution of true muon momenta produced by the decay
mode K + → µ+ νµ in the dinucleon decay Monte Carlo. The peak at 236 MeV/ccorresponds to the monochromatic momentum resulting from a kaon decaying whileat rest. The momentum values that differ from 236 MeV/c occurred when the parent
kaon decayed in flight before it could come to a stop through ionization energy loss.
These off-peak muons accounted for about 27% of all muons in the dinucleon decayMonte Carlo that were produced through kaon decay.
K + → π + π 0 Features
In the pion decay mode, K + → π + π 0 , the outgoing particles are again emitted
back-to-back, this time with a precise momentum of 203 MeV/c. However, the π 0
immediately (τ = 8.4 × 10−17 s) decays into two gammas, each of which is boosted

106

muons

104

3

10

102

10
0

100

200

300

400

500

600

700

800

true muon momentum (MeV/c)

Figure 8.5: True momentum distribution of muons produced by the kaon decay channel K + → µ+ νµ in the dinucleon decay Monte Carlo. The line represents
the Cherenkov threshold for µ+ , 120 MeV/c. The peak corresponds to themonochromatic momentum of 236 MeV/c resulting from a kaon decaying whileat rest. The off-peak values are from in-flight decays of the parent kaons.

along the direction of the parent pion. This gives the gammas momenta ranging from20 MeV/c to 227 MeV/c, as shown in Fig. 8.6, with a small tail resulting from in-flightkaon decays.
The π + , on the other hand, have a precisely determined momentum of 203 MeV/c,
which is barely above Cherenkov threshold (pthresh [π ± ] = 160 MeV/c). The small
number of Cherenkov photons produced by the π + , however, get drowned out by the
copious amount of light produced by all of the other visible particles in the dinucleon
decay event. The π + is effectively undetectable, and therefore no attempt was madeto search for any rings produced by π + in this analysis.

107

2500

gammas

2000

1500

1000

500

0
0

100

200

300

400

500

600

700

800

true gamma momentum (MeV/c)

Figure 8.6: True momentum distribution of gammas produced by the kaon decay channel
K + → π + π 0 in the dinucleon decay Monte Carlo. The range of momenta
produced by decays while the kaon is at rest is 20–227 MeV/c. The tail is fromin-flight decays of the parent kaons. Gammas are EM showering particles,
thus will always yield detectable Cherenkov light for any energy above a few
MeV/c.

Chapter 9

Dinucleon Decay Search
The underlying methodology of the search was to hypothesize that each event
analyzed was generated by dinucleon decay into kaons, and then compare how wellthe event’s reconstructed variables matched the expected values for a genuine signal
event, as determined by the dinucleon decay Monte Carlo.
The search was performed in four steps. The first (and most involved) step was
to classify all found rings in an event as either kaon, muon, or gamma candidaterings. The second step was to categorize the event based on the classification of
its constituent rings. In the third step, several precuts were applied which reducethe background and “cleaned up” the signal without reducing the search efficiency
too much. Finally, in the fourth and final step, thirty-seven reconstructed variables
were all combined into a single discriminatory variable through the use of a boosted
decision tree. A final cut was placed on the output of the boosted decision tree which
was used to define the final accepted signal-like sample.

108

109

9.1

Ring Classification

The ring classification stage depended upon the hypothesis that the observed
rings in an event arose from the reaction 16 O(pp) → 14 C K + K + , where the kaonsthen decayed via K + → µ+ νµ or K + → π + π 0 . Under this assumption, the topologyof the event may inform the choice of ring classification in a powerful way.
The term “ring classification” is used in this analysis, whereas “particle identification” (or “PID”) is not. This was a conscious choice to reflect that rings were being
matched up with the most appropriate particle type based on the hypothesis thatthe event is indeed a dinucleon decay event. This should yield correct results for thesignal events, however the particle type assignments should be incorrect almost by
design when considering background atmospheric neutrino events. Hence the term
“PID” was deemed somewhat inappropriate to use in this stage of the search.
In more conventional Super-Kamiokande analyses, PID is performed at the individual ring level, essentially only taking into account whether or not the particle
is of a showering or non-showering type. The former would indicate an electron,
positron, or gamma, and the latter would indicate a muon or a heavier charged particle. This level of distinction is nearly always sufficient for the studies performed at
Super-Kamiokande.
In the case of dinucleon decay into kaons, however, rings generated by charged
kaons would be identified as non-showering, but no further disambiguation could bemade between said kaon rings and rings generated by muons based on their respective
showering likelihoods alone. This was the motivating factor behind incorporating
topology into the ring classification process.
The next important assumption regarded which of the Cherenkov light producing
particles in a given dinucleon decay event were most likely to be found. Table 9.1, Ta-110

ble 9.2, and Table 9.3 show the frequency of finding a given number of rings producedby kaons, muons, and gammas in the dinucleon decay Monte Carlo, respectively. The
entries in the tables are organized by final state and the total number of rings foundin an event. Note that the γ’s were the by-products of the π 0 produced by the decayof the K + .
Final Stateµ+ νµ µ+ νµµ+ νµ π + π 0
π+π0 π+π0

Nring
3
4
3
4
5
3
4

0 K+ 1 K+ 2 K+ ≥ 1 K+
4%
60%
36%
96%
4%
17%
79%
96%
30%
61%
9%
70%
9%
63%
28%
91%
11%
35%
54%
89%
44%
51%
6%
56%
23%
62%
15%
77%

Table 9.1: Frequency of rings produced by kaons in p p → K + K + Monte Carlo events.

Final Stateµ+ νµ µ+ νµµ+ νµ π + π 0

Nring
3
4
3
4
5

0 µ+ 1 µ+ 2 µ+ ≥ 1 µ+
4% 42% 54%
96%
7% 14% 79%
93%
42% 58%–58%
29% 71%–
71%
22% 78%–78%

Table 9.2: Frequency of rings produced by muons in p p → K + K + Monte Carlo events.

Final Stateµ+ νµ π + π 0
π+π0 π+π0

Nring
3
4
5
3
4

0γ
6%
14%
37%
0%
0%

1γ
45%
22%
9%
9%
2%

2γ
49%
64%
54%
52%
22%

3γ 4γ ≥1γ
–
–
94%
–
–
86%
–
–
63%
39%
– 100%
59% 17% 100%

Table 9.3: Frequency of rings produced by gammas in p p → K + K + Monte Carlo events.

111

An important feature that can be seen in Table 9.1 is that there was at least one
ring produced by a kaon in the majority of all dinucleon decay events of the threechosen final states. The ring classification algorithm took advantage of this fact andused the kaon rings as the cornerstone of the full classification process.

9.1.1

Two Kaon Ring Classification

The key to ring classification and the overall event reconstruction in this search
was to first identify the kaon candidate ring(s) in the event. Accordingly, the firstpriority of the ring classification process was to determine whether or not there are
two rings that can be classified as kaon candidates. This was done by considering thepair of rings with the largest opening angle between their directions and checking to
see whether they satisfied the following criteria:(1) the opening angle must be > 154◦
(2) the vertex separation must be < 640 cm(3) the showering likelihood of both rings must be > −2.4(4) the Cherenkov angle of both rings must be < 40◦(5) the invariant mass must be < 1940 MeV/c2(6) the total momentum must be < 400 MeV/c
where the last three items were all calculated using the kaon-like values for both rings.
Figure 9.1 shows the distributions of the variables used in the two-kaon classification process for the dinucleon decay Monte Carlo. Each entry represents either thepair of rings with the largest opening angle between their directions in a given event

112(e.g., the reconstructed invariant mass), or one of those two rings counted individually(e.g., the showering likelihood).

9.1.2

Kaon Decay Product Ring Classification

After the two kaon candidates have been identified, the next step was to classify
all of the remaining rings in the event as one of the two possible kaon decay productring types, muon or gamma. This step of the classification relied heavily upon thegeometric relationship between the kaons and their decay products. See Fig. 8.1 foran illustration of the dinucleon decay event geometry.
The primary instrument in this step of the classification was the impact parameterfor each particle type hypothesis. The impact parameter was defined as the distance
between the vertex of a kaon decay product candidate ring and an infinite line drawn
along the directional vector of the most likely parent kaon candidate. See Fig. 9.2 for
an illustration of the impact parameter.
The particle type hypothesis which yielded the smallest impact parameter became
the classification of the ring. Figure 9.3 shows the impact parameters for all true µ+and γ rings in the dinucleon decay Monte Carlo as calculated using both the muonlike and gamma-like vertex of every ring. The difference between the incorrect andcorrect impact parameter can be seen in Fig. 9.4. The positive entries in this plot
represent rings where the correct classification would be assigned using the impact
parameter method.
To further refine the performance of the kaon decay product ring classification,
all muon and gamma candidate rings were also required to satisfy the type-specificcriteria listed below:• Muon candidates: 200 MeV/c < muon-like momentum < 300 MeV/c

113

+

+

All K K MC

1400

1000

1200

+

True K Pairs
ring pairs

ring pairs

800

Other Pairs
600

400

1000
800
600
400

200
200
0
0

20

40

60

80

100

120

140

160

0
0

180

200

400

opening angle between ring directions (deg)

600

800

1000

1200 1400 1600

1800 2000

vertex separation (cm)

5000

3500
3000

4000

rings

rings

2500
3000

2000
1500

2000

1000
1000
500
0
-8

-6

-4

-2

0

2

4

6

0
10

8

20

showering likelihood

30

40

50

60

Cherenkov angle (deg)

3000
1200
2500

ring pairs

ring pairs

1000
2000
1500
1000

600
400

500
0
1000

800

200

1200

1400

1600

1800

2000

invariant mass (MeV/c2)

2200

2400

0
0

200

400

600

800

1000

1200 1400 1600

1800 2000

total momentum (MeV/c)

Figure 9.1: Reconstructed variables used in the two-kaon classification process. Kaonlike variables were used in all applicable plots. The dashed lines show thesubset of events where the rings were actually generated by two kaons, andthe dotted lines show the subset of events where the rings were generated byother particles (including events where only one of the rings from the pair wasa true kaon ring). Vertical lines indicate the two-kaon classification criteria.

114

!"#
Figure 9.2: Illustration of the impact parameter variable used in the classification of thekaon decay product candidate rings. The illustration shows part of a reconstructed event where one kaon decayed via K + → µ+ νµ . The circles represent
reconstructed vertices and the arrows represent reconstructed ring directions.
The correct classification has been chosen for the true K + ring, and the impact parameter is shown for both possible classifications of the true µ+ ring.
The correct classification yields a smaller impact parameter.

3000

2200

Muon-Like Vertex

2000

Gamma-Like Vertex

Muon-Like Vertex
Gamma-Like Vertex

2500

1800
1600

2000

rings

rings

1400
1200
1000
800

1500

1000

600
400

500

200
0
0

100

200

300

400

500

impact of true µ + rings (cm)

600

700

0
0

100

200

300

400

500

impact of true γ rings (cm)

600

700

Figure 9.3: Impact parameter calculated using both the muon-like and gamma-like vertices of all true muon (left) and gamma (right) rings from the dinucleon decayMonte Carlo. Solid lines indicate the muon-like reconstructed vertex wasused in the calculation of the impact parameter, and dashed lines indicate the
gamma-like reconstructed vertex was used.

115

True µ + Rings

1400

True γ Rings

1200

rings

1000
800
600
400
200
0
-800

-600

-400

-200

0

200

400

600

800

incorrect impact - correct impact (cm)Figure 9.4: The difference between the impact parameter as calculated using the incorrect
and the correct reconstructed vertex for all true muon and gamma rings.
A positive value indicates that the correct particle type yielded the smaller
impact parameter, thus would be chosen by the ring classification algorithm.
The dashed line represents the true muon rings, and the solid line represents
the true gamma rings from the dinucleon decay Monte Carlo.

• Gamma candidates: gamma-like Cherenkov angle > 30◦
Figure 9.5 shows the distributions of the variables used in the criteria above.
The overall performance evaluation for kaon decay product ring classification was
estimated by assigning the correct classification to all kaon rings in every event (i.e.,
the kaons are assumed to have been properly classified), and then applying the impact
parameter method plus the secondary particle-specific criteria. The approximate misclassification rates under the given assumptions were about 20% for true muon ringsand 28% for true gamma rings. The reason γ rings were more often mis-classified isbecause the reconstructed vertex resolution was somewhat worse than that of muonrings (see Appendix A).

9.1.3

Pion Classification

If the kaon decay product classification process classified two rings from the same
parent kaon as gamma candidates, they were then taken to be from the same original

116

All K +K + MC
12000

True γ Rings

7000

Other Rings

10000

Other Rings
6000

8000

rings

rings

All K +K + MC

8000

True µ + Rings

6000

5000
4000
3000

4000
2000
2000

1000

0
100

200

300

400

500

600

muon-like momentum (MeV/c)

700

0
0

10

20

30

40

50

60

70

gamma-like Cherenkov angle (deg)Figure 9.5: Reconstructed variables used in the secondary criteria of the kaon decay product classification process. Left: Reconstructed muon-like momentum of allrings in the dinucleon decay Monte Carlo. Right: Reconstructed gamma-like
Cherenkov angle of all rings in the dinucleon decay Monte Carlo. In bothplots, the dashed line represents the true muon (left) or gamma (right) rings,
and the dotted line represents all other rings. The vertical lines indicate theclassification criteria.

π 0 . The π 0 was reconstructed in full and considered to be the kaon decay product. Itsvertex was taken to be the spatial average of the reconstructed gamma-like verticesof its two constituent gamma rings.

9.1.4

Low Gamma Momentum Correction

A final correctional step was taken to catch true muon rings that have been misclassified as gamma candidates. This was done when the total momentum of all
gamma candidate rings was less than 100 MeV/c. If there was only one gamma candidate ring, it was reassigned as a muon candidate, as long as there were not alreadytwo muon candidates in the event. If there were two gamma candidates, the most
non-showering ring was reassigned as a muon candidate first, and then the second.
Again, the rule was held where there can be no more than two muon candidates in asingle event. If there existed two muon candidates when the algorithm was attempting to reassign a current gamma candidate, the ring simply remained as a gamma
candidate.

117

9.1.5

Single Kaon Ring Classification

If an event failed to meet the two-kaon search criteria outlined in Section 9.1.1,
the single-kaon search criteria was then applied. About 23% of the kaons produced
in dinucleon decay events had momenta below Cherenkov threshold and about 27%
decay in flight, reducing the performance of the reconstruction software. For thesereasons, events with only a single found kaon ring represented a significant enoughportion of the signal that it entailed the creation of an additional algorithm designed
to classify the rings in such events.
To find the most appropriate kaon candidate in the event, each ring was consideredin turn as the “test” kaon candidate, as long as its kaon-like Cherenkov angle was< 40◦ (see Fig. 9.6). The other rings in the event were forced to conform to the mostappropriate kaon decay product test classifications under the assumption that thetest kaon candidate ring was a true kaon ring. The test decay product classifications
were determined by the algorithm set forth in Section 9.1.2. A variable called the“linearity” was then calculated for the test kaon candidate, effectively representing
the likelihood that a given ring was from a true dinucleon decay kaon.
7000

All K +K + MC

6000

True K+ Rings
Other Rings

rings

5000
4000
3000
2000
1000
0
10

20

30

40

50

60

70

kaon-like Cherenkov anlge (deg)Figure 9.6: Reconstructed kaon-like Cherenkov angles for all rings in the dinucleon decayMonte Carlo. The dashed line represents the true kaon rings, and the dottedline represents all other rings. The vertical line represents the single-kaonclassification criteria.

118

The linearity quantifies the linear geometry of the event, given the current test
ring classifications. The linearity was calculated by first taking the cosine of the angle
formed between a line drawn from the test kaon candidate’s kaon-like vertex point to
one of the test kaon decay product candidates’ appropriate reconstructed vertex point
and the directional vector of the test kaon candidate. This value was then multipliedby the cosine of the angle formed by the test kaon candidate’s direction vector andthe line between its vertex and the other decay product candidate ring’s appropriate
vertex. Figure 9.7 shows an illustration of how the linearity was determined. The
more well aligned all of the vertices were along the test kaon candidate’s direction,
the more negative this value would be. A perfectly linear configuration would yielda linearity value of −1.

θ$
θ%
µ+

!"#

Figure 9.7: Illustration of the linearity variable used in the single-kaon classification process. The illustration shows reconstructed K + K + → µ+ νµ µ+ νµ event where
the correct classification has been chosen for the three found rings (K + , µ+ ,µ+ ). The second kaon is here assumed to have been produced below Cherenkov
threshold. The circles represent the reconstructed vertices, and the arrows
represent the reconstructed directions. The linearity for the visible kaon candidate ring is defined as cos θ1 · cos θ2 .

Figure 9.8 shows the difference between the linearity calculated for every true
gamma and muon ring and the true kaon ring in the same event, for all events in the
dinucleon decay Monte Carlo which had only one true visible kaon ring. Entries withpositive values represent events where the true kaon ring would have been chosenover the other rings in the event as the kaon candidate. Fifty-seven percent of the119dinucleon decay Monte Carlo events which contained only one true kaon event werecorrectly classified as a single kaon event. Of those, 94% correctly identified the truekaon ring as the kaon candidate.

ring pairs+70µ + - K linearity60

γ - K linearity

+

50
40
30
20
10
0
-2

-1.5

-1

-0.5

0

0.5

1

1.5

2

true decay prod. linearity - true kaon linearity

Figure 9.8: Difference in linearity between true muon and gamma rings and true kaonrings. Positive values indicate the true kaon ring was favored over the otherring as the kaon candidate.

9.2

Event Categorization

The event categorization process is simply the process of labeling each event according to the classifications of its constituent rings. For example, an event in which
two rings were classified as kaon candidates, and one ring was classified as a muoncandidate would be categorized as a “K + K + µ+ ” event.
If at least one ring was not able to be classified as either a kaon, muon, or gamma
candidate by the classification algorithms described in Section 9.1, then the event is
simply categorized as an “other”-type event.

120

9.3

Precuts

After the ring classification and event categorization processes, several precuts
were applied with the goal of eliminating easily identifiable background and restricting
signal sample to events which were well reconstructed while retaining as large an
efficiency as possible. The precuts are listed below:(1) 1000 < P OT OT < 11000, where P OT OT is the total amount of light seen in the
ID in units of p.e.
(2) 3 ≤ NRING ≤ 5, where NRING is the number of found rings in the event.
(3) 0 ≤ NDCY E ≤ 2, where NDCY E is the number of found decay electrons.
(4) W ALL > 200 cm, where W ALL is the distance from the event vertex to thenearest wall; this is the fiducial volume (FV) cut.
(5) And finally, the event category must fall into one of the eight accepted dinucleon
decay event categories listed in Section 8.3.3.
Figure 9.9 shows the distributions of the reconstructed variables used in the precuts.
The first three precuts were intended to retain roughly all of the signal while
significantly reducing the background, which had a much wider spread in P OT OT
than the signal, and an event count distribution that fell rapidly with increasing ringcount. The decay electron count was limited to the maximum number we would
expect from the dinucleon decay signal events.
The fiducial volume cut was applied to ensure optimal and uniform reconstruction
resolution throughout the target region. Further, it created a well-defined boundarywith a volume that was simple to calculate. Both of these factors were crucial in the

121

10

SK1 Data

SK1 Data
25

Atm. ν MC
pp → K+K + MC

6

pp → K+K + MC20

eventsevents8

Atm. ν MC

4

15

10

2

5

0
1000 2000 3000

4000 5000 6000

07000 8000 9000 10000 11000

3

4

total charge in the ID (p.e.)5number of found rings

SK1 Data
600
500

SK1 Data

Atm. ν MC

14

Atm. ν MC

pp → K+K + MC

12pp → K+K + MCeventsevents

10
400
300

8
6

200

4

100
0

2
0

1

2

3

0
0

200

NDCYE (number of decay electrons)

400

600

800

1000

1200

1400

1600

WALL (cm)Figure 9.9: Reconstructed variables used in the precuts. From top-left to bottom-right:
P OT OT , N RIN G, N M U EDCY , and W ALL. The solid lines indicate thep p → K + K + Monte Carlo, the dashed lines indicate the atmospheric neutrinoMonte Carlo, and the markers indicate the SK-I data. Vertical lines indicate
the boundaries of the precuts.

estimation of the total exposure. The fiducial volume was defined as the region in
the ID that is at least 200 cm from the nearest wall.
Finally, the event category requirement ensures that the event was well reconstructed, and limits the search to categories that are most likely to come from a true
dinucleon decay event. See Section 8.3 for a more detailed explanation of why theseparticular event categories were chosen.
Table 9.4 shows signal efficiency, expected background, and data events passing
the precuts, broken down by event category. The total efficiency for the signal after
applying the precuts was 21.9%. The total expected background was 33.9 events after
normalizing to the SK-I livetime. The total number of events in the SK-I data passing

122

the precuts was 27.
Event Category
K + K + µ+
K + µ+ µ+
K+ K+ γK + K + µ+ µ+
K+ K+ γ γ
K + µ+ γ γ
K + K + µ+ γ
K + K + µ+ γ γ
total

p p → K + K + M.C.
Atmospheric ν M.C.
efficiency (%) events / SK-I livetime
6.1
0.6
5.6
6.1
1.1
0.4
2.7
0.0
0.6
0.3
4.6
26.1
1.1
0.3
0.2
0.1
21.9
33.9

SK-I Data
events
0
5
0
0
0
22
0
0
27

Table 9.4: Signal, background, and data after applying precuts.

The event category in which the SK-I data and the expected background rate differed the most was K + µ+ γ γ, where the atmospheric neutrino Monte Carlo predicted
26 events of background, and only 22 such events were found in the data. Using Poisson statistics, the probability of seeing 22 or fewer events given the expected value of
26 is about 20%.
Table 9.5 summarizes background Monte Carlo events which survived the precutsby neutrino flavor and weak current type, and Table 9.6 shows the breakdown ofneutrino interaction types of the surviving background. Charged current (CC) andneutral current (NC) muon neutrino events made up the majority of the backgroundafter the precuts. Single and multi-pion interactions were the dominant interaction
types.

9.4

Multi-Variate Analysis

For the final portion of the analysis, a multi-variate approach was taken to increase
the sensitivity of the search. Multi-variate techniques were considered appealing for123ν Mode
CC νµ
NC νµ
CC νe
NC νe
total

Events/200 Years Events/1489.2 Days%
1654
13.5 411303
10.6 32588
4.8 14517
4.2 134062
33.1 100

Table 9.5: Remaining background Monte Carlo events after precuts sorted by neutrino
flavor and weak current type. CC: Charged Current; NC: Neutral Current.
The first numerical column shows the unscaled number of events from the 200
year analysis sample. The second numerical column shows the number of events
scaled to the SK-I livetime.

Interaction
CC single pion delta resonanceCC multi-pion production
NC multi-pion production
NC diffractive pion production
NC single pion delta resonanceCC diffractive pion production
CC quasi-elastic
NC elastic
CC coherent pion production
NC coherent pion production
total

Events/200 Years Events/1489.2 Days
942
7.7
896
7.3
660
5.4
549
4.5
542
4.4
265
2.2
104
0.9
49
0.4
35
0.3
20
0.2
4062
33.1

%
23.2
22.0
16.2
13.5
13.3
6.5
2.6
1.2
0.9
0.5
100

Table 9.6: Remaining background Monte Carlo events after precuts sorted by neutrino
interaction type. CC: Charged Current; NC: Neutral Current. The first numerical column shows the unscaled number of events from the 200 year analysis
sample. The second numerical column shows the number of events scaled to
the SK-I livetime.

this study due to their ability to simultaneously analyze many different discriminatory
variables spanning many different event categories. This would be extremely difficult
to accomplish using a simple cuts-based analysis, such as the one used in the p → e+ π 0
search, which hinges largely upon the combined discriminating power of only twovariables: the invariant mass and the total momentum of the system.

124

After some consideration, a boosted decision tree was chosen to be used as themulti-variate tool for this analysis. Several factors played in this decision. First,
the performance of a boosted decision tree is not degraded by the addition of weakly
discriminating variables [74]. If a variable is so weakly discriminating as to be useless,
it is simply not used by the boosted decision tree.
Second, a boosted decision tree will be unhindered by variables that may havetrivial values for some events. This is accomplished through the use of a trivial bin
designated by the user. For example, an event which has only one kaon candidate ringmight be assigned a meaningless value of −10 for its reconstructed dinucleon mass
variable during the classification and categorization stage of the analysis. Thus if theboosted decision tree were to place a cut at 0 on this variable, it would effectively
be the same as cutting on whether or not an event had two kaon candidates. This
technique requires that the trivial bin be outside of the region of interest for the nontrivial entries, though it shouldn’t be too far from the region of interest, or else thebinning resolution of the variable will be degraded. In a multi-variate analysis, every
event has to have some value for every variable, so the technique of using a trivial
bin is particularly useful in this analysis because it allows events of all categories tobe processed by the same boosted decision tree.
Finally, boosted decision trees have a proven track record in background dominated studies in the world of particle physics [75–77].
The process of creating and applying the boosted decision tree was accomplished inthree stages: training, testing, and analysis. Each stage required its own independent
set of signal and background Monte Carlo. As such, the atmospheric neutrino anddinucleon decay Monte Carlo sets were both divided into three roughly equal portions
as described in Table 9.7. The actual SK-I data was only ever used in the final stageof performing the analysis, thus required no splitting.

125

Monte Carlop p → K+ K+
Atm. ν

Training
25,000 events
150 years

Testing
25,000 events
150 years

Analysis
25,000 events
200 years

Table 9.7: Division of the Monte Carlo in the multi-variate portion of the analysis.

9.4.1

Input Variables

In total, 37 variables were chosen to be used as inputs for the boosted decision
tree. The names of these variables along with a brief description of each follow:

Variables related to the reconstructed dinucleon• rec dinuc mass: The reconstructed dinucleon pair mass. Expected to peaknear the sum of the true mass of two nucleons, roughly 1900 MeV/c2 , for signal
events, and random for background events, which never contain two true kaons
produced by dinucleon decay.
• rec dinuc mom: The reconstructed dinucleon pair momentum. Expected tobe small (around 200 MeV/c) for signal events, and be random for background
events, which never contain two true kaons produced by dinucleon decay.
• k cand vtxsep: The separation of the vertices of the kaon candidates whichmake up the reconstructed dinucleon. Expected to be small for signal events,
and random for background events, which never contain two true kaons produced by dinucleon decay.

Variables related to the reconstructed π 0

126• rec pi0 mass: The reconstructed π 0 mass. Expected to peak at the true mass,
135 MeV/c2 , for signal events, and be somewhat random for background events,
which may contain true pions, but are reconstructed under a false hypothesis.
• rec pi0 mom: The reconstructed π 0 momentum. Expected to peak around
206 MeV/c for signal events, and have a random distribution for background
events, which do not contain monochromatic pions produced by kaon decay.
• g cand vtxsep: The separation of the vertices of the gamma candidates whichmake up the reconstructed π 0 . Expected to be small for signal events, andsomewhat random for background events, which are reconstructed under a false
hypothesis.
• rec pi0 impact: The impact parameter of the reconstructed π 0 . Expected tobe small for signal events and larger for background events, which are reconstructed under a false hypothesis.
• g cand totmom: The total momentum of all gamma candidates. Expected
to peak around 206 MeV/c for signal events, where it is assumed that bothtrue gamma rings are reconstructed as a single ring in events with only one
gamma candidate. Background events, which do not contain monochromaticpions produced by kaon decay, are expected to have a random distribution.

Variables related to the reconstructed momentum• k cand mom[0,1]: The reconstructed momentum of the kaon candidates. Expected to peak roughly around 800 MeV/c for signal events, and have a random
distribution for background events, which almost never even contain true kaons.

127• mu cand mom[0,1]: The reconstructed momentum of the muon candidates.
Expected to peak at 236 MeV/c for signal events, and have a random distribution
for background events, which do not contain monochromatic muons producedby kaon decay.
• g cand mom[0,1]: The reconstructed momentum of the gamma candidates.
Expected to have a distribution consistent with a monochromatic π 0 with amomentum of 207 MeV/c (about 20–230 MeV/c) for signal events, and have arandom distribution for background events, which do not contain monochromatic pions produced by kaon decay.

Variables related to the reconstructed Cherenkov angle• k cand cang[0,1]: The reconstructed Cherenkov angle of the kaon candidates.
Expected to be roughly around 30◦ for signal events, and have a random distribution for background events, which almost never even contain true kaons, letalone kaons produced by dinucleon decay.
• mu cand cang[0,1]: The reconstructed Cherenkov angle of the muon candidates. Expected to be roughly around 35◦ for signal events, and have a random
distribution for background events, which do not contain monochromatic muons
produced by kaon decay.
• g cand cang[0,1]: The reconstructed Cherenkov angle of the gamma candidates. Expected to be roughly around 42◦ for signal events, and have a somewhat more random distribution for background events, which may contain pions,but are reconstructed under a false hypothesis.

128

Variables related to the showering likelihood• k cand show[0,1]: The showering likelihood of the kaon candidates. Expected
to have a clear peak in the positive (non-showering) region for signal events,and have a somewhat more random distribution for background events, which
are reconstructed under a false hypothesis.
• mu cand show[0,1]: The showering likelihood of the muon candidates. Expected to have a clear peak in the positive (non-showering) region for signal
events, and have a somewhat more random distribution for background events,
which are reconstructed under a false hypothesis.
• g cand show[0,1]: The showering likelihood of the gamma candidates. Expected to have a clear peak in the negative (showering) region for signal events,and have a somewhat more random distribution for background events, which
are reconstructed under a false hypothesis.

Variables related to the event geometry• mu cand impact[0,1]: The impact parameter of the muon candidates. Expected to be small for signal events and larger for background events, which arereconstructed under a false hypothesis.
• g cand impact[0,1]: The impact parameter of the gamma candidates. Expected to be small for signal events and larger for background events, which arereconstructed under a false hypothesis.
• decay vtxsep: The separation between the reconstructed vertices of the kaon
decay product vertices. Expected to be large (around 2.6 m) for signal events,

129

and smaller in background events, which have one true vertex. This variable isfairly accurate for background events despite their being reconstructed undera false hypothesis. Their accuracy is due to the fact that the decay productcandidates are more likely to be correctly labeled as gammas and muons, thus
yielding accurate reconstructed information.

Variables related to the light pattern and the single vertex hypothesis• total rtot mag: The magnitude of the vector sum of the light contained ineach ring from the event vertex (RT OT ). Expected to be small in signal events
due to the symmetry of the light pattern and larger in background events, which
tend to be less symmetric.
• total rtot dot k: The magnitude of the vector sum of the light containedin each ring from the event vertex (called “rtot”) along the direction of thebrightest kaon candidate ring. Expected to be small in signal events due tothe back-to-back kaon rings of very similar momenta, and larger in background
events, which lack this feature.
• agood: The goodness of the standard vertex fitter, which assumes the event has
only one true vertex. Expected to be small for signal events, which have threetrue vertices, and larger for background events, which have one true vertex.
• off vtx sep: The separation between the standard (“official”) reconstructed
event vertex and the dinucleon decay hypothesis event vertex. Expected to besmall for signal events due to the symmetry of the underlying event geometry,and larger in background events, which are less symmetric and reconstructedunder a false hypothesis.

130

Variables related to the decay electrons• n dcy e: The number of found decay electrons. Expected to be larger for signal
events, which produce two true decay electrons, and smaller for background
events, which tend to produce fewer true decay electrons.
• dcy e vtxsep: The separation between the vertices of the found decay electrons. Expected to be large (around 2.6 m) for signal events, and small (∼0 m)
for background events, which do not have true separated pion or muon vertices.
Figure 9.10 shows the input variable distributions for the analysis sample of the
p p → K + K + signal Monte Carlo, the analysis sample of the atmospheric neutrino
background Monte Carlo, and the SK-I data. The data and atmospheric neutrinoMonte Carlo are in good agreement overall.

131

SK1 Data

0.8

SK1 Data

1.4

Atm. ν MC

0.7

Atm. ν MC

1.2

pp → K+K + MC

0.6

pp → K+K + MC1eventsevents

1.6

0.8

0.5
0.4

0.6

0.3

0.4

0.2

0.2

0.1

0
1600

1650

1700

1750

1800

1850

1900

1950

reconstructed dinucleon mass (MeV/c2)

0
0

2000

50

100

150

200

250

300

350

400

450

500

reconstructed dinucleon momentum (MeV/c)

SK1 Data
1

Atm. ν MC
pp → K+K + MCevents

0.8

0.6

0.4

0.2

0
0

100

200

300

400

+

500

600

700

vertex separation of K candidates (cm)Figure 9.10: Distributions of the boosted decision tree input variables used to perform thefinal analysis. Trivial bins have not been plotted. The solid line represents
the p p → K + K + Monte Carlo, the dashed line represents the atmospheric
neutrino Monte Carlo, and the crosses represent the SK-I data. This set ofplots shows variables related to the reconstructed dinucleon pair. No data
events were found to have two kaon candidates, thus there are no data entries
in this set of plots.

132

SK1 Data

SK1 Data

7

Atm. ν MC

7

Atm. ν MC

6

pp → K+K + MC

6

pp → K+K + MC5eventsevents

5
4

4

3

3

2

2

1

1

0
0

0
0

50

100

150

200

250

300

350

400

reconstructed π0 mass (MeV/c 2)

450

500

100

200

300

400

500

600

700

800

reconstructed π0 momentum (MeV/c)

900

1000

8
SK1 Data

SK1 Data

7

Atm. ν MC

6

pp → K+K + MC

10Atm. ν MC
pp → K+K + MC

8eventsevents

5
4
3

6

4

2
2
1
0
0

100

200

300

400

500

600

vertex separation of γ candidates (cm)

0
0

700

100

200

300

400

reconstructed π0 impact parameter (cm)500

SK1 Data
6

Atm. ν MC
pp → K+K + MCevents

5
4
3
2
1
0
0

100

200

300

400

500

600

700

800

total momentum of γ candidates (MeV/c)

900

1000

Figure 9.10: (Continued...) This set of plots shows variables related to the reconstructed
π0 .

133

SK1 Data
7

Atm. ν MC

pp → K+K + MC

6

pp → K+K + MC

µ + candidate rings

8

SK1 Data

Atm. ν MC

6

4+

K candidate rings

10

5
4
3
2

2
1
0
600

700

800

900

1000

1100

1200

1300

1400

0

1500

160

180

200

momentum (MeV/c)

220

240

260

280

300

320

340

momentum (MeV/c)

SK1 Data
14

Atm. ν MC

γ candidate rings

12

pp → K+K + MC

10
8
6
4
2
0
0

100

200

300

400

500

600

700

800

900

momentum (MeV/c)

Figure 9.10: (Continued...) This set of plots shows the reconstructed momenta.

134

14

SK1 Data

12

Atm. ν MC6µ + candidate rings

pp → K+K + MC
10
8
6

+

K candidate rings

SK1 Data

7

Atm. ν MC

4

4
3
2

2

1

0
10

0
10

20

30

40

50

60

pp → K+K + MC

520

Cherenkov angle (deg)

30

40

50

60

Cherenkov angle (deg)SK1 Data

14

Atm. ν MC

γ candidate rings

12

pp → K+K + MC

10
8
6
4
2
0
10

20

30

40

50

60

Cherenkov angle (deg)Figure 9.10: (Continued...) This set of plots shows the reconstructed Cherenkov angles.

135

SK1 Data

10

SK1 Data

10

Atm. ν MCµ + candidate rings

pp → K+K + MC

8

6

4

+

K candidate rings

Atm. ν MC

20-6

pp → K+K + MC

8

6

4

2

-4

-2

0

2

4

0
-6

6

-4

-2

showering likelihood

0

2

4

6

showering likelihood

SK1 Data
12

γ candidate rings

Atm. ν MC
pp → K+K + MC

10
8
6
4
2
0
-6

-4

-2

0

2

4

6

showering likelihood

Figure 9.10: (Continued...) This set of plots shows the showering likelihoods.

SK1 Data

20

SK1 Data

14

Atm. ν MC

18

Atm. ν MC

12

pp → K+K + MC

16

pp → K+K + MCγ candidate ringsµ + candidate rings

136

10
8
6
4

14
12
10
8
6
4

2
0
0

2
100

200

300

400

500

600

0
0

700

100

200

impact parameter (cm)

300

400

500

600

700

impact parameter (cm)

SK1 Data
7

Atm. ν MC

6

pp → K+K + MCevents

5
4
3
2
1
0
0

100

200

300

400

500

600

700

800

900

1000

vertex separation of kaon decay product candidates (cm)Figure 9.10: (Continued...) This set of plots shows variables related to the geometry of
the event.

137

16

SK1 Data

14

Atm. ν MC

12

pp → K+K + MC

pp → K+K + MC

1210

eventseventsAtm. ν MC

14

10
8

8

6

6

4

4

2

2

0
0

SK1 Data

16

500

1000

1500 2000 2500

3000 3500 4000

0
-2000

4500 5000

total rtot vector magnitude (p.e.)

-1000

0

1000

2000

3000

SK1 Data
12

pp → K+K + MC

Atm. ν MC
pp → K+K + MC88eventsevents5000

SK1 Data
10

Atm. ν MC

10

4000

total rtot vector along brightest kaon direction (p.e.)

6

6

4
4
2

2
0
0.4

0.5

0.6

0.7

AFIT goodness

0.8

0.9

1

0
0

200

400

600

800

1000

1200

separation between standard and dinucleon decay vertex (cm)Figure 9.10: (Continued...) This set of plots shows variables related to the symmetry of
the light pattern of the event, and the fitting goodness using a single vertex
hypothesis.

138

2.2
SK1 Data

30

SK1 Data
2

Atm. ν MC
pp → K+K + MC

pp → K+K + MC

1.6
1.4

20

eventsevents25

Atm. ν MC

1.8

15

1.2
1
0.8

10

0.6
0.4

5

0.2
0

0

1

number of found decay electrons

2

0
0

100

200

300

400

500

600

vertex separation of decay electrons (cm)Figure 9.10: (Continued...) This set of plots shows variables related to the found decay
electrons in the event.

139

9.4.2

Training and Testing of the Boosted Decision Tree

A boosted decision tree is collection of individual decision trees that are createdin an iterative process one after another, each using reweighted training events basedon the performance of the previous tree in the chain. The complete collection ofindividual decision trees is referred to as the forest. The details of the boosted decision
tree and how it is trained are described in further detail in Appendix B.
Figure 9.11 shows one of the decision trees from the forest of the boosted decision
tree used in this analysis. In this diagram, events start at the top node and travel
down the diagram until reaching one of the “signal-like” or “background-like” terminal“leaf” nodes. An event which passes the cut at a given node is considered signal-like,and travels to the right to the next node. Likewise, an event which fails the cut at
a given node is considered background-like, and travels to the left to the next node.
The output of each individual decision tree is +1 or −1, depending on whether thepath traversed by an event ends on a signal-like or background-like terminal leaf node,
respectively.
When using the boosted decision tree to perform an analysis, each event fed into the boosted decision tree is evaluated by every tree in the forest, and the final
output of the boosted decision tree is given by the sum of each of the outputs of theindividual decision trees, weighted by a “boost-weight” that is associated with eachtree. The expression for the output, which is also described in Section B.1, is:yBDT (x) =Xi∈f orestln(αi ) · hi (x),(9.1)

where x is the tuple of input variables, and hi (x) and αi are the output and weight of
the individual trees, respectively. The boost-weight of an individual tree is determined
during the training phase, and essentially measures the tree’s overall performance in140

separating signal from background. It is described in detail in Section B.3.

!"#"$%&'($)"*"#"%(&)$"+,-./,0"1"%(%)$2"
34"-0.5/"0"1"%6"!"#"2%7A('%"*"#"%(2)B"C?C5D,9C?C,-?C,E">"&)B"
3F;GF"=/<<0C9/6"!"#"A%''($'"*"#"%(A2B"!"#$%&'()*+!"#"7$$8(%8"*"#"%(72$"
90.,:;%,<5==">"7(8'"3%"90.?+@0-":;?+=6"!"#"$82('A"*"#"%()&8"G,.5+-,=F?HI%J">"KB(&B"
3G5<<5"L09/"=F?H09;+G6"!"#"888(&8"*"#"%(2$7"!"#$%&'()*+!"#"A$&(2B"*"#"%(&$$",-%)".+!"#")B28(7B"*"#"%(B8$"E,.5+-,<?<I%J">"$B)"
3E5?+"<?<"!'%%"M0NO.6"!"#"AA7$(A%"*"#"%(2'7"!"#$%&'()*+!"#"22&'(&&"*"#"%(&'%",-%)".+!"#")&7$(&2"*"#"%(882"
E,.5+-,LCP=0:">"B88"
3=<5DD"E5?+"LCP"=0:(6"!"#")%8('&"*"#"%()B&"!"#$%&'()*+!"#")27A(7'"*"#"%('A)",-%)".+

Figure 9.11: One of the decision trees from the boosted decision tree forest that was usedin this analysis. The variable and cut value used at each branching node is
depicted, along with an interpretation of the cut. N is the total weighted
number of events at a given node. P is the signal purity of the node, definedas P = S/(S+B), where S is the weighted number of true signal events at thenode and B is the weighted number of true background events at the node.
Signal and Background indicate a terminal signal-like and background-likeleaf node, respectively.

The goal of the testing stage was to adjust the tuning parameters of the boosted
decision tree to produce output for the signal and background Monte Carlo thatsatisfied certain criteria that were decided upon beforehand. The criteria chosen forthis analysis required that the output distributions have the following qualities: (1)good overall separation between signal and background; (2) overall smoothness withno large and/or irregularly exaggerated jaggedness in either the signal or background141

distribution; (3) relative smoothness in the background and signal distributions near
the signal-like tail region of the background distribution. The emphasis on the signallike tail region of the background is due to the fact that this is approximately where
the final cut will be placed, as this search was intended to be a zero backgroundsearch. Thus it was considered desirable to have well-behaved output distributionsin that region in order to reduce systematic errors that could arise from fluctuationsin the data.
The tuning parameters were adjusted through an iterative trial process until thecriteria listed above were met. The three parameters which were tuned in the process
were the total number of trees in the forest, the minimum number of events in aterminal node, and the pruning strength. The final configuration of the boosted
decision tree is presented in Table 9.8.
Number of trees:
Min. events per leaf:
Boost type:
Separation type:
Number of bins:
Pruning method:
Pruning strength:500
1Ada-BoostGini Index
50
Cost Complexity
25.0

Table 9.8: Configuration of the boosted decision tree.

The outputs of the training, testing, and analysis stages are compared in Fig. 9.12.
The shape difference between the training and testing samples indicated that the
boosted decision tree was slightly overtrained on the training sample, meaning thatthe separation between the training signal and background distributions was somewhat better than that of the testing distributions. This does not indicate that the
boosted decision tree would have inconsistent results from sample to sample, however, which is the primary concern regarding potential bias in the performance. A142very large degree of overtraining would indicate that the performance of the boosted
decision tree will not be optimal when applied to samples that were not used in the
training itself [74].
To check for consistent performance of the boosted decision tree on samples thatwere not used in the training process, the output of the testing sample and the outputof the analysis sample were compared. Figure 9.12 demonstrates that the boosted
decision tree performance was consistent across statistically similar samples of events.
Test Bkg
2.5

Test Bkg
2.5

Test Sig
Train Bkg

2

Anal. Bkg

Train Sig

2

1.5

1.5

1

1

0.5

0.5

0

-0.3

-0.2

Test Sig

-0.1

0

0.1

0.2

Boosted Decision Tree Output

0.3

0Anal. Sig

-0.3

-0.2

-0.1

0

0.1

0.2

0.3

Boosted Decision Tree Output

Figure 9.12: Comparisons of training, testing, and analysis outputs from the boosted decision tree. Left: boosted decision tree output for signal and backgroundtraining samples (markers) and testing samples (hatched). Disagreement indicates slight (though non-problematic) overtraining of the boosted decision
tree. Right: Boosted decision tree output for analysis samples (markers) andtesting samples (hatched). Agreement indicates consistent performance ofthe boosted decision tree on statistically similar datasets.

9.4.3

Final Cut Placement for Boosted Decision Tree Output

The final stage of the multi-variate analysis was to assess the boosted decision
tree output for the analysis sample of the signal and background Monte Carlo, andto choose the placement of the final cut that would be used to determine whether or
not there was evidence for any signal events in the actual SK-I dataset.
The starting point for determining the final cut placement was to determine the√point which maximized the significance, defined as S/ S + B, where S is the signal

143

efficiency and B is the background efficiency, for 1 signal event over 26 background
events using the testing Monte Carlo samples, which would correspond to seeing 1
signal event in the 27 SK-I data events that passed the precuts. The significancecurve for these parameters is shown in Fig. 9.13. The significance was found to be
maximal for a cut placed at about 0.12.

Signal purity
Signal efficiency*purity
S / S+B

Signal efficiency

Efficiency (Purity)

Background efficiency

0.7

1

0.6
0.8

0.5

0.6

SignificanceCut efficiencies and optimal cut value

0.4
0.3

0.4

0.2
0.2

0

For 1 signal and 26 background
events the maximum S / S+B is0.6689 when cutting at 0.1254

-0.3

-0.2

-0.1

0

0.1
0.1

0.2

0.3

0

BDT output

Figure 9.13: Significance curve for the boosted decision tree output determined by testingMonte Carlo for √1 signal event over 26 background events. The significanceis defined as S/ S + B. The maximum significance is attained for a cut
placed at about 0.12.

The signal efficiency and expected background values corresponding to a cut
placed at 0.12 were calculated using the analysis Monte Carlo sample. The cut value
was then incremented by steps of 0.01, the chosen granularity for the final cut placement. Table 9.9 shows the signal efficiency and expected background for a number ofpotential cut placements that were considered. The final cut placement was chosen
to be 0.12, yielding a signal efficiency of 12.6% and an expected background of 0.28
events. It was decided to move the cut no lower than this due to the goal of having

144

zero background in the final analysis of the data. An expected background of 0.3
events or higher was deemed to be too large, given that the systematic errors on the
background estimate were expected to be somewhat large.
Cut Placement
0.10
0.11
0.12
0.13
0.14

Signal Efficiency
14.4%
13.5%
12.6%
11.5%
10.4%

Expected Background0.51 evt./1489.2 days
0.38 evt./1489.2 days
0.28 evt./1489.2 days
0.19 evt./1489.2 days
0.16 evt./1489.2 days

Table 9.9: Signal efficiency and expected background for different final cut placements on
the boosted decision tree output. Figures were calculated using the analysis
Monte Carlo sample. 0.12 was chosen for the final cut placement to be used in
the analysis of the data.

9.4.4

Boosted Decision Tree Monte Carlo Results

Figure 9.14 shows the boosted decision tree output of the atmospheric neutrino
and p p → K + K + analysis Monte Carlo samples. The results are color-coded byevent category. Figure 9.15 shows an enlarged view of the Monte Carlo output inthe region near the final cut, placed at 0.12. Events in the data with a final output
greater than the cut value would be considered dinucleon decay signal candidates.
The color legend in both figures list the event categories by decreasing abundance ineach respective plot.
The final signal efficiency and expected background after applying the boosted
decision tree cut are shown in Table 9.10. The final total efficiency for the signal
after applying the boosted decision tree cut was 12.6%. The final total expectedbackground was 0.28 events after normalizing to the SK-I livetime.

145

1.4
K+ µ + π0

K + K+ µ +1.2

K+ µ+ µ +2.5

total efficiency (%)

events / 1489.2 days

3

K+ K+ µ+2

K+ K+ γ1.5

K+ K+ µ+ γ
K K π
+

1

+

0

K K µ π
+

0.5

+

+

0

K K µ µ
+

0
-0.4

-0.3

-0.2

-0.1

0

0.1

+

0.2

+

0.3

K + µ + π0
K + K+ µ + µ+

0.8

K + K+ µ + γ

0.6

K + K+ γ
0.4

K + K+ π00.2+

boosted decision tree output: Atm. ν MC

K + µ + µ+1

K + K+ µ + π0

0
-0.4

0.4

-0.3

-0.2

-0.1

0

0.1

0.2

+

0.3

+

boosted decision tree output: pp → K K MC

0.4

Figure 9.14: Boosted decision tree output for the Monte Carlo. Left: Atmospheric neutrino Monte Carlo. Right: p p → K + K + Monte Carlo. The legend lists
the event categories by decreasing abundance in each respective plot. The
vertical line indicates the final cut, placed at 0.12.

1.4
K+ µ + π0

1.4

K µ+ µ +1.2

K K µ+++

1

K+ K+ γ0.8

K+ K+ µ+ γ

0.6

K+ K+ π00.4

K+ K+ µ+ π0

0.2

K+ K+ µ+ µ +0

K + K+ µ +1.2+0.02 0.04 0.06 0.08

0.1 0.12 0.14 0.16 0.18

0.2 0.22

boosted decision tree output: Atm. ν MC

total efficiency (%)

events / 1489.2 days

1.6

1
0.8
0.6

K + µ + µ+
K + µ + π0
K + K+ µ + µ+
K + K+ µ + γ
K + K+ γ

0.4
0.2
0

K + K+ π0
K + K+ µ + π00.02 0.04 0.06 0.08

0.1 0.12 0.14 0.16 0.18+

0.2 0.22+

boosted decision tree output: pp → K K MC

Figure 9.15: Boosted decision tree output for the Monte Carlo in the region near the final
cut. Left: Atmospheric neutrino Monte Carlo. Right: p p → K + K + Monte
Carlo. The vertical line indicates the final cut, placed at 0.12.

9.4.5

Final Background Characteristics

Table 9.11 shows the surviving background Monte Carlo events after applying
the boosted decision tree cut, broken down by neutrino flavor and weak current type.
Charged current (CC) and neutral current (NC) muon neutrino interactions dominate
the background after the precuts.
Table 9.12 shows the specific neutrino interaction types of the final surviving
background Monte Carlo. Single-pion and multi-pion production from delta resonance

146

Event Category
K + K + µ+
K + µ+ µ+
K+ K+ γK + K + µ+ µ+
K+ K+ γ γ
K + µ+ γ γ
K + K + µ+ γ
K + K + µ+ γ γ
total

p p → K + K + M.C.
Atmospheric ν M.C.
efficiency (%) events / SK-I livetime
5.0
0.07
1.8
0.08
0.8
0.02
2.5
0.00
0.3
0.00
1.0
0.08
1.0
0.03
0.2
0.00
12.6
0.28

Table 9.10: Signal and background after applying the boosted decision tree cut.

ν Mode
CC νµ
NC νµ
CC νe
NC νe
total

Events/200 Years Events/1489.2 Days%
8.73
0.18 643.00
0.06 222.00
0.04 15
0.00
0.00
0
13.73
0.28 100

Table 9.11: Final remaining background Monte Carlo events sorted by neutrino flavor
and weak current type. CC: Charged Current; NC: Neutral Current. The
first numerical column shows the unscaled number of events from the 200
year analysis sample. The second numerical column shows the number of
events scaled to the SK-I livetime.

are the dominant sources of background after the final cut on the boosted decision
tree output.
Table 9.13 shows the breakdown of the true particle types (PID) for each of thedifferent ring classifications from the final background sample. A plurality of the kaoncandidate rings were created by muons. The rest were mostly created by charged
pions and protons. Half of the muon candidate rings were created by charged pions.
A majority of the gamma candidate rings were created by actual gammas.

147Interaction
CC multi-pion production
CC single pion delta resonanceNC single pion delta resonanceNC diffractive pion production
CC quasi-elastic
NC multi-pion production
CC diffractive pion production
NC elastic
CC coherent pion production
NC coherent pion production
total

Events/200 Years Events/1489.2 Days%
5.34
0.11 38.9
4.89
0.10 35.6
2.00
0.04 14.6
1.00
0.02 7.3
0.50
0.01 3.6
0.00
0.00 0.0
0.00
0.00 0.0
0.00
0.00 0.0
0.00
0.00 0.0
0.00
0.00 0.0
13.73
0.28 100

Table 9.12: Final remaining background Monte Carlo events sorted by neutrino interaction type. CC: Charged Current; NC: Neutral Current. The first numericalcolumn shows the unscaled number of events from the 200 year analysis sample. The second numerical column shows the number of events scaled to the
SK-I livetime.

True PID
p
π±
µ±
e±
γ
total

K + Candidates µ+ Candidates γ Candidates
5
5
0
8
11
1
12
2
0
1
1
3
1
3
6
27
22
10

Table 9.13: True PID of all rings in final background sample.

9.4.6

Final Signal Characteristics

Table 9.14 shows the breakdown of the true particle types (PID) for each of thedifferent ring classifications from the final signal sample. The kaon candidates had
a purity of 92.6%, the muon candidates had a purity of 81.8%, and the gamma
candidates had a purity of 67.4%. These high purities show that the classification
algorithm performed quite well on the dinucleon decay signal Monte Carlo events.
Table 9.15 shows the purity of each of the event categories for the final signal

148

True PID
K+
π+
µ+
e+
γ
total
purityK + Candidates µ+ Candidates γ Candidates
4273
276
27
57
128
32
224
2668
200
19
51
62
42
138
663
4615
3261
984
92.6%
81.8%
67.4%

Table 9.14: True PID of all rings in final signal sample.

sample. The purity is defined here as the fraction of events of a given category inwhich every ring in the event was correctly classified. The final entry represents thetotal purity of the final signal sample across all categories. The K + K + µ+ and K +K + µ+ µ+ categories both had very high purity; these represent the majority of theevents in the final signal sample. The good performance in these categories is dueto the two found kaons. The rest of the categories had purities ranging from 3050%. The reason that these categories had somewhat worse purities than K + K + µ+
and K + K + µ+ µ+ is that they either were missing the second kaon ring, or because
they contained gamma candidate rings, which have a worse purity than either kaonor muon candidate rings. Having correctly classified every ring in about two-thirdsof the final signal sample is here considered a very good performance benchmark,given the highly complicated nature of the signal events and the ring classification
algorithm that it entailed.

9.5

Systematic Uncertainties

The systematic uncertainties were broken down into four main categories: Monte
Carlo simulation, detector knowledge, event reconstruction, and boosted decision tree

149

Event Category Purity
K + K + µ+
84.7%
+ + +
K µ µ
43.5%
+
+
K K γ
39.6%
K + K + µ+ µ+
81.7%
K+ K+ γ γ
41.8%
K + µ+ γ γ
31.1%
+
+ +
K K µ γ32.5%+
+ +K K µ γγ
48.7%
overall
65.5%
Table 9.15: Event categorization purity for the final signal sample.

related uncertainties. The first category, Monte Carlo simulation, had a unique setof uncertainties for the signal efficiency and expected background. The latter three
categories were common to both signal and background.

9.5.1

Monte Carlo Simulation Uncertainties

SignalThe four main sources of uncertainty from the Monte Carlo for the signal efficiency were from hadronic interactions of the kaons, the simulation of correlated
decay events, the Fermi momentum of the parent nucleons, and the branching ratiosof the kaon decay channels.
Based on cross-section comparisons between the dinucleon decay Monte Carlo andthose compiled by the PDG [12], the systematic uncertainty for hadronic interactionsof the kaons was estimated to be 25%. Figure 5.4 shows the cross-sections from thep p → K + K + Monte Carlo overlaid on the experimental data points taken from the
PDG.
The effect of correlated decay is described in Section 5.1. A 100% systematic

150

uncertainty was taken for this effect. The uncertainty in the Fermi momentum wastaken to be 20%, following [4]. The uncertainty on the kaon branching ratios was
assumed to be negligibly small.
Background
For the atmospheric neutrino background, there were also four main systematic
ucertainty terms taken into account. All four of these effects were taken to have a
direct effect on the final background uncertainty. The neutrino flux was taken to have
an 8% uncertainty, the neutrino interaction cross-section 15%, nuclear effects weretaken to have a 20% uncertainty, and pion propagation simulation uncertainties, which
are the dominant source overall in the background estimate uncertainty, were taken
to have a 50% uncertainty. All of these uncertainties follow the general perscriptions
laid out in previous Super-Kamiokande studies [4, 57].

9.5.2

Detector Related Uncertainties

There were two main systematic uncertainties related to detector knowledge: theuncertainty in the fiducial volume and the energy scale.
A 4.4% uncertainty in the fidcucial volume was taken for signal and background,based on comparisons between the data and atmospheric Monte Carlo, which corresponded with a shift in the fiducial volume boundary cut of 200 cm (+22.2 cm or-21.5 cm) from the inner wall of the detector.
The energy scale uncertainty was also estimated based on comparisons between
data and atmospheric neutrino Monte Carlo, and was estimated to be at the 1.1%
level. A corresponding 1.1% shift in the P OT OT precut (Section 9.3) was used toestimate the overall uncertainty in the signal and background based on the energy

151

scale uncertainty.

9.5.3

Reconstruction Uncertainties

Six major sources of uncertainty were identified from the event reconstruction
process: showering likelihood, ring counting, Cherenkov angle, vertex position, decay
electron counting, and ring momentum reconstruction. The estimated size of theuncertainties for each of these sources is summarized in Table 9.16. The same values
were used for both signal and background uncertainty estimates.
Sytematic Source
showering likelihoodring counting
Cherenkov angle
vertex position
decay-e counting
ring momentum

Uncertainty
5.8% (0.2)
< 1%
0.7◦
30 cm
∼ 1%
1.2%

Table 9.16: Systematic uncertainties in event reconstruction variables.

Each of the respective variables was shifted by its corresponding uncertainty, andthe entire search process was re-run on the shifted Monte Carlo, from the ring classification to the boosted decision tree. The boosted decision tree was not re-trainedfor each shifted sample; the original boosted decision tree trained on the unshifted
Monte Carlo was used to estimate all uncertainties.

9.5.4

Search Method Uncertainties

Finally, to attempt to estimate the systematic uncertainties introduced by the
boosted decision tree, several techniques were used, and the results of each wereaveraged independently for the signal and background. The first technique was to152

simply shift the final cut placement up and down by an increment of 0.01. This
was the granularity used when choosing the original cut placement, as described in
Section 9.4.3.
The second technique was to assess the final signal efficiency and expected background after applying the trained boosted decision tree to the testing Monte Carlosample (see Section 9.4.2), and then the combined testing and analysis samples, simultaneously.
The results of these tests are summarized in Table 9.17. The average uncertaintyof the two techniques was taken as the final uncertainty of the boosted decision tree
performance.

+

pp → K K
Atm. ν+Cut Variation Sample Variation Average
8.1%
1.4%
4.7%
33.9%
26.8%
30.4%

Table 9.17: Systematic uncertainty of boosted decision tree performance.

9.5.5

Total Systematic Uncertainty

The final systematic uncertainties are listed in Table 9.18. The total uncertainty
on the signal detection efficiency was calculated to be 25.2%. The total uncertainty
on the expected background was calculated to be 68.1%. The leading sources of
uncertainty for both the signal efficiency and the expected background came from therespective Monte Carlo simulations. For the signal, the uncertainty in the simulationof the Fermi momentum of the original nucleons dominated. For the background, theuncertainty in the simulation of the propagation of pions was the dominant source.

153

Systematic Source
σ[ǫsig ] σ[Nbkg ]++p p → K K Monte Carlohadronic interactions
1.3%
correlated decay
2.5%
Fermi momentum
24.1%
K + decay B.R.
< 1%
Atm. ν Monte Carloν flux8.0%
ν cross-section
15.0%
π nuclear effect
20.0%
pion propagation
50.0%
Detector Knowledge
fiducial volume
2.2%
5.4%
energy scale
< 1%
< 1%
Event Reconstruction
showering likelihood
2.4%
14.3%
ring counting
< 1%
< 1%
Cherenkov angle
1.0%
3.6%
vertex position
3.9%
16.1%
decay-e counting
< 1%
< 1%
ring momentum
0.7%
5.4%
Boosted Decision tree
sample bias
4.7%
30.4%
total
25.2%
68.1%
Table 9.18: Systematic uncertainties in event reconstruction variables.

Chapter 10

Results and Conclusions
The 1489.2 days of fully contained SK-I data, corresponding to 91.7 kiloton · yearsof exposure, were searched for evidence of 16 O(pp) → 14 C K + K + . Twenty-seven
events in the data passed the precuts described in Section 9.3. Figure 10.1 shows the
final results of the search. No signal candidate events were found in the data, definedby the final cut on the boosted decision tree output at 0.12. The distribution of theboosted decision tree output for the SK-I data was found to be consistent with theatmospheric neutrino Monte Carlo, representing the background of the search.
The amount of background expected to pass the final cut was 0.28 ± 0.19 events
after normalizing to the 1489.2 days of SK-I livetime. The p p → K + K + signal
detection efficiency was estimated to be 12.6% ± 3.2%.
Following the conventions described in Appendix C.1, the results of the searchwere used to calculate a lower limit on the partial lifetime of the dinucleon decay
channel p p → K + K + :τ≥ 1.7 × 1032 years.
B.R. p p→K + K +(10.1)

This number represents the lifetime limit per 16 O nucleus at a 90% confidence level.
154events / 1489.2 days of livetime

155

SK1 Data
7

Atm. ν MC

6

pp → K+K + MC

5
4
3
2
1
0
-0.4

-0.3

-0.2

-0.1

0

0.1

0.2

0.3

0.4

boosted decision tree output

Figure 10.1: Final output of boosted decision tree. SK-I data is overlaid on analysis Monte
Carlo samples of the dinucleon decay signal and the atmospheric neutrino
background. The final cut placement is at 0.12; events with an output above
this threshold are considered signal candidates. No signal candidate eventswere found in the data.

The limit presented in this dissertation is the result of the first direct experimental
search for dinucleon decay into kaons. This limit is two orders of magnitude better
than the dinucleon decay limits for pion and lepton modes set by the Frejus experiment [12, 15], and seven orders of magnitude better than the limits for neutrino
modes and invisible modes set by Frejus, DAMA, and the BOREXINO CTF [12–15].
Using the lower lifetime limit from the search, one can then calculate the corresponding upper limit on the magnitude of the ∆B = 1 R−Parity violating SUSY
parameter, λ′′uds . The details of the calculation are described in Appendix C.2. The
upper limit on the λ′′uds parameter based on the first experimental lower lifetime limit
for the dinucleon decay channel p p → K + K + was calculated to be:|λ′′uds | <0.64
τ [years]

1/4= 7.8 × 10−9 .

(10.2)

156

This new result is more restrictive than the limit of ∼ 10−7 estimated by Goity andSher [22] based on a typical nuclear matter lifetime limit of ∼ 1030 years taken fromthe limits observed in searches for single nucleon decay and the searches for dinucleon
decay into pions or leptons performed by the Frejus collaboration.
According to the literature [20–22], a search for dinucleon decay into kaons should
provide the best currently accessible experimental limit on λ′′uds . Because λ′′uds scaleswith τ −1/4 , and τ ∝ ǫ · A, where ǫ is the signal detection efficiency and A is the
exposure in kiloton · years, one would need to increase ǫ · A by a factor of 10,000 to
improve the limit calculated in this analysis by a single order of magnitude. Thus,
adding the SK-II and SK-III datasets to the search would do little to improve thelimit obtained here from the SK-I dataset alone.
In regard to the ability of future experiments to improve the λ′′uds limit presented
here by an order of magnitude, it seems that it would be very difficult to accomplishusing conventional nucleon decay detection methods, though perhaps not strictly
impossible. Certainly, the target mass of the detector (or set of detectors) would
have to be in the range of megatons. Until such a detector exists, however, the newlimit on λ′′uds presented in this dissertation is expected to remain the best limit onrecord.

Appendix A

MVFIT Resolutions
The resolution for each reconstructed variable was obtained separately for thep p → K + K + Monte Carlo and the atmospheric ν Monte Carlo. Three of the precuts
listed in Section 9.3 were applied to both Monte Carlo samples before calculating theresolutions:• 1000 < potot < 11000, where potot is the total amount of light seen in the ID
in units of p.e.
• 3 ≤ nring ≤ 5, where nring is the number of found rings in the event• wall > 200 cm, where wall is the distance from the event vertex to the nearest
wall; this is the fiducial volume (FV) cut
A correction of +4.5%, +0.5%, and +5.0% was applied to the kaon-like, muon-like,and gamma-like reconstructed momenta, respectively. These corrections compensate
for the bias arising from the RTOT-to-momentum look-up tables used in the momentum reconstruction process (see Section 7.1.4), which were generated for RTOT
values corresponding to a 70◦ half-opening angle. The RTOT value used in MVFIT

157

158

is instead calculated using the unmasked window, defined as the ring’s Cherenkov
angle +10◦ (see Section 7.2).
A correction of -5 cm and -40 cm was also applied to the reconstructed vertex along
the ring direction for muon-like and gamma-like vertices, respectively. This was again
to compensate for bias introduced by using reconstruction algorithms that were notoriginally developed to perform with the light masking techniques used by MVFIT.
The momentum and vertex correction were both applied before generating theresolution tables and plots in this appendix (but were themselves based on the sametables calculated before the corrections were applied). The directional, vertex, andperpendicular (⊥) vertex resolutions all represent the minimal point that includes
68% of the entries. The mean and variance for the momentum and parallel (||) vertex
resolutions were calculated using Gaussian fits.
In all of the plots in this chapter, the solid lines represent the p p → K + K +Monte Carlo and the dashed lines represent the atmospheric neutrino Monte Carlo.
The vertical lines represent the 68% resolution point. In all plots, “∆ variable name”
indicates the reconstructed value of the variable minus the true value of the variable.
∆ vertex || indicates the difference between the reconstructed and true vertex along
the direction parallel to the true direction of the particle. Likewise, ∆ vertex ⊥indicates the difference between the reconstructed and true vertex in the plane perpendicular to the true direction of the particle.

159

A.1

Charged Kaon Ring Resolutions

The MVFIT reconstructed variable resolutions for rings created by charged kaons
are listed in Table A.1, and shown in Fig. A.1. 409 kaon rings from the atmospheric
neutrino Monte Carlo and 46,097 kaon rings from the p p → K + K + Monte Carlo
were analyzed.
K ± rings
Atm. ν p p → K + K +direction
resolution
4◦
4◦
momentum
mean
-12%
0%
variance
11%
6%
vertexresolution
160 cm130 cm⊥ resolution
20 cm
40 cm
|| mean-68 cm-27 cm|| variance
114 cm88 cm
Table A.1: MVFIT reconstructed resolutions for true charged kaon rings. ⊥ indicatesthe direction perpendicular to the ring direction, and || indicates the direction
parallel to the ring direction.

160

20

Atm. ν MC × 100
pp → K+K + MC

16

pp → K+K + MCtrue K ± rings / 1000

true K ± rings / 1000

Atm. ν MC × 100

14

18

14
12
10
8
6

12
10
8
6
4

4
2

2
0
0

10

20

30

4.5

0
-1

50

-0.8

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

∆ momentum / true momentum

Atm. ν MC × 100

0.8

1Atm. ν MC × 100
124

pp → K+K + MC

pp → K+K + MC

3.5

true K ± rings / 1000

true K ± rings / 1000

40

∆ direction (deg)

3
2.5
2
1.5

10
8
6
4

1
2

0.5
50

100

150

200

250

300

∆ vertex (cm)true K ± rings / 1000

0
0

350

400

450

0
0

500

50

100

150

∆ vertex

4.5

Atm. ν MC × 100

4pp → K+K + MC

200

250

300

350

400

450

500

to ring direction (cm)

3.5
3
2.5
2
1.5
1
0.5
0
-500

-400

-300

-200

-100

0

100

∆ vertex || (cm)

200

300

400

500

Figure A.1: MVFIT reconstructed variable resolutions for charged kaons. In order fromleft to right, top to bottom: direction, momentum, vertex, vertex perpendicular to ring direction, vertex parallel to ring direction.

161

A.2

Muon Ring Resolutions

The MVFIT reconstructed variable resolutions for rings created by muons arelisted in Table A.2, and shown in Fig. A.2. 25,357 muon rings from the atmospheric
neutrino Monte Carlo and 33,775 muon rings from the p p → K + K + Monte Carlo
were analyzed.
µ± rings
Atm. ν p p → K + K +direction
resolution
3◦
6◦
momentum
mean
+2%
-1%
variance2%
9%vertexresolution
35 cm
140 cm⊥ resolution
15 cm
65 cm
|| mean+4 cm
0 cm|| variance
23 cm
45 cm
Table A.2: MVFIT reconstructed resolutions for true muon rings. ⊥ indicates the direction perpendicular to the ring direction, and || indicates the direction parallelto the ring direction.

162

60

Atm. ν MC

true µ ± rings / 1000

true µ ± rings / 1000

pp → K+K + MC × 10

40
30
20

10

50
40
30
20
10

0
0

10

20

30

40

∆ direction (deg)

0
-1

50

20

Atm. ν MC

18

pp → K+K + MC × 10

-0.8

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

∆ momentum / true momentum

0.8

1Atm. ν MC

35

pp → K+K + MC × 10

16true µ ± rings / 1000

true µ ± rings / 1000

Atm. ν MC

60

pp → K+K + MC × 10

50

14
12
10
8
6

30
25
20
15
10

4
5

2
0
0

50

100

150

200

250

300

∆ vertex (cm)

350

400

450

0
0

500

50

100

150

∆ vertex

200

250

300

350

400

450

500

to ring direction (cm)Atm. ν MC

true µ ± rings / 100030

pp → K+K + MC × 10

25
20
15
10
5
0
-500

-400

-300

-200

-100

0

100

∆ vertex || (cm)

200

300

400

500

Figure A.2: MVFIT reconstructed variable resolutions for muons. In order from left toright, top to bottom: direction, momentum, vertex, vertex perpendicular toring direction, vertex parallel to ring direction.

163

A.3

Gamma Ring Resolutions

The MVFIT reconstructed variable resolutions for rings created by gammas orelectrons are listed in Table A.3, and shown in Fig. A.3. 168,935 gamma and electronrings from the atmospheric neutrino Monte Carlo and 56,972 gamma and electronrings from the p p → K + K + Monte Carlo were analyzed.
γ, e± rings
Atm. ν p p → K + K +direction
resolution
6◦
8◦
momentum
mean
-5%
+5%
variance
25%
29%
vertexresolution
110 cm140 cm⊥ resolution
20 cm
40 cm
|| mean+67 cm+63 cm|| variance
68 cm
79 cm
Table A.3: MVFIT reconstructed resolutions for true gamma and electron rings. ⊥ indicates the direction perpendicular to the ring direction, and || indicates the
direction parallel to the ring direction.

164

Atm. ν MC

35

Atm. ν MC

pp → K K MC × 1030

pp → K+K + MC × 10

+

+

100

true γ rings / 1000

true γ rings / 1000

120

80
60
40
20

20
15
10
5

0
0

10

20

30

40

∆ dir (deg)

0
-1

50

-0.8

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

∆ momentum / true momentum

0.8

Atm. ν MC

160

Atm. ν MC

pp → K+K + MC × 10

140

pp → K+K + MC × 10

true γ rings / 1000

50

true γ rings / 1000

25

40

30

20

1

120
100
80
60
40

10
20
0
0

50

100

150

200

250

300

∆ vertex (cm)

350

400

450

0
0

500

60

100

150

∆ vertex

200

250

300

350

400

450

500

to ring direction (cm)Atm. ν MC
pp → K+K + MC × 10

50

true γ rings / 1000

50

40
30

20
10

0
-500

-400

-300

-200

-100

0

100

∆ vertex || (cm)

200

300

400

500

Figure A.3: MVFIT reconstructed variable resolutions for gammas and electrons. In order from left to right, top to bottom: direction, momentum, vertex, vertex
perpendicular to ring direction, vertex parallel to ring direction.

Appendix B

Boosted Decision Tree
The ROOT-based TMVA [74] multivariate analysis software package was used to
implement the boosted decision tree in this analysis.
A decision tree is a classifier with a top-down binary tree structure. Starting withthe first node in the tree, all variables are scanned and a one-dimensional cut is placedon the variable which allows for the maximal separation of signal and background.
The events are then split into two groups, one passing the cut and traveling down the“signal-like” branch to one node, and the other failing the cut and traveling down the“background-like” branch to another node. The same process is then repeated foreach of the two new daughter nodes. This pattern is repeated until a certain stopping
criterion is reached, at which point a node is no longer split into two branches, but isinstead declared either a signal-like or background-like terminal leaf node, dependingon the makeup of the events which it contains. The process of creating of a tree issynonymous with training a tree.
Boosting refers to the process of training a classifier using a reweighted, or “boosted,”
training event sample that was the output of a previously trained classifier. In thecase of a boosted decision tree, the output of tree n is reweighted based on the per165

166

formance of that tree. Then, those boosted events are used to create tree n + 1. Notethat the output of tree n includes all of the original training events; no events arethrown out or go unused in any tree created during the training process.
The forest is the name given to the full collection of decision trees that are createdduring the training process. The forest itself is simply referred to as the boosted
decision tree, and its output is a weighted majority of the decision trees from which
it is comprised.

B.1

Final Output

The calculation for the final boosted decision tree output, yBDT (x), is shown below:

yBDT (x) =Xi∈f orestln(αi ) · hi (x)(B.1)

where x is the tuple of input variables, and hi (x) and αi are the output and weightof the individual trees, respectively. The calculation of the weight for a given tree is
described in Section B.3.

B.2

Separation Index

The branch node splitting is a one-dimensional cut on a single variable from thetuple of input variables. The binning for the variables is specified in the boosted
decision tree tuning parameters, and the trial locations for the branching cut lie
between the bins. All options are considered, and the variable and cut position which
maximize the separation index between the current node and those of the hypotheticaldaughter nodes are chosen. The hypothetical daughter nodes are weighted by theirrelative fraction of events in this calculation.

167

The separation index, or Gini index, is defined below, where P is the purity of
the node, N is the number of training events in the node, and Wi is the weight of theindividual events:

P =P
Gini =P

SW
PS

S WS +

N
X
i=1

B WB

!

Wi P · (1 − P )

(B.2)

(B.3)

The quality of separation, Qsep , is maximized at each branching node, and is given
by:

Qsep = Gparent − (Gbkg + Gsig )(B.4)

where Gparent is the Gini index of the parent node, and Gbkg and Gsig are the Gini
indices of the background-like and signal-like daughter nodes, respectively.

B.3

Adaptive Boost Method

The boosting method that was used in this analysis is the adaptive boost method(“Ada-Boost”). In this method, misclassified events (i.e., signal events which end up
in a background leaf, or vice versa) are given a larger event weight in the trainingof the next tree. All of the misclassified events in a single tree are multiplied by acommon boost weight, α, where err is the mis-classification rate of the tree.

α=1 − err
err

(B.5)

168

B.4

PruningPruning is a process by which statistically insignificant nodes are removed froma tree after it has been grown to maximum size (i.e., all branches have ended interminal leaf nodes). In this analysis, the expected error pruning method was used,which removes all daughter nodes for which the statistical error estimate of the parent
node is smaller than the combined statistical error estimates of the daughter nodes.
Thus a branching node may become a leaf node.
The statistical error estimate of each node, ǫnode , is calculated using the binomial
error, as shown below:ǫnode = F ·r

P · (1 − P )
N

(B.6)

where P is the purity of the node, N is the number of training events in the node,
and F is the pruning strength factor. The pruning criterion is then given by:ǫparent < ǫbkg + ǫsig(B.7)

where ǫparent is the statistical error of the parent node, and ǫbkg and ǫrmsig are thestatistical errors of the background-like and signal-like daughter nodes, respectively.

B.5

Training Distributions

The distributions of each of the input variables used in the training of the boosted
decision tree can be seen in Fig. B.1. See Section 9.4.1 for a complete description ofeach of the input variables.

169

0.001

1600 1650 1700 1750 1800 1850 1900

-50

0

0.001
300

400

500

0.002
0.00150.001
0.0005
0

200

400

600

Input variable: g_cand_vtxsep

Input variable: k_cand_mom[0](1/N) dN / 15.1

U/O-flow (S,B): (0.0, 0.0)% / (0.0, 0.1)%

0.003
0.00250.002
0.00150.001
0.0005
0

200

400

600

800

400

0.01
0.008
0.007
0.006
0.005
0.004
0.003
0.002
0.001
0

100

200

300 400

0.006
0.005

0.003
0.002
0.001

0.16
0.14
0.12
0.1

0.06
0.04
0.02
0

10

20

30

40

Input variable: k_cand_mom[1]Input variable: k_cand_cang[1](1/N) dN / 0.243(1/N) dN / 5.77

(1/N) dN / 1.09

0.1
0.05
0

-4

-2

0

2

4

6

0.004
0.003
0.002
0.001

8

500

550

600

650

700

0.04
0.03
0.02
0.01

800

10

20

30

40

50

Input variable: k_cand_show[1]Input variable: mu_cand_mom[0]Input variable: mu_cand_cang[0](1/N) dN / 1.09

0.15
0.1
0.05
-10

-8

-6

-4

-2

0

2

4

0.016
0.014
0.012
0.01
0.008
0.006
0.004
0.002

6

0.08
0.07
0.06
0.05
0.04
0.03
0.02
0.01
0
10

120 140 160 180 200 220 240 260 280 300

20

30

40

50

Input variable: mu_cand_impact[0]Input variable: mu_cand_mom[1](1/N) dN / 3.56

0.25
0.2
0.15
0.1
0.05
0
-10

-8

-6

-4

-2

0

2

4

6

mu_cand_show[0]

0.007
0.006
0.005
0.004
0.003
0.002
0.001
0

0

200

400

600

800

70

0.009
U/O-flow (S,B): (0.0, 0.0)% / (0.0, 0.0)%

0.3

U/O-flow (S,B): (0.0, 0.0)% / (0.0, 0.0)%

U/O-flow (S,B): (0.0, 0.0)% / (0.0, 0.0)%

0.35

60

mu_cand_cang[0]Input variable: mu_cand_show[0](1/N) dN / 22

mu_cand_mom[0](1/N) dN / 0.297

k_cand_show[1]70U/O-flow (S,B): (0.0, 0.0)% / (0.0, 0.0)%

0.2

0.018
U/O-flow (S,B): (0.0, 0.0)% / (0.0, 0.0)%

U/O-flow (S,B): (0.0, 0.0)% / (0.0, 0.0)%

0.25

60

k_cand_cang[1](1/N) dN / 3.06

k_cand_mom[1](1/N) dN / 0.295

k_cand_show[0]

750

0.05
U/O-flow (S,B): (0.0, 0.0)% / (0.0, 0.0)%

0.2
0.15

U/O-flow (S,B): (0.0, 0.0)% / (0.0, 0.0)%

U/O-flow (S,B): (0.0, 0.0)% / (0.0, 0.0)%

0.3
0.25

0.006
0.005

50

k_cand_cang[0]Input variable: k_cand_show[0]

0.4

700

0.18

k_cand_mom[0]

0.45

600

0.08

700 800 900 1000 1100 1200 13001400 1500

g_cand_vtxsep

0.35

500

Input variable: k_cand_cang[0]

0.007

0

600

rec_pi0_impact

0.008

0.004

1000

500

0.009

1000

U/O-flow (S,B): (0.0, 0.0)% / (0.0, 0.0)%

0.005
0.00450.0035300

rec_pi0_mom(1/N) dN / 20.2

rec_pi0_mass

800(1/N) dN / 13.1

0.0025

600

0.004

200

Input variable: rec_pi0_impact

0.003
U/O-flow (S,B): (0.0, 0.0)% / (0.0, 0.0)%

0.002

100

k_cand_vtxsep(1/N) dN / 0.889

0.003(1/N) dN / 17.8

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

0.004

200

0

Input variable: rec_pi0_mom

0.005

100

0.001
0.0005rec_dinuc_mom

Input variable: rec_pi0_mass

0

0.002
0.001550 100 150 200 250 300 350 400

rec_dinuc_mass

U/O-flow (S,B): (0.0, 0.0)% / (0.0, 0.0)%

0.00150.003
0.0025

U/O-flow (S,B): (0.0, 0.0)% / (0.1, 0.1)%

0.003
0.00250.002

0.0035

U/O-flow (S,B): (0.0, 0.0)% / (0.0, 0.0)%

0.00350.0025(1/N) dN / 11.7

0.004

Input variable: k_cand_vtxsep

0.003
U/O-flow (S,B): (0.0, 0.0)% / (0.0, 0.0)%

0.0045(1/N) dN / 7.6

Input variable: rec_dinuc_mom

SignalBackground0.005

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

Input variable: rec_dinuc_mass

0.008
0.007
0.006
0.005
0.004
0.003

1000 1200

120 140 160 180 200 220 240 260 280 300 320

mu_cand_impact[0]mu_cand_mom[1]Figure B.1: Distributions of the input variables used to train the boosted decision tree.
Trivial bins have been truncated to emphasize shape difference in the nontrivial range. Filled blue histograms indicate the p p → K + K + signal Monte
Carlo. Hatched red histograms indicate the atmospheric neutrino backgroundMonte Carlo.

170

0.005
10

20

30

40

50

0.12
0.1
0.08
0.06
0.04
0.02

60

-10

-8

-6

-4

-2

0

2

4

6

0.003
0.00250.002
0.00150.001
0.0005

8

0

200

400

600

Input variable: g_cand_mom[0]Input variable: g_cand_cang[0]Input variable: g_cand_show[0](1/N) dN / 0.267

0.00150.001
0.0005
0

200

400

600

800

0.08
0.06
0.04
0.02

1000

20

25

30

35

40 45

50

Input variable: g_cand_impact[0]Input variable: g_cand_mom[1](1/N) dN / 11.5

0.008
0.006
0.004
0.002

0

200

400

600

60

0.1
0.05

65

-10

-8

-6

-4

-2

0

0.005
0.004
0.003
0.002
0.001
0

100

200

300

400

500

0.1
0.08
0.06
0.04
0.02

600

20

30

40

50

Input variable: g_cand_show[1]Input variable: g_cand_impact[1]Input variable: g_cand_totmom(1/N) dN / 17.8

-10

-8

-6

-4

-2

0

2

4

0.002
0.001

6

0

0.00150.001
0.0005
200

400

600

800

16
14
12
10
8
6
4
2
0
0

1000

0.2 0.4 0.6 0.8

0.005(1/N) dN / 0.0144

1 1.2 1.4 1.6 1.8

0.001
0

200 400 600 800100012001400160018002000

off_vtx_sep

6
5
4
3
2
1
0.2

0.3

0.4

1000

0.001
0.00080.00060.0004
0

100

200

300

400

500

600

dcy_e_vtxsep
Input variable: total_rtot_dot_k

7

0

800

0.0012

2

U/O-flow (S,B): (0.0, 0.0)% / (0.0, 0.0)%

U/O-flow (S,B): (0.0, 0.0)% / (0.1, 0.2)%(1/N) dN / 36.6

Input variable: agood

0.002

600

n_dcy_e

Input variable: off_vtx_sep

0.003

400

Input variable: dcy_e_vtxsep

18

decay_vtxsep

0.004

200

g_cand_totmom

U/O-flow (S,B): (0.0, 0.0)% / (0.0, 0.0)%

0.002(1/N) dN / 0.0339

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

0.002500Input variable: n_dcy_e

0.004

0.003

0.001
0.0005g_cand_impact[1]Input variable: decay_vtxsep

0.00350.002
0.0015100 200 300 400 500 600 700 800

g_cand_show[1]U/O-flow (S,B): (0.0, 0.0)% / (0.0, 0.0)%

0.003

0.0025

U/O-flow (S,B): (0.0, 0.0)% / (0.0, 0.0)%

0.05

0.004

70

0.5

0.6

0.7

0.8

0.9

agood

0.0014U/O-flow (S,B): (0.0, 0.0)% / (0.0, 0.0)%

0.1

0.005

(1/N) dN / 11

0.15

0.006

(1/N) dN / 130

0.2

0.007
U/O-flow (S,B): (0.0, 0.0)% / (0.0, 0.0)%

U/O-flow (S,B): (0.0, 0.0)% / (0.0, 0.0)%

0.25

60

g_cand_cang[1](1/N) dN / 15.7

g_cand_mom[1]

0.3

6

0.12

(1/N) dN / 0.269

g_cand_impact[0]4Input variable: g_cand_cang[1]

0.006

800

2

g_cand_show[0]U/O-flow (S,B): (0.0, 0.0)% / (0.0, 0.0)%

U/O-flow (S,B): (0.0, 0.0)% / (0.0, 0.0)%

0.01

55

0.2
0.15

g_cand_cang[0](1/N) dN / 16.9

g_cand_mom[0]0.25U/O-flow (S,B): (0.0, 0.0)% / (0.0, 0.0)%

0.002(1/N) dN / 0.915

0.00250.1

U/O-flow (S,B): (0.0, 0.0)% / (0.0, 0.0)%

U/O-flow (S,B): (0.0, 0.0)% / (0.0, 0.0)%

0.003

0.12
U/O-flow (S,B): (0.0, 0.0)% / (0.0, 0.0)%

0.0035800

mu_cand_impact[1](1/N) dN / 0.806

mu_cand_show[1](1/N) dN / 17.4

mu_cand_cang[1](1/N) dN / 17.1

0.14

U/O-flow (S,B): (0.0, 0.0)% / (0.0, 0.0)%

0.01

Input variable: mu_cand_impact[1]0.16
U/O-flow (S,B): (0.0, 0.0)% / (0.0, 0.0)%

0.02
0.015

(1/N) dN / 0.305

Input variable: mu_cand_show[1]0.025

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

Input variable: mu_cand_cang[1]0.00120.001
0.00080.00060.00040.00020-2000 -1000

01000 2000 3000 4000 5000

total_rtot_dot_k

Figure B.1: (Continued...) Distributions of the input variables used to train the boosted
decision tree. Trivial bins have been truncated to emphasize shape differencein the non-trivial range. Filled blue histograms indicate the p p → K + K +signal Monte Carlo. Hatched red histograms indicate the atmospheric neutrino background Monte Carlo.

171

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

Input variable: total_rtot_mag

0.00160.00140.00120.001
0.00080.00060.00040.0002
0

1000

2000

3000

4000

5000

total_rtot_mag

Figure B.1: (Continued...) Distributions of the input variables used to train the boosted
decision tree. Trivial bins have been truncated to emphasize shape differencein the non-trivial range. Filled blue histograms indicate the p p → K + K +signal Monte Carlo. Hatched red histograms indicate the atmospheric neutrino background Monte Carlo.

B.6

Variable Ranking

A ranking of the input variables used to train the boosted decision tree can be
estimated by counting how often each is used. Each time a variable is used at a
branching node, it is weighted by the square of the separation gain achieved by the
splitting, as well as the number of events in the node. The variable ranking for thisanalysis is shown in Table B.1.
The most important variables, total rtot mag and total rtot dot k, are measures of the symmetry of the overall light pattern of the event. The dinucleon decay
events have a very high degree of symmetry, largely due to the back-to-back kaons.
The atmospheric neutrino interactions are less likely to create very symmetric distributions of light due to the initial momentum of the incoming neutrino.
The least important variables, rec pi0 mom (the reconstructed π 0 momentum)
and k cand cang[1] (the reconstructed Cherenkov angle of the second kaon candidate), are both members of collections of variables that describe a pair of clas-

172

sified kaon candidates or gamma candidates, respectively. Thus, after the othervariables in their respective collections have already been exploited, there is little
discrimination power left in these two. k cand cang[0] is much more powerful than
k cand cang[1] because every event that enters the boosted decision tree is required
to have at least one kaon candidate, as dictated by the precuts in Section 9.3, whereas
not every event has a second kaon candidate. In fact, very few background events
have two kaon candidates.

173

Rank
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37

Variable Name
total rtot mag
total rtot dot k
g cand impact[0]rec pi0 mass
mu cand mom[1]off vtx sep
dcy e vtxsep
k cand mom[0]k cand cang[0]mu cand mom[0]k cand mom[1]g cand cang[0]mu cand cang[0]g cand vtxsep
g cand totmom
decay vtxsep
g cand show[1]rec dinuc mass
k cand vtxsep
k cand show[1]agood
k cand show[0]g cand cang[1]g cand mom[0]mu cand impact[0]rec pi0 impact
g cand show[0]mu cand show[0]mu cand cang[1]g cand impact[1]rec dinuc mom
mu cand impact[1]mu cand show[1]g cand mom[1]n dcy erec pi0 mom
k cand cang[1]Importance ×102
7.0
5.4
5.0
4.7
4.7
4.6
4.6
4.7
4.5
4.1
3.9
3.6
3.5
3.3
2.7
2.6
2.5
2.4
2.4
2.3
2.2
2.1
2.0
1.8
1.7
1.5
1.4
1.3
1.1
1.1
1.1
0.9
0.9
0.8
0.7
0.6
0.1

Table B.1: Input variables for the boosted decision tree.

Appendix C

Limit Calculations
C.1

Lower Lifetime Limit on Decay Mode

The following expression for the probability of detecting n signal events is obtainedby simple Poisson statistics:

P(n|Γλǫb) =e−(Γλǫ+b) (Γλǫ + b)n
,
n!

(C.1)

where Γ is the decay rate, λ is the exposure, ǫ is the detection efficiency, b is theamount of background, and P(A|B) is the conditional probability of A, given that
proposition B is true.
To be able to take into account the systematic errors in the limit calculation,
Bayes’ Theorem is used. Applying Bayes’ Theorem to the Poisson function gives thefollowing:P(Γλǫb|n) =1
P(n|Γλǫb)P(Γλǫb),
A

(C.2)

where A is a normalization factor.
As the decay rate, detection efficiency, exposure, and background are all indepen174

175

dent variables, the above expression can be rewritten to give:P(Γλǫb|n) =1
P(n|Γλǫb)P(Γ)P(ǫ)P(λ)P(b).
A

(C.3)

Now the probability density function of Γ can be defined as:P(Γ|n) =Z Z Z

P(Γǫλb|n)dǫdλdb,(C.4)

e−(Γλǫ+b) (Γλǫ + b)n
P(Γ)P(ǫ)P(λ)P(b)dǫdλdb.
n!

(C.5)

which can be re-expressed as:1
P(Γ|n) =A

Z Z Z

The normalization factor A is given by:

A=

Z ∞

P(Γ|n)dΓ.

(C.6)

0

The systematic uncertainties on the detection efficiency, exposure, and backgroundcan now be incorporated into the prior distributions, P(λ), P(ǫ) and P(b), the shapesof which are assumed to be truncated Gaussian distributions:

176

P(λ) ∝

P(ǫ) ∝

P(b) ∝

 exp − (λ−λ20 )2(λ > 0)
2σ

,

λ


 0
(λ ≤ 0) exp − (ǫ−ǫ20 )2(ǫ > 0)2σǫ 0(ǫ ≤ 0) exp − (b−b20 )2(b > 0)
2σ
b


 0

(C.7)

,

(C.8)

.

(C.9)

(b ≤ 0)where λ0 , ǫ0 , and b0 are the estimations of the efficiency, exposure, and background,and σλ , σǫ , and σb are their respective systematic uncertainties.
For asymmetric systematic errors, the prior distribution is defined as:(α−α0 )2
(α > α0 )
exp
−
2

2σα+




(α−α0 )2P(α) ∝(0 < α ≤ α0 ) ,
exp
−
2
2σα−




 0
(α ≤ 0)(C.10)

where α is the quantity described by the prior, and σα+ and σα− are its corresponding
positive and negative systematic uncertainties.
Because the expected background in this study is near zero, the statistical error of
the background estimation is taken into account by treating the prior as a convolutionof a Poisson and a Gaussian distribution:P(b) =Z ∞0

e−bMC (bMC )nb(bC − bMC )2
exp −
nb !
2σb2

!

dbMC ,(C.11)

where nb is the number of background events in the unnormalized Monte Carlo, bMC
is the number of true background events in the livetime normalized Monte Carlo, andC is the oversampling factor of the Monte Carlo.

177

The prior distribution of the decay rate, Γ, is taken to be uniform from 0 to acutoff value: 1 (0 < Γ < Γcut )P(Γ) =,
 0 (Γ ≤ 0 or Γ ≥ Γcut )(C.12)

where the cutoff value, Γcut , is the upper limit of the decay rate for the calculationof the normalization factor A needed to avoid divergence. Γcut is set to 10−31 years−1 ,which should be sufficiently large for this study.
The confidence level (CL) can be calculated by integrating the probability density
function:CL =Z Γlimit

P(Γ|n)dΓ.

(C.13)

0

And finally, after setting the confidence level to the desired value and solving theabove equation for Γlimit , one can obtain an expression for the lower lifetime limit:τlimit = 1/Γlimit.

C.2(C.14)

Upper Limit on R-Parity Violating Parameter

Fermi’s Golden Rule gives the expression for the average rate of dinucleon decayΓ̄ as:1Γ̄ = 32π ρN

Z

d3 k1 d3 k2 ρ(k1 )ρ(k2 )vrel (1 − v1 · v2 )σtot (NN → X),(C.15)

where ρN is the average nucleon density, ki is the momentum of the ith nucleon, ρ(ki )178

is the density in momentum space of the ith nucleon, vrel is the relative velocity of
the nucleons, vi is the velocity vector of the ith nucleon (taken to be small), and σtot
is the total cross-section for the decay. The cross-section can be approximated by:

σtot (NN → X) ∼ ρN

128παs2|λ′′uds |4 R10,
vrel M2N(C.16)

where αs is the strong coupling, λ′′uds is the R−Parity prohibited SUSY coupling
constant, MN is the nucleon mass, and R is the ratio between the hadronic and SUSY
scales, given by:

R=

Λ̃
,
(Mg̃ Mq̃4 )1/5

(C.17)

where Λ̃ is the QCD scale, Mg̃ is the gluino mass, Mq̃ is the squark mass.
The partial lifetime τ is simply the inverse of the decay rate:
M2N
1
τ= ∼
R−10 .
′′
2
4
ρN 128παs |λuds |Γ̄(C.18)

Here it clearly shown that the lifetime scales with |λ′′uds |−4 .
Taking ρN = 0.25 fm−3 , αs ∼ 0.12, and MN = 0.94 GeV/c2 , and assuming that τ
represents an upper limit on the decay process lifetime gives the upper limit on λ′′uds
as:|λ′′uds | <

1.65 × 10−30 −10Rτ [years]

1/4.

(C.19)

If R is taken to be ∼ 1.1 × 10−3 , as per Goity and Sher [22], then the upper limiton λ′′uds as a function of the lower lifetime limit on dinucleon decay into kaons is given
by:

179

|λ′′uds | <0.64
τ [years]

1/4.

(C.20)

Bibliography[1] S. Fukuda, et al. Nuclear Inst. and Methods in Physics Research, A 501, 2-3
418 (2003)[2] K. Hirata, et al. Physics Letters B 220, 1 (1989)[3] C. McGrew, et al. Physical Review D 59 052004 (1999)[4] Nishino, et al. Physical Review Letters 102, 141801 (2009)[5] K. Kobayashi, et al. Physical Review D 72 052007 (2005)[6] M. Shiozawa, et al. Physical Review Letters 81, 16 (1998)[7] Y. Hayato, et al. Physical Review Letters 83 1529 (1999)[8] H. Georgi, S. L. Glashow. Physical Review Letters 32 438 (1974)[9] P. Langacker. Physics Reports 72C (1981)[10] P. Langacker. Arxiv preprint hep-ph hep–ph/9411247 (1994)[11] Super-Kamiokande. http://www-sk.icrr.u-tokyo.ac.jp/whatsnew/new-20091125e.html (2010)[12] C. Amsler, et al. Physics Letters B 667, 1 (2008)
180

181

[13] R. Bernabei, et al. Physics Letters B 493, 12 (2000)[14] H. Back, et al. Physics Letters B 563, 23 (2003)[15] C. Berger, et al. Physics Letters B 269 227 (1991)[16] T. Lee, C. Yang. Physical Review 104 254 (1956)[17] C. Wu, et al. Physical Review 105 1413 (1957)[18] A. Sakharov. Physics-Uspekhi 34, 5 392 (1991)[19] L. Carpenter, D. Kaplan, E. Rhee. Arxiv preprint hep-ph/0607204 (2006)[20] A. Barbieri, A. Masiero. Nuclear Physics B 267 679 (1986)[21] S. Dimopoulos, L. Hall. Physics Letters B 196, 2 (1987)[22] J. Goity, M. Sher. Physics Letters B 346, 1-2 69 (1995)[23] K. Abe, et al. ArXiv preprint hep-ex/0607059v2 (2006)[24] G. ’t Hooft. Physical Review Letters 37, 8 (1976)[25] S. Martin. Arxiv preprint hep-ph/9709356 (1997)[26] J. Rosiek. Physical Review D 41, 11 3464 (1990)[27] M. Chemtob. Progress in Particle and Nuclear Physics 54, 1 71 (2005)[28] R. Barbier, et al. Arxiv preprint hep-ph/0406039 (2004)[29] H. Nishino. Ph.D. Thesis, University of Tokyo (2009)[30] M. Honda, et al. Physical Review D 70, 4 (2004)[31] W. G. Webber, S. Stephens. Proc. 20th ICRC (Moscow, 1987) (1987)

182[32] E. Seo, J. Ormes, R. Streitmatter. The Astrophysical Journal 378 7365 (1991)[33] P. Papini, C. Grimani. International Cosmic Ray Conference (1993)[34] M. Boezio, et al. Physical Review Letters 82 4757 (1999)[35] M. Hof, et al. Nuclear Instruments and Methods in Physics Research A 454, 1
180 (2000)[36] T. Sanuki, et al. Arxiv preprint astro-ph/0002481 (2000)[37] A. Alcaraz, et al. Physics Letters B 490 27 (2000)[38] M. Ryan, J. Ormes, V. Balasubrahmanyan. Physical Review Letters 28 985(1972)[39] K. Asakimori, et al. The Astrophysical Journal 28 985 (1998)[40] I. Ivanenko, V. Shestoperov. International Cosmic Ray Conference 2 (1993)[41] Y. Kawamura, et al. Physical Review D 40 729 (1989)[42] A. V. Apanasenko, et al. Astroparticle Physics Journal 16 13 (2001)[43] K. Hänssgen, J. Ranft. Computational Physics Communications 39, 11 53 (1986)[44] S. Roesler, R. Engel, J. Ranft. Physical Review D 57 2889 (1998)[45] G. D. Barr, et al. Physical Review D D70 023006 (2004)[46] G. Battistoni, et al. Arxiv preprint hep-ph/0305208v1 (2003)[47] Y. Hayato. Nuclear Physics Proceedings Supplements 112 171 (2002)[48] C. H. L. Smith. Physics Reptorts 3 261 (1972)

183[49] S. Barish, et al. Physical Review D 16 3103 (1977)[50] S. Bonetti, et al. Il Nuovo Cimento A (1977)[51] M. Pohl, et al. Il Nuovo Cimento A 26, 11 (1979)[52] N. Armenise, et al. Nuclear Physics B 152 365 (1979)[53] A. Vovenko, A. Volkov, V. Zhigunov. Soviet Journal of Nuclear Physics (1979)[54] S. Belikov, A. Bygorsky, L. Klimenko. Zeitschrift für Physik A (1985)[55] J. Brunner, H. Grabosch, H. Kaufmann. Zeitschrift für Physik C (1990)[56] D. Rein, L. M. Sehgal. Annual Physics 133 79 (1981)[57] Y. Ashie, et al. Physical Review D 71, 11 112005 (2005)[58] M. Nakahata, et al. Journal of the Physical Society of Japan 55 3786 (1986)[59] T. Sjostrand. Computer Physics Communications 82 74 (1994)[60] M. Glück, E. Reya, A. Vogt. The European Physical Journal C-Particles 5 461(1998)[61] A. Bodek, U. K. Yang. Arxiv preprint hep-ex/0301036 (2003)[62] M. Derrick, et al. Physical Review D17 1 (1978)[63] S. Barlag, et al. Zeitschrift fur Physik C11 283 (1982)[64] P. Musset, J. P. Vialle. Physics Reports 39 1 (1978)[65] J. E. Kim, et al. Reviews of Modern Physics 53 211 (1981)[66] G. Rowe, M. Salomon, R. H. Landau. Physical Review C18 584 (1978)

184[67] M. Goossens, C. Geneva. cdsweb.cern.ch (1993)[68] T. A. Gabriel, J. E. Brau, B. L. Bishop. IEEE Transactions on Nuclear Science36 14 (1989)[69] C. Zeitnitz, T. A. Gabriel. Nuclear Instruments and Methods in Physics Research
A349 106 (1994)[70] . H.-E. R. A. Group, E. Bracci. cdsweb.cern.ch (1972)[71] K. Nakamura, S. Hiramatsu, T. Kamae. Nuclear Physics A 268, 3 381 (1976)[72] H. Ejiri. Physical Review C 48, 3 1442 (1993)[73] L. Salcedo, et al. Nuclear Physics A 484, 3-4 557 (1988)[74] A. Hocker, et al. Arxiv preprint physics/0703039 (2007)[75] V. Abazov, et al. Physical review letters 98, 18 181802 (2007)[76] H. Yang, B. Roe, J. Zhu. Nuclear Instruments and Methods in Physics Research
A 555 370 (2005)[77] H. Yang, B. Roe, J. Zhu. Nuclear Instruments and Methods in Physics Research
A 574 342 (2007)

Curriculum Vitae
Michael D. Litos
590 Commonwealth Avenue
Department of Physics, Boston University
Boston, MA 02215, U.S.A.

Mobile: (617)224-7377
Office: (617)353-6040
E-mail: mlitos@bu.edu

EDUCATION• Boston University, Boston, Massachusetts– PhD in Experimental Particle Physics, expected in May 2010
• Michigan State University, East Lansing, Michigan– Bachelor of Science in Physics, June 2003
– Additional Major in East Asian Languages & Culture (Japanese), June
2003
RESEARCH EXPERIENCE• Super-Kamiokande / T2K– Hardware:∗ Design, construction, installation and interfacing software for slow control and monitor system for front end electronics∗ Installation of PMTs, wavelength shifter plates, reflective Tyvek wallcoverings∗ Testing of ethernet communication daughter cards for front end electronics upgrade∗ Testing of prototype Hamamatsu PMTs in a water Cherenkov detector– Simulation:∗ Creation of detector simulation in Geant4, including ROOT-based
event display∗ Detector design optimization through MC modeling and testing∗ PMT configuration optimization through MC modeling and testing∗ Creation of custom MC input vector generator for dinucleon decay∗ Update and maintenance of MC generator functionality validation tool

– Analysis:

185

186

∗ First ever search for dinucleon decay into kaons∗ First ever multiple vertex fitting algorithm for a water Cherenkov detector∗ First Super-K analysis using a Boosted Decision Tree
∗ First Super-K analysis incorporating event topology for PID∗ Upgrade of Cherenkov ring fitting software to fit rings produced bycharged kaons
TEACHING EXPERIENCE• Boston University– Discussion session teaching assistant for college physics, Summer 2005
• Michigan State University– Help room teaching assistant for college physics, academic year 2002-2003
PRESENTATIONS• A Search for Dinucleon Decay into Charged KaonsMichael Litos
APS April Meeting, Parallel Talk
St. Louis, MO, April 2008
• A Search for Dinucleon Decay into Kaons Using a Large Water Cherenkov
Detector
Michael Litos
Boston University, Departmental Seminar
Boston University, Boston, MA, November 2009
• Search For Dinucleon Decay Into Kaon Modes Using Multivariate
Techniques
Michael Litos
APS April Meeting, Parallel Talk
Washington D.C., February 2010
• Looking for Double Proton Decay at Super-Kamiokande
Michael Litos
Fermilab, Special Seminar
Fermi National Accelerator Laboratory, Batavia, IL, March 2010
• Looking for Double Proton Decay at Super-Kamiokande
Michael Litos
SLAC, Special Seminar
SLAC National Accelerator Laboratory, Menlo Park, CA, March 2010

187

• Looking for Double Proton Decay at Super-Kamiokande
Michael Litos
LBNL, Special Seminar
Lawrence Berkeley National Laboratory, Berkeley, CA, March 2010
SELECTED PUBLICATIONS• Search for Dinucleon Decay into Kaons via 16 O(pp) → 14 CK + K + in a
Large Water Cherenkov Detector
M. Litos et al. [Super-Kamiokande Collaboration]
In progress• A Proposal for a Detector 2 km Away From the T2K Neutrino Source
E. Kearns et al. [T2K Collaboration]
Proposed to the T2K Collaboration (2005)
OTHER PUBLICATIONS• Search for Proton Decay via p → e+ π 0 and p → µ+ π 0 in a Large Water
Cherenkov Detector
H. Nishino et al. [Super-Kamiokande Collaboration]
Phys. Rev. Lett. 102:141801 (2009)• Kinematic reconstruction of atmospheric neutrino events in a large
water Cherenkov detector with proton identification
M. Fechner et al. [Super-Kamiokande Collaboration]
Phys. Rev. D 79:112010 (2009)• First Study of Neutron Tagging with a Water Cherenkov Detector
H. Watanabe et al. [Super-Kamiokande Collaboration]
arXiv [hep-ex] 0811:0735 (2008)• Solar neutrino measurements in Super-Kamiokande-II
P. Cravens et al. [Super-Kamiokande Collaboration]
arXiv [hep-ex] 0803:4312 (2008)• Search for Matter-Dependent Atmospheric Neutrino Oscillations inSuper-Kamiokande
K. Abe et al. [Super-Kamiokande Collaboration]
Phys. Rev. D 77:052001 (2008)• Search for Supernova Neutrino Bursts at Super-Kamiokande
M. Ikeda et al. [Super-Kamiokande Collaboration]
Astrophys. J. 669:519-524 (2007)

188• Search for neutral Q-balls in Super-Kamiokande II
Y. Takenaga et al. [Super-Kamiokande Collaboration]
Phys. Rev. Lett. B. 647:18-22 (2006)• A Measurement of atmospheric neutrino flux consistent with tau neutrino appearance
K. Abe et al. [Super-Kamiokande Collaboration]
Phys. Rev. Lett. 97:171801 (2006)• High energy neutrino astronomy using upward-going muons in SuperKamiokande-I
K. Abe et al. [Super-Kamiokande Collaboration]
Astrophys. J. 652:198 (2006)• Search for Diffuse Astrophysical Neutrino Flux Using Ultrahigh Energy Upward-Going Muons in Super-Kamiokande I
Molly E. C. Swanson et al. [Super-Kamiokande Collaboration]
Astrophys. J. 652:206-215 (2006)
REFERENCES• Prof. Edward Kearns
Boston University
Department of Physics
590 Commonwealth Avenue
Boston, MA 02215
Office Phone: (617)353-3425
E-mail: kearns@bu.edu• Prof. Christopher Walter
Duke University
Department of Physics
Physics Bldg., Science Dr.
Box 90305, Duke University
Durham, NC 27708
Office Phone: (919)660-2535
E-mail: chris.walter@duke.edu• Prof. Henry Sobel
University of California, Irvine
Department of Physics & Astronomy
4129 Frederick Reines Hall
Irvine, CA 92697-4575
Office Phone: (949)824-6431
E-mail: hsobel@uci.edu

189• Prof. Takaaki Kajita
Institute for Cosmic Ray Research
University of Tokyo
5-1-5 Kashiwa-no-Ha, Kashiwa City
Chiba, 277-8582, Japan
Office Phone: (81)471-36-3102
E-mail: kajita@icrr.u-tokyo.ac.jp