Double beta decay and neutrino mass models

Neutrinoless double beta decay allows to constrain lepton number violating extensions of the standard model. If neutrinos are Majorana particles, the mass mechanism will always contribute to the decay rate, however, it is not a priori guaranteed to be the dominant contribution in all models. Here, we discuss whether the mass mechanism dominates or not from the theory point of view. We classify all possible (scalar-mediated) short-range contributions to the decay rate according to the loop level, at which the corresponding models will generate Majorana neutrino masses, and discuss the expected relative size of the different contributions to the decay rate in each class. Our discussion is general for models based on the SM group but does not cover models with an extended gauge. We also work out the phenomenology of one concrete 2-loop model in which both, mass mechanism and short-range diagram, might lead to competitive contributions, in some detail.

However, from the theoretical point of view it is not a priori clear, whether the mass mechanism gives indeed the dominant contribution to the double beta decay rate, and many models have been discussed in the literature where this might not be the case. The classical example appears in left-right (LR) symmetric extensions of the SM [9,10] and also in Rparity violating (R P / ) supersymmetric theories with both trilinear R P / [11,12] and bilinear R P / [13,14] terms. Furthermore, leptoquark models [15] and more recently models with colour octet scalars [16] or colour sextet diquarks [17][18][19] have been discussed.
1 m ν is defined as m ν = j U 2 ej m j , where the sum runs over all light neutrinos. This is equivalent to the (e, e) entry of the Majorana neutrino mass matrix in the basis where the charged lepton mass matrix is diagonal.
Given that there is such a large list of possible lepton number violating models, is it possible to determine which contribution to the 0νββ decay rate is the dominant one? -Perhaps, if double beta decay were to be observed in the next round of experiments and either KATRIN [20] or cosmological data [21][22][23] also find hints for a neutrino mass scale of the order of, say, somewhat larger than O(0.1) eV, one could claim on the basis of minimality that the mass mechanism m ν gives (at least) the most important contribution to the total decay rate. However, once upper limits on the total neutrino mass ( m ν ) placed from cosmology drop below the level of O(0.1) eV, the question becomes exceedingly difficult to answer.
In that case, from the experimental point of view, there remain only a few possibilities to make progress. For example, measurements of the angular correlation between the two electrons from 0νββ [10,24,25] or measuring double beta plus decays [26] offer the possibility to identify the Lorentz structure (equivalently the chiralities of the emitted electrons) of the LNV processes. However, realistically the SuperNEMO proposal [25] could only accumulate the necessary statistics to identify the Lorentz structure, if the half-live of 82 Se is below 10 26 ys, while there yet exists no experimental proposal with a sufficient sensitivity for 0νβ + /EC decays to make use of the ideas discussed in [26].
From the theoretical point of view, as discussed in [27,28], the amplitude of 0νββ decay can be divided into a long-range and a short-range part. Here, short-range means that all particles appearing in the diagrams for double-beta decay are heavier than the nuclear Fermi scale, i.e. O(0.1) GeV. Current limits from 0νββ decay then correspond to lower limits on the effective scale Λ LNV of lepton number violation, where g eff is some mean of the couplings appearing in the diagram and m S and m F are the masses of the fields that mediate the 0νββ process, see the next section for details. This mass scale is testable, at least in principle, at the LHC, and the combination of future LHC limits (or a possible discovery, to express it in a more optimistic way) and double beta decay data might allow to test many, but maybe not all, of the possible short-range diagrams that contribute to the decay rate [29,30]. In this paper we take a different approach and study the question, whether the mass mechanism is dominant or not, from a purely theoretical point of view. As described above, current and next generation 0νββ decay experiments will test LNV interactions at the TeV scale. Such TeV-scale LNV interactions, on the other hand, appear also in the context of radiative neutrino mass models. In other words, a new physics (short-range) contribution to 0νββ decay will always also produce a non-zero neutrino mass. In this paper we discuss the relation between possible models for short-range contribution to 0νββ decay and the neutrino mass-generation mechanism. Our study is based on the complete list of possible decompositions of the d = 9 (short-range) double beta decay operator given in [31,32]. The general decomposition list given in [31,32] is equivalent, in principle, to defining all models which can give a contribution to double beta decay, and the black box theorem, The theorem guarantees that, once 0νββ decay has been observed, Majorana neutrino masses will be generated, independent of the underlying model, at the latest at the four-loop level. By itself, this theorem does not guarantee that the mass mechanism is the dominant contribution to 0νββ decay.
see below, guarantees that all these models will produce Majorana neutrino masses at or below four-loop order. Our approach therefore basically consists in classifying all the possible models contributing to 0νββ decay with respect to the loop level at which they will generate Majorana neutrino masses. We can then discuss the expected size of the two contributions to 0νββ: (1) the d = 9 short-range contribution O d=9 and (2) the contribution from the neutrino mass mechanism m ν which is radiatively induced in the corresponding models and conclude model-by-model, which one of the two is expected to dominate. Given that the list of [31] is tree-level complete, our discussion is quite general and covers actually models of neutrino mass from tree-level models to 4-loop models. Before entering into the details of our work, let us briefly comment on the well-known relation between short-range 0νββ contributions and neutrino masses, i.e., the black box theorem [33]. 2 The theorem guarantees that, once 0νββ decay has been observed, neutrinos are Majorana particles. However, the black box theorem does not demonstrate that the mass mechanism dominates 0νββ, since it only guarantees neutrino masses at the level of four-loop. Obviously, a four-loop diagram, additionally suppressed by m 2 u m 2 d m 2 e /Λ 5 LNV , can produce only tiny neutrino masses, which are many orders of magnitude below of what is required to explain oscillation data [35]. Nevertheless, the black box theorem, together with the general decomposition of the d = 9 double beta decay operator published in [31], defines the basic idea of our current paper. Indeed we find that all "models" listed in [31] produce neutrino masses at or below the 4-loop order as demanded by the theorem.
We comment that our work has some overlap with [36,37] and [38]. Babu & Leung [36] have written down all SM invariant ∆L = 2 operators from dimension-5 (d = 5, the wellknown Weinberg operator [39]) to d = 11 and showed the relation between the effective operators and neutrino masses on the basis of black-box like loop diagrams. The authors of [36,37] discuss then possible ultra-violet completions for several example operators and give estimates for the scales Λ LNV , for which those operators can explain current neutrino data. The authors of [38] provide a systematic study of these operators, for one-loop and two-loop neutrino mass models, and discuss also which of these could possibly be tested at the LHC. However, our discussion differs from these papers in that we are mostly interested in double beta decay and its relation to neutrino mass.
We mention also the work of [40], which pursues the link between the short-range contribution to 0νββ and neutrino masses, but takes a different approach from ours. The authors focus on three types of LNV effective interactions which consist only of leptons and Higgs doublets, and list the models in which those LNV interactions simultaneously drive both the new physics contribution to 0νββ and neutrino masses. The main difference between [40] and our work is that it is assumed in [40] that new physics resides in the leptonic sector only.
In our classification, we also rediscover several models discussed in the literature previously, like for example leptoquark models [15], which can give 1-loop neutrino masses as discussed in [41] or 2-loop neutrino masses, as in the model of [19] or the one in [42]. We do not cover, however, the possible contributions from light sterile neutrinos. There exists a vast amount of papers on this subject in the literature already [43][44][45][46][47][48][49][50][51][52][53][54][55] and we have nothing new to add on this particular subject.
Since neutrino mass models must not only produce the correct absolute values of neutrino mass, but also reproduce the observed flavour structure of the neutrino mass matrix, one also has to pay attention to constraints from flavour physics observables. In [56] the authors applied the hypothesis of "minimal flavour violation" (MFV) to effective operators that contribute to 0νββ and found that the MFV assumption constrains the effective operators to be smaller than the detectable level. In this work, we do not adopt any such theoretical assumption on the flavour structures of the parameters in the models. Instead we simply consider bounds on lepton flavour violating observables as constraints on model parameters. We believe this to be the correct approach since any of the "exotic" contributions to 0νββ decay requires the introduction of new scalars, not present in the SM, with their own Yukawa interactions with SM fermions. Thus the whole concept of MFV is not very well funded in any of the models of interest for 0νββ decay.
A few disclaimers might also be in order here. Our analysis concentrates on the true d = 9 operator, i.e. it covers only the short range part of the 0νββ amplitude. Our results thus do not cover, for example, the long-range diagrams of R-parity violating SUSY [57,58] or leptoquark models [15]. Also, we limit ourselves to scalar exchange, thus models with a coupling between new scalars and the SM gauge bosons, such as [59,60] are not considered. Also, this restrictions implies that we do not cover models with an extended gauge group either, especially we do not discuss models with left-right symmetry. And, finally, the list of decompositions in [31] is complete only at tree-level. Thus, we do not consider cases in which the neutrino mass is generated at some higher loop level, while the 0νββ amplitude appears at one-loop order, as for example in the recent papers [61,62].
The rest of this paper is organised as follows. At the beginning of section II, as a preparation, we will summarise the main results of [31] and discuss some generalities useful for the latter parts of the paper. We will then discuss the classification of the different possible models and estimate in each case the relative size of the contribution from the mass mechanism and the short range part of the amplitude. In section III, we will discuss one concrete two-loop model of neutrino mass in more detail. Section IV summarises and discusses our main findings. Tables with lists of the different models, classified as described in section II are deferred to the appendix.

