Long-range contributions to double beta decay revisited

We discuss the systematic decomposition of all dimension-7 (d=7) lepton number violating operators. These d=7 operators produce momentum enhanced contributions to the long-range part of the neutrinoless double beta decay amplitude and thus are severely constrained by existing half-live limits. In our list of possible models one can find contributions to the long-range amplitude discussed previously in the literature, such as the left-right symmetric model or scalar leptoquarks, as well as some new models not considered before. The d=7 operators generate Majorana neutrino mass terms either at tree-level, 1-loop or 2-loop level. We systematically compare constraints derived from the mass mechanism to those derived from the long-range double beta decay amplitude and classify our list of models accordingly. We also study one particular example decomposition, which produces neutrino masses at 2-loop level, can fit oscillation data and yields a large contribution to the long-range double beta decay amplitude, in some detail.


I. INTRODUCTION
Majorana neutrino masses, lepton number violation and neutrinoless double beta decay (0νββ) are intimately related. It is therefore not surprising that many models contributing to 0νββ have been discussed in the literature, see for example the recent reviews [1,2]. However, the famous black-box theorem [3] guarantees only that -if 0νββ decay is observed -Majorana neutrino masses must appear at the 4-loop level, which is much too small [4] to explain current oscillation data [5]. Thus, a priori one does not know whether some "exotic" contribution or the mass mechanism dominates the 0νββ decay rate. Distinguishing the different contributions would not only be an important step towards determining the origin of neutrino masses, but would also have profound implications for leptogenesis [6,7].
In terms of only standard model (SM) fields, ∆L = 2 terms can be written as nonrenormalizable operators (NROs) of odd mass dimensions. At mass dimension d = 5, there is only one such operator, the famous Weinberg operator [8], O W = 1 Λ (LLHH). At tree-level the Weinberg operator can be understood as the low-energy limit of one of the three possible seesaw realizations [9][10][11][12][13]. All other ∆L = 2 operators up to d = 11 -excluding, however, possible operators containing derivatives -have been listed in [14]. When complemented with SM Yukawa interactions (and in some cases SM charged current interactions), these higher dimensional operators always also generate Majorana neutrino masses (at different loop-levels), leading again to the Weinberg operator 1 at low energies.
All ∆L = 2 operators also contribute to 0νββ decay. From the nuclear point of view, the amplitude for 0νββ decay contains two parts: the long-range part and the short-range part. The so-called long-range part [23] describes all contributions involving the exchange of a light, virtual neutrino between two nucleons. This category contains the mass mechanism, i.e. the Weinberg operator sandwiched between two SM charged current interactions, and also contributions due to d = 7 lepton number violating operators. 2 The short-range part of the 0νββ decay amplitude [24], on the other hand, contains all contributions from the exchange of heavy particles and can be described by a certain subset of the d = 9 ∆L = 2 operators in the list of [14]. In total there are six d = 9 operators contributing to the shortrange part of the amplitude at tree-level and the complete decomposition for the (scalar induced) operators has been given in [25]. The relation of all these decompositions with neutrino mass models has been studied recently in [26]. 3 The general conclusion of [26] is that for 2-loop and 3-loop neutrino mass models, the short-range part of the amplitude could be as important as the mass mechanism, while for tree-level and 1-loop models one expects that the mass mechanism gives the dominant contribution to 0νββ decay. 4 In this paper we study d = 7 ∆L = 2 operators, their relation to neutrino masses and the long-range part of the 0νββ decay amplitude. We decompose all d = 7 ∆L = 2 operators and determine the level of perturbation theory, at which the different decompositions (or "proto-models") will generate neutrino masses. Tree-level, 1-loop and 2-loop neutrino mass models are found in the list of the decompositions. We then compare the contribution from the mass mechanism to the 0νββ decay amplitude with the long-range d = 7 contribution. Depending on which particular nuclear operator is generated, limits on the new physics scale Λ > ∼ g eff (17 − 180) TeV can be derived from the d = 7 contribution. Here, g eff is the mean of the couplings entering the (decomposed) d = 7 operator. This should be compared to limits of the order of roughly Λ > ∼ √ Y eff 10 11 TeV and Λ > ∼ Y 2 eff 50 TeV, derived from the upper limit on m ν for tree-level and 2-loop (d = 7) neutrino masses. (Here, Y eff is again some mean of couplings entering the neutrino mass diagram. We use a different symbol, to remind that Y eff is not necessarily the same combination of couplings as g eff .) Thus, only for a certain, well-defined subset of models can the contribution from the long-range amplitude be expected to be similar to or dominate over the mass mechanism. Note that, conversely a sub-dominant contribution to the long-range amplitude always exists also in all models with mass mechanism dominance.
We then give the complete classification of all models contributing to the d = 7 operators in tabular form in the appendix of this paper. In this list all models giving long-range contributions to 0νββ decay can be found, such as, for example, supersymmetric models with R-parity violation [32,33] or scalar leptoquarks [34]. There are also models with non-SM vectors, which could fit into models with extended gauge sectors, such as the left-right symmetric model [35][36][37]. And, finally, there are new models in this list, not considered in the literature previously.
We mention that our paper has some overlap with the recent work [38]. The authors of this paper also studied d = 7 ∆L = 2 operators. 5 They discuss 1-loop neutrino masses induced by these operators, lepton flavour violating decays and, in particular, LHC phenomenology for one example operator in detail. The main differences between our work and theirs is that we (a) focus here on the relation of these operators with the long-range amplitude of 0νββ decay, which was not studied in [38] and (b) also discuss tree-level and 2-loop neutrino mass models. In particular, we find that 2-loop neutrino mass models are particularly interesting, because the d = 7 long-range contribution dominates 0νββ only in the class of models.
The rest of this paper is organized as follows. In the next section we lay the basis for the discussion, establishing the notation and recalling the main definitions for ∆L = 2 operators and 0νββ decay amplitude. In the following section we then discuss an example of each: tree-level, 1-loop and 2-loop neutrino mass models. In each case we estimate the contribution to the mass mechanism and the constraints from the long-range amplitude. We study a 2-loop d = 7 model in some more detail, comparing also to oscillation data and discuss the constraint from lepton flavour violating processes. In section IV we then discuss a special case, where a d = 9 operator can give an equally important contribution to the 0νββ decay amplitude as a d = 7 operator. The example we discuss is related to the left-right symmetric extension of the standard model and, thus, of particular interest. We then close the paper with a short summary. The complete list of decompositions for d = 7 operators is given as an appendix.

