Neutron-antineutron oscillation search using a 0.37 megaton-years exposure of Super-Kamiokande

As a baryon number violating process with Δ B ¼ 2 , neutron-antineutron oscillation ( n → ¯ n ) provides a unique test of baryon number conservation. We have performed a search for n → ¯ n oscillation with bound neutrons in Super-Kamiokande, with the full dataset from its first four run periods, representing an exposure of 0 . 37 Mton-years. The search used a multivariate analysis trained on simulated n → ¯ n events and atmospheric neutrino backgrounds and resulted in 11 candidate events with an expected background of 9.3 events. In the absence of statistically significant excess, we derived a lower limit on ¯ n appearance * Deceased.


I. INTRODUCTION
The present baryon asymmetry of the universe provides indirect evidence for baryon number violating (BNV) processes [1], which cannot be sufficiently explained by mechanisms within the Standard Model (SM) [2]. Searches for BNV processes probe physics beyond the reach of the SM can be classified based on the baryon number violation (ΔB) involved. Processes with ΔB ¼ 1 are tightly constrained by null observations from proton decay searches, and processes with ΔB ¼ 3 are expected to conflict with nucleosynthesis scenarios [3]. The Standard Model allows for non-perturbative processes involving sphalerons that would wash out any baryon number asymmetry from processes that conserve B − L, where L is lepton number, before the electroweak phase transition [4]. Therefore, as a BNV process violating both B and B − L, neutron-antineutron oscillation provides a unique probe of baryon number violation and essential insight into the baryon asymmetry and baryogenesis.
Several models predicting n −n oscillations have been proposed since the 1970s, including those employing an SUð2Þ L ×SUð2Þ R ×SUð4Þ c gauge group to generate a baryon asymmetry [5,6] and others that propagate SM fields into extra space-time dimensions [7]. The predicted oscillation times vary from 10 9 s [7] to 5 × 10 10 s [6] and correspond to energy scales of 10 2 -10 3 TeV, well above the scale that can currently be probed by accelerators [8]. Recent developments from a lattice QCD calculation on the postsphaleron baryogenesis model [9] and a left-rightsymmetric model with extra dimensions [10] indicate that n −n oscillation is more sensitive to baryon number violation than previously expected, providing further motivation for experimental searches.
The probability of a free neutron oscillating to an antineutron can be parametrized as a simple 2 × 2 Hamiltonian and can be written as where ΔE is the energy difference between the neutron and antineutron and δm ¼ 1=τ n−n , where τ is the neutronantineutron oscillation time. In the case of degenerate neutron and antineutron energies, Eq. (1) has a simplified form, P n→n ðtÞ ≈ ðδmtÞ 2 ¼ t τ n−n 2 : For bound neutrons in nuclei, the probability can be written as [11] P nuc ðn →nÞ¼ where T nuc is the observed neutron lifetime in neutronantineutron oscillation, and R is the so-called nuclear suppression factor that accounts for the suppression of oscillations due to differences in the nuclear potentials of neutrons and antineutrons. Theoretical calculations of R using effective field theories vary [12,13], but, in the following, we adopt R ¼ 0.517 × 10 23 s −1 for 16 O based calculations by Friedman et al. [11]. Experimental searches for n −n oscillation rely on observing particles (mostly pions) produced when a neutron oscillates into an antineutron and annihilates with a nearby nucleon. There have been a number of n −n searches using either free neutrons [14] or bound neutrons [15][16][17][18][19][20][21], none of which have yielded a positive signal. Accordingly, constraints on the n −n oscillation time have been set at τ n−n > 0.86 × 10 8 s for free neutron oscillation [14] and at τ n−n > 2.7 × 10 8 s for bound neutrons [15].
In this paper, we present a search for n −n oscillations using the full dataset from the first four running periods of Super-Kamiokande and update the result presented in Ref. [15], which used data from the first period. The current analysis includes an updated dataset, an updated hadron production model, final state interactions, and adopts a multivariate method to achieve better discrimination between the background and signal processes. This paper is organized as follows. After a short description of the Super-Kamiokande detector in Sec. II, we describe the simulation of both the n −n signal and atmospheric neutrino background in Sec. III. The selection algorithm and analysis cuts are explained in Sec. IV, followed by discussion of systematic uncertainties in Sec. V. Analysis results and concluding remarks are presented in Secs. VI and VII, respectively.

