# ΔB=2 neutron decay into antiproton mode n\to \bar pe^+ν(\barν)

∆∆ B = 2 neutron decay into antiproton mode n → ¯ pe + ν (¯ ν ) Xiao-Gang He

1, 2, 3, ∗ and Xiao-Dong Ma

2, 4, † Tsung-Dao Lee Institute, and School of Physics and Astronomy,Shanghai Jiao Tong university, Shanghai 20024, China Department of Physics, National Taiwan University, Taipei 106, Taiwan Physics Division, National Center for Theoretical Sciences, Hsinchu 300, Taiwan School of Nuclear Science and Technology, Lanzhou University, Lanzhou 730000, China

We discuss the unique baryon number violation by two units neutron decay mode n → ¯ pe + χ , with χ being the standard model (SM) neutrino ν or antineutrino ¯ ν or any beyond SM light fermion,in the framework of eﬀective ﬁeld theory. This mode is kinematically allowed but rarely discussedtheoretically or searched for experimentally. We estimate the lower bound on its partial lifetimefrom that of the dinucleon decay np → e + χ per oxygen nucleus O set by the Super-Kamiokandeexperiment, with a conservative bound Γ − n → ¯ pe + χ > . × yrs. We also discuss its characteristicsignature for the future experimental search and astrophysical implications. Introduction.