II. GENERAL SETUP
The 0νββ decay amplitude can be separated into two pieces: (a) the long-range part [23], including the well-known mass mechanism, and (b) the short-range part [24] of the decay rate describing heavy particle exchange. Here, we will concentrate exclusively on the long-range part of the amplitude.
The long-range part of the amplitude exchanges a light, virtual neutrino between two point-like vertices. The numerator of the neutrino propagator involves two pieces, (m ν i +p / ). If the interaction vertices contain standard model charged current interactions, the m ν i -term is projected out. This yields the "mass mechanism" of 0νββ decay. However, if one of the two vertices involved in the diagram produces a neutrino in the wrong helicity state, i.e. (ν L ) c , the p / -term is picked from the propagator. Since the momentum of the virtual neutrino is typically of the order of the Fermi momentum of the nucleons, p F ≃ 100 MeV, the 0νββ amplitude from the operators proportional to p / is enhanced by p F /m ν O(10 8 ) 5 Decompositions of d = 7 operators were also discussed in [39,40].
136 Xe 2.0 · 10 −9 3.9 · 10 −7 4.7 · 10 −9 4.7 · 10 −9 3.3 · 10 −10 5.6 · 10 −10 with respect to the amplitude of the standard mass mechanism. Consequently, any operator proportional to p / will be tightly constrained from non-observation of double beta decay. Following [23] we write the effective Lagrangian for 4-fermion interactions as The leptonic (hadronic) currents j β (J α ) are defined as: where γ µν is defined as γ µν = i 2 [γ µ , γ ν ]. The first term of Eq. (1) is the SM charged current interaction, the other terms contain all new physics contributions. We normalize the coefficients ǫ β α relative to the SM charged current strength G F / √ 2. Recall, P L/R = 1 2 (1 ∓ γ 5 ) and we will use the subscripts L and R for left-handed and right-handed fermions, respectively. Note also that all leptonic currents with (1 − γ 5 ) will pick m ν i from the propagator, leading to an amplitude proportional to ǫ β L α × m ν (β L ∈ {S − P, V − A, T L }), which is always smaller than the standard mass mechanism contribution and thus is not very interesting. Thus, only six particular ǫ β α can be constrained from 0νββ decay. For convenience, we repeat the currently best limits, all derived in [1], in Table I. Eq.(1) describes long-range 0νββ decay from the low-energy point of view. From the particle physics point of view, these ∆L = 2 currents can be described as being generated from d = 7 operators. Disregarding the d = 7 "Weinberg-like" operator O W × (H † H), there are four of these operators in the list of Babu & Leung [14]: Here, O 2 is included for completeness, although it is trivial that the mass mechanism will be the dominant contribution to 0νββ decay for this operator, since it does not involve any quark fields. We will therefore not discuss the detailed decomposition of O 2 , which can be found in [38]. The operators O 3b,4a,8 will contribute to the long-range amplitudes j β J α , and the coefficient of the amplitudes is described as where Λ 7 is the energy scale from which the d = 7 operators originate, and ǫ d=7 is one of (or a combination of two of) the ǫ β α of Table 1. The factor 1/4 is included to account for the fact that Eq. (2) is written in terms of (1 ± γ 5 ) while chiral fields are defined using P L/R . This leads to the numerical constraints on the scale Λ 7 mentioned in the introduction, taking the least/most stringent numbers from Table I. All ∆L = 2 operators generate Majorana neutrino masses. However, operators O 3a and O 4b will generate neutrino mass matrices without diagonal entries, since L i L j ǫ ij = 0 within a generation. Neutrino mass matrices with such a flavour structure result in very restricted neutrino spectra, and it was shown in [42] that such models necessarily predict sin 2 (2θ 12 ) = 1 − (1/16)(∆m 2 21 /∆m 2 31 ) 2 . This prediction is ruled out by current neutrino data at more than 8 σ c.l. [5]. Models that generate at low energies only O 3a or O 4b can therefore not be considered realistic explanation of neutrino data. 6 Flavour off-diagonality of O 3a and O 4b does also suppress strongly their contribution to long-range double beta decay, in case the resulting leptonic current is of type j S+P (see appendix 7 ). This is because the final state leptons are both electrons, while the virtual neutrino emitted from the L in O 3a,4b is necessarily either ν µ or ν τ . In the definition of the "effective" ǫ β α , then neutrino mixing matrices appear with the combination j U ej U * µj (or U ej U * τ j ), which is identically zero unless the mixing matrices are non-unitary when summed over the light neutrinos.
Departures from unitarity can occur in models with extra (sterile/right-handed) neutrinos heavier than about ∼ 1 GeV. While the propagation of the heavy neutrinos also contributes to 0νββ, the nuclear matrix element appearing in the amplitude of the heavy neutrino exchange is strongly suppressed, when their masses are larger than 1 GeV [44,45]. Consequently, the heavy neutrino contribution is suppressed with respect to the light neutrino one and the sum over j U ej U * µj is incomplete, appearing effectively as a sum over mixing matrix elements which is non-unitary. Current limits on this non-unitary piece of the mixing are of the order of very roughly percent [46][47][48][49], thus weakening limits on the coefficients for O 3a and O 4b (for j S+P ), compared to other operators, by at least two orders of magnitude. 6 However, models that produce these operators usually allow to add additional interactions that will generate O 5 (O 6 ) in addition to O 3a (O 4b ), as for example in the model discussed in [43]. These constructions then allow to correctly explain neutrino oscillation data, since O 5 /O 6 produce non-zero elements in the diagonal entries of the neutrino mass matrix. 7 Decomposition #8 of O 3a also generates j TR which can contribute to 0νββ without the need for a nonunitarity of the mixing matrix.
To the list in Eq. (3) one can add two more ∆L = 2 operators involving derivatives: We mention these operators for completeness. As shown in [50], tree-level decompositions of O Dµ 1 always involve one of the seesaw mediators, and thus one expects this operator to be always present in tree-level models of neutrino mass. As we will see, if neutrino masses are generated from tree-level, the mass mechanism contribution in general dominates 0νββ, and consequently the new physics effect from O Dµ 1 cannot make a measurable impact. The second type of the derivative operators, O Dµ 2 , has also been discussed in detail in [50] with an example of tree-level realization, we thus give only a brief summary for this operator in the appendix.

