Quasi Dirac neutrino oscillations

Dirac neutrino masses require two distinct neutral Weyl spinors per generation, with a special arrangement of masses and interactions with charged leptons. Once this arrangement is perturbed, lepton number is no longer conserved and neutrinos become Majorana particles. If these lepton number violating perturbations are small compared to the Dirac mass terms, neutrinos are quasi-Dirac particles. Alternatively, this scenario can be characterized by the existence of pairs of neutrinos with almost degenerate masses, and a lepton mixing matrix which has 12 angles and 12 phases. In this work we discuss the phenomenology of quasi-Dirac neutrino oscillations and derive limits on the relevant parameter space from various experiments. In one parameter perturbations of the Dirac limit, very stringent bounds can be derived on the mass splittings between the almost degenerate pairs of neutrinos. However, we also demonstrate that with suitable changes to the lepton mixing matrix, limits on such mass splittings are much weaker, or even completely absent. Finally, we consider the possibility that the mass splittings are too small to be measured and discuss bounds on the new, non-standard lepton mixing angles from current experiments for this case.


Introduction
Neutrino oscillation experiments cannot distinguish Dirac from Majorana neutrinos, hence it is still unknown whether or not lepton number is conserved. Other processes, such as neutrinoless double beta decay [1,2], need to be probed in order to answer this question. However, while the nature of neutrinos is often seen as a dichotomy, presenting two sharply distinct scenarios, the Dirac neutrino case can be seen as a limit of the more general Majorana case in which lepton number violating mass terms are zero, and this limit can be approached smoothly.
In practice, one can start with a model of 2n Majorana neutrinos and get a phenomenology arbitrarily close to the one of a model of n Dirac neutrinos. This can already be seen with only one generation of active (ν) and sterile neutrinos (N c ). In the basis (ν, N c ) T the most general mass matrix reads: If m L = m R = 0, lepton number is preserved and neutrinos are Dirac particles. This limit can alternatively be characterized by two exactly degenerate mass eigenstates composed in equal parts of ν and N c : ν 1 = 1/ √ 2 (ν + N c ) and ν 2 = i/ √ 2 (−ν + N c ). 1 Small deviations from the limit m L = m R = 0 lead to a quasi-Dirac scenario where lepton number is no longer exactly preserved.
Let us rewrite eq. (1) using: As long as ε and θ are much smaller than one, we obtain: Departures from the Dirac case therefore can manifest themselves as either new mass splittings or new mixing angles (or, in general, both). Moreover, as this simple example shows, mass splittings and mixing angles are completely independent of each other. Note that for small values of ε and θ, lepton number violation is naturally suppressed, as expected. This can be most easily seen in our one generation scenario for the double beta decay observable m ν : for θ = 0 (ε = 0) it is straightforwardly calculated to be m ν εm D ( m ν 2θm D ). We have therefore the following situation. Oscillation experiments cannot distinguish a model with n Majorana neutrinos (containing n Weyl spinors) from one with n Dirac neutrinos (containing 2n Weyl spinors) with matching masses and mixing angles. Nevertheless, once we add to a model with Dirac neutrinos small sources of lepton number violation, oscillation probabilities will change. Some illustrative examples are shown in fig. (1). We plot there the electron neutrino survival probability for low-energy (reactor) neutrinos at distances up to (and slightly larger than) the typical distances of the KamLAND experiment [3]. In all plots the black lines show the expectation for the current global best fit point [4] for the ordinary neutrino parameters in the standard three generation case, to which we have added either a non-zero mass splitting to a Dirac state (top row) or one particular new quasi-Dirac angle (bottom row).
In section (2) we will discuss the general parametrization of masses and mixing angles for scenarios with three generations of quasi-Dirac neutrinos. However, from the examples shown in fig. (1) one can read off already some basic facts about oscillations of quasi-Dirac neutrinos, which we will work out in greater detail in section (3). First, small non-zero values of ε's are equivalent to introducing new, large oscillation lengths. Thus, the best constraints on ε will come from oscillation experiments with the largest possible baselines. And secondly, even if mass splittings are negligibly small, the new, non-standard angles which appear in this setup (called θ above) may affect oscillation probabilities in a way similar to standard angles, hence creating parameter degeneracies. For example, as fig. (1) shows, from P ee alone one cannot provide limits on a single angle. (In this example variations of θ 14 can be compensated by varying θ 12 .) Even by combining more than one oscillation probability, constraints can only be derived for certain combinations of angles and phases of the mixing matrix. We will discuss this in detail in section (3.2). Constraints on mass splittings are discussed in section (3.1). Figure 1: Electron neutrino survival probability for quasi-Dirac neutrinos with a fixed energy E ν = 4 MeV as a function of distance (left), and for fixed distance L = 200 km as function of E ν (right). The standard 3-generation neutrino oscillation parameters have been fixed at their best fit point values [4], to which a small perturbation has been added. In the top row, we show the effect of mass splittings: ε 2 2 = 0 (black), ε 2 2 = 10 −5 eV 2 (orange) and ε 2 2 = 2 × 10 −5 eV 2 (red). In the bottom row, it is possible to see the effect of introducing a non-standard angle: θ 14 = 0 (black), θ 14 = π/8 (orange) and θ 14 = π/4 (red). The exact definition of ε 2 and θ 14 will be given latter in section (2).
A word on nomenclature. The terminologies quasi-Dirac and pseudo-Dirac neutrinos appear nearly interchangeably in the literature. We prefer to define quasi-Dirac (QD) neutrinos as being a mixture of active and sterile states, in contrast with pseudo-Dirac (PD) neutrinos 2 which are composed of active states only. In both cases, the structure of mass and mixing matrices must be such that the lepton sector is close to preserving one or more U (1) symmetries.
With this definition, quasi-Dirac and pseudo-Dirac neutrinos are then very different objects, both theoretically and phenomenologically. Let us briefly mention that various aspects of pseudo-Dirac neutrinos have been considered in the literature: Magnetic moments and double beta decay [7], possible mass textures [8][9][10], and oscillatory behavior [11][12][13][14]. We note in passing that models of pseudo-Dirac neutrinos require neutrino mass matrices which no longer fit the solar and atmospheric neutrino oscillation data [15][16][17]. 3 Many more papers discussed the phenomenology of quasi-Dirac neutrinos. For example, double beta decay was first discussed in this context in [6], while [19] and [20,21] consider quasi-Dirac neutrinos as a possible explanation of the atmospheric and solar neutrino problems, respectively. More ambitiously, explaining atmospheric, solar and LSND neutrino oscillations simultaneously was discussed in [22,23]. However, all these proposals are by now ruled out experimentally, since they predict too much oscillations into sterile neutrinos. Limits on quasi-Dirac neutrino parameters, on the other hand, have been derived from solar neutrino data [24] as well as from solar, atmospheric neutrino data and cosmology [25]. Furthermore, in [26][27][28][29] QD neutrinos have been discussed in the context of neutrino telescopes, such as IceCube.
Quasi-Dirac neutrino oscillations were also discussed in [30], where it was claimed that to leading order in m R,L /m D the flavor composition of the mass eigenstates does not change (only mass splittings appear), hence oscillations for n = 3 pairs of quasi-Dirac neutrinos are described by the standard mixing matrix. This assertion was taken to be true by others [26][27][28]31], yet we want to stress that this claim is not correct, as can be seen from the eqs. (4)- (6). Already for one generation, these expressions show that the mass splitting and the departure from maximal mixing are both linearly dependent on m R,L and, more importantly, they are controlled by orthogonal combinations of these two parameters. As such, it is even possible to have no mass splittings at all and at the same time have arbitrary mixing angles.
There are also a number of more theoretical papers discussing how quasi-Dirac neutrinos could arise. One possibility is the so-called "singular" seesaw where the mass matrix for the singlet neutrinos (N c ) has a determinant equal or close to zero [32]. Quasi-Dirac neutrinos from such a singular seesaw with additional type-II seesaw contributions have been discussed in [33]. Another possibility [31] involves introducing additional singlets (S), as it is done for the inverse seesaw mechanism [34]. A double seesaw is then responsible for producing very light S states which, together with the active states, form quasi-Dirac neutrinos [31]. The authors of [35] use a Dirac seesaw to explain the necessary smallness of the Dirac neutrino mass terms first, and then generate quasi-Dirac states by the addition of a very small seesaw type-II term. The "mirror world" model of [36] is another way to obtain these particles.
In models with extended gauge groups quasi-Dirac neutrinos can also appear. An example is the E 6 inspired 331 model of [37]. Here, several electroweak triplets of the gauge group SU (3) L are needed to accomodate the Standard Model leptons, and the observed active light neutrinos are automatically quasi-Dirac states [38]. A very different idea, based on supergravity has been discussed in [39]. There it was pointed out that if neutrino Dirac terms are generated from the Kähler potential (instead of the superpotential), neutrinos would be quasi-Dirac, since Majorana terms come from higher order Kähler potential terms and thus are expected to be suppressed. This idea [39] is particularly attractive, since it could, at least in principle, explain the observed smallness of the Dirac neutrino mass terms.
In addition to n active neutrinos, models of Dirac neutrinos require the introduction of n Weyl spinors transforming trivially under the electroweak gauge group. For this reason, the study of quasi-Dirac neutrinos necessarily has some overlap with the physics of sterile neutrinos. Many experiments have searched for sterile neutrinos. Most famously, the SNO neutral current measurement rules out dominant contributions of sterile neutrinos to the solar neutrino oscillations [40]. Super-Kamiokande searched for steriles in atmospheric neutrinos [41]. OPERA [42], MINOS and DayaBay [43], IceCube [44] and NOνA [45] published searches for sterile neutrinos. For a more complete list of references see the recent reviews [46,47]. Note, however, that constraints on sterile are usually derived assuming best fit point values for the standard oscillation parameters, to which two new parameters (one angle and one mass splitting) are added in the fit. This approach does not cover the general quasi-Dirac neutrino parameter space. In particular, keeping the standard neutrino parameters fixed can lead to misleading conclusions about limits for the new/extra parameters.
There are also some hints for the existence of sterile neutrinos. However, all these hints point to a new and much larger mass scale in oscillations, i.e. ∆m 2 O(1) eV 2 . Since these indications imply masses and mixings very different from those of the standard oscillations, they can not be explained by quasi-Dirac neutrinos. We thus do not discuss these hints any further and refer only to the recent review [47].
The rest of this paper is organized as follows. In section (2) we discuss the basics of quasi-Dirac oscillations, constructing general expressions for the mixing matrix for the three generation case. In section (3) we discuss constraints on the new, non-standard parameters from various neutrino experiments. Constraints on quasi-Dirac mass splittings are discussed in section (3.1), while in section (3.2) we discuss the constraints on angles, for the case in which mass splittings are negligible. We then close with a short summary and discussion.

