Probing nonstandard lepton number violating interactions in neutrino oscillations

We discuss lepton number violating processes in the context of long-baseline neutrino oscillations. We summarise and compare neutrino flavour oscillations in quantum mechanics and quantum field theory, both for standard oscillations and for those that violate lepton number. When the active neutrinos are Majorana in nature, the required helicity reversal gives a strong suppression by the neutrino mass over the energy, $(m_{\nu}/E_{\nu})^{2}$. Instead, the presence of non-standard lepton number violating interactions incorporating right-handed lepton currents at production or detection alleviate the mass suppression while also factorising the oscillation probability from the total rate. Such interactions arise from dimension-six operators in the low energy effective field theory of the Standard Model. We derive general and simplified expressions for the lepton number violating oscillation probabilities and use limits from MINOS and KamLAND to place bounds on the interaction strength in interplay with the unknown Majorana phases in neutrino mixing. We compare the bounds with those from neutrinoless double beta decay and other microscopic lepton number violating processes and outline the requirements for future short- and long-baseline neutrino oscillation experiments to improve on the existing bounds.


I. INTRODUCTION
The seminal confirmation of neutrino flavour oscillations by the Super-Kamiokande and SNO experiments in 1998 and 2001, respectively, initiated a golden era in the experimental and theoretical studies of massive neutrinos in the Standard Model (SM) [1,2]. Five of the six (or eight) independent mixing parameters describing the three Dirac (or Majorana) neutrinos have now been pinned down by solar, atmospheric, accelerator and reactor neutrino oscillation experiments [3].
There are still a few pieces of the neutrino mass puzzle that remain unknown. First, the value of the Dirac phase δ which controls the magnitude of CP-violation in the neutrino sector. Second, the sign of the atmospheric squared mass difference ∆m 2 23 which will decipher the normal ordering (NO) or inverted ordering (IO) of the neutrinos. Both quantities will be determined in the coming years by accelerator and reactor oscillation experiments such as NOνA, DUNE, Hyper-Kamiokande and JUNO [4][5][6][7]. Neutrino oscillations are insensitive to the absolute scale of the neutrino masses m 0beta decay experiments such as KATRIN aim to combine their results with precision measurements of the cosmic microwave background and large-scale galaxy clustering to place a stringent upper bound on this scale [8][9][10][11].
Last but by no means the least is the question of the fundamental nature of the light neutrinos.
For massive neutrinos -or generally any massive fermion that is a singlet under gauge transformations -it is possible to write two Lorentz invariant mass terms. The first, proposed by Dirac, is the same as that of the charged fermions. The second, proposed by Majorana, violates lepton number by two units [12,13]. On the grounds of naturalness the latter has been favoured for some time, with effective models using the high scale of new physics (NP) to generate the tiny neutrino masses -also known as the seesaw mechanism. Common UV-complete versions of this mechanism incorporate heavy Majorana right-handed (RH) neutrinos into the larger gauge structure of a grand unified theory [14][15][16][17].
While these theories generically predict lepton number violating (LNV) phenomena, none have yet been observed in nature. If the light neutrinos are indeed Majorana, a standard explanation of this non-observation is to promote the non-anomalous combination of baryon number (B) minus lepton number (L) to a gauge group U (1) B−L . This symmetry is spontaneously broken at low energies and thereby suppresses (B − L)-violating processes by the high symmetry-breaking scale [18]. Other models which predict additional particles such as the Majoron and supersymmetric partners have been excluded at high significance by astronomical and collider observations [19][20][21][22][23].
If the neutrinos are instead purely Dirac, the absence of an LNV signal could be explained if B − L remains a global symmetry after the spontaneous symmetry breaking of a left-right (LR) symmetric model [24]. There is also the intriguing possibility that neutrinos are quasi-Dirac [25][26][27].
The most promising current effort to detect an LNV process is via searches for neutrinoless double beta (0νββ) decay which is the nuclear decay process (Z, A) → (Z +2, A)+2e − . If observed, the absence of outgoing neutrinos implies there to be an internal neutrino propagator, which is possible only if the neutrino is its own antiparticle, i.e. if neutrinos are Majorana fermions. Even if unrelated new physics (NP) were to trigger 0νββ decay, a positive signal would automatically imply Majorana neutrinos via the Schechter-Valle black-box theorem [28]. Current and future 0νββ experiments such as EXO-200, KamLAND-Zen, CUORE, GERDA-II and SuperNEMO are aiming to push the lower limit on the 0νββ decay half-life T 0νββ 1/2 up to 10 26 years for various nuclear isotopes [29]. This lower limit can be converted to an upper limit on the neutrino mass to cancel the different contributions to m ββ [30].
Alternatively one can introduce Wilson coefficients for the relevant dimension-d combinations of SM fields, suppressed by the scale of NP raised to the power (4−d). A general analysis of LNV channels has been conducted in the effective field theory formalism for operators up to dimension eleven [51]. Experimental signatures of alternative LNV processes have been explored in the literature, for instance the decays K + → π − µ + µ + , τ − → π − π − µ + and the lepton flavour and number violating (LNFV) muon conversion process µ − + (Z, A) → e + + (Z − 2, A) [30,52,53]. While these channels are less sensitive to LNV compared to 0νββ decay, searches have been conducted by a variety of experiments [54][55][56]. The contributions of lepton number conserving (LNC) but potentially lepton flavour violating (LFV) NSI to standard processes such as meson decays, neutrino oscillations and neutrino scattering have also been considered exhaustively in the literature [50,. For processes which are insensitive to lepton number -for example, the outgoing neutrino is not detected in charged pion decays π ± → ± α (−) ν α -precision probes of SM predictions such as lepton universality and the unitarity of the Cabibbo-Kobayashi-Maskawa (CKM) matrix can be used to constrain both LNC and LNV NSI [64].
Much of the literature has so far only considered the effect of LNC NSI on neutrino oscillations. This is for good reason -neutrino oscillation experiments are typically only concerned with neutrino flavour at production and detection, inferring the process ν α → ν β (orν α →ν β ) from the accompanying charged lepton + α ( − α ) at production and − β ( + β ) at detection. In many cases there is no detector at the neutrino source to identify the initial composition of flavours -this and the associated energy distribution must be inferred from separate measurements and Monte Carlo simulations [83][84][85]. Often there is also no way to determine the charge of the outgoing lepton at the far detector and therefore to discern the incoming lepton as a neutrino or antineutrino. This probe of lepton number is not a priority for most oscillation experiments because 'ν α →ν β ' is heavily suppressed if the mechanism is the standard Majorana neutrino helicity reversal [86][87][88][89]. Like any other helicity suppression this introduces a factor of ∼ (m ν /E ν ) 2 to the rate of the process, where m ν is the neutrino mass scale and E ν is the neutrino energy. Hence the small neutrino masses (m ν ∼ 0.1 eV) and typical large oscillation experiment neutrino energies (E ν ∼ 5 MeV − 2 GeV) combine to suppress the magnitude of the 'ν α →ν β ' and 'ν α → ν β ' processes by ∼ 10 −21 − 10 −16 .
On the other hand, the highly suppressed amplitude of such a process can be exploited -any positive signal of LNV such as an excess of 'wrong'-signed charged leptons at the far detector would then strongly imply NP. Experiments that have been sensitive to the charge of outgoing leptons, such as the long-baseline (LBL) oscillation experiment MINOS and the LBL reactor/solar oscillation experiment KamLAND, can thus be used to constrain LNV NSI.
In this paper we begin Sec. II with a brief discussion and derivation of the ν α → ν β oscillation probability in quantum mechanics (QM) and quantum field theory (QFT). Using the latter we will study the 'ν α →ν β ' process for Majorana neutrinos, obtaining the expected ∼ (m ν /E ν ) 2 suppression of the total rate. We also show that the total rate cannot be factorised into an oscillation probability in a similar way to ν α → ν β oscillations. In Sec. III we consider the impact of an LNV NSI -rather than the well-studied LNC NSI -on neutrino oscillations. We will demonstrate that the total rate is no longer suppressed and factorises when the chirality of the production and detection vertices are opposite. We write down a general expression for the non-standard oscillation probability and a simplified expression in the two-neutrino (2ν) mixing approximation, specifically for the ν µ −ν τ sector. This allows us in Sec. IV A to use a limit from the MINOS experiment on the 'ν µ →ν µ ' process to place bounds on the simplified 2ν parameter space of this effective model. We generalise to the three-neutrino (3ν) mixing scheme and re-evaluate constraints from MINOS in Sec. IV B along with those from the KamLAND experiment in Sec.
IV C. We compare these constraints to those from microscopic LNV processes such as 0νββ decay, µ − − e + conversion and radiative neutrino masses in Sec. IV D. To conclude, we summarise our results and briefly outline the potential for future oscillation experiments with similar sensitivity to improve on these bounds.
where P L is the left-handed (LH) chirality projector and we have employed four-component spinor notation for the fields. The kinetic term is diagonal in the basis of definite neutrino mass, labelled by the index i, while the interaction term is diagonal in the basis in which the neutrino and associated charged lepton have the same flavour, labelled by the index α.
The Lagrangian in Eq. (1) is valid for both Dirac and Majorana neutrinos (up to the indicated factor of 1/2 in the Majorana kinetic term). The former and latter are defined by Here ν D has a new RH component ν R , whereas ν M has the charge conjugate of its LH component (up to an arbitrary intrinsic charge parity which is commonly set to unity). For Dirac neutrinos the SM CC interaction term shown in Eq. (1) creates negative helicity neutrinos |ν(q, −) and annihilates positive helicity antineutrinos |ν(q, +) . For Majorana neutrinos |ν(q, +) is equivalent to |ν(q, +) . The creation and annihilation of the other two degrees of freedom in the Dirac case, or the 'wrong' helicity degree of freedom in the Majorana case, are suppressed by ∼ (m ν /E ν ) at the amplitude level [89].
The following discussion is valid in both cases because the interaction Lagrangian remains the same. The fields in the flavour basis can be rotated to those in the mass basis using the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) mixing matrix U , While oscillation experiments produce neutrinos in a particular flavour eigenstate |ν α , for oscillations to take place this must be a coherent superposition of the mass eigenstates. The time evolution of each massive state is governed by the Schrödinger equation, resulting in where 2|q| is the energy of each ultra-relativistic massive neutrino with mass m i and shared three-momentum q. An oscillation probability can be derived by computing the square of the overlap between the time-evolved initial flavour state |ν α (T ) and an arbitrary final flavour state |ν β , where we have taken T L and L is the oscillation baseline.
This simple derivation was the first attempt to quantify oscillation phenomena in a quantum mechanical framework [90][91][92]. It was soon realised, however, that this description relies on numerous unphysical assumptions, namely the use of plane wave states. Assuming the propagating |ν i states to be plane waves of equal momenta q forces the external particles at production to have definite energies and momenta. Energy-momentum conservation at production is then in tension with the creation of three |ν i with different energies E q = |q| 2 + m 2 i [93]. While it is possible to derive Eq. (5) without the equal momentum assumption, an overall uncertainty in the energy and momentum of the neutrino mass eigenstates is still a necessary component [94,95].
Any rigorous treatment of oscillation in QM must therefore describe the |ν i states with wave packets, which is also known as the internal wave packet model [96][97][98]. As in Eq. (5), the oscillation probability is proportional to the overlap of the wave packets. The loss of coherence seen at long distance in oscillation experiments qualitatively arises from dispersion of the wave packets, which propagate at different group velocities v = ∂Eq ∂q q=Q with Q the mean momentum [99][100][101][102].
Despite these improvements there are still fundamental issues with this picture. Firstly, it is not apparent what shape the neutrino wave packets should take. Secondly, a finite energy-momentum uncertainty requires the localisation in space and time of the production and detection interaction processes, which is not accounted for [101]. It is also difficult to define a meaningful Fock space for the flavour eigenstate neutrinos [103]. Finally, the QM approach does not consider the possible entanglement between the outgoing ν i and ± α at production [104]. A fully consistent framework in which to study neutrino oscillations is QFT, also known as the external wave packet model [105]. In this formalism the entire production, propagation and detection process can be described by a macroscopic Feynman diagram, as shown in Fig. 1. The external interacting particles are wave packets centred on the production and detection points x P and x D , while the neutrinos are described by an internal propagator. The neutrino wave packets used in QM can be derived from the external particle wave packets in QFT; in QFT, however, the coherence conditions at each stage in the process are explicit [102,106,107].
To calculate an oscillation probability in QFT we must first compute an overall rate. To do this one constructs the S-matrix element at first order in the Fermi coupling G F / where P I , D I , P F and D F are the initial and final state particles at production and detection, respectively, and x 1 and x 2 are space-time points in the vicinity of x P and x D . This tree level process is depicted to the left of Fig. 1. The external asymptotic states ψ ∈ {P I , D I , P F , D F } are now described by the wave packets where f (ψ) P ψ (p) is the momentum distribution function of an external particle ψ with mean momentum P ψ . To create and annihilate an internal neutrino the production and detection Lagrangian terms must take the general form where L P (x) and L D (x) are the 'reduced' Lagrangian interaction terms for production and detection with the neutrino fields and PMNS matrix elements removed.
The total rate is proportional to the spin-averaged S-matrix element squared, which can be expanded as where Tr denotes the Dirac trace. The A i factors are where Φ P and Φ D are integrals quantifying the overlap of external wave packets at production and detection respectively, explicitly written as where x 1 = x 1 − x P , x 2 = x 2 − x D and M P , M D are 'reduced' matrix elements defined by The trace appearing in the interference term of Eq. (10) is now where the underline of Φ P,D is shorthand for Φ † P,D γ 0 . We now require the well-known limit of the d 3 q integrals as L → ∞, first given in Ref. [108]. If ψ(q) is a twice differentiable function, for large For A i < 0 the integral falls as L −2 and can be neglected. Applying Eq. (15) to Eq. (14) with which has effectively set the neutrinos to be on-shell. As has been stated before in the literature [109,110], if the production and detection processes are of the same chirality, the trace on the right-hand side of Eq. (16) can be factorised in the relativistic limit (m i ≈ 0) as Furthermore, writing ( q + m) as a spinor sum and expanding |q| = we can see that for ν α → ν β the propagation of positive helicity neutrinos is doubly suppressed by (m i /2E q ) 2 compared to the propagation of negative helicity neutrinos.
Neglecting the positive helicity neutrino contribution to the spinor sum in Eq. (18), the negative helicity spinors can be absorbed into the overlap integrals on either side in the traces of Eq. (17): where in the relativistic limit the mass eigenstate indices can be neglected [102]. The trace appearing in the interference term can now be expressed as and we see that the contributions from production, propagation and detection have factorised at the squared amplitude level. The total rate for the ν α → ν β process is now related to the differential production flux, oscillation probability and detection cross section as where we have neglected experimental considerations such as the detector efficiency and fiducial volume [111]. The oscillation probability can now be solved for through rearrangement of Eq.
(21). As Ref. [102] shows, in the case of continuous fluxes of incoming particles, the differential production flux and detection cross section have the proportionality Using that |Φ P | 2 and |Φ D | 2 are independent of the mass m i , while also taking |q i | ≈ |q j | in the relativistic or quasi-degenerate mass limit (equivalent to the inequality |q i | − |q j | |q i |, |q j |), any dependence P να→ν β has on the specific form of the wave packets cancels [102]. In this limit P να→ν β is given by Eq. (5), confirming the result of the naive QM approach under the above conditions.
For Majorana neutrinos the process 'ν α →ν β ' is possible and suppressed by (m i /2E q ) 2 . To see this one can construct an amplitude like Eq. (6) where the production and detection Lagrangian terms are the same which is non-zero for Majorana neutrinos becauseν i (x) both creates and annihilates |ν i . The Feynman diagram for this process is depicted in Fig. 1 (right). The amplitude can again be squared and expressed as a sum of amplitudes A i . Using the Majorana fermion Feynman rules of Ref. [112], the amplitudes A i are identical to Eq.
where iΓ is the vertex factor for the detection process and C is the charge conjugation matrix. For jD appearing in the squared amplitude now contains the factor squared, and so is singly suppressed by (m i /2E q ) 2 compared to the propagation of negative helicity neutrinos in Eq. (18) [113][114][115].
We also see that in the relativistic limit the trace vanishes instead of factorising into components corresponding to production, oscillation, and detection. Strictly speaking it is therefore impossible to define an oscillation probability if helicity reversal is the dominant mechanism contributing to 'ν α →ν β ', only a probability for the entire process [109].