III. CLASSIFICATION
In this section we will discuss a classification scheme for the decompositions of the ∆L = 2 operators of Eq. (3), based on the number of loops, at which they generate neutrino masses. We will discuss one typical example each for tree-level, 1-loop and 2-loop models. The complete list of decompositions for the different cases can be found in the appendix.

A. Tree level
If the neutrino mass is generated at tree-level, one expects m ν ∝ v 2 /Λ, which for coefficients of O(1) give Λ ∼ 10 14 GeV for neutrino masses order 0.1 eV. The amplitude of the mass mechanism of 0νββ decay is proportional to . The d = 7 contribution is therefore favoured by a factor p F / m ν , but suppressed by (v/Λ) 3 . Inserting Λ ∼ 10 14 , the d = 7 amplitude should be smaller than the mass mechanism amplitude by a huge factor of order O(10 −27 ). However, this naive estimate assumes all coefficients in the operators to be order O(1). Since these coefficients are usually products of Yukawa (and other) couplings in the UV complete models, this is not necessarily the case in general and much smaller scales Λ could occur.
To discuss this in a bit more detail, we consider a particular example based on O 3 , decomposition #4, where two new fields, (1) a Majorana fermion ψ with the SM charge (SU(3) c , SU(2) L , U(1) Y ) = (1, 1, 0) and (2) a scalar S with (3, 2, 1/6), are introduced to decompose the effective operator, see Table III and Fig. 1. The Lagrangian for this model contains the following terms: Here, we have suppressed generation indices for simplicity. The first term in Eq. (6) will generate Dirac masses for the neutrinos. The Majorana mass term for the neutral field ψ (equivalent to a right-handed neutrino) can not be forbidden in this model. We will discuss first the simplest case with only one copy of ψ and comment on the more complicated cases with two or three ψ below. The contribution to 0νββ decay can be read off directly from the diagram in Fig. 1 on the left. It is given by With only one copy of ψ, the effective mass term contributing to 0νββ decay is m ν = (Y ν ) 2 e v 2 /m ψ and we can replace (Y ν ) e by m ν to arrive at the rough estimate of the constraint derived from the d = 7 contribution to 0νββ: Eq. (8) shows that the upper limit on the Yukawa couplings disappears as m ν approaches zero. When the masses are greater than roughly m ψ ≃ m S ∼ 10 TeV, the Yukawa couplings must be non-perturbative to fulfil the equality in Eq. (8). This implies that the mass mechanism will always dominate the 0νββ contribution for scales Λ larger than roughly this value, independent of the exact choice of the couplings.
We briefly comment on models with more than one ψ. As is well-known, neutrino oscillation data require at least two non-zero neutrino masses, while a model with only one ψ leaves two of the three active neutrinos massless. Any realistic model based on Eq. (6) will therefore need at least two copies of ψ. In this case Eq. (7) has to be modified to include the summation over the different ψ i and In this case, one still expects in general that limits derived from the long-range part of the amplitude are proportional to m ν . However, there is a special region in parameter space, where the different contributions to m ν cancel nearly exactly, leaving the long-range contribution being the dominant part of the amplitude. Unless the model parameters are fine-tuned in this way, the mass mechanism should win over the d = 7 contribution for all tree-level neutrino mass models.
The tables in the appendix show, that all three types of seesaw mediators appear in the decompositions of O 3 , O 4 and O 8 : ψ 1,1,0 (type-I), ψ 1,3,0 (type-III) and S 1,3,1 (type-II). In order to generate a seesaw mechanism, for some of the decompositions one needs to introduce new interactions, such as S † 1,3,1 HH, not present in the corresponding decomposition itself. However, in all these cases, the additional interactions are allowed by the symmetries of the models and are thus expected to be present. One then expects for all tree-level decompositions that the mass mechanism dominates over the long-range part of the amplitude, unless (i) the new physics scale Λ is below a few TeV and (ii) some parameters are extremely fine-tuned to suppress light neutrino masses, as discussed above in our particular example decomposition.

