A model of antineutron annihilation in experimental searches for neutron-antineutron transformations

Searches for baryon number violation, including searches for proton decay and neutron-antineutron transformation, are expected to play an important role in the evolution of our understanding of beyond Standard Model physics. The neutron-antineutron transformation is a key prediction of certain popular theories of baryogenesis, and the experiments such as the Deep Underground Neutrino Experiment and the European Spallation Source plan to search for this process with bound- and free-neutron systems. Accurate simulation of this process in Monte Carlo will be important for the proper reconstruction and separation of these rare events from background. This article presents developments towards accurate simulation of the annihilation process for use in a cold, free neutron beam for neutron-antineutron searches from antineutron-carbon annihilation, as carbon-12 is the target of choice for the European Spallation Source's NNBar Collaboration. Initial efforts are also made in this paper to perform analogous studies for intra-nuclear transformation searches in argon-40 nuclei, although this latter work will be properly addressed in our next publication.


INTRODUCTION A. Background
As early as 1967, A. D. Sakharov pointed out [1] that for the explanation of the Baryon Asymmetry of the Universe (BAU) there should exist interactions in which baryonic charge is violated besides mere departures from thermal equilibrium and symmetry. Thus, experimental searches for baryon number ( ) violating processes, and in particular the baryon minus lepton ( − ) number violating process of neutron-antineutron oscillation ( → ̅), are of great importance due to their possible connections to the explanations of the observed matter-antimatter asymmetry of the universe-as first laid out by V. A. Kuzmin [2] and followed in developments by many authors, see e.g. recent reviews [3][4][5].
Thus, the search for → ̅, along with nucleon decay, remains one of the most important areas of modern physics, hopefully leading to an understanding of phenomena related to the BAU.
The best lower limit on a measurement of the oscillation period with free neutrons, → ̅ , was attained at a reactor at the Institut Laue-Langevin (ILL) [6] in Grenoble, France, with a cold neutron beam. These neutrons flew through an evacuated, magnetically shielded pipe of 76 in length (corresponding to a flight time of ~0.1 ), until being allowed to hit a target of carbon ( 6 12 ) foil (with a thickness of ~130 ). This foil would have absorbed antineutrons, resulting in matterantimatter annihilation which was expected to yield a signal with a star-like topology made of several pions. Particle detectors and calorimeters surrounded the target to record such annihilation events, and was capable of reconstructing the vertex of the pion-star within the central plane of the 6 12 foil along with the visible energy. In total, the target received ~3 × 10 18 neutrons, with no recorded annihilation events, i.e. with zero background. This was due to an analysis scheme requiring two or more tracks ( ̅-annihilation or background-produced mesons, or their decay products) to be reconstructed in the detector as emanating from the 6 12 foil. As a result, the oscillation limit for free neutrons was established to be → ̅ ≥ 0.86 × 10 8 . (1) In the last two decades since obtaining this result, there have been significant technological developments within the field which have permitted the planning of another transformation experiment, recently proposed at the currently under construction European Spallation Source (ESS) [5,7,8]. According to preliminary estimates, such an experiment could explore this process with 2 − 3 orders of magnitude higher sensitivity than in [6], leading next generation free neutron experiments to be sensitive to oscillation time range → ̅~1 0 9 − 10 10 .
Another way to detect → ̅ is through intranuclear searches, and discovery is tantalizing possible. Searches for → ̅ can be performed in experiments with large underground detectors looking for any hints of the instability of matter. Within the nucleus, spontaneous ̅ production would lead to annihilation with another neighboring nucleon, resulting in the release of ~2 of total energy. However, such intranuclear transformations are significantly suppressed compared to → ̅ in vacuum [5,[9][10][11][12][13]. The limit on the → ̅ intra-nuclear transformation time (in matter) is associated with the square of the free transformation time [5] through a dimensional suppression factor, : In the nucleus, this suppression is due to differences between the neutron and antineutron nuclear potentials; however, in high mass detectors, this suppression can be compensated by the large number of neutrons available for investigation within the large detector volume. A number of nucleon decay search collaborations have been involved in the search for → ̅ in nuclei, such as Frejus [14] and Soudan-2 [15] in 26 56 , and IMB [16], Kamiokande [17], and Super-Kamiokande (SK) [18] in 8 16 ; there has also been a deuteron search performed at SNO [19]. In the Soudan-2 experiment, there is a limit on the transformation time in iron nuclei of ≥ 7.2 × 10 31 [15], which is in line with the limit for the free transformation time of → ̅ ≥ 1.3 × 10 8 . In SK, which extracted 24 → ̅ candidate events while expecting a background count of 24.1 atmospheric neutrino events, these limits were ≥ 1.9 × 10 32 [18] and → ̅ ≥ 2.7 × 10 8 , respectively.
The prevalence of background within SK and other large underground detectors, possibly shrouding a true event, prioritizes the rigorous modeling of both signal and background within an intra-nuclear context. Without any significant improvement in the separation of signal to background in new experiments, it will be possible to improve the appearance limit, but impossible to claim any real discovery. This contrasts the tantalizing figure that future experiments in large underground detectors could improve the restrictions on processes where Δ = ±2 up to ~10 33 − 10 35 [5] in the absence of background. An experiment possibly capable of such a search for → ̅ within the 18 40 nucleus is currently under construction, using large liquid argon ( 18 40 ) time projection chamber: the Deep Underground Neutrino Experiment (DUNE) [20].
Whether or not → ̅ is definitively observed above background in intra-nuclear experiments depends critically upon the separability of signal from background and the energy scale at which the new BSM mechanism will appear. In the case of an observation in intra-nuclear experiments the results will be of great importance for the understanding of fundamental properties of matter, along with building a precise theoretical model describing these properties. Although in the free neutron search [6] no background was detected, the question of background separation might become essential with the planned increase in sensitivity in searches using both free neutrons produced by spallation and bound neutrons in underground experiments, meriting further study beyond this work.
Thus, one requires detailed information about the processes during the annihilation of slow antineutrons on nuclei. The purpose of this work is to create a model describing the annihilation of a slow antineutron incident upon a 6 12 nucleus for the upcoming transformation experiment using a free neutron beam at ESS. Also, the first steps have also been taken towards a full, realistic simulation of the annihilation resulting from → ̅ within 18 40 nuclei for DUNE.

B. Past simulation for free and bound → ̅ searches
In general, the experiment requires maximum efficiency for detection and reconstruction of incredibly rare antineutrons to be separated from background. The development of Monte Carlo (MC) generators for → ̅ searches is not new, and has been an integral part of all past experiments. Sadly, the descriptions of these MCs, as known, are not always complete or seemingly consistent, and are not easily accessible. Information about the generator developed for the ILL experiment [6] is few and far between, unavailable [21], and lacking [22] in detailed explanation.
Intra-nuclear searches have been completed far more times than free neutron experiments, and so their accompanying generators are similarly abundant. Never-the-less, many of their descriptions are scattered throughout a multitude of dissertations and are poorly defended within published works. Similarly, open access to these simulations is lacking. For instance, SK [18] cites only three works in reference to their generator, one of which is a previous work of this paper's lead author, and two of which contain rather ancient antiproton annihilation data; how exactly these are implemented within their model is not available.
The authors are also aware of Hewes' work in relation to → ̅ in DUNE [23]. However, there exist similar issues to those seen in [18], among them the assumption that the annihilation occurs along the density distribution of the nucleus despite the supposed use of work in [13]. In both [18] and [23], only ~10 of exclusive annihilation channels are used, whereas our model utilizes ~100 derived from experimental and theoretical techniques. No previous studies are known to have been tested on their ability to reproduce antinucleon annihilation data, which is a central feature of our work. Our present model is also the first published to incorporate a proper description of the annihilation's dependence on the interaction radius within carbon. Work is underway on a proper implementation of this concept within 18 40 .

C. This work and its goals
Our goal is to create an adequately accurate generator, one which can serve as a platform to be used within all free and intra-nuclear → ̅ experiments. In this article, we present the main framework and approaches underlying the model, wherein the annihilation of an antineutron on the target nucleus is considered to consist of several sequential and independent stages. We use the approach originally undertaken in [24,25].
In the first stage of this approach, one defines the absorption point of an antineutron by the nucleus in the framework of the optical model. Our modeling was performed for 10 antineutrons incident upon a 6 12 nucleus [24,25].
For 18 40 , → ̅ is assumed to occur within the nucleus, where the nucleons have some Fermi motion, and the present paper shows some first steps in this direction; the process of → ̅ within 18 40 will be the focus of our future work. After the point of these quite different initial conditions, all of the following stages of the process for both 6 12 and 18 40 do not differ and are considered within a unified approach.
The second stage in this approach is the actual annihilation of the antineutron with one of the constituent intra-nuclear nucleons. In contrast to [24,25], where a statistical model for nucleonantinucleon annihilation into pions was used, the present paper instead uses a combined approach first proposed in [26] and will be described in Section II D. In this paper we use a version of the annihilation model originating in 1992, utilizing corresponding experimental data available at that time. While there do exist more recent findings from later analyses of LEAR data, these are not many in number and will not greatly affect the conclusions reached for ESS from the model, as these can only slightly modify the probabilities of various annihilation channels within the database of our simulation; these can be updated at a later time when we seek even greater precision for 18 40 .The third stage is the intra-nuclear cascade (INC), initiated by the emergence and nuclear transport of mesons from the annihilation; decays of short-lived resonances are also handled. In this paper, we use the original version of the model which takes into account the nonlinear effect of decreasing the nuclear density, along with a time coordinate [27], which is necessary for the correct description of the passage of resonances through the nucleus.
The final stage is the de-excitation of the residual nucleus.
In this paper, we present a general description of the model and the first results obtained for ̅ in preparation for the forthcoming ESS experiment. In future developments, the description of individual stages of the process can be modified and improved, but the approach remains the same.
The outline of this paper is as follows: Section II provides a detailed description of the model for all successive stages of the process under consideration. In Section III, a comparison will be made between simulation and experiment to test the model against existing at-rest ̅ annihilation data. In Section IV, some validation tests of the ̅ annihilation event generator output data are shown. In Section V, we summarize our work, and briefly consider a future path toward simulation of intra-nuclear transformations in 18 40 .

II. THE MODEL OF THE ANTINEUTRON ANNIHILATION ON THE NUCLEUS A. Absorption of the slow antineutron by the nucleus
In this work, we simulate the annihilation of a cold (~10 ) ̅ on a 6 12 nucleus. The calculation of the total annihilation cross section of an ̅ on 6 12 is a separate problem that is not considered within the scope of this model,, and instead the annihilation event itself is the starting point.
The interaction of a slow antineutron with the nucleus cannot be considered within the framework of the intra-nuclear cascade (INC) model, as is usually done for antinucleon energies above several tens of . Such an interaction also cannot be legitimately modeled using antineutron-nucleon cross sections. The approach used here to describe the interaction between the nucleus and the incoming slow antineutron resulting from the transformation is based on the integration of optical and cascade models. In the optical-cascade model, the initial conditions for the INC are formulated within the optical model. This approach was first applied in [28] to describe the annihilation of stopping antiprotons on nuclei when the antiproton is absorbed from the bound state made by the antiproton orbiting the atom. The same approach was used for the antineutron by L. A. Kondratyuk [24,25] in the discussion of future → ̅ search experiments. The radial ( ) distribution of the absorption probability density ( ) is directly related to the radial nuclear density ( ) and the radial wave function ( ), and is derived from the wave equation for a slow antineutron: This solution for a slow, plane wave antineutron incident on a 6 12 nucleus was presented in great detail in [24,25]. In order to define the annihilation point in simulation, it is desirable to use a simple analytic function. Therefore, we approximate the solution ( ) obtained in [24,25] as a Gaussian function, with a maximum situated at = + 1.2 , where is the radius of half density (with ( 6 12 ) = 2.0403 ) with a width of = 1 . This approximated function is presented in Fig. 1 as the solid orange curve with arbitrary units to demonstrate the penetration depth of the antineutron inside the nucleus.
The model assumes that the proton density within the nucleus ( ) is described as an electrical charge distribution, as obtained in high-energy electron scattering experiments. The function ( ) obeys a Woods-Saxon distribution: where = 0.5227 is the diffuseness parameter of the nucleus, and the radius of half density [29]. For practical reasons within the modeling process, the nucleus is split into seven concentric zones, within which the nucleon density is considered to be constant. Fig. 1 shows the density distribution of the nucleons for 6 12 , calculated by equation (4), along with a step approximation which divides the nucleus into seven zones of constant density. It is seen that although an antineutron penetrates more deeply compared to an antiproton (the dotted line in the Fig. 1), the absorption of the antineutron still occurs about the periphery of the nucleus. Since an antineutron would be strongly absorbed even within the diffuse periphery of the nuclear substance, another eighth zone with density = 0.001 ⋅ (7) is added which extends far beyond the nuclear envelope.  line is a Woods-Saxon density distribution, while the blue step function is an approximation used to divide the nucleus into seven zones of constant density. Right: The radial dependence of the absorption probabilities for the 6 12 are shown for an antineutron (solid orange) and an antiproton (dashed grey) [28]. Note that the eighth zone extends from the end of zone seven at = 4.44 to = 10 .

B. Antineutron annihilation within the nucleus
For intranuclear → ̅ transformation, the conversion of a bound neutron into an associated antineutron is significantly suppressed within the nuclear environment. The reason for this is the large difference between the values of the effective nuclear potential for the neutron and antineutron [5,10]. The question of the magnitude of the suppression of → ̅ occurring within the nucleus has been the subject of in depth nuclear theoretical discussions for a number of years [10][11][12][13][30][31][32]. The transformation would be more probable for neutrons with lower binding energy, and the maximum of the antineutron wavefunction is located beyond the nuclear radius [32]. In contrast to [18,23], for the correct description of the absorption process of the antineutron produced by → ̅ within 18 40 , it is necessary to determine the radial dependence of the probability density of the transformation within the nucleus. This development will be included in our next publication, focused on 18 40 . However, for a first approximation of the annihilation process within 18 40 , the simulation outputs shown in this article are considered for the case when the transformation occurs with equal probability for all neutrons within only the peripheral zone of the nucleus. Thus, we plan to demonstrate the importance of the radial annihilation dependence. With regard to modeling the transformation in the 18 40 nucleus, as discussed in the previous section, the difference between the absorption of the slow antineutron by the 6 12 nucleus comes in the first stage only. As with the 6 12 nucleus, the nucleon density distribution of the 18 40 nucleus is described by expression (4) and is approximated as being divided into seven zones of constant density where it is assumed that the neutrons are distributed throughout the nucleus identically to protons.

C. The nuclear model and nucleon momentum distribution
Within the INC model, the nucleus is considered to be a degenerate, free Fermi gas of nucleons, enclosed within a spherical potential well with a radius equal to the nuclear radius. Nucleons fill all energy levels of the potential well, from the lowest, when a nucleon can have the largest negative potential energy and ~0 momentum, to the highest echelons of the Fermi level, where the nucleon moves with Fermi momentum , and is retained within the nucleus only because of the binding energy (where ≈ 7 per nucleon).
In the interval [0, ], the three-momentum of the nucleon can take all permissible values. The differential probability distribution of the nucleons with respect to the total momentum and kinetic energy [29] takes the form: Here, is the kinetic energy of a nucleon within the nucleus, and ), then their Fermi momentum and energy are easily expressed in terms of the radius. Because every cell in phase space 3 3 contains a number of states ( is the spin of the nucleon) and the total number of protons or neutrons in the nucleus being equal to , it then follows from the normalization condition that and one finally gets that 3 is the volume of the nucleus, and remains the nucleon mass.
If the nucleus is subdivided into concentric spherical zones of constant density, the values of and for each zone are calculated similarly to equations (9) and (10), but with an -th radius, and the density of the nucleons within this -th zone. Fig. 2 shows the spatial distribution of the potential = −( + ) for protons and neutrons in both 6 12 and 18 40 nuclei.