II. THE SUPER-KAMIOKANDE EXPERIMENT
Super-Kamiokande (SK) is a cylindrical 50 kiloton water Cherenkov detector located in Kamioka, Japan, that is shielded by a 2700 meter water-equivalent rock overburden [22]. The detector consists of an outer detector instrumented with 1885 outward-facing 8-inch Photomultiplier tubes (PMTs) mounted 2 m from the detector's outer wall on a structure that optically separates it from the inner detector (ID). This structure also supports the 11,129 inward-facing 20-inch PMTs that form the ID and view its 32 kton target volume. The outer detector is primarily used as a veto for charged particles entering from outside the detector or identifying particles that exit the ID, and the ID itself is used reconstruct the energies, vertexes, and particle types of most interest to the present work.
The experiment started data taking in 1996 and underwent four data-taking phases since then labeled as SK-I, II, III, and IV. The SK-I period ran from 1996 until the detector underwent maintenance in 2001. During that period, an accident destroyed more than half of the SK PMTs, reducing the photocathode coverage from ∼40% to ∼19% for the SK-II period in 2002-2005. After replacing the missing PMTs in 2005, the detector restarted operations as SK-III in 2006-2008. Following upgrades of the frontend electronics and water purification system, the SK-IV period ran from 2008 until May of 2018, when the data taking was paused and the detector tank was opened for further upgrades. The analysis in this work uses the full dataset from the SK-I-IV periods. Details of the detector and its calibration can be found in [23].

III. SIMULATION
Following the oscillation of a neutron into an antineutron, the subsequent annihilation of the antineutron with a nucleon in the oxygen nucleus is expected to produce many visible particles, most of which are pions. The simulation of this signal is broken into stages: oscillation, hadronization, final state interactions of particles before exiting the nucleus, and finally propagation and subsequent reinteraction of those particles with detector media. During the first stage, the position of the oscillated neutron within the nucleus is determined using the standard Woods-Saxon distribution [11,24] with a Fermi momentum simulation based on the spectral function measured in [25]. The effect of nuclear binding energy is taken into account by subtracting it from the nucleon masses when calculating the annihilation products, using 39.0 MeV for s-state and 15.5 MeV for p-state nucleons, respectively. Thereafter the oscillated antineutron is assumed to have an equal probability of annihilating with any remaining nucleons.
Modeling of thenn ornp annihilation products is done based on available accelerator data. Due to a lack of antineutron scattering data, the hadronization simulation uses results from antiproton scattering experiments instead. Assuming isospin symmetry, we used data from thepp annihilation experiment Crystal Barrel [26,27] to simulate thenn annihilation. For thenp channel, we used thepn annihilation branching ratio measurements from the OBELIX experiment [28] and bubble chamber data [29][30][31] and then flipped the signs of the charged pions to matchnp. Tables I and II show the branching ratios for nn andnp adopted in the simulation. The branching ratios of kaonic channels are artificially constructed due to lack of experimental data, and the kaonic production fornp is less than 1=2 fromnn, and thus is omitted. Corresponding uncertainty calculations can be found in Sec. V, and the efficiency calculation is explained in Sec. IV.
Hadronization products are mostly pions. The pion interaction probability within the oxygen nucleus is expected to be large, and these so-called final state interactions (FSI) include quasielastic scattering (e.g., π þ n → π þ n), absorption (π þ þ n → p, π − þ p → n), charge exchange (π þ þ n → p þ π 0 , π − þ p → n þ π 0 ), and pion production (π AE þ n → π AE þ n þ π 0 ) [32].F o r pions above 500 MeV=c, the surrounding nucleons are treated as quasifree particles, while for lower momentum pions the interaction probabilities are calculated according to the model of Salcedo and Oset [33] in consideration of the effect of Pauli blocking. More details can be found in Ref. [32]. Atmospheric neutrino interactions in water are the dominant background to the search for n −n oscillation at SK. The theoretical calculation from the HKKM model [34,35] predicts the atmospheric neutrino flux at Kamioka in the energy region from sub-GeV up to several TeV after oscillation. Using this flux prediction, we simulated atmospheric neutrino interactions, including the outgoing particles and their subsequent interactions with the nuclear medium in water, with NEUT version 5.3.6 [36]. Final state interactions for both the signal and background are simulated with NEUT. Particles escaping the nucleus are passed to a GEANT3based [37] detector simulation. The simulation tracks particles through the detector medium, simulating their interactions in water as well as the production of secondary particles and the response of the PMTs to Cherenkov radiation. Detailed tuning and calibration has been performed to provide a tailored simulation of photon propagation in SK [38]. The interaction of hadrons with water is simulated using the GCALOR package [39], except for pions below 500 MeV=c, which are simulated using a model based on NEUT's FSI simulation. The final background is reweighted to the result of the analysis in [40], adjusting its central value to the best fit oscillation and systematic error parameters favored by the SK data.