The baryon-antibaryon asymmetry ofthe universe requires the violation of baryon number ( B ),which is one of the three Sakharov conditions for a suc-cessful baryongenesis mechanism [1]. The baryon numberviolation (BNV) is also a general feature in scenarios ofphysics beyond the standard model like the grand uniﬁedtheories [2–4], in which the | ∆ B | = 1 single nucleon de-cays like proton decay p → e + π are predicted. However,the scale of new physics (NP) associating with | ∆ B | = 1nucleon decays is constrained, by the current experimen-tal data, to be around 10 − GeV, which is impossibleto directly produce those heavy particles at current andfuture high energy colliders.On the other hand, there is a class of scenarios sup-pressing | ∆ B | = 1 nucleon decays but contribute dom-inantly to | ∆ B | = 2 processes like the neutron - an-tineutron ( n − ¯ n ) oscillation, hydrogen - antihydrogen(H − H) oscillation and dinucleon to dimeson/dileptondecays in nuclei (like pp → π + π + , e + e + ) [5–9]. Insuch scenarios, the associated NP scale is lowered tobe around O (1 − ) TeV from the current experimen-tal results on n − ¯ n oscillation [10, 11] and dinucleon todimeon/dilepton decays [12–14], and this TeV scale NPis appealing since they may be directly tested at collid-ers like the LHC and the other proposed ones throughthe search of the mode pp → (cid:96) +1 (cid:96) +2 + 4 jets [15]. In thislatter case with | ∆ B | = 2 dinucleon decays, we noticedthat there is a unique | ∆ B | = 2 single free or bound neu-tron decay mode n → ¯ pe + χ with χ being the standardmodel (SM) neutrino ν or anti-neutrino ¯ ν or any newlight fermion beyond the minimal SM such as a lightsterile neutrino, but such a mode has not yet been con-sidered both theoretically and experimentally. In thisletter we will investigate this neutron decay mode in theframework of eﬀective ﬁeld theory and point out its dis-tinct experimental signature relative to other dinucleondecays for the guidance of future experimental searches. ∗ [email protected] † [email protected] From the perspective of standard model eﬀective ﬁeldtheory (SMEFT), the leading order interactions con-tributing to n → ¯ pe + χ depend on the types of χ . For the χ to be the electron neutrino ν e , n → ¯ pe + ν e can be gener-ated, at leading order, by the same dimension-9 (dim-9)operators mediating the n − ¯ n oscillation through an in-sertion of a SM charged weak current, which is similarto the nuetron-antineutron conversion through the same n − ¯ n oscillation operators did by Gardner and Yan [16–18]. For the muon and tau neutrinos ν µ,τ , the leadingorder interactions for n → ¯ pe + ν µ,τ appear at dim 13, this is because the SM has a good lepton ﬂavor symme-try and therefore n → ¯ pe + ν µ,τ cannot be realized like thedecay n → ¯ pe + ν e by dim-9 interactions.On the other hand, for the antineutrino case, the de-cay mode n → ¯ pe + ¯ ν e,µ,τ also violates lepton number bytwo units ( | ∆ L | = 2) and their leading order interactionsappear at dim 12. After sending the Higgs ﬁeld to itsvacuum expectation value and expanding the quark andlepton doublet, we obtain the relevant eﬀective interac-tions mediating n → ¯ pe + χ with χ (cid:54) = ν e in the so-calledlow energy eﬀective ﬁeld theory (LEFT) below the elec-troweak scale [20, 21], which consisting of two up-typeand four down-type quarks and a charged lepton cur-rent, i.e., having the conﬁguration ( uuddddeχ ). For thedecay n → ¯ pe + ν e , it is a little bit complicated since itcan be generated at leading order from the same dim-9 n − ¯ n oscillation operators with the structure ( uudddd )by inserting a SM four-fermion vertex in one of the fourdown-type quark legs.Once the relevant LEFT operators are obtained, onecan perform a non-perturbative QCD matching for thesix-quark sectors using the baryon chiral perturbationtheory (B χ PT) formalism [22–24]. In this way, one endsup with operators consisting of nucleons, mesons and They are not generated at dim-12 in that the SMEFT has theproperty that the operator’s dimension is even (odd) if | ∆ B − ∆ L | / | ∆ B | = 2 and∆ L = 0 means the dimension of relevant operator is odd. a r X i v : . [ h e p - ph ] J a n eptons and the decay rate can be readily calculated.Here we note that the relevant interactions for pn → e + χ transition is the same as that of the decay n → ¯ pe + χ because these two processes are related to each otherby crossing symmetry. Therefore constraints on thebranching ratio of n → ¯ pe + χ can be obtained by usingknown bound from pn → e + χ transition. If realizing theinteractions via SMEFT operators, these interactionsmay be related to other dinucleon to dilepton transitionssuch as pp → e + e + and nn → ¯ ν ¯ ν . An example of sucha case will be given below. We will present the detailselsewhere for the above procedures to our accompanyinglong paper concerning the | ∆ B = ∆ L | = 2 dinucleon todilepton decays ( pp → (cid:96) + (cid:96) (cid:48) + , pn → (cid:96) + ¯ ν (cid:48) , nn → ¯ ν ¯ ν (cid:48) ) in afull EFT analysis [25]. EFT analysis for n → ¯ pe + χ . In the following, we willtake the EFT procedures outlined in the above and startdirectly from the hadronic level interactions for a model-independent analysis. We take the minimal assumptionthat all the dominant interactions contributing to n → ¯ pe + χ are encoded in the leading order dim-6 operatorswith the ﬁeld conﬁguration ( npeχ ).The relevant dim-6eﬀective | ∆ B | = 2 interactions are classiﬁed into twosectors in terms of the net lepton number | ∆ L | = 0 (for χ = ν ) or | ∆ L | = 2 (for χ = ¯ ν ) as follows L ∆ B =2 = (cid:88) i (cid:101) C ( pn )∆ L =0 ,i (cid:101) O ( pn )∆ L =0 ,i + (cid:88) i C ( pn )∆ L =2 ,i O ( pn )∆ L =2 ,i + h . c . , (1)where all the independent dim-6 operators in each sectorare parametrized as follows n → ¯ pe + ν : (cid:101) O ( pn ) SR = ( p T Cn )( ν L e R ) , (cid:101) O ( pn ) S R = ( p T Cγ n )( ν L e R ) , (cid:101) O ( pn ) VL = ( p T Cγ µ n )( ν L γ µ e L ) , (cid:101) O ( pn ) V L = ( p T Cγ µ γ n )( ν L γ µ e L ) , (cid:101) O ( pn ) T = ( p T Cσ µν n )( ν L σ µν e R ) , (2) n → ¯ pe + ¯ ν : O ( pn ) SL = ( p T Cn )( e T L Cν L ) , O ( pn ) S L = ( p T Cγ n )( e T L Cν L ) , O ( pn ) VL = ( p T Cγ µ n )( e T R Cγ µ ν L ) , O ( pn ) V L = ( p T Cγ µ γ n )( e T R Cγ µ ν L ) , O ( pn ) T = ( p T Cσ µν n )( e T L Cσ µν ν L ) , (3)where C is the charge conjugation matrix satisfying C T = C † = − C and C = −