FIG. 2:
The spatial distribution of the potential = −( + ), with appropriate partitioning of the nucleus into seven zones for protons (solid histogram) and neutrons (dotted histogram) for both 6 12 and 18 40 nuclei (for 6 12 , the solid and dotted histograms lay atop one another). is the average nuclear binding energy of 7 per nucleon.
The momentum distribution of the nucleons in individual zones will be the same as for a degenerate Fermi gas, and the probability of a nucleon to have momentum in the th-zone will continue to be determined by (5), although corresponding to th-zone's boundary Fermi momentum value. Fig. 3 shows the momentum distributions of nucleons for both 6 12 and 18 7 nuclei, obtained by summing all the momentum distributions for all individual zones. From Figs. 2 and 3, we can see that the nucleons located in the central zone of the nucleus have the highest value of , and, accordingly, the maximum value of the Fermi momentum . Therefore, the contribution to the total momentum distribution from the nucleons located in the central ( = 1) zone gives the high-momentum part which extends up to 250 − 270 / . Conversely, the nucleons located within the peripheral zone of the nucleus ( = 7) have momenta up to 80 − 100 / . Moreover, the contribution to the overall momentum distribution of a particular zone is greater the more nucleons within it. Thus, in our model there is a correlation of the momentum with the density and, respectively, with the radius. nuclei, summed over all zones. The thin lines show histograms which correspond to contributions from individual zones of the nucleus to the total momentum distribution (only odd-numbered zone distributions are shown so that the picture is not indecipherable). The zone number is shown at the right of each histogram. Note the logarithmic scale of probability in arbitrary units.
Thus, for 6 12 nuclei in the first stage, the nucleons are distributed within the nucleus according to the step density function (see Fig. 1). Next, according to the radial distribution of the antineutron absorption probability in the 6 12 nucleus ( ) (also Fig. 1), the point of annihilation is taken randomly by Monte Carlo technique. The radius of this point determines the number of the zone in which the nucleon partner (neutron or proton) is located, and with which the antineutron annihilates.
The physics underlying the transformation within an 18 40 nucleus is different, since the antineutron generated from → ̅ is not extranuclear in nature, but instead depends entirely on the magnitude of the intra-nuclear binding energy as a function of radius. Generally, it is thought that the bound transformation should occur near the surface of the nucleus (possibly even outside the nuclear envelope [32]), we assume that an antineutron resulting from the transformation has kinetic energy ̅ = + , and annihilates on the nearest nucleon neighbor within only the peripheral zone. Further, the simulation is done within the same scheme for both 6 12 and 18 40 nuclei. The annihilation partner has Fermi momentum randomly selected from the momentum distribution for a particular zone; then, the annihilation occurs.