IV. EVENT RECONSTRUCTION AND SELECTION
The present analysis uses the full dataset from the SK-I through SK-IV periods, corresponding to 6050.0 live-days. Events are required to be fully contained (FC), meaning the number of PMTs in the highest charge cluster of outer detector hits is less than 10 in SK-I and less than 16 in SK-II-III-IV. Timing information in each event's hit PMTs in the ID is used to reconstruct an overall vertex from the event, from which an iterative search based on the Hough transform [41] is performed to identify Cherenkov rings. Each Cherenkov ring is classified according to its hit pattern and opening angle as either "showering" (e-like) for particles that create electromagnetic showers such as e and γ or as "nonshowering" (μ-like) for particles such as μ and π AE . The momentum of each ring is determined by the particle type and the charge among all hit PMTs within a 70°cone around the ring with consideration of charge shared between multiple rings. An additional search for delayed electrons from muon decays is performed from 1.2-20 μs after the primary event trigger.
This analysis starts with FC events more than 2.0 m from the ID wall, which defines a 22.5 kton fiducial volume. The reduction efficiency is 92% for n →n signal events in fiducial volume. This sample is then processed in two stages, first applying simple analysis cuts before applying a multivariate technique to extract the signal.

A. Analysis precuts
Based on the distinct features of n −n and atmospheric neutrino events, several preliminary cuts are applied to reduce background rates while maintaining high signal efficiency. The n −n oscillation signal is expected to have multiple pions, while a large number of atmospheric neutrino interactions are elastic scatters with only one Cherenkov ring from the outgoing charged lepton. Therefore, the number of reconstructed rings is required to be >1. This cut removes ∼75% of the background while keeping 89% of the signal. Unlike the wide range of energies covered by atmospheric neutrinos, the n −n signal is more kinetically constrained, and thus a set of kinematic cuts are also applied. Here, the total reconstructed momentum is required to be within ½35; 875 MeV=c, the visible energy in [30,1830] MeV, and the total reconstructed invariant mass in ½80; 1910 MeV=c 2 . After the cut on the number of rings, these kinematic cuts further remove ∼50% of the background with a relative signal efficiency of 98%.

