Searches for Baryon Number Violation in Neutrino Experiments: A White Paper

Baryon number conservation is not guaranteed by any fundamental symmetry within the Standard Model, and therefore has been a subject of experimental and theoretical scrutiny for decades. So far, no evidence for baryon number violation has been observed. Large underground detectors have long been used for both neutrino detection and searches for baryon number violating processes. The next generation of large neutrino detectors will seek to improve upon the limits set by past and current experiments and will cover a range of lifetimes predicted by several Grand Unified Theories. In this White Paper, we summarize theoretical motivations and experimental aspects of searches for baryon number violation in neutrino experiments.


Contents
1 Introduction and Executive Summary

Theory Summary
The stability of ordinary matter has long been a subject of both theoretical and experimental interests. The electron is stable because of electric charge conservation. On the other hand, in the Standard Model (SM), the stability of the proton is guaranteed by the accidental global symmetry of baryon number at the renormalizable level. In models of quark-lepton unification, such as the Grand Unified Theories (GUTs), baryon number is necessarily violated. As a result, the proton is not stable, and decays dominantly into e + π 0 (in non-supersymmetric theories) or K +ν (in supersymmetric theories). Another compelling reason for baryon non-conservation is that understanding the origin of matter in the universe requires it as one of the three fundamental Sakharov conditions in addition to CP violation and thermal non-equilibrium. This striking prediction of GUTs on proton decay, which are otherwise inaccessible to laboratory experiments, motivated the construction of large-scale water Cherenkov detectors like IMB, Kamiokande, and its subsequent upgrade, Super-Kamiokande. Although there is no direct evidence of proton decay so far, but only stringent lower limits on the proton lifetime, it is important to continue the searches for proton (and bound neutron) decay and other baryon number violating (BNV) processes in general. A large class of GUTs predict proton lifetime to be within an order of magnitude above the current experimental limit. It is also important to keep in mind that the same experiments originally constructed to search for proton decay have now become truly multi-purpose experiments. In particular, they have played a major role in neutrino physics, starting with the serendipitous detection of SN1987A neutrinos, as well as the discovery of neutrino oscillations in atmospheric and solar neutrinos.
Therefore, the significance of current and next-generation neutrino experiments simultaneously searching for baryon number violation and studying neutrino properties cannot be overemphasized. While the main focus of the BNV experiments is on proton decay searches, there also exist other equally important baryon and/or lepton number violating processes, such as dinucleon decays and neutron-antineutron oscillations which must be studied as well along with their experimental detection prospects. Possible connections of BNV observables to other beyond the Standard Model physics, such as neutrino mass, baryogenesis, dark matter, flavor physics, and gravitational waves are also being explored. The recent lattice developments for the relevant nucleon and nuclear matrix elements of effective BNV operators are also crucial for reducing the theoretical uncertainties in the BNV predictions.

Experimental Summary
Experiments have long sought evidence of the decay of the proton as proof of physics beyond the SM (BSM). The lower limit on the proton's lifetime is currently of order 10 34 years. Experimental searches that seek to probe beyond this limit therefore need a huge source of protons and years of exposure. The large mass and long operation times of detectors used for observation of neutrino oscillations make them well-suited for searches for baryon number violation, including nucleon decay and neutronantineutron oscillations. The Super-Kamiokande neutrino experiment, which uses a water Cherenkov detector with a fiducial mass of 22.5 ktons and began operation in 1996, has published leading limits on 30 baryon number violating processes. Next-generation neutrino detectors, such as DUNE (40 kton fiducial mass liquid argon TPC), Hyper-Kamiokande (190 kton fiducial mass water Cherenkov), and JUNO (20 kton fiducial mass liquid scintillator), all include searches for baryon number violation as a major component of their physics programs and hope to improve upon the limits set by Super-Kamiokande, if not observe baryon number violation for the first time.
Detector mass is a crucial characteristic in next-generation baryon number violation searches. For small detectors, the exposure required to improve upon limits already set by Super-Kamiokande can exceed the likely lifetime of the experiment. Clearly Hyper-Kamiokande has the advantage in this respect. That being said, detector technology is also extremely important; DUNE's excellent imaging capabilities and JUNO's superb timing resolution offer advantages in some channels over Hyper-Kamiokande's larger mass. NOvA, a currently-running neutrino experiment with a 14 kton segmented liquid scintillator detector, is developing a search for neutron-antineutron oscillations that could potentially have sensitivity comparable to current limits. THEIA is a proposed water-based liquid scintillator detector that would combine the advantages of the large mass of a water Cherenkov detector with the good resolution of a liquid scintillator detector. With this worldwide program, should a baryon number violating signal be observed by any of the detectors in the next generation, confirmation from other detectors using different technologies would provide powerful evidence of physics beyond the Standard Model.
In addition to detector mass and technology, simulation and analysis techniques can also affect the potential of these searches. As with neutrino interactions, the experimental community has come to understand how important nuclear effects are in predicting the characteristics of final-state particles. Final-state interactions in the nucleus alter the multiplicity and momenta of final-state particles. Uncertainties in modeling final state interactions therefore introduce uncertainties into the signal efficiency estimates and lifetime limits. Furthermore, analysis techniques are continually improving. For example, Super-Kamiokande made improvements to the search for proton decay via p → e + π 0 by reducing backgrounds via neutron tagging. Potential improvements to searches in a liquid argon TPC could come from tagging of nuclear de-excitations.
The experimental neutrino physics community has long been conducting searches for baryon number violation using neutrino detectors. The next generation of neutrino detectors, situated in underground laboratories with good shielding to reduce cosmic ray backgrounds, will allow the continued pursuit of this goal, with massive detectors and continually improving analysis techniques.

