Progress and Simulations for Intranuclear Neutron – Antineutron Transformations in 4018 Ar

With the imminent construction of the Deep Underground Neutrino Experiment (DUNE) and Hyper-Kamiokande, nucleon decay searches as a means to constrain beyond Standard Model (BSM) extensions are once again at the forefront of fundamental physics. Abundant neutrons within these large experimental volumes, along with future high-intensity neutron beams such as the European Spallation Source, offer a powerful, high-precision portal onto this physics through searches for B and B − L violating processes such as neutron–antineutron transformations (n → n̄), a key prediction of compelling theories of baryogenesis. With this in mind, this paper discusses a novel and self-consistent intranuclear simulation of this process within 40 18Ar, which plays the role of both detector and target within DUNE’s gigantic liquid argon time projection chambers. An accurate and independent simulation of the resulting intranuclear annihilation respecting important physical correlations and cascade dynamics for this large nucleus is necessary to understand the viability of such rare searches when contrasted against background sources such as atmospheric neutrinos. Recent theoretical improvements to our model, such as the first calculations of the 40 18Ar intranuclear radial annihilation probability distribution and the inclusion of a realistic n̄A potential, are discussed. A Monte Carlo simulation comparison to another publicly available n→ n̄ generator within GENIE is shown in some detail. The first calculation of 40 18Ar’s n → n̄-intranuclear suppression factor, an important quantity for future searches at DUNE, is also completed, finding T R ∼ 5.6× 10 s−1.


I. INTRODUCTION
A. Theoretical Background Baryon number (B) violation is the only remaining component of the Sakharov conditions [1] which has yet to be confirmed experimentally, an innate requirement for proper baryogenesis. While it was shown [2] that B itself is not an exact symmetry of the Standard Model (SM) and is only infinitesimally violated via B − L nonpertubative electroweak processes, this infraction does not appear to explain the excess of matter over antimatter we observe. Thus, it is possible that not only B must be violated, but that B − L must too be broken via some beyond Standard Model (BSM) extension.
Such explicit B violating extensions can be constructed as particular types of (di)nucleon decays, mainly within the structures of dimension d = 6 and d = 9 operators. While it is believed that many d = 6 operators (governing things such as proton decay) may be heavily (or altogether) suppressed [3], this may not be the case for d = 9 * jbarrow3@vols.utk.edu † golubeva@inr.ru ‡ paryev@inr.ru § j-m.richard@ipnl.in2p3.fr operators with the basic structure c · qqqqqq/M 5 , where the high-mass scale M can be rather low at perhaps ∼ 100 TeV. These operators allow for modes of dinucleon processes of disappearance or oscillation, a most important member of which is neutron-antineutron transformation (n →n). Popular (modern and minimal) BSM extensions [4][5][6][7] permitting such a process can dynamically create a proper baryon abundance in the early universe, even while predicting a reasonable and possibly reachable upper limit on the mean transformation time in vacuum, allowing the theory to be well constrained if not eventually entirely eliminated experimentally. This prospect seems stronger than ever given new studies from lattice quantum chromodynamics calculations [8] and other recent works [9,10].
B. DUNE and Intranuclear n →n in 40 18 Ar The Deep Underground Neutrino Experiment (DUNE) will be the future heart of American particle physics. It will contain 40 kt of liquid argon ( 40 18 Ar) inside its collective fiducial mass, acting both as prospective neutrino target and ionization detection medium within a set of four gargantuan time projection chambers. While its main goals of study are eponymous, some nuclear and particle physicists see DUNE as the next large-scale step arXiv:1906.02833v1 [hep-ex] 6 Jun 2019 toward answering arguably more fundamental and controversial questions about the nature of matter and the laws which govern it and its origins. Of course, BSM physics searches like n →n fit within this program well.
This oscillatory or transformational process can be best encapsulated within a single value: the mean vacuum (free) transformation period, τ n→n . With the inclusion of a so called intranuclear suppression factor (taking into account the bound nature of the n), this free transformation time becomes T = T R τ 2 n→n , where T R is quite large. The complications of this derivation, along with the first calculation of this factor for the 40 18 Ar nucleus, will be described in some detail in later sections.
Previous searches for n →n have been carried out using both free [11] and bound [12][13][14][15][16][17][18] n's. The most successful thus far has been Super-Kamiokande's [16] bound water-Cerenkov n →n search within 16 8 O, which, given 24 candidate events and an expected background of 24.1 atmospheric neutrino events over ∼ 4 years within 22.5 kT of fiducial mass, set a limit on the intranuclear n lifetime of 1.9 × 10 32 years, which, when converted to a free transformation time, became ≥ 2.7 × 10 8 s. The free search at the Institut Laue-Langevin [11] puts a lower limit at the same order of magnitude with no apparent backgrounds.

