A neutrinoless double beta decay master formula from effective field theory

We present a master formula describing the neutrinoless-double-beta decay ($0\nu\beta\beta$) rate induced by lepton-number-violating (LNV) operators up to dimension nine in the Standard Model Effective Field Theory. We provide an end-to-end framework connecting the possibly very high LNV scale to the nuclear scale, through a chain of effective field theories. Starting at the electroweak scale, we integrate out the heavy Standard Model degrees of freedom and we match to an $SU(3)_c\otimes U(1)_{\mathrm{em}}$ effective theory. After evolving the resulting effective Lagrangian to the QCD scale, we use chiral perturbation theory to derive the lepton-number-violating chiral Lagrangian. The chiral Lagrangian is used to derive the two-nucleon $0\nu\beta\beta$ transition operators to leading order in the chiral power counting. Based on renormalization arguments we show that in various cases short-range two-nucleon operators need to be enhanced to leading order. We show that all required nuclear matrix elements can be taken from existing calculations. Our final result is a master formula that describes the $0\nu\beta\beta$ rate in terms of phase-space factors, nuclear matrix elements, hadronic low-energy constants, QCD evolution factors, and high-energy LNV Wilson coefficients, including all the interference terms. Our master formula can be easily matched to any model where LNV originates at energy scales above the electroweak scale. As an explicit example, we match our formula to the minimal left-right-symmetric model in which contributions of operators of different dimension compete, and we discuss the resulting phenomenology.


Introduction
Neutrinoless double beta decay (0νββ) is the process where two neutrons inside an atomic nucleus are transmuted into two protons and two electrons without the emission of neutrinos. An observation of this process would indicate that lepton number (L) is not a good symmetry of nature and that the neutrino mass has a Majorana component, implying that the mass eigenfields are self-conjugate. Current experimental limits are very stringent [1][2][3][4][5][6][7][8][9][10][11][12], e.g. T 0ν 1/2 > 1.07 × 10 26 yr for 126 Xe [5], with next-generation ton-scale experiments aiming for one to two orders of magnitude improvement. The simplest interpretation of 0νββ experiments assumes that lepton-number violation (LNV) is due to the exchange of light Majorana neutrinos. However in various beyond-the-SM (BSM) scenarios other sources of LNV exist that can induce 0νββ. For example, in left-right symmetric models, apart from the exchange of a light Majorana neutrino, there appear contributions from the exchange of heavy neutrinos and charged scalars. While a single nonzero 0νββ measurement can be attributed to any LNV interaction, in principle various LNV sources can be disentangled by measurements of different isotopes, the angular or energy distributions of the outgoing electrons, or by correlating with collider observables, provided sufficient theoretical control can be achieved.
In most BSM scenarios the LNV source responsible for 0νββ is induced at an energy scale Λ well above the electroweak scale. This scale separation justifies an effective field theory (EFT) approach. Such an approach has the advantage that 0νββ and its correlation with collider observables can be described in a model-independent fashion. The Standard Model can be seen as the renormalizable part of an EFT that includes higher-dimensional operators which are suppressed by powers of the scale of BSM dynamics, Λ. Within this EFT, the ∆L = 2 operators have odd dimension [13]. The first of ∆L = 2 term therefore appears at dimension five and provides a contribution to the neutrino Majorana mass [14]. In the standard type-I see-saw mechanism this dimension-five operator arises from integrating out heavy right-handed neutrinos at the typical GUT-scale Λ ∼ 10 15 GeV. LNV operators with a dimension > 5 are then suppressed by powers of v/Λ 10 −13 and can be safely neglected. In various models, however, the scale of LNV new physics is much lower and the dimension-five operator may be suppressed by loop factors and/or small Yukawa couplings. For instance, in the above-mentioned left-right symmetric models the dimension-five operator scales as y 2 /Λ, where y is a Yukawa coupling scaling as y ∼ m e /v ∼ 10 −6 . While dimension-seven and -nine LNV operators are suppressed by additional powers of Λ, they can be suppressed by only one or even zero powers of y. As such, the dimension-seven and -nine operators can be competitive with the dimension-five operator, for Λ in the 1 − 10 TeV range. Higher-dimensional LNV operators at the multi-TeV scale can be generated in radiative neutrino models [15][16][17][18][19][20] and have also been studied in Refs. [20][21][22][23][24].
In order to describe 0νββ in a model-independent way, one should include all LNV operators up to dimension-nine in the SM-EFT 1 . In this work, we consider these operators and derive the form of the SU (3) c × U (1) em -invariant ∆L = 2 Lagrangian at a scale of a few GeV. We then construct the resulting chiral EFT Lagrangian that describes ∆L = 2 interactions between the low-energy degrees of freedom: pions, nucleons, electrons, and neutrinos. Armed with the chiral Lagrangian, we calculate the leading-order (LO) 0νββ transition operator or "0νββ potential" using a consistent power counting that was introduced in Refs. [25][26][27]. The 0νββ transition operator is the basis of many-body calculations of the 0νββ amplitude. We identify the nuclear matrix elements (NME) that need to be calculated in order to obtain the LO 0νββ decay rate for 0 + → 0 + transitions. Somewhat remarkably we find that the set of NMEs necessary to describe the LO 0νββ rate arising from LNV sources up to dimension nine is the same as the one required to describe the light and heavy Majorana neutrino exchange contributions. As such, the required NMEs have already been studied extensively in the literature and we can compare results from different many-body methods.
Our main result is a so-called "Master formula" for the 0νββ decay rate which includes the contributions from all LNV sources up to dimension nine and takes into account all possible interference terms. The formula incorporates the matching and evolution of the effective LNV operators that appear at the electroweak scale all the way to the final 0νββ decay rate by the application of several different EFTs, illustrated in Fig. 1. This Master formula can be used to directly calculate the 0νββ decay rates for different isotopes in any given particle physics model after matching to the SM-EFT at the electroweak scale. Our analysis also highlights the role played by uncertainties in hadronic and nuclear matrix elements, and we discuss how this is reflected on the limits we can set on BSM models.
To illustrate the usefulness of the Master formula we study a specific high-energy model, the left-right symmetric model. We show that, after performing a matching calculation to the SM-EFT operators, the master formula can be directly used to set constraints on the model parameters. In particular, we investigate for what values of the right-handed scale the dimensionseven and -nine operators are competitive with respect to the dimension-five contributions.
In order to keep this paper somewhat compact, we heavily borrow from previous work that described in detail various technical aspects of the derivation [25][26][27]. In that work, we used the same SM-EFT and chiral EFT techniques but limited ourselves to LNV dimension-five and -seven operators. Here we extend the analysis to dimension-nine and derive the 0νββ transition operators at LO in the chiral expansion. Nevertheless, the current paper is self-contained and all information to use the Master formula is given in the main text and accompanying appendices.
The organization of the paper is perhaps best described by referring to Fig. 1. The prediction for the 0νββ decay rate is obtained through a sequence of effective field theories, in order to separate short-distance particle physics effects which can be treated perturbatively, from longdistance hadronic and nuclear physics that must be evaluated non-perturbatively. Here are the key elements: • We assume throughout that lepton number violation occurs at a scale Λ m W . This scale is shown in green in Fig. 1. At this scale there are a number of SU (3) c × SU (2) L × U (1) Yinvariant ∆L = 2 operators, beginning at dimension five (the usual Weinberg operator). In this paper we also consider gauge-invariant ∆L = 2 operators in the SM-EFT with dimension seven and nine. We denote the dimension of SM-EFT operators by dim-n with n = 5, 7, 9.
• At scales below m W , after electroweak-symmetry breaking, we match to a new effective field theory in which ∆L = 2 manifests via operators of different dimension than the original dim-5, -7, -9 operators. This happens because once the Higgs field obtains its vacuum expectation value (vev) and the Higgs and W boson are integrated out of the EFT, we generate at tree-level new dimension-3, -6 ,-7, and -9 operators. We evolve these operators to the GeV scale by considering the one-loop QCD renormalization of such operators, which captures the leading-order perturbative operator mixing and renormalization. This step is shown in blue in Fig. 1.
• The green and blue steps are described in Section 2. The QCD renormalization group equations are given in Section 2 and solved in Appendix D.
• The next step, shown in red, occurs at the GeV scale, where the quark operators are matched to chiral EFT, describing the interactions of neutrons, protons, pions, and leptons. Here a number of effects occur: -The dimension- 3 Majorana neutrino mass generates the standard neutrino potential arising from ν → ν c (antineutrino conversion into neutrino) but also a contact operator nn → ppee. Both effects occur at leading order in the power counting [27].
-Dimension-6 operators generate new interactions leading to ∆L = 2 beta decay, n → peν, and pion decay, π → eν. When these are combined with the SM weak interactions, we find additional contributions to the nn → ppee transition operator. Renormalization also requires the introduction of contact nn → ppee operators, discussed in Appendix C.
-Dimension-9 operators can generate both short-range contact operators, nn → ppee, as well as pion-range, n → pπee and ππ → ee, interactions. When the pion-range interactions are combined with the strong interactions we generate additional contributions to the nn → ppee transition operator. The three processes are shown in Figure 2.
These effects are described in Section 3. The hadronic input (low energy constants, or "LECs") required for our analysis is provided in Table 1. An important hadronic effect is the non-perturbative renormalization of the Majorana neutrino mass [27] and higherdimension operators (see Figure 6), due to the short-range nature of the nucleon-nucleon force. These effects are discussed in Appendix B and the expectations for the values of the new LECs are summarized in Table 1.
• All these effects at the O(10 − 100 MeV) scale combine to form the nn → ppee transition operator, to be used in many-body nuclear calculations. This step is shown in purple in Fig. 1.
• The next step, shown in orange in Fig. 1, is to evaluate the 0νββ transition operators between the ground states of the relevant nuclei. The nuclear many-body matrix elements (see Table 2 and Appendix A.2), together with the phase space factors (see Table 4 and Appendix A.1), and short-distance Wilson coefficients, combine to give the 0 + → 0 + decay rate. This final step is described in Section 4.
• Section 5 then applies our Master formula to obtain bounds on the dim-3, -6,-7, and -9 operators, assuming only a single operator dominates at a time. In Section 6 we consider an explicit example of a UV-complete model, namely the left-right symmetric model [28][29][30], to illustrate the utility of our general formalism. In this model several operators of different dimension are generated. Further details on the left-right model (LRM) are provided in Appendix E.
• In Section 6.3 we discuss the case in which the right-handed neutrinos in the left-right symmetric model have masses close to the GeV scale. Although the assumption Λ v no longer applies in this case, we show that the EFT framework can be straightforwardly generalized to include LNV arising at intermediate scales.
• We summarize our results and discuss future directions in Section 7. Figure 1: A schematic overview of the effective field theory approach to evaluating the 0νββdecay amplitude starting from high-scale ∆L = 2 dynamics. The different colors represent various effective field theories at different scales. See the main text for more details.
Finally, before diving into the details of our analysis, let us comment on the relation of our work with prior literature. First of all, the framework discussed here is similar in spirit to Refs. [21,24], whose results we generalize and extend here. Secondly, a master formula for 0νββ has previously appeared in the literature [31,32], and it has been used in many subsequent studies. We note that we disagree with Refs. [31,32] in the treatment of the hadronization of quark operators. In particular, Ref. [31] misses important hadronic physics, thereby underestimating the contributions of certain Wilson coefficients to the 0νββ amplitude by O(16π 2 ) compared to what we find here. One main difference is that Ref. [31] ignores the couplings to ππ -which we now know with a fair amount of certainty [33][34][35][36] -and to πN . The LECs of certain four-nucleon operators are also underestimated by O(16π 2 ), because non-perturbative renormalization is not considered. We further discuss these and other differences with Ref. [31] in Appendix F.

