Search for invisible modes of nucleon decay in water with the SNO+ detector

This paper reports results from a search for nucleon decay through 'invisible' modes, where no visible energy is directly deposited during the decay itself, during the initial water phase of SNO+. However, such decays within the oxygen nucleus would produce an excited daughter that would subsequently de-excite, often emitting detectable gamma rays. A search for such gamma rays yields limits of $2.5 \times 10^{29}$ y at 90% Bayesian credibility level (with a prior uniform in rate) for the partial lifetime of the neutron, and $3.6 \times 10^{29}$ y for the partial lifetime of the proton, the latter a 70% improvement on the previous limit from SNO. We also present partial lifetime limits for invisible dinucleon modes of $1.3\times 10^{28}$ y for $nn$, $2.6\times 10^{28}$ y for $pn$ and $4.7\times 10^{28}$ y for $pp$, an improvement over existing limits by close to three orders of magnitude for the latter two.


I. INTRODUCTION
Violation of baryon number conservation is predicted by many Grand Unified Theories [1] potentially explaining the matter-antimatter asymmetry of the universe.Searches for the decay of protons or bound neutrons act as important constraints on our understanding of physics beyond the Standard Model.
Modes of nucleon decay involving visible energy deposition by decay products, such as positrons, pions or kaons, have so far not been observed by large scale detectors [2][3][4].As such, interest has turned to less wellconstrained decay modes that could have escaped detection.A potential set of channels are the invisible decay modes in which one or more nucleon decays to final states which go undetected, such as those with neutrinos, other exotic neutral, weakly-interacting particles or charged particles with velocities below the Cherenkov threshold.Although no prompt signal would be observed from the particles directly emitted in the decay, the remaining nucleus would be left in an excited state, and could then emit a detectable signal as it de-excites.The search for the de-excitation signal of the final state nucleus is model-independent, as it puts no requirements on the particles produced in the decay.Theoretical models include modes with decays to three neutrinos [5] or to non-Standard Model particles such as the unparticle [6] or dark fermions [7], the latter providing a possible solution to the neutron lifetime puzzle [8].
The Sudbury Neutrino Observatory (SNO) and Kam-LAND experiments have conducted searches for such model-independent modes with KamLAND setting the current best limit for the invisible neutron decay life-time of τ (n → inv.) > 5.8 × 10 29 y at 90% C.L. [9] and SNO setting the best limit for invisible proton decays of 2.1 × 10 29 y [10], improving on previous limits by the Borexino Counting Test Facility (CTF) [11] and Kamiokande [12].Limits also exist for the dinucleon modes of 1.4×10 30 y [9] for the nn mode from Kam-LAND, 5.0×10 25 y [11] for the pp mode by the Borexino CTF and 2.1×10 25 y [13] for the pn mode.
The SNO + experiment has been running since December 2016, taking commissioning data with the detector filled with ultrapure water.During this phase, a new search has been made for invisible nucleon decay via the decay of 16 O.After an invisible nucleon decay, the 16 O is left in either the 15 O* excited state, if the decaying nucleon was a neutron, or in the 15 N* state, if the decaying nucleon was a proton.44% of the time, 15 O* will de-excite to produce a 6.18 MeV gamma, and 2% of the time, the decay will produce a 7.03 MeV gamma.Similarly, 15 N* will produce a 6.32 MeV gamma in 41% of decays, with 7.01, 7.03 and 9.93 MeV gammas produced 2, 2 and 3% of the time, respectively [14].
This search has a unique sensitivity for two reasons: firstly, the branching fraction to produce a visible signal of a de-exciting oxygen nucleus is larger than the 5.8% for carbon [16] used by KamLAND.Secondly, the use of H 2 O rather than heavy-water (D 2 O) removes the solar neutrino charged-current and neutral-current signals, major backgrounds in the SNO search.
Dinucleon modes are also sought, based on the emission of de-excitation gamma rays from 14 O * , 14 N * and 14 C * for the nn, pn and pp invisible decay modes, respectively [15].The pp decay can lead to gammas of 6.09 MeV at 10.9% and 7.01 MeV at 20.1%, and the pn 1. Spectrum of de-excitation gamma rays emitted after the invisible decay of 16 O for single nucleon decay [14] (left) and dinucleon decay [15] (right) as a function of gamma ray energy Eγ between 5 and 9 MeV.decay has a 6.45 MeV gamma with 7.7% and 7.03 MeV gamma with 8.9% probability.The nn decay proceeds via many channels, with a summed branching ratio of 4.53% for gamma emission between 5 and 9 MeV.The branching ratio for single and dinucleon decay are shown in Fig. 1.