C. Past Intranuclear n →n Simulation Work
It is critical to recognize the interdependency of the computational modeling of BSM signals and backgrounds in the estimation of detector efficiencies and background rates within the context of large modern experiments. Thus, it becomes crucially important that one should not only take care to model the BSM signal and backgrounds as completely, consistently, and rigorously as possible, but also to employ limited approximations which still attempt to phenomenologically preserve as many innate physical correlations as achievable. This is especially true given nontrivial automated triggering schemes planned for future rare searches. Unfortunately, in the case of many previous n →n studies, this has not altogether been the case.
While previous work in 40 18 Ar by Hewes [19] and others from the DUNE Nucleon Decay Working Group using newly-constructed modules within the GENIE [20] Monte Carlo (MC) event generator has made excellent progress and developed technically fantastic analysis schemes using convolutional neural networks and boosted decision trees, many of the underlying physics assumptions of a hypothesized n →n signal within the GENIE default model are not entirely correct. Some approximate notions include: 1. The assumption that the annihilation occurs along the nuclear density distribution of the nucleus, even while [21] is openly discussed; in later sections, we shall compare this assumption with a more holistic quantum-mechanical model where there appears to be moderate disagreement.
2. Employing a single nucleon momentum distribution described by a non-local, relativistic Bodek-Ritchie Fermi gas (the radial dependence of the annihilating (anti)nucleons' momenta are ignored).
3. Only ∼ 10 annihilation channels (a la [16]) are assumed to be necessary to describe the annihilation products. This seems low, as ∼ 100 are known, many of them containing heavier resonances; these heavier species can be responsible for ∼ 40% of all pion (π 0,± ) production. 7. Only an rough estimation of the nuclear suppression factor of 40 18 Ar has thus far been used to approximate lower limits on the transformation period.
It is no doubt that some of these current technicalitie proceeds directly from the secondary nature of the GE-NIE n →n module's genesis, a consequence of GENIE's top-down structure and first-and-foremost focus on neutrino interactions. This being said, similar issues or inconsistencies are known to exist in other work [22], and detailed explanations of past simulation's internal processes are lacking [16,23], but we do not intend to litigate these here. Instead, we offer a new and antinucleondata-tested model which does not ignore these effects. Secondarily, we have taken the liberty to investigate and slightly perturb GENIE's publicly available n →n module in an attempt to understand and show commonalities and differences across simulations, allowing for discussion of how these may in turn effect the eventual feasibility of a true event's observation in DUNE.

A. Concepts and Pertinent Questions
Assuming that the n →n transformation does occur in a vacuum, an interesting question arises: what are the consequences for a nucleus? Of course, if an intranuclear n becomes ann and collides with another nucleon, an annihilation will take place, releasing ∼ 2 GeV of rest-mass energy, from which ∼ 4-5 mesons are emitted on average. After this point, the wounded nucleus will evaporate several nucleons and perhaps break into unstable daughter nuclei. Several questions are immediately raised: 1. When a n tentatively becomes ann, it ceases feeling a smooth potential of 50 MeV and instead experiences a (complex) potential whose magnitude is 100 MeV. How much is such a transformation suppressed by this change in potential? 2. A deep annihilation could produce multiple fragments with the primary mesons ultimately being absorbed. Alternatively, a peripheral annihilation would probably release a large fraction of the primary mesons and at most rip out only a few nucleons, albeit with a more asymmetric topology. So, where, preferentially, does the annihilation take place?
4. Some concerns have been expressed about the reliability of the estimate of the n →n lifetime inside nuclei within a straightforward nuclear-physics approach based on shell-model physics; see, e.g., [31]. Can one face these criticisms and demonstrate stable and consistent results?
In this section, we review these questions with methods based on the Sternheimer equation (see, e.g., [32], and refs. therein), as used by Dover, Gal and Richard [21,[33][34][35] throughout past discussions of these topics.
We then apply such methods to the 40 Ar isotope, which comprises the main component of the DUNE detector. For completeness, we briefly repeat the main steps of these calculations.

B. Lifetime of the Deuteron
As a warm-up, consider a simplified deuteron, consisting of a pure s-wave bound state of a proton (p) and n. We adopt the wave function by Hulthen which has been tuned to reproduce the correct binding energy and spatial extension. It reads (see, e.g., [38]) where N n is a normalizing factor, r is in GeV −1 , λ 1 = 0.2316 c and λ 2 = 5.98 λ 1 . It has been checked that using another wave function does not change the following results significantly, provided it fits the deuteron energy and radius.
In the presence of n →n transformations, the wave function becomes Assuming an arbitrary strength γ = 1/τ for the n →n transition, the inducedn component w is given by the Sternheimer equation (at first order in perturbation theory) where E is the unperturbed energy and V the antineutron-proton potential, resulting in a width This immediately implies the scaling law Γ ∝ γ 2 , or for the lifetime T of the deuteron where T R , sometimes called the "intranuclear suppression factor", is actually a reduced lifetime. For the optical potential V (r), there are some models that are tuned to reproduce the main features of low-energy antinucleonnucleon scattering and have predicted the shift and broadening of the low-lying levels of the protonium atom. A method to solve Eq. (3) is discussed in [21], which is similar to the one used in the earlier studies [33,34], and involves the matching of several independent solutions corresponding to various limiting conditions. The alternative below directly provides the desired solution. First, to get the neutron wave function from the neutron potential V n (r), one should solve the radial equation subject to u n (0) = u n (∞) = 0. A method adapted from aircraft engineering [36] consists of the change of variables r = r 0 x/(1 − x), where r 0 is a typical distance, and, for the wave function u n (r) =ũ n (x), of an expansionũ in which the coefficients a j are closely related to the values of the function at the points x i = i /(N + 1), for i = 1, 2, . . . N . This results in a N × N eigenvalue equation A U n = E U n , where A is the discretized Hamiltonian and U n is the vector of the (1 − x i )ũ n (x i ). See, for instance, [37], for an application to quarkonium in potential models. For the Sternheimer equation (3), one gets a simple matrix equation where the N ×N matrixĀ is the discretized Hamiltonian for then, E the n energy, and W = {(1−x i )w(x i )} is the vector containing then wave function. This calculation is fast and robust. Besides the deuteron energy E 0 = −0.0022 GeV and wave function u n , solving Eq. (3) requires a model for the antineutron-proton potential V . As shown by Fermi and Yang [39], the long-range part of theN N potential in isospin I is deduced from the N N interaction in isospin I by the G-parity rule: if a meson (or set of mesons) with G = +1 is exchanged, it gives the same contribution, while under G = −1 exchange, the sign is flipped. The complex short-range part of V is fitted as to reproduce the low-energy data onp scattering and protonium. We have adopted the so-called DR2 model [40], and it was checked that variants such as the Kohno-Weise potential [41] produce very similar results.
From these considerations, one can construct an example of then radial position distribution, as shown in Fig. 1. Note that all such distributions in this and later sections have been implicitly multiplied by r 2 , unless otherwise stated. However, an annihilation requires both the presence of a neighbor nucleon and ann, so the annihilation must take place at a slightly less peripheral site than then density itself. For the deuteron, the reduced lifetime is estimated to be about T D R 3 × 10 22 s −1 . This result shows great consistency with other recent calculations [42]. This estimate is remarkably stable. For instance, increasing the core of thenp interaction by a factor 10 results in only a 20% increase of T R . Even with a large |V |, then transformation is more suppressed, but it actually annihilates more efficiently. Remarkably, there is an almost exact cancellation between these two effects. It can also be seen that the calculation is sensitive mainly to the value of V (r) near 0.8 − 1 fm. This is fortunate, as low-energyp scattering on nucleons and the shift of the antiprotonic hydrogen atom probes essentially this region and so one cannot determine the interaction at closer distances. 1

C. Lifetime of the 40 Ar nucleus
In Refs. [21,33], there are estimates of the lifetimes of some nuclei important for past underground experiments, such as 16 O [16] or 56 Fe [15]. Here, we repeat the exercise for 40 Ar.
The detailed properties of atomic nuclei are well accounted for by sophisticated Hartree-Fock calculations. For many applications, it turns out to be rather convenient to use ad-hoc shell-model wave functions that are tuned to reproduce the main properties of the nuclei, in particular the spatial distribution of p's and n's. This was done in connection with the compilation of nuclear data, see, e.g., [44]. In the present study, we follow the strategy outlined in [45]. The p and n wave functions have been calculated for us by Karim Bennaceur, in the so-called "filling approximation": the nucleus is supposed to be spherical, implying that the states of each shell are populated with the same (integer or fractional) occupation number.
For each n shell, there is an inducedn wave-function governed by an equation analogous to (3), with a centrifugal term for p, d and f states, where V is now thē n-39 Ar potential, and m the corresponding reduced mass. In Fig. 2, some details are given for one of the external shells which contributes most to the instability, namely 1d 5/2 . It is seen thatn are produced in the tail of the n distribution, with subsequent annihilation at the surface of the matter distribution.
For comparison, the distribution of the 1f 7/2 shell are shown in Fig. 3: the peripheral character of the antineutron component is even more pronounced. The resulting radialn distributions for all shells are shown in Fig. 4.
If one now adds up the contributions of each neutron to the width, calculates the average width per n, and estimates the corresponding reduced lifetime, one gets a value of T Ar R ∼ 5.6 × 10 22 s −1 . As in the case of the deuteron, this value is remarkably stable against changes in the parameters of then-nucleon interaction. Thus, we similarly estimate the uncertainty to be about 20 %. This stability in the width can be understood: if one increases the absorptive potential, the n →n transition is more suppressed, but, on the other hand, the antineutron annihilates more efficiently. In Fig. 5 we show the factor, γ, by which the width of the 1d 5/2 level is modified when the realn-nucleus potential is multiplied by f r and its imaginary part by f i . If one changes these values by ±20%, far beyond what can be admitted to keep a good fit to the antinucleon-nucleus data, modifications to the width are less than 10%.
FIG. 5. The relative change to the width of the 1d 5/2 shell is shown when a factor fr is applied to the real potential, and fi to the imaginary part.
One can also illustrate how intranuclear n →n transformation and subsequent annihilation is a surface phenomenon. In Fig. 6, the absorptive potential is modified near r = r c by applying as a factor a "glitch" function with the form 1+0.4 exp(−20 (r−r c ) 2 ), where r c is varied. The width is seen to be modified only near the surface, and is insensitive to what happens in the center of the nucleus. To end this section, it is instructive to compare the method based on an exact solution of the perturbation equation (3) and the usual approximation where the n andn spatial distributions are not distinguished. If, fur-thermore, the n-andn-nucleus do not differ too much in their real part, then the width is given by (see, e.g., [43]): For the low-lying shells of 40 Ar, the difference is small with respect to our estimate. For the most external shell, the width is overestimated by a factor of about 3.5, overemphasizing the suppression of the process.

III. INTRANUCLEARn DYNAMICS
A. The need for a more precise accounting of nuclear effects As a rule, the Intranuclear Cascade (INC) model used for simulations of inelastic interactions of particles within the nucleus is assumed to hold at energies 30-40 MeV (the conditions of applicability of the model are considered in detail in [30]). The influence of other intranuclear nucleons on an incident particle is taken into account by the introduction of some averaged potential U (r), and so within the nucleus a cascade particle changes its energy by the amount of this potential. In the local Fermi gas approximation used within the INC model, the intranuclear nucleons are bound in the nucleus by a nuclear potential V N (r) = −T F (r) − N , where T F is the Fermi energy of the nucleon, and N the average binding energy per nucleon. Since the energies of the particles participating in the cascade are sufficiently large, a simplified account of the influence of the nuclear environment is justified and does not lead to distortion of the simulation results However, after modeling the interaction of a very slow n with a 12 C nucleus [30], some questions arose about the legitimacy of particular physical approximations which could no longer be ignored, requiring that we attain a more correct accounting of the nuclear environment in both cases of an extranuclear and intranuclear transformation and subsequent annihilation. We now discuss these pertinent changes in our model, which have been incorporated into all extranuclear and intranuclear simulations.

B. Then-nucleus intranuclear potential and hown modifies the nuclear medium
The influence of the nuclear environment on an incoming extranuclear (in the case of ESS using 12 C [30]) or intranuclearn leads to the modification of then's vacuum four-momentumpn = (En, pn) (p 2 n = E 2 n − p 2 n = m 2 n ) inside the target nucleus due to an effective scalar attractive nuclear potential of the form [54] where ρ(r) and V 0 together are the local nucleon density normalized to the atomic mass number A of the nucleus, then potential depth at saturation density is ρ 0 , and r is the distance between then and the center of the nucleus. We may assume that thenA andpA potentials are the same. With this potential, the totaln energy E n in the nuclear interior of ordinary nuclei can be expressed in terms of its in-medium mass m * n , defined as [54,55] 2) m * n (r) = m n + Un(r) (12) and its in-medium three-momentum p n , as in the free particle case, is [55] E n = m * 2 n + p 2 n .
Analysis [55] ofp production in proton-nucleus and nucleus-nucleus collisions at kinetic energies of a several GeV showed that thep potential at normal nuclear matter density is in the range of −100 to −150 MeV for outgoingp momenta below 2.5 GeV/c. Studies ofp production at AGS energies [57,58] suggestp potentials of −250 MeV and −170 MeV at density ρ 0 forp annihilation events at rest with respect to the nuclear matter and forp with momentum of 1 GeV/c, respectively. The real parts of anp optical potential in the center of the nucleus of −(150 ± 30) MeV and of −(220 ± 70) MeV were extracted in [59] from the data onp absorption cross sections on nuclei and on the annihilation spectra of π + 's and p's, correspondingly. Combined analysis [60] of data on antiprotonic X-rays and of radiochemical data showed that at the center of the nucleus thep potential is approximately −110 MeV in depth. So, in spite of various attempts to fix this potential, its depth at density ρ 0 still remains rather uncertain presently. For the sake of definiteness, in the subsequent calculations, we will adopt the realistic value V 0 = −150 MeV within Eq. (11).
The in-medium momentum p n is related to the vacuum momentum pn by the following expression: For example, with V 0 = −150 MeV, this shows that for n annihilation at rest, i.e. when pn = 0, then momentum |p n | in the center of the nucleus and at its periphery, corresponding to 10% of the central density, is equal to 510 and 167 MeV/c, respectively.
Within the non-interacting local Fermi-gas model used in our MC simulations [30], for the bound target nucleon total energy E N in the medium at the annihilation point r we will employ the formula [61,62] 2 ) The potential Un is not the usual Lorentz scalar potential UN S , determining along with the Lorentz vector potential UN V , the total in-medium antinucleon energy E N via the dispersion relation with Here, p N is the momentum of the nucleon N (N = {p, n}) in the Fermi sea, p F (r) is the nucleon Fermi momentum at the local point r, and the quantity N ≈ 7 MeV is the average binding energy per nucleon; at this point, 0 ≤ |p N | ≤ p F (r). Within the representation of Eqs. (13)- (15), the invariant collision energy s for the interaction of ann with a nucleon bound in the nucleus at the point r is The total collision energy En + E N entering into the second relation of Eq. (18) in the non-relativistic limit appropriate to our case can be calculated as Contrary to the on-shell interaction, forn annihilation at rest, from Eq. (19) it is seen that this energy is always less than m n + m N and its maximum value is m n + m N − N .
If a bound target n is transformed into ann (for the 40 18 Ar case), we assume that its total energy E n , defined by Eq. (15), is equal to that E n of then, determined by Eq. (13) above. Namely: It is interesting to note that, for p F (0) = 250 MeV/c, Eq. (20) gives for then momentum |p n | the values of 430 and 40 MeV/c if the transition n →n of the target n at rest occurs in the center of the nucleus and at its periphery, respectively (corresponding to 10% of the central density).
The value of s for the intranuclear 40 18 Ar case is given by the first relation of formula (18).

IV. MONTE CARLO SIMULATION METHODS
In the preceding sections, we have attempted to summarize new theoretical additions to the dynamics of our MC simulation. We will now review the stages of our generator, accompanied by some discussion of other recent changes given new antinucleon annihilation data, and show some differences from previously published outputs [30].

A. Fundamentals
Most all of the background concerning our MC generator can be found in our recent work [30], however, we will mention some improvements which have been implemented here. We can briefly summarize the content of the model as follows: 1. An annihilation point is taken from a corresponding probability distribution as shown later in this work for 40 18 Ar within Fig. 12. This point lies within a particular zone of a set of eight concentric spherical shells representing the volume of the nucleus. Each shell has its own uniform nuclear density and each their own single nucleon momentum distribution, thus acting as a zoned local Fermi gas.
2. The annihilation occurs, producing π's and higher mass resonances such as η, ω, and ρ's, according to tabulated channels with branching ratios taken from [30]. After the decay of all resonances, on average ∼ 4-5 mesons are produced from the initial annihilation process.
3. The annihilation products are then transported through the nuclear environment quasi-classically using a full intranuclear cascade model (a local nuclear density decrease is also included). Meson resonances are decayed according to their individual lifetimes inside and outside the nucleus, though particularly long-lived species are treated as stable inside the nucleus.
4. The products are ejected from the nucleus and the nuclear remnant(s) is allowed to de-excite, evaporating nucleons and fragments of higher mass.

B. Improvements from OBELIX and Crystal Barrel Data
The annihilation model was created in 1992 [26], originally using experimental branching ratios obtained before this time. Later analysis of experimental data from pp annihilation at rest obtained from the LEAR (CERN) p beam by the OBELIX [46] and Crystal Barrel [47] collaborations are now integrated into the internal annihilation model. Note that all tables in all proceeding sections are outputted from samples of 10,000 events. In Table I we see the absolute changes in the branching fractions of individual annihilation channels in accordance with this newer experimental data. The first column shows the annihilation channel, the second column shows the value of this channel's branching fraction in the corresponding table in [30], while the third shows the new branching fraction of this channel changed in accordance with the experimental data [46,47]. Once summed, all channels are then renormalized to unity with these new fractions taken into account. A comparison of the elementary processes forpp annihilation was carried out taking into account experimental data [46,47] and annihilation simulations from our recent work [30]. Table II shows the average multiplicities of mesons produced inpp annihilation at rest. The first column re-presents the simulation results from Table IV in [30], while the second column presents the results of modeling when taking into account data from [46,47]. The third column also re-presents the experimental data itself [30,[48][49][50][51] for ease of comparison. It follows from Table  II that the average multiplicities of annihilation mesons with changes in the branching ratios of some individual channels do not change significantly and the difference between the two calculation options is within the uncertainty of the experimental data. Nevertheless, in this and all future simulations, this new table will be used for modeling allpp annihilations.  [47] pp → π 0 π 0 π 0 ω 0.40 1.24 [47] a OBELIX gives data forpp → π + π + π − π − with a branching ratio of BR(pp → π + π + π − π − ) = 6.4 ± 0.09% (a gas target). This value includes both resonant and independent production of these four charged π's. To avoid double counting, we have subtracted from this value all the fractions of all of the channels that give rise to π + π + π − π − in the final state. Thus, in thepp annihilation table contained in [30], we have experimental data showing BR(pp → ρ 0 ρ 0 ) = 0.67%, BR(pp → π + π − ρ 0 ) = 2.02%, and importantly BR(pp → π + π − ω) = 3.03%, whose contribution to the full channel ofpp → π + π + π − π − (including BR(ω → π + π − ) = 2.3% is 3.03 × 0.023 = 0.07. From all of these considerations, we introduce the final ratio to the table as 6.4 − 0.67 − 2.02 − 0.07 = 3.64%. Table II shows the average multiplicities of mesons produced inpp annihilation at rest. The first column re-presents the simulation results from Table IV in [30], while the second column presents the results of modeling when taking into account data from [46,47]. The third column also re-presents the experimental data itself [30,[48][49][50][51] for ease of comparison. It follows from Table   II that the average multiplicities of annihilation mesons with changes in the branching ratios of some individual channels do not change significantly and the difference between the two calculation options is within the uncertainty of the experimental data. Nevertheless, in this and all future simulations, this new table will be used for modeling allpp annihilations.

C. Model modifications and a new comparison with experimentalpC annihilation at rest
The branching fraction modifications shown in Table I were implemented into the optical-cascade model, most recently discussed in [30]. Now, let's analyze how these changes effect the description of experimental data forpC annihilation at rest. The first line of Table III presents the experimental multiplicities of the emitted π's and the energy carried away by those π's and γ's. The second line of the table shows the results of a calculation with our original optical-cascade model from [30] before any modifications. The third line (Calculation #1) presents the results of a simulation taking into account the changes in the annihilation table described above. The last line (Calculation #2) presents the results a simulation accounting for both the changes in the annihilation table and modifications related to the nuclear environment. In both cases, an antinucleon potential of 150 MeV in the center of the nucleus, and varying in accordance with the nuclear density, is included. , and two new calculations taking into account the newest versions ofp annihilation branching ratios while also considering a new intranuclear antinucleon potential with an associated nuclear medium response. From analysis of the table, it follows that the average multiplicities vary slightly with the modification of the annihilation table and dynamics. Multiplicities of π's increase slightly with the modifications associated with the influence of the nuclear environment, and everywhere are within the experimental error. At the same time, the multiplicity of p's and n's emitted during the cascade development and emitted in the process of de-excitation was significantly reduced in the Calculation #2. This suggests that the dynamics of meson-nuclear interactions and energy dissipation in the residual nucleus change somewhat with the introduction of the influence of the environment. Clarification of this issue requires data on proton and neutron emission frompN annihilation experiments, along with further detailed study. Since the present work is devoted to experiments to search for transformations that are planned to be registered through the π-stars, we will focus on the characteristics of these topologies and their associated quantities. Fig. 7 shows the spectrum of π + emitted duringpC annihilation at rest in comparison to experimental data from [49] (triangles), and [53] (circles). Just as in the average multiplicities shown in Table III, there are no significant differences in the spectrum for these two variants of calculation, though a comparison to the old calculation reveals a slightly better fit to the data in the region around 250 MeV (the ∆-resonance).
In our calculations for the 40 18 Ar nucleus to be discussed in later sections, we will use a model that includes effects related to the influence of the nuclear environment.  Table  III, while the solid shows Calculation #2. All points are taken from experimental data in [49,53].

A. Changes in extranuclearnC annihilation simulations
Despite the absence of significant differences in the description of availablepC experimental data, there are certain aspects of the two variants of the models where differences are very noticeable. The changes discussed above lead to some rather important differences in annihilation stages as compared to previous work in [30], examples of which can be seen in Figs. 8-11. These can make possible different observable final states after the intranuclear cascade which are important for experiments planning to utilize smaller nuclei (such as the NNBar Collaboration at the European Spallation Source). Note that all plots in all proceeding sections are outputted from samples of 10,000 events.
For the intranuclear cascade, the operating energy conservation law for the annihilation process of an extranuclear neutron is written as where E n is the total energy of then inside the nucleus at the point of annihilation, E N is the total energy of the nucleon annihilation partner at the same point, and E * is the excitation energy of the nucleus after the annihilation. In the degenerate Fermi gas model, this is defined as varying from 0 to T F N i , where i is the zone number in which the annihilation takes place, T F N i is the boundary Fermi energy of the i-th zone, and T N i is the Fermi energy of the annihilation partner in the same zone. If one takes into account relations (14)- (16) and (19) from Sec. III, it follows that where N = 7 MeV/nucleon. In Fig. 8, we see the distributions (old and new) of the initial amount of total energy carried by the two annihilating nucleons. Due to the fact that the previous iteration of our calculations did not include a truenpotential or consider intranuclear nucleons to be off-shell, the skew of this distribution was always greater than the available rest-mass energy of an annihilating pair; this was discounted as a kind of virtual intranuclear effect, which later disappeared and showed conservation of energy and momentum once entering the final state. This is now changed and made internally consistent, showing a proper distribution less than the combined free rest masses of the pair of annihilating nucleons. Secondarily, the strength of then-potential can be seen to smooth out the zoned structure present in the previous calculation.
The dynamical (position-correlated) momentum of the initial annihilation pair (and, by conservation, their annihilation products) can be studied in Fig. 9. In the old variant of the model, then was assumed to come from a transformation down-range of a cold n source with a mean energy of only ∼meV, ignoring then-potential. Thus, the original momentum distribution of the annihilation products (dashed line) was effectively a direct observation of the non-interacting zoned local Fermi gas single nucleon momentum distribution folded with the radial annihilation probability distribution. In the new variant of the model described above, shown in the solid histogram, the mass of then is defined by expression (12) and the momentum of then follows from (14). The direction of the momentum of the nucleons are isotropically distributed, and thus the total momentum of the annihilation products varies in absolute value from |P n − P N | to |P n + P N |, smoothing and spreading out the structure associated with the presence of zones in the target nucleus. The peak in the histogram in the region just below 100 MeV/c corresponds to annihilations on the outside of the nucleus (within the diffuse eighth zone), where then-potential is taken to be 0 while no off-shell mass is accounted for, and so the momentum of the annihilating pairs is equal to the momentum of the nucleon partner within the seventh zone. Here we see the total available initial and final mesonic/pionic and photonic parameter space for an n →n signal, most commonly shown via total momentum versus invariant mass plots (similar to [16]). The results of new modeling are shown at top, and the old version of the model is shown below. In Figs. 10, the effect of the antineutron potential and off-shell nature of nucleon masses are clearly seen in the top plot, while the bottom plot shows a small momentum range due to the absence of the antineutron momentum; the right-ward inclination of the old parameter space is a consequence of the on-shell mass' effects on the overall kinematics of the annihilation. Thus, we have significantly different initial conditions for the transport of annihilation mesons through the nucleus. Figs. 11 show the same variables after transport, but re-scattering, ∆-resonance, and absorption of annihilation mesons leads to an overall decrease in the observed invariant mass, lessening the apparent final state differences between the two simulation variants. Note the two-lobes seen in the leftward-reaching blob-a sign of single initial meson absorption. While this structure is still somewhat apparent in the new plot, it is not as prevalent; the likely cause of this is the ability for thenpotential to comparatively give momentum to the annihilation products more effectively, leading to slightly fewer total meson absorptions while changing the expected momentum distribution. Overall, the new parameter space is slightly more spatially if not statistically constrained, which could lead to higher hypothetical experimental efficiencies for future experiments. The new (top) and old (bottom) total mesonic/pionic/photonic final state parameter space is shown for extranuclearnC annihilation.

B. IntranuclearnAr annihilation simulations
As mentioned previously, there now exist multiple generators for n →n within the particle physics community, notably developed in [19] using GENIE [20]. While we believe that the demonstration of our independent generator's capabilities [30] in the reproduction of antinucleon data is well-established for 12 C, such a complete set of physical observables does not readily exist for larger nuclei, especially not 40 18 Ar. Thus, out of a need for ample comparisons, we endeavor to show the commonalities and differences between each of these n →n generators for intranuclearn 39 Ar annihilation useful to DUNE. We do this by attempting to make some of the same assumptions (roughly) as GENIE, and vice-versa. For instance, we can and do generate events by simulating the annihilation position sourced from a Woods-Saxon distribution within our generator (alongside the more modern version as developed in Sec. II); similarly, with little work, we have perturbed the default settings of the GENIE n →n generator module to utilize a noninteracting local Fermi gas nuclear model along with a full intranuclear cascade. The inclusion of ann-potential within GENIE has not yet been investigated; implementation of the modern annihilation position probability distribution is currently underway. While none of these comparisons across generators are ever to be exact, their approximately equivalent formalism can eventually serve to inform us on stability of quantities which characterize the possible final state topologies of a true n →n event signal with respect to their associated backgrounds. This stability across models, and their interplay with model detector reconstruction, will be studied in detail in future work with DUNE collaborators.
Two of the probability distributions of intranuclear radial position upon annihilation for these generators are shown in Fig. 12 in orange and dashed blue, and are surprisingly similar even with quite different physical assumptions. The quantum mechanical, shell-by-shell distributions discussed in Fig. 4 are all taken in a weighted average to create the final orange curve, from which our generator can source its initial annihilation positions in a binned fashion. All GENIEv3.0.2 events source their positions from a realistic, smooth Woods-Saxon (nuclear density) distribution [20] incredibly similar to our parameterization for the nuclear density: When this curve is multiplied by r 2 , one generates the dashed blue annihilation probability distribution. These curves effectively demonstrate how even the most simple of assumptions, many only quasi-classical, can lead to quite good approximations; however, we note the preponderance of events even further toward or beyond the surface of the nucleus using a fully quantum-mechanical formalism. The increased likelihood of such surface annihilations, along with their associated correlation with lower momenta and higher final state meson multiplicity, will be shown in the coming figures to be an arguably critical part in the proper evaluation of future experimental efficiencies and possible lower limits on mean intranuclear n →n transformation time. For similar plots and discussions, see [26][27][28][29][30]. Some of the differences between this work and GE-NIEv3.0.2's generator pertain to the initial dynamics of the intranuclear annihilation. An example of this can be seen in Fig. 13, showing the initial annihilation meson total energy for this work (using the orange curve in Fig. 12) and GENIE (using the dashed blue curve in Fig. 12), each using a version of a local Fermi gas nuclear model. Like the extranuclearn annihilation described In orange, we see the modern, quantum-mechanically derived, shell-averaged, true intranuclear radial position of annihilation probability distribution as developed in Sec. II, present in our generator. Probability distributions are normalized to the same arbitrary integral for a direct comparison, and the scale is arbitrary. on 12 C above, energy balance in the annihilation point is given as E ann + E * = E n + E N ; taking into account that we have E n = E n , and that we have the total energy E ann + E * = m n + m N − 2 = 1.866 GeV .
It can be seen in Fig. 13 that our distribution of energies available to an annihilation is always less than 1.866 GeV, then dynamically changes with radius. However, the GENIE distribution (which assumes a similar binding energy per nucleon as our model) does not vary with radius; thus, the total range of the annihilation energy within the GENIE simulation can be explained simply by considering the minimum/maximum potential magnitudes of annihilation pair momentum (corresponding to an anti/co-parallel intranuclear collision) while assuming an approximately constant intranuclear defect nucleon mass of ∼ 910 MeV/c 2 . The sharp rise of the GENIE distribution around 1.82 GeV can be seen to correspond to the addition of momentum distribution shapes around zero momentum (see Fig. 14).
One can see the different initial single nucleon momentum assumptions in Figs. 14. The GENIE non-local Bodek-Ritchie or local Fermi gas nuclear models mentioned here serve effectively only as a set of initial conditions which enable certain nucleon momentum and radius correlations (or lack thereof). The non-local Bodek-Ritchie has been considered the default operating model for n →n simulations; however, for most of the rest of this article, we will compare local Fermi gas models to each other for a succinct simplicity. Note the different characteristic ranges of momenta; in general, the shapes and ranges of each nucleon local Fermi gas model (solid lines, top and bottom figure) are incredibly similar, while the non-local Bodek-Ritchie is quite unique (dashed lines), especially with its phenomenological, short-range "correlation" tail. For GENIE, we see that the shapes of all distributions are identical: p f (n) = p f (n) ≈ p f (p); this is not the case for our model, where p f (n) = p f (n) ≈ p f (p) due to then-potential. 3 The most important aspect of correlated behavior which has been previously unaccounted for in GENIEaffiliated work on n →n is that of initial (anti)nucleon momentum and radius. In Figs. 15, we show comparison between our and GENIE's local Fermi gas nuclear models, which by definition preserves these correlations. All figures assume a (zoned or smooth) Woods-Saxon annihilation probability distribution, and outputs are rather similar in our simulations using work from Sec. II. Thē n-potential is apparent in the top plot, which appears smoothed and lacking of any zoned discontinuities due to the strength of the interaction (as seen in the middle plot, showing correlations forn annihilation partners). All plots correctly predict a falling-off of nucleon momentum at higher radii, a key consideration for event reconstruction and background rejection. This behavior is not present within GENIE when using the default Bodek-Ritchie model.
The initial annihilation meson total momentum distribution is seen for our and GENIE models in Fig. 16; these distributions are equivalent to the initial two-body annihilation pair momentum distributions by conservation. Each histogram shares a Gaussian-like shape due to the randomized momentum selection from underlying distributions, though our output shows much higher available momentum due to the interaction of then with the modified nuclear medium.
All of this leads us to consider the available initial (and final) mesonic parameter space. As seen in Figs. 17, this serves as an initial condition of the annihilationgenerated mesons before intranuclear transport (thus, for GENIE, hA/hN2018 models are at this stage equivalent). Due to the non-dynamical off-shell masses of annihilation pairs, GENIE predicts higher invariant masses, while our model shows them decrease due to off-shell mass defects in correlation with radial position. Overall, the space is quite differently filled for each model, though considering this is before the intranuclear cascade, one cannot necessarily predict much about the final state.
The follow-up to these figures can be studied in the comparison of Figs. 18, where our and GENIE gener- ated events proceed through a full intranuclear cascade (hN2018 only is shown here). Note that our model includes photons in the final state from high-mass reso-nance decays, while GENIE does not. The disconnected regions toward the left of the plots are signs of single π emission after at least one or more meson absorptions. While overall the distribution of events is rather consistent, critically, the high density of events with large invariant mass and low momentum (bottom right of plots) among these local Fermi gas models shows the importance of modeling correlations between position and momentum as they imply a comparatively large number of escaping (and possibly visible) π's in the final state. It is with these areas that one may hope to find a significant rejection of background events, possibly allowing for an actual observation of an n →n event.
We end the comparison of these generators with Tab. IV and Figs. 19 and 20, which show many similarities and some differences across them. Multiplicities in Tab. IV are seen to be most different between our model and GENIE in the realm of outgoing nucleons (resulting from nucleon knock-out or evaporative de-excitation); this should not be surprising, as GENIE does not currently contain a public version of an evaporation model within either its hN2018 (full intranuclear cascade) or hA2018 (single effective interaction) model tunes. Exciting work by GENIE developers in this regime is expected to be completed soon, and we look forward to being able to compare our results. Small differences can also be seen in the π 0,± multiplicities, which are partially due to the fact that we predict more n →n events toward or beyond the surface of the nucleus (see the blue and orange radial annihilation position probability distributions of Fig. 12), but there are also nontrivial dependencies given the larger number of possible branching channels we simulate [30] compared to GENIE [20].  To give a more complete context to Table IV, we plot several (absolute magnitude) final state momentum spectra for π + and p species, two key constituents in the eventual experimental search for n →n in DUNE. In Fig. 19, we see that our models agree quite well with the full intranuclear cascade simulation from GENIE (using hN2018) in both multiplicity and shape; there is some lack of shelf-structure around the ∆-resonance within the hA2018 simulation (recall this models the cascade as a single effective interaction), a sign of the competition between cross-sections (or rates) of processes such as ∆ decay and pion absorption.
In Fig. 20, we see the outgoing proton spectra. Here, our model and GENIE differ greatly (and GENIE even among itself) across lower momenta. Though there is a full intranuclear cascade model within GENIE (hN2018), it can be directly seen that it does not yet include any nucleon evaporation currently. In some respect, these differences should be expected due to the novel nature of the phenomena we are modeling (transporting ∼ 4-5 mesons is no easy business), and the fact that our generator was comparatively purpose-built to reproduce antinucleon annihilation data. We will note, however, that if one takes a more experimental viewpoint, these are not actually so disparate; indeed, if we consider an approximate, conservative, minimum proton kinetic energy de-tectability threshold in liquid 40 18 Ar to be 100 MeV (i.e., 450 MeV/c) [63], we see that above this value much of the shape and magnitude of all distributions are quite similar. Thus, in some respect, we expect these models to appear rather degenerate for protons when taking detector response into account.

VI. CONCLUSIONS AND OUTLOOK
In this work, we have endeavored to give an update to the community on recent developments in the modeling ofnA annihilation events in service of future BSM n →n searches. Critical among the findings of this paper is the calculation of a new intranuclear suppression factor for 40 18 Ar, T Ar R ∼ 5.6 × 10 22 s −1 , along with a new and associated calculation of the 40 18 Ar radial position annihilation probability distribution. Also, efforts have been made to implement this and other important nN annihilation dynamics into an independently developed, antinucleon-data-driven MC generator to service both the ESS NNBar Collaboration and DUNE. Samples of 100,000 events for both communities are now being prepared and will be made available upon request to the authors. Comparisons and discussions of differences and similarities have been made to data where available, previous MC results, and other publicly available event generators such as GENIE.
Within this work, a kind of forward path has been illuminated for the (intranuclear) BSM (di)nucleon decay community, showing the importance of some initial physical correlations in the modeling of BSM signals, most importantly that of the event constituent's momenta and intranuclear position. However, the effects of final state interactions cannot be understated. It is with these findings, and our associated event generator, that we hope to empower future experiments to better understand probable signal topologies for rare decays. Collaboration is ongoing with DUNE community members within the Nucleon Decay Working Group to study the implications of this modeling work on possible efficiencies, atmospheric neutrino background rejection rates, and lower limits on the n →n mean transformation time inside simulated DUNE detector geometries. for providing a stimulating environment where these topics have been thoroughly discussed. We also wish to thank Yuri Kamyshkov for his constant stimulus and interests in this work, helping to drive us all along.