D. Annihilation model
Unlike papers [24,25], where a statistical model for nucleon-antinucleon annihilation into pions was used, the present paper uses a combined approach first proposed in [26]. The phenomena of ̅ annihilation can lead to the creation of many particles through many possible (at times ~200) exclusive reaction channels; many neutral particles may be present, which can make experimental study quite difficult. Experimental information for exclusive channels is known only for a small fraction of possible annihilation channels, and therefore a statistical model based on (3) symmetry [33] has been chosen to describe the ̅ annihilation. Work to generalize the unitary-symmetric model for ̅ annihilations, along with the development of methods for calculating the characteristics of mesons produced from the annihilation, was performed by I. A. Pshenichnov [34]. According 8 to the model, the ̅ annihilation allows for the production of between two and six intermediate particles. Given the estimates of the phase space volume at low momenta, the production of a larger number of intermediate particles is unlikely. Intermediate particles, such as , , and mesons, are all possible; the channels with strangeness production are not considered within this version of the model. This unitary-symmetric statistical model predicts 106 ̅ annihilation channels, and 88 ̅ annihilation channels, but this differs from experiments, which effectively measure only ~40 channels for ̅ and ~10 channels for ̅ annihilation. However, neither the statistical model, nor the experimental data, can provide a complete and exclusive description of the elementary nucleon-antinucleon annihilation processes. For this reason, semiempirical Tables I and II of annihilation channels are employed for use in annihilation modeling. These are obtained as follows: First, all experimentally measured channels were included in Tables I and II. Then, by using isotopic relations, probabilities were found for those channels that have the same configurations but different particle charges. Finally, the predictions of the statistical model with (3) symmetry were entered for the remaining intermediate channels.
Sometimes the probabilities of intermediate channels measured in different experiments differ significantly. In this case, the data in the semi-empirical tables were corrected within experimental accuracies in order to describe the topological cross section for ̅ and ̅ in a consistent way. In our approach, a substantively large collection of experimental data was used: multi-particle topologies, inclusive spectra, topological pion cross sections, and branching ratios of various resonance channels. The following pages show the semiempirical tables, with probabilities of various ̅ and ̅ annihilation channels included. In further modeling of ̅ interactions, it is considered that channels for ̅ are identical to ̅ channels, and that annihilation channels for ̅ are charge conjugated to ̅ channels.
Considering the laws of energy and momentum conservation for each annihilation, the procedure for simulating the characteristics of both the intermediate particles and their various decay products consists of the following: first, a single channel from the table is randomly selected via Monte Carlo technique as the initial state, with all necessary momenta of all annihilation particles determined according to the pertinent phasespace volume. This takes into account the Breit-Wigner mass distribution for meson resonances, while all pions have a mass value of 0.14 2 .
The subsequent disintegration of unstable mesons is modeled according to experimentally known branching ratios. All major decay modes for meson resonances have been considered, such as in Table III. All experimental data used for comparison with this annihilation model are described in great detail in [26]. However, there do exist more recent data obtained from LEAR by the Crystal Barrel [35] and OBELIX [36] Collaborations on some exclusive channels which show somewhat different branching ratios from those used by us. We plan to make a revision of the annihilation tables in the near future, taking into account all data. Below is a comparison of the simulation results and experimental data on ̅ annihilation at rest.  Fig. 4 shows the pion multiplicity distribution generated by ̅ annihilation, while Fig. 5 shows the charged pion momentum distribution. From considering Table  IV and Figs. 4 and 5, it follows that the Monte Carlo and available experimental data are in general agreement with the main features of ̅ annihilation. In all, our annihilation model utilizes a complex series of tables with a much larger number of predicted and pertinent channels than [18,23]. As this approach demonstrates a good description of the experimental data for ̅ annihilation at rest, we believe that it is also adequate for an accurate description of ̅ annihilation at rest, and so can be implemented within our ̅ annihilation simulation.  [26]. Note that 2) indicates that the probabilities are obtained from isotopic relations. The sum of all branching ratios is normalized to 100 percent.  . Similarly note that 1) indicates a probability attained from experiment; see references used in [26]. Note that 2) indicates that the probabilities are obtained from isotopic relations. The sum of all branching ratios is normalized to 100 percent.     [41].