1, and the ﬂavor of SM left-handed neutrino is suppressed. Except the neutron decay n → ¯ pe + χ , the above interactions can also contribute tothe dinucleon decay pn → e + χ in nuclei and the conver-sion e − p → nχ and e − n → ¯ pχ in the electron-deuteron(e-d) scattering. This correlation of diﬀerent processesdue to the crossing symmetry can help us to estimate the partial lifetime of n → ¯ pe + χ from the experimental lowerbound on that of pn → e + χ per oxygen nucleus O setby the water Cherenkov Super-Kamiokande (SK) exper-iment [12]. In the following, we assume each time thereis only one operator dominant for the decays n → ¯ pe + χ and np → e + χ , then the bound on the latter mode canbe used to extract a bound on the former.As mentioned earlier that the low energy eﬀective in-teractions may be related to NP at high energies, whichcan be parametrized through the SMEFT interactions.Here we give an example how a SMEFT operator cangenerate some of the low energy operators. We ﬁnd thedim-12 operator O S, ( S ) Q L = ( Q i T a CQ jb )( Q k T c CQ ld )( Q m T e CQ nf )( L T g CL (cid:48) h ) × (cid:15) ab (cid:15) cd (cid:15) eg (cid:15) fh T SAA { mn } [ kl ][ ij ] , (4)is of a very interesting case. Where C is again the chargeconjugation matrix, Q and L are the SM quark andlepton doublets, respectively. T SAA { mn } [ ij ][ kl ] = (cid:15) ijm (cid:15) kln + (cid:15) ijn (cid:15) klm is a color tensor responsible for the operator tobe color invariant, (cid:15) ab ( (cid:15) ijk ) is the second (third) ranktotally antisymmetric Levi-Civita symbol.The hadronic counterpart of the six-quark part of O S, ( S ) Q L can be done through the B χ PT matching, whichhas been used before for the dim-9 six-quark n − ¯ n oscil-lation operators in [24]. We ﬁrst identify that the six-quark part of O S, ( S ) Q L belongs to the irreducible repre-sentation ( L , R ) of the two-ﬂavor QCD chiral group SU (2) L × SU (2) R , then its leading oder hadronic coun-terpart is realized by appealing to the nucleon ﬁeld Ψ =( p, n ) T and pion matrix u in such a way that the resultantoperator has the same symmetry property as O S, ( S ) Q L in-cluding the Lorentz covariance, baryon number and chiralproperty. After carrying out the procedures described inthe above, we obtain the leading order hadronic operatortaking the form ( u † ) u L a ( u † ) v L b Ψ T a C [ g × + ˆ g × γ ] Ψ b with u L , v L = 1 , u, d quarks. Expand to zerothorder of pion ﬁelds, we ﬁnd the operator leads to O ( pn ) SL and O ( pn ) S L in Eq. (3) with the coeﬃcients C ( pn ) SL = − g × C S, ( S ) Q L , C ( pn ) S L = − g × C S, ( S ) Q L , (5)where C S, ( S ) Q L ≡ Λ − is the Wilson coeﬃcient of oper-ator O S, ( S ) Q L with Λ NP being the associated NP scale.Here the hadronic low energy constants g × and ˆ g × are pertinent to the non-perturbative QCD matching ofquark level operators in the B χ PT. g × can be relatedto the n − ¯ n oscillation matrix element via chiral sym-metry, the latter has been determined by using latticeQCD method [26], that implies g × ∼ × − GeV ,while the value for ˆ g × has not been calculated. Anaive dimensional analysis would give ˆ g × to be of order ∼ Λ χ / (4 π ) ∼ . × − GeV [27], where Λ χ ∼ . g × , the mass dimension ofˆ g × should be compensated by Λ QCD ∼

