Search for proton decay into three charged leptons in 0.37 megaton-years exposure of the Super-Kamiokande

A search for proton decay into three charged leptons has been performed by using 0.37$\,$Mton$\cdot$years of data collected in Super-Kamiokande. All possible combinations of electrons, muons and their anti-particles consistent with charge conservation were considered as decay modes. No significant excess of events has been found over the background, and lower limits on the proton lifetime divided by the branching ratio have been obtained. The limits range between $9.2\times10^{33}$ to $3.4\times10^{34}$ years at 90$\,$% confidence level, improving by more than an order of magnitude upon limits from previous experiments. A first limit has been set for the $p\rightarrow\mu^-e^+e^+$ mode.

A search for proton decay into three charged leptons has been performed by using 0.37 Mton·years of data collected in Super-Kamiokande. All possible combinations of electrons, muons and their antiparticles consistent with charge conservation were considered as decay modes. No significant excess of events has been found over the background, and lower limits on the proton lifetime divided by the branching ratio have been obtained. The limits range between 9.2 × 10 33 to 3.4 × 10 34 years at 90% confidence level, improving by more than an order of magnitude upon limits from previous experiments. A first limit has been set for the p → µ − e + e + mode.

I. INTRODUCTION
The Standard Model of elementary particles describes strong, weak and electromagnetic interactions based on gauge symmetries. Grand Unified Theories (GUTs) [1] can unify three interactions in the Standard Model in a single gauge group with one coupling constant. In most GUTs, grand unification is predicted typically at 10 15−16 GeV which is unreachable by accelerators, whereas the effects of the grand unification might be detected through rare phenomena beyond the Standard Model. The most promising such phenomenon is violation of baryon number, and the most sought after signature is proton decay [2]. The Super-Kamiokande (SK) experiment has been leading the search for proton decay and has set the most stringent limits on the lifetime for various channels predicted by GUT models. For example, the p → e + π 0 and p → νK + are favored decay modes in non-supersymmetric and TeV-scale supersymmetric GUT models, respectively, yet no significant signal was observed, excluding proton lifetimes up to 10 34 years [3,4]. The simplest unification model, minimal SU(5) [5], has been ruled out by SK and earlier experiments. Other nucleon decay channels motivated by unification, such as the charged antilepton plus meson channels were searched for in SK recently; no proton decay signal was found [6].
However, baryon number may be violated irrespective of GUTs and is expected in many scenarios beyond the Standard Model [7]. If the usual lower dimensional operators responsible for p → e + π 0 etc. are suppressed, then different proton decay channels can dominate from higher dimensional operators. This can naturally occur when considering lepton flavor symmetries [8]. In this case, trilepton nucleon decay channels such as p → µ − e + e + or p → e − µ + µ + can be dominant. As they are generated from effective dimension d = 10 operators, these processes probe scales of around 100 TeV. A minimal model for this based on leptoquarks has been put forward in Ref. [8], which also suggested that such processes could be connected to the recent anomalies observed in B-meson decays.
Some of these decay modes were already searched for by the Irvine-Michigan-Brookhaven-3 (IMB-3) [9] and Harvard-Purdue-Wisconsin (HPW) [10] experiments. The data were consistent with the expected background in both experiments and no significant signal was confirmed. Lifetime limits were set to be 10 30 -10 32 years for each decay mode. SK can significantly extend the search of the previous experiments.