II. THE SNO + DETECTOR
The SNO + detector is inherited from SNO [17], with several major upgrades to enable the use of liquid scintillator as the primary target rather than D 2 O.The detector consists of a 6 m-radius spherical acrylic vessel (AV) surrounded by 9394 inward-facing photomultiplier tubes (PMTs) mounted on a stainless steel support structure with an 8 m radius from the center of the AV.During the initial water phase, the AV was filled with approximately 905 tonnes of ultrapure water.The cavity where the detector is installed is also filled with ultrapure water to shield the innermost regions from radioactivity in the PMTs and the surrounding rock.Among the upgrades for SNO + was a new rope net to counteract the buoyancy of the scintillator and hold down the AV [18].The PMTs and front-end electronics have been reused, with work done to repair PMTs that failed over the lifetime of SNO.Aspects of the trigger system and data acquisition (DAQ) software were upgraded to handle the higher data rates and light yield expected in the scintillator phase.

III. DATA SELECTION
The results reported in this paper are based on the analysis of 235 days of data recorded between May 4th 2017 and December 25th 2017.During this time, the de-  tector was live for 95% of the time, with 16.9% of that spent performing calibration or maintenance.A series of data quality checks were made to select the data for analysis with specific selection criteria for the detector state, event rate, occupancy, and number of poorly calibrated channels, resulting in the rejection of 29.3% of the live time.The removal of time-correlated instrumental effects, cosmic ray muon events, and trigger dead time between events resulted in the loss of an additional 2.4% of live time.The final analyzed data set has a live time of 114.7 days with an uncertainty of 0.04%.
During the SNO + water phase, significant work was done on commissioning the water processing and recirculation systems.Changing background levels associated with these activities motivated a time-dependent analysis.The data were split into six data sets, during each of which the background levels were relatively stable, each with its own background estimate and set of analysis cuts.Table I details the live time of each data set.
Channels that failed calibration checks were excluded from the analysis, though they still contributed to the hardware trigger.The number of offline channels varied over time but on average was around 800 channels.A stable and well-calibrated channel can still register hits caused by electronic cross-talk and PMT dark noise.Hitcleaning algorithms, used to exclude cross-talk hits from the analysis, typically remove approximately 2% of hits in an event.The dark noise is measured and simulated on a run-by-run basis.

IV. EVENT RECONSTRUCTION
A. Vertex reconstruction SNO + uses a set of algorithms to reconstruct the position and direction of Cherenkov events based upon maximizing the likelihood of a distribution of PMT hit times that have been corrected for the time residuals, i.e., the time of flight relative to the assumed vertex position, and of the angle between the true direction and the line from the trial position to the PMT.These algorithms only consider hits on well-calibrated online channels, within 50 ns of the modal PMT hit time.
Three additional observables were used to classify events.The In Time Ratio (ITR) is the ratio of PMT hit time residuals within a prompt timing window of [−2.5 ns, 5 ns] to the total PMT hit time residuals in an event.β 14 is an isotropy parameterization based on the 1st and 4th Legendre polynomials of the distribution of PMT hits in the event [19], calculated using the angle between two PMTs with respect to the reconstructed position.The projection of a particle's reconstructed direction unit vector onto the corresponding event position unit vector u • r determines whether the particle appears inward-or outward-going relative to the center of the AV.For the physics and calibration analyses, it is required that ITR>0.55 and −0.12 <β 14 <0.95.