200 MeV sinceit is the scale for the non-perturbative QCD eﬀect, then g × ∼ Λ ∼ . × − GeV . We see both estima-tions agree within an O (10) uncertainty. For numericalestimations, we will use this later lower value for illus-trations. The complete details for the above analysis willbe presented in our accompanying long paper on LEFTand SMEFT operators [25].It is also interesting that the above operator O S, ( S ) Q L not only generate n → ¯ pe + ¯ ν and np → (cid:96) + ¯ ν (cid:48) , butalso pp → (cid:96) + (cid:96) (cid:48) + and nn → ¯ ν ¯ ν (cid:48) . If the future SKexperiments could ﬁnd the | ∆ B | = 2 signals from allthose three diﬀerent channels, it probably points to theNP incorporated in this operator. Decay rate calculations . From the above interactionsin Eqs. (1,2,3), it is straightforward to calculate the tran-sition amplitudes for n ( k ) → ¯ p ( p ) e + ( p ) χ ( p ) and hence-forth the decay rates. The decay rate can be written asΓ n → ¯ pe + χ = 12 m n π m n (cid:90) ds (cid:90) dt (cid:12)(cid:12) M n → ¯ pe + χ (cid:12)(cid:12) , (6)where m n is the neutron mass, the Mandelstam variablesare deﬁned as s = ( p + p ) , t = ( k − p ) and u =( k − p ) .Taking into the interactions in Eqs. (1,2,3) into con-sideration, after calculating the amplitudes and ﬁnishingthe phase space integration, and assuming one term dom-inates at a time, we obtain the following decay rate foreach type of the couplings,Γ n → ¯ pe + ν = m n π (cid:110) (cid:15) (cid:12)(cid:12) (cid:101) C ( pn ) SR (cid:12)(cid:12) , (cid:15) (cid:12)(cid:12) (cid:101) C ( pn ) S R (cid:12)(cid:12) ,(cid:15) (cid:12)(cid:12) (cid:101) C ( pn ) VL (cid:12)(cid:12) , (cid:15) (cid:12)(cid:12) (cid:101) C ( pn ) V L (cid:12)(cid:12) , (cid:15) (cid:12)(cid:12) (cid:101) C ( pn ) T (cid:12)(cid:12) (cid:111) , Γ n → ¯ pe + ¯ ν = m n π (cid:110) (cid:15) (cid:12)(cid:12) C ( pn ) SL (cid:12)(cid:12) , (cid:15) (cid:12)(cid:12) C ( pn ) S L (cid:12)(cid:12) ,(cid:15) (cid:12)(cid:12) C ( pn ) VL (cid:12)(cid:12) , (cid:15) (cid:12)(cid:12) C ( pn ) V L (cid:12)(cid:12) , (cid:15) (cid:12)(cid:12) C ( pn ) T (cid:12)(cid:12) (cid:111) , (7)If several terms exist simultaneously, there are in generalinterference terms. In our numerical estimate later, weassume one term dominates at a time and therefore theinterference terms are neglected. The dimensionless pa-rameters (cid:15) i stem from the phase space integration overthe kinematic variables and take the following numericalvalues (cid:15) = 5 . × − , (cid:15) = 7 . × − , (cid:15) = 5 . × − ,(cid:15) = 1 . × − , (cid:15) = 6 . × − . (8)where the suppression of (cid:15) i is because of the small phasespace for the process.We now analyse the dinucleon to dilepton decay pn → e + χ in nucleus to obtain a correlation to the pro-cesses discussed above. From the same interactions in Eqs. (1,2,3), we can calculate the transition rate for pn → e + χ in nucleus in the following manner [28]Γ pn → e + χ = 1(2 π ) √ ρ p ρ n (cid:90) d k d k ρ p ( k ) ρ n ( k ) × v rel . (1 − v · v ) σ ( pn → e + χ ) , (9)where ρ p ( k )( ρ n ( k )) is the proton (neutron) density dis-tribution in momentum space and ρ p ( ρ n ) is the aver-age proton (neutron) density. v ( v ) is the velocityof the nucleon p ( n ) and v rel . is their relative velocity v rel . = | v − v | . The total cross-section for the free nu-cleon scattering process p ( k ) n ( k ) → e + ( p ) χ ( p ) takesthe form σ ( pn → e + χ ) = 14 E E v rel . (cid:90) d Π (cid:12)(cid:12) M Npn → e + χ (cid:12)(cid:12) , (10)where E ( E ) is the energy of the initial state nucleon p ( n ). d Π is the relativistically invariant two-body phasespace. For the oxygen nucleus O, we neglect the smalleﬀect due to the eﬀects of nucleon Fermi motion andnuclear binding energy [29]. Under the quasi-static ap-proximation of nucleons, then the transition rate reducesinto Γ pn → e + χ = ρ N m N (cid:12)(cid:12) M pn → e + χ (cid:12)(cid:12) Π , (11)where the average nuclear matter density ρ N approxi-mately equals 0 .