II. SETUP AND CLASSIFICATION
In this section, we classify neutrino mass models based on the decomposition of the d = 9 0νββ effective operators, according to the number of loops in the resulting neutrino mass diagrams. In each class, we will compare the size of the two contributions to the 0νββ: (i) d = 9 operator itself, and (ii) the mass mechanism m ν induced by the d = 9 operator. This classification therefore allows to identify those models, for which one can expect nonstandard contributions (beyond the ordinary mass mechanism) to be important for 0νββ decay. Note, that m ν , is equivalent to the e-e entry in the neutrino mass matrix, (M ν ) ee , in the basis where the charged lepton mass matrix is diagonal. Thus, in this section we concentrate on the comparison of short-range contributions to the size of this entry in (M ν ). Of course, a complete fit to all neutrino data will need to take into account also all other entries in (M ν ). A specific example, how this can be done and the additional constraints from both, oscillation data and lepton flavour violating decays, is discussed in section 3. In the discussion in this section, we always keep generation indices in the unknown couplings of the different models. Other indices could be constrained combining double beta decay with, for example, oscillation data. A specific example for this is worked out in section 3.

A. Generalities
The short-range double beta decay operator can be generated at tree-level via only two topologies shown in Fig. 2. The bosons (depicted with dashed lines in the diagrams) in these topologies could be either scalars or vectors, but we will consider only scalar exchange here. Assigning the outer fermions with either the left ("L") or the right ("R") chirality in all possible permutations allows to derive the complete list of "decompositions" (or proto-models) that can contribute to the 0νββ decay amplitude at tree level [31]. The fermion propagator in topology-I contains two terms, but in the short-range part of the amplitude the first term is suppressed relative to the second by a factor of p / /m ψ ≃ p F /m ψ , where p F is the typical Fermi momentum in the nucleus and the mass m ψ is suppossed to be larger than O(100) GeV. Considering then only decompositions which pick the mass term from the propagator results in a total of 135 possible decompositions for Topology-I (T-I in the following), while there are 27 decompositions in Topology-II (T-II), if we limit ourselves to scalar exchange [31,32]. For tables showing the different decompositions see the appendix. Babu & Leung [36] have listed ∆L = 2 operators from d = 5 to d = 11. Among the d = 9 operators in their list, the following five are relevant for double beta decay: Here, we have suppressed the indices of generation, SU(3) c , SU(2) L , and Lorentz spinor, which are contracted appropriately. In the list of decompositions shown in [31], there appears also one operator not given in [36], which is 3 As the black box theorem demonstrates, one can obtain the Weinberg operator by connecting the quark legs in these effective operators with the SM Yukawa interactions.  mass models with respect to the number of loops, which we discussed in the introduction, is modified from this naive expectation, once the decomposition of the operators are specified.
In fact, as shown below, many decompositions of the operators in Eq. (4) contain automatically the particle content (and interactions) such that neutrino masses are generated at lower order, i.e., both tree-level and 1-loop neutrino mass models are found. And, surprisingly, also the opposite case exists: If we restrict ourselves to decomposing the operators only with scalar and fermion mediators, none of the decompositions of the operator O 14 generates a genuine neutrino mass diagram at the 2-loop level. We will come back to this important point later in more detail. Genuineness is one of the key concepts in our classification method. The term, genuine n-loop neutrino mass model, is defined as the model in which the neutrino mass term is generated at the n-loop level and for which, simultaneously, diagrams with loop level lower than n are guaranteed to be absent, see [63] for more details.

