PHYSICAL REVIEW D 91, 072009 (2015)

Search for dinucleon decay into pions at Super-Kamiokande
J. Gustafson,4 K. Abe,1,29 Y. Haga,1 Y. Hayato,1,29 M. Ikeda,1 K. Iyogi,1 J. Kameda,1,29 Y. Kishimoto,1,29 M. Miura,1,29
S. Moriyama,1,29 M. Nakahata,1,29 T. Nakajima,1 Y. Nakano,1 S. Nakayama,1,29 A. Orii,1 H. Sekiya,1,29 M. Shiozawa,1,29
A. Takeda,1,29 H. Tanaka,1 T. Tomura,1,29 R. A. Wendell,1,29 T. Irvine,2 T. Kajita,2,29 I. Kametani,2 K. Kaneyuki,2,29,*
Y. Nishimura,2 E. Richard,2 K. Okumura,2,29 L. Labarga,3 P. Fernandez,3 S. Berkman,5 H. A. Tanaka,5 S. Tobayama,5
E. Kearns,4,29 J. L. Raaf,4 J. L. Stone,4,29 L. R. Sulak,4 M. Goldhaber,6,* G. Carminati,7 W. R. Kropp,7 S. Mine,7P. Weatherly,7 A. Renshaw,7 M. B. Smy,7,29 H. W. Sobel,7,29 V. Takhistov,7 K. S. Ganezer,8 B. L. Hartfiel,8 J. Hill,8N. Hong,9 J. Y. Kim,9 I. T. Lim,9 T. Akiri,10 A. Himmel,10 K. Scholberg,10,29 C. W. Walter,10,29 T. Wongjirad,10 T. Ishizuka,11S. Tasaka,12 J. S. Jang,13 J. G. Learned,14 S. Matsuno,14 S. N. Smith,14 T. Hasegawa,15 T. Ishida,15 T. Ishii,15 T. Kobayashi,15
T. Nakadaira,15 K. Nakamura,15,29 Y. Oyama,15 K. Sakashita,15 T. Sekiguchi,15 T. Tsukamoto,15 A. T. Suzuki,16
Y. Takeuchi,16,29 T. Yano,16 S. Hirota,17 K. Huang,17 K. Ieki,17 T. Kikawa,17 A. Minamino,17 T. Nakaya,17,29 K. Suzuki,17
S. Takahashi,17 Y. Fukuda,18 K. Choi,19 Y. Itow,19 G. Mitsuka,19 T. Suzuki,19 P. Mijakowski,34 J. Hignight,20 J. Imber,20
C. K. Jung,20 J. L. Palomino,20 C. Yanagisawa,20 H. Ishino,21 T. Kayano,21 A. Kibayashi,21 Y. Koshio,21 T. Mori,21
M. Sakuda,21 Y. Kuno,22 R. Tacik,23,31 S. B. Kim,24 H. Okazawa,25 Y. Choi,26 K. Nishijima,27 M. Koshiba,28 Y. Suda,28
Y. Totsuka,28,* M. Yokoyama,28,29 C. Bronner,29 K. Martens,29 Ll. Marti,29 Y. Suzuki,29 M. R. Vagins,29,7 J. F. Martin,30
P. de Perio,30 A. Konaka,31 M. J. Wilking,31 S. Chen,32 Y. Zhang,32 and R. J. Wilkes33(Super-Kamiokande Collaboration)
1

Kamioka Observatory, Institute for Cosmic Ray Research, University of Tokyo,
Kamioka, Gifu 506-1205, Japan
2
Research Center for Cosmic Neutrinos, Institute for Cosmic Ray Research, University of Tokyo, Kashiwa,
Chiba 277-8582, Japan3
Department of Theoretical Physics, University Autonoma Madrid, 28049 Madrid, Spain
4
Department of Physics, Boston University, Boston, Massachusetts 02215, USA
5
Department of Physics and Astronomy, University of British Columbia, Vancouver,
British Columbia V6T1Z4, Canada
6
Physics Department, Brookhaven National Laboratory, Upton, New York 11973, USA
7
Department of Physics and Astronomy, University of California, Irvine, California 92697-4575, USA
8
Department of Physics, California State University, Dominguez Hills, Carson, California 90747, USA
9
Department of Physics, Chonnam National University, Kwangju 500-757, Korea
10
Department of Physics, Duke University, Durham North Carolina 27708, USA
11
Junior College, Fukuoka Institute of Technology, Fukuoka, Fukuoka 811-0295, Japan
12
Department of Physics, Gifu University, Gifu, Gifu 501-1193, Japan
13
GIST College, Gwangju Institute of Science and Technology, Gwangju 500-712, Korea
14
Department of Physics and Astronomy, University of Hawaii, Honolulu, Hawaii 96822, USA
15High Energy Accelerator Research Organization (KEK), Tsukuba, Ibaraki 305-0801, Japan
16
Department of Physics, Kobe University, Kobe, Hyogo 657-8501, Japan17
Department of Physics, Kyoto University, Kyoto, Kyoto 606-8502, Japan18
Department of Physics, Miyagi University of Education, Sendai, Miyagi 980-0845, Japan
19
Solar Terrestrial Environment Laboratory, Nagoya University, Nagoya, Aichi 464-8602, Japan20
Department of Physics and Astronomy, State University of New York at Stony Brook,
New York 11794-3800, USA
21
Department of Physics, Okayama University, Okayama, Okayama 700-8530, Japan
22
Department of Physics, Osaka University, Toyonaka, Osaka 560-0043, Japan
23
Department of Physics, University of Regina, 3737 Wascana Parkway, Regina,
Saskatchewan S4SOA2, Canada
24
Department of Physics, Seoul National University, Seoul 151-742, Korea
25
Department of Informatics in Social Welfare, Shizuoka University of Welfare,
Yaizu, Shizuoka 425-8611, Japan
26
Department of Physics, Sungkyunkwan University, Suwon 440-746, Korea
27
Department of Physics, Tokai University, Hiratsuka, Kanagawa 259-1292, Japan
28
The University of Tokyo, Bunkyo, Tokyo 113-0033, Japan
29
Kavli Institute for the Physics and Mathematics of the Universe (WPI), Todai Institutes for Advanced
Study, University of Tokyo, Kashiwa, Chiba 277-8582, Japan30
Department of Physics, University of Toronto, 60 St., Toronto, Ontario M5S1A7, Canada
31
TRIUMF, 4004 Wesbrook Mall, Vancouver, British Columbia V6T2A3, Canada

1550-7998=2015=91(7)=072009(20)072009-1© 2015 American Physical Society

J. GUSTAFSON et al.

PHYSICAL REVIEW D 91, 072009 (2015)32

Department of Engineering Physics, Tsinghua University, Beijing 100084, China
Department of Physics, University of Washington, Seattle, Washington 98195-1560, USA
34
National Centre For Nuclear Research, 00-681 Warsaw, Poland(Received 17 February 2015; published 28 April 2015)33

A search for dinucleon decay into pions with the Super-Kamiokande detector has been performed
with an exposure of 282.1 kiloton-years. Dinucleon decay is a process that violates baryon number by
two units. We present the first search for dinucleon decay to pions in a large water-Cherenkov detector.
The modes 16 OðppÞ → 14 Cπ þ π þ , 16 OðpnÞ → 14 Nπ þ π 0 , and 16 OðnnÞ → 14 Oπ 0 π 0 are investigated.
No significant excess in the Super-Kamiokande data has been found, so a lower limit on the lifetimeof the process per oxygen nucleus is determined. These limits are τpp→πþ πþ > 7.22 × 1031 years,
τpn→πþ π0 > 1.70 × 1032 years, and τnn→π0 π 0 > 4.04 × 1032 years. The lower limits on each mode areabout 2 orders of magnitude better than previous limits from searches for dinucleon decay in iron.
DOI: 10.1103/PhysRevD.91.072009

PACS numbers: 13.30.Eg, 14.20.Dh, 29.40.Ka

I. INTRODUCTION
Baryon number seems to be a conserved quantity inorder to explain the stability of matter, yet baryon number
must be violated to explain the predominance of matterover antimatter in the Universe [1]. Single nucleon decay
processes that violate baryon number by one unit but
conserve ΔjB − Lj (where B is baryon number and L is
lepton number) have been studied extensively at SuperKamiokande [2–6] in the context of grand unified theories(GUTs) that predict an unstable proton. Lower limits on
proton decay modes such as p → eþ π 0 of 1.2 × 1034 years[7] have ruled out the simplest GUTs, such as minimal
SU(5) [8]. Other processes, such as those that violate B by
two units, are of great theoretical interest although notnecessarily predicted by GUT-inspired models. Searches
for such processes constitute an important subset ofsearches for baryon number violation, both in their own
right and given the continued nonobservation of single
nucleon decay. The discovery of such a process would be aspectacular revelation of new physics, while nonobservation provides important constraints on the theories which
predict them. Dinucleon decay—the simultaneous decay of
two bound nucleons into leptons or mesons—constitutes
such a process. In this paper we search for the processes
16OðppÞ→ 14 Cπ þ π þ , 16 OðpnÞ→ 14 Nπ þ π 0 , and 16 OðnnÞ →14
Oπ 0 π 0 (abbreviated throughout to pp → π þ π þ ,pn → π þ π 0 , and nn → π 0 π 0 ).
Processes violating baryon number by two units canoccur in a broad class of models [9]. The experimentally
accessible processes of this kind are dinucleon decay andneutron-antineutron oscillation (n-nbar). The latter process
may also be pursued in experiments with free neutrons [10],whereas dinucleon decay requires the close overlap ofbound nucleons. The tree-level Feynman diagram for
n-nbar and dinucleon decay to pions is the same, with
two spectator quarks in the latter case. An example of a*

Deceased.