B. Multivariate analysis
Event displays of a simulated n →n signal event and a simulated background event are shown in Fig. 1.D u e to the high ring multiplicity, the performance of ring reconstruction for n →n signal events is not as satisfactory as typical sub-GeV neutrino events. To compensate for the limitation of ring reconstruction and to include more discriminant features, we applied a multivariate analysis (MVA) to events passing the precuts. Compared to a conventional box-cut analysis [15], this analysis significantly enhance the separation between n →n signal and background. An estimation using the same Monte Carlo (MC) set shows that the sensitivity of the MVA method is twice that of the box-cut method.
Compared to atmospheric neutrino backgrounds,nn or np annihilation within oxygen are generally expected to be more constrained kinematically and have more Cherenkov rings isotropically distributed in the detector. To exploit these features, we introduced 12 variables into the MVA, among which three are conventional kinematic quantities, including the visible energy, total momentum, and total invariant mass.
The remaining nine input variables are as follows. Since only a fraction of atmospheric neutrinos has sufficient energy to produce multiple charged particles, signal events are typically expected to have more visible Cherenkov rings. The number of such rings is used as a variable. However, the full reconstruction is limited to five rings, as in the case of Fig. 1. Therefore, an additional variable that counts ring fragments, or potential rings, is also introduced.
The total momentum of an n −n event is limited by the momenta of the interacting nucleons, while a background event can carry more momentum from the incident neutrino and is expected to be more forward going at the energies needed to produce multiple particles. Therefore, this search employs four variables to quantify the isotropy of candidate events. The energy ring ratio is defined as ðE tot − E max Þ=½E tot · ðn ring − 1Þ, where E max is the energy of the ring with highest energy in an event, E tot is the total energy of the event, and n ring is the number of rings. For the n −n signal, the annihilation energy is more uniformly distributed among the outgoing pions and therefore, the distribution of this variable is expected to have a sharper peak than of backgrounds. Signal events are also expected to have higher sphericity than backgrounds, so this analysis adopts a sphericity variable [42]. Fox-Wolfram moments, which are superpositions of spherical harmonics that measure correlations between particle momenta (see Ref. [43] for details) are also adopted to describe the correlation between rings. This analysis employs the first and second order Fox-Wolfram moments, since higher orders were found to provide little extra discrimination ability.
Finally, three variables related to particle identification are used: the number of e-like rings, the number of decay electrons, and the maximum distance to any decay electron from the primary vertex. Due to the large number of signal modes with one or more π 0 s in the final state, signal events are expected to have more e-like rings from their decays into photons. Corresponding distributions for signal and background MC set after the precuts are shown in Fig. 2.
These 12 variables are used in the construction of a multilayer perceptron (MLP) [44], which is trained on n −n signal and atmospheric neutrino background MC. The MLP consists of a network of layers of nodes that are weighted and interconnected in order to optimize the discrimination between event types. Input variables form the input layer nodes and are combined in the MVA into a single node at the output layer, which is the estimator describing how signal-or backgroundlike an event is. Between these layers there can be so-called hidden layers, whose structure and connectivity can be altered to optimize performance. In this analysis, a trial-and-error optimization for the hyperparameters of the MLP structure was performed and the final structure was determined to be 1 hidden layer with 18 hidden nodes.
The signal efficiency and background efficiency as a function of the estimator value is shown in Fig. 3, where 0 corresponds to backgroundlike and 1 is signal-like. A sensitivity analysis was performed assuming a 0.37 megaton-years exposure and realistic systematic errors (described below) using the Rolke method [45] to determine the optimal cut position in the output estimator. The optimized cut was found to be 0.789, where the signal (background) efficiency from the MVA alone is 5.0% (0.1%). Combined with the preselection efficiency, the total signal efficiency is 4.1% with an expected background of 0.56 events per year, or 9.3 events over the entire data period. Selection efficiencies for each of the signal channels can be found in the last column of Tables I and II. Among the multiple types of neutrino interactions, the dominant background in this analysis is from deep inelastic scattering (DIS), with secondary contributions from FIG. 1. The event display of a simulated n →n oscillation event (top) and a simulated background event (bottom). The top figure shows thenp annihilation producing six pions. In total, there are 13 rings of these pions and their decay products, among which five were reconstructed (colored rings). The beige ring is a successfully reconstructed π þ . The dashed small ring is a π − mistakenly reconstructed as a electron-like particle, and the dashed large ring is a similarly misreconstructed π þ . The two solid cyan rings are 2γ s from the decay of a π 0 . The bottom figure shows the deep inelastic scattering process of a ν e . In this event, there are four rings reconstructed by the algorithm, all of which are γs. charged current pion production (CC 1π), neutral current pion production (NC 1π), and charged current elastic scattering (CC EL). Figure 4 shows the remaining backgrounds before the MVA cut. After applying the MVA cut, the remaining backgrounds in the final sample are shown in Table III, with ν μ þν μ contributing 5.8 events, and ν e þν e contributing 3.5 events. The contribution from ν τ þν τ is less than 0.1 event and is not shown in Table III.