B. Tree-level neutrino mass models
Let us start with a rather trivial example, illustrated in Fig. 3, in which the d = 9 contribution to 0νββ is related to the Majorana neutrino mass at the tree level. Here, as everywhere else in this paper, subscripts on fields denote their transformation properties (or charge in case of U(1) Y ) under the SM group, SU(3) c × SU(2) L × U(1) Y . A new scalar field is denoted by the symbol S, and a fermion field, which is understood as either a vector-like fermion or Majorana fermion, is denoted as ψ. Thus, ψ 1,1,0 has the same quantum numbers as a right-handed neutrino, while S 1,2,1/2 is equivalent to a (copy of) the SM Higgs doublet.
Here, the singlet fermion field ψ 1,1,0 is allowed to have a Majorana mass m ψ . The dot (·) denotes the anti-symmetric tensor (iτ 2 ) for SU(2) L . This Lagrangian generates an effective operator for a short-range contribution to 0νββ, which corresponds to O SR 1 in the notation of [28]. (Here we use the notation O SR i for the five relevant short-range operators, defined in [28], to distinguish them from the lepton number violating operators O j , j = 11, 12, 14, 19, 20 and "−".) m S 1,2,1/2 is the mass of S 1,2,1/2 . Following the method adopted in [1, 31] and using the experimental bound [3][4][5] T 0νββ 1/2 ( 136 Xe) > 1.6 · 10 25 [ys] (8) one finds the bound on the coefficient of the d = 9 operator as which can be interpreted as a constraint on the Yukawa coupling: .
If the scalar mediator S 1,2,1/2 acquires a vacuum expectation value (vev), the right diagram in Fig. 3, which has the same topology as the type-I seesaw mechanism, will contribute to neutrino mass(es): The experimental bound on the effective neutrino mass m ν 0.3 eV, which is found from Eq. (8) under the assumption of the mass mechanism being dominant, gives where we have used v SM ≃ 174 GeV. If S 1,2,1/2 is identified as the SM Higgs doublet H, as in the ordinary type-I seesaw, the constraint on Y Sνψ shown in Eq. (12) is obviously much stronger than Eq. (10). 4 In other words: The mass mechanism dominates the contribution to 0νββ decay by far, if we consider SM Higgs exchange. However, S 1,2,1/2 is not necessarily the SM Higgs, it could be an additional new state, such as appear, for example, in multi-Higgs doublet models. In this case, neutrino masses would still be generated through the type-I seesaw mechanism, with the vev S of the scalar : Three different 1-loop diagrams for neutrino mass [78], which can appear in one of the 0νββ decay decompositions. For a discussion see text.
S 1,2,1/2 independent of the SM vev v SM , such as occurs, for example, in the neutrinophilic neutrino mass model of [64][65][66][67]. For this case one finds that if the relaxed constraint holds, Eq. (10) becomes more stringent than Eq. (12), i.e., the short-range diagram will be the dominant contribution to the 0νββ decay amplitude in this case. If S 1,2,1/2 has exactly zero vev, in the literature often called the "inert doublet", we can no longer directly relate the the relative size of the d = 9 operator with the mass mechanism. The only conclusion one can derive in this particular case is the trivial constraint that the standard model Higgs coupling with L and ψ 1,1,0 must obey Eq. (12).