dinucleon decay diagram from one of the models discussed
in Ref. [9] is shown in Fig. 1. The expected time scale for
dinucleon decay to pions can be related to the time scale for
n-nbar, if the same operator is dominant for each process.
This is discussed further in the Appendix.
One dinucleon decay search has been performed
before at Super-Kamiokande, namely a search for decay
to charged kaons [11]. The Frejus experiment, which usediron nuclei, performed an extensive search covering dinucleon decay final states [12]. The search in Ref. [12] included
the modes analyzed in this paper. The more exclusive limits
in this paper relative to those in Ref. [12] reflect the large
exposure of Super-Kamiokande. Our limit calculationaccounts for systematic uncertainties, particularly thosedue to intranuclear pion interactions. This is the first searchfor dinucleon decay into pions in a water-Cherenkov
detector.
This paper proceeds as follows. Section II describes
some basic features of the Super-Kamiokande experiment.
Section III describes relevant aspects of our Monte Carlosimulation, for both the dinucleon decay signal andatmospheric neutrino background. Section IV describes
standard data selection and reconstruction algorithms usedin Super-Kamiokande. Section V describes the analysis of
all three modes, first discussing some features of boosted

FIG. 1. Feynman diagram for dinucleon decay to pions fromone of the models discussed in Ref. [9], with spectator quarks
q1 ; q2 ∈ ðu; dÞ. The X 1 and X 2 scalar particles belong toð6̄; 1; − 13Þ and ð6̄; 1; 23Þ representations of the Standard ModelSUð3Þ ⊗ SUð2Þ ⊗ Uð1Þ gauge group.

072009-2

SEARCH FOR DINUCLEON DECAY INTO PIONS AT …

PHYSICAL REVIEW D 91, 072009 (2015)þ þ

decision trees that are relevant for the pp → π π and
pn → π þ π 0 searches. Sections VI and VII describe ourresults and conclusion.

neutrino MC is the same MC sample used in other SK
analyses, as described in detail in Ref. [16].
A. Dinucleon decay

II. THE SUPER-KAMIOKANDE DETECTOR

As dinucleon decay requires the close overlap of two
nucleons, we restrict attention to the oxygen nucleus for the
decaying nucleons. Fermi momentum, correlation between
nucleons, nuclear binding energy, and pion-nucleon interactions are taken into account. This last effect is especiallyprominent, given the dense nuclear material the pions
traverse.
Fermi momentum and nuclear binding energy are takenfrom measured electron-12 C scattering [17]. Nuclear binding energy is taken into account by modifying the nucleon
masses. Reference [18] estimated that 10% of single
nucleons that undergo decay have wave functions thatare correlated with other nucleons in the nucleus; weassume the same level of correlation for dinucleon decay
pairs. Such correlated decay events typically have a lower
total invariant mass, due to the momentum of the correlated
nucleon, which is below the Cherenkov threshold (Fig. 2).
The positions of the decaying nucleons in 16 O are approximated using the Woods-Saxon nuclear density model [19].
Dinucleon decay leaves the remaining nucleus in an excitedstate from which it deexcites by the emission of gamma
rays. However, the deexcitation gamma rays are on theorder of a few MeV, while the pions are on the order of a
GeV, so the extra light coming from the gamma rays is
negligible relative to the pions. The main difference indinucleon decay simulation from single-nucleon decay isthe available phase space. For dinucleon decay, the initial
total momentum and energy are summed over both nucleons (or three nucleons, if there is a correlated decay).
Pion-nucleon interactions in the nucleus are simulatedusing NEUT [20]. The interactions modeled include
quasielastic scattering (in which the incoming and outgoing pion have the same charge, but there is some energyloss), charge exchange (e.g. π þ þ n → π 0 þ p), absorption (e.g. π þ þ n → p), and hadron production (e.g.
104III. SIMULATION
Dinucleon decay signal and atmospheric neutrino background Monte Carlo (MC) simulations are generated toestimate signal efficiency, expected background rates, andsystematic uncertainties. Since each of the four SK datataking periods have different detector configurations, separate sets of MC for both the dinucleon decay signal and theatmospheric neutrino background are used for each period.
The dinucleon decay MC is based on single nucleon decay
MC used for previous SK analyses [2–6]. The atmosphericInvariant mass of pions

103

EventsSuper-Kamiokande (hereafter “SK”) is a large, cylindrical water-Cherenkov detector. Full details of the detector
can be found in Ref. [13]; here, we note only the most basic
features.
The detector contains two optically separated regions:the inner detector (ID) and outer detector (OD). Inwardfacing, 20-inch photomultiplier tubes (PMTs) collectCherenkov light in the ID, while light from penetrating
particles (typically cosmic-ray muons) is detected by
8-inch, outward-facing PMTs in the OD. The fiducial
volume of the ID is defined as the cylindrical volume withsurfaces 2 meters inward from the ID wall. This corresponds to 22.5 ktons, and 7.5 × 1033 oxygen nuclei.
SK has had four data-taking periods. It began taking data
in April 1996, with 11 146 PMTs covering 40% of the ID
surface. The first data-taking period, called SK-I, continued
until July 2001, totaling 1489.2 live-time days, or 91.7
kton-years. An accident occurred in November 2001,
which caused the implosion of about half of the ID
PMTs. The remaining 5182 PMTs were redistributed
uniformly on the surface of the ID, providing 19%
photo-coverage. To prevent such accidents in the future,each PMT has been enclosed in acrylic and fiber-reinforced
plastic cases. The second data-taking period, which used
this decreased photo-coverage, lasted from December 2002
until October 2005, and is called SK-II. This lasted 798.6
live-time days, or 49.2 kton-years. Replacement PMTs
were produced and the photo-coverage was returned to
40% in 2006. The third data-taking period, from July 2006
until September 2008, is called SK-III, comprising518.1 live-time days (31.9 kton-years). The detector was
upgraded in the summer of 2008, with improved electronics
that record all PMT hit information without dead time
[14,15]. This is the current configuration of the detector,
called SK-IV. In this paper, SK-IV data are used up to
March 2014, which comprises 1775.6 live-time days
(109.3 kton-years). The total SK-I–SK-IV data set is282.1 kton-years.

102
10
1
10-1
10-2

500

1000

1500

2000

Invariant mass (MeV/c 2)

FIG. 2 (color online). Invariant mass of the pions frompp → π þ π þ , before nuclear effects. Plots for pn → π þ π 0 and
nn → π 0 π 0 look similar. The dashed histogram shows the
spectrum for events with a correlated nucleon.

072009-3

J. GUSTAFSON et al.
þþ

PHYSICAL REVIEW D 91, 072009 (2015)0π þ p → π þ π þ p). Quasielastic scattering is relevant
at all pion momenta, as is charge exchange, although it has
a lower overall cross section. Absorption has a sizable crosssection in the pion momentum range 100–300 MeV=c, but
it is negligible elsewhere. Hadron production is relevant forpion momentum greater than 500 MeV=c, and is therefore
relevant for most pions immediately after the decay. With
these effects taken into account, less than half of thegenerated events in each mode escape the nucleus in theoriginal pion configuration (Table I).
After the pions leave the nucleus, their propagation
and Cherenkov light emission is passed through a custom
detector simulation based on GEANT3 [21]. Charged pions
are susceptible to further secondary interactions in thewater, while neutral pions decay before such interactions.
The propagation of charged pions through the water below
500 MeV=c is simulated by a custom code based on
Ref. [22], while charged pions above this momentum arehandled by GCALOR [23].
Monte Carlo events are generated in a broadened fiducial
volume (FV), greater than 1 meter from the ID wall as
opposed to 2 meters which defines the FV. This allows us toaccount for events that migrate across the FV boundary
when reconstructed. The final detection efficiency is

TABLE I. Number of pions for pp → π þ π þ (top), pn → π þ π 0(middle), and nn → π 0 π 0 (bottom) after nuclear effects. Percentages indicate the percentage of events with the specified numberof π þ ’s (columns) and π 0 ’s (rows). Entries with a “  ” indicate asmall contribution (< 0.1%). The number of pions prior to
nuclear effects is shown in bold typeface.
pp → π þ π þ
πþ

0

1

2

3


2.7%
2.8%

11.8%
21.8%
3.0%

45.1%
5.2%


2.4%



π0
0
1
2

þ 0

pn → π π
π

þ

0

1

2

3


7.1%
10.5%

6.8%
46.8%
4.1%

12.3%
3.7%


1.4%



π0
0
1
2

nn → π 0 π 0
πþ

0

1

2


11.5%
43.0%

2.6%
25.1%
4.2%

3.5%
2.8%


π0
0
1
2

defined as the ratio of events surviving all selection criteria
to the total number of events generated in the fiducial
volume. For each SK period, we generate a total of 25 000
events. Three times this number of events is generated inmodes for which we use a multivariate search method. Thisis because the multivariate method requires three separate
samples, as described later.
B. Atmospheric neutrinos
The atmospheric neutrino MC sample is based on the
Honda atmospheric neutrino flux [24], and the NEUT
neutrino interaction and nuclear effects simulator [20].
The dominant neutrino interaction background for the
dinucleon decay searches comes from charged-current
single pion production (CC1π; ν þ N → l þ N 0 þ π for
nucleons N; N 0 and charged lepton l), as well as bothcharged- and neutral-current deep inelastic scattering(CCDIS and NCDIS, respectively). The DIS interaction
has the form ν þ N → l þ N 0 þhadrons, where l is acharged lepton (CCDIS) or a neutrino (NCDIS). Pionsin CC1π interactions are produced following the model ofRein and Sehgal [25]. This is the dominant hadronproduction interaction for hadronic invariant mass W
below about 2 GeV=c2 . Deep inelastic scattering is treateddifferently depending on W. For 1.3 GeV=c2 ≤ W ≤2.0 GeV=c2 , only pions are considered as outgoing hadrons,
and the number of pions is two or greater to differentiate itfrom CC1π. For W > 2.0 GeV=c2, other mesons are also
considered, such as η mesons and kaons. After leaving thenucleus, particles are propagated through the water using thesame tools described above for dinucleon decay.
The atmospheric neutrino samples correspond to an
exposure of 500 years for each SK period. Event rates
are weighted to include the effect of neutrino oscillations.
Final event rates after all selection criteria have been
applied are normalized to the corresponding SK-I–SK-IV
detector live time.
IV. REDUCTION AND RECONSTRUCTION
This search uses the fully contained data set. Fully
contained events are defined as events in which all charged
particles start and stop in the ID. Exiting and entering
particles, particularly cosmic-ray muons, are excluded byrequiring minimal activity in the OD. The threshold for a
triggered event in the fully contained data set is 5.7 MeV ofvisible energy in SK-I, -III, and -IV, and 8 MeV of visibleenergy in SK-II. Visible energy is defined as the energy to
produce the observed light in the event if it were produced
by a single electron. The majority of triggered events arecosmic-ray muons and low-energy radioactivity frommaterial surrounding the detector. A dedicated set of event
selection algorithms is applied to eliminate the majority ofthese backgrounds, such that all but a fraction of remaining
events come from atmospheric neutrinos [16].