E. The Intra-nuclear Cascade (INC) Model
Inelastic nuclear interactions are clearly statistical in nature, as they can be realized in many possible states. A statistical approach is key to describing such systems, and replaces the evolution of a system's wave function with the description of the evolution of an ensemble of the many possible states of the system. There are two dramatically different stages of a deeply inelastic interaction: 1) a fast, out-of-equilibrium stage in which energy is redistributed between the various degrees of freedom within the nucleus as a finite open system, and 2) the slow equilibrium stage of the decay of the thermalized residual nuclei.
The INC model is a phenomenological model describing the out-of-equilibrium stage of inelastic interactions and operates with the notion of the probability of a nuclear system being in a given state. Transitions between different states are caused by two-body interactions, which lead to secondary particles exiting the nucleus, dissipating the excitation energy in the process. However, this phenomenological model is linked to fundamental microscopic theory. It was shown in [42] that it is possible to transform a nonstationary Schrödinger equation for a many body system into kinetic equations, if large energy (and so short time) wave packet formulations are used.
To explain, if the duration of the wave packet's individual collisions are shorter than the interval of time between consecutive collisions, then the amplitudes of these collisions will not interfere. This condition is essentially analogous to the condition of a free gas approximation: 0 < , where 0 is the duration time of the collision, and is the mean-free-path time. This condition allows for the consideration of a particle's motion as in a dilute gas with independent particle motion on free path trajectories perturbed by binary collisions. Under these conditions, in a quasi-classical way, one can use the local momentum approximation by assigning a particle a momentum ⃗ ( ) between consecutive collisions. In this case, the quantum kinetic equation is transformed into a kinetic equation of Boltzmann-type describing the transport of particles within nuclear media; this differs from the conventional Boltzmann equation only by accounting for the Pauli exclusion principle. Thus, the INC model is a numerical solution of the quasi-classical kinetic equation of motion for a multi-particle distribution function using the Monte Carlo method.
We will now focus our discussion on the scope of the INC model and the possibility of generalizing its use, such as in the event of the absorption of a slow antineutron. The principles underlying the model are altogether justified if the following conditions are met [29,42,43]: a) The wavelengths, , of the majority of moving particles should be less than the mean distance between nucleons within the nucleus, i.e. < Δ, where In this case, the system acquires quasiclassical characteristics, and one can speak of the trajectories of particles and two-body interactions within the nucleus. For individual nucleons, this corresponds to an energy of more than tens of . Of course, this condition cannot be met in the case of a slow antineutron, and therefore, its absorption is described in the framework of the optical model. b) The interaction time should be less than the time between successive interactions ≤ , where ≈ ≈ 10 −23 , and is the nucleon radius. The mean-free-path-length time is ( is the cross section in ). This requirement is equivalent to the condition of requiring sufficiently small cross sections of elementary interactions and proves problematic for pions produced from the annihilation and lying within the energy range of the Δ-resonance, where > 100 . However, it should be kept in mind that the effective mean-free-path-length within the nucleus is increased by the Pauli exclusion principle; secondarily, because the uptake of the antineutron is predominately on the periphery of the nucleus, where the nuclear density is low and the distance between the nucleons large, one can expect that the INC model would work in this case. Never-the-less, the comparison of the simulation results with experimental data is the main criterion for the applicability of the model. The standard INC model is based on a numerical solution of the kinetic equation using a linearized approximation, which implies that the density of the media does not change in the development of the cascade, i.e. ≪ (where is the number of cascading particles, and is the number of nucleons making up the target nucleus). Such an approximation is violated in the case of multi-pion production in and interactions at , ≥ 3 − 5 , and also in the case of annihilation, especially when considering light nuclei such as 6 12 . A version of the model, which takes into account the effect of a local reduction in nuclear density, was first proposed in [43]. This version of the model considers the nucleus as consisting of separate nucleons, the position of their centers computed by Monte Carlo method according to the prescribed density distribution ( ) such that the distance between their centers is no less than 2 , where = 0.2 is the nucleon core radius. A cascading particle may interact with any intra-nuclear nucleon which lies inside the cylinder of diameter 2 + extending along the particle's velocity vector (here, is the interaction radius, while is the deBroglie wavelength of the particle). The is a parameter of the model and is chosen for better agreement with the experimental data. The key point to understand here is the ability to determine the probability of the cascading particle interacting with another constituent nucleon. We now consider this process in more detail.
Within the standard cascade model, the randomly chosen interaction point is computed from a Poisson distribution for the mean-free-pathlength. In this case, the probability ( ) of the particle experiencing collisions along the pathlength in media with density , where the particle has a total cross section , is defined as: If on the path-length there lie individual particle centers, each has an equal collision probability for the particle to collide on of centers and = 1 − . This probability is described by a binomial distribution: From the Poisson distribution (11), it follows directly that the probability of a particle experiencing no collisions along is simply: The same probability for this process can be obtained from the binomial distribution in (20): If one takes (0) = (0, , ), and when considering that = ( + /2) 2 , then: An essential feature of present version of the INC model is the fact that after interactions occur inside the nucleus, the nucleon is considered to be cascade particle and not a constituent part of the nuclear system. Thus, a reduction in nuclear density takes place during the cascade development. In order to describe the evolution of the cascade and the decays of unstable meson resonances over time, an explicit time-coordinate has been incorporated into the model.
So, one can summarize the physical considerations that underlie the INC model as follows: • The nuclear target is a degenerate Fermi gas of protons and neutrons within a spherical potential well with a diffuse nuclear boundary. The real nuclear potentials for nucleons ( ), antinucleons ( ⃐⃗ ⃗ ), and mesons ( , , ) effectively takes into account the influence on the particle of all intranuclear nucleons. The depth of the potential well for the antinucleon and mesons within the nucleus remains a free-parameter of the model. Recognizing that the annihilation process usually occurs on the periphery of the nucleus, a good approximation for this is considered to be ⃐⃗ ⃗ ≈ 0 and , , ≈ 0. In the future, a detailed study is planned to focus on the influence of these potentials on the simulation output.
• Hadrons involved in collisions are treated as classical particles. A hadron can initiate a cascade of consecutive, independent collisions upon nucleons within the target nucleus. The interactions between cascading particles are not taken into account. • The cross sections of hadron-nucleon interactions are considered within the nucleus to be identical to those in vacuum, except that Pauli's exclusion principle explicitly prohibits transitions of cascade nucleons into states already occupied by other nucleons.
Elementary processes, such as those seen in the channels (16) shown above, are described by empirical approximations from analysis of experimental data on and interactions at kinetic energies < 20 [29,43] Now consider some of the features of the INC model related to the introduction of unstable meson resonances into the model. Modeling annihilation with meson resonances (i.e. ̅ → + + + ) was described in the preceding section. It is assumed within the model that -mesons produced by annihilation decay quickly enough to avoid interacting with any intra-nuclear nucleons. In contrast, -mesons produced by annihilation can both interact with other intra-nuclear nucleons and decay within or outside the nucleus. The competition between the decay of the -meson and its interaction with intra-nuclear nucleons is determined by the expression for the mean-free-path: where = ( ) −1 , = (ℎ ) −1 , is the nuclear density, and is the Lorentz factor. The mean lifetime of the -meson is large enough for the particle to be considered stable within the nucleus, which can then decay upon exit. The model uses the experimentally measured decay modes of the meson resonances described above. When the annihilation products are allowed to disintegrate, their three-body decay is simulated by evaluation of the permissible phase-space volume.
To accommodate the passage of -and--mesons through nuclear material, in addition to channels listed in (16), the following interactions are also considered: Along with the creation of -and -mesons by annihilation, the model also accounts for the creation of mesons through interactions between annihilation pions and nucleons, such as For cross sections of reactions in (18), estimates given in [26] were employed. For those few reactions shown in (19), experimental cross sections were taken from compilation [44]. As these interactions are considered at relatively low energy, the angular distributions for reactions shown in (18) and (19) are assumed to be isotropic in the center of mass of the system. Reactions with three particles in the final state are simulated via their pertinent phase-space volume.