B. One-loop level
We now turn to a discussion of one-loop neutrino mass models. For this class of neutrino mass models, naive estimates would put Λ at Λ ∼ O(10 12 ) GeV for coefficients of O(1) and neutrino masses of O(0.1) eV. Thus, in the same way as tree neutrino mass models, the mass mechanism dominates over the long-range amplitude, unless at least some of the couplings in the UV completion are significantly smaller than O(1), as discussed next.
As shown in [51], there are only three genuine 1-loop topologies for (d = 5) neutrino masses. Decompositions of O 3 , O 4 or O 8 produce only two of them, namely Tν-I-ii or Tν-Iiii. We will discuss one example for Tν-I-ii, based on O 3 decomposition #2, see Table III and Fig. 2. The underlying leptoquark model was first discussed in [34,52], and for accelerator phenomenology see, e.g., [53]. The model adds two scalar states to the SM particle content, S(3, 1, −1/3) and S ′ (3, 2, 1/6). The Lagrangian of the model contains interactions with SM and the scalar interactions and mass terms: Lepton number is violated by the simultaneous presence of the terms in Eq. (9) and the first term in Eq. (10) [52]. Electro-weak symmetry breaking generates the off-diagonal element of the mass matrix for the scalars with the electric charge −1/3. The mass matrix is expressed as in the basis of (S −1/3 , S ′−1/3 ), which is diagonalized by the rotation matrix with the mixing angle θ LQ that is given as The neutrino mass matrix, which arises from the 1-loop diagram shown in Fig. 2, is calculated to be where N c = 3 is the colour factor. The loop-integral function ∆B 0 is given as with the eigenvalues m 2 1,2 of the leptoquark mass matrix Eq. (11) and the mass m d k of the down-type quark of the k-th generation. Due to the hierarchy in the down-type quark masses, it is expected that the contribution from m b dominates the neutrino mass Eq. (13).
and this gives roughly The constraint on the effective neutrino mass m ν 0.2 eV is derived from the combined KamLAND-Zen and EXO data [41], which is T 1/2 ≥ 3.4 × 10 25 ys for 136 Xe. The same experimental results also constrain the coefficient of the d = 7 operator generated from the Lagrangians Eqs. (9) and (10) as ǫ T R T R < ∼ 5.6 × 10 −10 (cf. Table I), which gives Therefore, for (λ S ) e1 (λ D ) 1e ≃ (λ S ) e3 (λ D ) 3e , the mass mechanism and the d = 7 contribution are approximately of equal size withM ≃ 750 GeV. Since m ν ∝M −2 , while ǫ O 3,#2 ∝M −4 , the mass mechanism will dominate 0νββ decay forM larger thanM ≃ 750 GeV, unless the couplings (λ S ) e1 (λ D ) 1e are larger than (λ S ) e3 (λ D ) 3e . We note that, leptoquark searches by the ATLAS [54,55] and the CMS [56-58] collaborations have provided lower limits on the masses of the scalar leptoquarks, depending on the lepton generation they couple to and also on the decay branching ratios of the leptoquarks. The limits derived from the search for the pair-production of leptoquarks are roughly in the range 650 − 1000 GeV [54][55][56][57][58], depending on assumptions.
The other 1-loop models are qualitatively similar to the example discussed above. However, the numerical values for masses and couplings in the high-energy completions should be different, depending on the Lorentz structure of the d = 7 operators, see also the appendix.

