Large $\nu$-$\bar{\nu}$ Oscillations from High-Dimensional Lepton Number Violating Operator

It is usually believed that the observation of the neutrino-antineutrino (${\nu}$-$\bar{\nu}$) oscillations is almost impossible since the oscillation probabilities are expected to be greatly suppressed by the square of tiny ratio of neutrino masses to energies. Such an argument is applicable to most models for neutrino mass generation based on the Weinberg operator, including the seesaw models. However, in the present paper, we shall give a counterexample to this argument, and show that large $\nu$-$\bar{\nu}$ oscillation probabilities can be obtained in a class of models in which both neutrino masses and neutrinoless double beta ($0\nu\beta\beta$) decays are induced by the high-dimensional lepton number violating operator ${\cal O}_7 = \bar{u}_R l^c_R \bar{L}_L H^*d_R + {\rm H.c.}$ with $u$ and $d$ representing the first two generations of quarks. In particular, we find that the predicted $0\nu\beta\beta$ decay rates have already placed interesting constraints on the $\nu_e \leftrightarrow \bar{\nu}_e$ oscillation. Moreover, we provide an UV-complete model to realize this scenario, in which a dark matter candidate naturally appears due to the new $U(1)_d$ symmetry.


I. INTRODUCTION
It is well established that neutrinos have tiny masses and mixings by observing neutrino oscillations between different flavors [1][2][3][4][5][6][7]. Such phenomena cannot be interpreted within the Standard Model (SM), giving us one of the strong motivations towards new physics.
A traditional smoking gun for the LNV is the neutrinoless double beta (0νββ) decays [35].
However, the LNV effects are usually predicted to be very small in most models based on the conventional dimension-5 Weinberg operator, since the corresponding amplitudes are always induced by the tiny Majorana neutrino masses and thus greatly suppressed. One prominent example is the neutrino-antineutrino (ν →ν) oscillations [36][37][38][39][40][41][42][43][44][45][46]  However, in the present paper, we shall show that the restrictions above could be waived for a class of models with the neutrino masses generated from the following specific dimension-7 LNV operator [47][48][49][50][51][52]: which is pictorially shown in Fig. 2a 1 . Here, u R , d R , l R represent the right-handed up-type quarks, down-type quarks and charged leptons, while L L and H the left-handed lepton and Higgs SU(2) L doublets withH ≡ iσ 2 H, respectively. In such a kind of models with sizable couplings only to the first two generations of quarks, the neutrino-antineutrino oscillation probabilities can be large enough so that such phenomena might be observed with the usual (anti)neutrino sources in the foreseeable future.
where we have redefined the Wilson coefficients C ud ll ′ into the basis with u R , d R and l R being the right-handed particles of the mass eigenstates while ν Ll still the flavor eigenstates. With O 6 , it is straightforward to draw the two-loop Feynman diagrams as depicted in Fig. 2b, assumed to be the dominant contributions to the neutrino masses. The explicit calculation gives the neutrino mass matrix as follows: where m u,d denotes the three generations of up-and down-type quark masses, and V ud the corresponding CKM matrix elements. Since the loop integrals I(m 2 W , m 2 u , m 2 d ) are quadratic divergent, we have estimated that I ∼ Λ 2 /(16π 2 m W ) 2 = 4Λ 2 /(16π 2 g 2 v 0 ) 2 in the second equality. If we assume all Wilson coefficients to be of O(1), the neutrino masses are dominated by the top-bottom terms leading to the UV cutoff to be Λ ∼ 1 TeV. Such a low UV cutoff makes it possible to detect sizeable LNV effects, so we shall take it in the following discussions.
Note also that the predicted neutrino mass matrix elements are proportional to the charged lepton masses. By taking into account the charged lepton mass hierarchy: m τ > m µ ≫ m e , it is generically expected that the component (m ν ) ee should be much smaller than other elements in the neutrino mass matrix. Effectively, we can take the approximation (m ν ) ee ≈ 0, which reduces to the well-studied texture-zero matrix. According to Refs. [33,34,53,54], only the normal ordering can fit the current neutrino oscillation data, which can be regarded as one of predictions of this scenario. Moreover, the nearly vanishing (m ν ) ee even restricts the lightest neutrino mass to be located within the range 0.001 eV m 0 0.01 eV, and the Majorana phase α 21 in the standard parametrization of the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) mixing matrix [55,56] to be 0.8π α 21 1.2π [33,34].
It is known that the smoking gun for the LNV is the neutrinoless double beta decay processes. For Majorana neutrinos, there always exists the traditional long-range channel shown in Fig. 3a, in which the LNV is induced by the insertion of the neutrino mass (m ν ) ee so as to flip the chirality of the internal neutrinos. Therefore, the detection of the 0νββ processes could help determine |(m ν ) ee |. Nevertheless, it has already been noted that in some types of neutrino models [27-29, 31, 33, 50, 57-60], Fig. 3a does not give the main contribution, whereas some other modes would dominate, as just the case for the present scenario shown in Fig. 3b. In Fig. 3b, the LNV occurs at the vertices due to the insertions ofÕ 6 . In order to apply the formalism in Refs. [61][62][63], we need first to transform the scalar-scalar interactions in Eq. (2) into the desired form in terms of vector-vector ones: where we have used the Fierz identity and made the charge conjugations of the currents.
We can then extract the relevant part for the 0νββ decays since the quarks and charged leptons involved in the process are u, d and e, respectively.
Note that in the above expression, we have defined the interaction coupling as according to the conventional formalism in Refs. [61][62][63].
By taking into account both decay channels in Fig. 3, we obtain the general formula for the half lifetime T 0νββ 1/2 as follows [62]: where C mm , C mλ , and C λλ can be determined by the corresponding phase space integrations and nuclear matrix elements. By using the numerical values of C listed in Table 2 of Refs. [62], we find that the first two terms in Eq. (9) are much smaller than the last one due to the suppressions of the nearly vanishing neutrino mass (m ν ) ee . As a result, it is concluded that the new long-range mode induced by O 7 gives the dominant contribution to the 0νββ decays in the present scenario.
By comparing the current experimental limits on the 0νββ half-life of different nuclei [64][65][66][67][68][69][70], we can actually give the strong constraints to the Wilson coefficient C ud ee by assuming Λ ∼ 1 TeV. The numerical limits are collected in Table I, from which it is seen that the strongest constraint is given by the target 136 Xe [65,66] with |C ud ee | < 1.9 × 10 −4 . Finally, we should mention that, contrary to the general expectation, such a stringent upper bound on C ud ee does not place any constraint to the neutrino mass element (m ν ) ee , since the neutrino mass formula in Eq. (5) depends on the different Wilson coefficients C cs ll ′ related to the second-generation quarks.
More recently, there is some interest in the estimation of the sensitivity of the µ − − e + conversion in nuclei in the literature [71][72][73], given that the the sensitivity of this mode is expected to be increased greatly in the near future due to the tremendous experimental improvement in the similar µ − −e − conversion mode. Nowadays, the most stringent limit on this channel was set by the SINDRUM II Collaboration in 1998, with the 90% C.L. upper bound on the rate as follows [74]: which was obtained by assuming that the process occurs by the coherent scattering to the ground state of calcium. With the forthcoming next-generation µ − − e − conversion experiments such as Mu2e [75] at Fermilab and COMET [76] at J-PARC in Japan, the sensitivities to the µ − − e + conversion are expected to reach R Al µ − e + ∼ 10 −16 and R Al µ − e + ∼ 10 −14 [72], respectively. Following Ref. [72], we can estimate the µ − − e + conversion rate in our scenario to be where we have neglected the two-loop contribution to the µ − -e + conversion rate, which is subdominant in our considered cutoff scale Λ ∼ 1 TeV. The first factor in Eq. (11) comes from the coefficient before the LNV operator according to our conventions. The factor G F / √ 2 is contributed by the W -boson propagator and its couplings, while 1/q 2 is the dominant contribution from the neutrino propagator and estimated to be of O(1/(100 MeV) 2 ), which is the typical distance between nucleons in a nuclei. |ψ 100 (0)| 2 is the 1s ground state probability density function of the muon in the captured atom, which has a mass dimension of 3.
Finally, all the other quantities related to the phase-space and nuclear matrix elements are characterized by the energy scale Q, and the power of Q is determined by the requirement that the final expression has the mass dimension of a decay rate. By fitting the known nuclear matrix element of titanium for the long-range light neutrino exchange, the scale Q is estimated to be 15.6 MeV [72]. By a similar argument, the muon capture rate can be approximated as Thus, by taking the ratio between Eqs. (11) and (12), we can obtain If we take Λ ∼ 1 TeV and the Wilson coefficient C ud µe, eµ ∼ 1, this ratio is of order of 10 −24 , which is too small compared with the current bound in Eq. (10) and the sensitivity of next-generation experiments.
For other LNV channels, such as rare meson decays K ± → π ∓ µ ± µ ± [77], D + → K − e + µ + [78] and rare tau decays τ − → e + π − π − [79], the sensitivities are even lower, so that we expect that they do not give rise to strong constraints to the present scenario. Therefore, we do not consider them in our following discussion.