II. THE SUPER-KAMIOKANDE EXPERIMENT
SK is the largest pure water Cherenkov detector, located at the Kamioka mine in Gifu prefecture, Japan. The SK detector consists of a stainless steel tank (39.3 m diameter, 41.4 m height), 50 kton of ultra pure water and photomultiplier tubes (PMTs). In order to reduce cosmic ray muon background, the detector is located 1,000 m under the peak of Mt. Ikenoyama (2,700 m water equivalent). The water tank is optically separated into two concentric cylindrical volumes by the support frames equipped with inward-facing 20-inch and outward-facing 8-inch PMTs. The inner volume is 33.8 m in diameter and 36.2 m in height, and is called the Inner Detector (ID). The ID contains 32 kton of water and monitored by 11,129 inward facing PMTs (about half for certain period as explained later). Outside of the ID is a 2 m thick water called the Outer Detector (OD). The OD is monitored by outward-facing 1,885 PMTs and mainly serves as an active cosmic ray muon veto as well as a shield against gamma rays from surrounding rock. 20-inch and 8-inch PMTs are uniformly mounted on the ID and OD surfaces, respectively. The details of SK detector are described elsewhere [11,12].
The analysis in this paper uses data taken with four different detector configurations. SK-I started in 1996 and stopped in 2001 for maintenance. SK-II was operated from 2002 to 2005 with about half the number of ID PMTs compared to SK-I due to the accident during the maintenance work after SK-I. In order to avoid further such accidents, the PMTs were protected by covers made of fiber reinforced plastic and acrylic for the photocathode starting from SK-II. SK-III started in 2006 and stopped in 2008. The number of PMTs in SK-III was recovered to almost the same number as SK-I. Readout electronics and data acquisition system were upgraded for the SK-IV period. SK-IV continued from 2008 and ended for the upgrade of SK in May 2018. Photocoverage of the ID is 19% in SK-II and 40% in other periods.
In this analysis, all the detected particles must be fully contained (FC) in the ID with the reconstructed vertex inside the fiducial volume (FV). Such events are selected by preselection cuts [3,4,6]. The FV is defined as the volume 2 m inside the top, bottom and barrel walls of the ID and corresponds to 22.5 kton mass. Contamination of non-neutrino background events due to cosmic muons is negligible after the preselection cuts. Data around 1 ms of the expected neutrino beam arrival timing from the Japan Proton Accelerator Research Complex (J-PARC), which has a repetition rate of 2.48 s, have been removed. All the available data from SK-I to SK-IV are used in this analysis. We use the data of 372.6 kton·years in total by summing up 91.7, 49.2, 31.9 and 199.8 kton·years of SK-I, II, III and IV data, respectively.

III. SIMULATION
H 2 O molecules are the sources of proton decay in SK searches, with 2 protons from hydrogen and 8 protons from oxygen. In the simulation, only a uniform phase space is assumed for kinematics of outgoing charged leptons and any additional correlations are not taken into account. The protons in hydrogen (free protons) have an initial mass and momentum of 938.27 MeV/c 2 and 0 MeV/c, respectively. On the other hand, the protons in oxygen (bound protons) interact with other nucleons and have some initial momentum. In the simulation, three nuclear effects are taken into account: nuclear binding energies in 16 O, Fermi motion, and correlated decay. Two nuclear binding energies (p-state and s-state) are accounted for with Gaussian spreads and are subtracted from the initial proton mass [6]. Fermi motion is estimated based on electron-12 C scattering data [13]. The bound proton sometimes correlates with the surrounding nucleons during its decay. This effect is predicted to occur with 10% probability and produces a broad distribution at lower mass in proton mass distribution [14].
For the background, only atmospheric neutrino events are considered since other non-neutrino backgrounds are negligibly small. The simulation of this sample consists of three steps: neutrino production in the atmosphere (neutrino flux prediction), neutrino interaction with water and particle tracking in the detector. The flux of atmospheric neutrino is calculated by the model of M. Honda et al. [15,16]. The interactions of atmospheric neutrinos with hydrogen or oxygen nuclei in water are simulated by the NEUT program [17]. This simulation covers a wide neutrino energy range from several tens of MeV to hundreds of TeV. Hadrons generated by neutrino interactions in the oxygen nucleus often cause secondary interactions within the nucleus. The interactions of pions, kaons, etas and nucleons in the target nucleus are simulated in NEUT by using a cascade model [17,18]. Simulated data samples equivalent to 500 years of detector exposure are generated for each SK period.
The particle propagation, Cherenkov radiation, propagation of Cherenkov photons in water and PMTs, and electronics response are simulated by a Geant3 based package [19] with custom modifications for use in SK, such as pion interactions in water and wavelengthdependent water transparency.
The simulation scheme for the signal sample is the same as for the other recent SK nucleon decay searches [3,4,6], and the latest SK oscillation analysis [20] for background sample. The oscillation parameters are taken from the latest atmospheric neutrino oscillation analysis [20].

