Neutrinoless double beta decay in chiral effective field theory: lepton number violation at dimension seven

We analyze neutrinoless double beta decay ($0\nu\beta\beta$) within the framework of the Standard Model Effective Field Theory. Apart from the dimension-five Weinberg operator, the first contributions appear at dimension seven. We classify the operators and evolve them to the electroweak scale, where we match them to effective dimension-six, -seven, and -nine operators. In the next step, after renormalization group evolution to the QCD scale, we construct the chiral Lagrangian arising from these operators. We develop a power-counting scheme and derive the two-nucleon $0\nu\beta\beta$ currents up to leading order in the power counting for each lepton-number-violating operator. We argue that the leading-order contribution to the decay rate depends on a relatively small number of nuclear matrix elements. We test our power counting by comparing nuclear matrix elements obtained by various methods and by different groups. We find that the power counting works well for nuclear matrix elements calculated from a specific method, while, as in the case of light Majorana neutrino exchange, the overall magnitude of the matrix elements can differ by factors of two to three between methods. We calculate the constraints that can be set on dimension-seven lepton-number-violating operators from $0\nu\beta\beta$ experiments and study the interplay between dimension-five and -seven operators, discussing how dimension-seven contributions affect the interpretation of $0\nu\beta\beta$ in terms of the effective Majorana mass $m_{\beta \beta}$.

D Conversion of nuclear matrix elements 48 1