25 fm − for either proton or neutron,and m N = ( m n + m p ) /

2. The two-body ﬁnal statephase factor Π = (1 / π ) (cid:0) − m e / ( m n + m p ) (cid:1) . Therelevant amplitudes are calculated from the interactionsin Eqs. (1,2,3). After ﬁnishing the phase space integra-tion and combining all pieces together, we obtainΓ pn → e + ν = ρ N m N π (cid:110) η (cid:12)(cid:12) (cid:101) C ( pn ) SR (cid:12)(cid:12) , η (cid:12)(cid:12) (cid:101) C ( pn ) S R (cid:12)(cid:12) ,η (cid:12)(cid:12) (cid:101) C ( pn ) VL (cid:12)(cid:12) , η (cid:12)(cid:12) (cid:101) C ( pn ) V L (cid:12)(cid:12) , η (cid:12)(cid:12) (cid:101) C ( pn ) T (cid:12)(cid:12) (cid:111) , Γ pn → e + ¯ ν = ρ N m N π (cid:110) η (cid:12)(cid:12) C ( pn ) SL (cid:12)(cid:12) , η (cid:12)(cid:12) C ( pn ) S L (cid:12)(cid:12) ,η (cid:12)(cid:12) C ( pn ) VL (cid:12)(cid:12) , η (cid:12)(cid:12) C ( pn ) V L (cid:12)(cid:12) , η (cid:12)(cid:12) C ( pn ) T (cid:12)(cid:12) (cid:111) , (12)where again we neglect the interference terms betweenany pair of Wilson coeﬃcients under the assumption ofone operator dominant each time. η i are dimensionlessparameters and deﬁned as η = 0 , η = (cid:18) − δ (cid:19) = 1 ,η = δ (cid:18) − δ (cid:19) = 7 . × − ,η = (cid:18) δ (cid:19) (cid:18) − δ (cid:19) = 2 ,η = 4 (cid:18) δ (cid:19) (cid:18) − δ (cid:19) = 4 , (13)3here δ = m e /m N and the nucleon velocity eﬀect is ne-glected. If we keep the small nucleon velocity eﬀect, η , will be modiﬁed to be proportional to the squared veloc-ity, i.e., (cid:101) η , ∼ v N ∼ E/m N ∼ − . The nucleonvelocity v N is estimated to be around 0 . c by taking theaverage nucleon binding energy ∆ E ∼ c is the speed of light. For η , , , thenucleon velocity eﬀect is small relative to their leadingcontribution and can be safely neglected. In the follow-ing, we use the modiﬁed value (cid:101) η ∼ v N ∼ − for η butstill keep η ’s value as a conservative estimation for thebound on partial lifetime of n → ¯ pe + χ .In order to employ the experimental result for np → e + χ , we assume each time there is only one operatordominant. According to the decay rate in Eq. (7) for n → ¯ pe + χ and Eq. (12) for np → e + χ , after keepingone term each time and eliminating the common Wilsoncoeﬃcient, we obtain the following relationsΓ − n → ¯ pe + χ = 2 π m n ρ N m N π η i (cid:15) i Γ − pn → e + χ = 2 π ρ N m n η i (cid:15) i Γ − pn → e + χ (cid:38) π ρ N m n η i (cid:15) i Γ − , SK pn → e + χ , (14)where in the last step we take the SK lower bound onthe partial lifetime of pn → e + χ as our input, in whichΓ − , SK pn → e + χ (cid:38) . × yrs [12]. Then our main resultfor the prediction of the | ∆ B | = 2 neutron decay mode n → ¯ pe + χ is as followsΓ − n → ¯ pe + χ (cid:38) . × yrs , for (cid:101) O ( pn ) SR , O ( pn ) SL , (15)Γ − n → ¯ pe + χ (cid:38) . × yrs , for (cid:101) O ( pn ) S R , O ( pn ) S L , (16)Γ − n → ¯ pe + χ (cid:38) . × yrs , for (cid:101) O ( pn ) VL , O ( pn ) VL , (17)Γ − n → ¯ pe + χ (cid:38) . × yrs , for (cid:101) O ( pn ) V L , O ( pn ) V L , (18)Γ − n → ¯ pe + χ (cid:38) . × yrs , for (cid:101) O ( pn ) TR , O ( pn ) T , (19)we see that the extrapolated lower bound on the partiallifetime varies from O (10 ) to O (10 ) years dependingon which operator is dominant. Furthermore, if we take (cid:101) η ∼ − as did for (cid:101) η , the lower bound is improvedto be O (10 ) years. As a comparison with the exper-imental bound on the inclusive mode Γ − , PDG n → e + anything (cid:38) . × yrs quoted by the Particle Data Group [30],we ﬁnd the indirect bound obtained here are improvedby at least ten orders of magnitude.Even though we obtained the above bound on the SMneutrino case, the bound is also valid for beyond the SM light fermion with a mass m χ < m n − m p − m e , theanalysis is similar to the above SM neutrino case. If themass m χ (cid:54) = 0, we would expect a stronger bound on thepartial lifetime due to an even smaller phase space.We can also use the known constraints on dinucleonto dilepton transition rate to set bound on NP scale.Applying the SK experimental limits on the transitions( pp → (cid:96) + (cid:96) (cid:48) + , pn → (cid:96) + ¯ ν (cid:48) , nn → ¯ ν ¯ ν (cid:48) ) from our exam-ple operator O S, ( S ) Q L [12, 14], we ﬁnd the scale of NPΛ NP (cid:38) − pp → (cid:96) + (cid:96) (cid:48) + + 4 jets or the futureLHeC via e − p → (cid:96) + + 5 jets. Discussions.