Definitions for quasi-Dirac neutrino oscillations
Dirac neutrinos can be described either in the weak or in the mass basis. The two pictures are equivalent. We will choose the latter one. Consider then a lepton-number preserving model with three active and three sterile neutrinos (ν and N c ). 4 In the basis where the charged lepton mass matrix is diagonal, the relevant part of the Lagrangian reads In order to diagonalize the matrix m ν , both active and sterile neutrinos must be rotated, ν → V ν and Strictly speaking, the neutrino mass matrix is not yet diagonal since it is still mixing different states (active and sterile neutrinos). This can be solved by rewriting ν i and N c where the masses and the 3 × 6 mixing matrix Ω have a special form (V is a 3 × 3 square matrix): If the pattern of masses and mixing in eqs. (10) and (11) is perturbed, neutrinos are no longer Dirac particles and lepton number is violated. Note that this is equivalent to switching on the lepton number violating masses m L and m R in eq. (1). We shall now look into the possible departures from the Dirac limit as seen from the mass basis.
In the case of masses, it is possible to split the three pairs of m ψ i , m ψ i+3 , hence we may introduce three ε i such that with the understanding that, for quasi-Dirac neutrinos, the ε i are small in comparison to the atmospheric and solar mass scales. In total there are now five mass parameters relevant for oscillation experiments: the usual ∆m 2 Atm and ∆m 2 , plus three new ε i mass splittings. (As usual, the overall mass scale of neutrinos does not enter the oscillation probabilities.) Let us now turn our attention to a generic mixing matrix Ω with dimensions n × m. Such a matrix can be described by 2nm real numbers, yet orthonormality of rows (ΩΩ † = 1) imposes n 2 conditions on them, and furthermore it is possible to absorb n phases into the charged lepton fields, hence there is a total of n (2m − n − 1) real physical degrees of freedom in Ω. For a 3 × 6 matrix, this corresponds to 12 angles and 12 phases, but note that 5 of these phases cannot be observed in neutrino oscillation experiments (they correspond to column phases). The matrix Ω can be explicitly parametrized as follows [48] (called below the SV parametrization). First, consider an elementary rotation in the (i, j) entries given by the complex number θ ij ≡ θ ij exp iφ ij such that, in the (1,2) case, it has the form In the SV parametrization, the i-th row of Ω (≡ Ω i ) is then given by the expression where e (i) is a column vector with entries e (i) j = δ ij . We do not give here Ω SV in full because the expression is very lengthy.
For a particular arrangement of the 24 angles and phases in eq. (14), Ω takes the special form (11) which is associated with the Dirac limit. Note that, as usual, one can write the 3 × 3 square matrix V with three angles and one phase: Unfortunately, it is very complicated to describe the Dirac limit in the SV parametrization. Hence we make a small modification by introducing the following 6 × 6 rotation matrix: With this definition, the mixing matrix in eq. (11), with V parametrized as in eq. (15), corresponds to Ω (θ ij , φ ij ), as in (16), with θ i4 = θ i5 = θ i6 = 0 (i = 1, 2, 3), φ 12 = φ 23 = 0 and φ 13 = δ. In other words, with this definition the Dirac limit for Ω simply corresponds to keeping only the standard three generation neutrino mixing angles non-zero. We can then write the probability of neutrino oscillation from a flavor α to a flavor β for an energy E and after a length L as: 5 Note that this expression is insensitive to column rephasings Ω → Ω diag e iκ1 , e iκ2 , e iκ3 , e iκ4 , e iκ5 , e iκ6 . It is easy to show that eq. (17) reduces to the standard oscillation formula in the Dirac limit.