Historical Context
"Is ordinary matter stable?" This question has been a subject of experimental and theoretical interest over many decades. Ordinary matter made up of electrons and nucleons (protons and neutrons) is stable first because electrons are stable because of electric charge conservation; but what about nucleons? So far, experimental searches for the decay of protons and bound neutrons have not found evidence for nucleon decay and have only led to decay mode dependent upper limits on the nucleon lifetime. From a theoretical point of view, the stability of matter can be guaranteed by assigning a baryon number B = +1 to the proton (the lightest baryon), which is known as the principle of conservation of baryon number formulated by Weyl in 1929 [1] (see also Refs. [2][3][4]). Baryon number conservation, however, is not guaranteed by any fundamental symmetry within the SM. Instead, baryon number is an accidental classical (global) symmetry of the SM Lagrangian, which is violated by small amounts via non-perturbative effects (namely, SU(2) L sphalerons) [5][6][7][8]. These are however negligibly small at temperatures low compared with the electroweak scale, but are important in the early universe. be considered an appealing extensions of the SM. In order for SUSY to provide a solution to the fine-tuning problem with the Higgs mass, it was necessary that the SUSY breaking scale should not be much higher than the electroweak symmetry breaking scale of 250 GeV. Thus, it was widely expected that supersymmetric partners of SM particles would be observed at the LHC. However, no evidence of these superpartners (in particular, the squarks and gluinos, which interact strongly) has been seen at the LHC running at 14 TeV [13]. This subsection summarizes the state of theoretical knowledge of both non-SUSY and SUSY GUTs. In our discussion, we mostly focus on minimal models; however, nucleon decay predictions of a wide range of models are summarized in Table 2. For reviews on this subject, see Refs. [20,[48][49][50].
• SU (5) GUTs: The simplest GUT model is the GG model [16] with the gauge group G GUT = SU (5). In this model, the GUT group is spontaneously broken to the G SM in a single step by a Higgs field in the adjoint (24) representation. Finally, SM symmetry is broken down to SU (3) c ×U (1) em when a scalar field in the fundamental (5) representation acquires vacuum expectation value. The latter field contains the SM Higgs doublet that remains light, whereas its color-triplet partner needs to reside roughly above 10 11 GeV since it leads to proton decay predominantly through p →νK + channel [51]. The SM fermions of each generation are contained in a 5-and a 10-dimensional representations. Notably, the 5 contains the lepton doublet and the d c quark field, while the 10 contains the quark doublet, and the u c as well as the e c fields. The gauge multiplet belonging to the adjoint representation contains twelve SM and twelve BSM vector bosons (labeled by X and Y with electric charges 4/3 and 1/3, respectively). These new gauge boson interactions involve both quarks and leptons; consequently, the new interactions violate baryon number B, and lead to proton decay via dimension-6 operators of the form such as u c γ µ Qe c γ µ Q etc [48,[52][53][54]. An example diagram is presented in Fig. 1 (diagram on the left).
Non-observation of proton decay requires these gauge bosons to be superheavy, and this bound can be computed easily by approximating the left diagram in Fig. 1 by a four-fermion interaction; by doing so, one obtains, where g GUT is the unified gauge coupling and m p and M X are the proton and superheavy gauge boson masses, respectively. Then the current experimental bound of τ p (p → e + π 0 ) > 2.4 × 10 34 yrs from the Super-Kamiokande Collaboration [55] typically implies M X ∼ M GUT 5 × 10 15 GeV. In addition to p → e + π 0 , current (as well as future) experimental bounds (sensitivities) on other important proton decay modes are summarized in Table 1.
In fact, the minimal SU (5) GUT (GG model) in combination with imprecise unification of gauge couplings predicts the proton lifetime of order τ p ≈ 10 28 − 10 32 yrs and was already ruled out by early proton decay experiments. Moreover, there are additional flaws of the GG model: (i) gauge couplings do not unify, (ii) it predicts an incorrect fermion mass relation, m e /m µ = m d /m s , (iii) there is a doublet-triplet splitting problem [62,63] (a generic problem in most of the GUT models), i.e., the fine tuning needed to render the electroweak-doublet Higgs in the 5-dimensional SU (5) Higgs light while keeping the color-triplet components at the GUT scale, and (iii) neutrinos remain massless, as in the SM. Nevertheless, many proposals were made in the literature towards building a viable GUT by extending the GG model; for an incomplete list of realistic models that also incorporate non-zero neutrino masses, see, e.g., Refs. [20,[64][65][66][67][68][69][70][71][72][73][74][75][76][77][78][79][80][81][82][83]. One of the most minimal renormalizable extensions of the GG model was advanced recently in Ref. [82] and employs fields residing in the first five lowest dimensional representations of the SU (5) group. In this model, a vectorlike fermion in the 15-  Figure 1: Left-diagram: Dominant proton decay mode p → e + π 0 in non-SUSY GUTs (here (X µ ) a i = (X a µ , Y a µ ) T is the iso-doublet gauge field). Right-diagram: Dominant proton decay mode p → νK + in SUSY GUTs (example diagram with Higgsino dressing, see text for details). The blob here represents the dimension-5 operator induced by colored Higgs exchange. dimensional representation is assisted with a 35-dimensional scalar multiplet to serve three purposes: (a) the wrong fermion mass relations are corrected, (b) unification is achieved at a high enough scale to be consistent with the current experimental bounds, and (c) neutrinos receive non-zero masses via 1-loop quantum corrections. In addition to having the least number of parameters in the Yukawa sector, this model inextricably connects the origin of neutrino masses with the experimentally observed difference between the charged lepton and down-type quark masses. Due to its minimality, this model turns out to be very predictive. The outcomes are as follows: (a) neutrinos are Majorana particles and predicted to have normal mass ordering, (b) the lightest neutrino is massless, (c) there are four light scalar relics at or below a 100 TeV mass scale, and (d) the proton lifetime is expected to be within τ p (p → e + π 0 ) ≈ 10 34 − 10 36 yrs with an upper bound of τ p 2.3 × 10 36 yrs.
One can also attempt to correct the wrong fermion mass relations by adding a 45-dimensional Higgs field [84] or by considering non-renormalizable operators [85]. However, these modifications reduce the predictivity of the model. Higher-dimensional operators, however, are not enough to provide sufficient contributions to neutrino masses. In building realistic extensions along this line, the options include the following (i) a 15-dimensional scalar [64] or (ii) a 24-dimensional fermionic [65] representation. The former (latter) option generates correct neutrino mass via a type-II [86][87][88][89] [91][92][93][94][95]) seesaw mechanism. The analyses of both these theories typically suggest that proton lifetime has an upper bound of τ p (p → e + π 0 ) 10 36 yrs [68]. On the other hand, predictivity can be achieved by employing the single operator dominance [96,97], where a single higher-dimensional operator dominates each Yukawa entry. Group-theoretical factors from GUT symmetry breaking can lead to predictions for the ratios of quark and lepton masses at the unification scale [96,97] that can be utilized to correct the bad mass relations between the down-type quarks and the charged-leptons. For a recent analysis utilizing this concept in the context of non-SUSY SU (5) GUT that predicts nucleon decay rates from specific quark-lepton Yukawa ratios at the GUT scale, see, e.g., Ref. [98].
• SUSY SU (5) GUTs: If supersymmetry is realized in nature and the SUSY breaking scale is not too far above the electroweak scale, then gauge coupling unification takes place close to M GUT = 2 × 10 16 GeV, and therefore the gauge-boson-induced proton lifetime is predicted to be τ p ≈ 10 35 yrs [99], which is above the current experimental lower limit. However, SUSY GUTs predict new proton decay contributions arising from dangerous dimension-5 [52,100] operators of the type QQQL, among which the dominant decay mode is p → νK + [101,102] generated by the exchange of the colored Higgs multiplet (of mass M Hc ); these operators are proportional to 1/M Hc . Therefore the colored Higgs fields must be superheavy. By dressing these dimension-5 operators by gauginos or Higgsinos, dimension-6 operators are generated [103][104][105] that cause the proton to decay as shown in Fig. 1 (diagram on the right). Typically the Wino exchange dominates; however, for larger tanβ, the Higgsino exchange (the one shown in Fig. 1) also becomes important. Also, the Higgsino exchange dominates the Wino exchange if the Higgsino mass is much larger than the Wino mass-such a situation is realized in, e.g., mini-split SUSY. This is because the loop diagram like the right figure in Fig. 1 is accompanied with a chirality flip and thus is proportional to the mass of the exchanged particle. Interestingly, assuming SUSY particles of masses of order of the electroweak scale, the lifetime of the proton is typically found to be τ p (p → νK + ) 10 30 yrs [105,106]. On the other hand, the Super-Kamiokande experiment gives a stringent limit on this channel: τ p (p → νK + ) 6.6 × 10 33 yrs [58]. This apparent contradiction rules out minimal SUSY SU (5) GUT [42,43] with low-scale SUSY (typically for sfermion masses TeV).
The solution to this problem resides in the Higgs boson mass which is currently accurately measured at 125.35 ± 0.15 GeV. This is close to the upper limit of ∼ 131 GeV predicted in supergravity unified models [108,109] and implies that the weak SUSY scale is high lying in the several TeV region which in turn implies that the sfermion masses could be high, even as large as 10 TeV and above. Such a scale could still be natural on the hyperbolic branch of radiative breaking of the electroweak symmetry [110]. Thus the proton lifetime exhibits a sensitive dependence on the Higgs boson mass [111]. Using the precise value of the Higgs boson mass one finds for the proton decay mode νK + a lifetime range of 3 × 10 34 − 2 × 10 35 yrs for mSUGRA and 3 × 10 34 − 10 36 yrs for SUGRA models with nonuniversal soft breaking (NUSUGRA) consistent with analyses that a heavy SUSY spectrum might revive the minimal SU (5) GUT [99,105,[112][113][114][115][116]. Further, SUSY CP phases enter in proton decay and can lead to cancellations in the decay rate enhancing the lifetime [117]. Similar cancellations due to specific flavor mixings are possible and can suppress proton decay rate, see, e.g., Refs. [118,119]. The cancellation mechanism is quite general and can apply to wide class of unified models including SO(10) [120]. Additionally, higher-dimensional operators can affect proton decay textures which are in general different from the textures for the fermion masses and affect proton lifetime [121,122] as well as can easily increase the triplet mass and thus increase the proton lifetime significantly, see, e.g., Refs. [119,123,124]. Recently, proton decay in minimal SUSY SU (5) has been revisited in Ref. [107] (see also Ref. [125]) where various soft SUSY-breaking conditions are imposed and uncertainties asso- Figure 2: Theoretical predictions of proton lifetime for representative GUT models are presented (for the underlying assumptions, see text). (c here represents the coefficient of a Planck suppressed dimension-5 operator, for details, see Ref. [107].) Current Super-K data rule out the gray shaded regions. Future projections/sensitives from JUNO, DUNE, THEIA, and Hyper-K are also specified in the diagram (see Section 3 for details). ciated with numerous phenomenological inputs in the calculation of the proton lifetime are analyzed. For example, in constraint MSSM (cMSSM) case, assuming sparticle masses O(10) TeV, the proton lifetime is found to be τ p (p → νK + ) (2 − 6) × 10 34 yrs which can be tested in the near future.
We emphasize that unusual decay modes, such as p → µ + π 0 and p → µ + K 0 , can also be within the reach of the future experiments even in the minimal SUSY SU(5), if there exists flavor violation in sfermion mass matrices [113]. These decay modes can also be enhanced in flipped SU(5) GUTs with R symmetry, as recently discussed in Refs. [126,127].
• SO(10) GUTs: GUTs based on SO(10) gauge symmetry [18,19] are especially attractive since all SM fermions of each family are unified into a single irreducible 16-dimensional representation. Additionally, 16 contains a right-handed neutrino which naturally leads to non-zero neutrino mass via a type-I seesaw mechanism. Furthermore, SO(10) is free of gauge anomalies, whereas, in contrast, in SU (5), the gauge anomaly due to the 5 R cancels the anomaly due to the 10 L of fermions (separately for each generation). The unification of all fermions of each generation into a single multiplet suggests that SO(10) may serve as a fertile ground for addressing the flavor puzzle. In fact, it turns out that SO(10) GUTs are very predictive and have only a limited number of parameters to fit charged fermion as well as neutrino masses and mixings, which have been extensively analyzed in the literature . As shown in Ref. [144], the most economical Yukawa sector with only SO(10) gauge symmetry consists of a real 10 H , a real 120 H , and a complex 126 H Higgs fields. Another widely studied class of models [128], with minimal Yukawa sector, utilizes a complex 10 H and a complex 126 H Figure 3: Various symmetry breaking chains [171] in SO(10) GUTs (reproduced from Ref. [172]). The left (right) diagram shows breaking chains along the Pati-Salam (Georgi-Glashow) route.
fields, which relies on the additional Peccei-Quinn symmetry exterior to the SO(10) gauge symmetry. Moreover, from decays of the heavy right-handed neutrinos, the matter-antimatter asymmetry of the Universe [8,150,151] can also be generated [152] in SO(10) GUTs.
Since the SO(10) group has rank 5, whereas the SM has rank 4, there are various possible pathways for the GUT symmetry to break down to the SM as demonstrated in Fig. 3. Depending on which Higgs multiplet is employed to break the GUT symmetry, at the classical level, there are typically four possibilities with minimal Higgs content for the intermediate gauge symmetry [20,153] [154] is obtained if GUT symmetry is broken by a 54 [155,156] (210 [157][158][159]) Higgs representation, whereas left-right symmetry with (without) D-parity is achieved if a 210 [157][158][159] (45 + 54 [160][161][162]) Higgs is used. Depending on the intermediate gauge symmetry, these models with predominant decay mode p → e + π 0 predict a lifetime in a wide range that varies in between 10 28 and 10 40 yrs [153]. For a recent study along similar lines using a non-minimal Higgs sector in generic SO(10) models see Ref. [163,164]. In these works, amounts of threshold corrections required to be consistent with present bounds on proton decay for various intermediate breaking chains are quantified. In Ref. [165], a minimal renormalizable model with a symmetry breaking sector consisting of Higgs fields in the 10 + 54 + 126 is analyzed and shown to have an upper limit on the lifetime τ p (p → e + π 0 ) 5 × 10 35 yrs. A variant of this framework with the adjoint 45 scalar instead of 54 has been studied thoroughly in [166,167] due to its very special quantum structure [168][169][170] and particular robustness with respect to the leading Planck-scale effects. are found to be M GU T 10 16 GeV and τ p 10 35 yrs, respectively.
For another model that uses a combination of type I and type II seesaw for neutrino masses and leads to a prediction of proton lifetime in the accessible range is in Ref. [174]. The model solves also the baryogenesis and DM problems of the SM.
• SUSY SO(10) GUTs: In SUSY SO(10) GUTs, if the intermediate symmetry is broken by a 126 Higgs that breaks B − L by two units, in addition to giving large masses to the right-handed neutrinos (that in turn generates tiny neutrino masses), R-parity of the MSSM automatically emerges from within the SO(10) symmetry. Therefore, SUSY SO(10) GUTs are highly attractive as they naturally provide a DM candidate and prohibit d = 4 baryon number violating operators without any additional symmetry assumptions. The minimal renormalizable model [128,197,198] with a symmetry breaking sector consisting of 10 + 126 + 126 + 210 was found to be very successful in fitting fermion data [130][131][132][133][134][136][137][138][139][140][141]. Soon after, it was realized [137,[199][200][201] that the intermediate symmetry breaking scale that is required to be of order 10 12 − 10 13 GeV from fits to light neutrino masses, subsequently leads to certain colored particles from various Higgs fields with masses of order 10 10 GeV that spoil the successful perturbative gauge coupling unification. Besides that, current proton decay limits completely rule out this version of the model with TeV scale SUSY spectrum.
These caveats can be resolved by extending the minimal version of the SUSY SO(10) GUT. Without introducing any new parameter in the Yukawa sector, it was shown in Ref. [146] that a minimal realistic extension is to add a 54 multiplet to the symmetry breaking sector. This setup allows the intermediate breaking scale of order 10 12 − 10 13 GeV to fit neutrino data, yet preserves perturbative gauge coupling unification. The viability of this model requires a mini-split SUSY spectrum with sfermion masses of order 100 TeV using the current experimental lower bound τ p (p → νK + ) > 6.6 × 10 33 yrs. Even though squarks and sleptons are far beyond the reach of LHC, the model can be tested at colliders with the discovery of LSP (Wino) along with its charged partners. Improvement of the proton lifetime limits in the channel p → νK + by one order of magnitude will highly disfavor this version of the model (expected upper limit on the proton lifetime with O(100) TeV sfermions is τ p 6 × 10 34 yrs). On the other hand, implementation of Peccei-Quinn (PQ) symmetry [202] to solve the strong CP problem (for a recent review see Ref. [203]) in renormalizable SUSY SO(10) with minimal Yukawa sector allows TeV scale sfermion masses as shown in Ref. [147]. This is possible since Higgsino mediated proton decay rate is strongly suppressed by an additional factor of (M PQ /M GUT ) 2 [147,204] and the expected proton lifetime in this framework is τ p ≈ 10 33 − 10 36 yrs. This scenario is exciting since the proton decay suppression mechanism is inherently linked to the solution to the strong CP problem.
The axion solution to the strong CP problem is particularly interesting since axion can be a cold DM candidate. Models in which the axion field is embedded into a scalar representation responsible for breaking the GUT symmetry show an interesting correlation between the proton decay rate and the axion mass. Models of this type were first proposed in the context of SU (5) in Ref. [205] and SO(10) in Ref. [206]. For recent works along this line, see, e.g., Refs. [207][208][209][210].
Another option to accommodate low scale SUSY in the SO(10) context is to permit more parameters in the Yukawa sector by introducing a 120 multiplet on top of 10 + 126 [186,211]. With new parameters in the Yukawa sector, it is possible to realize some cancellations to save the theory from too rapid proton decay even with a TeV scale SUSY spectrum. As shown in Ref. [186], type-II seesaw case is highly disfavored by the current proton decay lower limit, whereas for type-I scenario, the proton lifetime for the channel p → νK + can reach up to 10 37 yrs.
SUSY SO(10) GUTs that adopt small Higgs representations, to be specific, 16 instead of 126, do not guarantee automatic R-parity, but they still provide quite simply matter-parity which keeps the lightest SUSY particle stable to serve as DM. Nevertheless, models of this class (consisting of 16 + 16, 10, and 45 Higgs fields) are particularly interesting, not only because doublet-triplet splitting occurs naturally but also an interesting correlation of the decay rates between the two major proton decay modes p → νK + and p → e 0 π + emerges [212]. This interdependence shows that the empirical lower limit on the lifetime for p → νK + provides a constrained upper limit on the lifetime for p → e + π 0 mode and vice versa. Based on this correlation, estimated upper limits for proton lifetimes have been obtained in the context of an updated version of Ref. [212] (to be published; private communications from the authors of Ref. [212]; see also Sec. 7.2 of Ref. [213]). This updated version takes two factors into account: (a) the current LHC constraints on SUSY spectrum by allowing for an inverted hierarchy of the sfermion masses with a light stop (∼ 800 GeV to 1 TeV) and lighter higgsino/bino (∼ 700-800GeV), together with heavy masses (∼ 18-20 TeV) for the sfermions of the first two generations on the one hand, and (b) the appropriate normalization factors in the definitions of certain effective couplings involving identical fields, which were missed in Ref. [212], on the other hand. With this updating, the correlation mentioned above yields conservative estimates for the upper limits of τ p (p → e + π 0 ) less than about (2 − 10) × 10 34 yes, and τ p (p → νK + ) less than about (1 − 8) × 10 34 yrs. These upper limits are within the reach of the Hyper-K and Dune detectors including their planned upgrades respectively.
Moreover, GUT models with family symmetries (for an incomplete list of references, see, e.g., Refs. [214][215][216][217][218][219][220][221][222][223][224][225][226][227]) are particularly interesting since they tend to explain the flavor puzzle. An important aspect of this class of models is that Yukawa couplings emerge dynamically from minimization of the flavon potential, thereby reducing the number of parameters considerably. For example, assuming S 4 flavor symmetry on top of SO(10) gauge symmetry, in Ref. [193], in addition to explaining the hierarchies in the charged fermion masses and mixings, neutrino observables are also predicted (such as θ 13 ∼ 9 • ). Furthermore, the proton decay mode p → νK + in this model is found to have a lifetime of ∼ 10 34 yrs, which is within reach of the upcoming experiments.
It should be noted that unified models such as SO(10) with appropriate symmetry breaking that produce the SM gauge group are also constrained by electroweak physics. It is then of interest to investigate the consistency of the unified models (such as Yukawa coupling unification) with the recent result from the Fermilab E989 experiment [12] on the muon anomalous magnetic moment. The Fermilab experiment has measured a µ = (g µ − 2)/2 with a significantly greater accuracy than the previous Brookhaven [228] experiment and the combined Fermilab and Brookhaven experimental data gives a 4.2σ deviation [229] from the SM. An investigation of the Yukawa coupling unification for the third generation in a class of SUSY SO(10) unified models [230,231] shows that the SO(10) model is fully consistent with the Fermilab result.
Finally we mention two new classes of SO(10) models not discussed thus far. As noted above, SO(10) models require several Higgs representations to break the GUT symmetry to the SM gauge group. Thus a 16 or a 126 of Higgs field is needed to change rank, and one of 45, 54 or a 210 is needed to break the symmetry down further, and to achieve electroweak symmetry breaking and to generate quark and lepton masses an additional 10 dimensional representation is also needed. Recently a class of SO(10) models have been proposed consisting of 144 + 144 of Higgs fields [232,233] which can break the symmetry down to SU (3) c × U (1) em . Proton decay analysis in this model exhibits a cancellation mechanism consistent with the current experimental constraints and the possibility of observation in future experiment [120]. Another class of SO(10) models relates to the doublet-triplet splitting problem which as noted earlier is quite generic in grand unified models. One way to resolve it is the so called missing VEV mechanism [234,235] (see also Refs. [189,212,236]) where the vacuum expectation value of a 45 Higgs field which breaks the SO(10) symmetry lies in the (B − L)-preserving direction, and generates masses for the Higgs triplets but not for the Higgs doublets from a 10-plet. This mechanism works in SO(10) and has no analog in SU (5). A second approach is the missing partner mechanism [230,237,238]. For SO(10) it leads to a variety of models discussed in [230]. B − L = −2 violating dimension-7 and dimension-9 operators have been computed in this class of models [239]. Thus previous analyses using a bottom up effective low energy operator approach can now be extended to a top down one [239] for GUT scale baryogenesis and for B − L = −2 proton decay such p → νπ + , n → e − π + , e − K + and n −n oscillations.
In short, well-motivated non-SUSY and SUSY GUTs generically predict rates of BNV processes that can be probed by next-generation experiments if not already ruled out by the current experimental data. A sketch of theoretical predictions for selected models and experimental reach of upcoming detectors are depicted in Fig. 2 for non-SUSY and SUSY GUTs. Furthermore, nucleon decay predictions for a wide range of models are summarized in Table 2. For details on theoretical assumptions associated with each model's predictions, the readers are referred to the original literature.
As a cautionary remark, it is worth noting that none of the predictions in Fig. 2 or in Table 2 is actually sharp; one typically encounters ranges stretching over several orders of magnitude. This has to do with a number of theoretical uncertainties affecting the precision of the calculations at various levels of significance. These can be loosely divided into three main classes corresponding to different ways the quantitative estimates based on diagrams in Fig. 1 are influenced: (i) uncertainties in the determination of the masses of the relevant leptoquark fields (i.e., the GUT scale), (ii) uncertainties in the couplings (gauge, Yukawa) governing the GUT-scale dynamics and (iii) uncertainties in the relevant hadronic or nuclear matrix elements. As for the first class, the most prominent of these effects are the uncertainties related to the generally unknown shape of the GUT-scale spectrum of the models at stakes, to the proximity of the GUT and the Planck scales enhancing the uncontrolled corrections from higher-dimensional operators (for instance those due to gravity effects) [240][241][242] and to the limited precision attainable in the perturbative accounts (see e.g. Ref. [170]), all inflicting uncertainties in the GUT-scale matching conditions. The second class corresponds to the limits in our understanding of the flavor structure of the relevant B-and L-violating charged currents (cf. [243,244]) which, indeed, gets only partially reflected in the flavor observables accessible at low energies (the quark and lepton masses and mixings). It is also worth mentioning that unknown CP phases in the GUT Yukawa couplings [245] can change the proton decay rate significantly (see, e.g., Ref. [246] and Ref. [107]). The third class concerns the general difficulty with ab initio QCD calculations in the strongly coupled low-energy regime. While the modern lattice techniques have recently brought enormous progress in point (iii), cf. Sect. 2.9, (i) and (ii) still represent a formidable challenge to any precision calculations with rather limited prospects for any significant near-future improvement.
Thus far, we have discussed gauge-mediated proton decay, which dominates if the mass scale of vector-bosons and scalar-bosons are of a similar order. This happens for the latter contribution because the first-generation Yukawa couplings suppress the scalar-mediated proton decay. Naively, color-triplet scalar contributions become significant only if M S O(10 −4 )M V . It should be emphasized that the scalar-mediated contributions depend highly on the Yukawa sector of the theory and are model-dependent to a large extent. For minimal SO(10) GUTs, using typical Yukawa coupling structures, color-triplet masses need to be heavier than about 10 10 − 10 11 GeV [247]. In the context of a minimal model based on 10 H and 126 H with Peccei-Quinn symmetry, when scalar-mediated contributions dominate, the proton decay spectrum is found to be quite different from the one typically anticipated [247]. For example, (i) proton dominantly decays into νK + or µ + K 0 for lighter (3, 1, 1/3) or (3, 1, −1/3), respectively, and (ii) BR(p → µ + π 0 ) BR(p → e + π 0 ) is expected. Both these features are distinct from gauge-mediated proton decays; hence, if scalar-mediated contributions dominate, proton decay can provide very useful insight into the Yukawa structure of the theory.
• GUTs in extra spatial dimensions: Extra spatial dimensions provide a useful avenue for GUT model building in which some of the usual problems of four-dimensional GUTs can be addressed. For example, some versions of orbifold GUTs [184,[248][249][250] use compactification to break the unified gauge symmetry avoiding the doublet-triplet splitting problem. Moreover, localizing fermion fields appropriately in the bulk of extra dimension, a natural understanding of hierarchical fermion mass spectrum can be achieved in certain versions of GUTs [251][252][253]. In 6D GUTs, the origin of multiple families of matter fields can be realized by quantization of flux in torus constructed from the compact extra two dimensions [254][255][256]. In these classes of GUTs, the Yukawa couplings become calculable parameters and the quarks and lepton spectrum can be obtained from a very small number of parameters. In the extra-dimensional GUT models, proton decay can have distinct features. For example, in a class of 5D models, the proton decay mediated by vector bosons has an additional suppression due to the wave-function profiles of first-generation fermions and the underlying vector bosons in the bulk [185,249,257,258]. In SUSY orbifold GUTs, the dimension-5 operators are often suppressed due to U (1) R symmetry [185,258]. In 6D models with flux compactification, the proton decay also receives a non-trivial contribution from the higher Kaluza-Klein modes of the vector bosons [259]. This, along with the special flavour structure of these theories, implies proton dominantly decaying into µ + π 0 which is otherwise suppressed in the 4D GUT models. The absence of indication of low energy supersymmetry in the experimental searches so far has inspired investigations of the GUTs in the regime of split [138] or high scale SUSY [260]. Gauge coupling unification in this class of models is achieved by keeping only a part of SUSY spectrum at the intermediate energy scale while the remaining MSSM fields stay close to the GUT scale. Despite having most of the spectrum at intermediate or at very high energies, this class of theories are still reasonably constrained from the perspectives of Higgs mass, stability of the electroweak vacuum, flavour transitions, dark matter and proton decay [261][262][263][264]. It has been shown that if TeV scale Higgssino forms almost all of the dark matter of the universe then it leads to particular ranges for the unified couplings and the scale of unification which in turn put an upper bound on the proton decay requiring the proton lifetime less than 7 × 10 35 years [263,264].
Finally, we point out that a public software package SusyTCProton for nucleon decay calculations in non-SUSY and SUSY GUTs is available, which includes, e.g., the full loop-dressing at the SUSY scale; for details, see Ref. [265].

Pati-Salam Partial Unification
The first step towards unification was made in Ref. [14] (see also Refs. [15,213,266]) that are based on partial unification with non-Abelian gauge group G 422 = SU (4) C × SU (2) L × SU (2) R . By proposing the G 422 symmetry, this work [14] introduced many new concepts into the literature, which include: (i) Quark-lepton unification; (ii) The first realistic model of matter and its three interactions which quantized electric charge, and thus (as was realized later) the first quantum-theoretically consistent model of magnetic monopoles; (iii) The first left-right symmetric gauge structure providing a compelling reason for the existence of the right-handed neutrino and the right-handed gauge-boson W R ; and (iv) The first model that questioned baryon number conservation and proton stability in the context of higher unification; while this was in the context of integer-charge quarks, nevertheless it served to remove partly the major conceptual barrier against BNV violation and proton instability that existed in the community in the early 1970s.
In this PS model, unlike SO(10), there is no gauge-mediated proton decay. In fact in the original PS model, proton decay appears only if the quarks are chosen to have integer charges. This is why PS symmetry can be realized at rather low energy scales. In minimal models, the PS breaking scale can be as low as v R 10 3 TeV [267,268], which is coming from experimental limits on the branching ratios for K 0 L → µ ± e ∓ mediated by the new gauge bosons X µ carrying 4/3 charge under B − L. However, in the extended models, where additional fermionic degrees of freedom are introduced, the PS symmetry can even break down close to the TeV scale that has the potential to be directly probed at colliders; for a recent study, see, e.g., Ref. [269]. In light of recent anomalies [270][271][272][273] in the beauty-meson decays, low scale PS models have gained a lot of attention. Within this setup, the vector leptoquark X µ is an attractive candidate to explain some of the flavor anomalies [274][275][276][277][278][279][280][281][282]. Figure 4: Nucleon decay modes in the SU (4) C × SU (2) L × SU (2) R model with inclusion of color sextet Higgs fields ∆ R a la Ref. [283]. Here the scalar sub-multiplets ∆ qq R , Σ 8 , and {Σ 3 , ∆ q R } are color sextet, di-quark, and leptoquarks, respectively; moreover, ∆ 0 R is SM singlet that acquires a VEV and breaks B − L by two units. The Feynman diagram on the left (right) leads to nucleon decay of the type 3q → qq (3q → ).
The difference between the left and the right Feynman diagrams in Fig. 4 is: the right-diagram is obtained from the left-diagram by replacing the di-quark Σ 8 by a leptoquark Σ 3 . A priori, the amplitudes for these alternative decay modes could compete with each other. It is interesting to note that, in a certain region of the parameter space of the theory, the three lepton final states can dominate over single lepton final states; moreover, all these BNV processes can be within the observable range. Stringent limits on three-body decays p → e + (µ + )νν have been set by the Super-Kamiokande experiment, with τ > 2 × 10 32 yrs [294]. Unlike two-body nucleon decay channels, phase space alone provides only a crude approximation for the resulting particle spectra from these processes. For the trilepton modes, dynamics and phase space can be approximately accounted for using a general effective Fermi theory formalism of electroweak muon decay µ → e + νν [295]. Searches for nucleon decay in other channels have also been conducted and typically have present lower bounds in the range of 10 31 − 10 33 yrs [13].

Neutron-antineutron Oscillation
Oscillations of electrically neutral particles are well-known phenomena; for example, neutrino oscillations and neutral meson ( established. Therefore, one would naturally expect to observe neutron-antineutron (n−n) oscillations. However, the conservation of baryon number forbids the transition from a neutron (with one unit of baryon number) to an antineutron (with negative one unit of baryon number). As already pointed out, in the SM there is a global symmetry which forbids baryon number violation. However, once we go beyond the SM, there is no obvious reason to expect baryon number conservation to hold. Indeed, baryon number violation is one of the necessary conditions noted by Sakharov for the generation of a baryon asymmetry in the universe [296]. It was observed by Kuzmin that n-n oscillations might provide a source of baryon number violation connected with this generation of baryon asymmetry in the universe [297]. Proton decay is mediated by four-fermion operators, which have coefficients of the dimensional form 1/M 2 , whereas n-n oscillations are mediated by six-quark operators, which have coefficients of the dimensional form 1/M 5 . Consequently, if there were only one high scale M responsible for baryon number violation, one might expect that n-n oscillations would be suppressed relative to proton (and bound neutron) decay. However, there are theories beyond the SM where BNV nucleon decays are absent or highly suppressed and n-n oscillations are the main manifestation of baryon number violation [283,298]. Observation of n − n transition would show that baryon number is violated by two units |∆B| = 2, which can be naturally expected within GUTs and other extensions of the SM. However, the selection rule |∆B| = 2 is again very different from p → e + π 0 that follows ∆B = −1 (and ∆L = −1, hence ∆(B − L) = 0). Nucleon decay with a selection rule ∆B = −1 are induced by dimension-6 (or dimension-5) operators and, therefore, would correspond to the existence of new physics at an energy scale of about 10 15 GeV, whereas n − n transitions with |∆B| = 2 are induced by dimension-9 operators, and therefore, would imply new physics around the 100 TeV scale. Some early studies of n-n oscillations include [283,[299][300][301][302][303]. Some recent reviews include [304,305]. Figure 5: Feynman diagram for n − n oscillation in the model of Ref. [283]. Here, sub-multiplets ∆ R,ab are color sextets and ∆ 0 R acquires a VEV that breaks SU (4) C × SU (2) L × SU (2) R symmetry to the SM.
Although, in general, a theory may violate L without violating B, in most of the well-motivated extensions of the SM, these two are connected via fundamental symmetries. As an example, within leftright symmetric models, U (1) B−L is promoted to a gauge symmetry [306]. As long as this symmetry is exact, ∆B = ∆L is realized. This is highly motivating since the existence of |∆L| = 2 generates the Majorana neutrino mass term, which is then connected to the existence of |∆B| = 2 operators. Therefore, n − n oscillation may be directly linked to neutrinoless double beta decay processes; see Section 2.8.1.
It was shown in Ref. [283] that if the Higgs multiplet that breaks the SU (4) C × SU (2) L × SU (2) R group to the SM group is ∆ R (10, 1, 3) (instead of the Higgs set chosen in the original PS model), it generates Majorana neutrino mass via |∆L| = 2 processes while at the same time yielding n − n transition via |∆B| = 2 processes at an observable rate, keeping the proton stable. A term of the form ∆ 4 R in the scalar potential gives rise to a cubic term among three color sextet scalars once the B − L violating vacuum expectation value of the neutral component of ∆ R is inserted, which causes n − n oscillations [283,304] as shown in Fig. 5. The n − n oscillations arising in this class of models with TeV scale color sextet scalars have an upper limit of τ nn 10 10−11 sec [307,308], within the reach of the next-generation experiments; see Fig. 19. This upper bound is a consequence of the model requirements to satisfy the neutrino oscillation data via type-II seesaw, observed baryon asymmetry via the post-sphaleron baryogenesis mechanism [309], as well as the low-energy flavor changing neutral current (FCNC) constraints. Besides, TeV scale color sextets have rich phenomenology and can be uniquely probed at the hadron colliders [310][311][312][313][314][315][316]. There exist other simplified models of n − n [317][318][319], which, however, do not have an upper limit on τ nn , although they can still be tested using the collider probes of colored particles.
An interesting class of BSM theories hypothesizes that our usual four spacetime dimensions are embedded in a space with higher (spatial) dimensions, such that SM fermions have strongly localized wave functions in the extra dimensions [320]. It has also been shown that in these extra-dimensional models, neutron-anti-neutron oscillations can occur at an observable rate while baryon-number-violating nucleon decays are suppressed far below experimental limits keeping the proton stable [298,321,322]. The BNV nucleon decays are suppressed by separating quark wave function centers sufficiently far from lepton wave function centers in the extra dimensions, but this does not suppress n-n oscillations since the operators that mediate these oscillations only involve quarks. These models can also produce a seesaw mechanism with naturally light neutrino masses and viable DM candidates [323,324]. Another study of baryon-number violation without proton decay is [325].
An experimental search for n-n oscillations was carried out using neutrons from a reactor at the Institut Laue-Langevin (ILL) and obtained a lower limit τ nn > 0.86 × 10 8 sec (90 % CL) [326]. Neutron-antineutron oscillations also lead to instability of matter with associated ∆B = −2 dinucleon decays yielding mainly multi-pion final states. These have been searched for in deep underground nucleon decay detectors, most recently in Super-Kamiokande [327,328], obtaining a lower limit on a lifetime characterizing such matter instability of τ > 3.6 × 10 32 yr, which, when converted to a corresponding free neutron oscillation time, is τ nn > 4.7 × 10 8 sec. This difference is due to the suppression of n − n oscillations in matter as a result of nuclear potential differences. The oscillation time is matter (τ m ) is related to free neutron oscillation time (τ n−n ) as τ m = R τ 2 n−n , with the conversion factor R for Oxygen evaluated to be about 0.5 × 10 23 sec −1 [329]. Upcoming experiments at the ESS (European Spallation Source) and DUNE (Deep Underground Neutrino Experiment) plan to improve these bounds to τ nn 10 9−10 sec [330][331][332]. These oscillation times can be translated into new physics scales of roughly τ nn Λ 6 QCD 1/5 ∼ O(100 − 1000) TeV, which are well above energies directly accessible at colliders. For recent reviews on this subject, see Ref. [304,305]. Existing upper bounds on the rates for the dinucleon decays nn → π 0 π 0 , nn → π + π − , and np → π + π 0 imply upper bounds on the rates for dinucleon decays to dileptons nn → e + e − , nn → µ + µ − , nn → νν, and np → + ν , where = e, µ. These have been calculated in [298,333].

Other B − L violating Processes
Nonzero neutrino masses, if these are Majorana fermions, are evidence for new physics that violates lepton number by two units. GUTs, as argued earlier, will ultimately lead to a connection between lepton-number and baryon-number violating new physics. While the details depend on both the mechanism behind nonzero neutrino masses and the underlying GUT, it is often the case that the lepton number violating physics is a subset of different phenomena correlated with ∆(B − L) = 2 phenomena. Some connections have been explored at the effective operator level, in a number of papers; e.g., Ref. [303,334].
The strongest bounds on these ∆(B − L) = 2 processes come from searches for nucleon decay. These are nicely summarized in Ref. [13], and some are very strong, of order 10 32 years. It is also the case that all of the strongest bounds on ∆(B − L) = 2 neutron and proton decay come from twentieth century experiments, including Frejus [342] and IMB [343,344], keeping in mind that it is not possible to distinguish a nucleon decay into a neutrino from that into an antineutrino. The DUNE experiment, thanks to its tracking, calorimetric, and particle-ID capabilities, is well positioned to be sensitive to ∆(B − L) = 2 with lifetimes that approach 10 34 years [57].

Effective Field Theory of B Violation
The most suitable framework for studying low-energy, below electroweak-scale phenomena is the Low Energy Effective Field Theory (LEFT) [345], which differs from the SM with respect to the internal symmetry as well as the degrees of freedom. In going from the SM to the LEFT, the heavy degrees of freedom namely, the Higgs boson h, the electroweak gauge bosons W ± µ and Z µ and the top-quark t are integrated out, whereas the internal symmetry is spontaneously broken from SU LEFT operators can describe a wide variety of fermion-fermion interactions, also encompassing the BNV ones. The best examples are operators describing the proton decay into mesons and leptons or entirely into three charged leptons [346,347]. In the SM baryon and lepton number conservation emerge as accidental symmetries. This suggests that the violation of these symmetries within the LEFT is a sign of beyond the SM (BSM) physics. The backdrop for conducting analyses on and to test the veracity of BSM models is the Standard Model Effective Field Theory (SMEFT) and it also acts as the necessary bridge between the low and the high energy sectors. A complete matching between a BSM model and the LEFT would involve an intricate and meticulous multi-step procedure of integrating out heavy fields and matching parameters. But a simpler way to catch a glimpse of possible BSM origins of the rare processes encapsulated by LEFT operators is to find their embedding within SMEFT operators and then based on symmetry arguments unfurl them into Feynman diagrams consisting of heavy propagators [348].
We have highlighted two examples in Fig. 6, where the first column presents LEFT contact operators of mass dimension 9 describing proton decay processes (i) p → K + e + e − ν and (ii) p → e + e + e − The second column shows the relevant SMEFT operators of the same or higher mass dimension that describe the same contact interaction as the LEFT ones [348]. Finally, the last column reveals the BSM heavy field propagators that can mediate such processes.
The ∆B = 0 operators of the SMEFT provide a model-independent framework for exploring Bviolation, both for B-violating nucleon decay [53,290,291,349] and for n-n oscillations [283,300,302,303,350]. The flavor structure of the operators at different dimensions allows one to establish nucleon decay-mediating processes that can dominate. This has been explored in further detail in more recent works such as Refs. [281,321,[351][352][353][354]. In SMEFT, ∆B = 1 appears at dimension-6. Any flavor ∆B = 1 term leads to nucleon decay, including particles heavier than the proton such as charm or the tau lepton that can induce proton decay through off-shell contributions (see e.g. Ref. [355]). These contributions are strongly constrained by two-body nucleon decays such as p → e + π 0 . At dimension > 6, non-trivial lepton number ∆L = 0 allows to enforce dominance of some operators (e.g. Ref. [291]). ∆B and ∆L are connected through the dimensionful operators (e.g. [303,351]). Higher dimensional operators can also often lead to multi-body channels, such as n → K + µ + e − e − at dimension-9, or multi-nucleon decays with ∆B > 1 [352]. One can use limits on rates for p → + M and n →νM , where M denotes a pseudoscalar or vector meson to obtain indirect limits on rates for p → + + − , n →ν + − , p → + νν, and n →ννν [356]. Fig. 7 displays characteristic examples of processes with distinct ∆B and ∆L structures. Figure 7: Process examples with baryon and lepton number violation by ∆B and ∆L units, respectively. "Instanton" refers to processes that break the same quantum numbers as non-perturbative electroweak instantons. 0ν2(4)β refers to neutrinoless double (quadruple) beta decay. The minimal mass dimension d of the underlying effective operator is shown. Operators also carry flavor. Reproduced from Ref. [352].

Discrete Symmetries and Supersymmetry
B-violation processes appear in many extensions of the SM, a notable example being supersymmetric (SUSY) theories. Already in the MSSM realization nucleon decay-mediating dimension-4 operators QLd c and u c d c d c appear, where Q, L are left-chiral quark and lepton doublets and u c , d c are the utype, d-type superfields, respectively. To forbid rapid proton decay through these interactions, models often impose a Z 2 symmetry called R-parity (matter parity) (e.g. Ref. [41]). However, at dimension-5, one encounters nucleon decay-mediating QQQL, which can be forbidden by "proton hexality" Z 6 symmetry that contains R-parity as a subgroup [357]. Since all global symmetries are expected to be violated at some level [358], it is appealing to consider discrete gauge symmetries. Such symmetries can appear as remnants of spontaneously broken local U (1) symmetries [359]. Thus, care must be taken to ensure anomaly cancellation. Favorable discrete symmetries that allow for rich phenomenology, resolve theoretical puzzles (e.g. µ-problem) and forbid dangerously rapid nucleon decay have been identified, such as those of Refs. [357,[360][361][362][363]. Discrete (gauge) symmetries can also be considered in the context of GUT models.
Unlike the SM, the baryon and lepton numbers are no longer accidental symmetries of the classical Lagrangian in the MSSM. In the most general supersymmetric theory, gauge-invariant terms that violate baryon and lepton numbers are allowed, which are In the low-energy MSSM model, one forbids these terms by imposing a new discrete symmetry, known as "R-parity" [364]. One can allow some of the R-parity-violating (RPV) terms while still avoiding excessively rapid proton decay, and these could lead to |∆L| = 2 decays such as K + → π − µ + µ + and K + → π − µ + e + ; thus upper bounds on the branching ratios for these decays can be used to obtain upper bounds on the coefficients of certain RPV terms [337].
Here we consider the possibility of only baryon number violation. Therefore, the only non-zero term in the above superpotential we allow is d i d j u k . The associated coupling λ ijk is antisymmetric in the first two flavor indices, leading to in total nine couplings λ dsu , λ dbu , and so on. Due to this antisymmetric nature, coupling to three quarks of the first generation is absent; this is why n − n oscillation is highly suppressed. However, there exist stringent bounds on these couplings coming from dinucleon decay [365][366][367][368][369][370][371]. For example, let us consider the coupling λ dsu , which violates baryon number by one unit and strangeness by one unit; however, it conserves B − S. Hence it is easy to understand that the best bound comes from dinucleon decay into two mesons of identical strangeness via dimension-9 operators. Using the current limits of dinucleon decay to two kaons, pp → K + K + for instance, Super-K provides a lower limit of > 1.7 × 10 32 yrs [372] (for pion final states, limits of similar order exist from Super-K [373]), and one obtains the following constraint on the R-parity-violating coupling, where m S represents a common mass scale of SUSY particles and Λ is the hadronic scale. Dinucleon decays also appear in non-SUSY models, e.g., dinucleon decay to two kaons and dinucleon decay to leptons occur in a class of models studied in Ref. [325]. As suggested in Ref. [375], baryon number violation in two processes with at least one obeying the selection rule ∆(B − L) = ±2 can infer Majorana nature of neutrinos. As an example, consider two BNV processes: (i) p → e + π 0 (|∆B| = 1) and (ii) n − n oscillation (|∆B| = 2). The former and the latter processes correspond to ∆(B−L) = 0 and ∆(B−L) = −2, respectively. Suppose both these BNV processes are observed in the experiments. In that case, neutrinoless double beta decay is guaranteed to exist just by combining these two vertices as depicted in Fig. 8. This would subsequently confirm the Majorana character of neutrinos. One can get to the same conclusion if instead (i) p → e + π 0 (∆(B − L) = 0) and (ii) n → e − π + (∆(B − L) = −2) processes are considered [375]. Since nucleon decay in well-motivated GUT models is within reach of ongoing and upcoming experiments, one can hope to resolve this outstanding puzzle (i.e., Majorana or Dirac) by observing BNV processes, even if 0νββ is not directly observed by then.
By following the arguments given above, the presence of Majorana mass for neutrinos can be directly inferred via the weak instanton effects. Non-perturbative instanton/sphaleron configurations of weak interactions that lead to solutions to an effective operator involving twelve SM doublet fermions can be written as: [uddudd]·[uude]·[νν]. Of this operator, the first, second, and third pieces correspond to n − n oscillation, proton decay p → e + π 0 , and Majorana neutrino mass, respectively. Therefore, if p → e + π 0 and n − n oscillation are observed, the sphaleron configuration would imply that neutrinos have Majorana masses; this non-trivial connection known as the "B − L triangle" [375] is demonstrated in Fig. 9.

Leptogenesis
In the SM, baryon number and lepton number are accidental symmetries of the classical Lagrangian. However, B and L are exactly conserved only in perturbation theory and are not respected by nonperturbative effects. Sphaleron processes which are effective non-perturbative interactions constructed out of twelve left-handed SM fermions are the key to leptogenesis [152]. These interactions change baryon number B and lepton number L by a multiple of three (∆B = ∆L = 3N CS , where the n − n p → e + π 0 ν = ν c sphaleron Figure 9: B − L triangle [375] demonstrating how Majorana mass for neutrinos can be confirmed if BNV processes are observed. integer N CS is the Chern-Simons number characterizing the sphaleron gauge field configuration), but on the contrary, preserve B − L. Sphaleron processes are in thermal equilibrium in a wide range of temperatures: T EW ∼ 10 2 GeV T T sph ∼ 10 12 GeV [8]; owing to these interactions, in Sakharov's conditions [296] for baryogenesis, lepton number violation can replace baryon number violation.
The seesaw mechanism is the most natural candidate to explain tiny neutrino masses within the GUT framework or BSM frameworks with a local B −L symmetry [306]. In this scenario, tiny neutrino masses arise due to the heaviness of right-handed partners. The heavy right handed neutrinos have Majorana masses breaking B − L by two units which then gives a Majorana masses of light neutrinos. Neutrino mass in type-I seesaw is given by, where, m D = v/ √ 2 Y ν with Y ν being the Dirac neutrino Yukawa coupling, and M R is the right-handed neutrino mass. For order one Yukawa coupling for the neutrinos, one then obtains: consistent with neutrino oscillation data. In the vanilla leptogenesis scenario when N 1 , the lightest of the right-handed neutrinos, decays into lepton-Higgs pairs, a lepton asymmetry in these CP-violating out-of-equilibrium decays are generated, which is then partially converted to baryon asymmetry by the sphaleron processes and the right amount [387] of matter-antimatter asymmetry of the Universe can be reproduced. Leptogenesis and its connection to absolute neutrino mass and neutrino mixing parameters in the context of SO(10) GUTs are studied in Refs. [140,[388][389][390].
The CP -asymmetry from the right-handed neutrino decays can be estimated to be ∼ m ν M R /(16πv 2 ). A sufficiently large CP -asymmetry 10 −7 to generate adequate baryon asymmetry imposes the socalled Davidson-Ibarra bound M R 10 9 GeV [391]. However, if a pair of right-handed neutrinos are quasi-degenerate in mass, the CP asymmetry can be resonantly enhanced [392][393][394] (see also Refs. [393,[395][396][397]). Subsequently, the lower bound on M R comes only from the requirement that sufficient baryon asymmetry is induced at T > T sph . In this resonant leptogenesis scenario, right-handed neutrinos can have masses as low as sub-GeV to TeV range and low-scale seesaw models of this type have the virtue of being directly probed in experiments [384]. For implications of low-scale resonant leptogenesis from GUT, see e.g., Ref. [398].

Baryogenesis
The two leading proton decay modes, p → e + π 0 and p → νK + as extensively discussed above, both conserve B − L symmetry (i.e., ∆(B − L) = 0). For this reason, the minimal SU (5) model does not have any way to explain the origin of matter even though it has both baryon number violation as well as CP violation. Sphaleron processes would wash out any baryon asymmetry generated in such models. However, going beyond these d = 6 operators, the next-to-leading BNV operators correspond to d = 7, an explicit computation of which for a class of SO(10) models is given in Ref. [239], which break this symmetry by ∆(B − L) = −2 and generate nucleon decay modes such as p → νπ + , n → e − π + , n → e − K + . As mentioned above, spontaneous breaking of B − L in SO(10) GUTs is directly related to the neutrino mass generation. Intriguingly, in a class of SO(10) GUTs (applicable to both non-SUSY and SUSY), these d = 7 nucleon decay modes for which the lifetime can be in the experimentally accessible range are inherently linked [353,354] to the matter-antimatter asymmetry in the Universe. Due to their B − L breaking nature, electroweak sphaleron processes would not wash out such a GUT-scale-induced baryon asymmetry; hence asymmetry generated would survive down to low temperatures. The vacuum expectation value (VEV) of 126 Higgs breaks the B − L generator of the SO(10) and provides a large Majorana mass for right-handed neutrinos as well as generates trilinear scalar couplings that induce d = 7 nucleon decay operators; an example diagram is presented in Fig. 10. In this figure, a trilinear coupling ωHρ * is induced via VEV insertion of a ∆(1, 1, 0) scalar that breaks the B − L generator. Here, the quantum numbers of these fields are: ω(3, 1, −1/3), H(1, 2, 1/2), ρ(3, 2, 1/6) and this term originates from a quartic coupling in the scalar potential of the form 126 4 , which is gauge invariant. From the left (right) vertex in the Feynman diagram Fig. 10, it can be understood that ω (ρ) has B − L = −2/3 (B − L = +4/3); this clearly dictates that the associated trilinear coupling leads to the decay ω → H * ρ which violates B − L by −2. Thus, when tree as well as loop diagrams for the decay ω → H * ρ are combined, a net B − L asymmetry is generated. This GUT-scale-induced baryon asymmetry in ω → H * ρ decay that has deep correlations with neutrino mass generation and nucleon decay is shown [353,354] to correctly reproduce the observed baryon-to-entropy ratio in the Universe.

Dark Matter
The natural appearance of DM in SUSY GUTs has already been mentioned. Here we briefly discuss some DM candidates in non-SUSY GUTs. There are a large number of DM candidates in the literature, and a comprehensive review of these models is beyond the scope of this white paper. It is remarkable that GUTs based on SO(10) gauge symmetry automatically contain matter parity P M = (−1) 3(B−L) , as a discrete subgroup that can naturally stabilize the DM without the need for imposing ad hoc symmetries by hand. To be specific, stability of the DM is guaranteed by the discrete subgroup Z (2.7) The smallest representation in SO(10) that can realize this symmetry breaking pattern (i.e., with a remnant Z 2 symmetry) is 126, which can also be used to generate Majorana masses for right-handed neutrinos as discussed above. Under the matter parity, only the 16-and 144-dimensional representations are odd, whereas the rest of them are even (here, we restrict ourselves to representations not bigger than 210). Given the fact that the SM Higgs doublet residing in a 10-dimensional representation is even while a 16dimensional representation containing fermions of each generation is odd, the lightest component of an additional even fermion or an odd scalar multiplet will be stable due to the residual Z 2 symmetry [172,[403][404][405][406][407][408]. However, this stability will be lost if a Higgs in the 16-dimensional or 144-dimensional representation acquires a VEV. This leads to the list of possible DM candidates in SO(10) GUTs given in Table 3. DM candidates have also been presented in the context of extra-dimensional models (e.g., [323,324].

Flavor Violation
It is widely known that low-energy SUSY models with arbitrary mixings in the soft breaking parameters would lead to unacceptably large flavor-violating effects (e.g., [409,410]). It is essential to assume flavor universality in the mechanism that breaks SUSY to be consistent with experimental observations. Starting with flavor-universal SUSY breaking boundary conditions at a high scale, running effects of the renormalization group equations (RGEs) can still induce sizable flavor mixings in soft breaking parameters that lead to extensive flavor violations at low energies. Many constraining flavor-violating processes arise in the leptonic sector and are inherently connected to neutrino parameters. In SO(10) GUTs, Dirac neutrino Yukawa couplings Y ν induce observable lepton flavor violations (LFVs) that put strong constraints on SUSY breaking parameters [411,412]. Present experimental bounds on some of the most important LFV processes are [13]: In the RG evolution of the left-handed slepton soft-masses (m 2 L ) ij , flavor violating off-diagonal entries are induced that are proportional to Y ν Y † ν : where M R k are the heavy right-handed neutrino masses. Assuming scalar mass universality and gaugino mass unification as is usually done in cMSSM, there exist only a few parameters in the Here m 0 and m 1/2 are the universal soft SUSYbreaking (SSB) scalar and gaugino masses, respectively; A 0 is the universal tri-linear coupling and µ the Higgs mass term, which is determined from electroweak symmetry breaking condition. The fit to fermion spectrum in a specific GUT fully determines tan β as well as Y ν and the right-handed neutrino masses M R k . Consequently, one can compute the rates for LFV processes as functions of the remaining SUSY breaking parameters: {m 0 , m 1/2 , A 0 } (in this way, the gaugino, slepton, and squark masses are expressed in terms of {m 0 , m 1/2 , A 0 }). Such correlations within a minimal SO(10) GUTs are presented in Fig. 11; for details see Ref. [147]. Quark flavor violating processes such as K 0 − K 0 mixing, b → sγ, etc. also put strong constraints on SUSY breaking parameters. For a systematic study of constraints on Yukawa couplings due to GUT relations and finding correlations between constraints on flavor violation between the lepton and quark sectors, see Ref. [412].

Gravitational Wave Production
After the great success of the direct observation of gravitational waves (GWs) by the LIGO collaboration [414], GW detection has been considered to be a powerful probe of new physics that complements experiments in particle physics. BNV effects are usually associated with new symmetries (e.g., those preserving B −L) at high scales. Spontaneous breaking of these symmetries may generate gravitational radiation in the early Universe, which appears as cosmic GW background today. The GW exploration will be very helpful to understand these processes and any new physics behind them. Here we discuss potential sources of GWs generated from BNV-related physics: GWs via cosmic strings generated from intermediate symmetry breaking in GUTs and those via first-order phase transition.
The production of cosmic strings is another important prediction for most GUTs in addition to proton decay. GUTs beyond SU (5) provide a series of intermediate symmetries before breaking down to the SM. The spontaneous breaking of the GUT symmetry to the SM gauge symmetry occurs in several steps, and topological defects are produced accompanied by the symmetry breaking. Those defects include monopoles, domain walls and cosmic strings [415]. The former two are problematic as they would come to dominate the energy density of the Universe. This problem is solved by including a period of inflation after their production. The production of cosmic strings is usually associated with the breaking of a U (1) symmetry, which appears as a subgroup of SO(10) or larger groups [416]. These strings, if produced after inflation, evolve to a network. The network follows a scaling solution during the Hubble expansion and the energy density does not overclose the Universe [417][418][419].
The measurement of the cosmic GW background provides a novel way to probe grand unification [437,438]. It is worth mentioning that the new physics scale probed by GW measurements is usually not the GUT scale M X , but an intermediate scale between M X and the electroweak scale [439]. This scale is not arbitrary, but correlated with M X via the unification of gauge couplings. Measuring proton decay in neutrino experiments can determine both M X and the intermediate scale as well as Gµ. For example, given the lower bound τ (p → π 0 e + ) 1.6×10 34 yrs from Super-K [440], restrictions on Gµ in the following typical breaking chains are obtained, where D in the upper case refers the left-right matter parity symmetry and the minimal particle content consistent with the SM and neutrino masses have been considered [441]. Once any proton decay signal is confirmed in next-generation neutrino experiments, GW measurement provides a good opportunity to study the details of GUT breaking. Note that the the spontaneous breaking of the lowest intermediate gauge symmetries for all chains in Eq. 2.12 also leads to the U B−L breaking. It provides a window to connect GW signals with the Majorana nature of neutrinos and to test leptogenesis via seesaw mechanism [442].
For the symmetry breaking scale below 10 8 GeV, a strong first-order phase transition might provide a compelling source of observable GW signals. During the phase transition, bubbles of the broken phase nucleate and expand in the Universe. The collision of the bubbles and the resulting motion of the ambient cosmic fluid provides a source of a stochastic background of GWs that can be observable at GW observatories [443].
The GW spectrum is in general given by the sum of three contributions: bubble wall collisions, sound waves and hydrodynamic turbulence (see Ref. [443,444] for recent reviews). Despite the energy budget, the spectrum has a generic feature: a peak in the middle and polynomial suppressions in both low and high frequencies. It is distinguishable with that from cosmic strings as there is no plateau in the high frequency band. The frequency at the peak depends on the phase transition temperature linearly. In addition to the temperature, two more parameters are crucial to determine the GW spectrum: α the ratio of the released latent heat during the phase transition to the total energy density and β/H the ratio of the inverse duration of the phase transition to the Hubble rate. α dominates the GW strength. Enhancing α by one order of magnitude can enhance the GW strength by several orders of magnitude. β/H influences both the strength and frequency of GW. Faster the phase transition proceeding, higher frequency of the GW signal is expected. Theoretical efforts have been made to enhance α → O(0.1) and reduce β/H → O(10 2 -10 3 ) such that part of the parameter space of Bor L-violating new physics can be reached by the exploration of next-generation GW interferometers (see, e.g., discussions in [445][446][447][448][449][450][451]).

Lattice Developments
Rates of proton decays and neutron oscillations depend on nucleon and nuclear matrix elements of effective baryon number-violating operators. Prior to development of lattice QCD methods, these matrix elements were initially calculated with various nucleon models. However, eliminating theory uncertainties completely requires ab initio QCD calculations on a lattice with physical quark parameters. Calculations of proton decay amplitudes have been pursued since simulating nucleons on a lattice became possible, with methodology gradually improving over the last three decades. Neutron oscillation amplitudes were studied only relatively recently, when realistic lattice QCD calculations were already feasible.
Apart from the usual systematic effects inherent in lattice QCD calculations (discretization, finite volume, etc), calculations of BNV amplitudes are complicated by the need to preserve chiral symmetry of the quarks. Since the effective proton decay and neutron oscillation operators contain chiral quark fields, using a simple fermion discretization such as Wilson action may lead to undesirable mixing of operators on a lattice. This mixing would complicate renormalization and introduce additional systematic effects. Fortunately, formidable progress has been made towards lattice calculations with chiral fermions at the physical point [452][453][454].
Amplitudes of p →¯ π, p →¯ K decays may be computed on a lattice directly as matrix elements between the proton and a meson. Such calculations have initially been performed in quenched QCD with Wilson valence quarks [455] and DWF quarks [456], and later in unitary QCD with N f = 2 + 1 flavors of DW fermions [457,458]. All these calculations have been performed with unphysical masses of light quarks corresponding to pion masses m π 300 MeV, and required some form of chiral extrapolation to the physical point. Although these amplitudes did not exhibit strong pion mass dependence, applying chiral perturbation theory (ChPT) to baryons is often questioned, and proton decay amplitudes have not been computed in ChPT beyond the leading order (LO). In the framework of the chiral-bag proton model, it has also been suggested that the proton decay matrix elements may depend dramatically on the light quark masses [459]. These deficiencies have been addressed in the recent work [460], in which the proton decay amplitudes have been calculated with 10-20% precision using chirally symmetric DW fermions at the physical point, albeit with relatively coarse lattice spacing a 0.14 fm. The physical-point results are largely in agreement with previous calculations, and no suppression at light-quark masses has been observed.
The proton-to-meson transition amplitudes can also be "indirectly" estimated from proton-tovacuum proton decay constants α, β [461] using the leading-order chiral Lagrangian. These decay constants have been computed using quenched QCD with Wilson quarks [462][463][464] and Domain Wall quarks [456], as well as in unitary QCD with dynamical DW fermions [458,460,465]. The "indirect" estimates of proton decay amplitudes are typically higher than the "direct" results, which is likely attributable to pion loop effects. A calculation with next-to-leading order ChPT is highly desirable to understand this discrepancy. The proton-to-vacuum decay constants are also highly important to analysis of p → µ − e + e + and p → e − µ + µ + proton decays.
A particular BSM theory yields predictions for the effective scale of the physics contributing to n-n oscillations and predictions for the coefficients of the various six-quark operators O i contributing to these oscillations. To calculate the expected rate of n-n oscillations in a given model, one then needs to compute the matrix elements n|O i |n . These matrix elements have dimensions of (mass) 6 , and since at the hadronic level, the only important mass scale in the problem is Λ QCD 250 MeV, it follows that the matrix elements are ∼ Λ 6 QCD . Early calculations of the n-n matrix elements were performed using the MIT bag model [302,303]. Recently, these matrix elements have been calculated using lattice QCD [466,467]. This calculation has been performed on a lattice with chirally symmetric action at the physical point, and the operators have been renormalized using two-loop perturbative anomalous dimension [468]. As examples, for one operator, denoted Q 5 in [466], the lattice QCD calculation (normalized at 2 GeV) yielded a value in agreement with the MIT (fit A) value, while for other operators, the lattice QCD calculation yielded values larger in magnitude by up to an order of magnitude. Since the general amplitude for n-n oscillations has contributions from several different operators and since the matrix elements of these operators have different signs, there can be destructive interference, so that one cannot make a statement about the overall size of the n-n transition amplitude from a knowledge of individual matrix elements of operators without knowing the details of a given UV model predicting which operators occur and the values of the coefficients of these operators. Nevertheless, the lattice QCD calculations of these matrix elements are valuable inputs for the analysis of the predictions of a given BSM model for n-n oscillations.
Finally, lattice calculations of proton decay and neutron oscillation amplitudes inside nuclei may soon become feasible [469,470]. Although lattice simulation of nuclei requires substantially more resources due to the sign problem, development of exascale computing resources over the next decade is expected to make ab initio calculations of nuclear matrix elements possible.

Experimental Overview
This section summarizes results and sensitivities for baryon number violation searches for currently running, planned, and proposed neutrino detectors, as well as some avenues for potential improvement in nuclear modeling and detection techniques for future detectors. A comparison of current limits and future sensitivities for important modes are shown in Table 1.

Super-K
The Super-Kamiokande neutrino detection experiment represents the 2nd generation of large water Cherenkov neutrino and nucleon decay experiments, following Kamiokande and IMB which ran in the 1980s and 1990s. Super-K (also SK) is currently running and has been in operation since 1996. The experiment addresses major neutrino topics such as atmospheric neutrino oscillation, solar neutrino oscillation, gravitational collapse supernova bursts in the Milky Way galaxy, searching for the diffuse background of supernova neutrinos, and indirectly searching for dark matter via annihilation or decay to neutrinos. The SK detector is also the far detector for the T2K long-baseline neutrino oscillation experiment. The interaction target is 22.5 to 27.2 ktons of ultrapure water in a cylindrical stainless steel tank located 1 km underground in the Kamioka mine in Japan. The target mass is viewed by 11, 000 50-cm photomultiplier tubes, and is surrounded by an optically isolated active veto of roughly 2-meters of water viewed by 1, 800 PMTs. Since 2020, the detector water system has been upgraded to allow gadolinium sulfate in solution to enhance the capture and detection of neutrons [471]. This Super-Kamiokande detector naturally enables the search for a variety of baryon number violation processes in the same target mass, with only atmospheric neutrinos as a competing background. Foremost are searches for nucleon decay, often cited as one of two possible meanings embedded in the Super-Kamiokande name: Nucleon Decay Experiment.
Unfortunately for fans of new physics, the SK experiment has not detected any signs of nucleon decay. The extensive exposure, 450 kt-yr in a recent publication, sets a high bar for future experiments to overcome. For example, the lifetime lower limit for p → e + π 0 is now at 2.4 × 10 34 y based on a recent analysis that uses several improvements over older publications [55]. The improvements include: background reduction by tagging events where neutron capture on hydrogen is detected nearby the interaction, expansion of the fiducial volume from 22.5 to 27.2 ktons, and dividing the search region into a section where the nearly background-free decay of the free proton in H 2 O is accentuated. The simplicity and clarity of this search is illustrated in the three-panels of Fig. 12, showing the low background region at high invariant mass and low net momentum (note the extremely low expected background of the free proton region), and the absence of any data events in that region.
Overall, the Super-K experiment has published leading limits on 30 baryon number violating processes including: an extensive survey of antilepton plus meson final states that conserve (B − L), dinucleon decay modes that violate baryon number by two units, and three-body decay modes including those with fully leptonic final states. Special attention has been given to decay channels favored by SUSY, which are distinguished by the presence of a kaon in the final state. The most recent published limit for the key decay mode p →νK + is 5.9 × 10 33 years based on an exposure of 260 kt-yr [472], but a preliminary result of 8 × 10 33 years has been reported in conferences. In addition, Super-K holds the leading limit on the interesting ∆B = 2 process of neutron oscillation into anti-neutron within the oxygen nucleus [328]. This may be converted to an effective free n −n lifetime of 4.7 × 10 8 s, which is comparable to (slightly exceeding) the leading free neutron experiment [473]. Table 4 lists a sample of published limits.
Major improvement in sensitivity to lifetimes much beyond the limits listed here await the next generation of more massive detectors. Until then, operation of the Super-K experiment with concomitant data analysis will provide incremental progress but valuable methodological developments. Ongoing studies of nucleon decay include continued development of reconstruction techniques, refinement of the intranuclear simulations relevant to both signal and background, first studies of algorithms for neutron capture on gadolinium, and new searches for novel modes.

Neutron-Antineutron transformations in NOvA
The NOvA Far Detector is sensitive to the spontaneous conversion of neutrons to antineutrons. Since all neutrons in NOvA are bound in nuclei, the resulting antineutron would immediately annihilate on a neutron or proton, typically yielding three to six pions. The experimental signature is therefore Figure 12: The search for p → e + π 0 in 450 kton-y of Super-Kamiokande data [55]. The leftmost panels show the results of simulated proton decay, where the nuclear effects that degrade the signal for proton decay in 16 O are distinct from the low momentum events from the decay of the free hydrogen nuclei. The center panels show the predicted background due to atmospheric neutrino interactions. The rightmost panels show the SK experimental data. The upper row represents the original fiducial volume, used in the majority of Super-K publications, requiring the interaction vertex to be 200 cm from the PMT plane. The lower row represents the additional fiducial volume that is between 100 cm and 200 cm from the PMT plane. Figure 13: Simulation of neutron-antineutron oscillation in a carbon nucleus in the NOvA Far Detector, yielding π + π − π 0 . The purple lines show the true trajectories of the charged pions. a star with approximately zero net momentum and a visible energy typically between 0.8 GeV and 1.5 GeV, depending on the mode (see Fig. 13).
Neutron-antineutron oscillations are suppressed in nuclei by a factor R such that t free = t bound /R, where t free is the free neutron oscillation time. The suppression factor varies by nuclide. About half of the neutrons in NOvA are bound in carbon nuclei, with the majority of the remainder in chlorine. Although no calculations exist for these elements, one calculation of R for oxygen gives 0.517 × 10 23 s −1 [327,329]. Likely the suppression factor for carbon is somewhat lower than for the larger oxygen nucleus, so there is a possibility that the effective suppression for NOvA as a whole is lower than for a water-based detector, allowing NOvA to set stronger limits on t free , all other things being equal. NOvA must confront a significant background due to its surface location with only 3.6 meters water-equivalent overburden. Since 2018, NOvA has run a dedicated trigger searching for neutronantineutron-like events. By the end of the planned NOvA run in 2026, 112 kt-years of data will have been collected. To reduce background from cosmic rays, the trigger requires candidates to be contained, to have a similar extent in z as they do in y (unlike downward-going cosmic rays), to have only short tracks that could be produced by pions of the expected energy, to have a total physical size and number of hits consistent with a neutron-antineutron oscillation event, and finally to be symmetric between the xz and yz views.
Backgrounds to this search include atmospheric neutrinos and cosmic rays. On the surface, cosmic rays are the more challenging background. We expect that the most difficult to exclude will be from neutrons and gammas produced in the overburden or rock berms adjacent the detector. Two approaches are being pursued to study these backgrounds: first, a data-driven approach using the energy sidebands above and below the expected signal visible energy, and second, a simulation using CORSIKA [477].
The analysis to search for neutron-antineutron events in triggered data is under development. Further information available to an offline event classifier includes more sophisticated measures of event symmetry, calibrated hit energies, event time duration, and various reconstructed track variables. Using a CNN to distinguish between neutron annihilation events and cosmic rays, as is planned for DUNE [478], is also a possibility. If successful in suppressing cosmic ray backgrounds below the level of atmospheric neutrino backgrounds, NOvA will achieve a sensitivity on the free neutron oscillation time of somewhat longer than 10 8 s at 90% C.L.

MicroBooNE
The MicroBooNE detector is an 89 ton active Liquid Argon Time Projection Chamber (LArTPC) detector located on-surface and exposed to neutrinos from the Booster Neutrino Beamline (BNB) and Neutrinos at Main Injector (NuMI) beamline at Fermilab. The excellent spatial and calorimetric resolution offered by MicroBooNE's LArTPC, also shared by the future Deep Underground Neutrino Experiment (DUNE), enables precise measurements of neutrino scattering as well as beyond-Standard Model (BSM) searches, including intranuclear neutron-antineutron transitions. In a LArTPC, this process is characterized by a unique and striking star-like signature containing multiple final-state pions from the antineutron annihilation with a nearby nucleon. The pions are visible as tracks (from charged pions) or showers (from neutral pion electromagnetic decays) pointing back to the annihilation vertex, with approximately zero net momentum and total energy of about twice the nucleon mass.
Although MicroBooNE is too small in size and too overwhelmed with cosmogenic background to perform a competitive search, it aims to perform the first ever search for neutron-antineutron (n −n) transitions in a LArTPC, which will serve as a proof-of-principle demonstration of LArTPC capabilities in searching for this process. The analysis MicroBooNE has developed [479] makes use of state-of-theart reconstruction tools, including deep learning methods developed for LArTPC image data analysis, to look for neutron-antineutron transition signatures in data collected during neutrino beam-off times.
The signal events are simulated using the GENIE Neutrino MC Generator (GENIE) v.3.00.04 with hA-Local Fermi Gas (hA-LFG) as a default model, while backgrounds are expected to be contributed by cosmogenic activity during multiple 2.3 millisecond intervals of exposure (corresponding to a total of 372 seconds) and are evaluated using a data-driven approach. The signal interaction is simulated per one exposure interval of 2.3 ms, which also contains multiple reconstructed clusters of cosmogenic activity.
Pre-selection is applied to all clusters reconstructed during 2.3 ms intervals ("events") based on a Boosted Decision Tree (BDT) score, to control the background rate and remove obvious backgrounds (cosmics). The BDT is trained using only the topological information derived from the extent of the reconstructed clusters in channel and time space. As shown in Tab. 5, pre-selection reduces the cosmic background clusters such that the number of n −n clusters and cosmic clusters per 2.3 ms interval become of the same order. The pre-selected clusters are subsequently used to train a sparse Convolutional Neural Network (CNN) to differentiate signal from background clusters. For this step, 2D projections of clusters from MicroBooNE's three planes (U, V, Y) are formatted in such a way so as to retain only the important pixels, out of the full image, in the form of a collection of spacepoints containing information of wire position, time-tick and charge deposition per spacepoint. Considering statistical-only sensitivity as a figure of merit, an optimized CNN selection provides 73.6±0.034 % signal efficiency and 8.77e-3±7.4e-4 % cosmic background efficiency. As shown in Tab. 5, the final selection highly suppresses the cosmic background clusters while the number of signal clusters remains high.
The preliminary sensitivity for MicroBooNE's n −n search is calculated assuming 372 seconds of exposure, corresponding to 3.13 neutron·years, and considering statistical uncertainties only, and is shown in Tab. 6. The sensitivity is calculated using the frequentist-based Rolke method [480]. Although the n −n search in MicroBooNE is not competitive compared to existing limits from SNO or Super-Kamiokande, this analysis serves as a demonstration for the capability of future, larger and well-shielded LArTPC's, including the future DUNE, to perform such searches for baryon number violation with significantly higher exposure and thus higher sensitivity. The MicroBooNE analysis is ongoing, to incorporate effects of systematic uncertainties.

Hyper-K
Continuing the neutrino and nucleon decay physics program in Kamioka, Japan, Hyper-Kamiokande (HK) is a third-generation water Cherenkov detector currently under construction roughly 8 km south of Super-Kamiokande. The detector will use a 187 kton water target, roughly 8 times that of its predecessor, and will observe natural neutrinos from the sun, core collapse supernovae, and atmospheric neutrinos as well as a 1.3 MW neutrino beam from an upgraded J-PARC accelerator. The experiment is expected to begin operations in spring of 2027 and is expected to improve on BNV searches at SK (c.f. Section 3.1) by an order of magnitude or more using similar analysis techniques.
Importantly, Hyper-Kamiokande will be instrumented with improved 50 cm photomultiplier tubes with increased quantum and collection efficiencies, resulting in twice the photon detection efficiency of the sensors used in SK. Further, HK's sensors will have roughly half the timing resolution for single photoelectron signals. Both of these features positively impact searches for nucleon decays at HK. Notably, atmospheric neutrino backgrounds can be reduced by 30% relative to SK's achievement [55] by the augmented ability to tag the faint light from the 2.2 MeV gamma ray produced by neutron capture on hydrogen. Focusing specifically on the p →νK + mode, the improved PMT timing resolution allows Hyper-Kamiokande to better separate light from the below-Cherenkov-threshold K + 's decay products and that from photons produced by the recoiling oxygen nucleus. This improves the detection efficiency of the K + → µ + ν µ mode from 9.1% [472] to 12.7% in HK [56]. Altogether Hyper-Kamiokande's discovery potential after 10 years will exceed 3σ if the effective lifetime of the proton decay into this channel is less than 2 × 10 34 years. For the other "flagship" mode, p → e + π 0 , a lifetime less than 6 × 10 34 years, i.e. slightly more than twice the current SK limit, will lead to a 3σ detection.
Unified theories predict branching ratios for the various possible proton decay modes, indicating that observations of multiple channels and therefore comprehensive coverage of the possibilities is important for determining the symmetries of the underlying model. Table 7 lists Hyper-Kamiokande's median sensitivity to several of these modes, including some |∆(B − L)| = 2| and ∆B = 2 modes. It should be noted that for modes in which the initial state cannot be fully reconstructed from the decay products, such as p → e + νν, the search becomes background limited. In such cases the sensitivity increases with only the root of the exposure, resulting in only a factor of a few improvement in existing limits. Accordingly, the challenge for HK going forward is to leverage its improved detector to reduce backgrounds to all modes. Currently efforts to reduce backgrounds include adoption of improved particle reconstruction algorithms, improved detector calibration and light collection with additional photosensors, and introducing enhanced neutron tagging methods, either via algorithms or future gadolinium doping.

DUNE
The Deep Underground Neutrino Experiment (DUNE) promises one of the largest highly instrumented fiducial detector masses of any future large underground facility [481,482]. With 40 kt of liquid argon (LAr) some 1500 m below Lead, South Dakota to shield against cosmic ray backgrounds, DUNE's immense wire readout particle ionization charge-collection system across four separate modules forms its three-dimensional LAr time projection chambers, allowing scientists to exploit bubble-chamber quality images [483][484][485] for world-leading precision physics studies of the SM and beyond. With potentially MeV-scale precision [486], the ability to distinguish γ and e species [487], low cosmic µ backgrounds, and very low LAr ionization kinetic energy thresholds for even heavy charged species such as protons (ps) [57,482], the overall physics potential of DUNE goes far beyond its initial purpose as a ν detector built to better constrain and measure oscillation parameters such as δ CP . Indeed, the bubble-chamber-like capabilities of DUNE allow for observation of complex event topologies with potentially high multiplicities. Combined with state of the art detector reconstruction and particle identification (PID) methodologies, as well as a gargantuan number of intranuclear nucleons, there is great potential for the DUNE far detector to unlock the secrets behind rare processes.
Sensitivity to several of these processes has been studied [59] using the full DUNE simulation and reconstruction analysis chain, including the impact of nuclear modeling and FSI on a Boosted Decision Tree (BDT)-based selection algorithm. With an expected 30% signal efficiency, including anticipated reconstruction advances, and an expected background of one event per Mt·yr, a 90% confidence level (CL) lower limit on the proton lifetime in the p →νK + channel of 1.3×10 34 years can be set, assuming no signal is observed for a 400 kt-year exposure. Another potential mode for a baryon number violation search is the decay of the neutron into a charged lepton plus meson, i.e., n → e − K + . The lifetime sensitivity for a 400 kt-year exposure is estimated to be 1.1 × 10 34 years. Neutron-antineutron (n →n) oscillation is a baryon number violating process that has never been observed but is predicted by a number of BSM theories. The expected limit for the oscillation time of free neutrons for a 400 kt-year exposure is calculated to be 5.53 × 10 8 s. * Reconstruction of these events, which have final state particle kinetic energy of order 100 MeV, is a significant challenge, made more difficult by final-state interactions (FSI), which generally reduce the energy of observable particles. The dominant background for these searches is from atmospheric neutrino interactions. For example, a muon from an atmospheric ν µ n → µ − p interaction may be indistinguishable from a muon from K → µ → e decay chain from p →νK + decay, such that identification of the event relies on the kaon-proton discrimination. Neutron-antineutron oscillations can be detected via the subsequent antineutron annihilation with a neutron or a proton. The annihilation event will have a distinct, roughly spherical signature of a vertex with several emitted light hadrons (a so-called "pion star"), with total energy of twice the nucleon mass and roughly zero net momentum. Reconstructing these hadrons correctly and measuring their energies is key to identifying the signal event. As with nucleon decay, nuclear effects and FSI make the picture more complicated, as FSI can reduce the multiplicity of pions and make the pions less energetic.

JUNO
JUNO is a next-generation liquid scintillator detector under construction in southern China [60,489]. It consists of a 20 kton liquid scintillator target inside an acrylic spherical vessel surrounded by 17, 612 20-inch and 25, 600 3-inch photomultiplier tubes (PMTs). Its choice of technology, which affords it a good timing resolution and a low energy threshold, combined with its large size, unprecedented for a detector of this type, give it unique capabilities in the search for nucleon decays.
JUNO will be particularly sensitive to the p →νK + decay channel. The main reason is that, as highlighted in Sec. 3.7.2, this decay unleashes a three-fold sequence of events, each of which is detectable in a liquid scintillator detector such as JUNO: a prompt signal from the K + 's loss of kinetic energy, a short delayed signal (τ = 12.4 ns) from its decay daughter (most commonly a µ + ), and a long delayed signal from the daughter's decay (most commonly into a Michel positron with τ = 2.2 µs). This threefold coincidence provides a powerful handle to suppress backgrounds, which in JUNO are dominated by atmospheric neutrino interactions.
The sensitivity to the p →νK + decay channel has been studied using JUNO's custom Monte Carlo simulation framework [490] including a realistic detector performance. A modified version of the GENIE 3.0 [491] generator that accounts for final state interactions and residual nucleus deexcitations is used to generate p →νK + decays and atmospheric neutrino backgrounds. The short K + lifetime causes the energy depositions of the K + and its subsequent decay daughter to overlap in time. However, the 3-inch PMTs, whose transit time spread is around 1.5 ns, allow in many cases to disentangle these two signals after the time-of-flight correction has been applied, as illustrated on the left panel of Fig. 14. Fitting the two pulses simultaneously [492], as illustrated on the right panel of Fig. 14, allows to reconstruct their time separation and energy deposition. The distribution of the best-fit time separations is broader for the signal than for the atmospheric neutrino background, allowing to discriminate between these two with high efficiency [493].
There are additional handles that further enhance the signal-to-background separation. Among them is the use of a fiducial volume cut, the consideration of a muon veto to suppress cosmogenic backgrounds, and the consideration of the visible energy of the candidate signals, among others. A preliminary selection using all these criteria yields an efficiency of 31% to p →νK + decays, with only 0.3 background events in a period of 10 years. A Feldman-Cousins [494] estimation of the sensitivity to this decay mode yields a lower limit for the proton lifetime of 8.34 × 10 33 years at 90% C.L. with 10 years of data, in absence of a signal.
It is worth noting that JUNO is also expected to have good sensitivity to other decay modes, namely p → µ + µ + µ − and the n → 3ν invisible mode. JUNO will have an energy resolution of 3%  Figure 14: Left: Example of hit time distribution in JUNO for a p →νK + decay as seen with the 3-inch PMTs. Despite the proximity between the first two signals, the threefold coincidence is clearly visible. Right: Example of multi-pulse fitting in a case where the time separation between the two first pulses is 11 ns. Both images are obtained from Ref. [489].
at 1 MeV, which is unprecedented for a liquid scintillator detector. This will be an important asset in those searches where the signal is a mono-energetic energy deposition, as explained in Sec. 3.7.3. JUNO's sensitivity to these modes is still under evaluation, and the results are expected to be released soon.

THEIA
THEIA is a detector concept that utilizes advances in photon detection technology (fast timing, chromatic separation) with water-based liquid scintillator to create a highly scalable detector (up to the 100 kton scale) with better energy resolution than a pure water detector. Due to the large mass and presence of scintillation light, THEIA would have sensitivity to several modes of nucleon decay that is either comparable with or better than next-generation detectors. Additional details are available in [495].
3.7.1 p → e + π 0 and related modes Decay modes which have pions in the final state are subject to inelastic intranuclear scattering for bound protons, which causes the pion to be reabsorbed about 60% of the time. This dominates the total efficiency to see this decay mode in water as well as water-based liquid scintillator. For THEIA the dominant background comes from atmospheric neutrino interactions, and is independent of depth. Compared with other water Cherenkov detectors, THEIA would have improved neutron tagging (90% efficiency), better energy resolution, and sensitivity to below-Cherenkov threshold charged particles. These features can be all used to better reject atmospheric neutrino backgrounds.

p →νK and related modes
Kaon decay modes are less affected by intranuclear effects and produce a three-fold coincidence signal in the detector through the subsequent decays of their daughter particles [13].
• K + → π + π + π − (5.58%) • K + → π 0 e + ν e (5.07%) In a pure water Cherenkov detector the kaon is below the energy Cherenkov threshold and does not produce a detectable signal; however, the multiple pions can produce identifiable multi-ring signals, and the subsequent decays of the pions and muon create a coincidence trigger. In THEIA, the kaon would produce scintillation light as well, though much of the kaon signal would overlap with it's decay products due to it's 12-ns lifetime [60]. Since these first two signals would be difficult to distinguish, it is safe to assume that THEIA would perform no better than JUNO (pure scintillator) in terms of signal efficiency, but could be built to a much larger scale.

n → 3ν and related modes
Finally, THEIA would have sensitivity to a class of low-energy modes known as "invisible" modes, where all of the final state particles cannot produce light emissions (such as all-neutrino final states). When these occur in a multi-nucleon nucleus then the nucleus will often be left in an excited state and will release deexcitation γs and nucleons. In an 16 O nucleus this manifests as a 6.18-MeV γ with a relatively high branching ratio (44%). The primary background for this signal comes from solar neutrinos, internal radioactivity, and cosmogenic activation of oxygen to 16 N. THEIA would be the only large-scale water Cherenkov detector available to look for this decay due to its large mass and great depth. Leading limits on invisible neutron and dineutron decay are set by KamLAND [496], and on invisible proton and diproton decay by SNO+ [497].

Detection Techniques in LArTPCs
The search for baryon number violation is a prime goal of particle physics and is being carried out in large underground detectors. Most of the current lifetime limits are affected by backgrounds, predominantly from the 100 events per kiloton year which arise from atmospheric neutrino interactions. The current best limits come from water detectors, but a promising way to reduce backgrounds is with Figure 16: Search for the decay products from nuclear de-excitation of the residual radioactive daughter nuclei left after the disappearance of proton from 40 Ar nucleus. . a large liquid argon time projection chamber (LArTPC), almost all 40 Ar. Large LArTPCs provide the capability to image charged particles with mm-scale resolution and thus explore signatures from particle and nuclear physics at the same time [498]. The concurrent detection of light in a photodetector system provides an opportunity to identify signatures from nucleon decay that are not present in neutrino interactions. Analysis strategies can then be tuned to both reduce background and increase efficiency. A photon-detection system can also measure energy calorimetrically, working as a crosscheck of the energy measured by the ionizing particles in a TPC and thus improving the energy resolution when both measurements are used together. There are about one hundred nucleon decay modes accessible to large underground liquid argon TPCs, similar to the DUNE Far Detector [481], or to the considered module of opportunity liquid argon TPCs, under consideration [499], such as Q-Pix based LArTPC [500]. We see a valuable opportunity to combine signatures from particle and nuclear physics and take advantage of the fast and precise timing capabilities of a photon-detection system and the detailed pattern recognition capabilities of a LArTPC [484,501] to improve the sensitivity of future nucleon decay searches. As an example, consider the mode p →νK + [472]. For this mode, one looks for the K + signature via a short ionization track, followed primarily by K + → µ + νµ decay which makes an observable µ + track. The µ + decays to an observable positron. There are also other less likely K + decay modes. The kaon and its decay products can be reconstructed as images and the decay chain could be tested for kinematic consistency. Three potential backgrounds include 1) the large number of quasi-elastic atmospheric ν µ interactions where the recoil proton is misidentified as a K + , 2) production of K + by atmospheric neutrinos, and 3) neutrino production of K L outside the detector which enter the LArTPC and charge exchange to a K + .
There are several additional signatures of p →νK + decay to investigate, mostly requiring LArTPC mm-scale resolution and a low-energy charge detection threshold ideally combined with a fast timing from the photon system: • The 40 Ar could often become an excited state of 39 Cl, with emission of de-excitation gamma rays. This nuclear process needs additional studies since any proton inside the 40 Ar can decay. The proton that decays might be deeply bound within the nucleus and the nucleus would be, at a minimum, left with an excitation energy much higher than the neutron separation energy and maybe even other separation energies. There is a possibility of even higher excitation if part of the decay energy is transferred to the nucleus. On the LArTPC response side, the detection of MeV-scale gamma-rays was demonstrated [486] with efficiency of 50% and energy resolution of 24% at 0.5 MeV, and with an efficiency of almost 100% and energy resolution of 14% at 0.8 MeV. If the gammas from 40 Ar → 39 Cl de-excitation are measured by the photon system, in addition to the TPC charge detection, we would have a new unique time-tag of the proton decay process, see Fig 16. • The product 39 Cl decays with a 56 minute half life to 39 Ar (feasible to tag in a deep detector). Then, 39 Ar decays 93% of the time to excited states and will give you characteristic gamma rays, roughly 54% of the time a 1.27 MeV gamma, 39% of the time a 1.52 MeV gamma and 46% of the time a 250 keV gamma (this one typically comes with the 1.27 MeV gamma).
• The K + lifetime is 12 ns, and a photon system with good timing could distinguish it from its decay products, as well as improve the identification of Michel electrons, from the 2 µs decay of µ's. In recently published results of the ProtoDUNE-SP liquid argon time projection chamber performance [486], a time resolution of 14 ns was achieved in the case of a single photon-detector channel. It is expected that photon-detector time resolution will improve as a square root of the number of independent channels, driving the resolution well-below the 12 ns lifetime of K + .
• Bound neutron decay to excited states of 39 Ar might also make de-excitation gamma.
Our intent is, therefore, to further study these signatures and to examine how the signal sensitivity will be further enhanced, and how backgrounds could be rejected more efficiently with additional information from nuclear de-excitations and with precise timing from photon-detectors in TPCs.

Effect of Different Nuclear Model Configurations on Sensitivity to Intranuclear Neutron-Antineutron Transformations
Sensitivity studies for baryon number violating processes are currently relegated to simulations using generated Monte Carlo (MC) event samples. For example, using the GENIE MC event generator [502] and DUNE detector simulation and reconstruction software packages such as LArSoft [503], a first foray into sensitivity investigations for separating intranuclear n →n from (predominately neutral-current) atmospheric neutrino backgrounds in the DUNE far detector was considered in [332]. However, there remains much work to be done [504]. The dependencies of convolutional neural networks (CNN's) and other automated (machine learning) methods such as multivariate boosted decision trees' (BDTs') responses to various topological inputs is not entirely clear. The nature of these algorithms' responses to underlying choices in what are generally considered to be broadly consistent nuclear model configurations (NMCs) must be further studied to assess effective uncertainties, including those originating from disparate models of nuclear Fermi motion and final state interactions (FSIs) from (non)stochastic intranuclear cascades [488,505]. Given that the event triggers for n →n will likely utilize the signal's expected region of interest (ROI) as informed primarily by the MC simulation of the process and its separation from background via automated methods, such dependencies are highly important, particularly when physical correlations are being ignored; this is but one important component of the larger, still developing paradigm, which should include improved reconstruction and PID. In blue is the naive intranuclear radial position ofn annihilation, a probability distribution generated by a Woods-Saxon nuclear density as presented in GENIE [502]. In orange is the modern, quantum-mechanically derived intranuclear radial position of annihilation probability distribution as developed in [488]. The vertical scale is arbitrary. Top Center: The initial (anti)nucleon momentum distributions are shown using a local Fermi gas (LFG) model with an additionaln potential [488]. Top Right: The same for the GENIEv3.0.6 [502], showing an LFG model and the default nonlocal Bodek-Ritchie (BR) model. Bottom Left: A two dimensional plot of intranuclearn momentum-radius correlation is shown using an LFG and the newly-derived annihilation position distribution [488] (top left, orange). Bottom Center: The same using GENIEv3.0.6's [502] LFG model of (anti)nucleon momentum and a Woods-Saxon nuclear density (top left, blue), showing good correlation. Bottom Right: The same using GENIEv3.0.6's [502] nonlocal BR nuclear model of (anti)nucleon momentum and a Woods-Saxon nuclear density (top left, blue), showing no positional correlation, and thus over-selecting high momenta.
To explain this particular aspect in some detail, let us begin by stating the obvious: there is actual importance in the maintaining of particularly relevant physical correlations in MC simulations which act to inform the detection of unknown (rare) processes. To illustrate this, consider Figs. 17 [488]. Some correlations which have gone previously unexplored include the expected position of the intranuclear n annihilation following n →n conversion, as shown at top left. Considering the {n,n} mass splitting which suppresses n →n lessens as one hypothetically decreases the binding energy, a radial dependence is expected in the transformation probability beyond the simplistic assumption of the assumed Woods-Saxon nuclear density-indeed, such transformations are expected to occur predominately near the nuclear surface where n binding is low. Further, and as illustrated in the top central figure, given the strength of the annihilation cross section, the optical potential describing the Fermi motion of a previously convertedn is expected to be deeper than that of normal nuclear matter, imbuing the convertedn with higher available momentum. Similarly, due to n →n occurring predominately on the surface of the nucleus, one may expect a reduction in the available Fermi motion-derived momentum; thus, simulations which are inherently nonlocal in their assumptions of Fermi motion (as is the case in the GENIE [502] default nuclear model) can lead to biases in any particular ROI. The (non)locality of a given nuclear model is shown by a decrease in (anti)nucleon momentum as one moves further out toward the nuclear envelope (r → R), as shown in the bottom figures for two generators; if a hemicircular region is occupied in this parameter space, then no momentum-radius correlations are preserved. Figure 18: The final state mesonic/pionic parameter space (total momentum versus invariant mass) [327] after stochastic intranuclear transport ofn annihilation generated mesons, compared for a few NMCs, not including detector effects. The ROI is generally considered to the be "hot-spot" in the lower right hand corner, implying the expected low Fermi momentum and high invariant mass derived from the annihilation of two nucleons creating a topologically spherical π-star; differences in these may lead to different detector signal efficiencies via automated methods. Left: an LFG model with an additionaln potential and a full intranuclear cascade [488,505]. Right: GENIEv3.0.6 [502] using the default nonlocal BR relativistic Fermi gas and a full intranuclear cascade via the 2018 hN Intranuke model (hN).
The investigation in maintaining these physically relevant correlations can go even further, as one may expect fewer FSIs to be experienced by annihilation-generated mesons due to reduced views on intranuclear scattering centers near the nuclear periphery. Also, when evaluating the n →n signal's ROI (before consideration of further skewing due detector effects), choices across models of these FSIs can have a critical role in determining signal efficiencies through topological selection of high multiplicity events involving knock-out ps, which may be overproduced [488]. When considering outgoing mesons only, the nature of the ROI can be seen to remain disparate through the comparison of various NMCs across multiple (and single) generators, as shown in Figs. 18.
The importance of these correlations goes beyond the signal simulation; indeed, the same NMCs and correlations should be respected consistently across atmospheric neutrino background simulations. Iterating across the available NMCs, between and within single generators (for instance, GENIE [502]), and intermixing these together allows for the evaluation of uncertainties in a "universe" style approach, though going beyond simple knob turning; a project of this scale requires massive automated analysis techniques to be successful, especially given the nonreweightable nature of some of the more theoretically well-motivated stochastic FSI models employed. Such uncertainty evaluations will permit a better understanding of the signal to background ratio, informing the final expected sensitivities to n →n: if this ratio remains stable across NMCs, then the minutia of certain physical correlations in simulation will be shown to be unnecessary; however, the opposite is the more likely case. One study showed that a change in the NMC reduced the sensitivity by roughly 40% all else being equal [504]. Thus, this greatly encourages not only the evaluation of dependencies of automated analyses' (CNNs, BDTs) responses, but also the necessity of improved reconstruction. DUNE has had successes in machine learning being applied to PID [59], which was not included in [332]; implementing a CNN score describing the probability of particular track's PID could better discriminate signal from background considering the unique high-multiplicity π-star expected to emanate from a nucleus aftern annihilation.
From these discussions, one may hypothesize as to the regions of τ nn parameter space probed by DUNE [59], Super-Kamiokande [328], and ESS NNBAR [305]. In the context of post-sphaleron baryogenesis [308,309], one can compare potential lower limits to model predictions, as shown in Figure 19. A spread in the expected values of τ nn given various nuclear model configurations is concerning, especially as the underlying nuclear physics, though well constrained [504], is not definitively known. Figure 19: A reproduction of Fig. 9 in [308] is shown with sensitivity estimates for DUNE-assuming a 25% efficiency-and ESS NNBAR-assuming an ILL-like efficiency [326]-each in a case of zero background compared to the Super-Kamiokande I-IV limit [328].
By considering the above in greater detail into the future, sensitivities to B − L-violating ∆B = 2 processes such as n →n and related dinucleon decay modes are expected to greatly improve through enhanced physical modeling, reconstruction, and PID. Uncertainties due to specific NMC choices remain under investigation, potentially lowering the sensitivities without other improvements. These steps are required in order to pursue truly complementary physics goals beyond proton decay and neutrino oscillation studies, enabling the broader particle physics community to fully utilize future large underground detectors to better search for BSM physics, potentially informing us of our universe's origins.

Inclusive nucleon decay searches
As summarized in Sec. 2, B-violating processes are associated with a broad variety of theories and can manifest through distinct nucleon decay channels. The strongest limits can be set on particular exclusive nucleon decay processes, motivated by specific theoretical models. However, as discussed in Sec. 2.5, higher order operators leading to complicated processes can readily dominate over typical two-body nucleon decay channels. With increased complexity of interactions, performing exclusive experimental searches for all specific nucleon decay channels becomes unfeasible beyond the simplest realizations.
These considerations highlight the importance of inclusive nucleon decay searches and set limits on classes of baryon number-violating processes simultaneously. Such searches can be classified as follows [352]: • model-independent and invisible mode searches -model-independent searches probe all nucleon decays simultaneously, unconstrained by specific final states. These searches also constitute the primary approach for studying invisible channels such as n → invisible. Invisible channels such as n → 3ν could become significant in models based on extra-dimensions (e.g. [506]). Analysis of p → invisible will also serve as a test of charge conservation. Prominent searches of this class include isotope fission products as well as signatures of nuclear de-excitation, with current limits being τ 10 30 yrs. The latter method is expected to be a promising target for future neutrino experiments such as JUNO, Hyper-Kamiokande and DUNE.
• N → X + anything -allowing for model-independence while also taking advantage of precise particle identification in current and upcoming large neutrino experiments, such as Super-Kamiokande (see related searches of Ref. [476]), motivates N → X + anything searches, where X is a SM particle of unknown energy. These include: A) N → (e ± , µ ± , π ± , K ± , ρ ± , K * ,± ) + anything: primary decay charged particles are often directly visible with high efficiency in experimental searches. From electric-charge conservation any proton decay will eventually result in at least one positron, albeit potentially space-timedelayed if it is originated from the decay of a heavier charged final-state particle. B) N → (π 0 , K 0 , η, ρ, ω, K * ,0 ) + anything: searches involving neutral mesons rely on the electromagnetically interacting daughter particles for detection (e.g. π 0 → γγ). C) N → γ + anything: assuming multi-body nucleon decays, the emission of a photon will typically be suppressed compared to the same multi-body channel without a photon. D) N → ν + anything: similar to other invisible searches, can also be studied through neutrino flux measurements from accumulation of nucleon decays inside the Earth.
• ∆B > 1 processes -inclusive searches also offer prospects for ∆B > 1 processes. Note that inclusive ∆B = 1 searches also constrain ∆B > 1. The increased available energy for final state particles in multi-nucleon decay processes modifies the search compared to single-body decays.

Conclusion and Outlook
We have discussed the current status and future prospects of baryon number violation searches with large underground detectors for neutrino experiments. These detectors used for next-generation neutrino experiments have the capability to improve the existing limits on nucleon lifetime by up to two orders of magnitude. A broad class of GUT models (both non-supersymmetric and supersymmetric) can be probed. Some of them predict an upper limit on the nucleon lifetime, which might be fully within reach of the future experiments. Apart from the |∆B| = 1 nucleon decay, there exist other interesting BNV observables as well, such as neutron-antineutron oscillation (|∆B| = 2) and B → LM decays (|∆(B − L)| = 2), which also have promising prospects for underground detectors to be used in future neutrino experiments. The discovery of baryon number violation will be an unambiguous signal of new physics, and therefore, it is important to search for as many BNV channels as possible.   Selection efficiency (%) 73.6 8.77e-3 Table 5: The number of entities for each reconstruction and selection stage. The numbers are evaluated using a simulation sample that corresponds to roughly 10 times the MicroBooNE 372 second exposure. The selection efficiency indicates the ratio between the "Events" and "Events (after finalselection)".