C. Two-loop level
We now turn to a discussion of 2-loop neutrino mass models. As shown in the appendix, in case of the operators O 3 and O 4 , 2-loop models appear only for the cases O 3a and O 4b . As explained in section II, these operators alone cannot give realistic neutrino mass models. We thus base our example model on O 8 . The 2-loop neutrino mass models based on O 8 are listed in Tab. V in the appendix. In this section, we will discuss decomposition #15 in detail, which has not been discussed in the literature before.
In this model, we add the following states to the SM particle content: With the new fields, we have the interactions which mediate O 8 operator, as shown in the left diagram of Fig. 3. Here, i runs over the three quark generations. While Y d i LαS k and Y u i ψH could be different for different i, for simplicity we will assume the couplings to quarks are the same for all i and drop the index i in the following. We will comment below, when we discuss the numerical results, on how this choice affects phenomenology. For simplicity, we introduce only one generation of the new fermion ψ, while we allow for more than one copy of the scalar S 3,2,1/6 . Note that, in principle, the model would work also for one copy of S 3,2,1/6 and more than one ψ, but as we will see later, the fit to neutrino data becomes simpler in our setup. The fermion ψ 2/3 mixes with the up-type quarks through the following mass term: Due to the strong hierarchy in up-type quark masses, we have assumed the sub-matrix for the up-type quarks in Eq. (21) is completely dominated by the contribution from top quarks. The mass matrix Eq. (21) is diagonalized with the unitary matrices V L and V R as and the mass eigenstates Ψ 2/3 i are give as where the index a for the interaction basis takes a ∈ {t, ψ}. The interactions are written in the mass eigenbasis as follows: The 2-loop neutrino mass diagram generated by this model is shown in Fig. 3. Using the formulas given in [59], one can express the neutrino mass matrix as . (26) Here N c = 3 is the colour factor and I(z k,i , r i , t i ) is the loop integral defined as The dimensionless parameters z k,i , r i , t i are defined as and loop momenta q and k are also defined dimensionless. Due to the strong hierarchy in down-type quark masses, we expect that neutrino mass given in Eq. (26) is dominated by the contribution from bottom quark. If we assume in Eq. (26) that all Yukawa couplings are of the same order, then the entries of the neutrino mass matrix will have a strong hierarchy: (m ν ) ee : (m ν ) µµ : (m ν ) τ τ = m e : m µ : m τ . Such a flavor structure is not consistent with neutrino oscillation data. Therefore, in order to reproduce the observed neutrino masses and mixings, our Yukawa couplings need to have a certain compensative hierarchy in their flavor structure.
Since the neutrino mass matrix, and thus the Yukawa couplings contained in the neutrino mass, have a non-trivial flavour pattern, these Yukawas will be also constrained by charged lepton flavour violation (LFV) searches. Here we discuss only µ → eγ which usually provides the most stringent constraints in many models. In order to calculate the process µ → eγ we adapt the general formulas shown in [60] for our particular case. The amplitude for µ → eγ decay is given by Here, ǫ α is the photon polarization vector and q β is the momentum of photon. Three different diagrams contribute to the amplitude for µ → eγ, which are finally summarized with the two coefficients σ R and σ L given by where Here, we have assumed that both the ψ −2/3 and the ψ −5/3 have the same mass M ψ . This neglects (small) mass shifts in the ψ −2/3 state, due to its mixing with the top quark. Due to the large value of M ψ , that we use in our numerical examples, this should be a good approximation. Note also, that the contribution from the top quark is negligible for those large values of M ψ used below. The functions F 1 (x) and F 2 (x) are defined in Eqs (40) and (41) in [60] as The branching ratio for µ → eγ can be expressed with the coefficients σ R and σ L as where Γ µ is the total decay width of muon. Later, we will numerically calculate the branching ratio to search for the parameter choices that are consistent with the oscillation data and the constraint from µ → eγ. Before discussing constraints from lepton flavour violation, we will compare the longrange contribution to 0νββ with the mass mechanism in this model. This model manifestly generates a d = 7 long-range contribution to 0νββ. The half-life of 0νββ induced by the long-range contribution is proportional to the coefficient ǫ V +A V +A which is expressed in terms of the model parameters as Here, we use the limit on ǫ V +A V +A from non-observation of 136 Xe 0νββ decay, see Table I. With one copy of the new scalar, the bound of Eq. (37) is directly related to the effective neutrino mass Eq. (26) and places the stringent constraint: where we have used the approximate relation with z k,1 = (m S k /m t ) 2 , r 1 = (m b /m t ) 2 , t 1 = (M W /m t ) 2 , and I(z k,1 , r 1 , t 1 ) ∼ 5 × 10 −2 for a scalar mass of m S = 10 TeV and M ψ ≃ 0.8 TeV. Note that this parameter choice is motivated by the fact that the model cannot fit neutrino data with perturbative Yukawa couplings with scalar masses larger than m S > ∼ 10 TeV. As one can see from Eq. (38), the long-range contribution to 0νββ clearly dominates over the mass mechanism in this setup. In short, this neutrino mass model predicts large decay rate of 0νββ but tiny m ν . This implies that, if future neutrino oscillation experiments determine that the neutrino mass pattern has normal hierarchy but 0νββ is discovered in the next round of experiments, the 0νββ decay rate is dominated by the long-range part of the amplitude. Recall that O 8 contains e c . This implies that the model predicts a different angular distribution than the mass mechanism, which in principle could be tested in an experiment such as Super-NEMO [61].
Note that, to satisfy the condition Eq. (38), cancellations among different contributions to m ν are necessary. This can be arranged only if we consider at least two generations of the new particles in the model (either the scalar S or the fermion ψ).
Here we discuss more on the consistency of our model with the neutrino masses and mixings observed at the oscillation experiments. Instead of scanning whole the parameter space, we illustrate the parameter choice that reproduces the neutrino properties and is simultaneously consistent with the bound from lepton flavour violation. To simplify the discussion we use the following ansatz in the flavour structure of the Yukawa couplings: with a dimensionless parameter y. With Eq. (40), the neutrino mass matrix Eq. (26) is reduced to where Λ is defined as and I is given as We introduce three copies of the new scalar S −1/3 k . The resulting mass matrix Eq. (41) has the same index structure as that of the type-I seesaw mechanism, and therefore, the matrix Λ can be expressed as following the parameterization developed by Casas and Ibarra [62].
Here,m ν is the neutrino mass matrix in the mass eigenbasis, and the mass matrix m ν is diagonalized with the lepton mixing matrix U ν as for which we use the following standard parametrization Here c ij = cos θ ij , s ij = sin θ ij with the mixing angles θ ij , δ is the Dirac phase and α 21 , α 31 are Majorana phases. The matrix R is a complex orthogonal matrix which can be parametrized in terms of three complex angles as Note that it is assumed in this procedure that the charged lepton mass matrix is diagonal. After fitting the neutrino oscillation data with the parametrization shown above, there remain y, Y uψH and the masses M ψ , m S k for k = 1, 2, 3 as free parameters. For simplicity, we assume a degenerate spectrum of the heavy scalars m S = m S k . In Fig. 4-(a), we plot the half-life T 0νββ 1/2 as a function of m ν 1 for fixed values of the coupling Y uψH = 0.6 and the masses M ψ = 800 GeV and m S = 10 TeV. The parameter y is taken to be 10 −3 , since this minimizes the decay rate of µ → eγ, as we will discuss below. We have used oscillation parameters for the case of normal hierarchy. The region enclosed by the red curves is d = 7 long-range contribution to 0νββ, and the blue curves correspond to the mass mechanism contribution only, which is shown for comparison. The gray region is already excluded by 0νββ searches, and for the model under consideration only the cyan region is allowed. As one can see from Fig. 4-(a), the total contribution to 0νββ is dominated by the d = 7 long-range contribution. Note that the mass mechanism and the long-range contribution are strictly related only under the assumption that Y uψH and Y dLαS k are independent of the quark generation i. This is so, because the 2-loop diagram is dominated by 3rd generation quarks, while in 0νββ decay only first generation quarks participate. If we were to drop this assumption and put the first generation couplings to Y u 1 ψH < ∼ 10 −2 ×Y u 3 ψH and Y d 1 LαS k < ∼ 10 −2 ×Y d 3 LαS k , the half-life for the long-range amplitude would become comparable to the mass mechanism, without changing the fit to oscillation data.
Note that non-zero Majorana phases are necessary to allow for cancellations among the mass mechanism contributions, so as to make m ν small as required by Eq. (38). In Fig. 4-(b), we plot the half-life T 0νββ 1/2 as a function of the scalar mass m S . Here we fixed the the decay rate versus m ν 1 (left) and m S (right). The gray region is the current lower limit in 0νββ decay half-life of 136 Xe. In the plot to the left the region between the red curves is the one allowed by the long-range contribution to the decay rate of 0νββ calculated scanning over oscillation parameters for the case of normal hierarchy and m S = 10 TeV. We also show the allowed region for the half-live for the mass mechanism as blue lines for comparison. The cyan region correspond to the parametric region where our model can be consistent with current 0νββ experimental data. In the plot to the right the red curve is the long-range contribution to the decay rate for the fixed oscillation parameters m ν 1 = 1.23 × 10 −3 eV , α 21 = 0, α 31 = π/2, s 2 23 = 1/2 and s 2 12 = 1/3 and the remaining oscillation parameters ∆m 2 31 and ∆m 2 21 fixed at their best-fit values for the case of normal hierarchy. oscillation parameters to m ν 1 = 1.23×10 −3 eV , α 21 = 0, α 31 = π/2, s 2 23 = 1/2 and s 2 12 = 1/3 and the remaining oscillation parameters ∆m 2 31 and ∆m 2 21 to their best-fit values for the case of normal hierarchy. The plot assumes that the matrix R is equal to the identity. The plot shows that the half-life increases to reach approximately T 0νββ 1/2 ∼ 10 26 yr for m S = 10 TeV. Now we discuss the constraint from lepton flavour violating process µ → eγ. In Fig. 5, we show Br(µ → eγ) as a function of the scalar m S and the parameter y for fixed values of the coupling Y uψH = 0.6 and the fermion mass M ψ = 800 GeV, which is the same parameter choice adopted in Fig. 4. These plots show that the current experimental limits on Br(µ → eγ) put strong constraints on the model under consideration. In Fig. 5-(a), we plot Br(µ → eγ) with different values of the parameter y = {10 −1 , 10 −2 , 10 −3 }. We have used again the parameters m ν 1 = 1.23 × 10 −3 eV, α 21 = 0, α 31 = π/2, s 2 23 = 1/2 and s 2 12 = 1/3 fixing the remaining oscillation parameters ∆m 2 31 and ∆m 2 21 at their best-fit values for the case of normal hierarchy. With the choice of y = 10 −1 , the entire region of m S is not consistent with the current experimental limits. On the other hand, we can easily avoid the constraint from µ → eγ by setting the parameter y to be roughly smaller than 10 −2 . Note that the curves with y = 10 −1 and y = 10 −3 do not cover the full range of m S . This is because the fit to neutrino data would require Yukawa couplings in the perturbative regime. (We define the boundary to perturbativity as at least one entry in the Yukawa matrix being smaller than √ 4π.) It is necessary to have smaller values of the parameter y to obey the experimental bound. This feature is also shown in Fig. 5-(b) where we plot the Br(µ → eγ) as a function of y with different values of the mass m S = {1, 5, 10} TeV. As shown, for y 10 −2 it is possible to fulfil the experimental limit, having the Br(µ → eγ) a minimum around y = 10 −3 . Because of the perturvative condition, the curves with m S = 5 TeV and m S = 10 TeV end in the middle of the y space. The reason for the strong dependence of Br(µ → eγ) on the parameter y can be understood as follows: As shown in Eq. (42) the Yukawa couplings Y dLαS k and Y eαψS k are related in the neutrino mass fit, but only up to an overall constant, 1 y . For values of y of the order of 10 −3 both Yukawas are of the same order and this minimizes Br(µ → eγ). If y is much larger (much smaller) than this value Y dLαS k (Y eαψS k ) becomes much larger than Y eαψS k (Y dLαS k ) and since the different diagrams contributing to Br(µ → eγ) are proportional to the individual Yukawas (and not their product) this leads to a much larger rate for Br(µ → eγ).
In summary, for all 2-loop d = 7 models of neutrino mass, which lead to O 8 , the longrange part of the amplitude will dominate over the mass mechanism by a large factor, unless there is a strong hierarchy between the non-SM Yukawa couplings to the first and third generation quarks. Such models are severely constrained by lepton flavour violation and 0νββ decay. We note again, that these models predict an angular correlation among the out-going electrons which is different from the mass mechanism.