Current experimental limits and future prospects
As discussed in the previous section, the full parameter space for a system of 3 pairs of QD neutrinos has 30 free parameters: Two independent ∆m 2 ij plus one overall mass scale, three ε 2 i , twelve angles and twelve phases. Even discounting the five Majorana phases and the overall mass scale, which can not be probed in oscillation experiments, the remaining number of parameters is much too large to fit simultaneously. Nearly all experimental data, on the other hand, is consistent with the standard picture of only three active neutrino species participating in oscillations [4], i.e. two mass squared differences (∆m 2 Atm and ∆m 2 ), three mixing angles (θ 23 , θ 12 and θ 13 ) plus one phase (δ) are sufficient to describe the data. As mentioned in the introduction, there are also some hints for sterile neutrinos with a mass scale of the order of ∆m 2 ∼ O(eV) [46,47]. However, all these hints are at most of the order of (2 − 3) σ, we will thus not take them into account in the following. Instead, since the standard three generation picture seems to describe the data well, we will consider "small" perturbations and derive limits on particular combinations of non-standard parameters.
In order to deal effectively with the large number of parameters controlling quasi-Dirac neutrino oscillations, we will consider two simplified scenarios: 1. First, we take one non-zero ε i at a time. In these one-parameter extensions, very stringent limits on ε i are found, in agreement with earlier analysis, see for example [24,25]. We then extend this analysis to two new parameters: One mass splitting plus one new angle. This second step allows us to identify "blind spots" in the oscillation experiments, i.e. degenerate minima in particular directions in parameter space, where limits on mass splittings are much worse than in the one parameter fits. We then discuss a particular parametrization of these degenerate directions in parameter space, where the effects of ε i can be decoupled from oscillation experiments nearly completely.
2. In the second setup, we discuss the limit where mass splittings are too small to be measured in oscillation experiments, hence there are just angles and phases of Ω to deal with. In this situation, it can be shown that from the 24 parameters in Ω only 13 combinations enter the oscillation probabilities of active neutrinos. Moreover, since there is only very limited information on oscillations involving ν τ , we can in practice restrict ourselves to experiments involving ν e 's and ν µ 's. There are then only 7 combinations of the 24 angles θ ij and phases φ ij which appear in the oscillation probabilities. We discuss the construction of these 7 quantities, the current constraints and possible tests for quasi-Dirac neutrinos in this limit.
In our analysis we do not take into account the data from every existing oscillation experiment. Given the scarcity of data on τ neutrinos, we ignored it altogether, concentrating instead on the available charged current data for e and µ neutrinos and anti-neutrinos. Also, we focus on those experiments, which should provide the most important constraints for quasi-Dirac neutrinos. First, we consider KamLAND [3] since it fixes most accurately the so-called solar mass splitting ∆m 2 . From the solar neutrino experiments we fit Super-K elastic scattering data [49] and from Borexino the measured pp [50] and 7 Be fluxes [51]. This choice is motivated by the fact that Super-K [49] has the most accurate data in the high energy range, while [50,51] fix the low energy part of the solar neutrino spectrum.
We take the atmospheric neutrino data from [52]. We concentrate in our fit on the sub-GeV sample, since the lowest energetic neutrinos will be most sensitive to small values of ε 2 i , as demonstrated in fig. (2). In addition, we use data from the MINOS collaboration (muon and anti-muon neutrino survival) [53] and T2K (muon neutrino survival and muon to electron neutrino transition) [54], since these two experiments determine ∆m 2 Atm better than atmospheric data. And, finally, we take into account data from DayaBay [55], since from the three current reactor neutrino experiments DayaBay determines θ 13 with the smallest error.
In our fits, we use a simple χ 2 method to determine the allowed ranges of model parameters. We take statistical and systematic errors from the experimental publications, to which we added a further (small) systematic error for the uncertainties in our theoretical calculations. This latter systematic error was chosen such that our simulations reproduce the allowed parameter ranges for the standard oscillation parameters, determined by the respective experiment, within typically (1-1.5) σ c.l. ranges. Note that we do not attempt to do a precision global fit for standard neutrino oscillation parameters. Rather, we consider reproducing the experimental results for the standard case as a test for the reliability of our derived limits.  Table 1: 95 % upper limits on ε 2 i derived from different experimental data sets. Two numbers are given for each case; the first one is the limit obtained marginalizing over two standard oscillation parameters (see text), the second (in brackets) is the limit obtained for the best fit point value of the standard oscillation parameters. For a discussion see text.