B. Vertex systematic uncertainties
Systematic uncertainties in the reconstructed vertex were evaluated using a 16 N calibration source [20], previously used in the SNO experiment, to produce tagged 6.1 MeV gammas.
A series of 16 N scans were taken during the data taking period in 2017.During a scan, 16 N data was collected in a series of runs at points spaced about 50 cm apart along the principal axes of the detector, typically through the center along the x, y and z axes, where the z axis points upwards through the neck of the AV.Additional scans were also taken off-axis in the xz and yz planes to evaluate asymmetries in the detector and reconstruction.

Position uncertainties
To evaluate uncertainties associated with the reconstructed event vertex position and direction, the measured detector response to the 16 N calibration source was compared with predictions from Monte Carlo simulation, shown in Fig. 2.
For events which were tagged by the source PMT and passed the β 14 and ITR cuts, the difference in the reconstructed vertex position and source position was taken in each of the x, y and z axes.The resulting one-dimensional distributions were fit with a distribution function representing the position of the first Compton electron, esti- mated from the Monte Carlo model, convolved with a Gaussian function and an exponential tail.
The uncertainties in reconstruction were characterized in terms of a constant offset between the position of the source and the mean reconstructed position, x→ x + δ i , a position-dependent scale factor in which the position offset scales linearly with its value, x → (1 + δi 100 )x, and a resolution describing the width of the distribution of reconstructed event positions, x → x + R(0, δ i ), where R is a random number drawn from a Gaussian distribution of width δ i and mean 0. The vertex offset was calculated as the volumeweighted mean of the difference between the Gaussian fitted means of data and Monte Carlo simulation while the scale was found as the slope of a linear fit to the differences, both listed in Table II.These were applied during signal extraction by shifting and scaling the position of each event according to the uncertainties along each axis independently and recomputing the timing probability density functions (PDFs) used for signal extraction.
The position resolution of the data events was found to be ∼200 mm.The difference in resolutions between the data and Monte Carlo events was modeled as a Gaussian of standard deviation δ i = σ 2 Data − σ 2 MC for each 16 N run.The results were combined in a volume-weighted average, independently for each detector axis.The resulting values for δ i are also listed in Table II.These were applied during the signal extraction, smearing the positions of all Monte Carlo events by a Gaussian distribution of the appropriate width to reproduce the data.

Angular resolution
Reconstructed direction is also evaluated using the 16 N source, by taking into account the high degree of colinearity between Compton scattered electrons and the initial gamma direction.The angle θ between the 'true' direction, taken to be the vector from the source position to the reconstructed event vertex, and the reconstructed event direction was calculated and the distribution of these angles were compared for data and Monte Carlo events.To reduce the effect of position reconstruction uncertainties, only events that reconstructed more than 1200 mm from the source position were used.The resulting distributions were fit with a functional form of two exponentials as employed by SNO [21]: where β m and β s are the slopes of the two exponential components, parameterizing the distribution for small and large angles respectively, and α m is the fraction of the events following the exponential with slope β m .The volume-weighted average of the differences in β was computed across all runs and the resulting value used as an estimate of the uncertainty in angular resolution, transforming cos θ → 1 + (cos θ − 1)(1 + δ i ).Events that were transformed beyond −1 were randomly assigned a value between −1 and 1.The resulting uncertainties are listed in Table II.

β14 uncertainties
The systematic uncertainty of β 14 , shown in Table II, was determined from the difference between the means of Gaussian fits to data and Monte Carlo simulations of a 16 N run with the source at the center of the detector.This found a shift of −0.031 ± 0.004, which was applied to the Monte Carlo predictions.