From our analysis in this paper, we ex-pect the partial lifetime for the unique | ∆ B | = 2 neutrondecay mode n → ¯ pe + χ is far longer than the currentexperimental sensitivity. Nevertheless, this decay modehas several distinctive features in experimental searchesand also strong astrophysical implications. First, in thebound nuclei, if a neutron decays into an antiproton,then the antiproton could annihilate with a neighbourproton to release two bunch of energetic mesons in theopposite direction through the QCD interactions, or two O (1 GeV) gamma photons via the QED interaction, suchsignals are clear evidence of this decay mode if they couldhappen. They may be searched for in the future neu-trino experiments like the Super-K, DUNE [31], JUNO,etc. Second, unlike the | ∆ B | = 2 dinucleon to dime-son/dilepton decays [12–14], this decay mode can takeplace for the free neutron and it opens up the possibilityto search this | ∆ B | = 2 signal from free neutron decayand oscillation experiments like the European SpallationSource [32]. This neutron decay mode could also haveimpact on the astrophysical processes like the cooling ofneutron star and the evolution of the universe, which canbe used to constrain the relevant NP complementarily. ACKNOWLEDGEMENT

This work was supported in part by NSFC (Grants11735010, 11975149, 12090064), by Key Laboratory forParticle Physics, Astrophysics and Cosmology, Ministryof Education, and Shanghai Key Laboratory for Parti-cle Physics and Cosmology (Grant No. 15DZ2272100),and in part by the MOST (Grants No.109-2112-M-002-017-MY3 and 109-2811-M-002-535). XDM would liketo thank Fu-Sheng Yu for his invitation as a visitor atLanzhou Uni. due to the inﬂuence of the COVID-19 ontraveling back to Taiwan, and also the valuable discus-sions with him on this work. [1] A. Sakharov, Sov. Phys. Usp. , 392 (1991). [2] H. Georgi and S. Glashow, Phys. Rev. Lett. , 438(1974).