IV. LEFT-RIGHT SYMMETRIC MODEL: d = 7 VERSUS d = 9 OPERATOR
Writing new physics contributions to the SM in a series of NROs assumes implicitly that higher order operators are suppressed with respect to lower order ones by additional inverse powers of the new physics scale Λ. However, there are some particular example decompositions for (formally) higher-order operators, where this naive power counting fails. We will discuss again one particular example in more detail. The example we choose describes the situation encountered in left-right symmetric extensions of the standard model.
Consider the following two Babu-Leung operators: O 8 can be decomposed in a variety of ways, decomposition #14 (see Table V) is shown in Fig. 6 to the left. The charged vector appearing in this diagram couples to a pair of right-handed quarks and, thus, can be interpreted as the charged component of the adjoint of the left-right symmetric (LR) extension of the SM, based on the gauge group SU(3) C × SU(2) L × SU(2) R × U(1) B−L . In LR right-handed quarks are doublets, Q c = Ψ3 ,1,2,−1/6 , the ψ 1,1,0 can be understood as the neutral member of L c , i.e. the right-handed neutrino, and the Higgs doublet is put into the bidoublet, Φ 1,2,2,0 . The resulting diagram for 0νββ decay is shown in Fig. 6 on the right. Fig. 6 gives a long-range contribution to 0νββ decay. We can estimate the size of ǫ O 8 from these diagrams: The first of these two equations shows ǫ O 8 for Fig. 6 on the left (notation for SM gauge group), the second for Fig. 6 on the right (notation for gauge group of the LR model). Here, g 1 and g 2 could be different, in principle, but are equal to g R in the LR model. v SM is the SM vev, fixed by the W -mass. In the LR model, the bi-doublet(s) contain in general two vevs. We call them v d and v u here and v 2 SM = v 2 d + v 2 u . In Eq. (49) only v u = v SM sin β, with tan β = v u /v d , appears. Note that we have suppressed again generation indices and summations in Eq. (49). We will come back to this important point below. Now, however, first consider O 7 . From the many different possible decompositions we concentrate on the one shown in Fig. 7. The diagram on the left shows the diagram in SM notation, the diagram on the right is the corresponding LR embedding. It is straightforward to estimate the size of these diagrams as: Arbitrarily we have called the 4-point coupling in the left diagram g 2 3 . In the LR model again the couplings are fixed to g L and g R . In the last relation in Eq. from a d = 9 operator. This a priori counter-intuitive result is a simple consequence of the decomposition containing the SM W L boson. Any higher-order operator which can be decomposed in such a way will behave similarly, i.e. 1/Λ 5 ⇒ 1/(Λ 3 v 2 SM ). 8 We note that in this particular example the contribution of O 7 is actually more stringently constrained than the one from O 8 . This is because O 8 leads to a low-energy current of the form (V + A) in both, the leptonic and the hadronic indices, i.e. the limit corresponds to ǫ V +A V +A . O 7 , on the other hand, leads to ǫ V +A V −A , which is much more tightly constraint due to contribution from the nuclear recoil matrix element [64], compare values in Table I.
We note that, one can identify the diagrams in Fig. 6 and Fig. 7 with the terms proportional to λ and η in the notation of [64], used by many authors in 0νββ decay. For the complete expressions for the long-range part of the amplitude, one then has to sum over the light neutrino mass eigenstates, taking into account that the leptonic vertices in the diagrams in Figs. 6 and 7 are right-handed. Defining the mixing matrices for light and heavy neutrinos as U αj and V αj , respectively, as in [64], the coefficients ǫ O 8 and ǫ O 7 of the d = 7 and d = 9 operators are then the effective couplings [64]: Orthogonality of U ej and V ej leads to 6 j=1 U ej V ej ≡ 0. However, the sum in Eq. (51) runs only over the light states, which does not vanish exactly, but rather is expected to be of the order of the light-heavy neutrino mixing. In left-right symmetric models with seesaw (type-I), one expects this mixing to be of order and in the LR model (right).
. In this case one expects the mass mechanism to dominate over both λ and η , given current limits on W L − W R mixing [65] and lower limits on the W R mass from LHC [66,67]. However, as in the LQ example model discussed previously in section III A, contributions to the neutrino mass matrix contain a sum over the three heavy right-handed neutrinos. In the case of severe fine-tuning of the parameters entering the neutrino mass matrix, the connection between the light-heavy neutrino mixing and m ν can be avoided, see section III A. In this particular part of parameter space, the incomplete 3 j=1 U ej V ej could in principle be larger than the naive expectation. Recall that the current bound on nonunitarity of U is of the order of 1 % [49]. For 3 j=1 U ej V ej as large as 3 j=1 U ej V ej ∼ O(10 −2 ) λ and/or η could dominate over the mass mechanism, even after taking into account all other existing limits. We stress again that this is not the natural expectation.
In summary, there are some particular decompositions of d = 9 operators containing the SM W or Higgs boson. In those cases the d = 9 operator scales as 1/(Λ 3 v 2 SM ) and can be as important as the corresponding decomposition of the d = 7 operator.