C. Optical calibration
The detector's optical parameters, including the attenuation and scattering of light in the water and acrylic, and the PMT angular response, are based on calibration measurements of the same materials from SNO [19].Further optical calibrations were carried out using the 'laserball' [22], a multi-wavelength laser pulse diffuser capable of running with 337, 365, 385, 420, 450 and 500 nm wavelengths, deployed within the detector.
Using the laserball data, the group velocity of light in the SNO + water was verified to be consistent with the values used in the SNO + simulation and reconstruction [23].
Laserball runs were taken along the z axis after major detector maintenance periods to re-measure the PMT gain, electronic time delays and their dependence with integrated pulse charge.The delays were compared for the different laserball runs, with time drifts typically less than 0.5 ns.Larger observed drifts are consistent with actual changes in the detector hardware, for example, the replacement of offline channels.
From visual observation during detector upgrades, it is known that the reflectivity of the PMT concentrators has degraded over time.The concentrator diffuse reflectivity was tuned in simulation by comparing the PMT hit time residual spectrum between a central 16 N run and its Monte Carlo simulation, with particular attention given to peaks in the late light (with residual times between 47 and 80 ns) due to reflections from the concentrators, PMT glass, and the AV.The total reflectivity was found to show no change but the diffuse reflectivity was increased from 2.0% to about 22% to provide a better match with the observed data.
The overall efficiency of the electronics and PMTs was calibrated by aligning the simulated spectrum of the number of prompt PMT hits to that from the 16 N calibration data at the center of the detector.

D. Energy reconstruction
The kinetic energy of an event T e is reconstructed by comparing the observed and expected numbers of prompt hits, defined as those with time residuals within [−10, 8] ns, based on simulation using the reconstructed position and direction.Since events are reconstructed under the hypothesis of an incoming electron, the observable energy of a gamma particle will reconstruct below its true energy due to the effects of Compton scattering and Cherenkov threshold, such that a 6 MeV gamma reconstructs around 5 MeV.Based upon an early 16 N calibration run, 6.4 prompt PMT hits are expected per MeV of electron-equivalent deposited energy.

E. Energy systematic uncertainties
The relative energy scale and detector energy resolution were determined by fitting the reconstructed energy spectrum from the 16 N calibration source, as shown in Fig. 3, with the predicted energy spectrum, generated from simulation, convolved with a Gaussian distribution [24].The fit is characterized in terms of a scale, as a linear correction to the energy, and resolution, relating to the width of the spectrum.
Events were associated with detector volume cells based on their reconstructed positions.The cells were distributed across z and ρ ≡ x 2 + y 2 , with dimensions 57 cm × 200 cm, determined based on statistics.Events were also separated into seven bins based on their u • r value.In addition to 82 distinct positions from scans within the AV, the 16 N source was also deployed at 19 positions along a vertical path between the AV and PMT array.Since the mean free path of the 6.1 MeV gamma is approximately 40 cm in water, all cells within the AV contained data but, due to limitations in the deployment positions of the source, some cells outside the AV were not populated.The energy scale values of the cells, mapped in z and ρ, were fit with a continuous function to provide a smooth calibration of energy across the detector.For each z-ρ cell, fits were performed within all seven u • r bins and then averaged to provide a solid angle-weighted energy scale in the cell.
The deployment of the source was simulated at the same positions and with the same detector conditions.Half of each of the measured and simulated data sets were used to calibrate the energy scale as a function of position.The resulting calibrations were applied to the remaining halves of the 16 N data sets, which were then fit to determine the relative energy scales and resolutions.Uncertainties were volume and solid angle-weighted according to the selection criteria for different analysis re-  gions.For the nucleon decay analysis, the relative energy scale and resolution uncertainties in the signal volume were 2.0% and 1.8× T e /MeV%, respectively, where T e is the reconstructed energy.The total energy uncertainties are dominated by the variation with position and the statistical uncertainty of the calibration data set.Contributions from source-related and time-related uncertainties were taken into account.