One parameter limits
We will first discuss limits derived on ε i assuming one ε i = 0 at a time and taking all non-standard angles to be zero. Table (1) shows limits on ε 2 i and the corresponding experimental data sets used to derive the limits.
For each case listed in table (1), we have calculated the upper limits on the ε 2 i twice: (a) marginalizing over two of the standard neutrino oscillation parameters and (b) for the best fit point value of the standard parameters. Marginalization over standard oscillation parameters leads to less stringent limits. However, the importance of this marginalization procedure differs widely for different experiments. For example, in the case of KamLAND, bounds on ε 2 1 of the order of roughly 10 −5 are derived marginalizing over the allowed ranges of ∆m 2 and sin 2 θ 12 , while for the best fit values of these last two parameters, the limits are more stringent by "only" roughly a factor 2.
As the table shows, the strongest constraints on ε 2 1 and ε 2 2 come from solar neutrino data. This is easily understood from fig. (3), which shows the electron neutrino survival probability as a function of neutrino energy for different values of ε 2 2 and for the best fit values of the standard solar oscillation parameters. For low values of neutrino energies, vacuum oscillations dominate and so very small ε 2 2 can be probed up to a scale essentially determined by the Earth-Sun distance (∼ 10 −12,−11 eV 2 ). Note that a non-zero ε 2 2 reduces P ee , but this reduction could be hidden in the relatively large error bar of the low-energy measurements. 6 Nevertheless, at higher neutrino energies, a similar reduction of P ee is produced due to matter effects in the sun. Since Super-K data provides a very accurate measurement of P ee at these higher energies, one can rule out values of ε 2 2 which can not be excluded by the Borexino measurements alone. The situation is very similar for ε 2 1 : limits on ε 2 1 and ε 2 2 are then of the order of 10 −11 eV 2 if we take the best fit values of the standard solar oscillation parameters; slightly less stringent numbers are obtained when marginalizing over the standard parameter uncertainties. Figure 3: Solar neutrino survival probability as a function of neutrino energy, for different choices of ε 2 2 . Solar angle and mass splitting have been fixed at their best fit values in this plot [4].
Solar neutrino experiments have essentially no sensitivity to ε 2 3 . This is simply due to the smallness of θ 13 (sin 2 θ 13 0.0215 [4]). Thus, we have to rely on experiments testing the atmospheric scale to derive limits on ε 2 3 . Table (1) quotes numbers for two cases. In the first scenario, we have combined data from DayaBay [55], T2K [54] and MINOS [53]: DayaBay fixes most accurately θ 13 , while both MINOS and T2K measure ∆m 2 Atm with rather small errors. Here, the limit on ε 2 2 is (not surprisingly) less stringent than the one derived from KamLAND (or solar). Depending on whether or not ∆m 2 Atm and sin 2 θ 23 are used at their best fit value or marginalized over, we get very different limits on ε 2 3 . This is due to the fact that when scanning over the standard oscillation parameters, the χ 2 function has two almost degenerate minima: one for small values of ε 2 3 and another for ε 2 3 of the order 10 −3 eV 2 . However, as the table also shows (second case), this non-standard solution is excluded, once we add Super-K atmospheric neutrino sub-GeV data to the fit. With the combination of these four experiments limits on ε 2 2 and ε 2 3 are again of order 10 −5 eV 2 . In the last line of table (1) we give our forecast of the sensitivity of the planned experiment JUNO [56]. JUNO will measure ∆m 2 and ∆m 2 Atm very precisely and thus it will also be able to derive limits on any ε 2 i . However, our results indicate that, despite being a very precise experiment, JUNO will not lead to a major improvement over existing limits on ε 2 i . Here, it is important to stress that limits using the best fit point and limits marginalizing over standard parameter uncertainties are very different. This can be traced back again to a near-degeneracy in the χ 2 -function: For 2 i of the order of ∆m 2 ij one has two only slightly different oscillation lengths contributing in the fit, which can give a better description than a single oscillation length.
In summary, strong limits on mass splittings can be derived from atmospheric and solar neutrino data (ε 2 3 ∼ O(10 −5 ) eV 2 and ε 2 1,2 ∼ O(10 −11 ) eV 2 , respectively) in the case where no other extra parameter is added to the standard neutrino oscillation picture.