C. 1-loop models
As shown in [78] (see also [79]), there exist a total of four genuine 1-loop diagrams, which can contribute to the d = 5 Weinberg operator at the renormalizable level. 5 Three diagrams, shown in Fig. 4, can be related to 0νββ decay decompositions. 6 Here, we have added a "ν" 5 There are also three more non-genuine diagrams, discussed in [78], which can be understood as one-loop generated vertices for one of the three tree-level seesaws. 6 The remaining diagram Tν-1-i can be understood as opening-up of the quartic scalar vertex in Tν-3, by inserting an additional scalar.
to the naming conventions of [78], in order not to confuse the 1-loop neutrino mass diagrams with the double beta decay topologies.
Let us discuss the relation between the neutrino mass diagram Tν-1-ii, which is shown as the left-most diagram in Fig. 4, and the decomposition of the d = 9 0νββ diagram , which is classified with the ID number T-II-2 O 11 in [31], as an example of this class of models. The decomposition leads to a Lagrangian, which contains the following terms: The lepton number violation can then be assigned to be due to the presence of the coupling µ SSS . With this Lagrangian, Eq. (14), the effective d = 9 Lagrangian that contributes to 0νββ process as short-range effects is given as: The experimental bound, see Eq. (8), on 0νββ decay can then be interpreted again as an upper limit on the new leptonic Yukawa interactions as: With the interactions shown in Eq. (14), neutrinos acquire Majorana masses through the diagram Tν-1-ii. This neutrino mass generation mechanism through the leptoquark-Higgs coupling was first proposed in [80] and discussed in detail in [41,81].
Assuming the coupling µ SSS is smaller than the average of the leptoquark quark masses m LQ = (m S 3,2,1/6 + m S 3,1,−1/3 )/2, we can roughly estimate the neutrino mass as where m d k is the mass of the down-type quark (of generation k), which enters this estimation as the vertex at the left-upper corner of Tν-1-ii of Fig. 4. Applying the bound on the effective Majorana mass m ν < 0.3 eV (which is obtained from Eq. (8) with the assumption of the mass mechanism dominance) to Eq. (17), we have Assuming that the flavour structure of the new Yukawa interactions in Eq. (14) is not strongly hierarchical, one concludes that third generation quarks give the largest contribution to the neutrino mass. If the vev S 1/2 is as large as the SM Higgs vev, this constraint is more than six (three) orders of magnitude more stringent than Eq. (16) for k = 3 (k = 1). Note again that if the vev S 1/2 vanishes then, as in the tree-level case, the d = 9 contribution to 0νββ and the mass mechanism in this class of models are independent of each other. If the leptoquark mass is set to O(1) TeV and Yukawa couplings are taken to be O(0.1), the trilinear coupling µ SSS must be O(100) keV to reproduce O(0.05) eV of neutrino masses, which is the minimum value necessary to reproduce data on atmospheric neutrino oscillations. Such a small value of the coupling µ SSS can (obviously) be probed only in a 0νββ process dominated by the mass mechanism, since there is no other LNV process, for which experiments have even remotely comparable sensitivity.
To finish this discussion of the one-loop neutrino mass case, recall the possible two-loop contributions to neutrino masses in this model. As shown in the next subsection (and the list in the appendix), the model described with Eq. (14) generates neutrino masses at the two-loop level, even if the value of S 1/2 is identical with zero. However, the one-loop contribution discussed in this subsection can easily (even with a value of S 1/2 smaller than MeV) dominate over the two-loop diagrams, as can be seen from Eqs. (17) - (18). Therefore, we classify this type of the models separately from the genuine two-loop models, which will be defined and explained in detail in the next subsection.
Very similar arguments can be applied to the other two topologies, Tν-1-iii and Tν-3, shown in Fig. 4. Such one-loop neutrino mass models appear in quite a large number of decompositions. We give the complete list of this class of models in Table II in the appendix, together with the additional interaction that is required (and is allowed by the SM gauge symmetries) to generate the corresponding one-loop diagram.

D. 2-loop models
The effective operators O 11 , O 12 , and O 14 contain two lepton doublets, and thus naively one expects them to generate neutrino masses at 2-loop level, 7 when their quark legs are connected with the SM Yukawa interactions. However, this naive picture has to be modified, once the possible decompositions of the effective operators are taken into account.
In this subsection, we will first demonstrate why some of the decompositions of the operators O 11 , O 12 , and O 14 do not generate neutrino masses genuinely at the 2-loop level. Here, we use the terminology genuine 2-loop diagrams, following [63], to imply that the 7 When the two lepton doublets are anti-symmetric in SU (2) L , an additional loop is necessary to obtain a neutrino mass term, i.e., the resulting neutrino mass diagram contains three loops, which correspond to O 11a , O 12b and O 14a in the list of [37]. In 0νββ decay only the operators O 11b , O 12a and O 14b , which contain symmetric pieces in SU (2) L , can appear. These are called O 11 , O 12 and O 14 for brevity here and in [31]. corresponding model (or decomposition in our case) generates neutrino masses at 2-loop order and that no lower order diagram exists. Let us first discuss the decomposition of O 14 . This operator contains: Here we explicitly wrote the 2-component spinor indices for the fermions: dotted for a righthanded field (complex conjugate of left-handed field), undotted for a left-handed field. As everywhere else in this paper, we restrict the discussion to fermions and scalars as mediators.
For O 14 this implies that the spinor indices on Q and u R must be contracted for the effective operator being a Lorentz scalar. There are only two choices to assign these two quarks to the outer legs of a d = 9 tree diagram: (i) They form a Yukawa interaction with a scalar mediator, i.e., (u R Q · S). This is shown as the upper diagram of Fig. 5. Or: (ii) Each of these quarks forms a Yukawa interaction with a fermion mediator ψ and one of the scalar mediators, i.e., (u R ψS)(ψQS ′ ). This is shown as the lower diagram of Fig. 5. When the loops are closed, as shown on the right of Fig. 5, via the SM Yukawa interaction y u (Q·H * u R ), the resulting quark loop is divergent. In other words, these loops are infinite corrections for (i) a mass term mixing the scalar mediator and the SM Higgs doublet m 2 SH SH * , or (ii) a corresponding term mixing two scalar mediators m 2 SS ′ SS ′ . Therefore, the original tree- level Lagrangian generating these operators must contain these scalar mass terms as counter terms for the infinities, and the quark loop appearing in neutrino mass diagrams must be substituted with those scalar mass terms. The inner loop in these classes of neutrino mass diagrams actually corresponds to an infinite (1-loop) correction to the scalar quartic interaction λ SS ′ HH SS ′ HH. In the diagram on the right, the "inner" loop involving S generates a Yukawa interaction y QψHψ Q · H. Again, this correction is infinite, thus requiring a tree-level counter term which must be contained in the original Lagrangian. Given these additional (but required) interactions, the models contain neutrino mass diagrams at 1-loop order. Note that, the left diagram in Fig. 6 corresponds to Diagram (A) of Fig. 14 in [38] (with appropriate Higgs insertions) and also Diagram (c) of Fig. 5 in [82]. The right diagram is Diagram (B) in [38] and (d) in [82]. Quite a number of possible decompositions of O 11 and O 12 follow this pattern and are thus listed in the tables in the appendix as 1-loop models, together with the additional-but-necessary interactions.
After filtering out all decompositions that result in non-genuine 2-loop neutrino mass diagrams, we have found that for all remaining decompositions there are only three types of genuine 2-loop diagrams, all of them based on O 11 . These are shown in Fig. 7. The naming scheme in this figure follows [63]. In the appendix, we present the complete list of the genuine 2-loop neutrino mass models and specify the class of neutrino mass diagrams, into which each model falls. In Table IV, two of the decompositions based on the BL operator O 19 are also listed. These appear in the table for two-loop models due to the fact that the intermediate fermion is of Majorana type, i.e., for these decompositions, the "asymmetric" following the classification of [63]. The left diagram (CLBZ-1) was discussed in the context of a neutrino mass model in [19], and corresponds to the diagram in Fig. 10 in [38] and (e) in [82].
The diagram (PTBM-1) in the middle was discussed in [83], and corresponds to Diagram (D2) in Fig. 14 in [38] and (f) in Fig. 5 in [82]. The diagram on the right corresponds to Diagram (C) in [38] and (f) in Fig. 5 [82]. A model based on this diagram will be discussed in Sec. III.
operator O 19 ∝Le R is always accompanied by the "symmetric" operator O 11 ∝LL, and the associated operator generates neutrino masses at the two-loop level. The catalogue of the effective operators appearing with their associated operators is given in the tables of [31]. Among these three genuine diagrams, neutrino mass models based on the CLBZ-1 and the PTBM-1 diagrams have already been studied in the context of the decomposition of the d = 9 operators [19,83]. Therefore, we will discuss a 2-loop neutrino mass model that is associated with the remaining possibility, i.e., the PTBM-4 diagram. Here, as in the previous subsections, we compare the d = 9 contribution to 0νββ with the mass mechanism contribution and postpone the detailed discussion on phenomenology of this model till Sec. III.
The example we choose is based on T-I-4-ii-b, O 11 . The Lagrangian for this decomposition contains the terms We use the notationˆ S 6,3,1/3 = ( S 6,3,1/3 ) X (T6) X IJ andψ 6,2,1/6 = (ψ 6,2,1/6 ) X (T6) X IJ . The tensors T 6 and T6 in the SU(3) c are given in [31]. Here, τ is the Pauli matrix vector for a triplet of SU(2) L . The effective d = 9 operator resulting from this Lagrangian can be written with the following linear combination of the basis operators O SR i∈{1-5} of the short-range contributions to 0νββ decay as: and the experimental bound Eq. (8) constrains a combination of the coefficients On the other hand, the neutrino mass generated from the 2-loop diagram based on the effective operator Eq. (21) also contributes to 0νββ through the mass mechanism. The size of 2-loop neutrino mass diagram can be roughly estimated as [36][37][38]63] Applying the experimental bound from 0νββ to Eq. (23) and substituting N c = 6 because of the colour sextet combination in the loop, we can place the bound on the couplings of the third generation quarks: .
As this rough estimation shows, the mass mechanism and the short-range part of the amplitude give similar contributions to 0νββ decay. Note that assuming flavour democratic Yukawa couplings, the constraints Eq. (22) and Eq. (24) become equally strong if the mass scale of the new particles is taken to be roughly ∼ 300 GeV. For larger mass values, the short range contribution can dominate only if Yukawas with index "3" are smaller than those with index "1", otherwise the mass mechanism dominates. A more detailed discussion using the full expression for the two-loop neutrino mass integral will be presented in Sec. III.