F. De-excitation of the residual nucleus
For inelastic nuclear reactions, after the rapid stage of the intra-nuclear cascade ( ≃ 0 ) and once statistical equilibrium ( ≅ (5 − 10) 0 ) is established inside the residual nucleus, a slow stage begins ( ≫ 0 ) involving the disintegration of the highly excited residual nucleus (note that 0 ≤ 10 −22 , which is the average time required for a particle to pass completely through the nucleus). The INC model is able to describe the dissipation of energy throughout the nucleus. At the end of the cascade stage, the nuclear degenerate Fermi gas contains a number of "holes" ℎ , which is equal to the number of collisions of cascade particles with nucleons within the nucleus. Also, there exists some number of excited particles , which is equal to the number of slow cascade nucleons trapped by the nuclear potential well. The excitation energy of the residual nucleus * , is the sum of the energy of all such quasiparticles calculated from the Fermi energies : The resulting residual nuclei have a broad distribution on the excitation energies * , momenta, masses, and charges. The INC model correctly accounts for the fluctuations of the cascade particles, and reliably defines the entire set of characteristics for residual nuclei.
The de-excitation mechanism for a residual nucleus is determined from the accumulated excitation energy of the nucleus [45]. Under low excitation energies (where * ≤ 2 − 3 ), the primary de-excitation mechanism is the consecutive emission (evaporation) of particles from the compound nucleus [46]. When the excitation energy of the nucleus is approaching the total binding energy (where * ≥ 5 . ), the prevalent mechanism is explosive decay [47]. For intermediate energies, both mechanisms coexist.