IV. EVENT RECONSTRUCTION
Events with charged particles are reconstructed by using charge and time information of the hit PMTs. A reconstruction scheme consists of vertex fitting, ring counting, particle identification (PID), momentum reconstruction, Michel electron search and neutron tagging. The reconstruction method is almost the same as the one used in other recent nucleon decay searches [3,4,6] or oscillation studies [20] in SK. Some updates for charge separation and neutron tagging improve the sensitivity of this search.
In the first step, the vertex is reconstructed by assuming that Cherenkov light comes from one point at the same time. Then the ring edge and the direction of the ring are estimated. Finally, the vertex is reconstructed more precisely by assuming that photons are emitted along the track of the charged particle.
Cherenkov rings are then reconstructed by using the pattern recognition algorithm known as the Hough transformation [21]. Ring candidates are evaluated by a likelihood method to determine if they are true or fake. In case more than one rings (multi-ring) are identified, the contribution of each ring to the detected photoelectrons in each PMT is estimated. The opening angle of the ring can be calculated from the reconstructed vertex position and the edge of the ring. The final stage of ring counting discards the candidate rings mostly caused by multiple Coulomb scattering of charged particles by a ring's angle relative to other rings and by visible energy.
Reconstructed rings can be classified as electron-like (e-like) or muon-like (µ-like) by using the pattern of PMT hits. Cherenkov rings of muons tend to have clear ring edges. In contrast, Cherenkov rings of electrons tend to be relatively diffuse due to electromagnetic showers and scattering. Expected PMT charge patterns for electrons and muons are compared with observed hit patterns using a likelihood function. Information about the opening angle is included in the likelihood function. With the emission of Cherenkov rings, gamma rays are usually identified as e-like and charged pions as µ-like.
The momentum of the charged particle is reconstructed from the total number of photoelectrons in a 70 degree cone around the ring, which is corrected by using a conversion table depending on the particle type. We correct for time drift of the gain, which varies according to the year of PMT manufacture. For the multi-ring case, the expected charge distribution for each ring is calculated. Then the momentum is assigned to each ring according to the fraction of expected charge. For this charge separation algorithm, the expected charge in the backward direction of a µ-like ring was tuned to reproduce the data. The precision of total mass reconstruction was improved with this tuning, especially for the p → µ + µ + µ − events. The energy scale of the detector is checked precisely by using Michel electrons, stopping muons and neutral pion samples [20].
Michel electrons are tagged by searching for PMT hit clusters after the primary event. Since about 20% of µ − are captured by nuclei and do not emit a decay electron, the tagging efficiency for µ − is lower than that for µ + . Free neutrons traveling in water are thermalized and captured by oxygen or hydrogen nuclei. Neutrons are predominantly captured by the interaction n + p → d + γ (2.2 MeV). This 2.2 MeV gamma ray signal is searched for to identify the neutron (neutron tagging). The capture signal occurs a few hundred microseconds after the initial neutrino interaction signal, and the tagging was possible only with the improved electronics introduced in SK-IV. The performance of the neutron tagging was recently improved by lowering the neutron tagging threshold (the maximum number of hit PMTs within 10 ns sliding time window) thanks to an additional cut on the continuous dark noise hits after one initial dark noise pulse and new parameters in the neural network. The tagging efficiency was improved from 22% [22] to 25%.

V. EVENT SELECTION
The following selection criteria are applied to separate proton decay signals from atmospheric neutrino background events. The selection criteria resemble those of other recent SK nucleon decay searches [3,6]. The same selections are applied to the data and MC simulation (both proton decay signal and atmospheric neutrino background).
C1: There must be three reconstructed Cherenkov rings.
C2: PID of Cherenkov rings must be consistent with the decay mode. For example, there must be one elike ring and two µ-like rings for the p → e + µ + µ − and p → e − µ + µ + decay modes. SK cannot tell the charge of the final state lepton, so that cuts and backgrounds for the p → µ + e + e − and p → µ − e + e + (p → e + µ + µ − and p → e − µ + µ + ) are essentially the same.
C4: Total mass (M tot ) and momentum (P tot ) of 3-ring events should satisfy 800 < M tot < 1050 MeV/c 2 and P tot < 250 MeV/c. In case of the p → µ + e + e − (p → µ − e + e + ) mode, one additional cut is used. The invariant mass of two e-like rings events should be above 185 MeV/c 2 .
C5: There should be no tagged neutron (only for SK-IV).
These cuts are applied to reduce atmospheric neutrino background, mainly deep inelastic scattering (DIS), based on the kinematics of the outgoing charged particles. Signal selection efficiencies of the C1 cut are lower for decay modes with more muons as shown in FIG.1. This is due to the higher Cherenkov threshold for muons compared to that for electrons. More than 90% of atmospheric neutrino background events are rejected by requiring three rings.
The C3 cut requires a number of Michel electrons depending on the number of muons in the final state. For the p → µ + µ + µ − mode, 2 or 3 decay electrons are required to keep a good signal efficiency. The signal efficiency of this cut depends on the charge and number of muons as the tagging efficiency of Michel electrons for µ + is higher than that for µ − . For example, efficiency for p → µ + e + e − is higher than that for the p → µ − e + e + .
After the C3 cut is applied, the main background for the p → µ + e + e − mode is ν µ charged-current (CC) π 0 production events, in which a π 0 decays to two gamma rays and is identified as two e-like rings. Such background can be reduced by requiring the invariant total mass of two e-like rings to be different from the π 0 mass (C4 cut for one-muon mode). CC π 0 production events are the dominant background to the p → e + e + e − as well (incoming neutrino is ν e in this case). However, since the background rejection is not beneficial to the proton decay signal efficiency, a cut on the invariant mass for two e-like rings is not applied for the p → e + e + e − mode.
Total mass and momentum cuts (C4) are the most effective cut in this analysis. They require that the kinematics of three charged particles is consistent with that from proton decay signals, i.e., their invariant mass should be consistent with the proton mass and the total momentum should be below the upper limit of Fermi motion of protons in oxygen nuclei (the momentum of the proton is 0 for free protons). The lower and higher tails of the proton mass and momentum distributions become larger due to effects of correlated decay.
The probability for a neutron to be generated by deexcitation of a nucleus after proton decay is estimated to be a few percent [23]. On the other hand, neutrons are often generated in the background process, dominated by DIS interactions of atmospheric neutrinos in water. By applying the C5 cut, about half of the background events are rejected while more than 90% of the signal events are kept. The