072009-4

SEARCH FOR DINUCLEON DECAY INTO PIONS AT …

PHYSICAL REVIEW D 91, 072009 (2015)

Events that remain are passed through a reconstruction
program. The reconstruction for each event determines anoverall event vertex and the number of Cherenkov rings, aswell as the direction, momentum, and particle classificationof each ring. Events are assumed to originate from a single
vertex, and the distribution of observed charge is used tofind the first Cherenkov ring. Additional rings are foundusing a Hough transform method [26]. Full details on our
event reconstruction can be found in Refs. [16] and [27].
Each Cherenkov ring is classified as “e-like” or “μ-like,”referring to electrons and muons. This is a conventional
way of denoting whether or not a ring displays thecharacteristics of electromagnetic showers. Showeringrings are produced by electrons and photons, but not
protons, muons, or charged mesons. The classification ofeach ring is determined by its light pattern. Electromagnetic
showers produce fuzzy rings due to scattering and pair
production, while heavier particles produce rings withsharper edges.
The momentum of each ring is determined based on theparticle type and the number of photoelectrons (referred toas charge) within a 70 degree half-angle cone of the ring
direction. The charge is corrected for light attenuation and
PMT angular acceptance. The overall charge in the cone iscorrected for variations in the PMT gains.
In addition to the standard ring-finding program, analgorithm described in detail in Ref. [28] is used to detect
particle tracks. Unlike the standard ring finder, this algorithm does not search for ring edges, but rather tests particle
hypotheses as either electrons or photons (specifically, twophotons from a π 0 decay) using only the light pattern. It has
been used in atmospheric neutrino [29] and long-baseline[30] oscillation analyses to reject neutral-current single π 0
events. It is used in the pn → π þ π 0 search to identify a π 0from the decay and use its reconstructed mass as adiscriminating variable for the analysis.
After the primary event, electrons from stoppingmuons (Michel electrons) are found by searching forclusters of hits detected at the same time. The time
window for such clusters extends to 20 μs after the main
event. In the SK-I, -II, and -III data-taking periods, there
was an impedance mismatch in the electronics whichcaused signal reflection around 1000 ns after the main
event. Thus, the time period 800–1200 ns after the main
event was excluded from the Michel electron search. The
improved electronics of SK-IV eliminated this signal
reflection, and as a result the tagging efficiency was
improved from 80% to 96% for μþ decays and 63% to
83% for μ− decays between SK-I–SK-III and SK-IV.
Hence, event selections that place an upper limit on thenumber of Michel electrons reduce the event rates morein SK-IV than in SK-I–SK-III.
After the reconstruction, a final set of selection criteria is
applied to isolate fully contained fiducial volume (FCFV)
events:(i) No significant OD activity: the largest OD hit
cluster contains less than ten hits (SK-I) or 16 hits
(SK-II–SK-IV).
(ii) The total visible energy is greater than 30 MeV.
(iii) The reconstructed event vertex must be at least 2 mfrom the inner detector wall. This is the fiducial
volume cut as defined in Sec. II.
After the FCFV selection, the fraction of events comingfrom sources other than atmospheric neutrinos is estimatedto be < 1%. In the next section, all analyses start with the
FCFV events.
V. ANALYSISThe experimental signature of a dinucleon decay to pions
event, for all three modes considered, is two back-to-backpions, each forming at least one Cherenkov ring of thecorrect type. The characterization of such events is madedifficult due to the interactions of charged pions in both theresidual oxygen nucleus and the water. Interactions in thewater lead to truncated Cherenkov light signatures, which
typically lead to poor momentum reconstruction. Thus,
for modes involving at least one π þ (pp → π þ π þ andpn → π þ π 0 ), simple selection criteria based on kinematic
variables (such as the total momentum and invariant massthat are used in other nucleon decay analyses) do not give
optimal sensitivity. For these modes, a multivariate method
is applied to statistically identify events with some, but
perhaps not all, of the expected characteristics that differentiate signal and background. Specifically, a boosted
decision tree (BDT) method is used. The method isimplemented in the ROOT-based TMVA [31] analysis library.
This was also used in the dinucleon decay analysis
in Ref. [11].
For nn → π 0 π 0, the situation is more straightforward, as
the π 0 ’s that escape the nucleus decay to energetic photons
before interaction in the water, and the electromagnetic
showers are well reconstructed. For this mode, simple
kinematic selection criteria are sufficient for the analysis.
A. Boosted decision trees
Here we describe some relevant details of the BDT,
applicable for both pp → π þ π þ and pn → π þ π 0 searches.
Similar to other multivariate methods, a BDT takes a setof input variables—in this case, reconstructed information
for both the dinucleon decay signal and atmospheric
neutrino background—with the goal of producing a single
discriminating output variable that separates signal andbackground. To get to this single variable, three stages are
required. They are referred to as training, testing, andanalysis, and each uses separate sets of signal and background MC. We describe the process briefly.
Boosted decision trees are sets of individual decision
trees, which are simple binary trees formed by selection
cuts on single variables. Trees are built using the training

072009-5

J. GUSTAFSON et al.

PHYSICAL REVIEW D 91, 072009 (2015)

sample as follows. Starting with the full training sample, aselection cut on a single variable is made that best separates
signal and background. For the analyses in this paper, thisoptimal separation is determined by the Gini index
g ¼ pð1 − pÞ, where p is the purity of the sample after
the selection. Other separation types exist, though they do
not lead to significantly different performance [31]. The
sample then splits off into two samples, or nodes, based on
this selection. Each of the daughter nodes does the samething, separately finding the variable that best separates
between signal and background, and so forth until some
stopping criterion is reached (such as a minimum numberof events remaining on a node). Terminal nodes are
declared signal-like or background-like based on the purityof signal events in the node.
The training process continues over many such trees(typically 1000 or so) using the same training sample,
with the prescription that events in subsequent trees are
reweighted, or “boosted.” Several boosting algorithms
exist, and are described in Ref. [31]. Upon reweighting,the best set of discriminating variables may change.
Variables can also be used more than once in a given tree.
Once the training sample has been processed, additional
data sets may pass through the trained set of trees. On an
event-by-event basis, the final BDT output is given by
Y BDT ¼

N tree1 X
w h ðxÞ:N tree i¼1 i ið1ÞHere N tree is the size of the ensemble of trees, wi is theweight of each tree in the ensemble, determined by the
boosting method, hi is the “performance” of each event for
the ith tree (whether it ended on a signal- or backgroundlike node), and x is the set of input variables. Variables are
given a ranking, or “importance,” which is a numericalvalue that roughly corresponds to the number of times thevariable is used in the training process relative to the othervariables in the set.
After the BDT is trained, two additional, statisticallyindependent sets of signal and background MC are necessary. The test samples are used to determine the optimal
boost method, as well as some additional parameters. The
analysis samples are used to obtain the final sensitivity.
Unlike in the training and testing process, no information is

given to the BDT regarding the composition of the analysis
samples as either signal or background. Both the test andanalysis samples are also important in detecting training
bias in the BDT. This bias is a situation in which the BDT
has been trained to model statistical fluctuations in the
training sample, as opposed to modeling from generaltrends in the input variables. Such overtraining can lead todifferent performance between MC samples, even if they
are statistically similar.
Table II quantifies the signal and background MC
samples used for each stage. Separate trees are trained,
tested, and analyzed for each SK period. Real data areprocessed through the trained BDTs as a final step.
1. pp → πþ πþThe interactions of π þ ’s in the nucleus and the watercomplicate the event signature for pp → π þ π þ. A set of
selection criteria is applied to keep pp → π þ π þ events, andreduce the atmospheric neutrino background. The key
characteristics are that the event has two “back-to-back”
rings, as the ring-finding algorithm is still good at reconstructing the ring direction even with the poor momentum
resolution, and that the rings are identified as μ-like. Due tointeractions in the water, the π þ ’s can simulate showering
effects and therefore be misidentified as e-like. However,the converse—in which rings produced by electrons orphotons are misidentified as μ-like—is rare. The following
selection criteria are applied.
(A1) There is more than one Cherenkov ring.
(A2) The two most energetic rings are μ-like. These areassumed to correspond to the π þ candidates.
(A3) The angle between the two most energetic rings isgreater than 120 degrees, a minimal requirement for“back-to-back” pions.
(A4) The total visible energy is less than 1600 MeV,
approximately the maximum visible energy expectedfrom a dinucleon decay event ð2Mp − 2Mπþ Þ, where
Mp is the proton and M πþ the pion mass.
The signal efficiency, expected background, and data
events remaining after the initial selection criteria areshown in Table III. After applying these criteria, we enter
the remaining events in the BDT. The following set of input
variables is used, ordered by their relative importance, as
described in the BDT section above.

TABLE II. Independent MC sample sizes used for each stage of the BDT for both the pp → π þ π þ and pn →π þ π 0 searches, prior to any event selection. Dinucleon decay sample sizes show only the number of events, while theatmospheric neutrino sample sizes show the approximate number of events and the equivalent number of years. The
same statistics are used for SK-I–SK-IV.