III. COMPARISON WITH EXPERIMENT
The optical-cascade model described throughout this work has been used to analyze experimental data taken from antiproton annihilation at rest on 6 12 target nuclei. Table V shows the average multiplicity of emitted pions and protons. Experimental data on average pion multiplicities (values of which are shown at the bottom of the ̅ row) are taken from [38]. The final column of the table indicates the average energy of pions and photons (resulting from the decay of -and -mesons) emitted from the nucleus. Calculated values for the average multiplicities of pions (values of which are shown at the top of the ̅ row) are within accuracies of the experimental data. Since the antiproton primarily annihilates on the surface of the nucleus, most of the mesons produced fly out of the nucleus without any interaction. In the case of a light nucleus such as 6 12 , the effect of absorption of annihilation mesons is not large and the average multiplicity of pions emitted appear to be quite similar to the multiplicity of pions in ̅ annihilation (4.910). For comparison, Table V also shows results which simulate the annihilation of a slow antineutron on a 6 12 nucleus. The comparison shows that the average pion multiplicity for an ̅ annihilation is somewhat lower, and that the average multiplicity for exiting nucleons slightly higher than the case of a stopped antiproton. This is due to the fact that the antineutron penetrates more deeply into the nucleus (seen in the solid line shown in Fig. 1) compared to an antiproton (seen in the dashed line shown in Fig. 1), and so there are more intra-nuclear interactions between annihilation mesons and constituent nucleons. Thus, the number of mesons emitted from the nucleus and their total energy are reduced, while instead the number of nucleons that were kicked from the original nucleus during the fast cascading stage (and then emitted from the nucleus during the de-excitation process) is increased. In the case of peripheral annihilation of an antineutron on 18 40 , the pions are almost entirely free to leave the nucleus, increasing the value of , and so the number of emitted nucleons is significantly lower than ̅ annihilation. Note here that the calculation completed for 18 40 is made in a very rough approximation with respect to the annihilation radius; this property requires further detailed investigation.  Fig. 6 shows the charged pion multiplicity distribution emitted from the nucleus due to ̅ annihilation (shown as the solid histogram with points), ̅ (shown as the dotted histogram), and ̅ (the dashed histogram). As was expected, the differences in these distributions, as with the mean number of emitted pions, are not significant, although there appears to be some bias towards a smaller number of pions for ̅ and a larger number for ̅ . Fig. 7 shows the distribution by number of events with the charge carried out by pions. For the ̅ annihilation the maxima of the distribution are = −1 and = 0, which practically corresponds to mesons exiting the nucleus without any interaction with nucleons. The optical-cascade model demonstrates good agreement with the experimental data. In the case of an annihilation with an ̅, the distribution has a maximum that is shifted to = 0 and = +1, respectively. In the case of a peripheral annihilation for ̅ , the distribution has a narrower maximum than ̅ .  [49], open squares- [50]. The dotted histogram shows an ̅ simulation. The dashed histogram shows an ̅ simulation. Fig. 9 shows the momentum distribution for + exiting the nucleus, which is rather similar to the momentum distribution of pions created by ̅ annihilation (as seen in Fig. 5). To understand the uncertainty of the model, calculations were done 1) without any nuclear potential for the antineutron, and, as an option, 2) with a model where the antineutron nuclear potential is introduced similarly to [26]. For mesons propagating inside the nucleus, we have not assumed any nuclear potentials. Both model calculations are presented in Fig. 9, and show rather good agreement with experimental data, although there is some exaggerated absorption behavior corresponding to the Δ-resonance region (~260 ). The difference between experimental measurements appears to be of the same order as the uncertainty in the calculation. Never-the-less, in the near future, we plan to study this question in more detail.
FIG. 9. The + momentum distribution is shown for antiproton annihilation at rest on 6 12 nuclei.
The points show experimental data from [38,52]. The histograms show calculations, where the solid line shows an option of the model with a nuclear potential for the antineutron, and the dashed line is the calculation done without this potential. The nuclear potential for annihilation mesons inside the nucleus is assumed to be zero. Fig. 10 shows the energy spectrum of protons exiting the nucleus from ̅ annihilation at rest. In the low energy regime (up to 50 ), evaporative protons provide a significant contribution to the spectrum. The model again shows good agreement with the available experimental data.
FIG. 10. The exiting proton kinetic energy spectrum due to antiproton annihilation at rest on 6 12 nuclei. The solid histogram shows the simulation result. The dotted histogram shows the contribution which evaporative protons impart to the whole distribution. The points show the experimental data taken in [53].
From the comparisons made above between the simulation results of the optical-cascade model and experimental data of ̅ annihilation at rest, it follows that the model as a whole describes experiments well, thus accurately reflecting the dynamics of the annihilation process and the propagation of annihilation mesons throughout the nucleus.