Two parameter case
While the discussion in the previous subsection seems to show that constraints on QD mass splittings are very stringent, we will now see that this conclusion is valid only under the assumption that no other non-standard parameter is different from zero.
As a simple example, consider the electron neutrino survival probability at distances short enough that the effects of ∆m 2 can be neglected. 7 We shall consider the particular example of a non-zero ε 2 3 and a non-zero θ 16 angle, defined in section 2. One finds that where c ij and s ij are short-hands for cos θ ij and sin θ ij and It is straightforward to see that the above expression for the neutrino survival probability remains (nearly) unchanged if we swap θ 13 by θ 16 . This is true up to very small terms proportional to ∆(P ee ) ∝ (∆m − ee − ∆m + ee )(c 13 − c 16 )s 13 s 16 . More specifically, in the limit where ∆m ± ee have the same value, P ee is only a function of the combination sin 2 θ 13 + sin 2 θ 16 . Thus, there will be a near-degeneracy of the relevant χ 2 function involving these two angles, and so values (or limits) derived for one of these parameters, without varying the other, will be misleading.
There is, however, another more interesting degeneracy associated to eq. (18). In calculating this expression we have used a certain parametrization for the mass splitting, which we may call the symmetric parametrization: Choosing sin θ 13 = tan θ 16 , the second term inside the bracket in eq. (18) vanishes (this choice corresponds to Ω e6 = 0). So, by adjusting ∆m 2 31 and ε 2 3 we can keep ∆m 2 31 − ε 2 3 /2 constant and equal to the best fit point value of ∆m 2 Atm , in which case there will be no upper limit on ε 2 3 itself coming from the electron neutrino survival experiments. Note that we could have defined m i , m i+3 → m i , m 2 i + ε 2 i . 8 We call this the asymmetric parametrization. Rewriting eq. (18) with this parametrization, the first term inside the bracket would not depend on ε 2 3 at all, so it becomes obvious that for the choice of sin θ 13 = tan θ 16 all dependence of P ee on ε 2 3 disappears. Fig. (4) shows these parameter degeneracies in the space (ε 2 3 , sin 2 θ 13 , sin 2 θ 16 ), using only the DayaBay data (on the left column). The underlying scan was done in the asymmetric parametrization, which is numerically simpler to implement. The plot in the upper and middle panel show clearly that there is no upper limit on ε 2 3 in this scan. The lower plot shows the degeneracy in parameter space under the exchange of θ 13 ↔ θ 16 .
We can break this particular degeneracy in parameter space, by adding more experiments. T2K measures two probabilities: (a) The muon neutrino survival probability, P µµ , and (b) the electron neutrino appearance probability, P µe , both at values of L/E which give access to the atmospheric neutrino mass scale, ∆m 2 Atm . If the only non-standard angle different from zero is θ 16 , then P µµ will not depend on θ 16 at all, while P µe will have θ 16 -dependence which is different from the one of P ee . Thus, adding T2K data to the scan is enough to break the degeneracy in θ 13 ↔ θ 16 hence an upper limit on ε 2 3 reappears. The middle column of fig. (4) illustrates this point; it shows the results of a combined scan over ε 2 3 , sin 2 θ 13 and sin 2 θ 16 for DayaBay plus T2K data. By comparison of the right with the middle column of fig. (4), one can clearly see that the addition of Super-K data generates a strong upper limit on ε 2 3 , for this particular choice of parameter subspace.
Given these results, one might wonder if there are particular directions in parameter space for which oscillation experiments become completely blind to QD neutrino mass splittings. Recall that the blind (or: degenerate) direction discussed above for P ee corresponds to the particular choice of Ω e6 = 0. In a similar way, for example, P µµ would loose any sensitivity to ε 2 3 if Ω µ6 = 0. Thus, with some special choice of θ 16 and θ 26 such that both Ω e6 and Ω µ6 are zero at the same time, one can indeed make DayaBay and T2K blind to variations of ε 2 3 .