V. SYSTEMATIC ESTIMATION
In this analysis systematic uncertainties are separated into two categories, those that arise from uncertainty in the physics modeling, such as the hadronization process and final state interaction, and those related to the detector response and event reconstruction.

Signal
Uncertainty in the momentum of the oxygen nucleons is expected to impact the resulting momentum of the n −n annihilation products. A systematic uncertainty is derived from the difference between the default spectral function model (described in Sec. III) and the Fermi gas model [46] used in the atmospheric neutrino simulation. It yields an uncertainty in the signal efficiency of 7%. Our model does not account for the recent measurement of higher momentum distribution, as summarized in Ref. [47]. With the dominance of FSI systematic uncertainties, the influence is estimated to be small.
Measured uncertainties in the branching fraction of each annihilation channel also introduce a systematic uncertainty in the hadronization process, resulting in uncertainties in the pion multiplicity of signal events. This uncertainty is accounted for by assigning uncertainties on the branching ratio of each channel listed in Table I  and Table II based on the statistical uncertainty in the results from the Crystal Barrel [26,27] and the OBELIX experiments [28]. They were then propagated to the analysis by reweighting the various final states accordingly and result in a 4% uncertainty on the signal efficiency.
Final state interaction modeling is the dominant systematic error on the signal efficiency. To estimate this uncertainty, we generated separate MC sets, each with different FSI model parameters that control the strength of the interaction cross sections and are allowed by fits to pion-nucleon scattering data [20]. These MC samples were processed through the same event selection, and the largest change in the signal efficiency is taken as the uncertainty. In this analysis, the largest deviation came from a variation with enhanced quasielastic scattering and absorption, but with decreased inelastic scattering, which produces fewer hadrons and thus lower efficiency. The assigned uncertainty is 31%.

Background
Uncertainties on the atmospheric neutrino background were calculated using the fit result from the SK atmospheric neutrino analysis [40]. A set of weights was constructed for each event, describing how it changes under a 1σ variation of each systematic error parameters used in that analysis. Applying these weights to the MC set allows the uncertainty to be conservatively propagated to the background prediction.
The overall atmospheric neutrino flux normalization has an uncertainty of 15% [40] in the dominant background energy range between 1 and 10 GeV, resulting in a 7% uncertainty in the background rate. In total, the uncertainty introduced by modeling of the flux was estimated to be 8%. Neutrino oscillation parameter uncertainties, particularly from θ 23 , also introduce a 3% uncertainty. Uncertainties from the neutrino interaction modeling are the most significant contribution to the error budget. The total uncertainty from neutrino interaction was estimated at 24%, among which the main contribution was found to originate from uncertainties in the deep inelastic scattering model and its cross section.

B. Detector systematics
Uncertainties in the detector's energy scale and the reconstruction's ability to accurately identify the number of and particle type of each ring introduce uncertainties in both the signal efficiency and background rate. The energy scale uncertainty is evaluated using calibration sources and control samples, such as cosmic ray muons and their decay electrons [40], and is 3.3% in SK-I, 2.8% in SK-II, 2.4% in SK-III, and 2.1% in SK-IV. It results in a 5% and 11% uncertainty on the signal efficiency and background rate, respectively. Similarly, differences in the water quality in the top and bottom regions of the SK tank introduces an asymmetry in the energy scale that introduces an additional 4% signal efficiency uncertainty and 6% background rate uncertainty.
Ring counting introduces 2% uncertainty in signal efficiency and 1% in background rate. This uncertainty is estimated by comparing the ring counting likelihood distribution of MC sets and a controlled sample of data [48].
For MVA variables besides ring counting, energy scale, and nonuniformity, we use an inclusive controlled sample (FC data after precuts, before MVA), and compare data and MC prediction, as shown in Fig. 2. The source uncertainties are assigned from the deviation of data and MC sets. These source uncertainties are then propagated to efficiency uncertainties. The individual systematic sources and their uncertainties are summarized in Table IV, while the total efficiency and uncertainty are presented in Table V.