V. SUMMARY
We have studied d = 7 ∆L = 2 operators and their relation with the long-range part of the amplitude for 0νββ decay. We have given the complete list of decompositions for the relevant operators and discussed a classification scheme for these decompositions based on the level of perturbation theory, at which the different models produce neutrino masses. For tree-level and 1-looop neutrino mass models we expect that the mass mechanism is more important than the long-range (p / -enhanced) amplitude. We have discussed how this conclusion may be avoided in highly fine-tuned regions in parameter space. For 2-loop neutrino mass models based on d = 7 operators, the long-range amplitude usually is more important than the mass mechanism. To demonstrate this, we have discussed in some detail We also discussed the connection of our work with previously considered long-range contributions in left-right symmetric models. This served to point out some particularities about the operator classification, that we rely on, in cases where higher order operators, such as d = 9 (O 9 ∝ Λ −5 LNV ), are effectively reduced to lower order operators, i.e. d = 7 Our main results are summarized in tabular form in the appendix, where we give the complete list of possible models, which lead to contributions to the long-range part of the amplitude for 0νββ decay. From this list one can deduce, which contractions can lead to interesting phenomenology, i.e. models that are testable also at the LHC. corresponding decompositions are listed at the "Mediators" column. The symbols S and ψ represents the Lorentz nature of the mediators: S ( ′ ) is a scalar field, and ψ L(R) is a left(right)handed fermion. The charges of the mediators under the SM gauge groups are identified and expressed with the format (SU(3) c , SU(2) L ) U (1) Y . It is easy to find the contributions of the effective operators to neutrinoless double beta decay processes at the "Projection to the basis ops." column. The basis operators are defined as Here we explicitly write all the indices: α, β for lepton flavour, the lower (upper) I for 3 (3) of SU(3) colour, i, j, k, l for 2 of SU(2) left, ρ, σ for Lorentz vector, and a, b, c, d (ȧ,ḃ) for left(right)-handed Lorentz spinor. The lowest-loop contributions (i.e., dominant contributions) to neutrino masses are found at the columns "m ν ". We are mainly interested in decompositions (=proto-models) where new physics contributions to 0νββ can compete with the mass mechanism contribution mediated by the effective neutrino mass m ν . An annotation "w. (additional interaction)" is given in the column of "m ν @1loop" for some decompositions. This shows that one can draw the 1-loop diagram, putting the interactions that appear in the decomposition and the additional interaction together. The additional interactions given in the tables are not included in the decomposition but are not forbidden by the SM gauge symmetries, nor can they be eliminated by any (abelian) discrete symmetry, without removing at least some of the interactions present in the decomposition. For example, using the interactions appear in decomposition #11 of Babu-Leung operator #8 (see Tab. V), one can construct two 2-loop neutrino mass diagrams mediated by the Nambu-Goldstone boson H + , whose topologies are T 2 B 2 and T 2 B 4 of [59]. This also corresponds to the 2-loop neutrino mass model labelled with O 1 8 in [38]. However, to regularize the divergence in diagram T 2 B 4 , the additional interaction (Q c ) a Ii (iτ 2 ) ij (L) ia S I is necessary, and this interaction generates a 1-loop neutrino mass diagram. Consequently, this decomposition should be regarded as a 1-loop neutrino mass model. 9 We also show the 1-loop neutrino mass models that require an additional interaction with an additional field (second Higgs doublet H ′ ) with bracket. 10 The two contributions to 0νββ are compared in Sec. III with some concrete examples. The comparison is summarized at Tab. II. In short, the mass mechanism dominates 0νββ if neutrino masses are generated at the tree or the 1-loop level. When neutrino masses are generated from 2-loop diagrams, new physics contributions to 0νββ become comparable with the mass mechanism contribution and can be large enough to be within reach of the sensitivities of next generation experiments. However, the 2-loop neutrino masses generated from the decompositions of the Babu-Leung operators of #3 and #4 are anti-symmetric with respect to the flavour indices, such as the original Zee model and, thus, are already excluded by oscillation experiments. Therefore, if we adopt those decompositions as neutrino mass models, we must extend the models to make the neutrino masses compatible with oscillation data. In such models, the extension part controls the mass mechanism contribution and also the new physics contribution to 0νββ, and consequently, we cannot compare the contributions without a full description of the models including the extension. Nonetheless, it might be interesting to point out that decomposition #8 of the Babu-Leung #3 contains the tensor operator O ten.
3a (e, e), which gives a contribution to 0νββ and generates neutrino masses with the (e, e) component at the two-loop level. On the other hand, 2-loop neutrino mass models inspired by decompositions of Babu-Leung #8 possess a favourable flavour structure. This possibility has been investigated in Sec. III C with a concrete example.
There is another category of lepton-number-violating effective operators, not contained in the catalogue by Babu and Leung: operators with covariant derivatives D ρ . These have been intensively studied in Refs. [39,40,50]. The derivative operators with mass dimension seven are classified into two types by their ingredient fields; One is D ρ D ρ LLHH and the other is D ρ Lγ ρ e R HHH. With the full decomposition, it is straightforward to show that the tree-level decompositions of the first type must contain one of the seesaw mediators. Therefore, the neutrino masses are generated at the tree level and the mass mechanism always dominate the contributions to 0νββ. The decompositions of the second type also require the scalar triplet of the type II seesaw mechanism when we do not employ vector fields as  4,6,7,9, 10,13,15 0νββ and that A MM of the mass mechanism. When the neutrino mass is generated at the tree and one-loop level, the new physics scale Λ 7 must be sufficiently high to reproduce the correct size of neutrino masses, consequently, the long-range contributions A LR are suppressed and the mass mechanism dominates the contribution to 0νββ. As usual in such operator analysis, these estimates do not take into account that some non-SM Yukawa couplings, appearing in the ultraviolet completion of the operators, could be sizably smaller than one, which would lead to lower scales Λ 7 . Also, for loop model the scales could be overestimated, since they neglect loop integrals.
The neutrino masses generated at the two-loop level from the decompositions of the Babu-Leung #8 operator should be estimated with d = 7 LLHHHH † operator (as illustrated in sect. III C). In addition, they receive additional suppression from the lepton Yukawa coupling y ℓ β , which further lowers the new physics scale Λ 7 . Note that in particular for the 2-loop d = 7 models, as the concrete example in sect. III C shows, the estimate for A MM /A LR can vary by several orders of magnitude, depending on parameters. However, both the estimate shown here and the explicit calculation in sect. III C give numbers A MM /A LR ≪ 1 , such that the long-range contribution dominates always over the mass mechanism for these decompositions.
mediators, and the new physics contributions to 0νββ become insignificant again compared to the mass mechanism. In Ref. [50], the authors successfully obtained the derivative operator (e R c γ ρ Liτ 2 τ W ρ H ′ )(Hiτ 2 H ′ ) at the tree level and simultaneously avoided the tree-level neutrino mass with the help of a second Higgs doublet H ′ (1, 2) +1/2 and a Z 2 parity which is broken spontaneously. Here we restrict ourselves to use the ingredients obtained from  #3 (L α L β )(Q)(d R H) S(1, 1) +1 ψ L,R (3, 2) − 5   to 0ν2β are given as the combinations of the basis operators in the "Projection to the basis ops." column. The tensor operators O ten. play an important role in the long-range contribution.

Mediators
Projection to the basis ops. m ν @tree m ν @1loop m ν @2loop #1 (L α L β )(H)(Qu R ) S(1, 1) +1 S ′ (1, 2) + 1  decompositions and do not discuss such extensions. Within our framework, the derivative operators are always associated with tree-level neutrino masses. In this study, we have mainly focused on the cases where the new physics contributions give a considerable impact on the 0νββ processes. Therefore, we do not go into the details of the decompositions of the derivative operators.