VI. SYSTEMATIC UNCERTAINTIES
Dominant systematic uncertainties for the signal efficiency are associated with the uncertainties in correlated decay and Fermi motion models. Since the mechanism of correlated decay is not well understood, variation of the signal efficiency was evaluated with 0% and 20% probabilities of correlated decay compared to our nominal estimate of 10%, and the spread was taken as the uncertainty. In the simulation of the signal, Fermi motion is simulated based on the electron-12 C scattering experiment [13]. On the other hand, the Fermi motion model for the background sample is based on the Fermi gas model. This model difference is considered as a source of systematic uncertainty.
For the background, systematic uncertainties on the neutrino flux and cross section models are taken into account in the estimated background rates. These uncertainties are estimated by an event-by-event weighting method based on the neutrino oscillation analysis in SK [20]. A pion generated by a neutrino interaction can interact with a nucleon in oxygen (final state interaction, FSI). It is also possible to interact with other nuclei in water after escaping the original nucleus (secondary interaction, SI). FSI/SI are simulated by a pion cascade model and their uncertainties are taken into account.
Systematic uncertainties for the detector performance and reconstruction are taken into account for both signal and background. In order to estimate these uncertainties, control sample data and MC are compared for each source of systematic uncertainties. We consider uncertainties for FV, detector non-uniformity, energy scale, ring counting, PID, decay electron tagging and neutron tagging. Systematic uncertainty for the detector exposure is negligible. We assigned a 1% error for the detector exposure to be conservative.
Systematic uncertainties for the signal and background are summarized in TABLE III and TABLE IV, respectively. The dominant uncertainties for the background due to the event reconstruction are energy scale and detector non-uniformity. Uncertainties of the energy scale [20] are taken into account for all the charge-related reconstruction parameters. The effect of detector nonuniformity of the energy scale [20] is taken into account for the total momentum reconstruction.
The dominant error for p → µ + µ + µ − (TABLE III) comes from uncertainty of the decay electron tagging. Since the number of candidate events with 3 µ-like rings and 2 or 3 decay electrons (selections for the p → µ + µ + µ − ) is smaller than for the other modes, the statistical error of the atmospheric neutrino control sample data used to estimate the systematic error is larger.