dinuc. decay
atm-ν

Training

Testing

Analysis25000 ev
∼750; 000 ev (150 yr)

25000 ev
∼750; 000 ev (150 yr)

25000 ev
∼1; 000; 000 ev (200 yr)

072009-6

SEARCH FOR DINUCLEON DECAY INTO PIONS AT …

PHYSICAL REVIEW D 91, 072009 (2015)

TABLE III. Efficiencies, expected backgrounds, and data
events after selection criteria are applied in the pp → π þ π þsearch. Background rates are scaled to the appropriate SK period
live time. Statistical errors are shown.

Eff. (%)
Bkg.
data

SK-I

SK-II

SK-III

SK-IV

11.2  0.2
33  0.9

10.5  0.2
17  0.5

12.0  0.2
13  0.4

12.1  0.2
45  1.2

27

14

8

43

(a1) The angle between the two most energetic μ-likerings. It is expected to peak around a perfect“back-to-back” angle of 180 degrees for the signal,and have a falling spectrum for increasing angle inthe background.
(a2) The ratio of charge carried by the highest-energy
ring to the total charge of all rings. This is expected
to be smaller for signal, since the pions are
expected to have roughly equal energy. By contrast,
we expect most of the energy from a chargedcurrent atmospheric neutrino event to be in the
charged lepton ring.
(a3) The total visible energy. This has a broad distribution
up to 1600 MeV for the signal, but peaks at a lower
value for the background around 200 MeV.
(a4) The maximum distance between a Michel electron
vertex and the primary event vertex. We expect thisto be larger for energetic muons than for pions, since
the muons penetrate water without sudden hadronicinteractions. Muons above the Cherenkov thresholdsometimes appear in the signal (about 9% of thetime), but are much more common in the background,
due to charged-current νμ interactions such as single
pion production.
(a5) The maximum angle between any μ-like ring direction vector, and any Michel electron vertex vector (aspointed to from the primary vertex). Similar to (a1),this is designed to look for “back-to-back” events,
and is an additional check on the topology of the
μ-like rings.
(a6) The magnitude of the vector sum of corrected charge
associated with all rings, where the corrections areapplied for various factors, as described in Sec. IV.
This is similar to the total momentum, since thecorrected charge associated with each ring is used todetermine its momentum. Due to typically poor
momentum resolution from hadronic energy losses,
it is preferable to use this variable instead of the total
momentum.
(a7) The number of Cherenkov rings. Though we aretrying to isolate two pions, they typically interact inthe water, often producing visible secondary tracks,
so that the distribution of the number of rings peaks
at three. By contrast, the Cherenkov ring spectrum

for the background falls with increasing ringnumber.
(a8) The number of Michel electrons. Naively we would
expect two Michel electrons from two charged pions.
However, the pions could undergo charge exchangeor absorption, leading to zero or one Michel electron.
Though the signal and background distributions bothpeak at one decay electron, this variable is relatively
independent of most of the other variables, which
depend on Cherenkov-ring counting.
(a9) The number of μ-like rings. Similar to the total
number of rings, this variable specifically seeks
scatters that produce a secondary charged pion.
The relative importance of each input variable is shownin Table IV. Figure 3 shows plots of the two mostdiscriminating variables: the angle between π þ candidates
and the ratio of charge carried by the most energetic ring.
The number of events in the data is small after the selection
criteria is applied, but the data and atmospheric neutrino
MC are consistent.
The final distributions of the BDT output for all
SK periods are shown in Fig. 4. For pp → π þ π þ, thegradient boost [31] method is used. The final efficiencies,
backgrounds, and candidate data events are shown in
Table V. SK-II has a slightly lower efficiency than theother SK periods, as fewer events pass selection
criterion (A2).
The final cut on the BDT output was determined using
only signal and background MC, not the real data. The cut
takes two basic
pﬃﬃﬃ considerations into account: 1) the cut
maximizes ϵ= b, where ϵ is the signal efficiency and b is
the background rate; and 2) the cut is not too close to thesignal peak, since variations are propagated through the
BDT to estimate systematics, and the part of the distribution near the peak changes the most under these variations.
For 2), 0.7 was the closest value to the peak that we
considered. As the efficiency (background) distribution
increases (decreases) nearly monotonically for all SK
periods, it was found that a cut at 0.7 satisfied 1) while
obeying the constraint of 2).

TABLE IV. Relative importance of each variable in the
pp → π þ π þ search, averaged across SK periods.
VariableAngle between μ-like rings
Ratio of charge carried by most energetic ring
Visible energy
Max. distance to Michel vertex
Max. angle between μ-like ring and Michel vertex
Magnitude of vector sum of corrected charge
Number of rings
Number of Michel electrons
Number of μ-like rings

072009-7

Importance
0.16
0.15
0.15
0.13
0.13
0.12
0.071
0.055
0.045

J. GUSTAFSON et al.

PHYSICAL REVIEW D 91, 072009 (2015)

SK-I
7

5

6

4

5

EventsEvents

SK-I
6

3
2

1
0

0

130

140

150

160

170

180

0

0.2

0.4

0.6

Angle (degrees)

Charge ratio

SK-II

SK-II

3

0.8

1

0.8

1

0.8

1

0.8

1

4

EventsEvents

3
2

1

120

4

2

3
2

1

1
0

0

120

130

140
150
160
Angle (degrees)

170

180

0

0.2

0.4

0.6

Charge ratio

SK-III

SK-III

EventsEvents

2

1

0

0

120

2

130

140

150

160

170

180

0

0.2

Angle (degrees)

SK-IV
8
6

EventsEvents

7
5
4
3
2
1
0

120

130

140

150

160

0.4

0.6

Charge ratio

170

180

Angle (degrees)

SK-IV

10
9
8
7
6
5
4
3
2
1
0
0

0.2

0.4

0.6

Charge ratio

FIG. 3 (color online). The angle between μ-like rings and the ratio of charge carried by the most energetic ring, for the pp → π þ π þsearch, for signal MC (solid histogram), atmospheric neutrino MC (dashed histogram), and data (crosses). Distributions are normalized
to the number of events in the data. Details of each variable are described in the text.

072009-8

SEARCH FOR DINUCLEON DECAY INTO PIONS AT …

PHYSICAL REVIEW D 91, 072009 (2015)

SK-II

12

6

10

5

8

4

EventsEvents

SK-I

6

3

4

2

2

1

0

0
-1

-0.5

0

0.5

1

-1

0
BDT Output

SK-III

SK-IV

3.5

0.5

1

0.5

1

16
14

3

Events2.5

Events-0.5

BDT Output

2

1.5

12
10
8

1
0.5
0
-1

-0.5

0

0.5

6
4
2
0

1

-1

-0.5

BDT Output

0
BDT Output

FIG. 4 (color online). BDT output for pp → π þ π þ signal MC (solid), atmospheric neutrino background MC (dashed), anddata (crosses). The vertical green line indicates the BDT cut value, and the arrow indicates that only events to the right of the cut
are kept.

The largest component (∼40%) of the remaining atmospheric background is from charged-current single pion
production (CC1π). In this mode, typically a νμ willproduce a muon and a charged pion at a large angle, thus
passing the π þ candidate angle requirement. Other contributing background interaction modes include CCDIS,NCDIS, charged-current quasielastic (CCQE), and neutralcurrent single pion production (NC1π). Table VI quantifies
the neutrino interaction modes of the remaining background. The percentages agree within statistical fluctuations across SK periods.
A pp → π þ π þ MC event is shown in Fig. 5. An
atmospheric neutrino event surviving all selection criteria,

including the final BDT cut, is shown in Fig. 6. Twocandidate events were found: one in SK-II (Fig. 7) andanother in SK-IV (Fig. 8). No events were observed in
SK-I or SK-III. Table VII shows values of the BDT input
variables for each candidate, compared to the means of thesignal and background distributions. While both candidate
events appear signal-like in some variables—particularly
the angle between μ-like rings (both candidates) and energy
variables (SK-IV candidate)—the candidates are consistent
with the expected background rate for each SK period(Table V), as well as the total background rate of 4.5
expected across SK-I–SK-IV. Thus, we conclude that there
is no evidence for pp → π þ π þ in the data.

TABLE V. Efficiency, expected background events, and candidate data events for the SK-I–SK-IV pp → π þ π þ search.
Statistical errors are quoted for efficiency and background.
Background is quoted both for the appropriate SK live timeand per megaton-year.

TABLE VI. Atmospheric neutrino interaction modes in theremaining background for the pp → π þ π þ search by percentage,as estimated by atmospheric neutrino MC.

SK-I

SK-II

SK-III

SK-IV

Eff. (%)
6.1  0.2 5.3  0.2 6.4  0.2 5.8  0.2
Bkg. (MT-yr) 17.8  1.8 14.3  1.6 17.4  1.7 14.2  1.6
Bkg. (SK live.)
1.6
0.70
0.56
1.6
Candidates
0
1
0
1

Mode

SK-I

SK-II

SK-III

SK-IV

CC1π
CCDIS
NCDIS
CCQE
NC1π

45  7%30  5%10  4%9  3%
5  3%

42  7%31  6%9  4%
9  3%
8  4%

51  7%24  5%6  3%
10  3%8  3%

42  7%34  7%13  5%10  4%0 þ 1%

072009-9

J. GUSTAFSON et al.

PHYSICAL REVIEW D 91, 072009 (2015)

FIG. 5 (color online). A pp → π þ π þ MC event. The intensity
scale indicates the amount of charge deposited on the PMTs. The
main display shows the ID, while the small display on the top
right shows the OD. The white þ’s indicate the horizontal andvertical location of the reconstructed vertex, projected onto thedetector wall. Orange rings are fit as μ-like, and cyan rings ase-like. Solid rings indicate that the ring’s opening angle is near theultrarelativistic limit in water of θc ¼ cos−1 ð1=nÞ ≈ 41 degrees,
while dashed or dotted rings are more collapsed. The two solidorange rings correspond to true π þ ’s, and have an angle of
164 degrees between them. The dotted cyan ring is fit as e-like,but is actually a hard scatter of one of the π þ ’s.