III. NON-STANDARD LEPTON NUMBER VIOLATING INTERACTIONS
We will now consider interactions at production and detection which are different from the usual SM CC interaction. CC-like and neutral current (NC)-like NSI which conserve lepton number, but introduce a new source of LFV to the SM have been studied extensively in the literature; see Refs.
If an experiment was indeed sensitive to the sign of the charged lepton ± β produced in the CC-like interaction of the incoming neutrino at detection, it would be able to distinguish between neutrinos and antineutrinos, or in other words negative and positive helicity Majorana neutrinos.
This has been possible in the past -a magnetised far detector was used to determine the charge from the curvature of tracks in the steel scintillator near and far detectors of MINOS. The prompt energy deposit from a + β and neutron capture on hydrogen, measurable for KamLAND, was also used as a distinct signal from the − β case [134]. It is therefore worth discussing the LNV equivalents of the CC-like and NC-like NSI discussed previously, as the non-observation of an excess of 'wrong'signed charged leptons by these experiments is able to say something about the possible size of an LNV equivalent ε coefficient.
The operators to the left of Table I trigger the standard LNC NSI studied in the literature, regardless of neutrinos being Dirac or Majorana. These are given in the notation of Refs. [133,135], which enumerate the effective operators of dimension-d ≤ 6 generated below the electroweak scale.  Table I are matched to ∆L = 0 operators in the 59 operator basis of the SM effective field theory (SMEFT).
The LNV operators to the right of Table I can also trigger an interaction at detection for an incoming neutrino. It is not possible to match these operators to ∆L = 0 dimension-six SMEFT operators above the electroweak scale, but instead to dimension-seven |∆L| = 2 operators with an additional Higgs field H to conserve hypercharge U (1) Y . These dimension-seven operators have been considered previously along with the complete list of LNV odd-dimensional effective operators up to dimension eleven in Ref. [51]. The operators in Table I correspond  also been considered more specifically in Ref. [136], but parametrised there in terms of a Wilson where Λ NP is the scale of NP. In this paper we will adopt the following notation. A general Lagrangian, normalised to the Fermi interaction, and which can trigger an LNV interaction at production or detection, can be written as where the second term includes all possible Lorentz contractions of the leptonic current j (ρ) = (νO (ρ) ) and hadronic current J (σ) = (ūO (σ) d), where (ρ, σ) run over the usual scalar (S), pseudoscalar (P ), vector (V ), axial vector (A) and tensor (T ) Dirac structures. We choose linear combinations of these proportional to the chirality projection operators P L and P R .
The ε (ρ,σ) coefficients parametrise the strength of the NSI compared to G F / √ 2. We can define these coefficients as matrices in the flavour basis (ε (ρ,σ) ) or mass basis (γ (ρ,σ) ) of the neutrino field, related by the rotation where U is the usual PMNS matrix. It is now straightforward to see the effect of these interactions on oscillations. If the production and detection leptonic currents are of the same chirality, ν α → ν β will be unsuppressed, while for Majorana neutrinos 'ν α →ν β ' will be suppressed by (m i /2E q ) 2 . If the production and detection are of opposite chirality, 'ν α →ν β ' will now instead be factorisable and suppressed by |ε (ρ,σ) | 2 . A commonly used example is a V − A ('L') leptonic current at production and a V + A ('R') leptonic current at detection arising from a LR symmetric model with an additional broken SU (2) R gauge symmetry [137].
An LBL oscillation experiment sensitive to the sign of the outgoing lepton at the far detector, while not detecting the ± α at the production process (occurring in the beam pipe), could measure the ratio where N ± β is the number of detected ± β at the far detector. If we assume NP to be the cause of the first term in the numerator of Eq. (27), and that this total rate is factorisable, the terms in R αβ can be decomposed as where the (ρ, σ) superscript denotes a NP effect. All probabilities are given this superscript because the process 'ν α →ν β ' reduces the number of neutrinos that undergo ν α → ν β , and thus the sum β P (ρ,σ) να→ν β does not equal unity. The standard oscillation probability must be scaled by a normalisation factor depending on the NP and that ensures β (P (ρ,σ) where the normalisation factor is Expanding the denominator in Eq. (29) it is straightforward to see that β P (ρ,σ) and β P (ρ,σ) , so taking ε → 0 recovers the SM prediction. We assume that the NP is a small effect, ε (ρ,σ) 1, and thus neglect this modification to the probabilities, i.e. set N (ρ,σ) α ≈ 1. The factor N (ρ,σ) α cancels in the ratio R αβ regardless. As another simplification we take the NSI coefficients to be real in this work.
Expanding the probability to the right of Eq. (29) gives where we have rotated the γ In the 3ν mixing scheme these are complicated functions of the three mixing angles, three squared splittings and two Majorana phases. If one were to consider atmospheric or accelerator 'ν µ →ν µ,τ ' oscillations, the 2ν mixing approximation is valid because the oscillation Hamiltonian is dominated by the atmospheric mass splitting. The mixing matrix in this case is where ϑ is the single mixing angle (approximately corresponding to θ 23 ) and η is a Majorana phase [138], while the squared mass splitting is δm 2 (corresponding to ∆m 2 23 ). F λλ take on the simplified forms, where ϕ = η − δm 2 L 4Eq . We will make use of these functions in the next section when studying the results from MINOS in the 2ν mixing approximation.

INTERACTIONS
We will now use the parametrisation discussed above to put constraints on the ε coefficient parameter space of these LNV NSI. We will first discuss constraints from the MINOS experiment in the 2ν mixing approximation in Sec. IV A, moving onto analysis of both MINOS and KamLAND in the 3ν mixing scheme in Secs. IV B and C. In Sec. IV D we will compare these constraints to the more common limits from microscopic LNV processes such as 0νββ decay, µ − − e + conversion, rare meson decays and the radiative generation of neutrino masses. We summarise all limits in Table III.

MIXING APPROXIMATION
The MINOS experiment initially took data from 2005 to 2012 and used the low energy NuMI beam to detect neutrinos with a near detector at Fermilab and a far detector at a baseline of L = 735 km from the source at the Soudan mine [139]. After a break the experiment continued from 2013 to 2016 as MINOS+, using the medium energy NuMI beam [140]. The experiment observed the disappearance of ν µ produced from π + decays (in the focusing beam configuration) andν µ from π − decays (defocusing), allowing the atmospheric mixing parameters dominating ν µ → ν µ disappearance to be probed. The experiment also confirmed ν e andν e appearance, constraining the reactor mixing angle θ 13 . Most importantly for our discussion, charged lepton sign identification was possible in the near and far detectors through the use of 1.3 T toroidal magnetic fields -ν µ ,ν µ , ν e andν e events could therefore be distinguished from the curvature of the outgoing µ − , µ + , e − and e + tracks, respectively.
Before MINOS began taking data, the expected fluxes of ν µ andν µ in the focusing and defocusing configurations were predicted to high precision from hadron production data and in situ measurements. These were more recently updated in Ref. [85]. In the focusing configuration the background ofν e produced from pion decays upstream of the target and avoiding deflection by the magnetic field is non-negligible and an important systematic error to correct [141]. There are alsō ν e produced downstream from secondary interactions in the beam pipe wall [142].
The ratio in Eq. (27) can be split into a signal part S µµ arising from the 'ν µ →ν µ ' process and a background part B µµ arising from the standard oscillation of the background antineutrinos ν µ →ν µ . The MINOS analysis of Ref. [142] removes the predicted energy-dependent value of B µµ from the total measured R µµ and derives the constraint S µµ 0.026. From we can put corresponding constraints on the range of possible ε (ρ,σ) µλ . The factorised form of Eq. (35) of course assumes the chirality of the production and detection processes to be opposite. For example, we could consider a V − A leptonic current at production (ρ = σ = L) and a V + A leptonic current at detection (ρ = R). The two possible Lorentz contractions with the hadronic current are σ = R and σ = L; in a LR symmetric model the former corresponds to the exchange of a W R boson, the latter to W L − W R mixing if the boson mass eigenstates are mismatched, as depicted in Fig. 2 [89]. To simplify this work we will only consider these two cases, setting ε (ρ,σ) µα for all other (ρ, σ) to zero and cleaning up the notation with ε noting that there is a subtle difference between the LNV NSI being at production and detection because the outgoing lepton ± α at production is not measured. Hence in theory one must sum over the different initial flavours as α S αµ 0.026, and MINOS is sensitive to ε eλ and ε τ λ . However, for the π ± energies in the NuMI beam it is kinematically impossible to produce a τ ± , and similar to standard π ± decays the production of an electron is chirality suppressed to that of a muon. We therefore neglect this subtlety and assume that an NSI at production is probed in a similar way to that at detection.
Working with the 2ν mixing expansion of P νµ→νµ given in Eqs. (31) and (34) be used to put a constraint on the (ε µµ , ε µτ ) parameter space. To do this we integrate the NuMI are assumed to be equal in the quasi-elastic scattering limit [143].
In Fig. 3 (left) we first plot the allowed regions in the (ε µµ , ε µτ ) parameter space for fixed L/E q = 735/3 km GeV −1 , using best fit values for δm 2 ≈ ∆m 2 23 , ϑ ≈ θ 23 (in the NO scheme), G F , V ud , g V , g A and four different values of the Majorana phase η [30]. This is equivalent to assuming the ν µ flux to be sharply peaked at 3 GeV and evaluating the oscillation probability and cross section at this energy. Appreciable constraints are possible for η = 0 and π, and are marginally better in the ε µτ direction. They are of order |ε µµ | 0.2 and |ε µτ | 0.1. When η = (n + 1/2)π for n ∈ Z a specific direction in the parameter space alleviates the constraints. The narrow ellipse arises because F (µ) τ 1 for the best fit parameters and these particular values of η. In Fig. 3 (right) we depict the allowed regions after the full numerical integration of the numerator and denominator of S µµ . For η = π/4 and π/2 the narrow bands of allowed values are reduced to ellipses more similar to the ellipses at η = 0 and π. The orientations of the ellipses also change marginally. Upper bounds are in the ranges |ε µµ | 0.2 − 0.5 and |ε µµ | 0.2 − 0.6.
While we have so far restricted our analysis to the MINOS experiment, it is interesting to consider what constraints could be made at different baselines for an experiment similar in design to MINOS. To the left of Fig. 4 we set η = 0 and examine the allowed regions in the ε µµε µτ space for different values of the baseline L, derived for illustrative purposes from the MINOS limit S µµ 0.026. We see that at zero distance this limit bounds |ε µµ | 0.16, while ε µτ remains unbounded. This is clear from the expansion of P νµ→νµ -the functions F  For η = 0 (and η = nπ, n ∈ Z) the oscillating part of P να→ν β exactly cancels the oscillating part of P να→ν β , and thus the bound on ε µµ does not depend on L. For η = π/2 the constraint at zero distance is far less stringent, but decreases appreciably from 0 km up to 1000 km. For η = π/4 the constraint worsens as L reaches ∼ 800 km but improves at larger baselines. For L 2000 km the constraints for non-zero η slowly oscillate but are roughly the same as for η = 0, i.e. |ε µµ | 0.15.
We show in Fig. 5 (right) a similar plot for ε µτ , setting ε µµ = 0 and shading the allowed regions as a function of the baseline. At zero baseline ε µτ is unbounded for η = 0, as discussed previously.
For large baseline the upper limits converge to |ε µµ | 0.16.
We summarise the constraints made in the 2ν mixing approximation on the two coefficients ε µµ and ε µτ in Table II. Here we allow one coefficient at a time to be non-zero, computing an upper bound in the fixed energy approximation (middle) and integrating over the energy (right).
The lower and upper values are the best and worst upper bounds, respectively, as the phase η is varied. One can see that ε µτ is unbounded for a specific value of η when the energy is fixed. In this analysis we have of course taken the best fit values of the standard mixing parameters to be fixed.
A rigorous fit to the data would need to let these parameters vary along with the ε coefficients, as taken into account for LNC NSI in Refs. [61,69,81,82].

MIXING SCHEME
We now generalise the analysis to the 3ν mixing scheme. We use the standard parametrisation of the PMNS mixing matrix where s ij = sin θ ij , c ij = cos θ ij and (α 2 , α 3 ) are Majorana phases. We now look to probe the 3 × 3 flavour structure of the NSI coefficient matrix ε βα in which all elements are taken to be real. We again expand the effective non-standard oscillation probability as in Eqs. (31) and (32), Comparing the numerical expression of F (µ) µ in the 3ν scheme to that in Eq. (34) for the 2ν scheme and taking the limits ∆m 2 12 → 0, ∆m 2 13 → ∆m 2 23 , we find that ϑ ≈ θ 23 , δm 2 ≈ ∆m 2 23 which we had already assumed in Sec. IV A, and that η ≈ (α 3 − α 2 )/2.
We now return to the interpretation of the MINOS bound S µµ 0.026. This ratio is defined in the same manner as Eq. (35) -multiplying Eq. (38) by the differential production flux and detection cross section and integrating over the NuMI beam energy between 0 and 20 GeV. Now in the 3ν mixing scheme, S µµ < 0.026 can be converted into an allowed region in the (ε µe , ε µµ , ε µτ ) parameter space which depends on the value of the Majorana phases α 2 and α 3 . Firstly, and in order to compare with bounds in the 2ν mixing approximation, we set ε µe = 0 and depict in Fig. 6 (left) the allowed regions in the (ε µµ , ε µτ ) parameter space for α 2 = 0 and three different values of α 3 . The ellipses are of similar size to those for the 2ν mixing scheme but have shapes and orientations (which we define as the ellipse eccentricity e and anticlockwise angle Θ from the positive ε µµ axis to the semi-major axis, respectively) with similar dependences on the Majorana phases. We show these dependences as contour plots in Fig. 6 (middle and right). We can see that the dependence of the angle Θ is approximately the same as η ≈ (α 3 − α 2 )/2, which is evident from the diagonal lines of roughly constant Θ along α 3 = α 2 + C. For α 2 = α 3 = 0 for example, we see that the constraint is better in the ε µµ direction -this is similar to η = 0 in the 2ν scheme.
We see that the largest eccentricity occurs at α 2 ≈ α 3 ≈ π -this coincides with the semi-major axis pointing in the ε µτ direction and therefore the upper bound on |ε µτ | can be slightly larger than the upper bound on |ε µµ |.
In Fig. 7  While we have so far considered only MINOS, OPERA was another LBL oscillation experiment to employ a magnetic field in the far detector [144,145]. Unlike MINOS, OPERA measured neutrinos from the CNGS beam at CERN that were above the threshold for τ ± production. The main physics goal of the experiment was to confirm ν τ appearance -around 10 τ ± events were recorded over four years of data taking [146]. Unfortunately the experiment was only able to distinguish the charge of one τ − event at 5σ significance, while the other charges were undetermined.
The statistics are therefore too low to comment on OPERA's sensitivity to lepton number. A future high-statistics LBL experiment above the τ ± threshold like OPERA would therefore be able to probe the τ sector NSI coefficients, ε τ e , ε τ µ and ε τ τ . If the initial ν e are produced from the beta decay of 8 B and propagate from the solar core to the solar surface, and then through the vacuum to the KamLAND detector. The oscillation probability must therefore take into account the resonant conversion of solar ν e to ν µ and ν τ through the Mikheyev-Smirnov-Wolfenstein (MSW) effect, a consequence of the slowly decreasing matter potential from the Sun's core to surface. We make the approximation that the conversion is adiabatic and utilise the ∆m 2 12 ∆m 2 23 hierarchy to write the standard ν α → ν β conversion probability in a similar form to that in Ref. [82], where R 23 and W 13 are Euler rotations making up the standard parametrisation of the PMNS matrix and U is a 2 × 2 unitary matrix satisfying where A CC = 2 √ 2G F E q N e and N e is the electron density in the Sun. In order to construct S ee we require the non-standard oscillation equivalent of Eq.   the scope of this work. The possibility that the NSI occurs at production would also complicate the analysis, becauseν e would experience a different matter effect while propagating through the Sun -we therefore concentrate on the NSI being at detection.
We can safely assume that, by the time the solar neutrinos reach Earth, they make up an incoherent mixture of flavour eigenstates. This effectively washes out any dependence on the Majorana phases, and we can approximate the non-standard oscillation probability as being the effective standard probability P eff νe→ν β (taking into account the MSW effect) multiplied by the NSI coefficient |ε eβ | 2 . This takes into account for example the oscillation ν e → ν µ,τ and then the LNV NSI process ν µ,τ p → e + n at detection. The signal ratio S ee and the limit set on it by KamLAND now becomes We now use Eq. (42) to set limits on NSI coefficients in the e sector, ε ee , ε eµ and ε eτ . Using the 8 B flux predicted by the solar model of Ref. [147], we integrate the differential flux, oscillation probability and low energy cross section (which we assume to scale as E 2 q ) over 0.02 MeV bins in the range 8.3 − 14.8 MeV, for both the numerator and the denominator.
Taking each LNV NSI coefficient to be non-zero at a time, we show the derived upper bounds in Table III Other processes capable of probing LNV NSI coefficients. Left: 0νββ decay. Middle: µ − − e + conversion. Right: kaon decay K + → π − µ + µ + .
of S ee contain different oscillation probabilities from the denominator.

VIOLATING PROCESSES
We will now briefly compare the constraints from the LBL experiments MINOS and KamLAND to those from conventional searches for LNV processes. Neutrinoless double beta (0νββ) decay continues to be the most promising method to verify the Majorana nature of the light neutrinos.
It is a highly useful probe because the black-box theorem ensures that any positive signal of 0νββ decay confirms the Majorana case even when the standard light neutrino exchange process is not the dominant mechanism [28,38,39]. An extensive part of the literature has studied non-standard mechanisms, including the LNV dimension-six operator considered in this work [34,40,45,136].
Here we simply extend this to the 3 × 3 flavour structure of the ε coefficient in order to compare with the MINOS and KamLAND upper bounds.
Using from Ref. [32] the general expression for the 0νββ decay inverse half-life when a RH interaction is present at one of the interaction vertices (as shown in Fig. 8): where C mm = 1.12 × 10 −13 , C γγ = 4.44 × 10 −9 and C mγ = 2.19 × 10 −11 are phase space and nuclear matrix element factors given by [32]. The inverse half-life, in a similar manner to the oscillation probability, can be expanded as where F  Best fit values for the mixing angles θ 12 , θ 13 and θ 23 , squared mass splittings ∆m 2 12 and ∆m 2 23 , Dirac CP phase δ are taken in the NO, with a lightest neutrino mass of m 1 = 0 eV. and lightest neutrino mass. Y λ (ζ, α 2 , α 3 , m 0 ) is the contribution from interference between the light neutrino exchange and the non-standard mechanism towards terms linear in the ε coefficients.
We now set T 0νββ 1/2 > 5.3 ×10 25 y derived from the 76 Ge experiment GERDA-II [148,149]. While KamLAND-Zen set a more stringent limit of T 0νββ 1/2 > 1.07 ×10 26 y with 136 Xe [150], the exact number we use is not crucial for the following constraints and discussion. Setting each ε eλ to be non-zero at a time, Eq. (44) can be solved to find an upper bound on the coefficient as a function of α 2 and α 3 . These are displayed in the contour plots of Fig. 9 for a smallest neutrino mass of m 1 = 0 eV in the NO scheme. The associated 'best' and 'worst' upper bounds are shown in Table   III. It can be seen that regardless of the value of the Majorana phases, the upper bound on |ε ee |, as found previously in the literature, is of order 10 −9 . For very finely tuned values of α 2 and α 3 , however, |ε eµ | and |ε eτ | are unbounded. Comparing all upper bounds in the e sector, we see that 0νββ decay is unequivocally the best method to probe ε ee . For a large portion of the (α 2 , α 3 ) parameter space it is also better at constraining ε eµ and ε eτ . However, for certain fine-tuned values of the phases these coefficients become unbounded and KamLAND can provide a better upper bound.
While 0νββ decay is certainly the most sensitive process to test for LNV as outlined above, it can only probe the e sector, whereas other observables may shed light on other flavour coefficients. An interesting process in this regard is the LNFV conversion of captured muons in nuclei, µ − +(Z, A) → e + + (Z − 2, A). Proposed many years ago by Pontecorvo [151], it has gained recent interest due to the upcoming searches for the LNC but LFV muon conversion µ − + (Z, A) → e − + (Z, A) by the COMET [152] and Mu2e [153] collaborations, which promise an increased experimental sensitivity by four orders of magnitude. While it is doubtful that the current limit R Ti µe 10 −11 [154] on the LNFV mode conversion rate can be improved in a similar fashion due to different background considerations [155], µ − − e + conversion is an important complementary probe to 0νββ decay.
To estimate the sensitivity of the LNFV µ − −e + conversion process on the LNV NSI coefficients considered in this paper, we follow the estimate in [136]. In this approach, and using our notation, the conversion rate is approximated as where the effective parameter ξ µe is defined by The two terms in the first line of Eq. (46) take into account that the LNV NSI can be at the interaction vertex of either the incoming µ − or the outgoing e + (the latter being shown in Fig. 8).
The only difference between the two diagrams is the exchange U * ei γ µi ↔ U * µi γ ei . In Eq. (45), q denotes the momentum scale of the intermediate neutrino in the process, q ≈ 100 MeV, and Q is the energy release of the emitted positron determining the size of the phase space with Q ≈ 15.6 MeV [156]. Due to the process being incoherent and partially going to excited final nuclear states, Q also approximately convolutes the nuclear matrix element of this transition.
Here, it should be emphasized that the relevant nuclear matrix elements have not been calculated in detail and Eq. (45) can only be regarded as a very rough estimate of the order of magnitude of the conversion rate. Nevertheless, using the experimental limit Eq. (45), we estimate a limit on |ξ µe | of order |ξ µe | 10 4 . This is barely stringent enough for the EFT assumptions to be consistent, with the limit corresponding to an effective operator scale Λ = (|ξ µe |G F ) −1/2 ≈ 3 GeV > q. Even with the most optimistic improvement of the future sensitivity to R µe 10 −16 [156], the coefficient |ξ µe | will only be probed at |ξ µe | ∼ 30, corresponding to a scale Λ NP ∼ 50 GeV.
We summarise the constraints on the individual coefficients in the µ sector (when each is taken to be non-zero) in Table III λ (ζ, α 2 , α 3 ), the dependences on the Majorana phases for the µ − − e + upper bounds are identical to those for 0νββ decay in Fig. 9, except being scaled by a factor of ∼ 10 12 . From Eq. (46) it is clear that µ − − e + can also probe the e sector coefficients, but sets bounds larger than 0νββ decay by a factor of 10 12 . However, it is of small interest to note that the dependence on the Majorana phases is the same as that for the MINOS upper bounds because the functions preceding . From the rare LNV meson decays such as K ± → π ∓ µ ± µ ± and B + → D − µ + µ + and rare τ decays such as τ − → π − π − µ + it is also possible to probe the LNV coefficients considered in this work as well as those in the τ sector [52,109,157]. However, at present these processes have sensitivities at or worse than that of µ − − e + conversion. Comparing the constraints in the µ sector in Table III, we see that the constraints from MINOS (and similar future LBL oscillations experiments) are currently far more stringent than microscopic LNV processes.
To conclude this section we will briefly mention the analyses in Refs. [51,136] on the contribution toward the Majorana neutrino mass M αβ from the operators in Table I. If the discussion is shifted back to the Wilson coefficients of the operators, i.e. we set and assuming that the Wilson coefficients C βα ∼ O(1), the minimum sum of neutrino masses m ν > 6 × 10 −2 eV sets a lower bound on Λ NP . For the scalar and tensor Dirac structures the loop mass contains one quark Yukawa coupling, and Λ NP > 1 × 10 6 TeV. For a vector Dirac structure the loop gets a contribution from the charged lepton Yukawa coupling and both quark Yukawa couplings if they are right-handed, giving Λ NP > 6 × 10 3 TeV. This can be decreased significantly to Λ NP ∼ 1 TeV if the couplings to third generation quarks are suppressed by some flavour symmetry [136].

V. CONCLUSIONS
In this paper, we have investigated the effect of lepton number violating non-standard interactions on long-baseline neutrino oscillations. If the light active neutrinos are of Majorana nature the 'ν α →ν β ' process become possible, either through a (m ν /E ν ) 2 suppressed light neutrino helicity reversal or an LNV charged current (CC)-like interaction at production or detection. The Majorana neutrino case is favoured from a model-building perspective and is currently being probed by a wide range of neutrinoless double beta (0νββ) decay experiments.
We first studied the different derivations of neutrino oscillations in quantum mechanics and quantum field theory. The QFT model is a more complete and physically consistent picture, taking into account the coherence of overlapping wave packets at production and detection. In this picture we derived the (m ν /E ν ) 2 suppression associated with the Majorana helicity reversal, showing that the total rate in this case cannot be factorised into a production flux, oscillation probability and detection cross section. In the QFT picture we next studied the LNV interaction at the detection process. It is immediately clear that if the chirality of the production and detection processes are opposite, an 'oscillation probability' P να→ν β can be factorised from the total rate of the process, justifying the use of a simple non-standard oscillation formula. The (m ν /E ν ) 2 suppression is replaced by a |ε| 2 suppression, where ε parametrises the strength of the LNV NSI compared to the Fermi coupling G F .
Using a bound made by the MINOS experiment on 'ν µ →ν µ ' oscillations, we put limits on the ε flavour coefficients in the case of a LH SM interaction at production and a RH leptonic current (connected in turn to a LH or RH hadronic current) at detection. The limits are also valid for the LNV NSI being at production. In the two-neutrino mixing scheme which is approximately valid for the ν µ − ν τ sector, we derived a simplified expression for the non-standard oscillation probability in terms of the coefficients ε µµ and ε µτ . Multiplying the non-standard oscillation probability by the NuMI beam flux and quasi-elastic CC cross section, and integrating over the energy range 0 − 20 GeV, we derived upper bounds on the absolute values of the two coefficients. While the value of the single Majorana CPV phase η alters the upper bounds somewhat, we conservatively obtained |ε µµ | 0.6 and |ε µτ | 0.7. If a future experiment like MINOS were to be at a smaller baseline we found that the value of η has a larger impact on the upper bound.
We subsequently generalised the constraints to the three-neutrino scheme, using the best fit values for the 3ν mixing parameters (including the most recent hinted value of the Dirac phase δ) and exploited the MINOS bound to constrain the µ sector parameter space, (ε µe , ε µµ , ε µτ ).
Likewise, we used a KamLAND measurement limiting the number of solarν e from the source of solar 8 B ν e to place constraints on the e sector parameter space, (ε ee , ε eµ , ε eτ ). For the latter we took into account the MSW effect and decoherence of the propagating neutrinos. We again found upper limits on the absolute values of the flavour coefficients but which do not vary as a function of the Majorana phases α 2 and α 3 . We briefly mentioned that a future high statistics experiment above the τ threshold (and with a magnetic field in the far detector like OPERA) could put constraints on the τ sector parameter space, (ε τ e , ε τ µ , ε τ τ ).
Of course, as has been considered before in the literature, these coefficients can be constrained by experimental searches of microscopic LNV processes. In order to compare with the KamLAND While the CPV properties of 'ν α →ν β ' oscillations were considered in Ref. [113] in the context of the light neutrino helicity reversal mechanism, future work could study the CPV properties of complex ε coefficients. It would also be interesting to consider phenomenology of LNV NClike matter oscillations, which also generically arise from dimension-six operators in the LEFT of the SM. Constraints on extensions of the 3ν scheme to additional sterile states may also be possible given the improving precision and possible charge-sign sensitivity of future SBL reactor experiments.