DayaBay
DayaBay+T2K DayaBay+T2K+SuperK Figure 4: Allowed parameter ranges for ε 2 3 , sin 2 θ 13 and sin 2 θ 16 for different experiments. The parameter planes always marginalize over the parameter not shown and all calculations used the best fit point value for ∆m 2 Atm . In the plots on the left only DayaBay data is taken into account; the middle panel combines DayaBay with T2K and the panel to the right shows the combination of DayaBay, T2K with Super-K atmospheric neutrino data. The different coloured regions present the 1, 2 and 3 σ c.l. allowed regions (cyan, blue and red). For discussion see text.
While it is possible, in principle, to calculate the combination of angles θ ij and phases φ ij (defined in section (2)) associated to these blind directions, in the following we will consider a simpler alternative. Consider a unitary rotation of the columns i and i + 3 of the mixing matrix Ω. Since we are not interested in column phases, such a rotation is governed by just two parameters (ϕ i and β i ): Now, recall that in the Dirac limit (see eq. (11)) the columns i and i + 3 of the mixing matrix Ω are proportional to each other, This means that applying the (ϕ i , β i ) transformation to the Dirac neutrino mixing matrix and, without loss of generality setting β i = 0, we obtain: .
From the last of these equations it can be seen that the i'th column (the (i + 3)'th column) of Ω vanishes, if one chooses ϕ i = 3π/4 (ϕ i = π/4).
DayaBay+SK+MINOS+T2K and Juno Figure 5: Allowed parameter space in the plane (ε 2 3 , ϕ 3 ) using DayaBay, T2K, MINOS and Super-K atmosheric neutrino data. The coloured plane shows the 2 and 3 σ c.l. allowed regions (blue and red). The dashed lines show the expected limits for JUNO. The asymmetric parametrization of the mass splitting (m 3 , m 2 3 + ε 2 3 ) was used, so for ϕ 3 = π/4 there is no sensitivity to ε 2 3 . Fig. (5) shows a scan over the allowed parameter space in the plane (ε 2 3 , ϕ 3 ) using DayaBay, T2K, MINOS and Super-K atmospheric neutrino data. In agreement with the above discussion, there is a blind spot where no limit on ε 2 3 exist. This blind direction corresponds to the choice of ϕ 3 = π/4. 9 Fig. (5) also shows that the addition of JUNO data can lead only to a marginally improved limit. 9 Shifting ε 2 3 one could alternatively define ( m 2 3 + ε 2 3 , m 3 ). In that case, the blind spot occurs at ϕ 3 = 3π/4 instead.
We now turn to a discussion of ε 2 1 and ε 2 2 . For these two parameters, again solar neutrino physics provides the most important constraints. As above for ε 2 3 , we can define a rotation angle ϕ 1 (ϕ 2 ) between the columns 1 and 4 (2 and 5) of the mixing angle which will mitigate the effects of a non-zero ε 2 1 (ε 2 2 ). Fig.  (6) shows the P ee probability for solar neutrinos as a function of neutrino energy for two different values of ε 2 2 and various values of ϕ 2 . 10 The results for ε 2 1 and ϕ 1 are completely analogous. As the figure shows, for ϕ 2 = π/4 again the effects of ε 2 2 completely decouple from the oscillation probability.   (7) shows the allowed parameter space in the two planes (ϕ 1 , ε 2 1 ) and (ϕ 2 , ε 2 2 ) using solar data and combining solar data with KamLAND. These plots have been calculated using the best fit point values for ∆m 2 and sin 2 θ 12 from the global fit [4]. The plots show in all cases that there exists a slight preference, between (2-2.5) σ in all cases, for non-zero values of ε 2 i . Note that the preferred solution of the solar data in the region of ϕ 1 ∼ 3/4π and ε 2 1 ∼ (10 −4.5 − 10 −4 ) is ruled out by KamLAND. However, even combining solar and KamLAND data some preference for non-zero ε 2 i of the order of very roughly 10 −10.5 eV 2 remains. We have traced back this preference for non-zero mass splittings in solar data to the well-known difference in the best fit points from ∆m 2 in solar and KamLAND data. As can be seen also in the latest global fits [4], solar data prefers a ∆m 2 around (4 − 5) × 10 −5 eV 2 , while KamLAND prefers ∆m 2 7.6 × 10 −5 eV 2 . This tension between the two data sets is roughly of the order of 2 σ, with the error bar dominated by the larger error on ∆m 2 in the solar data set. We have therefore recalculated the constraints on (ϕ 1 , ε 2 1 ) from solar data for a value of ∆m 2 = 4 × 10 −5 eV 2 . Fig. (8) shows the results of such a scan. As can be seen, in this calculation there is no longer any preference for a non-zero value of ε 2 1 . Note that such a low value of ∆m 2 is ruled out by many σ from the KamLAND data. Thus, a small non-zero mass splitting could provide, in principle, a solution for the observed tension between solar and KamLAND data.
In summary, by introducing one mass splitting at a time, we extracted bounds for the pairs of parameters (ε i , ϕ i ), i = 1, 2, 3. In the limit where ϕ i is (2 ± 1) /4π, one column of the mixing matrix vanishes and therefore the mass splitting ε i becomes unobservable. For this reason, one expects that for reasonably large values of ε i there must be tight limits on |ϕ i − (2 ± 1) /4π|, meaning that ϕ i has to be quite far from the Dirac limit (ϕ i = 0). On the other hand, for a small enough value of the mass splitting ε i , the associated oscillation length eventually become larger than the baseline of the relevant experiments, and in that case ϕ i becomes unconstrained. Figure 7: Allowed parameter range in the space ϕ 1 , ε 2 1 (top) and ϕ 2 , ε 2 2 (bottom). To the left: Solar data, to the right solar data + KamLAND. This plot uses the best fit point values for ∆m 2 and sin 2 θ 12 from the global fit. This combination of data shows a slight preference for a non-zero value of the mass splitting, for a discussion see text. Figure 8: Allowed parameter range in the space ϕ 1 , ε 2 1 (left) and ϕ 2 , ε 2 2 (right) from solar data, using ∆m 2 = 4 × 10 −5 eV 2 .