FIG. 6 (color online). Surviving atmospheric neutrino MC
event for the pp → π þ π þ search. The interaction type is νμcharged-current single π þ production. The uppermost μ-like ringcorresponds to a μ− , while the smaller μ-like ring is a π þ . The
e-like ring is mis-fit, as it comes from a hard scatter of the π þ .
The angle between the two μ-like rings is 170 degrees, thereby
mimicking the main event signature. The momentum of the muon
μ-like ring is 402 MeV=c, and the momentum of the π þ μ-likering is 301.3 MeV=c.

FIG. 7 (color online). pp → π þ π þ candidate event in SK-II.
The two μ-like rings are strikingly back-to-back, at almost exactly
180 degrees to each other. The e-like ring could be due to a
hard scatter. There is one Michel electron, not shown. The μ-likerings have momenta of 399.0 MeV=c (solid) and 280.5 MeV=c(dashed), and the e-like ring has momentum 42.0 MeV=c. The
visible energy is 387.0 MeV, the total momentum 124.2 MeV=c,
and the invariant mass 710.9 MeV=c2 .

2. pn → πþ π0There are many similarities between the pn → π þ π 0
search and the pp → π þ π þ search. The main difference isthat the π 0 provides a new set of discriminatory variables,

FIG. 8 (color online). pp → π þ π þ candidate event in SK-IV.
The two μ-like rings are at an angle of 151 degrees to each other,with momenta of 541.4 and 389.5 MeV=c. The dashed e-likering could be due to a hard scatter. There are two Michel electrons(not shown). The visible energy is 742.8 MeV, the total
momentum is 445.8 MeV=c, and the invariant mass is1022.3 MeV=c2 .

072009-10

SEARCH FOR DINUCLEON DECAY INTO PIONS AT …

PHYSICAL REVIEW D 91, 072009 (2015)

TABLE VII. Candidate variable values, and mean variable values for signal and background distributions, for each of the variables
used in the BDT for the pp → π þ π þ search. Values are given for SK-II and SK-IV, for which single candidate events were found.
Candidate

SignalBackground

Variable

SK-II

SK-IV

SK-II

SK-IV

SK-II

SK-IV

Angle between μ-like rings (degrees)
Ratio of charge carried by most energetic ring (dimensionless)
Visible energy (MeV)Max. distance to Michel vertex (cm)Max. angle between μ-like ring and Michel vertex (degrees)
Magnitude of vector sum of corrected charge (corrected p.e.)
Number of rings (dimensionless)Number of Michel electrons (dimensionless)Number of μ-like rings (dimensionless)

177.7
0.60
400.0
1245.6
166.3
306.9
3
1
2

151.8
0.46
745.6
1420.1
111.9
1799.5
3
2
2

152.3
0.55
680.5
440.7
94.2
601.3
2.89
0.91
2.14

152.1
0.55
693.7
424.8
109.0
1310.0
2.95
1.09
2.16

138.3
0.66
440.2
523.0
93.2
448.7
2.39
1.00
2.05

138.1
0.66
407.3
459.2
110.9
797.5
2.40
1.21
2.05

and some of the multivariate input variables only apply tothe reconstructed π 0 .
First, we apply a set of selection criteria:(B1) There is more than one Cherenkov ring.
(B2) There is at least one e-like and one μ-like ring. Theπ þ is assumed to correspond to the most energetic
μ-like ring, while the π 0 is assumed to correspondto the e-like ring (if one e-like ring is found) or thetwo most energetic e-like rings (if more than one
is found).
(B3) The number of Michel electrons is no more than one.
(B4) The “reduced” visible energy, defined as the total
visible energy minus the energy of the π 0 as definedin the previous criterion, is less than 800 MeV. Thisis similar to criterion (A4) used in the pp → π þ π þsearch; the maximum visible energy from the π þ
is Mproton − Mπ þ ≈ 800 MeV.
(B5) The angle between the reconstructed π þ and π 0 is
greater than 120 degrees.
The signal efficiency, expected background, and data
events remaining after the selection criteria are shownin Table VIII. About a third of the primary π 0 ’s in the
pn → π þ π 0 events surviving the selection criteria arereconstructed as a single e-like ring. This is mainly because
the typically high momentum of the π 0 ’s leads to a small
opening angle between the photons in the subsequent
π 0 → γγ decay, so that the rings are not distinct.

Events surviving the selection criteria are passed into a
BDT. The following set of input variables is used, orderedaccording to variable importance.
(b1) The momentum of the π 0 candidate. The π 0 momentum peaks near 1 GeV for the signal, which istypically well reconstructed. The background peaks
at a much lower value of ∼200 MeV.
(b2) The angle between the reconstructed π þ and π 0 . Thispeaks near 180 degrees for the signal.
(b3) The momentum of the π þ candidate. The momentum
resolution for charged pions from dinucleon decay istypically poor. Despite this, it still peaks at a higher
value in the signal than the background (470 MeV=c
vs 390 MeV=c), and is relatively independent of thebetter reconstructed π 0 momentum.
(b4) The invariant mass of the π 0 . This can be estimatedfor both one- and two-ring π 0 hypotheses, by meansof the specialized algorithm referred to in Sec. IV.
(b5) The ratio of charge carried by the highest-energy ringto the total charge of all rings.
(b6) The total visible energy.
(b7) The number of Michel electrons.
Variables (b5), (b6), and (b7) were described for the
pp → π þ π þ search. The ranking of these variables isshown in Table IX. The two most discriminating input
variables are shown in Fig. 9. The π 0 -π þ angle is

TABLE IX. Relative importance of each variable in the
pn → π þ π 0 search, averaged across SK periods.
TABLE VIII. Efficiencies, expected backgrounds, and data
events after selection criteria in the pn → π þ π 0 search. Background rates are scaled to the appropriate SK period live time.
Statistical errors are shown.

Eff. (%)
Bkg.
data

SK-I

SK-II

SK-III

SK-IV

21.0  0.3
132  1.8
136

21.9  0.3
69  1.0
66

21.6  0.3
48  0.6
45

21.1  0.3
147  2.0
171

Variableπ 0 candidate momentum
Angle between π 0 and π þ candidates
π þ candidate momentum
π 0 candidate invariant mass
Ratio of charge carried by most energetic ring
Visible energy
Number of Michel electrons072009-11

Importance
0.19
0.17
0.16
0.15
0.14
0.14
0.058

J. GUSTAFSON et al.

PHYSICAL REVIEW D 91, 072009 (2015)

8
7
6
5
4
3
2
1
0
120

SK-I
25

Events

20
15
10
5
130

140

150

160

170

0
0

180

200

400

SK-II

SK-II

130

140

150

160

170

180

18
16
14
12
10
8
6
4
2
0
0

200

400

600

800 1000 1200 1400

π0 momentum (MeV/c)

SK-III

SK-III
12

7
6

10

5

8

EventsEvents800 1000 1200 1400

π0 momentum (MeV/c)

π0-π+ angle (degrees)

4
3

6
4

2

2

1
0
120

600

π0-π+ angle (degrees)

EventsEventsEvents

SK-I
18
16
14
12
10
8
6
4
2
0
120

130

140

150

160

170

0
0

180

200

π -π angle (degrees)
0

400

600

800 1000 1200 1400

π momentum (MeV/c)

+

0

SK-IV

SK-IV

18
16
14
12
10
8
6
4
2
0
120

30
25

EventsEvents

35

20
15
10
5

130

140

150

160

170

180

π0-π+ angle (degrees)

0
0

200

400

600

800 1000 1200 1400

π0 momentum (MeV/c)

FIG. 9 (color online). The π 0 -π þ angle and π 0 momentum for the pn → π þ π 0 search. Signal and atmospheric neutrino MC (solid anddashed histograms) are normalized to the number of events in the data (crosses) for each SK period. Details of each variable aredescribed in the text.

072009-12

SEARCH FOR DINUCLEON DECAY INTO PIONS AT …

PHYSICAL REVIEW D 91, 072009 (2015)

SK-II

45
40
35
30
25
20
15
10
5
-0.6

EventsEvents

SK-I

-0.4

-0.2

0

0.2

0.4

0.6

-0.6

-0.2

0

BDT Output

SK-III

SK-IV

0.2

0.4

0.6

0.2

0.4

0.6

60
50

12

EventsEvents-0.4

BDT Output

16
14
10
8
6

40
30
20

4
2
-0.6

22
20
18
16
14
12
10
8
6
4
2

10
-0.4

-0.2

0

0.2

0.4

0.6

-0.6

-0.4

-0.2

BDT Output

0

BDT Output

FIG. 10 (color online). BDT output for pn → π þ π 0 signal MC (solid), atmospheric neutrino background MC (dashed), and data(crosses). The green vertical line indicates the BDT cut value, and the arrow indicates that only events to the right of the cut are kept. One
candidate event can be seen in the SK-I data distribution, just to the right of the cut at 0.19.

highly discriminating for the same reason as the angle
between π þ candidates in the pp → π þ π þ search. The
π 0 momentum is also a highly discriminating variable,
since the peak near 1 GeV=c for the signal π 0 momentum (which is expected for back-to-back pions emittedwith initial energy ∼2 GeV) is absent in the background.
The data and atmospheric neutrino MC are in good
agreement.
The final BDT outputs are shown in Fig. 10. For the
pn → π þ π 0 search, an adaptive boost method [31] waschosen.

TABLE X. Final BDT cut value, efficiency, expected background, and candidate data events for the SK-I–SK-IV pn →π þ π 0 search, with statistical uncertainties shown for efficiency
and expected background. Background rates are quoted for theappropriate SK live time and per megaton-year.
SK-I

SK-II

SK-III

SK-IV