V. BACKGROUND MODEL
Several sources of background were considered for this search, mainly from naturally occurring radioactive contamination from the 238 U and 232 Th chains in the various detector components.Other sources include interactions from solar neutrinos, atmospheric neutrinos and reactor neutrinos.The analysis cuts, shown in Table III, were chosen to limit the impact of these backgrounds in the region of interest (ROI) for the nucleon decay study.

A. Instrumental backgrounds
Backgrounds from detector instrumentation consisting of light emitted from within the PMTs, electrical breakdowns in the PMT base or connectors, and electronic pickup, were present during data taking.A suite of data reduction cuts were developed, based on those used by SNO, to remove these instrumental backgrounds.The cuts rely on low-level PMT information such as charge, hit time, and hit location.The total sacrifice within the nucleon decay ROI due to these cuts was measured to be 1.7% by applying the reduction cuts to 16 N calibration source data.The remaining instrumental contamination is evaluated using a bifurcated analysis method [19], in which two sets of data-cleaning cuts (bifurcation branches) were used simultaneously on the analysis data; the instrumental contamination is then calculated using the number of events that pass or fail one or both bifurcation branches.Using the 16 N sacrifice estimates, a correction to the contamination estimate is also applied for the estimated signal events flagged by the bifurcation branches.The number of events expected within each data set are included in Table IV.

B. Radioactive backgrounds
Three radioactive background analyses were performed in order to estimate the contribution of 214 Bi (U-chain) and 208 Tl (Th-chain) decays from the detector components and the detection medium in the nucleon decay ROI.One analysis was dedicated to the radioactivity from U and Th chains in the the water inside the AV, while two were dedicated to the radioactivity in the acrylic itself, the hold-down rope system, the water shielding, and the PMTs.Note that a new, sealed cover gas system in SNO + to suppress radon ingress, from the headspace volume above and into the water in the AV below, had not yet been brought online during the data taking period reported here.This resulted in somewhat elevated and variable levels of 214 Bi in the AV water due to the lack of a radon-free cover gas.

Internal radioactivity
214 Bi and 208 Tl decays within the AV water were distinguished by fitting to the β 14 distribution in a background analysis region defined by a radius of 4.3 m and energy above 4 MeV, to minimize the contamination from decays from the AV and water shielding.Monte Carlo simulations of 214 Bi and 208 Tl decays were used to construct PDFs of the data within the background analysis region.In each of the data sets, the rates of 214 Bi and 208 Tl were fit simultaneously to account for possible correlations, with the resulting concentrations shown in Table V.This is demonstrated for data set 3 in Fig. 4. The rate is then extrapolated from the background analysis region into the nucleon decay ROI based on the relative expected event rate between the two regions from MC simulations.

External radioactivity
Two independent analyses were performed to measure the radioactivity from the AV, hold-down ropes and water shielding.A two-dimensional likelihood fit measured the rates of U and Th components based on β 14 as a function of R 3 -where R is the reconstructed radius of the event as measured from the center of the AV-within a background analysis region chosen to preferentially contain events from the AV, ropes and water shielding.Results from the fit analysis are shown in Table V. 214  208 Tl were fit simultaneously, taking into account correlations between the two.The second analysis counted the events within four analysis regions, defined by R 3 and u•r to contain events from the water within the AV, AV and ropes, water shielding and PMTs respectively.Monte Carlo simulations were used to translate the number of events observed in each into total rates and to extrapolate to the expected number of events in the ROI, shown in Table IV.The analysis took into account the correlation of events between the different regions.Bias testing of the fits ensured that asymmetries in the spatial distribution of events were properly handled.The analyses agree with each other within uncertainties.