IV. ̅ ANNIHILATION GENERATOR VALIDATION
Colleagues within the ESS NNBar Collaboration have tested the model through its corresponding event generator comprehensively. The event generator outputs an annihilation point within 6 12 , annihilation product particle identities, energies, momenta, etc. These particles and variables are tracked as outputs and saved to file in three successive stages: 1) after the primary annihilation, 2) after all cascading ( , , , , , ) and evaporation particles have left the nucleus, and 3) after all decays of the meson resonances emitted from the nucleus ( , , ) are modeled.
As discussed in the previous sections, multiplicity, charge, momentum and energy distributions of particles show good agreement with antiproton annihilation experimental data, and all simulated variables quantitatively satisfy the fundamental tenets of the required physics. Specifically, the generator has been shown to conserve charge, energy, momentum, baryon number, etc., through all three stages of simulation. The output file type is .txt, and formatted in such a way as to easily separate the particle content and their respective physical variables through the stages. Analysis of the output has been completed by ESS colleagues using C++ and the CERN ROOT 5.34 scientific software framework [54]. 100,000 simulated ̅ annihilation events are available upon request from the authors. The currently available event files and the following plots are created from the completed simulation file data without either antineutron or meson potentials.
An important characteristic for any relativistic many particle system is the invariant mass. One may analyze the invariant mass distribution for annihilation mesons at the annihilation point, and then see how it distorts due to interactions throughout the nucleus. Detector performance might affect the invariant mass further, but this study in not the focus of this paper. Fig. 11 shows how the distribution of invariant mass changes for all outgoing pions and photons generated by ̅ annihilation products (solid), a result of interactions with nuclear media. The dotted line shows the original distribution of invariant mass of the initial ̅ annihilation products. Fig. 11 shows that the intra-nuclear interactions of annihilation mesons with nucleons have resulted in a significant redistribution of energy between mesons and other nuclear constituents, shifting and smearing the initial distribution of down to values of ~1.2 2 ⁄ . Note that the higher the initial value of , or the deeper the penetration of the antineutron into the nucleus, the larger the number of mesons which will interact with the nuclear environment, quickly devouring this particular part of the distribution.
Similarly, for Fig. 12, we see that the momentum distribution reconstructed from initial annihilation mesons is perturbed and expanded by transport through the nuclear environment. The structure shown in the dotted histogram shows a similar distribution as seen in Fig. 3, though implicitly convolved with Fig. 1, and considerate of different scales. After transport, this distribution distorts as particles cascade through the nucleus, shifting values up to as far as ~0. 8 . Following discussion of Fig. 2 within [18], a highly relevant plot of total momentum versus invariant mass output variables is shown below in Fig. 13 for outgoing pions and photons. The projection of the -axis is precisely Fig. 11, while the -axis is Fig. 12. Across Figs. 11-13, all bin widths are identical (10 2 ⁄ or ⁄ ), and all counting scales are logarithmic. Note the bright spot at ~1.9 2 in invariant mass and ~0.1 in total momentum; this shape curves slightly upward and to the right, and contains the ~35% of all exiting pions and photons which go through the nucleus without interaction.
FIG. 13: Stage 3 total momentum vs. invariant mass of pions and photons. Note the double lobe structure; this is due to the absorption of a single pion during transport. Also recognize the thin, sickle-like shape in the lower right-hand corner of the figure; this is due to invariant mass of the initial Stage 1 mesons which did not interact with the nucleus.