VI. RESULT
This full SK-I-IV dataset corresponds to an exposure of 0.37 megaton-years. After applying the cuts above, 11 events are found in data, which is consistent with the expected background of 9.3 AE 2.7 events. Furthermore, data and MC sets are in good agreement both before (Fig. 5) and after (Fig. 6) the MVA cut. The input variables to the MVA show a similar agreement (Fig. 2). Accordingly, we find no evidence for neutron-antineutron oscillations. Figure 7 shows the 11 candidate events within the detector. The spatial distribution is uniform as expected. Figure 8 shows the distribution in time. The dependence of  events after precut on time is due to the live time of SK. Performing a K-S test on the distrbution yields a maximum distance between data and MC sets at 0.33. To determine the likelihood of this result, this procedure was repeated on simulated datasets with the same size as the observation and assuming a constant rate. Among these pseudoexperiments, 14% had a K-S distance larger than 0.33, indicating no significant deviation from the assumed uniform distribution and is consistent with the expectation from atmospheric backgrounds. A comparison of the expected atmospheric neutrino background, signal efficiency and observed data in each SK run period is shown in Table VI. Signal efficiencies and background rates are slightly different across these run periods, and the majority of candidate events are found in SK-IV, which has the longest live time. The Poisson probability of observing nine or more events in SK-IV with an expectation of 5.5 events (ignoring systematic uncertainties) is 10.6%. This observation is similarly consistent with the background expectation.
In the absence of a statistically significant excess in data, a lower limit is established. To account for both statistical and systematic uncertainties, we used Rolke method in confidence interval calculation. The observation limit on neutron lifetime is set at 3.6 × 10 32 years (90% C.L.). Equation (3) and R ¼ 0.517 × 10 23 s −1 [11] are used to derive the corresponding limit on the n −n oscillation time, τ n→n > 4.7 × 10 8 s. A comparison between the expected sensitivity and this result is shown in Table VII. Alternative calculations of the nuclear suppression factor R can be found in Refs. [12,13]. Table VIII compares the present results with those from other bound neutron experiments and free neutron oscillation experiments. Papers before the year 2000 typically report τ n→n assuming R ¼ 1 × 10 23 =s, and the previous SK result considered uncertainty in the theoretical prediction of R. For better comparison and easier conversion, τ n→n is presented as ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi T n−n =R p with corresponding nuclear suppression factor R listed in Table VIII. The limit from free neutron experiments is also listed for reference, but should not be directly compared with the limits from bound neutron experiments due to the medium modifications of the six-quark operator which is still a subject for theoretical interpretation [49]. Among the bound neutron experimental searches so far, this analysis gives the most stringent limit on n →n oscillation.

VII. CONCLUSION
We performed a n −n oscillation search with SK-I-IV data using a multivariate analysis. Compared to previous results [15], the updated final state interaction model predicts fewer pions and less separation between signal and neutrino backgrounds. With the advanced MVA method and the inclusion of multiple new variables, the sensitivity of this analysis is still greatly enhanced.
For the 0.37 megaton-year exposure at SK, we observed 11 events with an expected background of 9.3 AE 2.7 events. There is no statistically significant excess of data events, so a lower limit on the neutron lifetime is set at 3.6 × 10 32 years at 90% C.L., corresponding to a lower limit on the neutron-antineutron oscillation time in 16 Oo f τ n→n > 4.7 × 10 8 s. This is the world's most stringent limit on neutron-antineutron oscillation so far, with 90% improvement from the previous best limit [15], and is reaching the predicted parameter space of some theoretical models [6].