III. NEUTRINO-ANTINEUTRINO OSCILLATIONS VIA O 7
Now we consider the neutrino-antineutrino oscillations when the effective operator O 7 gives rise to the measured Majorana neutrino masses. Besides of the conventional mechanism with the neutrino mass insertion shown in Fig. 1, O 7 can generate the new contributions to the phenomena of ν l →ν l ′ , with the Feynman diagrams presented in Fig. 4. Note that the LNV only occurs at each of two vertices in these new amplitudes, so that it is possible to induce large oscillation probabilities by avoiding the huge suppression from the mass insertions. Since the neutrino mass eigenstates propagate in space, we need first tō transform the neutrino fields in Eq. (6) into their mass states: in which the relation ν l = j U lj ν j has been used to definẽ where U represents the PMNS matrix. The amplitude for this new contribution to ν l →ν l ′ is thus given by where we have approximated the neutrino wave-functions as by assuming neutrinos propagate nearly at the speed of light so that their momenta can be estimated as p j ≈ E j − (m 2 j /2E j ). For simplicity, we have also defined Γ l ′ l j as the coefficients involving the products of C ud lj and U lj , and used K νν to denote the relevant nuclear form and kinematic factors, which are usually chosen to be real. Consequently, the probabilities for the neutrino-antineutrino oscillations are given by where we define the strong phases as ϕ jk = L(m 2 j − m 2 k )/(2E). In general, the matrix Γ l ′ l j is complex, so that it is expected that CP violation can be observed in the ν l →ν l ′ oscillations.
The CP conjugate process isν l → ν l ′ , for which the oscillation probability can be obtained as follows By taking the difference of the above two formula, we can obtain the CP asymmetry for the neutrino-anti-neutrino oscillation channels between flavors l and l ′ as follows: We now estimate the typical size of the probabilities for ν l →ν l ′ by taking the ratio of Eq. (18) with the corresponding neutrino-neutrino oscillations of the same flavor dependences: where K ν denotes the form factors for the conventional neutrino-neutrino oscillations, assuming to be K ν ≈ K νν . If Λ ∼ 1 TeV and C ud l ′ l ∼ 1 from neutrino mass calculations, the neutrino-antineutrino oscillations can be only mildly suppressed with a factor of O(10 −6 ) compared with the neutrino-neutrino counterparts. However, as shown in the previous section, the 0νββ decay experiments have already presented strong limits on the elements C ud ee < 1.9 × 10 −4 . Thus, we expect that the ν e ↔ν e channel should be much smaller than other channels. Except for ν e ↔ν e , other modes are not much constrained currently, so that their amplitudes can be large, and provide interesting signatures for this neutrino mass generation mechanism. Especially, opposed to the conventional channels in Fig. 1 where V (σ i , φ i , S, H) represents the scalar potential and is assumed to be stable so that the  which have already been discussed previously. In particular, by replacing the blob in Fig. 2b with the 1-loop box diagram in Fig. 5, the Majorana neutrino masses are generated at the three-loop level. We remark that to generate a realistic neutrino mass matrix with three different nonzero eigenvalues, at least three σ's and φ's are needed, which are labelled by the subscript i in Eq. (22).
It is quite useful to express the cutoff scale Λ and the Wilson coefficients C ud ll ′ in O 7 in terms of the parameters in the present model. Direct computations of Fig. 5 give the following relation: where I i 1 is the 1-loop integral involving the fields φ i and σ i . We have also assumed that the leptoquark mass m S is the largest mass scale in the loop, so that we can extract correct powers of m S to make the 1-loop integral I i 1 dimensionless and of O(1). If we further identify the cutoff scale in Eq. (1) as Λ ≡ m S , then the Wilson coefficients correspond to Such an identification can be further justified by the computation of the three-loop neutrino masses, given by where I i,u,d  [80,81], the µ − − e − conversions in nuclei [82][83][84][85], and threebody LFV decays l → l 1 l 2l3 [86,87], which might spoil our previous arguments based on the effective operator O 7 . But the LNV observables usually involve different coupling dependences from the LFV observables, so it is easy to evade such LFV limits by some level of tuning of model parameters. For example, the process of µ ± → e ± γ is among the most sensitive LFV probes since it usually constrains the model most stringently. In our model, there are two kinds of one-loop diagrams contributing to this process, the ones with φ ± -N loops and those with σ-χ ± loops. In order to simplify our discussion, we work with the assumption that the φ ± -N loops always dominate over the σ-χ ± ones. As a result, the main contribution to µ ∓ → e ∓ γ is where we have defined the loop integral By comparing with the upper bound B(µ + → e + γ) < 5.7 × 10 −13 from the MEG Collaboration [80] and assuming all of the φ's have the same mass of order of 1 TeV, we can obtain the following constraint on f µi,ei : With essentially the same argument, the σ-χ ± loops for µ ∓ → e ∓ γ can constrain the combination i g µi g * ei to a similar order. Furthermore, we can obtain a slightly less stringent constraint with the limit on the µ − − e − conversion in gold nuclei [82].
In face of the strong limits from the LFV processes, one may worry about their constraints on the neutrino masses and proposed neutrino-antineutrino oscillations, especially for the ν µ ↔ν e channel. However, the mismatch of the coupling dependences in LNV and LFV processes makes the advertised LNV phenomena compatible with these LFV constraints. We would like to emphasize the importance of the U(1) d symmetry in this model. Firstly, without it, the same fields may create the Weinberg operator at tree or one-loop level, spoiling the arguments above. For example, if U(1) d is replaced by Z 2 , the neutral fermions N L,R could obtain their Majorana masses so that the dominant neutrino masses come from the Ma's one-loop diagram enclosed by N and φ as in Ref. [25], which is a simple realization of the Weinberg operator. Furthermore, the presence of the leptoquark S usually involves baryon number violations, such as proton decays, via the following vertices: whereQ L ≡ iσ 2 Q L . However, the presence of the U(1) d symmetry forbids the existence of these vertices, so that the baryon number symmetry is still preserved in the present model.
A further interesting aspect of this U(1) d symmetry is that it guarantees the lightest neutral particle as a dark matter candidate, which could be the scalar from the mixing of σ and the electromagnetic neutral component of φ or the Dirac fermion N. But it is well beyond the scope of the present paper to discuss in detail dark matter physics and other aspects of this model, which we would like to present elsewhere.