PMT backgrounds
Events from radioactive decays in PMTs have a very low probability to reconstruct within the nucleon decay ROI, but occur at a very high rate, making it difficult to estimate their contribution to the total background rate using only Monte Carlo simulations.A data-driven method was used instead to put a constraint on the rate of PMT background events in the ROI.The 208 Tl decay produces a beta particle alongside a gamma cascade that can occasionally be detected in the PMT itself.The vertex reconstruction of PMT events is dominated by the Compton scatter of the gamma while the signal from the electron will appear early and concentrated in a single PMT.Events with early hits and high charge were tagged as PMT events, with simulation studies showing a tag- 86.9 ± 1.1   The values for the AV water come from the internal background analysis while the water shielding, AV and ropes numbers come from the external fit analysis, with the AV and ropes scaled together.For the external fit analysis, the statistical uncertainty is obtained from the fit, while the systematic uncertainty includes the contamination from PMT events.
ging efficiency of close to 30 %, while the misidentification rate is 1.4 %.The number of tagged events was used to estimate the number of expected events originating from the PMTs.

(alpha, n) interactions
Another potential source of background is the set of (alpha, n) reactions on 13 C atoms within the AV itself, induced by alpha particles from the decay of 210 Po, a daughter of 222 Rn.In about 10% of the cases, a high energy gamma or electron-positron pair is produced together with the neutron, which can act as a background for the nucleon decay search.Using predictions based on the leaching coefficient of 210 Po in water, temperature and the surface contamination on the AV [25], less than 600 (alpha, n) decays were expected during the period of data taking.Monte Carlo simulations of these events predict less than 1 event reconstructing in the nucleon decay ROI during the whole period under analysis.

C. Neutrino induced backgrounds
A dominant background for the nucleon decay search is the elastic scattering of 8 B solar neutrino interactions.Such events were largely excluded by a cut on cos θ , the reconstructed direction relative to the direction to the Sun.Monte Carlo simulations of 8 B solar neutrinos were constrained based on recent measurements from Super-Kamiokande [26], with oscillations applied using the BS2005-OP solar model [27].
Antineutrinos produced by nuclear reactors also contribute to the background.The expected number of events is estimated based on Monte Carlo simulations using the world reactor power data [28] with oscillations applied based on a global best fit [29].
Atmospheric neutrino interactions can create backgrounds for the nucleon decay search, particularly through neutral-current interactions with the oxygen nuclei.The interactions can liberate a nucleon from the 16 O atom, leaving either 15 N* or 15 O*, identical to the nucleon decay signal.However, a large fraction of these events can be tagged by looking for neutron followers appearing after the initial event.In order to estimate the size of this background, GENIE [30,31] was used to simulate high-energy atmospheric neutrino interactions.The GENIE Monte Carlo was verified with two studies.One study counted nucleon decay-like events with neutron followers to probe the size of the neutral-current background, and a second compared the energy, time and multiplicity of Michel electron followers directly to data.Both searches found good agreement between GE-NIE and data and the difference between the two is used as part of the atmospheric background uncertainty estimate.
Due to the location of SNO + at a depth of 2092 m, equivalent to close to 6000 m of water, the rate of cosmicray muons entering the detector is approximately three per hour.Spallation products from these cosmic-ray muons are removed by cutting all events within 20 s of a muon event, as was used during SNO.This was shown to reduce the remaining number of events from spallation products to less than one event during the data taking period [32].

VI. ANALYSIS METHODS
A blind analysis was carried out, removing events with the number of PMT hits approximately corresponding to between 5 and 15 MeV from the data set.After the analyses were finalized, blindness constraints were lifted and the whole of the data was made available for analysis.
The expected signal was simulated using Monte Carlo techniques to develop analyses that maximize the signal acceptance while minimizing the effect of backgrounds.The observables for each event used in the search were the reconstructed kinetic energy T e , the cube of the position radius R 3 , normalized by the cube of the radius of the acrylic vessel R 3  AV , position on the z-axis of the detector z, and direction relative to the solar direction cos θ , as well as the event classifiers β 14 , ITR, and u • r.
Two analysis methods, a spectral analysis and, as a cross-check, a counting analysis, were used to set a limit at 90% C.I. (credibility interval) on the number of signal decays (with a prior uniform in the decay rate), S 90% .This is then converted into a lifetime on the invisible nucleon decay modes using where N is the number of targets.For the neutron and proton modes, this is defined as the number of neutrons (or protons) in 16 O atoms inside a 6 m radius in the AV, 2.41 × 10 32 .The difference between this 6 m radius and the fiducial radius for a particular data set is accounted for in the selection efficiency of that mode.For the dinucleon modes, N is defined as the number of 16 O atoms within the same volume, 3.02 × 10 31 .To calculate the limit on the number of decays from the limit on the observed signal, an acceptance efficiency is calculated for each mode as the product of the theoretical branching ratios [14,15] and the selection efficiency from cuts on the observables with the total shown in Table VI.

