Progress and Simulations for Intranuclear Neutron--Antineutron Transformations in ${}^{40}_{18} Ar$

TL;DR

This paper presents a novel intranuclear simulation of neutron-antineutron transformations in argon nuclei, crucial for rare decay searches in upcoming large-scale neutrino experiments, with new theoretical calculations and simulation comparisons.

Contribution

It introduces a self-consistent intranuclear simulation for neutron-antineutron transformations in argon, including the first calculation of the suppression factor and comparisons with existing models.

Findings

New simulation model for intranuclear $n ightarrow\bar{n}$ processes in argon.

First calculation of argon's $ar{n}$-intranuclear suppression factor.

Comparison with existing Monte Carlo generator models.

Abstract

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 $\mathcal{B}$ and $\mathcal{B}-\mathcal{L}$ violating processes such as neutron--antineutron transformations ($n\rightarrow\bar{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}_{18} Ar$, 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, including the first calculation of ${}^{40}_{18} Ar$'s $\bar{n}$-intranuclear suppression factor, and Monte Carlo simulation comparisons to another publicly available $n\rightarrow\bar{n}$ generator within GENIE, are also discussed.

Tables

Table 1. Table 1: A list of several 𝒑 ¯ ​ 𝒑 bold-¯ 𝒑 𝒑 \bar{p}p annihilation products and their respective old ( Golubeva:2018mrz, ) and new branching fractions with inclusion of data from ( Salvini:2004gz, ; Amsler:2003bq, ) .

Channel	Probability (%) from (Golubeva:2018mrz, )	Probability (%) Used for This Work
$\overline{𝒑}𝒑\to {𝝅}^{+}{𝝅}^{+}{𝝅}^{-}{𝝅}^{-}$	$2.74$	$3.64$ (Salvini:2004gz, ) 555OBELIX gives data for $\overline{p}p\to {\pi }^{+}{\pi }^{+}{\pi }^{-}{\pi }^{-}$ with a branching ratio of $BR\left(\overline{p}p\to {\pi }^{+}{\pi }^{+}{\pi }^{-}{\pi }^{-}\right)=6.4\pm 0.09\%$ (a gas target). This value includes both resonant and independent production of these four charged $\pi$’s. To avoid double counting, we have subtracted from this value all the fractions of all of the channels that give rise to ${\pi }^{+}{\pi }^{+}{\pi }^{-}{\pi }^{-}$ in the final state. Thus, in the $\overline{p}p$ annihilation table contained in (Golubeva:2018mrz, ), we have experimental data showing $BR\left(\overline{p}p\to {\rho }^0{\rho }^0\right)=0.67\%$, $BR\left(\overline{p}p\to {\pi }^{+}{\pi }^{-}{\rho }^0\right)=2.02\%$, and importantly $BR\left(\overline{p}p\to {\pi }^{+}{\pi }^{-}\omega \right)=3.03\%$, whose contribution to the full channel of $\overline{p}p\to {\pi }^{+}{\pi }^{+}{\pi }^{-}{\pi }^{-}$ (including $BR\left(\omega \to {\pi }^{+}{\pi }^{-}\right)=2.3\%$ is $3.03\times 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\%$.
$\overline{𝒑}𝒑\to {𝝅}^{\mathrm{𝟎}}{𝝅}^{\mathrm{𝟎}}$	$0.02$	$0.154$ (Amsler:2003bq, )
$\overline{𝒑}𝒑\to 𝜼𝜼$	$0.01$	$0.0312$ (Amsler:2003bq, )
$\overline{𝒑}𝒑\to {𝝅}^{\mathrm{𝟎}}𝝎$	$0.58$	$0.460$ (Amsler:2003bq, )
$\overline{𝒑}𝒑\to 𝜼𝝎$	$0.34$	$0.960$ (Amsler:2003bq, )
$\overline{𝒑}𝒑\to {𝝅}^{\mathrm{𝟎}}{𝝅}^{\mathrm{𝟎}}{𝝅}^{\mathrm{𝟎}}$	$1.12$	$0.610$ (Amsler:2003bq, )
$\overline{𝒑}𝒑\to {𝝅}^{\mathrm{𝟎}}{𝝅}^{\mathrm{𝟎}}𝜼$	$0.54$	$0.514$ (Amsler:2003bq, )
$\overline{𝒑}𝒑\to 𝝎𝝎$	$1.57$	$0.358$ (Amsler:2003bq, )
$\overline{𝒑}𝒑\to {𝝅}^{\mathrm{𝟎}}{𝝅}^{\mathrm{𝟎}}𝝎$	$0.79$	$1.53$ (Amsler:2003bq, )
$\overline{𝒑}𝒑\to 𝜼{𝝅}^{\mathrm{𝟎}}𝝎$	$0.30$	$0.60$ (Amsler:2003bq, )
$\overline{𝒑}𝒑\to {𝝅}^{\mathrm{𝟎}}𝝎𝝎$	$0.37$	$0.344$ (Amsler:2003bq, )
$\overline{𝒑}𝒑\to {𝝅}^{\mathrm{𝟎}}{𝝅}^{\mathrm{𝟎}}{𝝅}^{\mathrm{𝟎}}𝝎$	$0.40$	$1.24$ (Amsler:2003bq, )

Table 2. Table 2: Meson multiplicity comparisons from previous work ( Golubeva:2018mrz, ) and the new generator as described here, which include experimental integration of more recent OBELIX ( Salvini:2004gz, ) and Crystal Barrel ( Amsler:2003bq, ) data sets.

	$\overline{𝒑}𝒑$ Simulation from (Golubeva:2018mrz, )	$\overline{𝒑}𝒑$ Simulation with New Model	$\overline{𝒑}𝒑$ Experiment
$𝑴\left(𝝅\right)$	$4.91$	$4.95$	$4.98\pm 0.35$ (Klempt:2005pp, ), $4.94\pm 0.14$ (Minor, )
$𝑴\left({𝝅}^{\pm }\right)$	$3.11$	$3.09$	$3.14\pm 0.28$ (Klempt:2005pp, ), $3.05\pm 0.04$ (Klempt:2005pp, ), $3.04\pm 0.08$ (Minor, )
$𝑴\left({𝝅}^{\mathrm{𝟎}}\right)$	$1.80$	$1.86$	$1.83\pm 0.21$ (Klempt:2005pp, ), $1.93\pm 0.12$ (Klempt:2005pp, ), $1.90\pm 0.12$ (Minor, )
$𝑴\left(𝜼\right)$	$0.09$	$0.09$	$0.10\pm 0.09$ (Levman:1979gg, ), $0.07\pm 0.01$ (Klempt:2005pp, )
$𝑴\left(𝝎\right)$	$0.20$	$0.27$	$0.28\pm 0.16$ (Levman:1979gg, ), $0.22\pm 0.01$ (Hamatsu:1976qz, )
$𝑴\left({𝝆}^{+}\right)$	$0.19$	$0.19$	—–
$𝑴\left({𝝆}^{-}\right)$	$0.19$	$0.18$	—–
$𝑴\left({𝝆}^{\mathrm{𝟎}}\right)$	$0.19$	$0.18$	$0.26\pm 0.01$ (Hamatsu:1976qz, )

Table 3. Table 3: A list of updated multiplicities from experimental data, our original work concerning \isotope [12][6]C ( Golubeva:2018mrz, ) , and two new calculations taking into account the newest versions of 𝒑 ¯ bold-¯ 𝒑 \bar{p} annihilation branching ratios while also considering a new intranuclear antinucleon potential with an associated nuclear medium response.

	$𝑴\left(𝝅\right)$	$𝑴\left({𝝅}^{+}\right)$	$𝑴\left({𝝅}^{-}\right)$	$𝑴\left({𝝅}^{\mathrm{𝟎}}\right)$	${𝑬}_{𝒕𝒐𝒕}$ (MeV)	$𝑴\left(𝒑\right)$	$𝑴\left(𝒏\right)$
$\overline{𝒑}𝐂$ Experiment	$4.57\pm 0.15$	$1.25\pm 0.06$	$1.59\pm 0.09$	$1.73\pm 0.10$	$\mathrm{𝟏𝟕𝟓𝟖}\pm \mathrm{𝟓𝟗}$	—–	—–
$\overline{𝒑}𝐂$ Old Calculation	$4.56$	$1.21$	$1.63$	$1.72$	$\mathrm{𝟏𝟕𝟑𝟔}$	$1.14$	$1.21$
$\overline{𝒑}𝐂$ Calculation $\#\mathrm{𝟏}$	$4.56$	$1.21$	$1.63$	$1.72$	$\mathrm{𝟏𝟕𝟑𝟖}$	$1.11$	$1.20$
$\overline{𝒑}𝐂$ Calculation $\#\mathrm{𝟐}$	$4.60$	$1.22$	$1.65$	$1.73$	$\mathrm{𝟏𝟕𝟔𝟐}$	$0.96$	$1.03$

Table 4. Table 4: Final state average stable particle multiplicities for several 10,000 event samples across multiple 𝒏 ¯ ​ 𝑨 ​ 𝒓 bold-¯ 𝒏 𝑨 𝒓 \bar{n}Ar annihilation MCs. Also included is the initial total annihilation energy. See ( Golubeva:2018mrz, ) for a full description of the zoned local Fermi gas and the intranuclear cascade used in this work. See the GENIE v3.0 manual ( Andreopoulos:2015wxa, ) for discussions of nuclear models and intranuclear cascades.

	$𝑴\left(𝝅\right)$	$𝑴\left({𝝅}^{+}\right)$	$𝑴\left({𝝅}^{-}\right)$	$𝑴\left({𝝅}^{\mathrm{𝟎}}\right)$	$𝑴\left(𝒑\right)$	$𝑴\left(𝒏\right)$	${𝑬}_{𝒐}^{𝒕𝒐𝒕}$ (MeV)
$\overline{𝒏}{}_{\mathrm{𝟏𝟖}}{}^{\mathrm{𝟑𝟗}}\mathrm{𝐀𝐫}$ Golubeva-Richard-Paryev (Zoned Local Fermi Gas, INC)	$3.813$	$1.239$	$1.008$	$1.566$	$3.459$	$4.823$	$1.846$
$\overline{𝒏}{}_{\mathrm{𝟏𝟖}}{}^{\mathrm{𝟑𝟗}}\mathrm{𝐀𝐫}$ Golubeva-Woods-Saxon-Paryev (Zoned Local Fermi Gas, INC)	$3.781$	$1.22$	$0.998$	$1.563$	$3.63$	$4.896$	$1.845$
$\overline{𝒏}{}_{\mathrm{𝟏𝟖}}{}^{\mathrm{𝟑𝟗}}\mathrm{𝐀𝐫}$ GENIEv3.0.6 (Default Bodek-Ritchie, hA2018)	$3.610$	$1.183$	$0.991$	$1.436$	$3.021$	$3.151$	$1.925$
$\overline{𝒏}{}_{\mathrm{𝟏𝟖}}{}^{\mathrm{𝟑𝟗}}\mathrm{𝐀𝐫}$ GENIEv3.0.6 (Default Bodek-Ritchie, hN2018)	$3.280$	$1.159$	$0.968$	$1.153$	$6.192$	$6.654$	$1.922$
$\overline{𝒏}{}_{\mathrm{𝟏𝟖}}{}^{\mathrm{𝟑𝟗}}\mathrm{𝐀𝐫}$ GENIEv3.0.6 (Local Fermi Gas, hA2018)	$3.594$	$1.173$	$0.9776$	$1.444$	$3.045$	$3.174$	$1.908$
$\overline{𝒏}{}_{\mathrm{𝟏𝟖}}{}^{\mathrm{𝟑𝟗}}\mathrm{𝐀𝐫}$ GENIEv3.0.6 (Local Fermi Gas, hN2018)	$3.24$	$1.155$	$0.956$	$1.129$	$6.269$	$6.718$	$1.905$

Equations

$$\begin{array}{c}{\Psi }_n=\frac{1}{\sqrt{4\pi }}\frac{u_n(r)}{r},{\int }_0^{\infty }{u}_n(r{)}^2\mathrm{d}r=1,\\ u_n(r)=N_n\left[\exp (-{\lambda }_{1}r)-\exp (-{\lambda }_{2}r)\right],\end{array}$$

$$\frac{1}{\sqrt{4\pi }}\left[\frac{u_n(r)}{r}|pn⟩+\frac{w(r)}{r}|p\bar{n}⟩\right].$$

$$-\frac{w^{\prime \prime }(r)}{m}+Vw(r)-Ew(r)=-\gamma u_n(r).$$

$$\Gamma =$$

$$=$$

$$T=T_R{\tau }^2,$$

$$-{u_n}^{\prime \prime }/\mu +V_n(r)u_n(r)=E{u}_n(r),$$

$${\tilde{u}}_n(x)=\sum _{j=1}^N{a}_j\sin (j\pi x)$$

$$(\bar{A}-E\text{mathbbm}1)W=\gamma U_n,$$

$$\Gamma \approx -2/\bar{W},\bar{W}={\int }_0^{\infty }|u_n(r){|}^{2}W(r)\mathrm{d}r.$$

$$U_{\bar{n}}(r)=V_0\frac{\rho (r)}{{\rho }_0},$$

$$m_{\bar{n}}^{\ast }(r)=m_{\bar{n}}+U_{\bar{n}}(r)$$

$$E_{\bar{n}}^{\prime }=\sqrt{m_{\bar{n}}^{\ast 2}+{\mathbf{p}}_{\bar{n}}^{\prime 2}}.$$

$$E_{\bar{n}}^{\prime }=\sqrt{m_{\bar{n}}^{\ast 2}+{\mathbf{p}}_{\bar{n}}^{\prime 2}}=\sqrt{m_{\bar{n}}^2+{\mathbf{p}}_{\bar{n}}^2}=E_{\bar{n}}.$$

$$E_N^{\prime }=\sqrt{m_N^2+{\mathbf{p}}_N^{\prime 2}}+V_N(r)\approx m_N+\frac{{\mathbf{p}}_N^{\prime 2}}{2m_N}+V_N(r),$$

$$V_N^i(r)=-\frac{P_{FN}^i{}^2(r)}{2m_N}-{ϵ}_N,$$

$$P_{FN}^i(r)={\left[3{\pi }^2\rho (r)/2\right]}^{1/3}.$$

$$s=(E_{\bar{n}}^{\prime }+E_N^{\prime }{)}^2-({\mathbf{p}}_{\bar{n}}^{\prime }+{\mathbf{p}}_N^{\prime }{)}^2=(E_{\bar{n}}+E_N^{\prime }{)}^2-({\mathbf{p}}_{\bar{n}}^{\prime }+{\mathbf{p}}_N^{\prime }{)}^2.$$

$$E_{\bar{n}}+E_N^{\prime }=m_{\bar{n}}+m_N+\frac{{\mathbf{p}}_{\bar{n}}^2}{2m_{\bar{n}}}+\frac{{\mathbf{p}}_N^{\prime 2}}{2m_N}+V_N(r).$$

$$E_n^{\prime }=\sqrt{m_n^2+{\mathbf{p}}_n^{\prime 2}}+V_N(r)=\sqrt{m_{\bar{n}}^{\ast 2}+{\mathbf{p}}_{\bar{n}}^{\prime 2}}=E_{\bar{n}}^{\prime }.$$

$$E_{ann}+E^{\ast }=E_{\bar{n}}^{{}^{\prime }}+E_N^{{}^{\prime }},$$

$$E^{\ast }=T_{FN}^i-T_N^i,$$

$$E_{ann}+E^{\ast }=m_{\bar{n}}+m_N-{ϵ}_N,$$

$${\rho }_{WS}^{Ar}=(1+\exp (\frac{r-3.6894}{0.5227}){)}^{-1};$$

$$E^{\ast }=T_{Fn}^i-T_{Fn}^i+T_{FN}^i-T_N^i$$

$$E_{ann}+E^{\ast }=m_{\bar{n}}+m_N-2ϵ=1.866\mathrm{GeV}.$$

Full text

I Introduction

I.1 Theoretical Background

Baryon number ($\mathcal{B}$) violation is the only remaining component of the Sakharov conditions (Sakharov:1967dj, ) which has yet to be confirmed experimentally, an innate requirement for proper baryogenesis. While it was shown (tHooft:1976rip, ) that $\mathcal{B}$ itself is not an exact symmetry of the Standard Model (SM) and is only infinitesimally violated via $\mathcal{B}-\mathcal{L}$-conserving non-pertubative electroweak instanton processes, this marginal infraction does not appear to naturally explain the observed excess of matter over antimatter. Thus, it is possible that not only must $\mathcal{B}$ be violated, but that $\mathcal{B}-\mathcal{L}$ must too be broken via some beyond Standard Model (BSM) mechanism.

Such explicit $\mathcal{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 (Nussinov:2001rb, ; Arnold:2012sd, ), this may not be the case for $d=9$ operators with the basic structure $c\cdot\,qqq\bar{q}\bar{q}\bar{q}/M^{5}$, where the high-mass scale $M$ can be rather low at perhaps $\sim 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\rightarrow\bar{n}$). Popular (modern and minimal) BSM extensions (Kuzmin:1985mm, ; Babu:2006xc, ; Babu:2013yca, ; Grojean:2018fus, ) permitting such a process can dynamically create a proper baryon abundance in the early universe, even while predicting a reasonable and possibly observable 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 (Rinaldi:2018osy, ) and other interesting, highly general recent works (Anca, ; Grojean:2018fus, ; Nesvizhevsky:2018tft, ; Berezhiani:2018pcp, ).

I.2 DUNE and Intranuclear $n\rightarrow\bar{n}$ in ${\isotope[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 (${\isotope[40][18]{Ar}}$) inside its 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 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\rightarrow\bar{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, $\tau_{n\rightarrow\bar{n}}$. Taking account of the bound nature of the intranuclear $n$, the inclusion of a so called intranuclear suppression factor $T_{R}$ converts this free transformation time into a bound transformation time $T=T_{R}\,\tau_{n\rightarrow\bar{n}}^{2}~{}$, where $T_{R}$ is quite large. Some details of this derivation, along with the first calculation of this factor for the ${\isotope[40][18]{Ar}}$ nucleus, will be described in some detail in later sections.

Previous searches for $n\rightarrow\bar{n}$ have been carried out using both free (BaldoCeolin:1994jz, ) and bound (Jones84, ; Takita86, ; Berger90, ; Chung02, ; Abe15, ; Bergevin10, ; Aharmim:2017jna, ) $n$’s. The most successful thus far has been Super-Kamiokande’s (Abe15, ) bound water-Cherenkov $n\rightarrow\bar{n}$ search within \isotope[16][8]O, which, given $24$ candidate events and an expected background of $24.1$ atmospheric neutrino events over $\sim 4$ years within $22.5$ kT of fiducial mass, set a limit on the intranuclear $n$ lifetime of $1.9\times 10^{32}$ years, which, when converted to a free transformation time, became $\geq 2.7\times 10^{8}s$. The best free $n$ search occurred at the Institut Laue-Langevin (BaldoCeolin:1994jz, ), and puts a lower limit at roughly the same order of magnitude with no apparent backgrounds.

I.3 Past Intranuclear $n\rightarrow\bar{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 take care to model BSM signals and backgrounds as completely, consistently, and rigorously as necessary by employing limited approximations which attempt to preserve as much physics as possible. This is especially true given nontrivial automated triggering schemes planned for future rare event searches. Unfortunately, in the case of many previous $n\rightarrow\bar{n}$ studies, this has not altogether been the case.

While previous work in ${\isotope[40][18]{Ar}}$ by Hewes (Hewes:2017xtr, ) and others from the DUNE Nucleon Decay Working Group using newly-constructed modules within the GENIE (Andreopoulos:2015wxa, ) 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\rightarrow\bar{n}$ signal within the GENIE default model are not entirely correct. Some approximate notions include:

The assumption that the annihilation occurs along the nuclear density distribution of the nucleus, as discussed in (Hewes:2017xtr, ) and other works, even while the outcomes of (Friedman:2008es, ) are referenced openly and often used; this ignores the surface dominance of the transformation and subsequent annihilation, a key prediction of (Dover:1996ee, ; Friedman:2008es, ). Using two models (including a more holistic quantum-mechanical one similar in spirit to (Friedman:2008es, )), these two assumptions are tested, showing moderate disagreement. 2. 2.

Employing a single nucleon momentum distribution described by a nonlocal, relativistic Bodek-Ritchie (Bodek:1980ar, ; Andreopoulos:2015wxa, ) Fermi gas (the radial dependence of the annihilating (anti)nucleons’ momenta are ignored). 3. 3.

Only $\sim 10$ annihilation channels (a la (Abe15, )) are assumed to be necessary to describe the annihilation products. This seems low, as $\sim 100$ are known, many of them containing heavier resonances; these heavier species can be responsible for $\sim 40\%$ of all pion ($\pi^{0,\pm}$) production. 4. 4.

A true cascade model has not yet been employed, and instead has been approximated as a single effective interaction (GENIE’s Intranuke hA2015 (Andreopoulos:2015wxa, ) was used for previous results) 5. 5.

There also exists no de-excitation model(s) (nucleon evaporation, etc.) within current publicly available builds of GENIE for ${\isotope[40][18]{Ar}}$ 6. 6.

No comparison tests against antinucleon annihilation data have yet been considered. 7. 7.

Only an rough estimation of the nuclear suppression factor of ${\isotope[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 technicalities proceed directly from the secondary nature of the GENIE $n\rightarrow\bar{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 (Gustafson16, ), and detailed explanations of past simulation’s internal processes are quite lacking (Gibin:1996ig, ; Abe15, ), but to individually contend these here is not the goal of this article.

I.4 This Work

Instead, a new and antinucleon-data-tested model is offered here which does not ignore or grossly approximate these effects. In this work, first, the intranuclear oscillation and subsequent annihilation of an $\bar{n}$ is considered using the Sternheimer equation and a proper handling of neutron wave functions. For each neutron shell of ${\isotope[40][18]{Ar}}$, we calculate the induced $\bar{n}$-component by solving the Sternheimer equation with a realistic $\bar{n}A$ interaction, as done in previous works by Dover, Gal, and Richard. From this, the shell-by-shell and average radial annihilation probability distributions are derived, the first calculation of its kind for ${\isotope[40][18]{Ar}}$ and now an integral part of the Monte Carlo simulations discussed later in this work. As a consequence of this study, the first calculation of the intranuclear $n\rightarrow\bar{n}$ transformation suppression factor for ${\isotope[40][18]{Ar}}$ has been completed, finding a value of $T_{R}^{Ar}\sim 5.6\times 10^{22}\,\mathrm{s}^{-1}$, some $\sim 12\%$ smaller than previous coarse estimates taken from \isotope[56]Fe (Hewes:2017xtr, ). Secondly, the first discussion and implementation of novel annihilation pair dynamics in the presence of a realistic $\bar{n}$ potential are completed, where mass defects and momentum modification have now been taken into account within these simulations in a self-consistent way. Third, we consider how changes to our model compare to experiment and past simulation work in extranuclear $\bar{p}$\isotope[12]C annihilation at rest (Golubeva:2018mrz, ) before moving on to consider extranuclear $\bar{n}C$, and then intranuclear $n\rightarrow\bar{n}$ and $\bar{n}N$ annihilation within ${\isotope[40][18]{Ar}}$. Finally, we analyze our investigation and slight perturbation of GENIE’s publicly available $n\rightarrow\bar{n}$ module, attempting to understand the commonalities and differences across both simulations; critically, a discussion is begun as to how ignorance of some particularly physically relevant correlations may affect the eventual feasibility of a true event’s observation in DUNE.

II The Intranuclear $\bar{n}$ Lifetime of Deuterium and $\isotope[40]{Ar}$

II.1 Concepts and Pertinent Questions

Assuming that the $n\rightarrow\bar{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 an $\bar{n}$, an annihilation will eventually take place within the nucleus, releasing $\sim 2$ GeV of rest-mass energy, from which $\sim 4\mbox{-}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:

When a $n$ tentatively becomes an $\bar{n}$, it ceases feeling a smooth potential of $\lesssim 50$ MeV and instead experiences a (complex) potential whose magnitude is $\gtrsim 100$ MeV. How much is such a transformation suppressed by this change in potential? 2. 2.

A deep annihilation could produce multiple fragments with the primary mesons ultimately being absorbed. Alternatively, a peripheral annihilation (Dover:1996ee, ) 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? 3. 3.

Many measurements were accumulated at Brookhaven throughout the 1960-70s, and later at CERN thanks to the LEAR facility (1982-1996), which benefited from a pure, intense, and cooled antiproton ($\bar{p}$) beam. For a review, see, e.g., (Amsler:1991cf, ; Klempt:2002ap, ). However, LEAR was shared by many experiments with various aims dealing with fundamental symmetries, strangeness physics, exotic mesons, etc., and experiments providing knowledge of vital systematics for $\bar{N}N$ and $\bar{N}A$ measurements were not given top priority. With this in mind, is our knowledge of antinucleon-nucleon ($\bar{N}N$) and antinucleon-nucleus $(\bar{N}A)$ interactions sufficient to carry out such an investigation? 4. 4.

Some concerns have been expressed about the reliability of the estimate of the $n\rightarrow\bar{n}$ lifetime inside nuclei within a straightforward nuclear-physics approach; see, e.g., (Costa:1983wc, ; Oosterhof:2019dlo, ). Can one face these criticisms and demonstrate stable and consistent results?

In this section, a review of these questions is discussed with methods based on the Sternheimer equation (see, e.g., (1970PhRvA…1..321S, ), and refs. therein), as used by Dover et al.Friedman et al.  (Dover:1982wv, ; Dover:1985hk, ; Dover:1994ha, ; Friedman:2008es, ) throughout past discussions of these topics. Such methods are then applied to the $\isotope[40]{Ar}$ isotope, which comprises the main component of the DUNE detector. For completeness, the main steps of these calculations are repeated.

II.2 Lifetime of the Deuteron

As a brief warm-up, consider a simplified deuteron, consisting of a pure $s$-wave bound state of a proton ($p$) and neutron ($n$). One may adopt the wave function by Hulthen, which has been tuned to reproduce the correct binding energy and spatial extension. It reads (see, e.g., (Hulthen:220876, ; Adler:1975ga, ))

[TABLE]

where $N_{n}$ is a normalizing factor, $r$ is in GeV**-1**, $\lambda_{1}=0.2316\,\hbar c$ and $\lambda_{2}=5.98\,\lambda_{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\rightarrow\bar{n}$ transformations, the wave function becomes

[TABLE]

Assuming an arbitrary strength $\gamma=1/\tau$ for the $n\rightarrow\bar{n}$ transition, the induced $\bar{n}$ component $w$ is given by the Sternheimer equation

[TABLE]

This gives an exact estimate of the first-order correction to the wave function and, hence, of the second order correction to the energy without involving a summation over the unperturbed states. Here, $E$ is the unperturbed energy and $V$ the antineutron-proton optical potential, resulting in a width

[TABLE]

This immediately implies the scaling law $\Gamma\propto\gamma^{2}$, or for the lifetime $T=\Gamma^{-1}$ of the deuteron

[TABLE]

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 antinucleon-nucleon 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 (Friedman:2008es, ), which is similar to the one used in the earlier studies (Dover:1982wv, ; Dover:1985hk, ), 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

[TABLE]

subject to $u_{n}(0)=u_{n}(\infty)=0$. A method adapted from aircraft engineering (Multhopp, ) 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)=\tilde{u}_{n}(x)$, of an expansion

[TABLE]

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,\dots N$. This results in a $N\times 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})\,\tilde{u}_{n}(x_{i})$. See, for instance, (Richard:1992uk, ), for an application to quarkonium in potential models. For the Sternheimer equation (3), one gets a simple matrix equation

[TABLE]

where the $N\times N$ matrix $\bar{A}$ is the discretized Hamiltonian for the $\bar{n}$, $E$ the $n$ energy, and $W=\{(1-x_{i})\,\tilde{w}(x_{i})\}$ is the vector containing the $\bar{n}$ 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 (Fermi:1959sa, ), the long-range part of the $\bar{N}N$ potential in isospin $I$ is deduced from the $NN$ 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 on $\bar{p}$ scattering and protonium. The so-called DR2 model (Dover:1980pd, ; Richard:1982zr, ) has been adopted, and it was checked that variants such as the Kohno-Weise potential (Kohno:1986fk, ) produce very similar results. Of course, the isospin $I=1$ of the potential is used for the full calculation of the deuteron lifetime.

However, an annihilation requires both the presence of a neighbor nucleon and an $\bar{n}$, so the annihilation must take place at a slightly less peripheral site than the $\bar{n}$ density itself. For the deuteron, the reduced lifetime is estimated to be about $T_{R}^{D}\simeq 3\times 10^{22}\,\mathrm{s}^{-1}$. This result shows great consistency with other recent calculations such as (Kopeliovich:2011aa, ). A calculation by Oosterhof et al. (Oosterhof:2019dlo, ), based of chiral effective field theory, is in some disagreement, but it has been revisited very recently by Haidenbauer et al. (Haidenbauer:2019fyd, ) who found a perfect agreement with the old estimate by Dover et al. (Dover:1982wv, ).

This estimate is remarkably stable. For instance, increasing the core of the $\bar{n}p$ interaction by a factor 10 results in only a $20$% increase of $T_{R}$. Even with a large $|V|$, the $\bar{n}$ 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-energy $\bar{p}$ 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.111JMR thanks Femke Oosterhof for discussions on this point.

In Ref. Dover:1982wv , it is shown that adding a realistic $D$-wave component, and a Sternheimer equation attached to it, does not modify the result significantly.

II.3 Lifetime of the $\isotope[40]{Ar}$ nucleus

In Refs. (Dover:1982wv, ; Friedman:2008es, ), there are estimates of the lifetime of nuclei that were important to analyze some past underground experiments, in particular \isotope[16]O (Abe15, ) and \isotope[56]Fe (Chung02, ). In the present paper, we study the case of \isotope[40]Ar which is relevant for the DUNE experiment.

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., (Ajzenberg-Selove:1977gcr, ). In the present study, besides a more efficient handling of the Sternheimer equation, our variant consists in using the strategy outlined within and neutron wave functions from (Bolsterli:1972zz, ); these wave functions correspond to a fit of the main static properties of the nucleus of interest. 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.

The second ingredient of our work is the $\bar{n}$-nucleus potential. With the noticeable exception of the OBELIX collaboration having studied antineutron scattering Ableev:1994ac ; Iazzi:2000rk , most data deals with the $\bar{p}$-nucleus interaction, either via $\bar{p}$ scattering or antiprotonic atom formation. The question is whether one can reasonably assume that the $\bar{n}$- and $\bar{p}$-nucleus potentials are nearly equal.

The most striking feature of $\bar{N}N$ cross-sections is the smallness of the charge-exchange component, already stressed by the team having discovered the antiproton (Chamberlain:1957, ), and confirmed in further measurements (Klempt:2002ap, ). Indeed, a one-pion-exchange alone would make the charge-exchange the largest contribution to the total cross section. This implies a large cancellation of the isospin $I=0$ and $I=1$ amplitudes, resulting in a somewhat isospin-independent $\bar{N}N$ interaction, and is confirmed by the available experimental data. If one compares the $\bar{p}p$ and $\bar{n}p$ total cross sections, they almost coincide, and differ only in tentative extrapolations towards lower energies; see Iazzi:2000rk ; Bugg:1987nq . This is confirmed by a comparison of the angular distribution of $\bar{p}$ scattering on the isostopes \isotope[16]O and \isotope[18]O at the same energy by the PS184 collaboration, with practically no differences Bruge:1986fd .

Considering this, we shall soon see in the text which follows that changes in the $\bar{n}$-nucleus interaction, such as its departure from understood $\bar{p}$-nucleus interactions, results in only very small modifications of the estimated lifetimes. Of course, the $\bar{p}$-nucleus potential fitting scattering experiments and antiprotonic-atom data is the strong-interaction part, and Coulomb effects have been adequately removed.

For each $n$ shell, there is an induced $\bar{n}$ wave function governed by an equation analogous to (3), with a centrifugal term for $p$, $d$ and $f$ states, where $V$ is now the $\bar{n}$-\isotope[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 that $\bar{n}$ are produced in the tail of the $n$ distribution, with subsequent annihilation at the surface of the matter distribution; the pattern is similar for all other shells.

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, though still decays exponentially as in all other shells.222Some questions have arisen from colleagues regarding the possibility of an annihilation occurring on another nucleus within the larger (detector) medium; this is not permitted due to the exponential decay of the $\bar{n}$ density distribution over internuclear distances, just as a $n$ is not found inside another nucleus.

The resulting radial $\bar{n}$ 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_{R}^{Ar}\sim 5.6\times 10^{22}\,\mathrm{s}^{-1}$. As in the case of the deuteron, this value is remarkably stable against changes in the parameters of the $\bar{n}$-nucleon interaction. Thus, a similar uncertainty of $\sim 20\,\%$ is estimated.

This stability in the width can be understood: if one increases the absorptive potential, the $n\to\bar{n}$ transition is more suppressed, but, on the other hand, the antineutron annihilates more efficiently. In Fig. 5 the factor $\gamma$ is shown, by which the width of the $1d_{5/2}$ level is modified when the real $\bar{n}$-nucleus potential is multiplied by $f_{r}$ and its imaginary part by $f_{i}$. If one changes these values by $\pm 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\%$. The same pattern is observed for the other levels. A consequence of this stability is that the estimated suppression factor $T_{R}$ keeps the same order of magnitude from one nucleus to the other, even from the deuteron to ${\isotope[40][18]{Ar}}$.

One can also illustrate how intranuclear $n\rightarrow\bar{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. Again, the example shown is $1d_{5/2}$, but the other levels exhibit a similar behavior.

To end this section, it is instructive to compare the method based on an exact solution of the perturbation equation (3) in the approximation where the $n$ and $\bar{n}$ spatial distributions are not distinguished. If, furthermore, the $n$- and $\bar{n}$-nucleus do not differ too much in their real part, then the width is given by (see, e.g., (Chetyrkin:1980ta, )):

[TABLE]

For the low-lying shells of \isotope[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 Intranuclear $\bar{n}$ Dynamics

III.1 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 $\gtrsim 30\mbox{-}40$ MeV (the conditions of applicability of the model are considered in detail in (Golubeva:2018mrz, )). 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)-\epsilon_{N}$, where $T_{F}$ is the Fermi energy of the nucleon, and $\epsilon_{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 $\bar{n}$ with a $\isotope[12]{C}$ nucleus (Golubeva:2018mrz, ), some questions arose about the legitimacy of particular physical approximations which could no longer be ignored, requiring that a more correct accounting of the nuclear environment be attained 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.

III.2 The $\bar{n}$-nucleus intranuclear potential and how $\bar{n}$ modifies the nuclear medium

The influence of the nuclear environment on an incoming extranuclear (in the case of ESS using \isotope[12]C (Golubeva:2018mrz, )) or intranuclear ${\bar{n}}$ leads to the modification of the $\bar{n}$’s vacuum four-momentum ${\tilde{p}}_{\bar{n}}=(E_{\bar{n}},{\bf p}_{\bar{n}})$ (${\tilde{p}}^{2}_{\bar{n}}=E^{2}_{\bar{n}}-{\bf p}^{2}_{\bar{n}}=m_{n}^{2}=m_{\bar{n}}^{2}$) inside the target nucleus due to an effective scalar attractive nuclear potential of the form (Cassing NPA700, )

[TABLE]

where $\rho(r)$ the local nucleon density normalized to the atomic mass number $A$ of the nucleus, $\rho_{0}$ is the saturation density, $V_{0}$ is the ${\bar{n}}$ potential depth at this density, and $r$ is the distance between the ${\bar{n}}$ and the center of the nucleus.

Using inferences from Fig. 5, one may assume that the $\bar{n}A$ and $\bar{p}A$ nuclear potentials are effectively the same. With this potential, a parameter of the model, the total ${\bar{n}}$ energy $E_{\bar{n}}^{\prime}$ in the nuclear interior of ordinary nuclei can be expressed in terms of its in-medium mass $m^{*}_{\bar{n}}$, defined as (Cassing NPA700, ; Sibirtsev NPA632, ) 333)The potential $U_{\bar{n}}$ is the effective $\bar{n}A$ scalar potential. The value of this potential in the center of the nucleus, $V_{0}$, is actually a free parameter of the model. This determines the total in-medium $\bar{n}$ energy through the “free space” dispersion relation 13 and is not the usual Lorentz scalar potential $U_{S}^{\bar{N}}$, determining along with the Lorentz vector potential $U_{V}^{\bar{N}}$, the total in-medium antinucleon energy $E_{\bar{N}}^{\prime}$ via the dispersion relation $E_{\bar{N}}^{\prime}=\sqrt{(m_{N}-U_{S}^{\bar{N}})^{2}+{\bf p}_{\bar{N}}^{{\prime}2}}-U_{V}^{\bar{N}}$ (Mishustin PRC71, ).)

[TABLE]

and its in-medium three-momentum ${\bf p}_{\bar{n}}^{\prime}$, as in the free particle case, is (Sibirtsev NPA632, )

[TABLE]

Analysis (Sibirtsev NPA632, ) of $\bar{p}$ production in proton–nucleus and nucleus–nucleus collisions at kinetic energies of a several GeV showed that the $\bar{p}$ potential at normal nuclear matter density is in the range of $-100$ to $-150$ MeV for outgoing $\bar{p}$ momenta below 2.5 GeV$/c$. Studies of $\bar{p}$ production at AGS energies (Koch PLB256, ; Spieles PRC53, ) suggest ${\bar{p}}$ potentials of $\simeq-250$ MeV and $\simeq-170$ MeV at density $\rho_{0}$ for $\bar{p}$ annihilation events at rest with respect to the nuclear matter and for ${\bar{p}}$ with momentum of 1 GeV$/c$, respectively. The real parts of an $\bar{p}$ optical potential in the center of the nucleus of $-(150\pm 30)$ MeV and of $-(220\pm 70)$ MeV were extracted in (Larionov PRC80, ) from the data on ${\bar{p}}$ absorption cross sections on nuclei and on the annihilation spectra of $\pi^{+}$’s and $p$’s, correspondingly. Combined analysis (Friedman NPA761, ) of data on antiprotonic $X$-rays and of radiochemical data showed that at the center of the nucleus the ${\bar{p}}$ potential is approximately $-110$ MeV in depth. So, in spite of various attempts to fix this potential, its depth at density $\rho_{0}$ still remains rather uncertain presently. For the sake of definiteness, in the subsequent calculations, a realistic value of $V_{0}=-150$ MeV will be used within Eq. (11).

The in-medium momentum ${\bf p}_{\bar{n}}^{\prime}$ is related to the vacuum momentum ${\bf p}_{\bar{n}}$ by the following expression:

[TABLE]

For example, with $V_{0}=-150$ MeV, this shows that for ${\bar{n}}$ annihilation at rest, i.e. when ${\bf p}_{\bar{n}}=0$, the $\bar{n}$ momentum $|{\bf p}_{\bar{n}}^{\prime}|$ 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 model444The usual model of a degenerate Fermi gas is one of free nucleons enclosed within a spherical potential well with a sharp border. The intranuclear nucleons fill all energy levels of this well and have momentum from [math] to the boundary value $P_{FN}$, which is itself a function of the nuclear density. In this case, the momentum of the nucleon does not depend on its position in the nucleus, and so we are thus considering a nonlocal (degenerate) Fermi gas model. Usually, in order to take better account of the influence of the diffuse nuclear boundary, the nucleus is divided into concentric layers (zones), each with a constant density and each with boundary values $P_{FN}^{~{}i}$. Thus, there now exists a correlation between the position of the nucleon in the nucleus and its momentum, and so we are now considering a local or zoned (degenerate) Fermi gas model. Ample description of this is contained in (Golubeva:2018mrz, ). used in our MC simulations (Golubeva:2018mrz, ), for the bound target nucleon total energy $E^{\prime}_{N}$ in the medium at the annihilation point $r$ the formula (Magas PRC71, ; Efremov EPJA1, )

[TABLE]

is used, where for every $i$th concentric zone of a spherical nucleus we have

[TABLE]

with

[TABLE]

Here, ${\bf p}_{N}^{\prime}$ is the momentum of the nucleon $N$ ($N=\{p,n\}$) in the Fermi sea, $P_{FN}(r)$ is the boundary Fermi momentum at the local point $r$, and the quantity $\epsilon_{N}\approx 7$ MeV is the average binding energy per nucleon; at this point, $0\leq|{\bf p}_{N}^{\prime}|\leq P_{FN}(r)$. Within the representation of Eqs. (13)–(15), the invariant collision energy $s$ for the interaction of an $\bar{n}$ with a nucleon bound in the nucleus at the point $r$ is

[TABLE]

The total collision energy $E_{\bar{n}}+E^{\prime}_{N}$ entering into the second relation of Eq. (18) in the non-relativistic limit appropriate to our case can be calculated as

[TABLE]

Contrary to the on-shell interaction, for ${\bar{n}}$ 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}-\epsilon_{N}$.

If a bound target $n$ is transformed into an $\bar{n}$ (for the ${\isotope[40][18]{Ar}}$ case), we assume that its total energy $E_{n}^{\prime}$, defined by Eq. (15), is equal to that $E_{\bar{n}}^{\prime}$ of the $\bar{n}$, determined by Eq. (13) above. Namely:

[TABLE]

It is interesting to note that, for $p_{F}(0)=250$ MeV$/c$, Eq. (20) gives for the ${\bar{n}}$ momentum $|{\bf p}_{\bar{n}}^{\prime}|$ the values of 430 and 40 MeV$/c$ if the transition $n\to{\bar{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 ${\isotope[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 (Golubeva:2018mrz, ).

IV.1 Fundamentals

Most all of the background concerning our MC generator can be found in our recent work (Golubeva:2018mrz, ), however, we will mention some improvements which have been implemented here. We can briefly summarize the content of the model as follows:

Unlike the optical approach used in our previous work (Golubeva:2018mrz, ) to calculate the radial annihilation position distribution of a slow external antineutron annihilating on \isotope[12]C, in this work, an intranuclear annihilation point is taken from a corresponding probability distribution shown later for ${\isotope[40][18]{Ar}}$ in 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. 2.

The annihilation occurs, producing $\pi$’s and higher mass resonances such as $\eta,\omega$, and $\rho$’s, according to tabulated channels with branching ratios taken from (Golubeva:2018mrz, ). After the decay of all resonances, on average $\sim 4\mbox{-}5$ mesons are produced from the initial annihilation process. 3. 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. 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.

IV.2 Improvements from OBELIX and Crystal Barrel Data

The annihilation model was created in 1992 (Golubeva:1992tr, ), originally using experimental branching ratios obtained before this time. Later analysis of experimental data from $\bar{p}p$ annihilation at rest obtained from the LEAR (CERN) $\bar{p}$ beam by the OBELIX (Salvini:2004gz, ) and Crystal Barrel (Amsler:2003bq, ) 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 1 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 (Golubeva:2018mrz, ), while the third shows the new branching fraction of this channel changed in accordance with the experimental data (Salvini:2004gz, ; Amsler:2003bq, ). Once summed, all channels are then renormalized to unity with these new fractions taken into account. A comparison of the elementary processes for $\bar{p}p$ annihilation was carried out taking into account experimental data (Salvini:2004gz, ; Amsler:2003bq, ) and annihilation simulations from our recent work (Golubeva:2018mrz, ).

Table 2 shows the average multiplicities of mesons produced in $\bar{p}p$ annihilation at rest. The first column re-presents the simulation results from Table IV in (Golubeva:2018mrz, ), while the second column presents the results of modeling when taking into account data from (Salvini:2004gz, ; Amsler:2003bq, ). The third column also re-presents the experimental data itself (Klempt:2005pp, ; Minor, ; Levman:1979gg, ; Hamatsu:1976qz, ; Golubeva:2018mrz, ) for ease of comparison. It follows from Table 2 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 all $\bar{p}p$ annihilations.

IV.3 Model modifications and a new comparison with experimental $\bar{p}C$ annihilation at rest

The branching fraction modifications shown in Table 1 were implemented into the optical-cascade model, most recently discussed in (Golubeva:2018mrz, ). Now, let’s analyze how these changes effect the description of experimental data for $\bar{p}C$ annihilation at rest. The first line of Table 3 presents the experimental multiplicities of the emitted $\pi$’s and the energy carried away by those $\pi$’s and $\gamma$’s. The second line of the table shows the results of a calculation with our original optical-cascade model from (Golubeva:2018mrz, ) 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 of 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.

From analysis of the table, it follows that the average multiplicities vary slightly with the modification of the annihilation table and dynamics. Multiplicities of $\pi$’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 from $\bar{p}A$ 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 observation of predominately spherically symmetric, multipionic topologies (so-called $\pi$-stars), we will focus on the characteristics of these topologies and their associated quantities.

Fig. 7 shows the spectrum of $\pi^{+}$ emitted during $\bar{p}{\rm C}$ annihilation at rest in comparison to experimental data from (Minor, ) (green triangles), and (Mcgaughey:1986kz, ) (blue squares). Just as in the average multiplicities shown in Table 3, 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 $\Delta$-resonance).

In our proceeding simulations for the ${\isotope[40][18]{Ar}}$ nucleus to be discussed in later sections, we will use a model analogous to Calculation $\#2$ that includes effects related to the influence of the nuclear environment.

V New results for simulations of $\bar{n}A$ annihilation

V.1 Changes in extranuclear $\bar{n}C$ annihilation simulations

Despite the absence of significant differences in the description of available $\bar{p}{\rm C}$ 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 (Golubeva:2018mrz, ), 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

[TABLE]

where $E^{{}^{\prime}}_{\bar{n}}$ is the total energy of the $\bar{n}$ inside the nucleus at the point of annihilation, $E^{{}^{\prime}}_{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

[TABLE]

varying from [math] to $T^{i}_{FN}$, where $i$ is the zone number in which the annihilation takes place, $T^{i}_{FN}$ is the boundary Fermi energy of the $i$-th zone, and $T^{i}_{N}$ 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

[TABLE]

where $\epsilon_{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 true $\bar{n}$-potential 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 the $\bar{n}$-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, the $\bar{n}$ was assumed to come from a transformation down-range of a cold $n$ source with a mean energy of only $\sim$meV, ignoring the $\bar{n}$-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 the $\bar{n}$ is defined by expression (12) and the momentum of the $\bar{n}$ 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^{{}^{\prime}}_{\bar{n}}-P^{{}^{\prime}}_{N}|$ to $|P^{{}^{\prime}}_{\bar{n}}+P^{{}^{\prime}}_{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 the $\bar{n}$-potential is taken to be [math] 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.

The most impressive new figures in consideration of $\bar{n}{\rm C}$ extranuclear annihilation are most likely Figs. 10 and 11. Here we see the total available initial and final mesonic/pionic and photonic parameter space for an $n\rightarrow\bar{n}$ signal, most commonly shown via total momentum versus invariant mass plots (similar to Fig. 2 in (Abe15, )) at the annihilation point before any nuclear transport is completed. 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, $\Delta$-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. Overall, the new parameter space is slightly more spatially if not statistically constrained, which could lead to higher hypothetical experimental efficiencies for future experiments.

V.2 Intranuclear $\bar{n}{\rm Ar}$ annihilation simulations

As mentioned previously, there now exist multiple generators for $n\rightarrow\bar{n}$ within the particle physics community, notably developed in (Hewes:2017xtr, ) using GENIE (Andreopoulos:2015wxa, ). While we believe that the demonstration of our independent generator’s capabilities (Golubeva:2018mrz, ) in the reproduction of antinucleon data is well-established for $\isotope[12]{C}$ (Golubeva:2018mrz, ), such a complete set of physical observables does not readily exist to constrain the model for larger nuclei, especially not for ${\isotope[40][18]{Ar}}$. Thus, out of a need for ample comparisons, we endeavor to show the commonalities and differences between each of these $n\rightarrow\bar{n}$ generators for intranuclear $\bar{n}\,\isotope[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\rightarrow\bar{n}$ generator module to utilize a noninteracting local Fermi gas nuclear model along with a full intranuclear cascade. The inclusion of an $\bar{n}$-potential within GENIE has not yet been investigated; implementation of the modern annihilation position probability distribution (Sec. II.3) is currently underway. While none of these comparisons across generators are ever to be exact, their approximately equivalent formalism can serve to inform us as to the stability of quantities which characterize the possible final state topologies of a true $n\rightarrow\bar{n}$ signal event 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 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; this position is thrown, a nearest neighbor nucleon found within a given zone (and then ”moved” to the annihilation position), a Fermi momentum computed for each of the pair given their initial positions, and then a total phase space calculated. All GENIEv3.0.6 events source their annihilation positions from a smooth Woods-Saxon (nuclear density) distribution (Andreopoulos:2015wxa, ) incredibly similar to our own continuous parameterization for the nuclear density, from which we derive our local zone densities:

[TABLE]

GENIE also similarly throws momenta from (non)local single nucleon momentum distributions. When this curve is multiplied by $r^{2}$, one generates the blue annihilation probability distribution. These curves effectively demonstrate how even the most simple of assumptions, some 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 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\rightarrow\bar{n}$ transformation time; the interplay of these quantities and any changes in the final state $\pi$-star topology observable in the DUNE detectors has not yet been completely investigated. For similar plots and discussions, see (Golubeva:1992tr, ; Golubeva:1996ij, ; Golubeva:1997mc, ; Golubeva:1997fs, ; Golubeva:2018mrz, ).

Some of the differences between this work and GENIEv3.0.6’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 mesons’ total energy for this work (using the orange curve in Fig. 12) and GENIE (using the blue curve in Fig. 12), each using a version of a local Fermi gas nuclear model. Like the extranuclear $\bar{n}$ annihilation described on $\isotope[12]{C}$ above, energy balance in the annihilation point is given as $E_{ann}+E^{*}=E^{{}^{\prime}}_{\bar{n}}+E^{{}^{\prime}}_{N}$; taking into account that we have $E^{{}^{\prime}}_{\bar{n}}=E^{{}^{\prime}}_{n}$, and that

[TABLE]

we have the total energy

[TABLE]

It can be seen in Fig. 13 that our distribution of energies available to an annihilation is always less than $1.866$ GeV. GENIE’s distribution (which assumes a similar binding energy per nucleon as our model) 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 $\sim 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 nonlocal Bodek-Ritchie or local Fermi gas nuclear models mentioned here serve effectively as a set of initial conditions which enable certain nucleon momentum and radial position correlations (or lack thereof). The nonlocal Bodek-Ritchie has been considered the default operating model for $n\rightarrow\bar{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 nonlocal 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}(\bar{n})=p_{f}(n)\approx p_{f}(p)$; this is not the case for our model, where $p_{f}(\bar{n})\neq p_{f}(n)\approx p_{f}(p)$ due to the $\bar{n}$-potential. 666For simplicity, throughout this section we have labeled certain plots with “Golubeva-Richard-Paryev” (for original work done with the modern shell model-derived annihilation position probability distribution and modification of the nuclear medium due to $\bar{n}$ interactions), “Golubeva-WS-Paryev” (for original work done with a Woods-Saxon-derived annihilation position probability distribution and nuclear medium modification), and then the various GENIEv3.0.6 model identifiers (Andreopoulos:2015wxa, ).

The most important aspect of correlated behavior which has been previously unaccounted for in GENIE-affiliated work on $n\rightarrow\bar{n}$ is that of initial (anti)nucleon momentum and radius. In Figs. 15, we show comparisons between our and GENIE’s local Fermi gas nuclear models, which by definition preserves these correlations, alongside the inherently nonlocal Bodek-Ritchie nuclear model. All figures assume a (zoned or smooth) Woods-Saxon-like annihilation position probability distribution for easier direct comparisons; all outputs from our simulations using the modern annihilation probability distribution from Sec. II appear somewhat similar when graphed in these coordinates, and so are not shown here to conserve space. The $\bar{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 for $\bar{n}$ annihilation partners). All local models 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 nonlocal Bodek-Ritchie nuclear model, shown at the bottom; the asymmetry present in this plot is due to the Fermi momentum cutoff, above which only phenomenologically ’short-range correlated’ $\bar{n}$’s populate.

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 the $\bar{n}$ with the modified nuclear medium.

All of this leads us to consider the available initial (and final) mesonic parameter space, again à la Fig. 2 of (Abe15, ). As seen in Figs. 17, this serves as an initial condition of the annihilation-generated mesons before intranuclear transport; thus, for GENIE, hA/hN2018 models (see Secs. 2.5.4 and 2.5.5 of the GENIE v3.0 manual (Andreopoulos:2015wxa, ) for full discussions) 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 generated 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 resonance decays, while GENIE does not. The disconnected regions toward the left of the plots are signs of single $\pi$ 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) $\pi$’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\rightarrow\bar{n}$ event. A full characterization of these effects across many nuclear model configurations within the DUNE detectors is forthcoming within the DUNE Nucleon Decay Group.

We end the comparison of these generators with Tab. 4 and Figs. 19 and 20, which show many similarities and some differences across them. Multiplicities in Tab. 4 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) models. 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 $\pi^{0,\pm}$ multiplicities, which are partially due to the fact that we predict more $n\rightarrow\bar{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 (Golubeva:2018mrz, ) compared to GENIE (Andreopoulos:2015wxa, ).

To give a more complete context to Table 4, we plot several (absolute magnitude) final state momentum spectra for $\pi^{+}$ and $p$ species, two key constituents in the eventual experimental search for $n\rightarrow\bar{n}$ in DUNE; note that the $\pi^{0,-}$ distributions are quite similar. 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 structure around the $\Delta$-resonance within the hA2018 simulation (recall this models the cascade as a single effective interaction using tabulated reaction rates), a sign of the competition between cross-sections (or rates) of processes such as $\Delta$ 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 $\sim 4\mbox{-}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 detectability threshold in liquid ${\isotope[40][18]{Ar}}$ to be $\gtrsim 100$ MeV (i.e., $\gtrsim 450$ MeV$/c$) (MicroBooNE1025, ), 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 of $\bar{n}A$ annihilation events in service of future BSM $n\rightarrow\bar{n}$ searches. Critical among the findings of this paper is the calculation of a new intranuclear suppression factor for ${\isotope[40][18]{Ar}}$, $T_{R}^{Ar}\sim 5.6\times 10^{22}\,\mathrm{s}^{-1}$, along with a new and associated calculation of the ${\isotope[40][18]{Ar}}$ intranuclear radial annihilation probability distribution. Also, efforts have been made to implement this and other important $\bar{n}N$ 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 reasonable 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\rightarrow\bar{n}$ mean transformation time inside simulated DUNE detector geometries.

Another future step elucidated by this work is the critical nature of current and future collaborations’ endeavors to holistically evaluate and compare various $\bar{n}A$ interactions using common formalisms for the calculation of (anti)nucleon wave functions, radial annihilation probability distributions, and intranuclear suppression factors for all pertinent nuclei. This is a rather monumental task, but one which should be completed in the same way as for other rare decay searches, such as within the $0\nu\beta\beta$ community.

Acknowledgements.

JMR would like to thank Michael Bender and Karim Bennaceur for useful advice about shell-model wave functions, Joseph Cugnon for discussions, and Dario Auterio and his collaborators for interesting discussions about the DUNE experiment. JLB’s work is supported by the U.S. Department of Energy, Office of Science, Office of Workforce Development for Teachers and Scientists, Office of Science Graduate Student Research (SCGSR) program. The SCGSR program is administered by the Oak Ridge Institute for Science and Education for the DOE under contract number DE‐SC0014664. JLB is grateful to the US DOE SCGSR Program and Fermilab for their generous support of this and other work. ESG would like to thank INR, Moscow, for the support of her travel. We all wish to thank the organizers of the INT Neutron Oscillations Workshop in Seattle during November 2017, the ILL-LPSC PPNS Workshop at Grenoble in May 2018, and the ESS Workshop at NORDITA in December 2018 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.