V. CONCLUSION
This work has endeavored to demonstrate the detail of the optical-cascade model for describing antineutron annihilation on 6 12 nuclei. It is quite important that the absorption of a slow antineutron is described within the framework of the optical model and that radial dependence of the annihilation probability is used within the initial stage of the simulation. A combination of experimental data with the results of a statistical model employing (3) symmetry is used to describe the annihilation process. The propagation of annihilation-produced pions and heavier meson resonances within the nucleus is described by the intra-nuclear cascade model, which takes into account the nonlinear effect of decreasing nuclear density. The process of deexcitation of residual nuclei is described by a combination of the evaporation model and the Fermi model of explosive disintegration. This combined approach shows good agreement with experimental data in the modeling of antiproton annihilation at rest on 6 12 nuclei, and provides a reliable predictor of the characteristics of slow antineutron annihilation on 6 12 . This model can thus be used as an event generator in the design of detector systems for planned experimental searches for neutron-antineutron transformation at the European Spallation Source, employing a free beam of cold neutrons.
This approach is universal, and can be used for simulating antineutron annihilation on many different nuclei. However, to search for the transformation occurring within nuclei (for example, within 18 40 , with no external source of antineutrons), a valid model can only be created when a proper definition of radial annihilation probability density is incorporated, allowing for the derivation of intra-nuclear → ̅ transformation constraints. The model proposed in this work can be thought of as a first step in the preliminary modeling of this full process.
Expressing an aside into future developments, it is planned that more precise modeling will be rendered for the intra-nuclear → ̅ transformation thought possibly to take place with 18 40 . Proper simulation of such a signal will be important to assess the feasibility of significantly improved transformation searches in the DUNE experiment. The goal of accurate and precise simulation is sought for the absolute suppression of the atmospheric neutrino background in the DUNE experiment, which will allow for an improvement in the search limits by several orders of magnitude. Inevitably, this will also require careful study and accurate simulation of atmospheric neutrino events in DUNE. The basic physics, model, file types, and analysis techniques presented in this article will continue to be employed, though special care must be taken in the integration of proper intra-nuclear transformation and radial annihilation probability distributions into the simulation. One other → ̅ generator, already developed internally to the DUNE collaboration [23], is currently being studied by multiple colleagues in both the DUNE and ESS collaborations; complex techniques have been developed using neural networks and multivariate boosted decision trees to study the separability of supposed signal from atmospheric neutrino background with promising results. However, when considering the subtleties of the simulation assumptions and techniques, some room for improvements within the generator are thought to be possible. Thus, the independent generator development described in this work, along with its future 18 40 extension, will help in understanding the potential and limits of exploration of rare processes like → ̅ where separation from background plays a major role. Altogether, this will hopefully lead to a fruitful collaboration and collective meta-analysis between generators and groups, which could reputably assess the probability of separating signal from background sources in such a large, underground experiment.