A. Spectral analysis
A spectral analysis was performed, fitting for the signal in the measured distribution of the observables η = [T e , R 3 , cos θ , β 14 , u • r], within the limits defined in Table III but with energy considered over the range 5-9 MeV for all data sets.The backgrounds in the fit included solar neutrinos, reactor neutrinos and atmospheric neutrinos as well as radioactivity from U and Th chain decays in the water, AV, ropes, water shielding and PMTs.Probability distributions for the signal and backgrounds, P s (η i ) and P b (η i ), were generated using Monte Carlo simulations with constraints on the radioactive backgrounds provided by the likelihood-fit external analysis.
To allow for the multiple data sets, the analysis simultaneously maximized the sum of the log likelihoods of each individual data set k, as described by: where n obs is the number of observed events in each data set, s is the signal decay rate, s,k is the acceptance efficiency of the signal in data set k, β b is the rate of background component b whose expectation βb is constrained by σ b , b is the acceptance efficiency for the background b and t k is the live time of data set k. Fits were bias tested using a sampling of fake data sets based on Monte Carlo simulations.To find S 90% , a profile likelihood [33] distribution is calculated by taking the value of the maximum likelihood for a given value of s.The upper limit at 90% C.I. is then found by integrating along this distribution.

B. Counting analysis
A counting experiment with a set of rigid cuts, shown in Table III, is also used, where the number of background events is calculated directly from the background analyses and is shown in Table IV.Due to changes in the level of backgrounds, candidate events were selected using different cuts during different periods of running.The signal acceptance within each data set is shown in Table VI.Using a Bayesian method [34], an upper limit on the number of signal decays that could have occurred is found by numerically solving where S 90% is the upper limit on the number of signal decays at 90% credibility level and, for each data set k, b k is the number of expected background events, combined from internal and external radioactivity, solar, reactor, atmospheric and instrumental backgrounds, n k is the number of observed events while k and t k are the signal efficiency after cuts and the live time of the data set.A is a normalization factor such that the integral tends to 1 as S 90% tends to infinity.

VII. RESULTS
The results of the spectral analysis for the neutron decay mode are shown in Fig. 5 with the fitted energy spectrum of the neutron decay signal at its maximum likelihood value plotted alongside the fitted backgrounds and data.Figure 6 shows the normalized and cumulative likelihood distribution for the neutron mode.The resulting limits on each mode of invisible nucleon decay are shown in Table VII alongside the existing limits.A breakdown of systematic uncertainties is given in Table VIII.
For the counting analysis, Table IX shows the observed events compared to the predictions for each data set.The results of the counting analysis are shown alongside those of the spectral analysis in Table VII