Cut
0.19
0.24
0.20
0.17
Eff. (%)
10.2  0.2 10.0  0.2 9.4  0.2 10.4  0.2
Bkg. (MT-yr)
2.7  0.7 2.3  0.7 2.2  0.7 2.9  0.8
Bkg. (SK live.)
0.25
0.11
0.07
0.32
Candidates
10
0
0

The final cut on the BDT output was designed so that lessthan one background event was expected across all SK
periods, with contributions from each SK period ordered bytheir live time. It is determined with MC only. The total
expected background is 0.75 events, consistent with the one
candidate event found in SK-I. Unlike the pp → π þ π þsearch, the signal distributions are not sharply peaked, sothat the efficiency does not change drastically as the cut
value is shifted.
The cut on the final BDT output, efficiencies, background rates, and candidates are shown in Table X. CCDIS
and CC1π are the main components of the remaining
background, detailed in Table XI.

TABLE XI. Neutrino interaction mode percentages for remaining background in the pn → π þ π 0 search for SK-I–SK-IV. The
percentages agree across SK periods within statistical fluctuations, which are large due to the low number of remaining events.
Mode

SK-I

SK-II

SK-III

SK-IV

CC1π
CCDIS
NCDIS
CCQE

33  17%38  17%8  8%
8  8%

29  17%46  19%20  14%5  5%

43  19%31  18%
21  15%5  5%

45  18%40  16%8  8%
8  8%

072009-13

J. GUSTAFSON et al.

PHYSICAL REVIEW D 91, 072009 (2015)

1200
960
720
480
240
0
0

500

1000

1500

2000

Times (ns)FIG. 11 (color online). A pn → π þ π 0 MC event. The μ-likering corresponds to a true π þ . The e-like ring contains twooverlapping γ’s from a π 0 . The angle between the μ- and e-likerings is 165 degrees. The e-like ring has a reconstructed
momentum of 893 MeV=c, while the original π 0 had momentum
876 MeV=c.

FIG. 13 (color online). The pn → π þ π 0 candidate in SK-I. The
e-like and μ-like rings are at an angle of 140 degrees. The e-likering has a high momentum of 987 MeV=c, but its fit π 0 mass islow, at 10 MeV=c2 . The μ-like ring has a momentum of
460.0 MeV=c. The total momentum and invariant mass are
700.1 MeV=c and 1266.4 MeV=c2 , respectively.

An example pn → π þ π 0 MC event display is shown in
Fig. 11, and an example background MC event is shownin Fig. 12. One candidate is found in SK-I, which is shownin Fig. 13. Table XII shows the values of the BDT input
variables for the candidate event, compared with the meansof the signal and background distributions. While thecandidate event appears closer to the signal mean values
for energy-related variables, it appears more backgroundlike for the π 0 − π þ angle and the reconstructed π 0 mass.
There are no candidate events in SK-II–SK-IV data.
Although the expected background in SK-I is much lower
than one event, it is consistent within about 1σ withthe expected background of 0.75 events in the entire
SK-I–SK-IV data set. The probability of observing one
or more background events in SK-I, without incorporating
systematic uncertainties, is 22%, and therefore the event is
consistent with background expectations. Further, the low
π 0 mass of the candidate event, as reconstructed by the
specialized algorithm discussed in Sec. IV, suggests that it
is not a pn → π þ π 0 signal. The e-like ring is more likely tobe an electron from, for instance, a CCνe 1π event. Thus,

TABLE XII. Candidate variable values, and mean variablevalues for signal and background distributions, for each of thevariables used in the BDT for the pn → π þ π 0 search. Letters foreach variable are assigned in the text. Values are given for SK-I,
for which a single candidate event was found.
Variable

FIG. 12 (color online). A surviving atmospheric neutrino MC
event for the pn → π þ π 0 search. The interaction type is CC1π 0 ,νμ n → μ− π 0 p. The two overlapping e-like rings are photons froma π 0 , with total momentum 939 MeV=c. The μ-like ringcorresponds to the μ− , with momentum 410 MeV=c. The angle
between them is 151 degrees, thereby appearing back to back.

Candidate Signal Background

π 0 momentum (MeV=c)986.5 756.0
π 0 -π þ angle (degrees)
142.8 152.9
π þ momentum (MeV=c)459.5 470.3
π 0 mass (MeV=c2 )10.0 143.0
Charge ratio (dimensionless)
0.77
0.56
Visible energy (MeV)
1171.6 1098.6
Michel electrons (dimensionless)
00.40

072009-14

233.7
140.4
394.5
102.2
0.58
495.0
0.51

SEARCH FOR DINUCLEON DECAY INTO PIONS AT …

PHYSICAL REVIEW D 91, 072009 (2015)þ 0

we conclude that there is no evidence for pn → π π in the
SK-I–SK-IV data set.
3. nn → π0 π0The nn → π 0 π 0 search is distinct from the other two
searches, in that no multivariate method is used. Similar
to other nucleon decay searches at SK [4–6], the total
momentum and invariant mass are highly discriminatory
variables, so that the search proceeds with a simple set of
selection criteria. The criteria are similar to those used forthe p → eþ π 0 analysis [4].
(C1) 2 ≤ number of Cherenkov rings ≤ 4.
(C2) All rings are e-like.
(C3) No Michel electrons.
(C4) Ptot ≤ 600 MeV=c.
(C5) 1600 MeV=c2 ≤ Minv ≤ 2000 MeV=c2 .
Pnring~ i j,The total momentum Ptot is defined as Ptot ¼ j i¼1
p~ i is thewhere nring is the number of reconstructed rings, p
3-momentum of thepith
ring,
and
the
invariantmass isﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
Pnring2
2
defined as M inv ¼ Etot − Ptot , where Etot ¼ i¼1 j~
pi j(the rings are assumed to come from photons, so there is norest-mass energy). A nn → π 0 π 0 MC event is shown in
Fig. 14, and a background MC event surviving all selection
criteria is shown in Fig. 15.
Criterion (C1) allows for both one- and two-ring π 0 fits.
Criteria (C2) and (C3) are designed to eliminate muonsand charged pions. Criterion (C4) is chosen to ensure a low
total momentum, and to reject cases where only one π 0 is

FIG. 14 (color online). A nn → π 0 π 0 MC event. The twooverlapping e-like rings on the top left correspond to two γ’s froma π 0 decay, and the other ring also corresponds to a π 0 . The true
momentum of the π 0 corresponding to the two-ring fit is872.7 MeV=c, and the true momentum of the π 0 corresponding
to the one-ring fit is 936.7 MeV=c. The total reconstructed
momentum of this event is 219 MeV=c, and the reconstructedinvariant mass is 1793 MeV=c2 .

FIG. 15 (color online). Surviving atmospheric neutrino MC
event for the nn → π 0 π 0 search. The interaction is νμ neutralcurrent deep inelastic scattering. In particular, the interaction
produced a π 0 and a η meson, from which four γ’s are produced.
All but the left-most ring correspond to true photons from the π 0and η meson. Many additional charged pions are also produced,
concentrated near the center of the display, thereby complicatingreconstruction. The left-most ring corresponds to a charged pion.
There is a ring near the center from multiple, closely overlapping
charged pion tracks, which is not found by the reconstruction
program. The total reconstructed momentum is 553 MeV=c, andthe reconstructed invariant mass is 1763 MeV=c2 .

reconstructed, as each π 0 from nn → π 0 π 0 has a momentum that peaks closer to 1 GeV=c. Criterion (C5) is chosen
to cover the ∼2Mn region as tightly as possible without
getting too close to the peak (Fig. 16). This eliminatesnearly all background, which peaks sharply near M π0 .
The signal efficiency, background, and data events after
each selection criteria are shown in Table XIII. The final
efficiencies and background rates are summarized inTable XIV. The remaining background is mostly neutral
current deep inelastic scattering (Table XV). The interaction mode percentages agree across SK periods within
statistical fluctuations, which are large due to the very low
number of remaining events. The lower background in
SK-IV is due to improved Michel electron tagging relative
to SK-I, -II, and -III.
Figure 17 shows the final total momentum and invariant
mass distributions for signal and background MC, and thedata. No candidate events are found for this mode, which isconsistent with the near-zero expected background for all
SK periods.
B. Systematic uncertainties
Systematic uncertainties for the dinucleon decay
searches are comprised of MC simulation uncertainties
specific to the signal and background MC, and common
signal and background uncertainties related to event
reconstruction and BDT bias (for modes which use a BDT).

072009-15

J. GUSTAFSON et al.

PHYSICAL REVIEW D 91, 072009 (2015)

SK-I

12

500

10

400

EventsEvents

SK-I

8
6
4

200
100

2
0

1000 2000 3000 4000 5000
Total Momentum (MeV/c)

6000

0

500

SK-II

250

5

200

Events

6

4
3

2000

2500

150
100
50

1
0

1000 2000 3000 4000 5000
Total Momentum (MeV/c)

6000

0

500

SK-III

Events

2.5
2
1.5
1
0.5
1000 2000 3000 4000 5000
Total Momentum (MeV/c)

6000

0

2500

500
1000
1500
2000
Invariant Mass (MeV/c 2)

2500

SK-IV
600

10

500

Events

12

8
6

400
300

4

200

2

100
1000 2000 3000 4000 5000
Total Momentum (MeV/c)

2000

200
180
160
140
120
100
80
60
40
20

SK-IV

0

1500

SK-III

301000

Invariant Mass (MeV/c 2 )3.5

Events

1500

SK-II

300

2

Events1000

Invariant Mass (MeV/c 2)

7

Events

300

6000

0

500

1000

1500

2000

2500

Invariant Mass (MeV/c 2 )FIG. 16 (color online). Total momentum (left) and invariant mass (right) for data (crosses), nn → π 0 π 0 MC (solid histogram), andatmospheric neutrino MC (dashed histogram). Criteria (C1)–(C3) have been applied to all distributions. The invariant mass criterion(C5) has been applied to the total momentum, and the total momentum criterion (C4) to the invariant mass. The thin green vertical lines
indicate the locations of the final cuts for each variable.