E. 3-loop models
From the point-of-view of effective operators, the Babu-Leung operators O 19 and O 20 require three SM Yukawa interactions to generate neutrino masses: two quark Yukawa interactions and one charged-lepton Yukawa interaction, to convert e R in the effective operators to L for a neutrino mass. This fact leads us to three-loop neutrino mass models. However, some of the possible decompositions of O 19 and O 20 contain the ingredients to generate neutrino masses at a level lower than three-loop. In such a case, the lower loop contributions can easily dominate neutrino masses and make the contribution from a three-loop diagram sub-dominant. This can happen for two reasons. First, there are decompositions based on O 19 or O 20 , in which O 11 necessarily also appears. We call this "associated" operators and classify those decompositions in the class corresponding to those lower loop levels, which usually will dominate over the 3-loop contribution. corresponds to the 1-loop generation of a certain vertex. We will discuss this case in a bit more detail.
Take the examples of the decompositions T-I-4-ii-a and T-I-4-ii-b, both O 20 . The Feynman diagrams are given in Fig. 8, and show T- graphically. As we will see, despite the similarity between these two cases, T-I-4-ii-b will lead to a genuine 3-loop model, while T-I-4-ii-a will not. Consider the examples of 3-loop diagrams for these two decompositions shown in Fig. 9. First of all, note that the loop diagrams shown in Fig. 9 should be understood as examples only, because there might be more than one diagram contributing to the full neutrino mass matrix for each decomposition. In the diagram for T-I-4-ii-a (left), one sees that the innermost loop effectively generates the vertex u R ψ 3,2,7/6 H † at 1-loop order. This one-loop sub-diagram is infinite and, therefore, a tree-level counter term is necessarily to be included in the La- grangian to absorb the infinity. In fact, the quantum numbers of the particles involved in the loop are such that the necessary vertex actually cannot be forbidden at tree-level by the SM gauge symmetry. This tree-level coupling has a value that is not fixed by the 0νββ decay amplitude, but the 2-loop (d = 7) diagram that results from this coupling, see Fig. 10, can easily dominate over the 3-loop diagram. 8 We have classified therefore all these cases of O 19 and O 20 , where such a tree-level vertex is allowed, as "2-loop d = 7" (Tab. V) in Appendix. The loop diagram based on the decomposition T-I-4-ii-b, which is shown as the right diagram of Figure 9, does not contain such an inner loop and, thus, such a construction is not possible for this decomposition. Decompositions of this type are therefore classified as genuine 3-loop models in the appendix.
As above, we use the notationŜ 6,1,4/3 = (S 6,1,4/3 ) X (T6) X IJ andψ 6,1,1/3 = (ψ 6,1,1/3 ) X (T6) X IJ . Together, the Y eψS , Y QψS , and Y LdS terms necessarily violate lepton number by two units. All generation indices in the couplings in Eq. (25) have been suppressed for simplicity. The contribution to the neutrino mass matrix can be roughly estimated as where Λ LNV ≃ m S 6,1,4/3 ≃ m S 3,2,1/6 ≃ m ψ is the mass scale of the heavy states, which is typically taken to be TeV. N c is a colour factor. Here, we assumed that all the SM fermion masses are much smaller than Λ LNV . Putting all the Yukawa couplings in Eq. (26) equal to unity and Λ LNV = 1 TeV and N c = 6 (for a colour sextet combination), one finds 9 This implies that the mass mechanism contribution to 0νββ is guaranteed to be subdominant in this class of models. Also, Eq. (27) shows that 3-loop models can potentially explain neutrino oscillation data only if all of the involved Yukawa couplings are set to be O(1). Thus, we expect such models to be quite constrained from upper limits on flavour violating decays of charged leptons. We will not discuss this class of models in more detail here, since their detailed phenomenology is outside the scope of this paper. The effective d = 9 operator resulting from the Lagrangian Eq. (25) can be written with the following linear combination of the basis operators O SR i∈{1-5} of the short-range contributions to 0νββ decay as: and the experimental bound Eq. (8) constrains a combination of the coefficients to be: The difference in the short-range bounds, Eq. (22) and Eq. (29), is due to the different values of nuclear matrix elements entering the transition operator. All other three-loop models will have constraints similar to the ones discussed here. They are listed in Table VI in the appendix. 9 Using m ν ≤ 0.3 eV, we can formally write the constraint on the Yukawa couplings in the form of: which is much worse than even the trivial constraint derived from perturbativity.
. To the right: four-loop d = 9 neutrino mass, see text.