3] K. Babu and S. M. Barr, Phys. Rev. D , 5354 (1993),arXiv:hep-ph/9306242.[4] P. Nath and P. Fileviez Perez, Phys. Rept. , 191(2007), arXiv:hep-ph/0601023.[5] S. Nussinov and R. Shrock, Phys. Rev. Lett. , 171601(2002), arXiv:hep-ph/0112337.[6] J. M. Arnold, B. Fornal, and M. B. Wise, Phys. Rev. D , 075004 (2013), arXiv:1212.4556 [hep-ph].[7] P. B. Dev and R. N. Mohapatra, Phys. Rev. D , 016007(2015), arXiv:1504.07196 [hep-ph].[8] S. Gardner and X. Yan, Phys. Lett. B , 421 (2019),arXiv:1808.05288 [hep-ph].[9] S. Girmohanta and R. Shrock, Phys. Rev. D , 015017(2020), arXiv:1911.05102 [hep-ph].[10] M. Baldo-Ceolin et al. , Z. Phys. C , 409 (1994).[11] K. Abe et al. (Super-Kamiokande), (2020),arXiv:2012.02607 [hep-ex].[12] V. Takhistov et al. (Super-Kamiokande), Phys. Rev.Lett. , 121803 (2015), arXiv:1508.05530 [hep-ex].[13] V. Takhistov (Super-Kamiokande), in (2016) pp. 437–444, arXiv:1605.03235 [hep-ex].[14] S. Sussman et al. (Super-Kamiokande), (2018),arXiv:1811.12430 [hep-ex].[15] J. Bramante, J. Kumar, and J. Learned, Phys. Rev. D , 035012 (2015), arXiv:1412.2140 [hep-ph].[16] S. Rao and R. E. Shrock, Nucl. Phys. B , 143 (1984).[17] W. E. Caswell, J. Milutinovic, and G. Senjanovic, Phys.Lett. B , 373 (1983).[18] S. Gardner and X. Yan, Phys. Rev. D , 056008 (2018),arXiv:1710.09292 [hep-ph]. [19] A. Kobach, Phys. Lett. B , 455 (2016),arXiv:1604.05726 [hep-ph].[20] E. E. Jenkins, A. V. Manohar, and P. Stoﬀer, JHEP ,016 (2018), arXiv:1709.04486 [hep-ph].[21] Y. Liao, X.-D. Ma, and Q.-Y. Wang, JHEP , 162(2020), arXiv:2005.08013 [hep-ph].[22] E. E. Jenkins and A. V. Manohar, Phys. Lett. B ,558 (1991).[23] V. Bernard, N. Kaiser, and U.-G. Meissner, Int. J. Mod.Phys. E , 193 (1995), arXiv:hep-ph/9501384.[24] J. Bijnens and E. Kofoed, Eur. Phys. J. C , 867 (2017),arXiv:1710.04383 [hep-ph].[25] X.-G. He, X.-D. Ma, and G. Valencia, .[26] E. Rinaldi, S. Syritsyn, M. L. Wagman, M. I. Buchoﬀ,C. Schroeder, and J. Wasem, Phys. Rev. D , 074510(2019), arXiv:1901.07519 [hep-lat].[27] S. Weinberg, Phys. Rev. Lett. , 2333 (1989).[28] J. Goity and M. Sher, Phys. Lett. B , 69 (1995),[Erratum: Phys.Lett.B 385, 500 (1996)], arXiv:hep-ph/9412208.[29] H. Nishino et al. (Super-Kamiokande), Phys. Rev. D ,112001 (2012), arXiv:1203.4030 [hep-ex].[30] P. Zyla et al. (Particle Data Group), PTEP ,083C01 (2020).[31] B. Abi et al. (DUNE), (2020), arXiv:2008.12769 [hep-ex].[32] A. Addazi et al. , (2020), arXiv:2006.04907 [physics.ins-det]., (2020), arXiv:2006.04907 [physics.ins-det].