072009-16

SEARCH FOR DINUCLEON DECAY INTO PIONS AT …
TABLE XIII. Efficiency, expected background, and data events
for the SK-I–SK-IV nn → π 0 π 0 search, for each analysis selection criterion. Statistical errors are shown. The criteria aredescribed in the text.
SK-I

SK-II

SK-III

SK-IV

FCFV
Eff. (%) 94.1  0.7 96.6  0.7 94.4  0.7 94.1  0.7
Bkg.
11846  11.1 6389  5.8 4187  3.7 14214  12.7
Data
12299
6610
4355
14444
(C1)
Eff. (%)
Bkg.
Data

76.7  0.6
3434  5.6
3558

76.4  0.6 76.6  0.6
1846  3.0 1204  2.0
1969
1239

76.8  0.6
4135  6.7
4215

(C2)
Eff. (%)
Bkg.
Data

61.0  0.5
2131  4.3
2249

59.1  0.5 59.6  0.5
1175  2.3 728  1.5
1268
782

60.1  0.5
2502  5.0
2623

(C3)
Eff. (%)
Bkg.
Data

54.4  0.5
1441  3.4
1492

52.6  0.5
799  1.9
837

52.8  0.5
497  1.2
529

52.1  0.5
1543  1.1
1627

(C4)
Eff. (%)
Bkg.
Data

30.1  0.4
837  2.6853

27.0  0.4
470  1.4
503

29.1  0.4
290  0.9
310

29.3  0.4
940  3.0
1021

(C5)
Eff. (%) 22.1  0.3 18.8  0.3 20.9  0.3 21.4  0.3
Bkg.
0.05  0.02 0.04  0.01 0.03  0.01 0.02  0.01
Data
0
0
0
0

The largest systematic uncertainties for the dinucleon
decay signal come from the simulation stage. The dominantsystematic uncertainty for all modes comes from pion finalstate interaction (FSI) uncertainties. To estimate this, we
generated separate sets of dinucleon decay MC for 24
parameter variations consisting of 1-sigma changes in pionnucleon cross sections, corresponding to uncertainties
determined from fits to pion scattering data. This is thesame set of variations used in the T2K νe appearance result[30]. The variations affect the cross section for both lowenergy (charge exchange, absorption, quasielastic scattering) and high-energy (hadron production) interactions.
These variations were propagated through all stages of
TABLE XIV. Efficiency, expected background, and candidate
data events for the SK-I–SK-IV nn → π 0 π 0 search, with statistical uncertainties. Background is quoted both for each SK live
time and in megaton-years.
SK-I

SK-II

SK-III

SK-IV

Efficiency (%) 22.1  0.3 18.8  0.3 20.9  0.3 21.4  0.3
bkg. (MT-yr)
0.6  0.2 0.9  0.3 1.0  0.3 0.2  0.1
bkg. (SK live.)
0.05
0.04
0.03
0.02
candidates
0000

PHYSICAL REVIEW D 91, 072009 (2015)
TABLE XV. Neutrino interaction mode percentages for theremaining background in the nn → π 0 π 0 search for SK-I–SK-IV.
Mode

SK-I

SK-II

SK-III

SK-IV

NCDIS
CCDIS
CC1π

63  32%
15  16%21  15%30  18%50  22%20  14%67  25%
24  14%9  9%

50  50%0 þ 50%51  51%

the analysis. For the pp → π þ π þ and pn → π þ π 0 modes,
each set of MC was processed with the trained, unvaried
BDT used for the data and analysis MC. The systematic
error is conservatively taken to be the largest percentagechange in efficiency across all variations upon applying allof the analysis selection criteria. Other significant simulation systematic uncertainties come from Fermi motion(∼20%) and correlated decay (∼10%). These are not passed
through the BDT, but rather are estimated through a simple
reweighting method of the final event sample. For Fermi
motion, the reweighting comes from differences between
the Fermi momentum simulation for dinucleon decay andthe Fermi gas model. For correlated decay, the effect is notwell understood, and so a conservative 100% uncertainty
is placed on the detection efficiency of events with acorrelated decay, which is 10% at the simulation stage andstill about 10% in the final event sample.
The largest overall contributions to background systematic uncertainties are related to neutrino interactions. In
particular, for pp → π þ π þ and pn → π þ π 0 , the uncertainties associated with charged-current single pion production
are dominant. This is mainly due to a 40% uncertainty inthe π 0 =π þ production ratio. Other simulation systematic
uncertainties due to pion FSI and neutrino flux are also
considered.
The most important reconstruction systematic uncertainty for both dinucleon decay and atmospheric neutrinos
comes from the energy scale. The energy-scale uncertainty
is determined from comparing the absolute momenta ofdata and MC for several control samples, such as cosmicray muons, their resulting decay electrons, and π 0 ’s fromneutral-current interactions. Time variation in the reconstructed momenta of the control samples is also taken into
account in the uncertainty calculation. The energy scale
uncertainty is estimated to be 1.1%, 1.7%, 2.7%, and2.3% for SK-I, -II, -III, and -IV, respectively [5,30]. This
uncertainty is propagated into all energy-dependent variables, and is propagated through the trained, unvaried BDT
used for the data and analysis MC for pp → π þ π þ andpn → π þ π 0 . The energy-scale uncertainty is the dominantsystematic uncertainty for the atmospheric neutrino background in the nn → π 0 π 0 search. Other reconstructionrelated uncertainties taken into account are ring separation,
particle identification, and the ring directional resolution.
Finally, for the multivariate method, an additional
systematic uncertainty for the training bias of the BDT

072009-17

J. GUSTAFSON et al.

PHYSICAL REVIEW D 91, 072009 (2015)

Total Momentum (MeV/c)

nn → π0π0 SK-I-IV MC

atm. ν SK-I-IV MC

SK-I-IV Data

1800
1600
1400
1200
1000
800
600
400
200
0

500

1000

1500

2000

500

1000

1500

2000

500

1000

1500

2000

Total Invariant Mass (MeV/c 2)

FIG. 17 (color online). The total momentum and invariant mass distributions for nn → π 0 π 0 MC (left), atmospheric neutrino MC(center), and data (right) for SK-I–SK-IV. As seen in Fig. 16, there are no major structural differences between the SK-I–SK-IV
distributions, so events from all periods are combined. Events have passed criteria (C1)–(C3). The solid box indicates the total
momentum and invariant mass cuts [(C4) and (C5) in the text].

is considered. This is estimated by assessing the differencein the testing and analysis distributions of the BDT output
for both signal and background. Tables XVI–XVIII break
down the systematic uncertainties for each mode and each
SK period.

VI. RESULTS AND DISCUSSION
Since no significant excess was observed in the data, weset a lower limit on the lifetime of each dinucleon decay
mode. This is done using a Bayesian method, which
incorporates systematic uncertainties as follows [32]:ZZZ

TABLE XVI. Systematic uncertainties for the pp → π þ π þmode.

PðΓjnÞ ¼ A

pp → π þ π þð2Þ

× PðΓÞPðλÞPðϵÞPðbÞdλdϵdb;

Signal (%)

SK-I

SK-II

SK-III

SK-IV

Simulation
Reconstruction
BDT
Total

35.2
6.0
3.6
35.9

35.1
8.6
2.2
36.2

33.6
4.0
4.4
34.1

38.5
3.2
2.0
38.7

Background (%)

SK-I

SK-II

SK-III

SK-IV

Simulation
Reconstruction
BDT
Total

29.1
6.1
6.8
30.5

29.1
8.1
1.0
30.3

35.8
4.1
4.3
36.4

26.5
3.2
1.4
26.8

TABLE XVII.
mode.

e−ðΓλϵþbÞ ðΓλϵ þ bÞnn!

Systematic uncertainties for the pn → π þ π 0
pn → π þ π 0

Signal (%)

SK-I

SK-II

SK-III

SK-IV

Simulation
Reconstruction
BDT
Total

33.3
3.3
<1
33.4

32.2
1.7
1.6
32.3

28.4
5.6
<1
28.9

35.0
5.6
<1
35.3

Background (%)

SK-I

SK-II

SK-III

SK-IV

Simulation
Reconstruction
BDT
Total

22.1
1.8
6.3
23.1

19.9
1.8
7.4
21.3

24.0
3.3
10.3
26.3

27.8
3.8
11.3
28.6

where Γ is the true dinucleon decay rate, n is the number ofcandidate events, λ is the exposure, ϵ is the efficiency, b is
TABLE XVIII.
mode.

Systematic uncertainties for the nn → π 0 π 0nn → π 0 π 0Signal (%)

SK-I

SK-II

SK-III

SK-IV

Simulation
Reconstruction
Total

31.1
1.5
31.2

34.4
1.7
34.4

37.3
4.0
37.6

33.1
3.6
33.3

Background (%)

SK-I

SK-II

SK-III

SK-IV

Simulation
Reconstruction
Total

13.6
10.9
17.5

15.5
18.1
24.0

14.5
28.9
32.3

13.9
24.3
28.0

TABLE XIX. Limits on each mode, compared with previous
results in Ref. [12].
Modepp → π þ π þ
pn → π þ π 0
nn → π 0 π 0072009-18

Frejus limit (56 Fe)

This analysis (16 O)

7.0 × 1029 yrs
2.0 × 1030 yrs
3.4 × 1030 yrs

7.22 × 1031 yrs
1.70 × 1032 yrs
4.04 × 1032 yrs

SEARCH FOR DINUCLEON DECAY INTO PIONS AT …

the backgroundrate, and A is a normalization constant thatR
ensures 0∞ PðΓjnÞdΓ ¼ 1.
The distributions PðλÞ, PðϵÞ, and PðbÞ are taken to be
Gaussian, while PðΓÞ is uniform for Γ > 0. The uncertainty on λ is assumed to be negligible. The exposure is
16
O-years, not nucleon-years as in most proton decay
analyses. This is consistent with the exposure used inRefs. [11] and [12], which used total 16 O and 56 Fe nuclei,
respectively.
Finally, the confidence level is calculated as
C:L: ¼