Introduction
The neutrino oscillation experiments of the last two decades have shown that neutrinos are massive particles, requiring an extension of the minimal version of the Standard Model (SM) of particle physics. Neutrinos could have a Dirac mass term, as all other fermions in the SM. This would require including sterile, right-handed neutrinos in the SM Lagrangian, whose only purpose is to generate a neutrino mass. Yet neutrinos are the only observed fundamental and charge-neutral fermions, so they could instead have a Majorana mass. In the SM, a Majorana mass term is forbidden by the neutrino SU (2) L ×U (1) Y quantum numbers, making it impossible to construct a gauge-invariant, renormalizable mass operator in terms of left-handed ν L fields. Thus, in the SM one can distinguish neutrinos from antineutrinos, and define a quantum number, lepton number (L), which is conserved at the classical level. L is, however, an accidental symmetry of the SM. As soon as one introduces non-renormalizable operators, which parameterize physics at energy scales much larger than the electroweak scale, L is broken [1], and neutrinos acquire a Majorana mass, inversely proportional to the scale of new physics Λ. The smallness of the neutrino mass might therefore offer a unique window on high-energy physics. Neutrinoless double beta decay (0νββ) experiments are the most sensitive probe of lepton number violation (LNV). In this process two neutrons in a nucleus turn into two protons, with the emission of two electrons and no neutrinos, violating L by two units. The observation of 0νββ would have far reaching implications: it would demonstrate that neutrinos are Majorana fermions [2], shed light on the mechanism of neutrino mass generation, and give insight on leptogenesis scenarios for the generation of the matter-antimatter asymmetry in the universe [3]. The current experimental limits on the half-lives are already impressive [4][5][6][7][8][9][10][11][12][13], at the level of T 0ν 1/2 > 5.3 × 10 25 y for 76 Ge [12] and T 0ν 1/2 > 1.07 × 10 26 y for 136 Xe [13], with next generation ton-scale experiments aiming at a sensitivity of T 0ν 1/2 ∼ 10 27−28 y. By itself, the observation of 0νββ would not immediately point to the underlying physical origin of LNV. While 0νββ searches are commonly interpreted in terms of the exchange of a light Majorana neutrino, in generic beyond-the-SM (BSM) models, 0νββ receives contributions from several competing mechanisms (for a review see Ref. [14]). Well-studied examples are left-right symmetric models [15][16][17], which contain an extended gauge and Higgs sector, as well as heavy right-handed Majorana neutrinos. In these models light Majorana neutrinos acquire mass via the type-I see-saw (via right-handed neutrinos) and / or the type-II see-saw (Higgs triplet) and can mediate 0νββ. In addition, however, 0νββ receives contributions from the exchange of heavy right-handed neutrinos, mediated by the gauge boson of the additional SU (2) R gauge group, from the mixing of light-and -heavy neutrinos or from the exchange of Higgs triplets [14,[18][19][20]. Depending on the masses of the right-handed neutrinos and gauge boson, and on the Yukawa couplings of the left-and right-handed neutrinos to the Higgs, 0νββ can be dominated by lightneutrino exchange, heavy-neutrino exchange, or receive several contributions of similar size.
Keeping explicit model realizations in mind, in this paper we investigate 0νββ in the framework of the SM Effective Field Theory (SM-EFT) [1,21]. In this framework, the SM is complemented by higher-dimensional operators, expressed in terms of SM fields and invariant under the SM gauge group. The coefficients of these operators are suppressed by powers of the scale Λ at which new physics arises. There is a single gauge-invariant dimension-five operator [1]. This operator violates L by two units, and, as already mentioned, provides the first contribution to the neutrino Majorana mass. Going further, there are no ∆L = 2 dimension-six operators [21,22], but there are several at dimension-seven [23], and -nine [24,25], and higher [26]. 1 Notice that here we are not extending the SM field content with a light right-handed neutrino, but the construction of the effective operators can be generalized to include it [28].
We systematically study the constraints on SU (2) L × U (1) Y -invariant dimension-seven operators from 0νββ. After defining the operator basis in Sec. 2, in Sec. 3 we integrate out heavy SM degrees of freedom, such as the Higgs and the W boson, and match onto a low-energy ∆L = 2 Lagrangian that only contains leptons and light quarks, suitable for the descriptions of low-energy processes such as double-beta decay. The resulting Lagrangian contains the neutrino Majorana mass and transition magnetic moments, dimension-six and -seven semileptonic fourfermion operators, as well as dimension-nine six-fermion operators. Of these operators, those of dimension-six and -seven give rise to non-standard ∆L = 2 single beta decay and to longrange neutrino-exchange contributions to 0νββ not proportional to the neutrino mass. Instead, the dimension-nine operators, which involve four quarks and two electrons, induce new 0νββ contributions without the exchange of a neutrino.
In Sec. 4 we match the quark-level ∆L = 2 Lagrangian onto Chiral Perturbation Theory (χPT), the low-energy EFT of QCD, and we discuss the hadronic input needed to constrain dimension-seven operators. In Sec. 5 we introduce a power counting and derive the neutrino potentials in χPT up to the first non-vanishing orders. The power counting reduces the number of matrix elements that are relevant at leading order in the chiral counting. The contribution of dimension-six ∆L = 2 operators to 0νββ was considered in Refs. [18,[29][30][31][32], while six-fermion dimension-nine were studied in Refs. [24,[30][31][32][33][34][35][36]. In Sec. 5 we discuss similarities and differences between the neutrino potentials we obtain and the existing literature.
In Sec. 6 we obtain our main result which is the derivation of the master formula for 0νββ half-life up to dimension-seven in the SM-EFT expansion and the first non-vanishing order in χPT. For earlier versions of such formula see, for example, Refs. [29,35]. The master formula includes the following important effects: • QCD renormalization group evolution of the dimension-seven operators from the highenergy scale to the weak scale, followed by the QCD evolution of the induced dimension-six, -seven, and -nine operators from the weak scale to the QCD scale.
• Up-to-date hadronic input for the low-energy constants, which are becoming increasingly under control. We find that nine low-energy constants are needed. Six of these are wellknown from either experimental or lattice QCD (LQCD) input, while we estimate the remaining three with naive dimensional analysis. The reader is referred to Table 2 as well as Fig. 5 which illustrates the impact of the uncertainty on the unknown low-energy constants on the constraints on a particular ∆L = 2 Wilson coefficient.
• Consistent power-counting in the chiral effective theory for the neutrino potentials induced by the dimension-seven operators, see Table 4. For some operators we find the first nonzero contributions in 0 + → 0 + transitions to arise at next-to-or next-to-next-to-leading order in the chiral expansion.
• Long-distance contributions arising from either neutrino or pion exchange. When the latter is chirally suppressed, subleading short-range pion-nucleon and contact 4-nucleon contributions are considered. The full interference of all effects is included. 1 All L = 2, B = 0 operators have odd dimension [27].
We find the master formula to depend on only a handful of nuclear matrix elements, a smaller set than typically considered, and we perform comparisons of calculations of the nuclear matrix elements elements already existing in the literature (see Table 5 and Figs. 3 and 4). We test our power counting explicitly by comparing the sizes of different matrix elements and by comparing matrix elements related by symmetry. Bounds on the induced dimension-six, -seven, and -nine operators, as well as the original dimension-seven operators, are obtained in Sect. 7 and presented in Tables 6 and 7 and range from tens to hundreds of TeV, assuming a single dimension-seven operator (Tables 7 and 6) or single induced operator (Table 6) turned on at a time. In Sect. 8 we discuss scenarios in which both a light Majorana neutrino mass and a dimension-seven operator contribute to the 0νββ rate. We study what additional experimental input can be used to disentangle the various ∆L = 2 contributions to 0νββ . We summarize, conclude, and give an outlook in Sect. 9.

Dimension-seven SM-EFT operators
The complete list of dimension-seven ∆L = 2 operators, invariant under the gauge group of the Standard Model, was built in Ref. [23], and it is summarized in Table 1. A subset of the operators was published in Refs. [37,38], and a few redundancies were eliminated in Ref. [39]. At the scale of new physics, Λ, we have the following ∆L = 2 Lagrangian where the first term is the dimension-five Weinberg operator, with C (5) a 3 × 3 matrix in flavor space. Furthermore, i runs over the labels of the operators defined in Table 1. In Table 1, L and Q denote the left-handed quark and lepton doublets, L = (ν L , e L ) T , Q = (u L , d L ) T , while u R and d R are right-handed quarks, singlet under SU (2) L . H denote the scalar doublet where v = 246 GeV is the scalar field vacuum expectation value (vev), h(x) is the Higgs field, and U (x) is a SU (2) matrix that encodes the three Goldstone bosons. The covariant derivative D µ is defined as D µ = ∂ µ − ig s t a G a µ − g τ I 2 W I µ − g Y B µ , where t a and τ I /2 are SU (3) and SU (2) generators, in the representation of the field on which the derivative acts. Y is the hypercharge quantum number, Y = −1/2 for L and Y = 1/2 for H. is a completely antisymmetric tensor, with 12 = +1. C is the charge conjugation matrix, C = iγ 2 γ 0 , which, in this basis, satisfies C = −C T = −C † = −C −1 .
All the couplings C i have lepton flavor indices, which we omit unless explicitly needed, while the couplings of the four-fermion operators in Classes 5 and 6 also carry indices for the quark flavors. Here we are only concerned with couplings to the first generation of quarks.
There are a few special cases in the above operator basis. Firstly, the dimension-five operator and O LH trivially contribute to 0νββ as they simply gives rise to a Majorana mass term below the electroweak scale, The operator O LHB , and the component of O LHW that is antisymmetric with respect to the lepton flavor indices, do not give rise to 0νββ at tree level, but are well constrained by the transition magnetic moments of the neutrinos, as we discuss further in Section 7.1.2. Also, both O (2) LHD and O LLēH do not induce 0νββ at tree level. For these two operators, in Section 7.1.1 we consider radiative corrections, such as the one-loop mixing onto the neutrino mass (O LH ) and magnetic moment (O LHB and O LHW ) operators. The effects of O LLēH are however suppressed by three and one power of the electron Yukawa coupling, respectively. Alternatively, one can study ∆L = 2 decays such as µ + → e +ν eνµ [40]. We briefly discuss bounds on C LLēH arising from muon decay in Sec. 7.1.3. The remaining operators in Table 1 -namely, the following 8 operators O LLQdH , O LLQūH and O LeudH -induce tree-level corrections to 0νββ. Before discussing the effects generated by these operators at the electroweak scale, we briefly comment on the QCD running between the scale Λ and µ ∼ m W . As the majority of the dimensionseven operators do not involve quarks, or only involve a quark vector or axial current, most of these operators do not run under QCD at one loop. The only exceptions are O (1,2) LLQdH and O LLQuH . The latter runs like a scalar current, while the former two operators can be written as combinations of tensor and scalar currents, with O (1) can be obtained by replacing ij mn → im jn . The couplings of these operators are given by, Here the i and j indicate the generation of the left-and right-most lepton fields, respectively. The running is then given by where C F = (N 2 c − 1)/2N c , and N c = 3 is the number of colors. The analytic solutions to these equations are discussed in Appendix B, where we also give numerical relations between C i (Λ) and C i (m W ).
Note that Eq. (5) only takes into account the QCD running, which should be the dominant contribution to the RG up to scales, µ ∼ 10 TeV. For larger renormalization scales, which one is sensitive to if Λ is significantly above the electroweak scale, electroweak contributions could become relevant as well (since α 2 (µ) 1 2 α s (µ) for µ 10 TeV). However, as the largest RG effects result from relatively low scales, µ < TeV, and the electroweak RGEs are currently not known in the literature, we neglect their effects here.

Low-energy Lagrangian
After the breaking of electroweak symmetry, the low-energy ∆L = 2 Lagrangian contains neutrino Majorana masses and transition magnetic moments. In addition, there appear several dimension-six and -seven four-fermion operators as well as dimension-nine six-fermion operators, which give long-and short-distance contributions to 0νββ decay, respectively. We write ∆L=2 + L ∆L=2 .
We choose to work in the mass basis of the charged leptons, but the flavor basis of the neutrinos. This implies that the charged-current interaction and the charged-lepton Yukawa matrix are flavor diagonal, while the neutrino Majorana mass matrix in Eq. (6) is not. Thus the flavor indices i, j in Eq. (6), and in what follows, run over the charged leptons, i, j ∈ {e, µ, τ }. The neutrino mass and magnetic moment terms are discussed in Sec. 7, and here we focus on the operators that mediate 0νββ. Below the electroweak scale the gauge-invariant dimensionseven operators of Table 1 induce the following dimension-six, -seven, and -nine operators +C (6) SR,ijū L d RēL,i Cν T L,j + C The coefficients C (6,7,9) ij are all defined to be dimensionless. Keeping the lepton flavor structure, the matching coefficients for the dimension-six operators at the electroweak scale are given by 2 For the dimension-seven operators we have while the matching conditions for the dimension-nine operators are 5,ij = 0 .
Although we explicitly kept the lepton flavors in the matching coefficients, only one of the elements will actually contribute to 0νββ. This is due to the fact that we require two electrons in the final state, which for the dimension-nine operators implies only the C (9) i, ee element can contribute. In addition, this means that the long-range contributions of the dimension-six and -seven operators have to be mediated by ν e (since the SM weak current has to produce an electron), implying that only the C (6), (7) i, ee component can contribute as well. In the following we therefore drop the flavor indices and use the shorthand, The coefficients in Eqs. (10), (11), and (12) need to be evolved from the matching scale µ ∼ m W to scales µ ∼ 2 GeV, where the matching to chiral perturbation theory and LQCD calculations is performed. The vector operators, C VL, VR and C VL, VR , consisting of quark nonsinglet axial and vector currents, do not run in QCD 3 . The renormalization group equations (RGEs) of the scalar and tensor operators below µ = m W are given by T .
2 Note that the operators in Eqs. (7), (8), and (9) are defined to give rise to d → u transitions, whereas the opposite convention is used for the dimension-seven operators in Table 1. 3 In the MS scheme, the renormalization factor of the non-singlet axial current Z A MS receives non-vanishing contributions starting at two loops [41]. It is however always possible to introduce a finite renormalization that restores the non-renormalization of the flavor non-singlet current [42].

7
Here we have suppressed the flavor indices as the QCD running is independent of them. The above RGEs correct the anomalous dimensions derived in Ref. [43]. The RGEs of the dimensionnine operators are given by [44,45] The analytic (and numerical) relations between C i (m W ) and C i (2 GeV) that result from the above RGEs are discussed in Appendix B.

Chiral Perturbation Theory
Having obtained the relevant ∆L = 2 interactions around 2 GeV, we want to study their manifestation at even lower energies. We do so by applying the framework of chiral perturbation theory (χPT) [46][47][48], and its generalization to multi-nucleon systems, chiral EFT (χEFT) [49][50][51][52]. χPT is the low-energy EFT of QCD and consists of the interactions among the relevant low-energy degrees of freedom (mesons, baryons, photons, and leptons) that incorporate the symmetries of the underlying microscopic theory: QCD supplemented by electroweak fourfermion interactions and, in our case, ∆L = 2 operators. A particularly important role at low energy is played by the approximate symmetry of QCD under the chiral group SU (2) L ×SU (2) R . Since it is not manifest in the spectrum, which instead exhibits an approximate isospin symmetry, chiral symmetry must be spontaneously broken down to the isospin subgroup SU (2) I . The corresponding Goldstone bosons can be identified with the pions. Chiral symmetry and its spontaneous breaking strongly constrain the form of the interactions among nucleons and pions. In particular, in the limit of vanishing quark masses and charges, when chiral symmetry is exact, pion interactions are derivative, allowing for an expansion in p/Λ χ , where p is the typical momentum scale in a process and Λ χ ∼ m N ∼ 1 GeV is the intrinsic mass scale of QCD. These constraints are captured by χPT.
The χPT Lagrangian is obtained by constructing all chiral-invariant interactions between nucleons and pions. In principle, an infinite number of interactions exist, but they can be ordered by a power-counting scheme. We use the chiral index ∆ = d + n/2 − 2, where d counts the number of derivatives and n counts the number of nucleon fields [46]. The higher the chiral index, the more suppressed the effects of a coupling are by factors of p/Λ χ ∼ m π /Λ χ ∼ χ , where we introduced χ = m π /Λ χ . Chiral symmetry is explicitly broken by the quark masses and charges, and, in our case, by electroweak and ∆L = 2 operators, but the explicit breaking is small, and can be systematically included in the power counting by considering m q ∼ m 2 π ∼ p 2 . Because the ∆L = 2 interactions are associated with very small parameters, we only consider operators linear in the ∆L = 2 couplings.
The coupling constants of the effective interactions in χPT, usually called low-energy constants (LECs), are not fixed by symmetry, and they capture the nonperturbative nature of low-energy QCD. In principle these LECs are unknown but their sizes can be estimated from naive dimensional analysis (NDA) [53], or, preferably, fitted to data or calculated from QCD directly for instance by using lattice simulations. As we discuss below, for 0νββ processes most LECs are relatively well known although there are some exceptions.  Figure 1: Examples of irreducible (diagrams (a) and (b)) and reducible (diagram (c)) twonucleon LNV diagrams. Double and single lines denote, respectively, nucleon and lepton fields. The black square denotes an insertion of the neutrino Majorana mass. Notice that diagram (c) is non planar, i.e. the pions "go around" the neutrino line. The first two diagrams respect the χPT power counting, and their scaling is determined by the chiral index ∆ of the vertices and by the number of loops. The sum of two-nucleon irreducible diagrams defines the 0νββ two-nucleon transition operator, or "neutrino potential". In the third diagram the nucleons can be close to their mass shell, and the diagram is enhanced by m N /p with respect to the χPT power counting. This diagram is included by taking the matrix element of the neutrino potential between the nuclear bound-state wavefunctions.
In the mesonic and single-nucleon sector, all momenta and energies are typically ∼ p. The perturbative expansion of the χPT Lagrangian then implies that the scattering amplitudes can also be expanded in p/Λ χ , with every loop (using 4πF π ∼ Λ χ , where F π is the pion decay constant) or insertions of subleading terms in the χPT Lagrangian causing further suppression.
For system with two or more nucleons, in addition to the momentum p, the energy scale p 2 /2m N becomes relevant. Nucleon-nucleon amplitudes therefore do not have an homogeneous scaling in p, and the perturbative expansion of the χPT interactions does not guarantee a perturbative expansion of the amplitudes [49,50]. In Fig. 1 we show two types of contributions to the amplitude. Diagram (c) represents the so-called "reducible" diagrams, in which the intermediate state consists purely of propagating nucleons. In these diagrams the contour of integration for integrals over the 0th components of loop momenta cannot be deformed in way to avoid the poles of the nucleon propagators, thus picking up energies ∼ p 2 /m N from nucleon recoil, no longer a subleading effect, rather than ∼ p. These diagrams are therefore enhanced by factors of m N /p with respect to the χPT power counting and need to be resummed, typically by solving a Schrödinger equation. The resummation leads to the appearance of shallow bound states in systems with two or more nucleons.
Diagrams (a) and (b) exemplify "irreducible" diagrams, whose intermediate states contain interacting nucleons and pions. These diagrams do not suffer from this infrared enhancement, and here nucleon recoil remains a small effect. Irreducible diagrams involving pions and nucleons follow the χPT power counting [49,50] (commonly called "Weinberg power counting"), while the situation is more complicated for contact interactions, where different schemes exist such as "KSW" [54] or pionless EFT [55], where the NN interactions become relatively enhanced.
Reducible diagrams are then obtained by patching together irreducible diagrams with intermediate states consisting of A free-nucleon propagators. This is equivalent to solving the Schrödinger equation with a potential V defined by the sum of irreducible diagrams. Notice, in particular, that the potential is only sensitive to the scale p, and does not depend on properties of the bound states such as the binding energy. For external currents, such as the electromagnetic and weak currents, one can similarly identify irreducible contributions, that can be organized in an expansion in p/Λ χ , and separate them from the effects that arise from the iteration of the strong-interaction potential. For example, diagrams such as Fig. 1(c) are taken into account by taking the matrix element of the neutrino-exchange potential, induced by the irreducible diagrams, between the wavefunctions of the nuclear bound states.
In the following subsections we construct the chiral Lagrangian relevant for 0νββ processes, and discuss the hadronic input needed to determine its couplings. The Lagrangian contains charged-current operators with an electron and an explicit neutrino, which is later exchanged between two nucleons (see Fig. 2(b)) to give rise to long-range neutrino-exchange contributions to 0νββ. For these operators the hadronic input consists of the vector, axial, scalar, pseudoscalar, and tensor nucleon form factors, which, with the exception of a subleading LEC in the tensor form factor, are well determined either experimentally or via LQCD calculations.
In addition, the Lagrangian has operators with pions, nucleons, and two electrons, but no neutrinos (see Fig. 2(c)), which give pion-exchange and short-range contact contributions to 0νββ. In this case new LECs arise from the hadronization of four-quark operators. In the case of purely mesonic operators, these LECs are well determined [56,57]. For pion-nucleon and nucleon-nucleon operators at the moment they can only be estimated with NDA.
In Sec. 5 we then use the Lagrangian constructed in Sec. 4 to derive the two-nucleon operators (the so-called "neutrino potentials") that mediate 0νββ.

The ∆L = 2 chiral Lagrangian
After evolving the ∆L = 2 operators to low energies, µ ∼ 2 GeV, we match them to χPT. The construction of the chiral Lagrangian closely follows that of the standard χPT Lagrangians [47]. We describe the pions by where τ i are the Pauli matrices, F 0 is the pion decay constant in the chiral limit, and we use F π = 92.2 MeV for the physical decay constant. We also introduce the nucleon doublet N = (p n) T in terms of the proton (p) and neutron (n) fields. The pions transform as U → LU R † and u → LuK † = KuR † under SU (2) L × SU (2) R transformations, while the nucleon doublet transforms as N → KN . Additional ingredients are external scalar, vector, and tensor sources in the quark-level Lagrangian, which, for our purposes, take the following form VLē R γ µ Cν T L + C VLē L Ci VRē R γ µ Cν T L + C VRē L Ci where τ ± = (τ 1 ± iτ 2 )/2. The chiral Lagrangian is then given by chiral invariants constructed from the meson and baryon fields and the above spurions, which transform as follows, The dimension-9 operators, C (9) 1 and C (9) 4,5 , can not be written in terms of the above sources and additional chiral constructions are required. The former transforms as 5 L × 1 R while C (9) 4,5 transform as 3 L × 3 R . We will discuss their chiral representations separately below.

Mesonic sector
In the meson sector the interactions that are responsible for long-range neutrino-exchange contributions arise from the standard leading-order (LO) chiral Lagrangian where B is the quark condensate, related to the pion mass by m 2 π = B(m u +m d ). We use (m u +m d )/2 = (3.5 +0.7 −0.3 ) MeV [58], such that B 2.8 GeV. The dimension-six and -seven operators enter through the external sources, l µ , r µ , s, and p. Contributions from the dimension-six tensor operator require two additional derivatives which increase the chiral index by two. As such, the dominant contribution from C (6) T comes from the pion-nucleon sector which is discussed below. One of the advantages of the chiral notation is its compactness, which, however, has the downside of making it more difficult to see to which processes the operators contribute. Here we expand the ∆L = 2 interactions in Eq. (16) up to terms linear in the pion field which provide the main contribution to 0νββ processes In addition, the dimension-nine operators give rise to contributions that do not involve the exchange of a neutrino. In this case, the higher-dimensional operators induce interactions that convert two pions (π − ) into two electrons. Following Refs. [24,56,59] we write the chiral representations of these interactions as where L µ ij = i U ∂ µ U † ij and the dots stand for terms involving more than two pions. By dimensional analysis the low-energy constants g (mix) 8×8 scale as O(Λ 2 χ ), while g 27×1 = O(1). We follow the notation of Ref. [56], in which these three low-energy constants were estimated using SU (3)-χPT relations and LQCD calculations. The values of the LECs we use are given in Table 2, and are in reasonable agreement with naive dimensional analysis.  Table 2: Hadronic input for the LECs g S , g T , g 27×1 , g 8×8 , and g mix 8×8 , at the scale µ = 2 GeV. Currently we lack experimental or LQCD input for the LECs g N π 27×1 , g N N 27×1 , and g T , and we follow naive dimensional analysis.

Nucleon sector
The LO nucleon Lagrangian responsible for long-range neutrino exchange is given by Here v µ and S µ are the nucleon velocity and spin, v µ = (1, 0) and S µ = (0, σ/2) in the nucleon rest frame, andχ + = χ + − Tr(χ + )/2 where χ ± is defined below. We have applied the heavybaryon framework to remove the nucleon mass from the LO Lagrangian [61]. The values of the couplings g A and g T are given in Table 2. The LEC c 5 is related to the strong protonneutron mass splitting and we give its value below. The chiral covariant derivative is defined as The first two terms in Eq. (20) involve contributions from the vector operators C VL(VR) , while the last two terms involve contributions from the scalar couplings C VL(VR) . However, for both the dimension-six vector and tensor operators, the LO terms do not contribute to the 0νββ 0 + → 0 + transitions and non-vanishing interactions only appear at next-to-leading order (NLO).
The relevant NLO corrections can be written as follows where the coefficients of the first two and fourth operators are fixed by reparametrization invariance [62] in terms of the LO nucleon Lagrangian, g M = 1 + κ 1 with κ 1 3.7 the anomalous isovector nucleon magnetic moment, and g T is the only unknown LEC at this chiral order 4 , which by NDA scales as This is the most general chiral-invariant Lagrangian at this order, that is also hermitian, as well as reparametrization, parity, and time-reversal invariant. Apart from long-range neutrino-exchange contributions, the nucleon sector mediates shortrange contributions induced by the dimension-nine operators. These can involve a single pion exchange, through vertices of the formpn π − ee, or through nucleon-nucleon interactions of the formpnpn ee. For the C (9) 4,5 couplings, the short-range contributions to 0νββ are suppressed in the chiral power counting with respect to the long-range pion-exchange terms from Eq. (19). However, for the C (9) 1 coupling, the πN and N N interactions contribute at the same level as the ππ terms of Eq. (19) [24,25]. Thus, for C (9) 1 all three mechanisms have to be considered. Starting with the chiral realization of the pion-nucleon couplings there is one relevant operator, where the dots stand for terms involving additional pions and g πN 27×1 is a LEC of O(1). For later convenience we have pulled out a factor of g A in our definition of g πN 27×1 . For the nucleon-nucleon interactions we also find a single relevant operator where the dots again stand for terms involving additional pions, and g N N

27×1
O(1) is another unknown LEC. As for the previous LEC, we have pulled out a factor of g 2 V in our definition of g N N 27×1 . Additional structures, such aspS µ npS µ n, can be eliminated using Fierz identities and are not independent.
We note here that the distinction between long-and short-distance contributions loses its meaning as one goes to sufficient high order in the construction of the χPT Lagrangian. For example, the operators in Eqs. (24) and (25) receive a contribution from the neutrino Majorana mass, proportional to m ββ /Λ 2 χ , induced by the exchange of hard neutrinos, with momentum |q| > Λ χ , which are integrated out in χPT [65]. Similarly, the operators C

One-body currents for β decays
We now summarize the single β decay amplitude, which provides the building blocks necessary to construct the full 0νββ amplitude. The single β decay amplitude involves two types of diagrams, general relativistic expression for the matrix element p|uσ µν d|n , which depends on four form factors. However, one of these form factors vanishes in the isospin limit and the other involves two derivatives and appears at N 2 LO in the chiral expansion. In the notation of Ref. [63], which is commonly used in the literature [29,30,64], we can identify g T = 2T = −4.54, which, however, does not appear in Ref. [63].
which either involve a single vertex or a single pion exchange between the lepton and nucleon line. Using the Lagrangians constructed in the previous sections, we write the amplitude as with the sources given in Eq. (15). As discussed in Section 4.3, for some operators we will need expressions through NLO in the chiral expansion. Up to NLO, the currents become Here p and p stand for the momentum of the incoming neutron and outgoing proton, respectively, and q µ = (q 0 , q) = p µ − p µ . Furthermore, ε µναβ is the totally antisymmetric tensor, with ε 0123 = +1. At LO in χPT the form factors are given by where we followed the normalization of Ref. [66]. Vector current conservation enforces g V (0) = 1, up to small isospin-breaking corrections. For g A and κ 1 we used the experimental values [58]. There is some disagreement in the literature on the value of g M (0), with some authors using g M (0) = κ 1 = 3.7, rather than the correct g M (0) = 1 + κ 1 = 4.7. The error appears to stem from one of the first papers that studied the contribution of weak magnetism [67], which did not account for the non-anomalous contribution to the isovector nucleon magnetic moment in the non-relativistic limit. We notice that earlier papers, such as [18,68], correctly use g M (0) = 4.7. The isovector scalar charge g S (0) is related to the quark mass contribution to the neutron-proton mass splitting [69]. Using (m n −m p )| str = 2.32 MeV [70] and m d − m u = 2.5 MeV [58] gives g S (0) = 0.93, at the renormalization scale µ = 2 GeV, in very good agreement with the direct LQCD calculation of Ref. [60]. For the isovector tensor charge g T (0) we use the results of Ref. [60,71]. The numerical input we use is listed in Table 2.
The expression of the currents in Eq. (27) in terms of the form factors g V,A,M,S,P,T (q 2 ), while traditional, somewhat blurs the χPT expansion of the various contributions. For instance, at LO in χPT only the pseudoscalar form factor g P (q 2 ) has non-trivial momentum dependence, due to the pion propagator, while all other form factors are purely static. Furthermore, the standard notation in Eq. (27) makes the power counting less apparent by artificially hiding a factor of m N in g P . This means q 2 g P (q 2 )/m N = O(1) is actually a LO contribution, while the magnetic contribution, g M /m N , is suppressed by 1/Λ χ , such that pieces proportional to g M are higher order in the chiral counting. Thus, at LO in χPT, we could drop the magnetic contributions in Eq. (27) and use Eq. (28) for g V,A,P (q 2 ).
The form factors g V,A (q 2 ) and g A (q 2 ) acquire momentum dependence at N 2 LO in χPT. At this order this momentum dependence is encoded in the nucleon isovector charge and axial radii, respectively, r V = 0.76 fm [58] and r A = 0.49 fm [72], corresponding to vector and axial masses Λ V = 0.9 GeV and Λ A = 1.4 GeV in a dipole parameterization of the form factors. This subset of N 2 LO corrections is usually taken into account in the calculation of 0νββ matrix elements by including a dipole form factor for g V and g A , with different vector and axial masses [66]. While including such corrections does not formally improve the accuracy of the calculation, as other N 2 LO contributions, such as pion-neutrino loops or short-range nucleon-nucleon contributions, are not considered, the numerical impact of the axial and vector form factors is not negligible, giving an O(10 − 20%) correction [67,73,74]. This suggests that it might be important to consistently include all N 2 LO corrections to 0νββ .
While the momentum dependence of the g V,A,S,T form factors only enters at N 2 LO in the chiral expansion, the magnetic form factor has a correction at NLO with respect to Eq. (27), due to pion loops 5 [48]. The treatment of the magnetic form factor g M (q 2 ) in the 0νββ decay literature is at odds with this result, as it is often assumed g M (q 2 ) = g M (0)g V (q 2 ), which is not justified in χPT [48].
To conclude this section, we stress that while most of the currents in Eq. (27) have been studied up to N 2 LO, here we do not include these corrections in the construction of the two-nucleon operators that mediate 0νββ, as consistency requires the inclusion of other, unknown, contributions, such as the pion-neutrino loops mentioned above. Thus, even when we use calculations that include partial N 2 LO corrections, our results are formally valid at LO in χPT.

0νββ operators
The ingredients derived in the previous section allow us to construct the two-nucleon operators that mediate 0νββ decays. Fig. 2 shows three possible contributions. The first diagram depicts the standard contribution proportional to the neutrino Majorana mass. The second diagram depicts long-range neutrino-exchange contributions that arise from the ∆L = 2 charged current interactions in Eqs. (7) and (8). These contributions are obtained by combining the one-body currents of the previous section. Finally, operators such as O (1) LHD and O LLduD induce six-fermion dimension-nine operators at the GeV scale, whose contribution to 0νββ decays is represented by the third diagram in Fig. 2. These diagrams do not involve the exchange of a neutrino.
For each operator, we will construct the dominant contribution to 0 + → 0 + transitions, within the framework of chiral EFT. The application of chiral EFT is justified by the separation of the scales involved in 0νββ where the typical momentum exchange between the nucleons is of similar size as the Fermi momentum within nuclei q ∼ k F ∼ m π = O(100 MeV), which is much larger than the reaction Q value, typically around a few MeV.
For the diagrams in Fig. 2(a) and (b), the LO neutrino potential is obtained by tree-level neutrino exchange. This involves the single-nucleon currents, represented by the gray circle and square in Fig. 2, at the lowest order that yields non-vanishing results. Analogously to the stronginteraction potential, the two-body transition operators in chiral EFT are only sensitive to the momentum scale q ∼ k F , and are therefore independent of the properties of the bound states. In particular this implies that the transition operators do not depend on the often used "closure energy"Ē, which encodes the average energy difference between intermediate and initial states. This can be understood from Fig. 2. An insertion of the strong-interaction potential between the emission and absorption of the neutrino in Fig. 2(a) or (b) would generate a diagram which, in the language of Sec. 4, is irreducible. That is, it is always possible to choose the contour of integration such that the energy and momentum of the nucleons in the loop have to be ∼ k F , and the nucleon is far from on-shell. Insertions of the strong interaction potential between the emission and absorption of the neutrino, which would give rise to intermediate nuclear states, are therefore suppressed and can be ignored at LO. Instead, in chiral EFT the dependence on the intermediate states arises from the region where the neutrino momentum is very soft q 0 ∼ |q| k F . The exchange of soft neutrinos gives rise to effects that are suppressed bȳ E/k F [65]. Notice that the situation is different from 2νββ decay, where insertions of the strong interaction potential between the two points where the neutrinos are emitted are not suppressed (in between the first and second neutrino emission, there are only propagating nucleons and the diagrams are "reducible"), and the intermediate states do need to be considered.
For neutrino-exchange contributions, the LO chiral EFT potential is very similar to standard results. In fact, as we will see, the chiral EFT potential reduces to results in the literature in the limit where the closure energy vanishes,Ē → 0. The advantage of chiral EFT is that it is possible to systematically consider subleading corrections. These consist of corrections to single-nucleon currents, which are often included in the literature via momentum-dependent form factors, but also genuine two-body effects, such as loop corrections to Fig. 2(a) and (b), which induce short-range neutrino potentials even for the standard mechanism [65], and three-body effects [75].
Diagram 2(c) does not involve the exchange of a neutrino. In this case the resulting LO potential is of pion range, ∼ 1/m π , or shorter range, ∼ 1/Λ χ . We work at LO in this case as well, but it is straightforward to include subleading corrections in chiral EFT.
In deriving the neutrino potential we take advantage of the fact that the Q value and the electron energies E 1,2 have typical size O(5 MeV) and are thus much smaller than k F . We assign the scaling Q ∼ E 1,2 ∼ m π 2 χ such that these scales can be incorporated in the standard χEFT power counting. The assigned counting generally allows us to neglect the lepton momenta, nuclear recoil, and soft-neutrino exchange, except in a few cases where the matrix element of the LO operator vanishes for 0 + → 0 + transitions. In these cases we consider subleading contributions in the χPT power counting.
Before discussing the contributions in Fig. 2(b) and (c) from the dimension-six, -seven, and -nine operators, we first recall the potential generated by light Majorana-neutrino exchange to establish our notation. For definiteness, we define the neutrino potentials as −A, where A is the amplitude for the process nn → ppe − e − .

Light Majorana-neutrino exchange
In momentum space, the neutrino potential induced by light Majorana-neutrino exchange is where k 1,2 ∼ Q are the electron momenta,q = q/|q|, and the tensor operator is given by where m ν i are the neutrino mass eigenvalues and U is the PMNS matrix. The Fermi (F) function only receives contributions from the vector currents at leading order. In contrast, the Gamow-Teller (GT) and tensor (T) functions receive contributions from the nucleon axial current, including the induced pseudoscalar contribution dominated by the pion pole, and, at higher order, from the nucleon magnetic moment. Here we follow Refs. [67,73,74,76] and separate the direct axial, induced pseudoscalar, and magnetic contributions. We then have the following expressions for For the GT and T functions, we have and In order to compare with the 0νββ literature, we express the long-range neutrino-exchange potentials in terms of g V,A,P,M (q 2 ) where it is implied that they follow the χPT relations in Eq. (28).

Neutrino exchange without mass insertion
The dimension-six scalar operators C (6) SL and C (6) SR , and dimension-seven vector operators, C (7) VL and C (7) VR , give a potential that is very similar to the one that is induced by light Majorananeutrino exchange. At LO in χPT Here we used Eq. (28) to rewrite the potential that is induced by the dimension-seven operators, h GT,7 (q 2 ), as follows which is equal to VR does not contribute at LO because of vector current conservation. The scalar current C SR , combined with the standard model axial current, gives a contribution that is suppressed by q/Λ χ , and, in addition, is parity odd and does not contribute to 0 + → 0 + transitions. The first non-vanishing contributions from the scalar current appear at O( 2 χ ). The pseudoscalar contribution in Eq. (32) has been considered in the literature [29,30,32,64], while the C V L,V R terms have not, even though they appear at the same chiral order. In these works, the neutrino potential is derived by considering the pseudoscalar form factor at q = 0, and by neglecting the induced pseudoscalar component of the axial current. For the pseudoscalar density at zero momentum the value F (3) P = 4.4 is used, which is obtained from a quark-model calculation [63]. These approximations have two consequences. First of all, as pointed out already in Ref. [63], the value F (3) P = 4.4 fails to reproduce the pion pole dominance of the pseudoscalar density, which in χPT gives the much larger F The value of Ref. [63] thus corresponds to using a pion mass of 1100 MeV such that m π ∼ Λ χ . Secondly, neglecting the momentum dependence of the pion propagator in Eqs. (32) and (33) implies that the neutrino potential is of much shorter range than the typical pion range, affecting the value of the nuclear matrix elements.  T induces two operators whose matrix elements vanish in 0 + → 0 + transitions. Including the NLO corrections to the tensor, axial, and vector currents outlined in Section 4.3, we obtain (12) .

O
In addition we find a recoil piece (see Appendix C), which we neglect in our results below. These contributions involve 0νββ operators that depend on the nucleon momenta and whose nuclear matrix elements are unknown. We expect these unknown contributions to be small, however, with respect to Eq. (34) because they are not enhanced by the large isovector nucleon magnetic moment.
Our expressions for the neutrino potentials induced by tensor currents disagree with the literature in two respects. First of all, together with O (6) T , another tensor structure is commonly considered, O [29,30,36,64]. This operator however is identically zero (see Appendix A). This is in disagreement with Refs. [29,30,64] that find a non-zero neutrino potential for this tensor structure. Secondly, the first term in Eq. (34) is sometimes erroneously associated with O

O
The LO operators induced by C (6) VR and C (6) VL also turn out to give vanishing contributions to 0 + → 0 + transitions. By employing the NLO vector and axial currents in Eq. (27) and taking into account the electron momenta and the equations of motion for the electrons, we obtain These expressions agree with Ref. [18,77], in the limit |q| Ē , whereĒ is the closure energy, E = O(10 MeV). In principle there is an additional recoil contribution for the left-handed current C (6) VL , see Appendix C. We neglected this term in the above as it turns out to be suppressed with respect to the magnetic-moment contributions contained in the h M M GT,T terms [77].
For C (6) VR , the first tree-level two-body contribution is proportional to the electron mass or energy and thus of order O( 2 χ ) in the power counting. At the same order one should consider pion-neutrino loops, i.e. the contributions of C (6) VR to short-range ∆L = 2 operators without neutrinos, and three-body operators. While we leave a more detailed study for future work, we stress that the limits we obtain on C VR , and, consequently, on C LeudH , should be taken as order-of-magnitude estimates, rather than rigorous bounds.

Dimension-nine operators
Finally, we discuss the contributions from the dimension-nine operators. In the case of C (9) 4,5 , the most important operators are the pionic ones, while the pion-nucleon and nucleon-nucleon interactions are suppressed by two powers of χ . In contrast, the pionic, pion-nucleon, and nucleon-nucleon couplings all enter at the same order for the operator C  (25), which give rise to the following potential 4 g 8×8 + C The above potential disagrees with parts of the existing literature in several aspects. In Refs. [30,35] the dimension-nine operators defined in Eq. (9) appear as a subset of the most general set of dimension-nine four-quark two-electron operators. The conversion between C (9) 1,4,5 and the coefficients ε defined in Refs. [30,35] is given in App. A. When considering the lowenergy manifestations of these quark-level operators, the authors of Refs. [30,35] only take into account four-nucleon operators, which are of the same form as the one in Eq. (25), and estimate their coefficients by assuming factorization. This approach should provide a reasonable estimate for the bounds on ε LLR 3 as this coupling is related to C 1 , whose neutrino potential receives contributions of similar size from ππ, πN , andN N operators. On the other hand, the contributions of the operators O , are severely underestimated. In these cases, the neutrino potential is dominated by the ππ contribution, given in Eq. (37), and theN N pieces are suppressed by 2 χ . Thus, for O 4,5 , the neutrino potentials of Refs. [30,35] miss the dominant contributions to 0νββ.
The importance of the pion-exchange contributions for certain BSM mechanisms has long been recognized [34,78]. Usually, however, pion exchange is included for the scalar-pseudoscalar operators ε LLR 1 and ε RRR 1 [34,78], while its contribution for vector and axial operators, ε LRR 3 and ε RLR 1 , has been largely ignored [30,32,34,35,78]. These issues were already addressed in Refs. [24,25], which performed a systematic power counting in χPT. The above expression is in agreement with the results of Refs. [24,25]. 20 Finally we comment that in the literature the low-energy constants that describe the hadronization of the four-quark operators have often been estimated using the vacuum insertion approximation [30,34,78]. While, in those cases in which all the relevant hadronic channels are included, this leads to acceptable results, we remark that for the ππ channel more rigorous estimates exist, based on direct LQCD calculations [57] and on SU (3) χPT and LQCD [56]. 6 Master formula for decay rate and nuclear matrix elements Using the potentials in the previous sections we can write down an expression for the inverse half-life for 0 + → 0 + transitions [18,79] T 0ν where E 1,2 are the energies of the electrons, and E f and M i are the energy and mass of the final and initial nuclei in the rest frame of the decaying nucleus. The functions F (Z, E i ) take into account the fact that the emitted electrons feel the Coulomb potential of the daughter nucleus and are therefore not plane waves. They take the following form where R A = 1.2 A 1/3 fm and Z are, respectively, the radius and atomic number of the daughter nucleus. This procedure of calculating the Coulomb corrections assumes a uniform charge distribution in the nucleus and only the lowest-order terms in the expansion in r, the electron position, factor, is taken into account. More precise calculations of the phase space factors apply exact Dirac wave functions [80] and the effect of electron screening [81]. The use of exact wave functions leads to somewhat smaller phase space factors (up to 30% for the heaviest nuclei) while the effects of electron screening are at the percent level [80]. In what follows we do not use Eq. (39) to calculate the phase space factors but instead use the more accurate results of Ref. [32] (see Table 3) which were found to be close to those of Ref. [80]. We only use Eq. (39) when calculating differential decay rates in Sect. 8.1. The Fourier-transformed amplitude is given by 6 where V (q 2 ) is the sum of the potentials discussed in Section 5, and r = r n − r m is the distance between the m th and n th nucleon.
Organizing the amplitude in Eq. (40) according to the different leptonic structures, the contributions of a light Majorana neutrino mass and dimension-seven operators are given by Here we factored out the leptonic structures such that the A i only depend on nuclear (and hadronic) matrix elements and the Wilson coefficients of the ∆L = 2 operators. These are discussed in much more detail below.
With the definitions in Eq. (41), the final form of the inverse half-life can be written as where the G 0i are phase space factors given by Here θ is the angle between the electron momenta and we followed the standard normalization of Ref. [18]. The b 0k factors are obtained from the electron traces that result from taking the square of Eq. (41). They are given by Here we kept terms proportional to k 1 ·k 2 , which are odd in cos θ and therefore do not contribute to the total decay rate, but can potentially be observed in measurements of angular distributions. The definitions in Eq. (44) follow for the most part the existing literature [18]. For G 06 and G 09 , in order not to cloud the chiral scaling of the matrix element, we did not extract a factor of (R A m e ) −1 from A M , as commonly done in the literature [18]. The phase space factors G 06 and G 09 defined in Eqs. (43) and (44) are obtained by multiplying the results in Ref. [18,32] by (m e R A )/2 and (m e R A /2) 2 , respectively. In addition, we removed a factor of 2/9 from the definition of G 04 in order to avoid small dimensionless factors. The phase space factors are summarized in Table 3. These are extracted from the calculation of Ref. [32], with the trivial rescalings discussed above. With the definitions of Eq. (44), the different phase space factors for a given isotope are all of similar size, with no parametric enhancements or suppressions, such that the relative importance of different contributions is determined by the matching coefficients and by the nuclear matrix elements. With the modified phase space factors, we can now apply the χPT power counting purely on the level of nuclear matrix elements.

Nuclear matrix elements
To describe the nuclear parts of this amplitude, we follow standard conventions, e.g. those of Ref. [76], and define the following neutrino potentials 7 where the tensor in position space is defined by S (mn) (r) = 3 σ (m) ·r σ (n) ·r − σ (m) · σ (n) . In the χPT power counting, the matrix elements defined in Eq. . The latter suppression, however, is softened by the large isovector magnetic moment of the nucleon which numerically scales as (1 + κ 1 ) χ O(1).
The A i that appear in Eq. (41) can be obtained from the potentials in Section 5, and, for completeness, we give them explicitly in this section. A ν has the same leptonic structure as the amplitude induced by light Majorana-neutrino exchange. We can divide it in a component which is proportional to the Majorana mass m ββ , a long-distance component M ν, ld arising from the dimension-six and -seven operators in Eqs. (7) and (8), and a short-distance component M ν, sd , proportional to the coefficients of low-energy dimension-nine operators The nuclear matrix element for light Majorana-neutrino exchange has the well-known form T , and the dimension-seven vector operators C VL, VR . The contributions of these operators are not proportional to the neutrino mass, which is replaced by a nuclear scale. We take this into account by factoring one power of the nucleon mass out of the nuclear matrix element in Eq. (47). Combining the results of Secs. 5.2.1 and 5.2.2, we obtain where We see that C SL, SR give the largest contributions to M ν, ld , followed by the tensor operator C (6) T whose effects are formally suppressed by m 2 π /Λ 2 χ , but again this suppression is somewhat mitigated by the large value of g M . The dimension-seven operators are severely suppressed by the Yukawa couplings of the light quarks (since the relative factor can be written as m 2 π /Bv = (m u + m d )/v).
The short-distance component arises from the dimension-nine operators in Eq. (9), which always involve an additional power of 1/v with respect to the contribution from light Majorananeutrino exchange. To compensate for this factor, and for the absence of the neutrino mass, we factored two powers of m N out of the short-distance nuclear matrix element in Eq. (47). We then have where we defined In Eq. (54) we factored the LEC g N N 27×1 out of M sd, 2 as to make the NME independent of the renormalization scale. With the scaling of the LECs discussed in Sec. 4.1, the left-right operators C (9) 4,5 give the largest contribution to M ν, sd , while contributions from the purely lefthanded operator C (9) 1 are suppressed by 2 χ . The dimension-six vector and axial operators C (6) VL, VR induce the additional leptonic structures in Eq. (41). A M is generated through the nucleon magnetic moment and is proportional to C The terms proportional to the electron energies and to the electron mass receive contributions from both C (6) VL and C (6) VR , and are given by where One of the NME combinations is redundant as we can write M me,R = −(M E,L + M E,R + 2M me,L )/2. We choose to eliminate M me,R in the sections below.

Chiral power counting
With these definitions we have introduced nine independent combinations of nuclear matrix elements that determine the 0νββ rate at LO in χPT arising from dimension-5 and -7 operators  Table 4. As discussed in Sec. 5, the smallness of the electron's mass and energy is accounted for in the power counting by assigning the scaling E 1 ∼ E 2 ∼ m e ∼ m π 2 χ = Λ χ 3 χ . The power counting suggests that C give the largest contribution to the inverse half-life, and thus are the most constrained from 0νββ experiments. This expectation is verified in Sect. 7. C (6) T and C (6) VL give contributions of similar size, suppressed by two powers of χ . In both cases, the large nucleon isovector magnetic moment enhances the matrix elements leading to somewhat stronger bounds than expected. C  Table 4: Power-counting estimates of the contribution of low-energy dimension-six, -seven, and -nine operators, as well as m ββ to the amplitudes in Eq. (41). Here ν stands for the contribution of the light Majorana-neutrino exchange mechanism. Furthermore, χ ≡ m π /Λ χ , where Λ χ ∼ m N ∼ 1 GeV is the symmetry-breaking scale. For the power counting, we consider the electron mass and energies and to scale as compared to A M . This expectation is very well confirmed when using realistic values of the nuclear matrix elements. In the case of C VR , there is no contribution to A M , and thus the first correction to the half-life is of O( 3 χ ). As a consequence, the bound on this coefficient, which is particularly interesting for left-right symmetric models, is weaker than for the remaining dimension-six operators as is explicitly found in Sect. 7.
Dimension-seven and -nine operators are further suppressed due to inverse powers of the electroweak scale. Contributions from C  [73,74,76] 8 . In Appendix D we discuss how to convert the nuclear matrix elements of the original references to the notation of Eqs. (45) and (46) (see Table 9). The NME M AA T does not contribute to the light Majorana exchange mechanism, and thus requires a dedicated calculation. This matrix element is important only for C (6) VR and, as we argue in Appendix D, even in this case its contribution is numerically small. Therefore, in Sec. 7 we set M AA T to zero. A few comments are in order. First of all, the neutrino potentials derived in χPT are not sensitive to the closure energyĒ, whereĒ ∼ 1 − 10 MeV is much smaller than the typical Fermi momentum. The relations in Table 9 are valid in the limitĒ → 0, which should be a good approximation if the bulk of the nuclear matrix elements comes from the region r ∼ 1/k F . Secondly, the momentum dependence of the axial and vector form factors is an O( 2 χ ) effect in χPT, and some of the relations in Table 9 neglect the difference between the axial and vector dipole masses, which is justified at leading order. Refs. [76], [32], and [83] computed NMEs that, with these assumptions, should be equal, up to higher-order corrections. By comparing these NMEs we can thus explicitly test the validity of the chiral power counting.     Table 9.
As a first example, if the momentum dependence of g V (q 2 ) and g A (q 2 ) is neglected, the short-distance matrix elements M F,sd and M AA GT,sd are related by a Fierz identity Table 5 shows that the results from Ref. [76] obey Eq. (58) up to corrections that range from ∼ 10% for 76 Ge and 82 Se to ∼ 20% for 136 Xe, while in Refs. [32] and [83] the corrections are roughly 15% and 10% for all the nuclei that were considered. The results of Ref. [85] for 76 Ge also respect Eq. (58) at the 10% level. Once the momentum dependence of g V (q 2 ) and g A (q 2 ) is no longer neglected, the relation in Eq. (58) receives corrections at O( 2 χ ) in χPT, a size consistent with these numerical results.
Furthermore, using the identity q 2 = (q 2 + m 2 π ) − m 2 π , and again neglecting the momentum dependence of g V (q 2 ) and g A (q 2 ), we can derive the following relations between short-and long-distance matrix elements, that are valid through NLO in the chiral counting. The NMEs of Refs. [76], [83] and [85] respect the first three relations to 5% accuracy, the fourth to 10%. For Ref. [32], M AP,P P GT and M AP,P P GT,sd were constructed from pion-range NMEs using the relations of Table 9, which make the first two and the fourth equations in Eq. (59) trivial identities. The third relation in Eq. (59) is non-trivial, and it is well respected by the NMEs in Ref. [32]. These numerical results confirm that the relations in Eq. (59) are accurate up to (5-10)% corrections, which is of the same size as the expected O( 2 χ ) χPT effects. The large number of NMEs computed in Ref. [32] allows for additional consistency checks, which we discuss in Appendix D. In general, for the consistency checks performed in Appendix D, we observe that various relations between NMEs are respected up to 20%-30% corrections, the level one would expect from LO χPT. We conclude that the power counting is working satisfactory although stronger conclusions would require the explicit inclusion of NLO corrections.

Matrix elements from different many-body methods
In Figs. 3 and 4 we show results for the nine combinations of NMEs that determine the contribution of SM-EFT dimension-seven operators to 0νββ, obtained by combining the results of Refs. [76] (blue triangles), [32] (red squares), [83] (green circles), and [84,85] (orange diamonds). The calculation of Ref. [76] is based on the quasiparticle random phase approximation (QRPA) method. Refs. [32] and [83] are shell model calculations. Refs. [84,85] use the interacting boson model. Note that Refs. [32,76,83] include short-range correlations in various ways using CD-Bonn or AV-18 parameterizations. The choice of parameterization has a non-negligible effect for the sd NMEs. In Table 9 we have used results using the CD-Bonn parameterization for [32,76,83].
In order to generate the results presented in Figs. 3 and 4 we made a few assumptions. M T 6 and M sd, 2 depend on the ratios of LEC g T /g T and g ππ, πN 27×1 /g N N 27×1 . In Fig. 3, we assumed the unknown LECs to follow NDA, g T = g πN 27×1 = g N N 27×1 = 1, while g T and g ππ 27×1 are given in Table   28 R(M ν ) Horoi et al. [32] Menéndez et al. [83] Barea et al. [84,85] Figure 3: Comparison of the NMEs obtained using the calculations of Refs. [76] (blue triangles), [32] (red squares), [83] (green circles) and [84,85] [32,76,[83][84][85]. Notation is the same as in Fig. 3. 2. Varying the size of g T has a limited effect on M T 6 , while M sd,2 is quite sensitive to the precise values of the LECs. We discuss this in more detail below. In addition M E,L , M E,R , and M me,L depend on the matrix element M AA T , which is not evaluated in any of the references we use for the NMEs. Fortunately, this matrix element was computed in Ref. [77],  Figs. 3 and 4 show that the nonstandard NMEs computed with different many-body methods differ by at most a factor of 2-to-3. This level of agreement is similar to the one observed for the light-neutrino-exchange mechanism [66] -see the spread in M ν -and leads to an uncertainty in the 0νββ rate of about one order of magnitude. The calculation of Ref. [32] yields values of M P S which have very similar size, but opposite sign with respect to Refs. [76,83,85]. The sign difference has no impact in the single-coupling scenario explored in Sec. 7. It will affect scenarios in which several operators are turned on at the same time, but in this case the effect is mitigated by the ignorance of the relative phase between the coefficients. A similar argument applies to M me,L and M me,R , which, using the results of Ref. [32] are found to have similar size, but different sign with respect to the other calculations. The uncertainty on the shortdistance NME M sd,2 is somewhat larger than for the other NMEs. This is not unexpected as such matrix elements depend on short-distance details of nuclear wave functions which are more model dependent then long-range aspects. The relative sizes of the NME combination M sd,2 for various isotopes vary strongly between Refs. [32,76,83]. Although we do not understand this behaviour in detail, it might be related to possible accidental cancellations between the various contributions to M sd,2 . In the next section we explore the consequences of these uncertainties on the constraints on the scale of BSM lepton-number-violating physics.
It  Table 3). The above considerations imply that seven combinations of NMEs dominate 0νββ in the SM-EFT.

Single-coupling constraints
In this section we discuss the constraints on the low-energy operators, as well as the fundamental dimension-seven operators that arise at the scale Λ. We start by considering the bounds from 0νββ experiments and discuss other relevant observables in Sect. 7.1. Throughout this section we will assume that only one operator is present at a time. We study scenarios involving multiple 30 76 [9,11], CUORE [7], and KamLAND-Zen [13] experiments, assuming C i (µ = 2 GeV) = v 3 /Λ 3 i . The left, middle, and right tables correspond to the matrix elements of Refs. [76], [32], and [83], respectively. The lower limits on Λ are shown in units of TeV. couplings in Sect. 8. We apply the following experimental limits [7,12,13,86] (all at 90% c.l.) By inserting the phase-space factors of Table 3 and the NMEs in Table 5 into Eq. (42), we obtain limits on the coefficients of the ∆L = 2 operators. In Table 6 we show bounds on m ββ and the low-energy dimension-six, -seven, and -nine operators of Eq. (9), which were derived using the NMEs of Refs. [76], [32], and [83] in the left, middle, and right panels, respectively.
Using NMEs from Ref. [76] we find an upper bound m ββ < 0.084 eV, and slightly weaker bounds for the other NMEs. The limits we obtain are in agreement with, for example, Ref. [32]. All bounds are somewhat weaker than the most stringent bound reported in Ref. [13], m ββ < 0.061 eV which is based on different NMEs than considered here.
For the non-standard operators, Table 6 shows the constraints on the scale of new physics, Λ, assuming that C i (µ = 2 GeV) = v 3 /Λ 3 and only one coupling is turned on at a time. In addition, we assumed natural values for the unknown LECs, g T = g πN 27×1 = g N N 27×1 = 1. As expected from the discussion of the previous section, the most stringent constraints arise in the case of C  Table 4 would predict the limit on C (6) T,VL to be weaker by 2/3 χ , the actual constraints are somewhat stronger than expected due to the large isovector magnetic moment. For most of the remaining couplings the limits closely follow what one would expect from the power counting. For example, the limits on C VL,VR and C (6) VR are weaker than the limits on C  Table 4. Finally, we would expect the limits on C (6) VR to be weaker than the limit on C  agrees fairly well with the actual results. Here we took into account by hand the large nucleon magnetic moment.
The case of C 1 requires additional explanation. From the power counting we would expect this coupling to contribute at the same order as C (7) VL,VR . However, the matrix element M sd, 2 receives several contributions proportional to unknown LECs, g πN 27×1 and g N N 27×1 . As a result, the contribution of C can vary substantially depending on the values and signs of these LECs. This is illustrated in Fig. 5 where we show the constraint on C disappears. Although such a near-exact cancellation is not expected, and is sensitive to higherorder corrections, the limits on the scale Λ for C (9) 1 appearing in Table 6 should be taken as an order-of-magnitude estimate, at least until the values of g πN,N N 27×1 are further constrained. In contrast, varying the sign of the only other unknown LEC, g T , only leads to O(10%) effects in the limits on Λ for C (6) T . Although the above constraints are useful to test the power counting, the fundamental ∆L = 2 operators of interest are the dimension-seven operators of Table 1. We present the limits on these couplings in Table 7, where the left, middle, and right panels again employ the NMEs 32 76 [9,11], CUORE [7], and KamLAND-Zen [13] experiments, assuming C i (µ = Λ) = 1/Λ 3 . The left, middle, and right tables correspond to the matrix elements of Refs. [76], [32], and [83], respectively. The limits on Λ are shown in units of TeV.
of [76], [32], and [83], respectively. The bounds on the scale of new physics are obtained by assuming a single coupling is present at the high scale, and C i (µ = Λ) = 1/Λ 3 . The strongest limits are derived in the case of C LLQdH and C LLQuH because these operators mainly induce the stringently constrained C (6) SL,SR . Instead, the weakest limits are obtained in cases where only the low-energy dimension-seven and -nine operators are induced. This is the case, for example, for C (1) LHD and C LHW , which both mainly contribute to C VL and C (9) 1 . Since these operators induce C (9) 1 , the corresponding limits are sensitive to the values of the unknown LECs, g πN,N N 27×1 . In Fig. 6 we present the same information in a different format, focusing on the bounds on the dimension-7 operators arising from the KamLAND-Zen experiment [13].
It should be noted that the Wilson coefficients will in general depend on a dimensionless coupling, c i , in addition the scale Λ, i.e. C i = c i /Λ 3 . The presence of these c i implies that the limits on Λ in Table 7 (where we assumed c i = 1) do not necessarily correspond to constraints on particle masses in any given BSM theory. In particular, in weakly coupled BSM theories, c i < 1, the limits on the masses of particles could be significantly weaker than those on Λ given in Table 7. Thus, the stringent bounds on Λ derived above do not necessarily imply that the responsible BSM physics is out of reach of collider searches. Apart from a simple rescaling of the limits in Fig. 6, dimensionless couplings, c i = 1, would change the starting point of the RGEs. However, the numerical impact of such a change in Λ is rather minimal. For example, changing the starting point of the RG from Λ = 50 TeV to Λ = 100 TeV, changes the running of the C i by no more than 10%.
An alternative way to present the limits is shown in Table 8, where we show the bounds on the dimensionless couplings, c i = Λ 3 C i (Λ). Here we picked the scale Λ to be 10 TeV, and derived constraints using several calculations for the NMEs [32,76,[83][84][85]. The bounds in Table 8 are inversely proportional to these NMEs, c i ∝ M −1 i , while the limits on the scales have a much weaker dependence, Λ ∝ M 1/3 i . As a result, the variation between different nuclear calculations is more pronounced in Table 8 than in Table 7.   [87] and KamLAND-Zen [13] experiments. Here we assume c i (µ = Λ) = C i (Λ) Λ 3 and choose the scale of BSM physics to be Λ = 10 TeV. The columns from left to right, correspond to the matrix elements of Refs. [76], [32], [83], and, in the case of 76 Ge, [84,85], respectively. Hyvärinen et al. [76] Horoi et al. [32] Menéndez et al. [83] Figure 6: Constraints from the KamLAND-Zen experiment [13] on the scale of the dimensionseven operators. We assume C i (µ = Λ) = 1/Λ 3 and only turn on one operator at a time.

Other constraints
Although 0νββ leads to stringent constraints on the C i couplings, reaching scales of O(100 TeV), it is interesting to see how these compare to constraints from other probes. In particular, all operators in Table 1 induce radiative corrections to the neutrino masses. In Sec. 7.1.1 we therefore discuss the naturalness bounds that can extracted from the neutrino masses. We find that for several operators they are stronger than the bounds from 0νββ. Considering additional probes is particularly important for the operators C LHB and C LLēH , which do not induce 0νββ at tree level, and C (2) LHD , whose contribution to 0νββ is suppressed by the electron energy, and was not considered in Secs. 5 and 6. We address the contributions of these operators to the neutrino masses in Sec. 7.1.1, and take into account bounds from the neutrino transition magnetic moments in Sec. 7.1.2, and from non-standard muon decays in Sec. 7.1.3.

Neutrino mass
The operators in Table 1 can generate neutrino masses. The tree-level contribution is The other C i do no contribute at tree level, but can contribute to C LH through RG effects between µ = Λ and µ = m W . The complete neutrino mass is a combination of the contributions of the dimension-seven operators and the Weinberg operator. In total we have m ν = m ν is the contribution from the Weinberg operator. Since m (0) ν is unknown we can only set constraints if we assume that the dimension-five and -seven contributions are not unnaturally large compared to the total neutrino mass. That is, we assume there is no large cancellation between m (0) ν and δm ν . To get an idea of these naturalness limits we will, somewhat arbitrarily, impose |δm ν | 1 eV.
From Eq. (61), we can already estimate the constraint on C LH . Assuming C LH (µ = Λ) = 1/Λ 3 , we get Λ > 1200 TeV. For the other dimension-seven operators that contribute at loop level, we require the evolution between µ = Λ and µ = m W . The relevant one-loop RGE is given by The above expression provides us with C LH (µ = m W ), which together with Eq. (61) and |δm ν | 1 eV, leads to the constraints where we again assumed C i = 1/Λ 3 . Contributions of the operators appearing in the second line of Eq. (62) are severely suppressed by three powers of small Yukawa couplings. The corresponding limits are well below the electroweak scale such that we do not obtain sensible constraints.
Here we only considered contributions to the neutrino masses through corrections to the dimension-seven coupling C LH . In principle, one could consider corrections directly to the dimension-five coupling, C (5) , in Eq. (1) as well. Below the scale Λ, the SU (2)-invariant dimension-seven operators do not mix with this dimension-five operator. However, assuming the dimension-five term is not protected by symmetry considerations, one might expect the BSM interactions that induce the C i appearing in Eq. (62), to contribute to C (5) as well. These contributions would result from matching the BSM theory to the EFT and, if they arise from loop diagrams, could in principle scale as C (5) ∼ 1 (4π) 2 1 Λ , in which case they would dominate over those in Eq. (62) by a factor of Λ 2 /v 2 . Such contributions would lead to more stringent limits than those in Eq. (63). On the other hand, it is possible to realize smaller contributions to the neutrino masses than those induced by Eq. (62) if there is a fine-tuned cancellation at work. Which of these scenarios is realized, as well as the mentioned matching contributions, depend strongly on the specific BSM theory above the scale Λ. Here we refrain from estimating such model-dependent effects and only consider the terms that are calculable within the EFT framework. Nevertheless, one should keep in mind that specific BSM theories could give larger contributions to the neutrino masses than those captured by Eq. (62).
It is certainly possible to avoid the above naturalness limits by allowing for some amount of fine-tuning between, for example, dimension-five and -seven contributions to the neutrino mass. Nevertheless, taken at face value, the contributions to δm ν can lead to very stringent constraints. This is certainly true for C LH and C (2) LHD , for which the limits reach O(100 TeV) or more, while these couplings would be left unconstrained by 0νββ. Note that these naturalness limits even exceed the 0νββ constraints for C LHW and C (1) LHD , while 0νββ is more constraining for C LHDe (as well as for C LLQuH and C (1) LLQdH ). Of the remaining operators, C LHB does not contribute at one loop as it is anti-symmetric in flavor space, while C (2) LLQdH , C LLduD , and C LeudH , mix with C LH at two loops and require, respectively, one, two, and three Yukawa insertions. The 0νββ limits are more stringent in these cases, and we do not consider the contributions to δm ν .

Magnetic moments
Apart from neutrino masses, the operators in Table 1 also induce contributions to the magnetic moment of the neutrinos. These magnetic moments can be constrained by neutrino-electron scattering in solar and reactor experiments [38,88,89], or through astrophysical limits from globular clusters [90]. As we are mainly interested in an order-of-magnitude estimate, here we will employ the limits of Ref. [89] from the scattering of solar neutrinos.
Tree-level contributions of the dimension-seven operators to the magnetic moments are where µ and C LHB are anti-symmetric in flavor space. Following the notation of Ref. [89], the transition magnetic moments can be parametrized by three complex parameters, Λ i , as follows, where the PMNS matrix, U , appears due to the rotation to the mass basis. The constraints derived in Ref. [89] are In principle, a detailed analysis should take into account the flavor structure of C LHB,LHW as well as the unknown phases in U . As we are mainly interested the order-of-magnitude of the limits, we take the following estimate For C LHW this limit is weaker than both the limit from 0νββ as well as the naturalness constraint from the neutrino mass. However, the neutrino magnetic moments do provide the most stringent limit on C LHB , whose contributions to 0νββ and the neutrino mass are suppressed.

Muon decay
The operator O LLēH does not contribute to 0νββ at tree level, and its contribution to the neutrino mass in Eq. (62) is suppressed by three powers of the electron Yukawa coupling, leaving the coefficient C LLēH poorly constrained. In this section we discuss the constraints on C LLēH from non-standard muon decays. After electroweak symmetry breaking, the ∆L = 2 Lagrangian relevant for muon decay is where the coefficients C S and C T are C µe S,T , and its hermitian C µe * S,T , mediate, respectively, the ∆L = 2 decays µ + → e +ν eνµ and µ − → e − ν e ν µ , while C eµ S,T and C eµ * S,T induce µ − → e −ν eνµ and µ + → e + ν e ν µ . The experimental analysis of Ref. [40] searched forν e in the decay products of a µ + at rest, by looking for the charged current processes pν e → e + n and 12 Cν e → e + n 11 B following the decay of the muon. The muonic neutrino is not identified, and thus the experiment constrains µ + → e +ν e (ν + ν). The experimental setup is such that the contribution of neutrino oscillations, ν µ →ν e , is negligible [40]. If, in addition, we assume that there are no ∆L = 0 lepton-flavor violating operators, which would for example induce µ + → e +ν e ν µ , the limits on the branching ratio can be used to put bounds on C µ e S, T . In terms of C µe S, T , the branching ratio is The dependence of the decay rate on theν e energy is determined by the Michel parameterρ, which, at tree level, isρ = 3/4 for the scalar, andρ = 1/4 for the tensor operator.

Two-coupling analysis
The single-coupling limits of section 7 clearly show the constraining power of the 0νββ experiments, as they reach scales of O(100 TeV). However, in realistic lepton-number-violating scenarios one would generally expect to generate multiple ∆L = 2 couplings at the scale of new physics. In this section, we discuss scenarios in which both m ββ and a dimension-seven operator are turned on simultaneously. We study how such scenarios differ from the well-known light-Majorana neutrino case. Finally, in section 8.1, we briefly consider the possibility of distinguishing different ∆L = 2 operators using the energy and/or angular distributions of the electrons emitted in 0νββ. We begin with showing the limits in the |m ββ | − Λ 3 C LLQuH plane in Fig. 7. Here we assumed Λ = 600 TeV and used the NMEs of Ref. [76]. In the left panel we take a specific value for the relative phase between the dimension-seven coupling and m ββ , namely, Arg (C LLQuH m * ββ ) = 3/4π. As one can see, in this case the experimental limits form ellipses in the m ββ − Λ 3 C LLQuH plane. For a generic relative phase the picture is qualitatively the same. However, specific values of the relative phase, namely, 0 and π, allow for cancellations between the dimension-seven and m ββ contributions. As a result, free directions appear once we marginalize over the relative phase. This is clearly shown in the right panel of Fig. 7.
These free directions appear in part because C LLQuH contributes to the same leptonic structure as m ββ (see e.g. Eq. (47)). As such, we also consider operators that generate different leptonic structures. We show the m ββ -Λ 3 C LeudH plane in Fig. 8, now assuming Λ = 40 TeV. Although we marginalized over the relative phase, no free directions appear because the different leptonic structure prohibit a (complete) cancellation between m ββ and the dimension-seven contribution. Finally, both Fig. 7 and 8 illustrate that the different nuclei considered here do not have very different sensitivities, i.e. the ellipses and bands all have roughly the same slope. This is a generic feature that does not depend on the dimension-seven coupling under consideration. Unfortunately this implies that it will be difficult to unravel the underlying ∆L = 2 mechanism from just nonzero 0νββ total decay rates.
It is interesting to consider the impact of the dimension-seven operators on the interpretation ββ as a function of the lightest neutrino mass in Fig. 9. The allowed areas are obtained by using the standard parametrization in terms of the neutrino masses, m ν i , the sines (cosines) of the neutrino mixing angles, s ij (c ij ), and the Dirac phase δ 13 , We then marginalize over the Majorana phases, λ 1,2 , and the experimentally allowed values of the Dirac phase, while setting the mixing angles to their central values [58]. The top-left (-right) panel of Fig. 9 depicts the normal (inverted) hierarchy for several values of C LLQuH . Blue, gray, and green bands assume C LLQuH = {−1, 0, 1} · Λ −3 , respectively, with Λ = 600 TeV. The current limit on m    Figure 9: The left (right) column shows the allowed values for the effective parameter |m eff ββ | (defined in Eq. (72)) as a function of m lightest ν for the normal (inverted) hierarchy. The gray bands depict the case with all dimension-seven operators set to zero, while the red horizontal line shows the 0νββ limit from 136 Xe. In the top panels, the green and blue bands show the allowed values for the case that C LLQuH = 1/Λ 3 and C LLQuH = −1/Λ 3 , respectively, assuming Λ = 600 TeV. The middle panels show the same scenarios after marginalizing over the possible phase of C LLQuH . I.e. we take C LLQuH = e iα /Λ 3 and marginalize over α. Finally, the bottom panels show C LeudH = e iα /Λ 3 marginalized over α, and assuming Λ = 40 TeV. ββ to be nonzero and a finite 0νββ must exist at some level. We show similar plots in the middle row of Fig. 9, where the green band is obtained from marginalizing over the phase of C LLQuH . For a wide range of m lightest ν , the effective parameter m (eff) ββ and thus the 0νββ rate, can go to zero even for an inverted hierarchy.

Inverted Hierarchy
C LLQuH generates the same leptonic structures as m ββ and it is interesting to look at a coupling that induces a different phase-space factor. In the bottom row of Fig. 9, we depict the allowed region for m (eff) ββ assuming that m ββ and C LeudH are both turned on. In this case the effective parameter m (eff) ββ is always nonzero and the allowed m (eff) ββ region simply shifts upwards for the normal and inverted hierarchies (left and right panels, respectively).

Pinpointing the ∆L = 2 mechanism
In the best-case scenario in which a 0νββ signal is measured, it would be crucial to identify the underlying ∆L = 2 mechanism. Of course, a nonzero value of T 0ν 1/2 could be generated by any of the dimension-five or -seven couplings and additional information is required to disentangle them. In principle, one could think of using measurements of T 0ν 1/2 in different nuclei. Although the NMEs generally show similar patterns for different nuclei, leading to degenerate sensitivities, this is not always the case for the phase space factors. In particular, G 02 has an increased sensitivivity to the Q value compared to the other phase space factors (see Eq. (44)). This means that, 128 Te, which has a rather small Q value, will have a significantly smaller value of G 02 than 76 Ge. As C (6) VR contributes proportional to G 02 , this in turn implies that 128 Te is less sensitive to C (6) VR compared to 76 Ge [91,92]. This in principle provides a way to disentangle C (6) VR from the other operators, by measuring the decay rates in several isotopes. However, as discussed above, the nuclei considered here ( 76 Ge, 82 Se, 130 Te, and 136 Xe) have very similar sensitivities to the dimension-seven couplings, something which is even worsened once nuclear and hadronic uncertainties are taken into account. It would therefore be difficult to pinpoint the underlying ∆L = 2 mechanism from just 0νββ total rates of the nuclei under consideration here. Similar conclusions were reached in Refs. [93,94].
Additional information could come from ∆L = 2 signals at colliders such as the LHC. There are certainly scenarios in which colliders can compete with the 0νββ measurements [95]. While the limits derived in Sect. 7 already put some of the operators at very high scales of O(100 TeV), two effects, in combination, may mitigate these bounds and make collider searches competitive with 0νββ experiments. First, in specific models the Wilson coefficients C i may naturally be suppressed by small Yukawa couplings, allowing for a smaller scale Λ consistent with the 0νββ bounds obtained here. 9 In addition, for a fixed mass scale, we saw in Sect. 6.3 that the uncertainty in the values of the nuclear matrix elements can lead to an order-of-magnitude variation in the predicted 0νββ rate. This is the appropriate measure for comparison, since in the contact limit , the production rate at a collider experiment has the same scaling with Λ as the 0νββ rate, yet is unaffected by uncertainties in the nuclear matrix elements. The rate at a collider may be even higher if intermediate particles can be produced on-shell. It therefore remains an open question whether direct searches at the LHC or a future collider would be able to see a signal from the fundamental ∆L = 2 operators.
As such, here we focus on additional observables that can be measured by the 0νββ experiments [96], namely, the angular and energy distributions of the electrons produced in 0νββ. These distributions are determined by the leptonic structures in Eq. (41). Several dimensionseven operators generate different leptonic structures such that the angular and energy distributions carry information about the C i . Unfortunately, only the low-energy couplings C (6) VL and C (6) VR induce leptonic structures different from the one generated by m ββ . These vector couplings are induced by the high-energy dimension-seven couplings C LHDe and C LeudH . Consequently, all other dimension-seven couplings induce the same lepton structure as m ββ and will be degenerate with m ββ and each other.
Thus, the angular and energy distributions can in principle be used to disentangle C LHDe and C LeudH from the remaining couplings. These two couplings induce a dependence on cos θ whose slope has the opposite sign of the one induced by m ββ . In addition, although C LHDe gives rise to an energy dependence that is very similar to m ββ , the energy distribution of C LeudH is significantly different. This is illustrated in Fig. 10 which shows the angular and energy dependence in the left and right panels, respectively. The different lines correspond to the case of nonzero m ββ (dashed black), nonzero C LeudH (dashed red), and a scenario where both couplings are turned on (orange band). In the latter scenario we set |m ββ | = 0.05 eV and C LeudH = e iα /Λ 3 with Λ = 40 TeV, while we varied over the relative phase α. As can be seen from the left panel, the slope of the cos θ dependence does indeed differ by a sign between m ββ and C LeudH . Once both couplings are turned on the resulting slope lies somewhere in between the two extremes. Although many couplings could induce the same cos θ dependence as m ββ , the opposite slope can only point to either C LeudH or C LHDe .
The energy dependence is shown in the right panel of Fig. 10. Again there is a clear difference between the case in which only m ββ is turned on (dashed black) or only C LeudH is nonzero (dashed red). As one would expect, including both couplings (orange band) gives a combination of the two dashed lines. It should be noted that only C LeudH is able to induce an energy dependence that significantly differs from the m ββ case, while the C LHDe case looks very similar to that of m ββ .

Summary, conclusions, and outlook
In this work we have investigated neutrinoless double beta decay in the framework of the Standard Model effective field theory. In principle, the dominant contribution to 0νββ arises from the dimension-five Weinberg operator which is only suppressed by one power of the scale of beyond-the-SM physics. However, in several models competing contributions arise from higherdimensional operators and we therefore extended the analysis to include all ∆L = 2 operators of dimension seven.
In the first part of this work we classified the different dimension-seven operators and studied how they manifest at a relatively low-energy scale of a few GeV. We studied the evolution of the operators to lower energies by considering renormalization-group running and threshold effects from integrating out relatively heavy SM fields such as the Higgs and electroweak gauge bosons. This analysis gives rise to a set of effective dimension-six, -seven, and -nine ∆L = 2 operators that we evolve to slightly above the QCD scale using their renormalization group equations. All operators scale as 1/Λ 3 , where Λ is the scale of BSM physics, and their effective dimension is determined by powers of the electroweak scale.
In the second part we applied the framework of chiral effective field theory to construct the effective ∆L = 2 hadronic Lagrangian. For each effective operator at the quark-gluon level we build the chiral Lagrangian up to the order where we find the first non-vanishing contribution to the 0νββ decay rate. Depending on the effective operator under consideration, the chiral Lagrangian consists of pionic, pion-nucleon, and/or nucleon-nucleon interactions. Armed with the chiral Lagrangian we calculated effective two-nucleon 0νββ operators in a consistent powercounting scheme, and derived, within the same scheme, a Master formula for the 0νββ decay rate. Our results contain several new aspects • We used up-to-date hadronic input for several low-energy constants that connect ∆L = 2 quark-gluon operators to ∆L = 2 chiral operators. While remarkable progress has been made in recent years on several of the LECs, others, in particular those associated to ∆L = 2 pion-nucleon and nucleon-nucleon interactions, are still unknown. In the future it will be important to further constrain or compute these LECs. For illustrative purposes, we show in Fig. 5 how the current bound on the Wilson coefficient C 1 is affected by the uncertainty on the unknown LECs.
• We introduced a power-counting-scheme for 0νββ operators which includes, apart from the standard χEFT counting rules, the additional scales associated with 0νββ : the so-called "closure energy" and the Q value of the reaction. We showed that up to leading order in the power counting, the rate does not depend on the closure energy. In addition, we find that the leading-order rate only depends on several nuclear moments (scalar, vector, axial, and tensor) and not on the associated radii which are often included. These considerations greatly reduce the number of nuclear matrix elements that needs to be calculated. We confirmed these power-counting predictions by explicit comparison with several sets of nuclear matrix elements calculated in the literature.
• Based on the extended χEFT power counting we identified nine combinations of nuclear matrix elements, which determine the leading-order 0νββ rate up-to-and-including dimension-seven operators in the SM-EFT. Two combinations of nuclear matrix elements turned out to be numerically suppressed due to factors beyond the power-counting scheme (the large size of the nucleon isovector magnetic moment and the smallness of the electron mass with respect to the reaction Q values.) As such, the 0νββ rate is dominated by a relatively small set of nuclear matrix elements.
• We find that the nuclear matrix elements that are needed to constrain the contributions of dimension-seven operators can be lifted from existing calculations of 0νββ. With the exception of M AA T , the required matrix elements can be deduced from calculations of lightand heavy-Majorana-neutrino exchange, provided that the various components, M ij GT,T (sd) in Eq. (48), are listed separately and the calculations include the contributions from weak magnetism and induced pseudoscalar form factor.
• The matrix element M AA T is important in constraining C (6) V R , but is not evaluated in any of the recent nuclear matrix element literature. Here we used the value computed in Ref. [77]. It would be preferable if in the future this matrix element is reported along with the other M ij F,GT,M nuclear matrix elements such that all nuclear physics input to the 0νββ rate is internally consistent.
• We have compared different sets of nuclear matrix elements obtained with various manybody methods. We find that uncertainties on the non-standard matrix elements, based on the spread of the results, are of similar size as the uncertainty on the light-Majorananeutrino-exchange matrix elements. Typically the matrix elements vary at most by factors of two-to-three (and several are in much better agreement) depending on the chosen nuclear method. However, the sign and relative sizes of the matrix elements are in good agreement with each other and the chiral power counting.
In the final phenomenological part of this work, we studied the constraints on the fundamental ∆L = 2 operators. The above-described framework provides essentially a dictionary between high-scale ∆L = 2 physics and low-energy 0νββ measurements such that constraints on the scale of BSM physics can be immediately obtained. We obtain several interesting conclusions: • Depending on the ∆L = 2 operator under consideration, the limits on the BSM scale varies from Λ > 10 TeV to Λ > 400 TeV. For most operators these limits on the scale Λ are not too much affected by hadronic and nuclear uncertainties, except for operators which mainly induce so-called short-distance contributions to 0νββ which depend on unknown LECs associated to ∆L = 2 pion-nucleon and nucleon-nucleon interactions. LQCD calculations of these LECs, along the lines of Refs. [57,97], could improve this situation. Several dimension-seven SM-EFT operators do not contribute to 0νββ at a significant level. We studied complementary observables, such as the neutrino mass and magnetic moment, and muon decay, that can be used to probe such couplings.
• We find that 0νββ experiments with different isotopes (we studied 76 Ge, 82 Se, 130 Te, and 136 Xe) are rather degenerate with respect to the different ∆L = 2 mechanism they are sensitive to. We have illustrated this in Figs. 7 and 8 where it can be seen that different isotopes probe roughly the same combination of ∆L = 2 operators.
• The inclusion of non-zero dimension-seven ∆L = 2 couplings can affect the standard interpretation of (the absence of) 0νββ signals in terms of light Majorana-neutrino exchange. In this framework, it is possible to rule out the inverted ordering of the neutrino mass spectrum with sufficiently sensitive 0νββ experiments. The upper panels of Fig. 9 illustrates that this is no longer necessarily true once dimension-seven operators are included in the analysis, although some fine-tuning is required to suppress the 0νββ rate. At the same time, the inclusion of dimension-seven operators can lead to a non-zero 0νββ rate for all values of the lightest neutrino mass even for a normal hierarchy.
• While total 0νββ rates of different isotopes have little discriminating power with respect to the underlying source, additional information could be obtained by angular and energy differential rates. As shown in Fig. 10, the differential rates can potentially separate several ∆L = 2 dimension-seven operators from the dimension-five and other dimensionseven operators. This is particularly relevant for BSM models, such as left-right symmetric models, that induce low-energy vector-like ∆L = 2 operators.
Our work can be extended in several ways. First of all, in several models also ∆L = 2 dimension-nine operators provide relevant 0νββ contributions. We aim to extend the framework to include these operators in future work. This will enable one to match specific UV-complete models to the effective field theory framework. In particular, this would allow for a global analysis of Standard Model extensions involving lepton-number violation, including 0νββ and high-energy probes at the LHC or future high-energy colliders.

A Comparison with other operator bases
In this Appendix we compare our operator basis and Wilson coefficients to the one previously used in the literature. The basis introduced in Refs. [29,30,35] contains at the hadronic scale operators of dimension six (long range part), related to the ones in (7) and dimension nine (short range part), related to the ones in (9). They do not consider operators of dimension seven (see (8)) which naturally arise in our analysis based on SU (2) × U (1) gauge invariance.
For the short-range effective couplings associated to dimension-nine six-fermion operators, [30,35] the effective couplings xyz i (with x, y, z labeling the chirality of the two hadronic densities and the leptonic current, in that order) the mapping goes as follows: 5 .

B RG evolution
In this appendix we briefly discuss the scale dependence of the couplings mentioned in sections 2 and 3. The running of the dimension-seven operators between the high scale, Λ, and the electroweak scale is given by C(µ) = U (µ, Λ) · C(Λ), C = (C LLQuH , C LLQdH , C T (m W ) .
The couplings C SL(SR) decrease in the UV while the tensor coupling C T increases. The dimension-seven operators do not run, while for the dimension-nine operators we have, 1 , C 4 , C where λ = 1/(2C F + 1/N c ) = 1/N c . Here the couplings C 1 and C (9) 4 increase in the UV, and the behavior of C T (2 GeV) = 0.87 C The remaining operators in Eqs. (7), (8), and (9) are scale independent at one loop in QCD.
The neutrino potentials in Eqs. (81) and (82) 34), they are not enhanced by the large nucleon isovector magnetic moment. Therefore we expect their contribution to be numerically somewhat smaller. In the case of Eq. (82), the second term was included in the analysis of Refs. [32,77], where it was found to be much smaller than the magnetic term in Eq. (35). For this reason, we neglected it in our formulae for the decay rate in Section 6.

D Conversion of nuclear matrix elements
In this appendix, we provide the conversion between the NMEs defined in Sec. 6.1 and those of the original papers [32,76,77,[83][84][85].
For the matrix elements involving the exchange of a light neutrino, our definitions match those in Refs. [76,[83][84][85]. The only exceptions are M M M GT,T , for which Refs. [76,84,85] used g M (0) = κ 1 = 3.7 rather than g M (0) = 1 + κ 1 . In Section 6.1, we thus rescaled these matrix elements by powers of r M = (1 + κ 1 )/κ 1 . For the Gamow-Teller and tensor matrix elements, Ref. [32] does not separately provide the AA, AP , P P and M M components. However, we can reconstruct the needed NMEs from linear combinations of other matrix elements computed in Ref. [32], as detailed in Table 9. The definitions of the NMEs in the third column of Table 9 are given in Ref. [32] 10 .
The relations we use are valid at LO in the chiral expansion, when one can takeĒ → 0 and neglect subleading effects as the difference between the axial and vector form factors. We discussed some checks of these assumptions in Sec. 6.1. Additional consistency checks can be performed with the NMEs of Ref. [32]. In the limitĒ → 0, one would expect M F = M F ω = M F q and M GT ω = M GT q . These relations are respected to a few percent for M F and M F ω , while M F q appears to be ∼ 50% smaller than M F ω . The relation between the GT elements holds to 10 The relation between M M M GT and MR given in Table 9 takes into account a factor of 1/3 that is missing from the definition of HR in Eq. (21v) of (the first arXiv version of) Ref. [32]. We thank M. Horoi for clarification on this issue.