V. CONCLUSIONS
Neutrino-antineutrino oscillations are one of the generic LNV phenomena. However, the conventional theories, including the seesaw neutrino models, are based on the dimension-5 Weinberg operator and predict that such effects are extremely small due to the great suppression from neutrino masses. In the present paper, we have provided a counterexample to such an expectation, in which the neutrino masses originate from the high-dimensional LNV operator O 7 in Eq. (1) with sizable couplings only to the first two generations of quarks.
In this class of models, the LNV in the ν-ν oscillations occurs at the interaction vertices with nuclear targets, rather than through the insertion of the Majorana neutrino masses, so that it escapes the suppressions in the conventional mechanism. Based on our calculations, the ν-ν oscillation probabilities are only mildly suppressed by a factor of O(10 −6 ) compared with the neutrino-neutrino oscillations, while the suppression factor for the conventional channels is of O(10 −16 ∼ 10 −22 ). Moreover, such large oscillation probabilities make it possible to observe the CP asymmetries in the ν ↔ν channels which have been shown to be generic in models with complex Wilson coefficients. Clearly, our findings reopen the hope to measure these interesting effects by using the conventional neutrino sources such as reactor and accelerator neutrinos.
There are several other interesting features in this class of models characterized by the high-dimensional operator O 7 . Due to the specific dependence on the charged lepton mass hierarchy, the neutrino mass matrix is predicted to be of normal hierarchy. Also, by fitting the neutrino mass data, the cutoff for O 7 is found to be Λ ∼ 1 TeV by assuming O(1) Wilson coefficients related to the second-generation quarks. Furthermore, the neutrinoless double beta decays are expected to be large in this kind of models, and have already imposed stringent limits to the Wilson coefficient C ud ee . In particular, the ν e ↔ν e mode is restricted to be too small to be tested experimentally. However, other modes do not suffer such strong constraints, and can be still large enough to be of phenomenological interest.
Many aspects of the present scenario are worthwhile to be investigated further. Besides of the 0νββ decays and ν →ν oscillations studied in this paper, there are other LNV effects which are also expected to be large due to new contributions from O 7 . One kind of the promising processes is the LNV rare meson decays [41,[88][89][90][91][92][93][94][95][96][97], such as D ± s → µ ± µ ± π ∓ , which might be tested by the LHCb experiments. Another kind of interesting observables involves the LNV channels at the LHC, e.g., the like-sign lepton signature, pp → l ± l ± jj [98][99][100][101][102][103]. Finally, our UV-complete model realizing O 7 as the dominant LNV operator provides a concrete connection between neutrino and dark matter physics, thus deserving detailed studies.

ACKNOWLEDGMENTS
We would like to thank Yue-Liang Wu and Bohdan Grzadkowski for helpful discussions.