Z Γ
0

limit

PðΓjnÞdΓand the lower lifetime limit is just the reciprocal of Γlimit .
Table XIX shows the lifetime limits at the 90% confidence
level, compared with the previous limits set by the Frejus
Collaboration in Ref. [12].
The interpretation of our results depends on the underlying physics model. For instance, our limit for nn → π 0 π 0is the most stringent experimental limit for dinucleon
decay to light mesons, including that from Ref. [11],which measured τpp→Kþ Kþ > 1.77 × 1032 years. However,Ref. [11] was motivated by supersymmetric models thatforbid couplings between quarks of the same generation.
If (as with the models considered in Ref. [33]) suchcouplings are not forbidden, and no preference is given
for couplings to first- or second-generation quarks, ourresult for nn → π 0 π 0 provides the strongest general constraint on dinucleon decay to mesons.
As stated in Sec. I, there is a connection between
n-nbar and dinucleon decay to pions, if the same operator
is dominant for both processes. This is discussed further
in the Appendix. Our results provide a useful cross-checkto n-nbar searches. By isolating a single final statefor each mode, as opposed to the many non-negligible
final states of n-nbar [34], these analyses are simple androbust.

PHYSICAL REVIEW D 91, 072009 (2015)

ACKNOWLEDGMENTS
We gratefully acknowledge the cooperation of the
Kamioka Mining and Smelting Company. The SuperKamiokande experiment has been built and operated fromfunding by the Japanese Ministry of Education, Culture,
Sports, Science and Technology, the United States
Department of Energy, and the U.S. National Science
Foundation. This work was partially supported by the
Research Foundation of Korea (BK21 and KNRC), theKorean Ministry of Science and Technology, the National
Science Foundation of China, the European Union FP7 (DSlaguna-lbno PN-284518 and ITN invisibles GA-2011289442), the National Science and Engineering ResearchCouncil (NSERC) of Canada, and the Scinet and West-gridconsortia of Compute Canada.
APPENDIX: DINUCLEON DECAY LIFETIMES
Here we discuss the lifetime formula for dinucleon decay
to pions. We assume that the same operator which mediates
n-nbar also mediates dinucleon decay to pions, which
follows from crossing symmetry.
Reference [33] discussed both n-nbar and dinucleon
decay with the same operator. This operator takes thegeneral form
OΔB¼2 ¼

1
O6Q ;
Λ5sup

where O6Q is a six-quark operator containing two up-type
and four down-type quarks, and Λsup is a suppression scale,with the fifth power following from dimensional analysis.
Indices related to quark generation, color, etc., are suppressed. For all first-generation quarks, this operator leads
to both n-nbar and dinucleon decay. The (free) n-nbar
lifetime is given byτnn̄ ≈

Λ5sup
;
Λ6QCD

ðA1Þ

where ΛQCD is a hadronic matrix element near the QCD
scale, estimated to be 200 MeV. The dinucleon decay
lifetime τNN is given by [33,35]

VII. CONCLUSION
We have performed searches for dinucleon decay topions via the processes pp → π þ π þ , pn → π þ π 0 , andnn → π 0 π 0 with 282.1 kiloton-years of data with theSuper-Kamiokande detector. This is the first search fordinucleon decay to pions done in a water-Cherenkov
detector. We do not observe any significant event excess
for any of the modes, and the data remain consistent
with the atmospheric neutrino background. Therefore,
partial lifetime limits of τpp→πþ πþ > 7.22 × 1031 years,
τpn→πþ π0 > 1.70 × 1032 years, and τnn→π0 π0 > 4.04 ×
1032 years per 16 O nucleus have been set. The limits per
nucleus are about 2 orders of magnitude more stringent
than the previous limits set in Ref. [12], and constrainmodels that predict ΔB ¼ 2 processes.

τNN ≈πM 2N Λ10
sup;
8ρN Λ10QCD

ðA2Þ

where M N is the nucleon mass and ρN ≈ 0.25 fm−3 is theaverage nuclear density. From Eqs. (A1) and (A2),τNN ≈πM 2N 2Λ
τ2 ¼ T nuc τ2nn̄ :8ρN QCD nn̄ðA3Þ

Equation (A3) is mathematically similar to the formalism
for comparing free and bound n-nbar lifetimes, although inthe latter case the factor T nuc comes from different physics
related to the difference in neutron and antineutron nuclear
potentials. For dinucleon decay, T nuc ≈ 1.1 × 1025 s−1 ,072009-19

J. GUSTAFSON et al.

PHYSICAL REVIEW D 91, 072009 (2015)

while for n-nbar T nuc is 2 orders of magnitude less,with Ref. [36] calculating it to be 0.517 × 1023 s−1 . The
current lower limit on the n-nbar lifetime, from SuperKamiokande, is τnn̄ > 2.7 × 108 s [34]. This translates

to a lower lifetime for dinucleon decay to pions ofτNN→ππ ≈ 2.5 × 1034 years, about 2 orders of magnitude
beyond the limits obtained in this paper. Thus, our results
do not place any new constraints on n-nbar.

[1] A. D. Sakharov, Pis’ma Zh. Eksp. Teor. Fiz. 5, 32 (1967)[JETP Lett. 5, 24 (1967)]; Usp. Fiz. Nauk 161, 61 (1991)[Sov. Phys. Usp. 34, 392 (1991)].
[2] K. Abe et al. (Super-Kamiokande Collaboration), Phys.
Rev. D 90, 072005 (2014).
[3] K. Abe et al. (Super-Kamiokande Collaboration), Phys.
Rev. Lett. 113, 121802 (2014).
[4] H. Nishino et al. (Super-Kamiokande Collaboration), Phys.
Rev. D 85, 112001 (2012).
[5] C. Regis et al. (Super-Kamiokande Collaboration), Phys.
Rev. D 86, 012006 (2012).
[6] K. Kobayashi et al. (Super-Kamiokande Collaboration),
Phys. Rev. D 72, 052007 (2005).
[7] J. L. Raaf (Super-Kamiokande Collaboration), Nucl. Phys.
B, Proc. Suppl. 229-232, 559 (2012).
[8] H. Georgi and S. L. Glashow, Phys. Rev. Lett. 32, 438(1974).
[9] J. M. Arnold, B. Fornal, and M. B. Wise, Phys. Rev. D 87,
075004 (2013).
[10] K. Babu et al., arXiv:1310.8593.
[11] M. Litos et al. (Super-Kamiokande Collaboration), Phys.
Rev. Lett. 112, 131803 (2014).
[12] C. Berger et al. (Frejus Collaboration), Phys. Lett. B 269,
227 (1991).
[13] Y. Fukuda et al. (Super-Kamiokande Collaboration), Nucl.
Instrum. Methods Phys. Res., Sect. A 501, 418 (2003).
[14] H. Nishino, K. Awai, Y. Hayato, S. Nakayama, K. Okumura,
M. Shiozawa, A. Takeda, K. Ishikawa, A. Minegishi, and
Y. Arai, Nucl. Instrum. Methods Phys. Res., Sect. A 610,
710 (2009).
[15] K. Abe et al., Nucl. Instrum. Methods Phys. Res., Sect. A
737, 253 (2014).
[16] Y. Ashie et al. (Super-Kamiokande Collaboration), Phys.
Rev. D 71, 112005 (2005).
[17] K. Nakamura, S. Hiramatsu, T. Kamae, H. Muramatsu,
N. Izutsu, and Y. Watase, Nucl. Phys. A268, 381 (1976).

[18] T. Yamazaki and Y. Akaishi, Phys. Lett. B 453, 1 (199\9).
[19] R. D. Woods and D. S. Saxon, Phys. Rev. 95, 577(1954).
[20] Y. Hayato, Nucl. Phys. B, Proc. Suppl. 112, 171 (2002).
[21] CERN program library W5013 (1994).
[22] M. Nakahata et al., J. Phys. Soc. Jpn. 55, 3786 (1986);
A. S. Clough et al., Nucl. Phys. B76, 15 (1974).
[23] J. P. Albanese, J. Arvieux, J. Bolger, E. Boschitz, C. H. Q.
Ingram, J. Jansen, and J. Zichy, Nucl. Phys. A350, 301(1980).
[24] M. Honda, T. Kajita, K. Kasahara, S. Midorikawa, and
T. Sanuki, Phys. Rev. D 75, 043006 (2007).
[25] D. Rein and L. M. Sehgal, Ann. Phys. (Berlin) 133, 79(1981).
[26] M. Davies, Machine Vision: Theory, Algorithms, Practicalities (Academic Press, San Diego, 1997).
[27] M. Shiozawa (Super-Kamiokande Collaboration), Nucl.
Instrum. Methods Phys. Res., Sect. A 433, 240 (1999).
[28] T. Barszczak, Report No. UMI-31-71221. Ph.D thesis,
2005.
[29] R. Wendell et al. (Super-Kamiokande Collaboration),
Phys. Rev. D 81, 092004 (2010).
[30] K. Abe et al. (T2K Collaboration), Phys. Rev. D 88, 032002(2013).
[31] A. Hoecker et al., Proc. Sci., ACAT (2007) 040, arXiv:
physics/0703039. See especially Section 7 and Section 8.12.
[32] R. M. Barnett et al. (Particle Data Group Collaboration),
Phys. Rev. D 54, 1 (1996).
[33] C. Csaki, E. Kuflik, and T. Volansky, Phys. Rev. Lett. 112,131801 (2014).
[34] K. Abe et al. (Super-Kamiokande Collaboration), arXiv:
1109.4227.
[35] J. L. Goity and M. Sher, Phys. Lett. B 346, 69 (1995); 385,
500(E) (1996).
[36] E. Friedman and A. Gal, Phys. Rev. D 78, 016002(2008).

072009-20