Quasi-Dirac neutrinos in the limit ε i → 0
If masses are degenerate as in eq. (10), then the oscillation probability formula will not change under unitary rotations of the columns i and i + 3 of the mixing matrix, see eq. (20). In other words, for unitary matrices U (i) (i = 1, 2, 3) leaves P (ν α → ν β ) unchanged. Hence, there is a U (2) 3 redundancy in our description of Ω, and in turn this means that out of the 24 parameters describing the mixing matrix, oscillation experiments are only sensitive to 13. 11 This number is further reduced to 7 if we ignore tau neutrinos. In this latter case the oscillation probabilities can be written as: with the oscillating factors A ij ≡ −4 sin 2 m 2 i − m 2 j L/ (4E) and B ij ≡ 2 sin m 2 i − m 2 j L/ (2E) and the 7 parameters X i defined as follows: As a side remark, we would like to point out here that a similar approach could, in principle, be used in the presence of one ε i : in this case, instead of 7, there would be 9 combinations of angles and phases to take into account. 12 Note that the X i defined above can take any value in our framework, provided that the following constraints are obeyed: 1. the first six X i are non-negative real numbers; 2. neither X 1 + X 2 nor X 3 + X 4 can be larger than 1; 3. X 5 ≤ X 1 X 3 and X 6 ≤ X 2 X 4 ; 4. the norm of X 7 is fixed by X 5 and X 6 (|X 7 | 2 = X 5 X 6 ), so even though X 7 is a complex parameter, only arg (X 7 ) is an independent degree of freedom; These conditions are a consequence of the definitions of the X i and the fact that the rows of the mixing matrix Ω are orthonormal (ΩΩ † = 1). Taking them into account, we are able to pick out all valid points in the X i parameter space, without ever referencing back to specific entries of the mixing matrix. For reference, the values of these X i parameters in the Dirac limit as a function of the standard θ 12 , θ 13 , θ 23 and δ parameters (see eq. (15)) are the following: X 1 = sin 2 θ 13 , X 2 = sin 2 θ 12 cos 2 θ 13 , X 3 = cos 2 θ 13 sin 2 θ 23 , X 4 = sin 2 θ 12 sin 2 θ 13 sin 2 θ 23 + cos 2 θ 12 cos 2 θ 23 − 1 2 sin 2θ 12 sin 2θ 23 sin θ 13 cos δ , X 7 = sin θ 12 sin θ 13 cos 2 θ 13 sin θ 23 − sin θ 12 sin θ 13 sin θ 23 + e −iδ cos θ 12 cos θ 23 .
Overall, the bounds on these 7 parameters are broadly consistent with the standard three neutrino oscillation picture. In other words, by substituting in expressions (34)-(37) the numbers obtained for θ 12 , θ 13 , θ 23 and δ from global fits [4], we get values for the X i roughly in agreement with figs. (9) and (10). In order to see clearly that current data is consistent with the Dirac limit, note that in this latter case there are only 4 independent parameters. Thus, it follows that the standard three neutrino oscillation picture must correspond to three relations among the 7 X i . These are and from fig. (11) one can see that oscillation data is compatible with each of these equalities within ∼ 1σ. The three together are disfavored only at min χ 2 Dirac − min χ 2 = 1.9 so, assuming no mass splittings ε i , there is currently no significant indication for quasi-Dirac neutrinos.
arg X 9 ≡ arg Ω e1 Ω * µ1 (Ω * e2 Ω µ2 + Ω * e5 Ω µ5 ) . Figure 9: One, two and three σ regions in the (X 1 , X 2 ), (X 3 , X 4 ) and (X 5 , X 6 ) which are allowed by electron neutrino survival data from at KamLAND, DayaBay, SuperK and Borexino, muon neutrino survival data at MINOS and T2K and muon to electron transition data from T2K. This plots were obtained from a 7-dimensional scan of the X i defined in eqs. (26)- (29) and marginalizing over 5 variables.  (29). The χ 2 function for X 7 is essentially flat. These parameters are a function of the entries of the mixing matrix only; no mass splittings ε i were considered. Figure 11: ∆χ 2 functions for the three combinations of parameters which, when equal to 0 simultaneously, signal the Dirac limit (see eq. (38)).
We would like to point out that even in the absence of new mass scales, quasi-Dirac neutrinos can, in principle, be distinguished from some other scenarios through oscillation experiments. In particular, consider 3 active neutrinos and a non-unitary 3 × 3 mixing matrix V . The oscillation probabilities are then given by the expressions P (ν e → ν e ) = F J 12 ee , J 13 ee , J 23 ee + J 12 where The exact form of the functions F and F which control the 0-distance neutrino behavior is not important for the present discussion. The more important point is that with a non-unitary V one can have an oscillatory behavior which is impossible to reproduce with quasi-Dirac neutrinos, and vice-versa.
For example, the B ij coefficients in P (ν e → ν µ ) do not need to be related for a non-unitary V , while for quasi-Dirac neutrinos they must be the same (up to a minus sign -see eq. (25)). On the other hand, note that from the J ij ee and J ij µµ one can extract the modulus of the absolute value of all J ij eµ , hence by measuring P (ν e → ν e ) and P (ν µ → ν µ ), as well as the coefficients B ij in P (ν e → ν µ ), the coefficients of the oscillatory factors A ij in P (ν e → ν µ ) are fixed for a non-unitary V (up to ± signs). However, for quasi-Dirac neutrinos no such constraint exists. So, with this short theoretical argument, one can conclude that in principle these two non-standard neutrino scenarios can be distinguished through oscillation experiments.