F. 4-loop models
Finally, all operators O − = e R e R u R d R u R d R , with exception of decomposition T-I-5-i (see Table I in the appendix), will lead to four-loop neutrino mass models. The simplest possibility to construct a four-loop diagram for these operators is to use a SM charged current interaction. We estimate that this gives the dominant contribution to the neutrino mass. Here we show an example of the decompositions of the 0νββ decay operator O − in Fig. 11, which is based on decomposition T-I-3-ii (u R u R )(d R )(d R )(e R e R ). The four-loop neutrino mass diagram based on this decomposition is also shown on the right. Taking the limit m ψ 3,1,5/3 ∼ m S 6,1,4/3 ∼ m S 1,1,2 ≫ m W , m t one can estimate the order of magnitude of this four-loop diagram, which is, The expression Eq. (30) shows that this four-loop contribution would yield only (m ν ) τ τ ∼ O(10 −10 ) eV for m ψ 3,1,5/3 ∼ m S 6,1,4/3 ∼ m S 1,1,2 ∼ 1 TeV, even when choosing all SM fermion masses to be third generation. Since this is obviously many orders of magnitude below the values of neutrino masses required from oscillation experiments, models of this category by themselves cannot be considered realistic. Of course, neutrinos could be quasi-Dirac particles, explaining oscillation data by Dirac mass terms (using additionally introduced right-handed neutrinos), while 0νββ decay is dominated by the short-range diagrams such as the one shown in Fig. 11. However, constraints on Yukawa couplings will be similar to those derived in the previous subsections in Eq. (21) and Eq. (28), with the exact value depending on the decomposition under consideration. All four-loop cases are listed in Table VII in the appendix.

III. A CONCRETE 2-LOOP EXAMPLE
In this section we will discuss one concrete genuine 2-loop neutrino mass model in some more detail. The example we choose is based on the decomposition T-I-4-ii-b of the Babu- Leung operator O 11 , which has not been discussed in the literature before. However, all 0νββ decompositions that generate 2-loop neutrino masses behave quite similarly, in what concerns fits for neutrino oscillation data and constraints from lepton flavour violation searches. Thus, most of the discussion presented below can be applied qualitatively also to all other 2-loop decompositions. Any model of neutrino mass must not only generate the correct neutrino mass scale, but also be able to explain the observed neutrino mixing angles. For a recent update of all oscillation data, see, for example, [84]. In addition, since the neutrino mass matrix has a non-trivial flavour pattern, one also expects that low-energy models 10 of neutrino mass are constrained by charged lepton flavour violation (LFV) searches. Here we will discuss only µ → eγ, since the experimental upper limit on this process provides usually the most stringent constraints in many models. We note that the authors of [83] present a 2-loop model, which corresponds to the decomposition T-I-5-i and discuss also the constraints from other LFV searches, which we expect are very similar in our example.
Below we will discuss two variations of the model based on T-I-4-ii-b. First (in Sec. III B), we introduce only one copy of the exotic fermion ψ 6,2,−1/6 for simplicity. Next (in Sec. III C), we will allow to have three copies of these fermions, which allows to fit also quasi-degenerate neutrinos.