VII. RESULT
No events are found in the signal box region for the p → e + e + e − and p → µ + e + e − (p → µ − e + e + ) modes. In upper column, from left to right for the p → e + e + e − , p → µ + e + e − and p → µ − e + e + modes, respectively. In lower column, from left to right for the p → e + µ + µ − , p → e − µ + µ + and p → µ + µ + µ − modes, respectively. The background MC is normalized by livetime. SK-I-IV are combined in signal MC, background MC and data. Note that both the data and background MC plots for the p → µ + e + e − and p → µ − e + e + modes are the same, since SK cannot identify the change sign of the leptons. For the same reason, the data and background MC plots are the same for the p → e + µ + µ − and p → e − µ + µ + modes as well.
One event is observed in the upper signal box for both the p → e + µ + µ − (p → e − µ + µ + ) and p → µ + µ + µ − modes (all the event displays in [24]). Observed candidates are summarized in TABLE I and TABLE V Two dimensional plots of total mass and momentum for signal (left), background (center) and measured data (right) after the selections C1-C3 and C5 are applied. From the top to the bottom, for the p → e + e + e − , p → µ + e + e − , p → e − µ + µ + and p → µ + µ + µ − modes, respectively. Light blue shows free proton and dark blue shows bound proton in the signal plot. Two black squares show the lower and upper signal boxes. The dot size is enlarged in the signal box only for background. SK-I-IV are combined in signal MC, background MC (4 × 500 years), and data. For the p → µ − e + e + mode the number of signal MC points is lower than that of p → µ + e + e − by 19%, due to the different effective lifetimes (and therefore different decay electron production probabilities) of the µ − and µ + in water. Similarly, the signal MC for p → e + µ + µ − has 20% fewer events than that for the p → e − µ + µ + mode. As in FIG.1, the data and background MC figures are the same for modes that differ only by the charge sign of the leptons.
TABLE I. Summary of signal efficiency, expected background events and data candidates for each decay mode and each data taking period (SK-I to SK-IV). The error values correspond to the statistical uncertainty of the MC sample. Lower and Upper stand for Ptot < 100 MeV/c and 100 < Ptot < 250 MeV/c, respectively. The data events remained in the p → e + µ + µ − and p → e − µ + µ + modes are the same event.
Efficiency (%) Background (events) Candidate (events )  Modes  I  II III IV  I  II  III  IV  I II III     The candidate for the p → e + µ + µ − (p → e − µ + µ + ) modes is observed in the upper signal box of the SK-IV period. Reconstructed proton mass and momentum for this candidate are 882 MeV/c 2 and 160 MeV/c, respectively. Another candidate for the p → µ + µ + µ − mode is found in the upper signal box of the SK-I period. There are two decay electrons, and the total mass and momentum are 835 MeV/c 2 and 170 MeV/c, respectively. These events were visually inspected and they appear not to be mis-reconstruction events.
The expected background events in the p → e + µ + µ − (p → e − µ + µ + ) mode is 0.27±0.04 (stat.) events. Assuming a Poisson distribution of mean 0.27, the probability to observe ≥ 1 events is 18.4%. Considering the expected background events in the p → µ + µ + µ − to be 0.40±0.07 (stat.) events, the probability of observing ≥ 1 events is 25.8%. The observed events are consistent with expected backgrounds; therefore lower proton lifetime limits at 90% confidence level (CL) with respect to each proton decay Additionally a Ptot < 250 MeV/c cut is applied on the total mass plot and a 800 < Mtot < 1050 MeV/c 2 cut is applied on the total momentum plot. From the top to the bottom, the p → e + e + e − , p → µ + e + e − , p → e − µ + µ + and p → µ + µ + µ − modes are shown, respectively. Dotted black lines show the boundary of the signal box. SK-I-IV are combined in signal, background MC and data. Background MC is normalized by atmospheric neutrino flux, oscillation probability and livetime. Signal MC is normalized to the partial proton lifetime limit calculated in section VIII. The number of signal MC events for the p → µ − e + e + mode is lower than that of p → µ + e + e − by 19%, due to the different effective lifetimes (and therefore different decay electron production probabilities) of the µ − and µ + in water. Similarly, the signal MC for p → e + µ + µ − has 20% fewer events than that for the p → e − µ + µ + mode. As in FIG.1, the data and background MC figures are the same for modes that differ only by the charge sign of the leptons. mode are calculated by using a Bayesian method [25,26]. The limit calculation is the same as for recent nucleon decay analyses [3,4,6]. We have 8 signal regions (4 periods × 2 boxes) for each decay mode. The probability density function (PDF) is defined for each region as below.
Here, i is the index of each signal region, A i is a normalization factor, Γ is the decay rate, n i is the observed events, λ i is the exposure, i is the signal efficiency and b i is the expected background events. P (Γ) represents the probability distribution for the decay rate, assumed to be uniform. P (λ i ) and P ( i ) are the probabilities for the exposure and signal efficiency, respectively, described by a Gaussian. P (b i ) is the probability for the expected background defined by the convolution of Gaussian and Poisson distributions. All PDFs are combined and the upper limits of the decay rate Γ limit at 90% CL are estimated as follows.
Finally the lower limit of partial proton lifetime is calculated according to Here B is the branching ratio of each proton decay mode. By using these functions, the partial lifetime limits at 90% CL for each mode of proton decay into three charged leptons are calculated as summarized in TABLE V.

IX. CONCLUSION
Proton decay into three charged leptons was searched for by using 0.37 Mton·years of data collected by SK. The observed data were consistent with the atmospheric neutrino background prediction and no clear indication of proton decay was observed. According to the observation, the model [8] for these decay modes at an energy scale below 100 TeV was excluded by this analysis. The lower partial lifetime limits at 90% CL were calculated for each mode as summarized in TABLE V. Compared with the previous limits by IMB-3 and HPW experiments, each limit was improved by 15-1800 times in this analysis as shown in FIG.4. A first limit has been set for the p → µ − e + e + mode.