VIII. CONCLUSION
The results shown by the spectral and counting analyses in Table VII are in good agreement.In the case of the dineutron decay mode, the spectral analysis performed significantly better due to the difference in the spectral shape of the dineutron signal, which is not taken into account within the counting analysis.
The limit set in this work on the lifetime of the proton decay mode of 3.6×10 29 y is an improvement on the existing limit from SNO, however the neutron mode limit of 2.5×10 29 y is weaker than the current limit from Kam-LAND.
For the dinucleon modes, the nn limit of 1.3×10 28 y does not improve upon the existing limit set by Kam-LAND, but the pn and pp mode limits of 2.6×10 28 y and 4.7×10 28 y improve upon the existing limits by close to three orders of magnitude.For each of the decay modes, the best fit is given as well as the difference between the best fit and the shift-and-refit value for each source of uncertainty.The total is the sum in quadrature of each of the separate systematic uncertainties assuming no correlation between the components.The statistical uncertainty, found by integrating the likelihood function prior to convolution with the systematic uncertainties using a Feldman-Cousins confidence interval to 1σ, is shown for comparison.Finally, the convolved likelihood function is integrated to 90% to get the final limit shown on the last line.
Data Observed Expected set events events 1 1 1.17 FIG.1.Spectrum of de-excitation gamma rays emitted after the invisible decay of16 O for single nucleon decay[14] (left) and dinucleon decay[15] (right) as a function of gamma ray energy Eγ between 5 and 9 MeV.

FIG. 2 .
FIG. 2. Comparison of the x-component of reconstructed position between16 N data acquired at a central source position and Monte Carlo simulation.

FIG. 3 .
FIG. 3. Comparison of reconstructed energy between16 N data acquired at a central source position and Monte Carlo simulation.

TABLE I .
Live time within each data set, in days.

TABLE II .
Systematic uncertainties in the reconstructed position and direction of events based on analysis of16N data.

TABLE III .
Summary of cuts that define the ROI used in the counting analysis for each data set, determined from Monte Carlo simulations based on expected background levels.Data set 2 is split into different cuts for the top and bottom halves of the detector as necessitated by backgrounds.Cuts of −0.12 < β14 < 0.95 and ITR > 0.55 were used in every data set.No cuts were used on u•r.The spectral analysis shares these cuts except for considering a broader energy range of 5-9 MeV.
Bi and   FIG. 4. Data and Monte Carlo distributions of events in β14 for 214 Bi and 208 Tl decays within the AV water analysis region for data set 3. The normalizations of the Monte Carlo spectra were fit to the observed data.

TABLE V
. Measured U and Th contamination in the AV water, water shielding, AV and ropes in g/gH 2 O , g/gAV and g/grope, respectively.

TABLE VI .
The acceptance efficiency of the signal after applying counting analysis cuts (for both single nucleon and dinucleon decay modes) of each data set.The single nucleon decay efficiencies are given per nucleon while the dinucleon modes are given for the nucleus as a whole.Combined systematic uncertainties are given, statistical uncertainties were negligible.
. A breakdown of the systematic uncertainties is shown in TableX.Fitted energy spectrum across all data sets for neutron decay and backgrounds.Backgrounds were fit for each data set and the spectrum shown is the sum over all of the time-bins for the individual components.The contributions of Bi and Tl were fit independently for internal backgrounds but were merged in this plot.The signal is shown at its maximum likelihood, not at 90% C.I.The errors around the full fit include the MC statistical uncertainties summed in quadrature with the individual systematic uncertainties.These errors are bin-by-bin correlated and are dominated by the energy resolution and the energy scale systematic uncertainties.The likelihood ratio (the maximum likelihood as a function of the signal divided by the likelihood at the best fit value) for the neutron decay mode of the spectral analysis versus the limit of the number of signal decays per day.The associated cumulative distribution indicates the point at which the limit on the signal at 90% C.I. is drawn.

TABLE IX .
The observed and predicted number of events passing the counting analysis cuts for each data set.The uncertainty on the total number is found by treating each data set as independent and combining in quadrature.The first uncertainty is the statistical uncertainty, the second is the total systematic uncertainty.The statistical uncertainty is almost completely dominated by the background from the PMTs.Systematicn (events/day) p (events/day) pp (events/day) pn (events/day) nn (events/day)

TABLE X .
Systematic uncertainties for the counting analysis on the upper limit (at 90% C.I.) on the signal decays per day, shown as the difference to the unshifted value.