A. General formulas for neutrino masses and µ → eγ
The Yukawa part of the Lagrangian describing the interactions between the exotic diquark, S 6,3,1/3 , the leptoquark, S 3,2,1/6 , and the coloured vector-like fermion, ψ 6,2,−1/6 can be written as: Here, i, j are generation indices for quarks, we use Greek indices for lepton generations and k runs over the number of copies of ψ 6,2,−1/6 . This Lagrangian generates a 2-loop diagram which corresponds to PTBM-4 according to classification by [63]. Following the general formulas from [63], the neutrino mass matrix can be expressed as: where N c is a colour factor, with N c = 6 for this model. Summation over all flavour indices i, j, k is implied. F (m ψ k , m S 3,2,1/6 , m d i , m S 6,3,1/3 , m d j ) is a loop integral defined as: .
Due to the strong hierarchy in down-type quark masses, the integral in Eq. (33) is completely dominated by the contributions from bottom quarks, unless the couplings Y ij QQS , Y ki ψdS and Y βj LdS follow an equally strong inverse hierarchy. We have thus taken into account only the contributions from bottom quark exchange in our numerical evaluation. Since it is convenient to rewrite Eq. (33) in terms of dimensionless parameters, we define z, r and t k as Rescaling the loop momenta, the integral can then be written as: This integral has been analytically calculated several times in literature. We follow the procedure outlined in [63], based on the calculations of [83]. We will fit the neutrino mass calculated with Eq. (32) to neutrino oscillation data. The discussion depends on the number of copies of the fermion mediator ψ 6,2,−1/6 ; as mentioned above we will discuss two different scenarios in the following subsections. The rate of the LFV process µ → eγ has also been calculated several times in literature. We adapt the general formulas shown in [85] for our particular case. The amplitude for µ → eγ decay is given by where e is the electric charge, ǫ α is the photon polarization vector, q β is the momentum of photon, and σ αβ ≡ (i/2)[γ α , γ β ]. There are two contributions to the coefficient σ R in the model we are discussing; one is the one-loop diagram with the diquark and the exotic fermion, the other is that with a bottom quark and the leptoquark. The total σ R is given by where x k ≡ . The functions F 1 (x) and F 2 (x) are defined as which are presented in Eqs. (40) and (41) in [85]. The branching ratio for the µ → eγ process, neglecting the electron mass, can then be expressed with the coefficient σ R as where α is the fine-structure constant.
B. One generation of ψ 6,2,1/6 The analysis presented in this section uses very similar methods to the one in ref. [16], where double beta decay and LFV is discussed in a 1-loop neutrino mass model containing colour octets. We will first consider a variant of the model, in which there is only one copy of the fermion mediator ψ 6,2,1/6 . The expression for the neutrino mass matrix in this case is given by suppressing the index for ψ in Eq. (32), which gives where Since det(m ν ) = 0 in this case, this version of the model can fit only to the hierarchical neutrino mass spectra (both of the normal and the inverse type), but not to the degenerate spectrum. 11 The eigenvalues of Eq. (41) can be easily found to be: for normal hierarchy (inverted hierarchy). In Fig. 13, we give typical values for the common factor F , which are calculated with the assumption of a nearly degenerate spectrum of heavy particles with the mass scale M eff ≡ m ψ ≃ m S 3,2,1/6 ≃ m S 6,3,1/3 and (Y QQS ) 33 = (Y ψdS ) 3 = 1. From Eq. (43) and Fig. 13, one can estimate the constraints from neutrino masses on the size of the Yukawa couplings. In order to reproduce the neutrino mass suggested by atmospheric neutrino oscillation (m ν 3 ∼ 0.05 eV), keeping the common mass scale M eff at 1 TeV, the Yukawa couplings Y α LψS and Y β LdS must be set typically to O(10 −2 ).
Although the eigenvectors of Eq. (41) can be calculated analytically, numerical exercises might be more helpful to grasp phenomenological aspects of the model. In the following, we will generate random sets of Yukawa couplings (Y LψS ) α and (Y LdS ) β3 under the condition that they reproduce the latest neutrino oscillation data [84] within 3 σ C.L. We will only show plots with the Yukawa couplings that fit the normal hierarchical neutrino spectrum, because plots for the inverse hierarchical case look qualitatively similar. Let us start the discussion with double beta decay. The half-life of 0νββ induced by the Majorana mass of neutrino is proportional to the inverse-square of the effective neutrino mass: T 0νββ For the normal hierarchy case, the effective mass is roughly given as (m ν ) ee ∼ s 2 12 ∆m 2 21 ∼ 3 × 10 −3 eV, which results in half-lives of the order of T 0νββ  The short range-contribution due to the d = 9 operator (cf. Eq. (21)) is proportional to the following combinations of the parameters: i.e., while the neutrino mass matrix is dominated by Yukawa couplings of the third quark generation, double beta decay is sensitive only to the Yukawa couplings that couple to the first generation quarks. To discuss the relation between these two contributions to 0νββ, we introduce a scaling factor i.e., η 31 = 1 corresponds to quark flavour universality in the Yukawa couplings. In Fig. 14, we calculate half-lives induced from the short-range contribution with randomly generated Yukawa couplings, assuming different values of η 31 ∈ {1, 5, 10, 50}. Taking η 31 = 1, we find quite long half-lives, too large to be measured in realistic experiments. On the other hand, with η 31 = 10, we find a lower limit on M eff , which is approximately M eff 400 GeV. This is still not competitive with leptoquark searches at the LHC, which places constraints on the masses of leptoquarks at m S 3,2,1/6 ∼ (600 − 1000) GeV (depending on generation) already in the first run [87-89] Thus, it is reasonable to conclude that, as for the mass mechanism contribution, the short-range contribution to the half-life is also expected to be too long to be measured in the near future in this variant of the model, unless η 31 is very large (i.e., for highly inverse hierarchical Yukawa couplings in terms of the quark generations). Finally, we discuss briefly the LFV process µ → eγ. While the neutrino mass matrix is proportional to the combination of the Yukawa couplings ( Fig. 15, we show Br(µ → eγ) for two different choices of the set of (Y QQS ) 33 and (Y ψdS ) 3 , as a function of M eff , assuming again that the mass spectra of heavy particles are nearly degenerate for simplicity. With the choice (Y QQS ) 33 = (Y ψdS ) 3 = 10 −2 the LFV process can place a bound on M eff of roughly M eff > ∼ TeV. However, the bound depends strongly on the exact choice of the remaining Yukawa couplings Y LψS and Y LdS . On the other hand, the LFV process can exclude only few parameter points in the case of (Y QQS ) 33 = (Y ψdS ) 3 = 10 −1 , and no useful limit on M eff can be derived.
As we have seen in this subsection, this variant of the model can reproduce oscillation data without running into conflict with LFV searches. However, it is interesting to note that the (Y QQS ) 33 interaction with the size required for reproducing neutrino masses will result in sizeable decay rates of the diquark into third generation quarks (both tops and bottoms), which should be testable at the LHC.
C. Three generations of ψ 6,2,1/6 Next, we examine the model with more than one copy of the fermion mediator ψ 6,2,−1/6 , which can fit not only hierarchical neutrino mass spectra, but also a quasi-degenerate spectrum. Here, we introduce three copies of ψ 6,2,−1/6 , motivated by the observed generations of SM fermions.
To simplify the following discussion, we adopt the following ansatz in the flavour structure of the Yukawa couplings: 12 With this ansatz, all the flavour structure relevant to phenomenology can be represented with only one vector (apart from a possible normalization factor y). The neutrino mass matrix can then be cast into the form: where the Λ is defined as andÎ is given asÎ Comparing Eq. (48) with the neutrino mass and mixing matrix, we can find the direct relation between Λ and the measured neutrino data. Following the procedure originally developed by Casas and Ibarra for seesaw type-I [91], we parametrize Λ as Here,m ν is the matrix of eigenvalues of m ν , which is diagonalized with the neutrino mixing matrix U ν via for which we use the following standard parametrization c ij = cos θ ij , s ij = sin θ ij with the mixing angles θ ij , δ is the Dirac phase and α 1 , α 2 are Majorana phases. Finally, R is a complex orthogonal matrix which satisfies the condition R T R = 1. We use the following parametrization for the R matrix in terms of three complex angles θ 1 , θ 2 , and θ 3 as After fitting the neutrino oscillation data with the parametrization shown above, there remain y, (Y QQS ) 33 and the masses m S 6,3,1/3 , m S 3,2,1/6 , m ψ k as free parameters. For the calculation of the short-range contribution to the 0νββ decay, we also have the parameter η 31 . For simplicity, we set y = 1 and assume again a nearly degenerate spectrum for heavy particles, which is parameterized with M eff . We can then calculate half-lives T 0νββ 1/2 for both neutrino mass mechanism and the short-range contribution, as a function of m ν 1 , M eff and η 31 . In Fig. 16, we fixed the oscillation parameters s 2 13 , ∆m 2 31 and ∆m 2 21 at their best-fit values, while s 2 23 = 1/2 and s 2 12 = 1/3 and set δ as well as the Majorana phases to zero, just to sketch out some phenomenological aspects of this example. Each panel shows the half-life of the short-range contribution for 0νββ decay of 136 Xe as a function of η 31 (top panel), M eff (middle panel) and m ν 1 (bottom panel). In each panel, we examine several choices for the remaining parameters, which are explained in the figure caption. The corresponding halflives induced from the mass mechanism are also indicated. As shown in Fig. 16, half-lives can vary over many orders of magnitude with the choice of parameters. The amplitudes induced from the mass mechanism becomes the same order as that from the short-range, when η 31 ∼ 2.7 (6.5) for M eff = 0.5 TeV (1 TeV). As in the case with only one generation of ψ 6,2,1/6 , the mass mechanism dominates the 0νββ, if the ratio η 31 is taken to be unity and the heavy mass scale M eff is given at the typical LHC search sensitivities. However, since the three-generation case can fit the quasi-degenerate neutrino spectrum, 0νββ decay half-lives can be much shorter than in the one generation case and can saturate the experimental bound.
We now turn to Br(µ → eγ). Again, as in the one generation case, the neutrino mass matrix depends on Yukawa couplings, but is not directly related to Br(µ → eγ). Therefore, we have always the freedom to adjust (Y ψdS ) k3 and (Y QQS ) 33 so as to fit the neutrino masses. The other Yukawa couplings are then fixed by the neutrino data (and the choice of M eff ), and we can use them to calculate Br(µ → eγ). Fig. 17 shows some examples with a value of M eff = 1 TeV. The plots show that constraints from Br(µ → eγ) can be easily fulfilled. For this choice of M eff , only if both (Y QQS ) 33 and (Y ψdS ) 13 are set to order O(10 −2 ) or lower, the predicted Br(µ → eγ) can saturate the experimental bound.

IV. CONCLUSIONS AND DISCUSSION
We have discussed the relation between the d = 9 short-range contributions to the 0νββ decay amplitude with neutrino mass models. All contributions to 0νββ decay violate lepton number and, therefore, generate also Majorana neutrino masses. We have classified all possible (scalar-mediated) short-range contributions to the decay rate according to the loop level, at which the corresponding models will generate Majorana neutrino masses. Possibilities range from tree-level to 4-loop neutrino masses. For each case we have discussed one example briefly and given estimates of the typical constraints imposed by both the short-range contribution and the mass mechanism. Generally, one expects that for models with treeor 1-loop neutrino masses, the short-range 0νββ decay amplitude will be sub-dominant to the mass mechanism. For 2-loop models short-range 0νββ decay amplitude and mass mechanism can be comparable, while for 3-loop and 4-loop models the short-range part of the amplitude will dominate.
We have also discussed one particular example of a 2-loop model in more detail. Here, we have shown different parts of parameter space where mass mechanism or short-range amplitude dominant can each be dominant. In the study, we have taken recent neutrino oscillation data and constraints from LFV experiments into consideration.
In the appendix we give the full list of decompositions, classified according to our scheme, in tabular form.

V. APPENDIX
Here, we give tables in which all possible scalar short-range decompositions are classified according to the loop level at which they will generate neutrino masses. The decompositions which generate neutrino masses at tree, 1-loop, 2-loop, 3-loop and 4-loop level are listed in Tables I to VII. The identification number given to each decomposition is defined in [31,32]. The notation T-I and T-II refers to the two possible topologies of the decompositions of d = 9 0νββ effective operators, and the BL number is the Babu-Leung classification of the effective neutrino mass operator, given in [36]. The columns "Add. Int." specify additional interactions, with respect to those appearing in the decomposition. While they do not appear directly in the 0νββ diagram, these additional interaction can not be forbidden by any symmetry, without forbidding the corresponding 0νββ decay decomposition at the same time. Once present, they generate a neutrino mass diagrams at the quoted loop level. The columns "Diagram" specify the topology of the neutrino mass diagram, for example "type-I" for seesaw type-I and so forth. The identification numbers for the 1-loop and 2-loop neutrino mass diagrams are taken from the general topology classification given in [78] and [63], respectively.