Summary
In general, neutrinos can have lepton number violating (Majorana) and lepton number conserving (Dirac) mass terms. If the lepton number violating mass terms are smaller than the lepton number preserving ones, neutrinos are quasi-Dirac particles. Phenomenologically, this corresponds to the existence of three pairs of neutrinos with slightly different masses, hence oscillation experiments are sensitive not only to the usual solar and atmospheric mass scales, but also to three small mass splittings ε i . Furthermore, for quasi-Dirac neutrinos there are more than 3+1 angles and phases to be considered. In this work, we have analyzed the constraints on these quasi-Dirac neutrino parameters imposed by current neutrino oscillation data and also briefly discussed the potential of the future JUNO experiment to improve upon existing constraints.
In section (2) we have discussed a fully general parametrization of the lepton sector for three generations of quasi-Dirac neutrinos. In addition to the charged lepton masses, there is a total of 6 masses, 12 angles and 12 phases. Oscillation experiments are not sensitive to the overall neutrino mass scale nor to 5 of the phases (which are of the Majorana type). Hence we are left with a 24-dimensional model space, compared to the six-dimensional space for an ordinary three generation case (∆m 2 , ∆m 2 Atm , θ 12 , θ 13 , θ 23 and δ). It is numerically too costly to handle such a large number of parameters at the same time, hence we analyzed several different special cases. First, we took a single mass splitting ε 2 i = 0. If we split two neutrinos with mass m i into a quasi-degenerate pair of particles with masses m 2 i − ε 2 i /2 and m 2 i + ε 2 i /2 then a new oscillation length L ∝ 1/ε 2 i appears which is associated to the conversion of active to sterile neutrinos. Very stringent limits on ε 2 i in such one parameter extensions can be derived, of the order of 10 −11 eV 2 for ε 2 1,2 (from solar neutrino data) and 10 −5 eV 2 for ε 2 3 (dominated by Super-K atmospheric neutrino data).
Next, we considered the case when one mass splitting and one of the non-standard angles are allowed to take non-zero values at the same time. As we have shown, in this situation degeneracies of the χ 2 function can occur, implying that from a single experiment in many cases it will no longer be possible to derive meaningful limits on individual parameters. These degeneracies can be resolved by considering data from more than one experiment, accessing different P (ν α → ν β ).
We then considered the possibility of nullifying the effects of the ε i completely by changing some particular combinations of the angles θ ij of our parametrization. Instead of pursuing the exact form of these rather complex parameter combinations, we discussed a simpler definition, describing 3 angles ϕ i associated to rotations between the columns i and i + 3 of the quasi-Dirac mixing matrix, such that in the limit where these angles are equal to π/4 (3π/4) the i + 3(i) column of the mixing matrix vanishes and hence the associated neutrino mass disappears from the oscillation probability formula. We stress that for these particular parameter combinations no limits on the ε i can be derived from oscillation experiments. The regions in the planes (ε i , ϕ i ) which are allowed by various experiments are shown in figs. (5) and (7). In this context, it is interesting to note that the tension between the value of the solar mass scale preferred by global fits (∼ 7.6 × 10 −3 eV 2 ) and the lower one preferred by solar data (∼ 4 × 10 −3 eV 2 ) might be resolved by a non-zero value for either ε 1 or ε 2 .
Lastly, we considered the possibility that the mass splittings ε i are too small to be measured in oscillation experiments. Even in this scenario, one can have departures from the lepton-number-conserving Dirac scenario due to the new angles θ ij (and phases φ ij ). As mentioned above, there is a large number of such parameters. However, it can be shown that with 3 pairs of neutrinos with the same mass, oscillations will only depend on a total of 13 combinations of angles and phases. Additionally, if we focus just on electron and muon neutrinos, this number is further reduced to 7, corresponding to 6 angles and 1 phase. In the text we called these parameter combinations X 1···7 and stressed that they can not be identified with θ 12 , θ 13 , θ 23 nor δ, as these quantities by themselves are not physical. Instead, the 7 X i correspond to combinations of these and additional θ ij angles and φ ij phases.
In section (3.2) we made a 7-dimensional scan of these X i parameters in the absence of mass splittings. Their exact definitions, as well as the limits imposed on them by current data can be found there. Crucially, for Dirac neutrinos there are only 4 parameters. Hence the Dirac limit corresponds to 3 relations among the 7 X i . By testing these relations, we find that min χ 2 Dirac − min χ 2 = 1.9, i.e. current data is compatible with the Dirac scenario. Progress on tests for quasi-Diracness can be made in the future with a more precise measurement of P (ν e → ν µ ) and P (ν µ → ν µ ). Thus, more statistics taken in T2K, MINOS+ or NOνA and, in particular, the future precise measurements possible at DUNE should provide more sensitive probes for this particular setup of quasi-Dirac neutrinos without new mass scales.