Lepton number violation in the SM-EFT
Lepton number is an accidental symmetry of the renormalizable part of the SM, which is violated by higher-dimensional operators. The ∆L = 2 operators relevant for 0νββ all have odd dimension [13] and we focus on dimension-five, -seven, and -nine operators that, respectively, scale as Λ −1 , Λ −3 , and Λ −5 , where Λ is the scale at which lepton number violation arises. At lower energies, after electroweak symmetry breaking (EWSB) and integrating out heavy SM fields (top, Higgs-, W-, and Z-bosons) the arising effective operators can have a different canonical dimension due to positive powers of the Higgs vacuum expectation value, v 246 GeV (the SM-EFT approach assumes Λ v). In particular, at energies around a few GeV the most important ∆L = 2 operators have canonical dimension three, six, seven, and nine. To avoid confusion, when discussing the original gauge-invariant SM-EFT ∆L = 2 operators, we denote their dimensions by dim-n with n = 5, 7, 9. When discussing the operators after EWSB, which are only SU (3) c × U (1) em invariant, we refer to them as dim-n operators (without the overline) where n = 3, 6, 7, 9.
The structure of the gauge-invariant ∆L = 2 operators has been discussed in great detail in the literature [14,[20][21][22][23][24]. The only dim-5 operator is the Weinberg operator [14] which, after EWSB, gives rise to the neutrino Majorana mass. The 12 operators that appear at dim-7 have been classified in Ref. [23] and were studied in the context of 0νββ in detail in Ref. [25]. The complete set of dim-9 operators is currently unknown, but certain subclasses with particular field content have been constructed [22,24]. For instance, Ref. [24] identified the dim-9 operators consisting of four quark fields and two electron fields, finding that only eleven such operators exist. However, additional classes of dim-9 operators can be constructed that give rise to unsuppressed 0νββ. Examples are operators involving two quark fields, two electron fields, and the combination of Higgs fields and derivativesφ † D µ ϕ. While these operators require the exchange of a W boson to induce 0νββ, the associated factor of G F is compensated by two powers of v arising from the two ϕ fields after EWSB.
Here we do not list the gauge-invariant dim-n operators but refer to Refs. [14,[20][21][22][23][24][25] for more details. Instead we focus on the ∆L = 2 Lagrangian after EWSB and integrating out the heavy SM degrees of freedom. At a scale slightly below the electroweak scale, the Lagrangian consists of SU (3) c × U (1) em operators of increasing dimension. For applications to 0νββ it is convenient to organize the Lagrangian in operators that violate the number of charged leptons by zero, one, or two units L ∆e=0 contains operators that violate lepton number in the neutrino sector, starting from the dim-3 Majorana mass of left-handed neutrinos The SU (2) L × U (1) Y invariance of the SM implies that the first contribution to m ν arises from a dim-5 operator, such that m ν ∼ v 2 /Λ. The dots in Eq. (2) denote operators of higher dimension, such as the dim-5 neutrino magnetic moment or dim-6 LNV neutral-current semileptonic operators. In order to induce 0νββ, the two neutrinos in the operators in L ∆e=0 need to be converted into electrons via the SM weak interaction. The contributions to 0νββ from higherdimensional operators in Eq. (2) are thus suppressed at least by powers of Λ 2 χ /v 2 (if not m 2 π /v 2 ), where Λ χ ∼ 1 GeV is the chiral-symmetry-breaking scale [37], with respect to m ν . We therefore neglect these effects in this work.
A richer set of contributions arises from L ∆e=1 . This Lagrangian contains LNV operators with one charged lepton and one neutrino field. In order to compensate the charge of the electron field, one needs at least an additional quark or lepton bilinear, making dim-6 the minimal dimension of these operators: ∆L=2 + L ∆L=2 + . . .
The operators most relevant to 0νββ are semileptonic four-fermion operators. At dim-6 we have ∆L=2 contains all possible ∆L = 2 dim-6 charged-current operators. At tree level, all operators in Eq. (4) receive their first contributions from dim-7 operators [25], so that C Beyond tree level, the operators in Eq. (4) might also receive contributions of O(v/Λ) from the neutrino Majorana mass, but we neglect these loop corrections here. Dim-7 operators in L ∆e=1 give rise to corrections that are suppressed by Λ χ /v with respect to the dim-6 terms of Eq. (4).
Here we consider only the subset of SU (3) c × U (1) em invariant dim-7 operators that receive tree-level matching coefficients at the EW scale from dim-7 operators [25] The coefficients of these operators scale as C Operators of higher dimension in Eq. (3), such as dim-8 dipole operators or dim-9 charged-current six-fermion operators, give rise to contributions that are more and more suppressed by powers of Λ χ /v and v/Λ.
The final class of operators are LNV operators with two electrons, which can directly contribute to 0νββ without additional SM weak interactions. U (1) em invariance forces these operators to be at least dim-9 L ∆e=2 = L ∆L=2 + . . .
The set of SU (3) c × U (1) em invariant four-quark two-lepton operators can be written as [21,24] L (9) where O i and O µ i are four-quark operators that are Lorentz scalars and vectors, respectively. The scalar operators have been discussed in Refs. [21,24] and can be written as where τ ± = (τ 1 ± iτ 2 )/2 with τ i the Pauli matrices and α, β are color indices. The O i operators are related to the O i by parity. The vector operators take the form [24] O µ where the second column of operators is related to the first column by a parity transformation. While twenty-four dim-9 operators appear in Eq. (7), not all of them are necessarily induced by the gauge-invariant dim-n (for n ≤ 9) operators. As discussed in Ref. [25] only three dim-9 operators, those with Wilson coefficients C (9) 1L , C (9) 4L and C (9) 5L , receive tree-level matching from dim-7 operators. This implies that these coefficients can scale as C Ref. [24] showed that the eleven dim-9 operators with Wilson coefficients C (9) 2L, 3L, 4L, 5L , C 2L, 3L , C 1R , and C (9) 6,7,8,9 receive contributions from four-quark two-lepton dim-9 operators. As discussed above, by replacingd R γ µ u R →φ † D µ ϕ in the dim-9 operators of Ref. [24] we identify several additional dim-9 operators that can be induced by dim-9 operators, namely C (9) 1L , C 1R, 4R, 5R and C (9) 6,7,8,9 . To summarize, out of the twenty-four dim-9 operators in Eq. (7) we find that nineteen operators are actually induced by dim-n (for n ≤ 9) operators. U (1) Y invariance implies that there are no dim-9 operators that contribute to C (9) 1L , C 2R, 3R , and C (9) 2R, 3R . These operators are then induced only at the dim-11 or higher level. Although this means these operators should be further suppressed, we include them for completeness, thereby keeping all operators in Eq. (7) in our analysis. Finally, the effective SU (3) c × U (1) em operators we consider may not cover the contributions of all dim-9 terms as, unlike for dim-5 and dim-7 operators, a complete dim-9 basis is currently unavailable. Nevertheless, we expect the included SU (3) c × U (1) em operators to capture the dominant ∆L = 2 effects in most, if not all, models of LNV.
The coefficients in Eqs. (4), (5), and (7) 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 corresponding to C (6) VL, VR and C (7) VL, VR involve quark nonsinglet axial and vector currents and therefore do not run in QCD. The renormalization group equations (RGEs) of the scalar and tensor operators below µ = m W are given by where C F = (N 2 c − 1)/(2N c ) and N c the number of colors. The RGEs of the scalar dim-9 operators are given by in agreement with Refs. [38,39]. The RGEs do not depend on the lepton chirality, and we therefore omitted the subscripts L, R in Eq. (15). The equations for the C 1,2,3 coefficients are equivalent to those in Eq. (15), while the RGEs for the vector operators can be written as, The evolution of C 8 and C 9 , as well as that of the primed coefficients, are governed by RGEs of the same form. The (numerical) solutions to the above RGEs are given in Appendix D and a comparison to the literature is given in Appendix F.

The chiral Lagrangian
After obtaining the Lagrangian at the scale of a few GeV, the next step is to match the quarklevel theory onto chiral perturbation theory (χPT) [40,41]. χPT is the effective theory of QCD which describes interactions at low energy in terms of baryons, mesons, photons, and leptons. It relies on an expansion in χ ≡ p/Λ χ that is provided by the approximate chiral symmetry of QCD. Here p ∼ m π is the typical momentum and Λ χ 1 GeV. The interactions that are induced in the chiral Lagrangian can be derived by writing down all possible terms that have the same chiral symmetry properties under SU (2) L × SU (2) R as the quark-level operators. All such operators come with an unknown coefficient that parametrizes the non-perturbative nature of QCD. In the mesonic and single nucleon sector of the theory, these low-energy constants (LECs) can be estimated by naive dimensional analysis (NDA), and, as we will discuss, several of them can be extracted from experimental data or existing lattice QCD calculations. The power counting of LNV nucleon-nucleon operators is complicated by the non-perturbative nature of the nuclear force, which causes NDA estimates of the coefficients of nucleon-nucleon operators to be unreliable [27,[42][43][44]. We will determine the scaling of nucleon-nucleon operators by requiring that the LNV nn → ppee scattering amplitude is properly renormalized [27].
The chiral Lagrangian induced by the neutrino Majorana mass and by the operators in Eqs. (4) and (5) has been discussed in Refs. [25][26][27]. In Appendix C we summarize the main results, and discuss several short-distance effects that were overlooked in Ref. [25]. We focus here on the dim-9 operators in Eq. (7), for which partial results are given in Refs. [21,24]. For the dim-9 operators under discussion three types of interactions will appear: purely pionic interactions, pion-nucleon couplings, and nucleon-nucleon (NN ) terms. In what follows we construct the corresponding operators that give the LO contributions to the nn → ppee amplitude.
Unlike standard χPT applications, in the operators constructed below there appear explicit lepton fields. As such, operators can be constructed that depend on the electron mass and/or the momenta of the outgoing lepton fields. The typical Q-values, and the corresponding electron momenta, of experimentally measurable 0νββ-transitions are small, Q ∼ O(5 MeV), compared to the typical momentum exchange between nucleons q ∼ k F ∼ m π = O(100 MeV). In order to incorporate this new scale into the χPT power counting we assign Q ∼ m e ∼ m π 2 χ [25].

Scalar dim-9 operators
The scalar operators O 1 -O 5 generate the ππee, πNN ee, and NN NN ee LNV couplings shown in Fig. 2. The operators O 2,3,4,5 induce non-derivative pionic operators, while the first pionic operators induced by O 1 contain two derivatives and are therefore relatively suppressed. The mesonic chiral Lagrangian 2 for O 1,2,3,4,5 is where U = u 2 = exp (iπ · τ /F 0 ) is the matrix of pseudo-Goldstone boson fields, F 0 is the pion decay constant in the chiral limit, and L µ = iU D µ U † . We use F π = 92.2 MeV for the physical pion decay constant. By NDA the LECs of the non-derivative pion operators are expected to be g ππ 2,3,4,5 = O(Λ 2 χ ), while g ππ 1 = O(1). These expectations are very well respected by the extractions of Ref. [34][35][36] based on chiral symmetry and lattice data. In Table 1 we give the value of the LECs at µ = 2 GeV in the MS scheme, obtained in Ref. [36]. The physical amplitudes are scale and scheme independent provided one uses Wilson coefficients C (9) i evaluated at the same scale and in the same scheme as used for the g ππ i . The πN terms are only relevant for the O 1 operator and can be written as 27, N = (p, n) T , and S µ and v µ are the nucleon spin and velocity. In the nucleon restframe we have S α = (0, σ/2) and v µ = (1, 0). The LEC g πN 1 is unknown, but expected to be O(1) by NDA. In a power counting based on NDA, LNV four-nucleon interactions are relevant only for O 1 , in which case they would compete with the ππ and πN interactions g ππ 1 and g πN 1 . However, the LNV potential induced by the non-derivative ππ operators in Eq. (17) has the same short-distance behavior as the neutrino potential mediated by the neutrino Majorana mass, 2 The ππ couplings defined here are related to those of Refs. [25,35,45] by g ππ 1 = g27×1, g ππ 2 = g 6×6 , g ππ 3 = g mix 6×6 , g ππ 4 = g8×8, g ππ 5 = g mix 8×8 , while for the πN and N N couplings we have g πN 1 = g πN 27×1 and g N N 1 = g N N 27×1 . The notation of Refs. [25,35,45]  V (q) ∼ 1/q 2 at large |q|. Ref. [27] showed that for these potentials the nn → ppee scattering amplitude has a logarithmic UV divergence, which can be absorbed by promoting the NN operators stemming from O 2,3,4,5 to leading order. The relevant NN interactions are The couplings g N N i = O(1) in the Weinberg power counting, but need to be promoted to O((4π) 2 ) in the case of O 2,3,4,5 . The renormalization of the scattering amplitude does not require such enhancement for g N N 1 . The ππ, πN , and NN Lagrangians for the O 1,2,3 operators, which are related by parity to O 1,2,3 , can be obtained by replacing C (9) 1L, 2L, 3L → C (9) 1L, 2L, 3L , C 1R, 2R, 3R → C (19). 3 This leads to ππ, πN , and N N Lagrangians of the same form (with C (L,R) 1,2,3 → C (L,R) 1,2,3 ) after expanding the meson matrices u (and U ) to two, one, and zero pions, respectively.
From Eq. (19) we see that all scalar operators in Eq. (8) induce a LNV four-nucleon operator that contributes to the nn → ppee amplitude at the same order as the pion-range contributions from the ππee operators. This happens either because the ππee interaction is suppressed by two powers in the chiral counting (as it is for O 1 and O 1 ), or because of non-perturbative renormalization which enhances the four-nucleon operator to leading order (O 2,3,4,5 and O 2,3 ). In other words, for all scalar operators the ππee and N N interactions appear at the same order. For O 1 and O 1 there appear additional contributions from the πN interaction. More details on the renormalization of the nn → ppee scattering amplitude and the non-perturbative RGE of g N N i are given in Appendix B.

Vector operators
The vector operators induce mesonic interactions involving a derivative on the pion fields, which, up to a total derivative, give rise to contributions proportional to m e [21,24]. Instead, the contributions to the nn → ppee amplitude from the πN and NN interactions are proportional to |q| ∼ k F and therefore larger than the contributions from the purely mesonic terms by one power of 1/ χ . The πN Lagrangian induced by O µ 6,7,8,9 and O µ 6,7,8,9 can be written as Tr iD α U † τ + − g πN 8 C 8 + g πN 9 C (9) 9 Tr where g πN i are LECs of O(1). The chiral representations of the operators O µ 6,7,8,9 are related to the above Lagrangian by parity. The corresponding πN terms can be obtained by replacing Combining all contributions and expanding up to one pion we have where C V ≡ C 6 + C 8 + C (9) 6 + C 8 ,C V ≡ C and g πN V = g πN For the NN interactions induced by O µ 6,7,8,9 and O µ 6,7,8,9 we obtain where g N N 6,7 are LECs of O(1) within the Weinberg power counting. The relevant Lagrangian for the O µ 6,7,8,9 operators is again related by parity, and can be obtained by replacing L ↔ R and u → u † in Eq. (22). Expanding the NN interactions gives for terms without pions The non-perturbative renormalization of the scattering amplitude relates the πN and NN couplings through an RGE, see Appendix B, but this does not change the power-counting expectations.
To summarize, the chiral Lagrangian from the dim-9 vector is rather simple and at LO consists of only the interactions in Eqs. (21) and (23). Unfortunately, very little is known about the corresponding LECs.

The 0νββ transition operator
With the above chiral operators we can construct the LO two-nucleon operators that induce 0νββ, to which we refer to as the "0νββ transition operators" or "0νββ potentials", as commonly done in the literature. The calculation can be directly lifted from Ref. [25] as all chiral Table 1: The low-energy constants relevant for the dim-3, dim-6, dim-7 and dim-9 operators. Whenever known, we quote the values of the LECs at µ = 2 GeV in the MS scheme. The LECs g ππ 1,...,5 were first extracted in Refs. [33,35] using SU (3) relations between π + -π − , K-K and K → ππ matrix elements. Ref. [35] found g ππ in good agreement with the more precise results of Ref. [36].
interactions constructed above also appear in the dim-7 chiral Lagragian but, of course, accompanied by different Wilson coefficients and LECs. We refer to Ref. [25] for details and give here the results. We define the ∆L = 2 transition operator or potential as V = −A where A is the Born-level amplitude of nn → pp ee.
The LO ∆L = 2 0νββ potential from the scalar and vector dim-9 operators is Here the combinations C (9) ππ, πN, N N of scalar couplings are defined as and similarly for C (9) {ππ, πN, N N } R . In the above expressions we use q µ = (q 0 , q) = (p − p ) µ , where 2p and 2p are the relative momenta of the in-and outgoing nucleon pairs. We defined The potential in Eq. (24) can be expressed in terms of the Fermi (F), Gamow-Teller (GT), and Tensor (T) functions h F (q 2 ), h AP,P P GT (q 2 ), and h AP,P P T (q 2 ) defined in Appendix A, and, consequently, their matrix elements for experimentally interesting nuclei can be extracted from existing calculations [49][50][51][52][53].
4 0νββ amplitudes and master formula for the decay rate 4.1 0νββ amplitudes Now that we have derived the momentum-space 0νββ potential, we can obtain an expression for the inverse half-life for 0 + → 0 + transitions. This requires one to insert the two-nucleon potential between the initial-and final-state nuclear wave functions and sum over all nucleons. In the literature, many-body matrix elements are usually calculated in coordinate space and we therefore define the amplitude where we sum over all nucleon pairs, q is the momentum transfer, and r = r n −r m is the distance between the m th and n th nucleon. The 0νββ potential V (q 2 ) is the sum of potentials induced by the low-energy ∆L = 2 dim-3, -6, -7, and -9 operators and we write V 3 (q 2 ), often called the "neutrino potential" in the literature, is the contribution from light Majorana-neutrino exchange. We give its expression in Eq. (81). The potentials corresponding to higher-dimensional operators are V 6,7 (q 2 ), which were derived in Ref. [25] and are for convenience reproduced in Appendix C, and V 9 (q 2 ), given in Eq. (24). It is useful to separate the total amplitude in Eq. (26) in various terms which differ by their leptonic structure. While V 9 (q 2 ) only consists of three different leptonic structures, two additional structures appear in V 3,6,7 (q 2 ) and we define where E 1,2 (k 1,2 ) are the energies (momenta) of the electrons. Here we extracted an overall factor from the various sub-amplitudes A i . In particular, a factor of m e /R A is extracted, where m e is the electron mass and R A = 1.2 A 1/3 fm in terms of A, the number of nucleons of the daughter nucleus. This normalization was chosen in order to align the definition of the various nuclear matrix elements with those appearing in the literature, but stress that in the final decay rate all the factors of m e /R A will drop out.
The subamplitudes A i depend on the Wilson coefficients of the ∆L = 2 operators, on hadronic matrix elements, and nuclear matrix elements. The required LECs encoding hadronic matrix elements are listed in Table 1. It turns out that all nuclear input that appears in Eq. (28)  For the exact definitions we refer to Appendix A.2. All NMEs, apart from one (M AA T ), can be extracted from existing calculations of light-and heavy Majorana-neutrino exchange contributions. Furthermore, at LO in χPT the fifteen NMEs are related by five identities that can be used to further reduce the number of required many-body calculations or as a consistency check of the results [25]. In Table 2 we summarize several recent calculations of the NMEs, obtained by different groups applying different many-body methods. The NMEs often appear in certain linear combinations M i that are defined below.
It is useful to further decompose the sub-amplitudes in terms of contributions from LNV operators of different dimension R , M .
The subamplitude A ν multiplies the leptonic structure that arises from light Majorananeutrino exchange, from several long-range dim-6 and dim-7 contributions, and from shortrange dim-9 contributions. We have therefore decomposed it in a component proportional to the electron-neutrino Majorana mass m ββ , and the additional terms M (6) ν and M (9) ν , generated, respectively, by dim-6 and -7, and by dim-9 LNV operators. The short-distance component M (9) ν arises from V 9 and always involves an additional power of 1/v with respect to the contribution from light Majorana-neutrino exchange. To compensate for this factor and for the absence of the neutrino mass, we have factored out two powers of m N in Eq. (29). In terms of the standard building blocks defined in Appendix A.2, the combination of NMEs M i are defined as where is a new leading-order low-energy constant [27], defined in Eq. (80), and B ≡ − qq /F 2 π 2.7 GeV at µ = 2 GeV in the MS scheme. Only C (9) 1L and C (9) 4L, 5L in M (9) ν receive matching contributions from dim-7 operators [25], while the remaining terms are at least dim-9 [24]. In the above expressions we have defined in terms of matrix elements defined in Appendix A.2. g ππ T , g πN T and g N N T are the LECs of ππee, πN N ee and N N N N ee short-range operators induced by C (6) T , defined in Appendix C. The subamplitude A R only receives contributions from the dim-9 scalar operators involving right-handed electrons. It is only induced by dim-9 operators and is proportional to The subamplitudes A E and A me are not affected by dim-9 operators and their expressions are therefore the same as in Ref. [25], apart from additional short-range contributions that were not included there. They depend on the NME combinations where g E,me V L,V R are defined in Appendix C. The final subamplitude, A M , receives a contribution from the dim-6 operator C (6) VL and from the dim-9 vector operators in Eq. (7) which are, respectively, proportional to The πN N ee and N N N N ee couplings g πN VL and g N N VL are defined in Appendix C. The contributions from the various LNV operators to the subamplitudes A i can be organized according to their scaling in the χPT expansion parameter χ = m π /Λ χ , where Λ χ ∼ 4πF π . The scaling is summarized in Table 3. The Table indicates that low-energy dim-7 and -9 contributions are suppressed by at least one power of Λ χ /v compared to the dim-6 contributions. Because the dim-9 operators C (9)( ) 2,3,4,5 can induce the pionic operators in Eq. (17), their contribution is enhanced by two powers of χ with respect to the dim-7 operators and the dim-9 scalar operators C (9)( ) 1 and all vector operators C (9) V . In Table 2 we report recent evaluations of the NMEs for experimentally interesting 0νββ emitters, obtained with three different many-body methods: the quasi particle random phase approximation (QRPA) [49], the shell model [50], and the interacting boson model (IBM) [51,52].
Refs. [49][50][51][52] Table 2 we see that these expectations are well respected by the Fermi and Gamow-Teller NMEs. M M M GT is larger than expected, which can be understood by taking into account that the 2 χ suppression is partially compensated by the large isovector nucleon magnetic moment, g M 4.7. The tensor matrix elements are usually smaller because the tensor operator S (12) vanishes between nn pairs in the 1 S 0 channel, which is the dominant two-nucleon component [54]. Part of this suppression might be an artifact of the applied many-body methods, as Variational Monte Carlo calculations in lighter nuclei, such as 12 Be and 12 C, show the ratio M AP T /M AP GT to be roughly 25% [45].

Master formula for the 0νββ decay rate
Using the amplitude in Eq. (26), the expression for the inverse half-life becomes [55,56], Here M i is the mass of the decaying nucleus, while E 1,2 and E f are the energies of the electrons and final daughter nucleus in the rest frame of the decaying nucleus. The functions F (Z, E i ) are defined in Appendix A.1 and take into account the fact that the emitted electrons feel the Coulomb potential of the daughter nucleus and are therefore not plane waves.
Using the decomposition of the amplitude in Eq. (28) to separate the different leptonic structures, we obtain the final expression T 0ν This 'Master-formula' describes the 0νββ decay rate up to dim-9 operators in the SM-EFT. It includes all contributions from the low-energy ∆L = 2 operators in Eq. (1) and takes into account all interference terms. It should provide a useful tool to constrain any model of highscale LNV, using the most up-to-date hadronic and nuclear input. A differential version of Eq. (38) is given in Appendix A.1. The various components in Eq. (38) can be obtained as follows: • G 0i are phase space factors defined in Appendix A.1 and their numerical values are given in Table 4.
•   Table 2: Comparison of NMEs computed in the quasi particle random phase approximation [49], shell model [50], and interacting boson model [51,52] Table 3: Power-counting estimates of the contribution of low-energy dim-3, -6, -7, and -9 operators to the amplitudes in Eq. (28), in terms of m ββ , the Higgs vev v, and χ ≡ m π /Λ χ , where Λ χ ∼ m N ∼ 1 GeV. We take the electron mass and energies to scale as χ . This Table assumes the NMEs to follow the chiral EFT power counting. C (9) vector indicates any of the vector operators in Eq. (7). Finally, note that to estimate the overall scaling of the amplitudes one needs to take into account that, up to insertions of dimensionless couplings, the Wilson coefficients scale as follows: on the underlying LNV model; (ii) hadronic LECs, whose current knowledge is summarized in Tab. 1; (iii) nuclear matrix elements defined in Appendix A.2, whose numerical values from different many-body methods can be found in Tab. 2.
• Several hadronic LECs are at the moment unknown. An assessment of the ensuing theoretical uncertainty can be obtained by varying the LECs in a range around the values of Table 1. We stress that for all the operators in Eqs. (4), (5), and (7) On a more technical note, it should be stressed that the decay-rate formula is expressed in terms of the Wilson coefficients at a low-energy scale µ 2 GeV. In order to match the formula to specific BSM theories, some additional steps are required. At the high-energy scale Λ where any beyond-the-SM fields are integrated out, we need to perform a matching calculation to gaugeinvariant dim-5, dim-7, and dim-9 operators. The resulting operators need to be evolved down to the electroweak scale where they are matched to the operators in Eq. (1). This procedure was completed in Ref. [25] for the dim-7 operators, and below we study it for a particular BSM model where also relevant dim-9 operators are induced. Finally, the low-energy EFT operators are evolved to the low-energy scale using the RGEs in Eqs. (14)- (16). The numerical factors of the last step are given in Appendix D. All the steps leading from a generic LNV Lagrangian at scale Λ to Eq. (38) are illustrated in Fig. 1

Single-coupling constraints
We now investigate the constraints from 0νββ limits on the low-energy ∆L = 2 operators in Eq. (1). In particular, we apply the experimental limits [5,10,11] (all at 90% c.l.) For operators of dimension six and higher, we interpret the limit as a lower bound on the scale of new ∆L = 2 physics, Λ. The operators in the left column of Table 5 can be induced by dim-7 operators and we assume C The operators in the right column can only be induced by dim-9 operators and here we assume C As such, the probed scale in the left column is typically significantly higher, O(100 TeV), than in the right column, O(5 TeV). However, in cases where the dim-7 operators predominantly induce dim-7 operators, the additional suppression of 2 χ Λ χ /v in Table 3 makes the probed scale much lower and comparable to that in the right column.
We give the bounds in two cases. The top panel of Table 5 shows the limits obtained assuming that only one operator is active at the scale µ = 2 GeV. To highlight the impact of the QCD evolution, in the lower panel of Table 5 we show the limits in the assumption that the operators are turned on at µ = m W = 80.4 GeV. We can see that the QCD running gives O(1) corrections to the bounds. We thus find that the RGEs have a far milder effect than the O(10 3 ) effects that were found for some operators in Ref. [58]. The origin of this discrepancy is discussed in more detail in Appendix F.
In order to set these limits we had to make several assumptions. Firstly, we used the NMEs from Ref. [50]. Results from other groups and many-body methods roughly differ by factor of 2 to 3, depending on the NME under consideration. In particular, for the light Majorana-neutrino exchange the relevant NME is M (3) ν which differs by roughly a factor 2 between Refs. [49][50][51]53] and this impacts the limit on m ββ by the same amount. For the operators that scale as v 3 /Λ 3 or v 5 /Λ 5 the NME uncertainties give an uncertainty on Λ of roughly a factor 3 1/3 1.5 and 3 1/5 1.25, respectively.
The remaining uncertainty arises from the size of the LECs, in particular those associated with the ∆L = 2 pion-nucleon and nucleon-nucleon couplings. For the light Majorana-neutrino exchange in Weinberg's counting there appear no LO nucleon-nucleon LECs and there would only be a small uncertainty from higher-order chiral corrections. However, as demonstrated in Ref. [27], renormalization requires that a ∆L = 2 nucleon-nucleon operator is promoted to LO. Currently the contribution from this term has not been incorporated consistently in calculations for the heavy nuclei under considerations but estimates for light nuclei show that the nucleon-nucleon terms can alter the total amplitude by O(1) corrections [27], and the associated uncertainty is as large as the NME uncertainty.
Similar uncertainties affect the limits on the dim-7 and dim-9 operators. For the limits in Table 5 we assumed g N N We furthermore assumed that the pion-exchange contributions saturate the amplitude for the operators which induce the pionic ∆L = 2 operators in Eq. (17). That is, we assumed Weinberg's counting and neglected the nucleon-nucleon contributions, g N N ν = g N N 2,3,4,5 = g E,me V L,V R = 0, even though these terms are enhanced to LO by renormalization arguments. In addition, for the short-distance operators induced by the dimension-six operators, we assumed g ππ,πN,N N T = g πN,N N VL = 0. Our limits are therefore affected by O(1) uncertainties and should be revisited once more is known about the LECs.

An explicit example: the minimal left-right symmetric model
The formula for the total decay rate in Eq. (38) is given in terms of Wilson coefficients of effective operators. This formula can be matched to any model of heavy BSM physics that contributes to 0νββ. In what follows we demonstrate this by considering an explicit BSM model: the minimal left-right symmetric model (mLRSM) [28][29][30]. The mLRSM has been studied in the context of 0νββ in great detail [21,[59][60][61][62]. This scenario is interesting because it provides an elegant explanation of P-or C-violation through spontaneous symmetry breaking and it allows for the generation of neutrino masses at a relatively low scale, within reach of the LHC or possible future colliders. Furthermore, within the model there appear dim-5, dim-7, and dim-9 contributions to 0νββ. As such, the mLRSM is particularly well suited to demonstrate that Eq. (38) is able to capture all of these effects. We start by giving a brief overview of the model, and refer to Refs. [28-30, 63, 64] and Appendix E for a more detailed discussion. We stress that we do not aim to perform a full study of the model. Our main goal here is to illustrate the EFT framework and in particular the use of Eq. (38).
The model is based on the gauge group SU (3) c × SU (2) L × SU (2) R × U (1) B−L and we adopt the version of the model in which charge conjugation is conserved at high energies 4 , see e.g. Ref. [65]. The fermions are assigned to representations of the above gauge group as follows The introduction of right-handed neutrinos, ν R , is required by the symmetries of the model. In addition to the fermions, the mLRSM involves several scalar fields, namely a bidoublet transforming as φ ∈ (1, 2, 2 * , 0), as well as two triplet fields, ∆ L,R assigned to (1, 3, 1, 2) and 76 i , related to the dim-6, dim-7, and dim-9 operators from the GERDA [11], CUORE [10], and KamLAND-Zen [5] experiments. In the left column, we assume C 3 for the dim-6, dim-7, and dim-9 operators that receive contributions from gauge-invariant dim-7 operators at the electroweak scale. For the remaining dim-9 operators we assume C (9) i (µ) = (v/Λ (9) i ) 5 in the right column. We use the nuclear matrix elements of Ref. [50] and the values of the LECs described in the text. To illustrate the effect of the QCD running, we show the limits in two cases. In the upper panel of the Table we assume that the operators are turned on one at a time at the scale µ = 2 GeV, while in the lower panel we take µ = m W .
(1, 1, 3, 2), respectively. These fields can be written as and transform as φ → U L φU † R , ∆ L,R → U L,R ∆ L,R U † L,R under SU (2) L,R transformations. The neutral components of these fields obtain the following vacuum expectation values (vevs) In the first step of symmetry breaking the vev of the right-handed triplet, v R , breaks the chargeconjugation symmetry or parity and SU (2 This vev defines the high scale of the model and gives the main contribution to the masses of the right-handed gauge bosons, the right-handed neutrinos, and the heavy Higgs fields. The vevs of the bidoublet, κ and κ e iα , break SU (2) L × U (1) Y to U (1) em , and are of the order of the electroweak scale. v L contributes to the masses of the light neutrinos and to the ρ parameter, and is therefore required to be much smaller than the other vevs. Apart from the kinetic terms and the scalar potential, the Lagrangian contains interactions between fermions and scalars that give rise to the ∆L = 2 operators of interest in terms ofφ = τ 2 φ * τ 2 (which also transforms asφ → U Lφ U † R ), two quark Yukawa matrices, Γ andΓ, two lepton Yukawa matrices, Γ l andΓ l , and the symmetric 3 × 3 matrices M L and M R . Charge-conjugation (parity) invariance implies that the latter two matrices are related by , while the Yukawa matrices are forced to be symmetric (hermitian). After the scalar fields obtain their vevs, the Yukawa terms induce the Dirac masses for the quarks and leptons, while the M L,R terms induce neutrino Majorana masses and ∆L = 2 interactions. In the case of C-invariant Yukawa couplings, there is a relation between the left-and righthanded CKM matrices, V L and V R , namely, V L = K u V * R K d , where K u,d are diagonal matrices of phases. Instead, if one assumes parity to be preserved at high energies, V R can be expressed in terms of quark masses, V L , κ /κ, and sin α [66].

Matching to SM-EFT operators
We integrate out the heavy fields (with masses ∼ v R ) at a scale just below the one at which ∆ R obtains a vev, so that SU (2) R is broken while SU (2) L is intact. The heavy fields that contribute to ∆L = 2 interactions are the right-handed neutrinos and gauge fields, ν R and W R , as well as the heavy scalar fields, ϕ H , δ ++ R , and ∆ L 5 . At the scale µ m W R , the mLRSM then gives matching contributions to the following gauge-invariant effective Lagrangian Leudϕd R γ µ u R + C LϕDe mn ϕ m (D µ ϕ) n eeϕD ((iD µ ϕ) †φ ) 2 , 5 Here ϕH refers to the combination of SU (2)L doublets in φ that obtains a mass of O(vR), while we identify the remaining doublet with the SM Higgs doublet, denoted by ϕ, see Appendix E for details. Note that the scalar fields δ + R and δ 0 R do not mediate ∆L = 2 effects since they are absorbed in the now massive SU (2)R gauge bosons.
where the first, second, and third lines correspond to ∆L = 2 dim-5, dim-7, and dim-9 operators, respectively. We followed Ref. [25] for the definition of the dim-7 operators. The matching condition for the dim-5 operator is where M ν R = √ 2v R M † R is the mass matrix for the right-handed neutrinos and M D = is the Dirac Yukawa matrix, with ξ = κ /κ. Here the first term arises from integrating out the heavy right-handed neutrinos, corresponding to the usual type-I see-saw mechanism. The second term arises from a type-II see-saw mechanism and is induced by integrating out the ∆ L fields. This contribution is proportional to v R β i /m 2 ∆ L , where β i are the parameters in the scalar potential that couple ∆ L to the SM Higgs doublet. By use of the minimum equations this can be written in terms of the vev v L , as was done in Eq. (45), see Appendix E for details.
The dim-7 operators are induced by combining the right-handed neutrino Majorana mass term with the right-handed charged-current interaction. The W R boson can either couple to right-handed quarks or to (φ † D µ ϕ), which leads to mixing with W L proportional to ξ after EWSB. In total, we have The dim-9 operators result from diagrams in which either ν R or δ ++ R couples to two righthanded electrons and two W R bosons. Again, the W R bosons can couple to right-handed quarks, or mix with the W L boson and we obtain The above matching conditions hold at the scale m W R , so that the QCD running between the right-handed and the electroweak scale has to be taken into account before integrating out the heavy SM fields. Of the induced operators in Eq. (44), only C eeud is affected by QCD evolution. The other operators do not involve colored particles or consist of currents such that their QCD anomalous dimensions vanish. C (9) eeud follows the same RGEs as C 1 in Eq. (15), which leads to, This gives C   4 , with respect to the dim-5 operator. The smallness of M D can therefore compensate for the powers of (v/v R ) 2 . These estimates suggest that the dim-7 and dim-9 cannot be neglected for right-handed scales within reach of collider experiments.
6.2 Matching to SU (3) c × U (1) em -invariant operators and the 0νββ decay rate After integrating out the heavy SM fields at the electroweak scale the SU (2) L -invariant operators induce the following effective neutrino mass At dim-6 we have the following contributions [25], Finally, at dim-9 we obtain These operators still have to be evolved to µ 2 GeV. The running of the dim-9 operators is discussed in Section 2 and induces C

Matching in case of light right-handed neutrinos
So far we have been assuming that the right-handed neutrinos have O(v R ) masses and can be integrated out simultaneously with the right-handed W boson. In principle, it is possible for the right-handed neutrinos to have masses well below the right-handed scale if the M R couplings in Eq. (43) are small. In this case the right-handed neutrinos need to be integrated out separately from the BSM particles with O(v R ) masses. Scenarios with m W R > m ν R > Λ χ can straightforwardly be included in our framework by slightly modifying the matching discussed above.
At the scale µ = m W R we now match to an EFT in which we have integrated out the W R , ϕ H , δ ++ R and ∆ L fields, but not ν R . The resulting EFT contains apart from the SM-EFT Lagrangian, additional gauge-invariant operators that involve the right-handed neutrinos which are now relevant low-energy degrees of freedom. This induces several additional operators with respect to those in the SM-EFT. Those relevant to our analysis are [67] L ⊃ē R γ µ ν R C The matching at µ = m W R is then given by Leudϕ = 0, C Apart from the contributions to C L,R , this is equivalent to the matching in Eqs. (45), (46), and (47), without the contributions induced by the right-handed neutrinos. Note that the C There are now two cases we can consider, namely, m W < m ν R < m W R and Λ χ < m ν R < m W . Starting with the former case, one can use the RGEs in Eq. (15) to evolve the C (9) eeud coefficient from m W R to m ν R , at which point we move from a theory involving right-handed neutrinos to the SM-EFT. The corresponding matching equations are given by where m + ν R and m − ν R indicate scales just above and below the m ν R threshold. Below this threshold one can employ the RGEs in Eq. (15) to evolve the couplings to the electroweak scale, where the matching conditions in Eqs. (49), (50), and (51) still apply, and use the RGEs discussed in Section 2 to evolve the dim-9 operators to Λ χ .
In the case that m ν R lies below the electroweak scale the matching is again slightly different. One now evolves the couplings in Eq. (53) directly to the electroweak scale. Here, the same matching conditions as in Eqs. (49), (50), and (51) apply. One can then use the RGEs discussed in Section 2 to evolve the dim-9 operators to m ν R , where the right-handed neutrinos are finally integrated out. At this scale we obtain the following matching conditions Although the matching conditions are somewhat different for the three cases, m ν R ∼ m W R , m ν R < m W R , and m ν R < m W , in practice these differences only affect the running factors associated with the contributions from the right-handed neutrinos. Since the RGEs lead to factors of O(1), the numerical impact of the differences in the matching are also O(1).
The above matching shows that it is in principle straightforward to include new light degrees of freedom, such as right-handed neutrinos, in the EFT framework. In case of the right-handed neutrinos in the left-right model (LRM) this requires the inclusion of the two new operators in Eq. (52), involving ν R , with respect to those that appear in the SM-EFT. Compared to the case of heavy right-handed neutrinos only several RG factors need to be changed. It should be stressed that the above matching does not allow one to describe even lighter fields, below Λ χ . In this case one would have to keep the new light fields as degrees of freedom in the chiral EFT. Although the construction of the chiral Lagrangian will be analogous to the case discussed here, we leave the m ν R < Λ χ scenario to future work and only discuss the m ν R > Λ χ possibility in what follows.

Discussion
Irrespective of the mass of the right-handed neutrinos, after accounting for the QCD evolution and matching, we can apply the formula for the inverse half-life: Eq. (38). Looking at the amplitudes discussed in Section 4, we infer that the dim-5 term induces the standard lightneutrino exchange contribution which is proportional to A ν . The dim-7 terms give rise to A E and A me (via C The hierarchy between the dim-5, dim-7, and dim-9 contributions can be estimated by noting that C (7) /C (5) Table 3. It is instructive to go through this analysis at various levels of increasing complexity: • A first naive estimate can be obtained by neglecting any chiral suppression (i.e. taking the largest possible contributions of the dim-5,-7,-9 couplings in Table 3). This leads to a value of (v/v R ) 2 (Λ χ /v)(1/M D ) for the ratio of dim-7 to dim-5 contributions to 0νββ decay, and [(v/v R ) 2 (Λ χ /v)(1/M D )] 2 for the ratio of dim-9 to dim-5 contributions. These ratios are both O(1) for reasonable lepton Yukawa values M D ∼ m e /v and right-handed scales v R ∼ 10 TeV. This would imply that for such relatively low right-handed scales, the dim-5, dim-7, and dim-9 contributions are all of the same order.
• However, the contributions of the induced dim-7 and dim-9 operators are chirally suppressed by up to three and two powers of χ , respectively, see Table 3. As a result, the dim-9 contributions only compete with the dim-5 terms for relatively small values of m ββ , and they generally dominate the dim-7 terms for v R in the 1 − 10 TeV range.
• The above counting holds in the limit of no mixing between left-and right-handed Wbosons, ξ = 0, and assumes m ν R ∼ m W R . In this case, the dim-7 and dim-9 operators only give rise to interactions between right-handed fields at low-energies, C VR and C 1R , whose contributions are suppressed by chiral symmetry, see Table 3. On the other hand, the terms that involve left-handed fields, C VL and C (9) 4R , which depend linearly on ξ are suppressed by fewer powers of χ and can therefore be important even for relatively small values of ξ.
• The suppression of the dim-7 contributions compared to the dim-9 terms can in principle be avoided by taking M D larger than ∼ m e /v. However, for large M D the Type-I seesaw term in Eq. (45) gives a large contribution to the neutrino masses. This implies that a cancellation, and a certain measure of fine tuning, between the Type-I and Type-II seesaw mechanisms is needed to reproduce the light neutrino masses m ν ∼ 0.1 eV.
• In the case of light m ν R , the contributions to the dim-5, -7, and -9 operators induced by the exchange of right-handed neutrinos become enhanced by a factor of m W R /m ν R . For the dim-5 (dim-7) operators this enhancement is (in part) mitigated by the fact that the dim-5 terms have to reproduce the usual expression for m ββ , which has an upper bound of m ββ < 0.1 eV for reasonable values of m lightest ν < 0.1 eV. As a result, the enhancement of the type-I seesaw contribution has to be compensated; either by a smaller Dirac mass matrix or by a cancellation with the type-II seesaw term. In the case of the C-symmetric LRM it is the Dirac mass matrix that compensates, see Eq. (57), which now roughly scales as M D ∼ √ m ν R . This leads to dim-5 terms that reproduce the usual expression for m ββ , while the dim-7 terms scale as C (7) ∼ M D M ν R ∼ m −1/2 ν R and therefore are enhanced by a factor of m W R /m ν R compared to the heavy right-handed neutrino case. In contrast, the elements of the M R matrix become smaller since they scale like m ν R /m W R . The same holds for M L in the case of a C (or P) symmetry. These couplings appear in the δ ++ R contributions to the dim-9 operators, see Eq. (47), as well as in the type-II seesaw contributions to the dim-5 operator in Eq. (45). These two types of contributions are therefore suppressed by a factor of m ν R /m W R , compared to the heavy right-handed neutrino scenario. And to reiterate, dim-9 operators are enhanced by a factor of m W R /m ν R .
To see whether the above expectations hold up and to assess the implications of the TeV-scale mLRSM on 0νββ, we discuss the inverse half-life for representative regions of parameter space in the next subsection.

Phenomenology
In order to make the number of free parameters manageable we will assume a certain flavor structure for M ν R . In particular, we assume that the mixing matrix of the right-handed neutrinos is equal to that of the left-handed neutrinos. We diagonalize the left-and right-handed neutrino mass matrices as follows The assumption U = U PMNS combined with the charge-conjugation symmetry, allows us to express the Dirac mass matrix as [68] v The inverse half-life can then be expressed in terms of the mixing angles and masses of the light neutrinos and several model parameters To show the impact of the higher-dimension operators on 0νββ within the mLRSM, we plot the effective parameter as a function of the lightest neutrino mass in Fig. 3. Here we picked the following values for the model parameters These values correspond to m W R 4.5 TeV, see Appendix E, and are consistent with direct collider searches for W R bosons [69,70], heavy Majorana neutrinos [71,72], and doubly-charged scalars [73,74]. In addition, we take the central values for the mixing angles in the PMNS matrix [47], and marginalize over the Majorana phases as well as θ L and α. Finally, to show the impact of the parameter ξ, we show results with two values, namely, ξ = 0 and ξ = m b /m t . The latter value is inspired by LR models with a P symmetry which requires ξ sin α m b /m t , in order to reproduce the hierarchy between the top and bottom masses [65].
The upper panels in Fig. 3 show that the dim-7 and dim-9 contributions are subdominant for the inverted hierarchy, while they can have a significant impact in the normal hierarchy. By comparing the solid (ξ = m b /m t ) and dashed lines (ξ = 0), it is confirmed that the contributions proportional to ξ can be dominant even for modest values of ξ due to the chiral suppression of the terms independent of ξ.
The lower panels of Fig. 3 show the value of r ≡ |m eff ββ /m ββ | − 1. The blue (red) area corresponds to r if only dim-5 and dim-7 (dim-9) contributions to m eff ββ are considered. The bottom-left panel illustrates that, for the chosen values of the model parameters, the dim-9 contributions dominate over the dim-5 and dim-7 terms for light neutrino masses in the normal hierarchy. This confirms the power-counting expectations from the previous section. In the inverted hierarchy the dim-9 contributions are also larger than the dim-7 contributions, but both play a marginal role.
To study a case where the dim-7 contributions dominate over the dim-9 terms, we perform a similar analysis with the following values for the model parameters v L = 100 eV , v R = 10 TeV , V ud R = V ud L , m ∆ R = 10 TeV , m ν R 1 = 10 TeV , m ν R 2 = 12 TeV , m ν R 3 = 13 TeV . The main difference with Eq. (59) is the larger value of v L . Such large values require a significant cancellation between the type-I and type-II seesaw mechanisms in order to keep the neutrinos light enough, making these regions of parameter space less attractive. Nevertheless, we show the resulting plots in Fig. 4 and find that the dim-7 contributions have significant impact in the normal hierarchy. Finally, we investigate the impact of light right-handed neutrinos. As an example we take The difference with respect to Eq. (59) is the mass scale of the right-handed neutrinos which we now take to be O(10) GeV. The results are shown in Fig. 5. Clearly, in the case of light right-handed neutrino masses the impact of higher-dimensional contributions on m eff ββ can be significantly enhanced with respect to the case of heavy m ν R . In fact, in this case the dim-9 A similarly large effect for light right-handed neutrinos was found in Ref. [61]. That work however assumed a type-II mechanism for generating neutrino masses, which together with the assumed charge-conjugation invariance, implies M L = M † R . The right-and left-handed neutrinos are therefore diagonalized by the same matrix and, up to an overall scaling, have the same spectrum. In this situation, as the lightest active neutrino becomes lighter, the lightest right-handed neutrino mass decreases, leading to a large effect in both the normal and inverted hierarchies. By contrast, in our scenario type-I seesaw dominates the neutrino mass mechanism, and the right-handed neutrino masses remain fixed as the lightest active neutrino decreases in mass. Our results for this parameter set generalize the conclusions of Ref. [61], showing that even for type-I neutrino mass mechanisms the effects of the left-right-symmetric model on 0νββ can be significant. At the level of the amplitude, these new effects can easily dominate over the active Majorana mass contribution by a factor of 10 or 100, in the inverted or normal hierarchies, respectively.
To our knowledge such light right-handed neutrinos are not obviously excluded by other considerations. At the LHC, the right-handed neutrino can be produced with an electron, through an off-shell W * R . This leads to different signatures depending on whether the righthanded neutrino decays inside the detector or not. Assuming first that it escapes the detector, the signature is e + X + M ET where M ET refers to the missing transverse momentum. It is well known that here the signal produces an edge in the transverse mass m T distribution at high values for m T , and this can be used to obtain strong constraints on the masses of such resonances or the size of the effective higher dimension operator when below threshold. Here we update the analysis of Ref. [67], and use the m T search to place bounds on the Wilson coefficient C (6) R of the effective Lagrangian given by Eq. (52). The CMS search at 13 TeV with 35.9 fb −1 [75] considers various values for cutting on m T . We find that in the electron channel, a transverse mass cut of m T > 2 TeV gives the strongest bound, where n = 2 events are observed and 5 SM backgrounds expected. Using Bayesian statistics, this leads to a 90% credibility level upper bound of 3 signal events. Assuming a total acceptance of around 90%, this gives the bound |˜ R | < 6.5 × 10 −4 (62) in the notation of [67], where C A similar limit can be obtained from an ATLAS search [76] using 79.8 fb −1 , where no events are seen above m T > 3 TeV. We note in passing that this new LHC bound is roughly a factor of 8 stronger than that previously obtained using 7 TeV LHC data, and by a factor of 100 from neutron decay [67]. In our case, the left-right model gives C TeV. This is just below the current LHC limit, but recall we assumed all produced ν R 's escape the detector.
The right-handed neutrinos are unstable as they can decay through off-shell right-handed gauge boson exchange, or by off-shell W and Z exchange through mixing with the active neutrinos. To simplify the discussion, we consider only decays mediated by the W R . Then by comparing to muon decay, we find cτ ∼ O(cm)(10 GeV/m ν R ) 5 (m W R /4.5 TeV) 4 (9/N ef f ), where N ef f counts the number of open decay channels including color factors. Thus the decay of ν R , if it occurs in the detector, produces displaced vertices comprising either two jets and a charged lepton or neutrinos, or three leptons [77,78]. This is a complicated signal and we will not attempt to obtain any rigorous bounds, as this would require an involved simulation of the signal acceptance and recasting of existing displaced vertex searches, which is well beyond the scope of this paper 6 . Here we survey several displaced vertex searches and find none of them out-right exclude light right-handed neutrinos with masses of O(10 GeV), though further investigation is justified.
CMS has a roughly model-independent search for directly produced long-lived neutral particles decaying to displaced jets, using 2.6 fb −1 of data [81]. The CMS analysis considers pp → N N where N is long-lived and its decays produce displaced jets. We find this search to not be constraining as we now briefly discuss. In our case, we have two dominant production modes. For pp → Z * R → ν R ν R → displaced jets, this production channel directly maps onto the topology of the CMS search. However, the obtained limit for m N = 50 GeV is 1000 fb and not constraining. The other dominant production process is pp → W * R → e R ν R → e R + displaced jets/leptons, but here we cannot directly translate the CMS limit, since the kinematics are not exactly the same, though we expect it not to deviate too strongly from the previous limit. CMS also has a search for displaced jets with a strong limit of σ · BR 2 O(0.1 fb) for cτ 0 ∼ O(cm) and 38.5 fb −1 [82]. This limit is sensitive to the trigger efficiency, as they require H T > 1000 GeV. To crudely estimate H T , we observe that it is correlated with m T , and we find that for m T > 1000 GeV and m W R = 4.5 TeV, σ(pp → ν R ν R ) ∼ O(0.5 fb) which is close to the observed limit. However, the CMS analysis additionally requires 4 jets each with p T > 40 GeV. Right-handed neutrinos having masses O(10 GeV) will be highly boosted, which further reduces the signal efficiency, for the two displaced jets produced in each ν R decay will tend to merge into a single jet. While without a simulation we cannot say what the acceptance is, we expect a further suppression of at least 10%. This would make this CMS search possibly constraining.
It would be interesting to further explore LHC limits on light right-handed neutrinos that lead to displaced vertices.

Discussion and conclusion
In many well-motivated scenarios of physics beyond the Standard Model, lepton number is violated at an energy scale, Λ, well above the electroweak scale 7 . Yet, the probes of LNV with the broadest (in terms of mechanisms) and strongest sensitivity are searches for neutrinoless double beta decay, which are associated with typical nuclear scales. This large scale separation suggests that the phenomenology of 0νββ is best tackled by EFT methods, describing in a systematically improvable way the LNV dynamics both at high energy and at hadronic and nuclear scales. We stress here that an EFT approach, in conjunction with improved many-body methods, is the only path towards reaching controlled uncertainties in 0νββ calculations. This end-to-end EFT framework has been developed in several recent papers, namely Refs. [26,27] for the dimension-five Weinberg operator and Ref. [25] for dimension-seven LNV operators. The current paper summarizes and finalizes this effort by including the effects of LNV operators of dimension nine and correcting several omission in the previous papers with respect to shortdistance 0νββ mechanisms.
Our work distinguishes itself from the previous literature on 0νββ through the "end-to-end" EFT treatment, starting from the high-scale SM-EFT all the way to chiral EFT for the hadronic and nuclear aspects of the problem. For similar approaches to dimension-nine LNV operators we point to Refs. [21,24]. Other approaches often use EFT methods only to classify high-scale LNV operators and at the hadronic scale employ models and approximations that in certain cases lead to the wrong scaling of the 0νββ amplitude. In Appendix F we provide a comparison with previous approaches, highlighting the discrepancies with the EFT approach.
We now summarize the EFT framework and the main results of this paper: • Assuming a high-scale origin, the effects of LNV dynamics at scales below Λ can be expressed in an expansion in 1/Λ by considering gauge-invariant LNV operators of various dimension. After electroweak symmetry breaking via the Higgs vacuum expectation value, v, and integrating out heavy SM fields, such as W-bosons and Higgs fields, the operators can be further expanded in powers of 1/v. Finally, at energies where QCD becomes nonperturbative the various effective LNV operators can be matched to chiral EFT, the lowenergy EFT of QCD, resulting in an expansion in m π /Λ χ where Λ χ ∼ 1 GeV is the chiralsymmetry-breaking scale. In total, the 0νββ rate is expanded as where the exponents α, β, and γ depend on the LNV source and the required accuracy of the final expression. In this work, we have identified the 0νββ rate for the gauge-invariant dim-5 (α = 1), dim-7 (α = 3), and a subset of dim-9 (α = 5) operators up to leading order in the (Λ χ /v) and (m π /Λ χ ) expansions.
• The main outcome of this work is the master formula in Eq. (38) which is graphically depicted in Fig. 1. This formula allows one to calculate the 0νββ decay rate for 0 + → 0 + transitions in various isotopes as a function of effective ∆L = 2 Wilson coefficients. Hadronic, nuclear, and atomic information is captured by LECs, NMEs, and phase space factors given in Tables 1, 2, and 4, respectively. Perturbative corrections due to QCD renormalizationgroup evolution are discussed in Section 2 and Appendix D. Surprisingly we find that at leading order in the chiral expansion, all nuclear matrix elements that are necessary to describe 0νββ arising from higher-dimensional ∆L = 2 operators already appear for 0νββ induced by light Majorana neutrino exchange (corresponding to the ∆L = 2 dim-5 operator) and heavy Majorana neutrino exchange. That is, all matrix elements (apart from M AA T ) can be lifted from the existing literature and we strongly encourage future many-body calculations of light and heavy Majorana neutrino exchange to be organized in terms of the NMEs in Table 2. While different calculations vary by factors of 2-3, recent and future progress in nuclear many-body theory will allow to reduce these uncertainties [49][50][51][52][53][84][85][86][87][88][89][90][91][92].
• While all the nuclear input can be lifted from the literature, the same cannot unfortunately be said for the LECs connecting ∆L = 2 operators at the quark-gluon level to the hadronic level. Many of the LECs associated with the dim-9 operators are unknown and here we estimated them by naive dimensional analysis, or through arguments based on renormalization, leading to significant uncertainties. The number of unknown LECs grows once non-perturbative renormalization due to the nucleon-nucleon interaction of the LNV nucleon-nucleon operators is considered [27]. In particular, 4 LECs -namely g N N 2,3,4,5 -that are induced by O 2,3,4,5 and O 2,3 and are O(1) in Weinberg's counting are promoted to O((4π) 2 ) after non-perturbative renormalization of the nucleon-nucleon operators. This enhancement implies that for these quark-level operators, the induced nucleon-nucleon operators contribute to the transition amplitude at the same order as the pion-range contributions -i.e., at leading order. However, little is known about the exact values of the nucleon-nucleon LECs and in our numerical analyses we have set these effects to zero by hand, though we expect that they can affect the amplitude at the O(1) level. Lattice QCD calculations of the nn → ppee transition are direly needed to improve this situation. We have summarized the state-of-the-art values of the LECs in Table 1.
• An advantage of our EFT approach is that our work can be extended to higher orders in the various expansions. The largest corrections are due to higher-order chiral corrections that are suppressed by powers of m π /Λ χ . Such corrections were calculated for the dim-5 operator in Ref. [26] 8 , and up to next-to-next-to-leading order (O( 2 χ )) include (i) factorizable corrections, encoded by form factors in the nucleon currents; (ii) corrections due to genuinely new two-body operators, involving additional (unknown) LECs; and (iii) socalled "closure'" corrections arising from the exchange of ultrasoft neutrinos that depend on the nuclear excited states. Similar higher-order corrections can be calculated for the higher-dimensional LNV operators.
• Another advantage of the EFT approach is that it can be matched to any specific UVcomplete model of LNV. The power counting (see Table 3) allows one to isolate the important contributions and 0νββ decay rates can be immediately expressed in terms of model parameters. This allows for easy derivation of bounds on the specific high-energy model. For instance, our bounds on operators in the simplified circumstance that a single operator dominates the decay rate are given in Table 5. As a more realistic example, we have studied the minimal left-right-symmetric model in which various operators of different dimension contribute to 0νββ, with the results of our analysis shown in Figs. 3 and 4 for particular choices of model parameters. We showed in Section 6.3 and illustrated in Fig. 5 how to easily adapt our effective theory framework to include states with masses between the GeV scale and the scale of LNV physics.
• Our work essentially provides a connection between high-scale sources of LNV and lowscale experiments. As such, the framework can be used in future work to study the 0νββ phenomenology in other models for neutrino masses, as well as to compare 0νββ and collider processes such as pp → eejj at the LHC as probes of ∆L = 2 dynamics. Other extensions are the connection between models of leptogenesis and 0νββ data [93]. Finally, our framework is not applicable for non-SM states with masses below the GeV scale that occur, for example, in models with additional light sterile neutrinos. It would be interesting to extend the derivation of the 0νββ potential to include such states.

A Phase space factors and nuclear matrix elements A.1 Phase space factors
The definitions of the phase space factors appearing in Eq. (38) are given by, Here θ is the angle between the momenta of the outgoing electrons and we followed the standard normalization of Ref. [55]. The Fermi functions F (Z, E 1,2 ) take into account the fact that the outgoing electrons interact with the Coulomb potential of the daughter nucleus, and their wavefunctions are not plane waves. Their expressions are given by where R A = 1.2 A 1/3 fm and Z are, respectively, the radius and atomic number of the daughter nucleus. The Fermi functions describe the Coulomb corrections in the assumption of a uniform charge distribution in the nucleus and only account for the lowest-order terms in an expansion in the electron position. It is possible to go beyond these approximations by using exact Dirac wave functions [94] and including the effect of electron screening [95]. The use of exact wave functions leads to smaller phase space factors, with a reduction of up to 30% for the heaviest nuclei. The effects of electron screening are at the percent level [94]. The phase space factors in Table 4 do not rely on Eq. (63), but reflect the more accurate results of Refs. [53,94]. Eq. (63) can be used to get a quick estimate of the half-life, and of the differential decay rates we discuss below.
The b 0k factors are obtained from the electron traces that result from taking the square of Eq. (28). Here we follow the notation of Ref. [25], in which these factors take the following form These definitions agree with those commonly used in the literature [55], up to the trivial rescalings discussed in Ref. [25]. With the definitions of Eq. (65), the different phase space factors for a given isotope are all of similar size, with no parametric enhancement or suppression, such that the relative importance of different contributions is determined by the matching coefficients and by the nuclear matrix elements. We list the phase space factors for 76 Ge, 82 Se, 130 Te and 136 Xe in Table 4, for which we use the calculation of Ref. [53].
As discussed in Ref. [25], the measurement of the half-life in one or several isotopes will not by itself allow to disentangle the effects of dim-5, dim-7 or dim-9 operators. Some additional information can in principle be extracted from the differential decay rate with respect to the energy difference of the two electrons, y = (E 1 − E 2 )/Q, and the angle between the electron momenta, cos θ. The differential version of the master formula (38) is where all the dependence on y and cos θ is encoded in the unintegrated phase space factors g 0k The variable y = (E 1 − E 2 )/Q ∈ [−1, 1], and the dimensionless factorsb 0k are related to Eq. (65) byb 0k (y, cos θ) = 4b 0k /Q 2 , and E 1,2 = 1±y 2 Q + m e . As discussed in Ref. [25], g 02 is the phase space factor whose y dependence is distinct from the standard light neutrino exchange. A measurement of the electrons energy difference could therefore be used to single out operators, like C (6) VR , which mainly contribute to A E . Furthermore, the phase space factors in Eq. (65) exhibit a dependence on cos θ which is at most linear. The slope of the cos θ dependence of the decay rate could distinguish between the standard light neutrino exchange or the contributions of dim-9 scalar operators, for which one expects a negative slope, and C (6) VR,VL or dim-9 vector operators, which should produce a positive slope.

A.2 Nuclear Matrix Elements
To describe the nuclear parts of the amplitude, we follow standard conventions, e.g. those of Ref. [49], and introduce the following definitions where K ∈ {F, GT, T }, while j λ (|q|r) are spherical Bessel functions, with λ = 0 for F and GT, and λ = 2 for the tensor. The h ij K (r) functions describe long-range contributions, while the h ij K,sd (r) indicate short-range contributions. The factors of R A and m π have been inserted so that the h ij K,(sd) are dimensionless. The h ij K (q 2 ) are defined as follows where h AP T (q 2 ) = −h AP GT (q 2 ), h P P T (q 2 ) = −h P P GT (q 2 ), and h M M T (q 2 ) = h M M GT (q 2 )/2. In addition, at LO in χPT, g V (q 2 ) = 1, g A (q 2 ) = g A 1.27, g M (q 2 ) = 1 + κ 1 4.7, and g P (q 2 ) = −g A 2m N q 2 +m 2 π . The NMEs computed in the literature [49][50][51][52] adopt the dipole parameterization of the vector and axial form factors with vector and axial masses Λ V = 850 MeV and Λ A = 1040 MeV. The magnetic and induced pseudoscalar form factors are then assumed to be given by Using the above definitions, we express the NMEs as where the position-space tensor is defined by S (mn) (r) = 3 σ (m) ·r σ (n) ·r − σ (m) · σ (n) . The matrix elements defined in Eq. (72)  The NME of the potential V 9 (q 2 ) can be expressed in terms of the h ij K defined in Eq. (69). At LO in chiral EFT we can write Eq. (73) differs from Eq. (24) only by the momentum dependence of the axial and vector form factors g A,V (q 2 ), which is a subleading effect in χPT. The blue ellipse denotes iterations of the Yukawa potential V π , defined in Eq. (74).
are now well determined [34][35][36]. However, the ππee couplings induce ultraviolet divergences in nn → ppee scattering amplitudes and, consequently, in the 0νββ nuclear matrix elements. These divergences are analogous to those induced by the exchange of a light-Majorana neutrino, discussed in Ref. [27], and can be absorbed by promoting the NN counterterms to LO. This can be seen by repeating the analysis of Ref. [27] in the presence of dim-9 LNV operators. At LO in chiral EFT the strong interaction potential in the 1 S 0 channel has a short-range and a Yukawa component C is the LO NN coupling in the 1 S 0 channel,C ∼ {1/F 2 π , m 2 π /F 4 π }. The LNV nn → ppee scattering amplitude in the presence of dim-9 operators is obtained by sandwiching the ∆L = 2 potential V 9 (q), defined in Eq. (24), between the incoming and outgoing scattering wavefunctions, ψ ± p (r), which are obtained by solving the Schrödinger equation with the potential V 0 . The ∆L = 2 amplitude is thus given by where V 9 (r) is the Fourier transform of the potential in Eq. (24). In the case of the operators O 2,3,4,5 , in Weinberg's power counting V 9 reduces to the pionexchange contribution where the dots include the tensor component S (12) , which does not contribute to scattering in the 1 S 0 channel. This potential leads to a matrix element in Eq. (75) that is regulator dependent. Because of the short-range componentC, the scattering wavefunctions for the potential V 0 go as 1/r for r → 0, implying that the integral in Eq. (76) is logarithmically divergent. In momentum space, the UV sensitivity arises from diagrams such as the first in Fig. 6, which, when the blue ellipse is replaced by free nucleon propagators, are UV divergent. By computing the two-loop diagrams in dimensional regularization, we find an RGE that links the ππ and N N couplings induced by O 2,3,4,5 where g N N i = (m NC /(4π)) 2gN N i . Since the ππ couplings on the r.h.s. of these equations scale as Λ 2 χ , the non-perturbative renormalization of the scattering amplitude implies g N N i ∼ Λ 2 χ /F 2 π ∼ (4π) 2 . Using regulators that are more suitable for few-body nuclear physics calculations does not change this conclusion [27].
The ππ and πN couplings induced by O 1 and the πN couplings that arise from the dim-9 vector operators also cause divergences in the scattering amplitude that need to be absorbed by a contact operator. Defining In this case, however, the scaling of the nucleon-nucleon operators is not affected since both sides of Eq. (79) scale as F 2 π ∼ m 2 π , so that the size of g N N 1,6,7,8,9 agrees with Weinberg's power counting.

C ∆L = 2 potentials induced by dim-3, -6 and -7 operators
The chiral Lagrangians and the 0νββ potentials induced by the light neutrino Majorana mass and by dim-6 and -7 operators are given in Refs. [26,27] and [25], respectively. In this appendix we recall the main ingredients, and include a few missing operators.
The neutrino Majorana mass in Eq. (2) induces 0νββ by coupling to nucleons and pions via the vector and axial weak currents. The pion and single nucleon currents can be straightforwardly built in χPT [40,41], and they give rise to the long-range component of the neutrino potential, as reviewed in Refs. [25,26]. The exchange of hard neutrinos induces short-range contributions to the 0νββ potential, parametrized by the operator g N N The renormalizability of nn → ppee scattering amplitude requires the LEC g N N ν to scale as g N N ν ∼ 1/F 2 π [27], implying that g N N ν contributes to the neutrino potential at LO. The neutrino potential in momentum space from light Majorana-neutrino exchange is The LO neutrino potential is obtained by setting the single nucleon form factors in h F,GT,T to their LO values. As discussed in Ref. [25,26], including the q 2 dependence of the vector and axial form factors accounts for a subset of the N 2 LO corrections. g N N ν is at the moment unknown, but can be obtained by matching chiral EFT and lattice QCD calculations of the nn → ppee amplitude performed at the same kinematic point. The value of g N N ν thus depends on the details of the strong interaction potential, and its use in many-body calculations requires a consistent treatment of short-range effects in the strong and weak sector.
The dim-6 and dim-7 operators in Eqs. (4) and (5) contain neutrinos, which need to be exchanged between nucleon lines to cause 0νββ. In Ref. [25] we considered the long-range ∆L = 2 potentials induced by these operators, which are mediated by the pion and nucleon vector, axial, scalar, pseudoscalar, and tensor currents.
The 0νββ transition operators arising from dim-6 and dim-7 can be divided in four components V a is the largest piece and is dominated by the contribution of the pseudoscalar operators C Note that this potential goes like 1/(q 2 ) 2 for large q 2 , ensuring that this potential does not induce divergences similar to those produced by V 9 as discussed in Sect. B. The long-range neutrino exchange contributions of the tensor operator, C T , and the rightand left-handed currents operators, C VR and C (6) VL , are chirally suppressed. For this reason, short-range effects, arising from the exchange of neutrinos with momentum ∼ Λ χ , could become important, even in Weinberg's power counting. These effects can be accounted for by building hadronic operators with the same transformation properties as the tensor products of one nonstandard current and the SM weak interaction. Carrying out the lepton contractions while neglecting the lepton momenta, we find that the operators induced by C (6) T,VL,VR transform like the following non-local terms VL : VR : where L µνρσ = g µρ g νσ + g νρ g µσ − g µν g ρσ + iε µρνσ . From Eq. (84) we see that C T induces operators that transform as a Lorentz scalar, but have the same chiral properties as the vector operators O 6,7,8,9 in Eq. (13). C (6) VL generates operators that are Lorentz vectors, but with the same chiral properties as O 1 . Similarly, short-range operators induced by C (6) VR are Lorentz vectors with the same chiral properties as O 4 , but with negative parity.
The form of short-range operators follows from the Lorentz, chiral, and parity transformation properties of the operators in Eq. (84). We construct VR , on the other hand, induces short-range operators that are parityodd, if one neglects the lepton momenta, and thus does not contribute to Eq. (85) or 0 + → 0 + transitions.
From the discussion in Ref. [25] and Eq. (85), it follows that the 0νββ potential induced by the tensor operator C (12) .
The contributions from g N N, πN T and g ππ T were neglected in Ref. [25]. The leading potential induced by C The short-range terms proportional to g N N VL and g πN VL are formally of the same order as first contribution in Eq. (87), mediated by the nucleon magnetic moment.
Finally, C VL,VR induce contributions proportional to the lepton momenta. In the case of C (6) VL , these contributions are suppressed by χ with respect to Eq. (87), while they are the leading-order contribution in the case of C (6) VR . Also in this case additional contact interactions need to be introduced in order to make the amplitude regulator independent. These short-range interactions take the form, where g E,me V L,V R are LECs of O(1). The combination of these contributions and the long-range terms give, where M (1) L,R = 1 3

D Renormalization group equations
In this appendix we briefly discuss the running between the electroweak scale and Λ χ , which follows from the RGEs discussed in section 2. The solutions of the RGEs for the dimension-six operators are where η ≡ αs(m W ) αs(µ) and β 0 = 1 3 (11N c − 2n f ), with n f the number of active flavors. From the RGEs of the dimension-nine operators we obtain with C F = (N 2 c − 1)/(2N c ), while the solutions to the RGEs for C 2,3 and C 6,7 are more conveniently expressed as where γ (i) 1,2 are the eigenvalues of the anomalous-dimension matrices and R i the matrices that diagonalize them. They can be written as whereγ = 4N 2 c − 11 + 16/N 2 c . As discussed in Section 2, the remaining couplings either do not run, or have an equivalent evolution as the above couplings.

E The left-right model
In this appendix we discuss a few more details of the mLRSM that are needed to obtain the matching conditions in Eq. (44). The main ingredients that are missing from Section 6 are the masses of the heavy BSM fields. After the right-handed triplet field, ∆ R , acquires a vev, several fields obtain masses proportional to v R and therefore become heavy. These fields are the righthanded gauge bosons, W R and Z R , the right-handed neutrinos, ν R , and the left-handed triplet fields, ∆ L . In addition, part of the right-handed triplet fields, namely δ ++ R and Re δ 0 R , as well as part of φ become heavy. The δ ± R and Im δ 0 R fields are 'eaten' by the W R and Z R gauge bosons, while the remaining part of φ stays light and can be interpreted as the SM Higgs doublet.
To integrate out these heavy fields we first have to rotate to the mass basis. To do so we write φ in terms of two SU (2) L doublets, φ ≡ (φ 1 , φ 2 ), and use the relation to the mass eigenstates φ 1 which follows from the Higgs potential, see e.g. Refs. [96][97][98][99]. Hereφ i = iτ 2 φ * i and ϕ is the SM Higgs doublet, while ϕ H obtains a mass of O(v R ). The contributions to ∆L = 2 interactions arise from exchanges of ν R , W R , δ ++ R , and ∆ L , so that the relevant masses are the following where g R is the SU (2) R gauge coupling, U is the rotation matrix between the mass and flavor bases, ν (mass) R = U ν R , and ρ 1,2,3 are parameters in the Higgs potential [98]. After integrating out these heavy fields, the induced ∆L = 2 interactions will depend on the masses of these fields. In our analysis we choose values for m ∆ R , v R (which determines m W R ), and m ν R = diag(m ν R 1 , m ν R 2 , m ν R 3 ).
Similarly, integrating out ∆ L gives a contribution to C (5) which depends on its mass. Explicitly, we have ξβ 1 e i(δ β 1 +α) + β 2 e iδ β 2 + ξ 2 β 3 e i(δ β 3 +2α) M L , where β i and δ β i are parameters of the Higgs potential in the notation of Ref. [99]. However, we do not have to choose values for β i , δ β i , and m ∆ L as we can trade the combination of these parameters for the vev of ∆ L , v L . The reason for this is that the β i couplings represent the terms in the Higgs potential that are linear in ∆ L . The same terms also appear in the minimum condition, ∂V H ∂∆ L S i = S i = 0, where S i stand for the scalars in the mLRSM. The only other (nonnegligible) terms appearing in this condition are the terms quadratic in ∆ L , which determine the mass term, since m 2 . This leads to the following relation v L e iθ L v 2 = − v R 2(1 + ξ 2 ) ξβ 1 e i(δ β 1 +α) + β 2 e iδ β 2 + ξ 2 β 3 e i(δ β 3 +2α) m 2 where we used v 2 = κ 2 (1 + ξ 2 ), and the right-hand side features the same combination that appear in Eq. (97). After substituting the above relation in Eq. (97), one obtains the second term in Eq. (45). Finally, as mentioned, we assume a charge-conjugation symmetry which acts on the fields as This symmetry enforces a relation between the gauge couplings and Majorana mass matrices, g R = g L and M L = M † R , ensures the Dirac mass matrices are symmetric, and leads to the relation in Eq. (57).

F Comparison with the literature
The dim-9 operators of Eq. (7) have been discussed in the context of 0νββ in the literature before [21,24,31,58,100]. Since not all references apply the same framework based on χPT that is employed here, there are several instances were our results are significantly different with respect to older findings. Here we give a brief overview of the most significant discrepancies and focus on the comparison with Refs. [31,58].
The dim-9 SU (3) c × U (1) em -invariant four-quark two-lepton operators were first cataloged in Ref. [31], while Refs. [21,24] removed several redundancies. To perform a quantitative comparison, we note that the Wilson coefficients of Ref. [31] are related to the ones in Eq. (7) by The operators with subscript R are obtained by sending L to R in the last index of the coefficients in Eq. (100). We start by noting that, after using the translation of Eq. (100), our results for the RGEs of the dimension-nine operators do not agree with parts of the literature. In particular, the RGEs for C (9) 2,3 in Eq. (15) as well as those for the vector operators, C (9) 6,7,8,9 in Eq. (16) differ from those derived in Ref. [58]. The RGEs in Eq. (15) are in agreement with Refs. [38,39], while the vector operators were not considered in those references. However, since the running gives rise to O(1) effects, as illustrated by Table 5, these discrepancies have a relatively mild effect on resulting constraints.
More significant discrepancies arise from the different approaches to the hadronization of the quark-level operators. In this work we perform the matching of the quark-level operators to interactions in the chiral Lagrangian, in which the non-perturbative nature of QCD is captured by the LECs in Table 1. Instead, Ref. [31] relies on factorization to estimate the required matrix elements and neglects the effects of ππ and πN operators, effectively matching onto N N interactions only. This approach gives comparable results to the ones obtained here for several couplings, namely C (9) 1, 6,7,8,9 . The reason for this is that the corresponding operators induce N N interactions at leading order, which do not require enhancement with respect to Weinberg's counting (i.e. they follow NDA). However, there are several operators, namely C (9) 2,3,4,5 , for which the leading contributions come from ππ interactions and the N N terms are enhanced compared to Weinberg's counting. Both of these effects are missed by the approach of Ref. [31] and lead to underestimates of the contributions to the amplitude of 0νββ by a factor of O(16π 2 ) compared to the results found here. These significant differences can be made explicit by comparing the different estimates for the LECs. As mentioned Ref. [31] does not include πN and ππ interactions, effectively setting the corresponding LECs to zero. The factorization estimates for the N N interactions give in our notation where g S and g T are the isovector scalar and tensor densities. Ref. [31] employs the MIT-bag model estimates g S = 0.48 and g T = 1.38 of Ref. [101], for which lattice QCD determinations are available, see Table 1.
Note that the values used in Ref. [31], together with Eq. (101), imply O(1) values, or slightly smaller, for the g N N i . As mentioned above, this leads to similar, though slightly weaker, limits for the C (9) 1,6,7,8,9 operators. Instead, our estimates for g N N 2,3,4,5 are a factor O(16π 2 ) larger than those used in Ref. [31]. In addition, the ππ interactions, for which the LECs are known, give rise to contributions at the same order in these cases. As a result, we find O(10−100) stronger limits on C (9) 2,4,5 compared to Ref. [31]. These different values for the LECs also led to the conclusion in Ref. [58] that the limits on the scalar and tensor couplings LLR 1,2 , equivalent to our C 2L,3L , are significantly affected by the RGEs. The reason for this is that Refs. [31,58] employ a matrix element for LLR 1 which is smaller than that for LLR 2 by a factor of O(100). This implies that the largest contribution from LLR 1 actually arises from the mixing into LLR 2 , making the running a large effect. However, since the two operators transform in the same way under the symmetries of QCD we estimate O (9) 2L,3L to have matrix elements of similar size, so that the running is not more than an O(1) effect (as can be seen from Table 5). We stress that this expectation is borne out by the lattice determinations of the ππ matrix elements in Table 1. This significant discrepancy in the limits is in part due to the fact that Ref. [31] neglects the ππ interactions. The importance of these terms was already pointed out in Refs. [100] and [21,24]. In fact, the latter two references work within the framework of χPT, and therefore obtain very similar results as we do here. The main difference between this work and Refs. [21,24] is that the latter references assumed Weinberg's power counting and therefore did not take into account the enhancement of g N N 2,3,4,5 . However, the consequences of this assumption are much less severe since the ππ interactions contribute at the same order, and the LECs g N N 2,3,4,5 are currently unknown. An additional difference is that Ref. [21] neglected the color-mixed operators O 3 , O 5 , O 7 and O 9 and thus considered an incomplete low-energy basis [24].
Finally, we stress that although the 0νββ expressions in this work depend on several unknown LECs, the χPT framework allows for systematic improvement upon the current situation. In particular, the LECs which can only be estimated through NDA and RG arguments at the moment, namely g πN i and g N N i in Table 1, are in principle amenable to lattice QCD determinations. Such lattice calculations would provide control over the hadronic uncertainties of the contributions to 0νββ and are a necessary ingredient to obtain reliable limits on the dim-9 operators.