Neutrinoless double beta decay in effective field theory: the light Majorana neutrino exchange mechanism

We present the first chiral effective theory derivation of the neutrinoless double beta-decay $nn\rightarrow pp$ potential induced by light Majorana neutrino exchange. The effective-field-theory framework has allowed us to identify and parameterize short- and long-range contributions previously missed in the literature. These contributions can not be absorbed into parameterizations of the single nucleon form factors. Starting from the quark and gluon level, we perform the matching onto chiral effective field theory and subsequently onto the nuclear potential. To derive the nuclear potential mediating neutrinoless double beta-decay, the hard, soft and potential neutrino modes must be integrated out. This is performed through next-to-next-to-leading order in the chiral power counting, in both the Weinberg and pionless schemes. At next-to-next-to-leading order, the amplitude receives additional contributions from the exchange of ultrasoft neutrinos, which can be expressed in terms of nuclear matrix elements of the weak current and excitation energies of the intermediate nucleus. These quantities also control the two-neutrino double beta-decay amplitude. Finally, we outline strategies to determine the low-energy constants that appear in the potentials, by relating them to electromagnetic couplings and/or by matching to lattice QCD calculations.

By itself, the observation of 0νββ would not immediately point to the underlying mechanism of LNV. In an effective theory approach to new physics, LNV arises from ∆L = 2 operators of odd dimension, starting at dimension-five [11][12][13][14]. As discussed in detail in Ref. [15], if the scale of lepton number violation, Λ LNV , is in the range 1-100 TeV, short-distance effects encoded in local operators of dimension seven and nine provide contributions to 0νββ within reach of next generation experiments. However, whenever Λ LNV is much higher than the electroweak scale, the only low-energy manifestation of this new physics is a Majorana mass for light neutrinos, encoded in a single gauge-invariant dimension-5 operator [11], which induces 0νββ through light Majorananeutrino exchange [16][17][18].
To interpret positive or null 0νββ results in the context of this minimal extension of the SM (the three light Majorana neutrinos paradigm), it is critical to have good control over the relevant hadronic and nuclear matrix elements. Current knowledge of these is somewhat unsatisfactory [19], as (i) few of the current calculations are based on a modern effective field theory (EFT) analysis, and (ii) various approaches lead to estimates that differ by a factor of two to three. In this paper we present the first end-to-end EFT analysis of 0νββ induced by light Majorana-neutrino exchange, describing the physics from the scale Λ LNV all the way down to the nuclear energy scale. The EFT framework has allowed us to identify long-and short-range contributions to 0νββ previously missed in the literature, that are, by power counting, as large as corrections usually included. The main results of our work are expressions for the leading and next-to-next-to-leading order (N 2 LO) chiral potentials mediating 0νββ, and the amplitude induced by the exchange of ultrasoft neutrinos, with momenta much smaller than the Fermi momentum.

II. EFFECTIVE THEORY FRAMEWORK
The starting point of our analysis is the weak scale effective Lagrangian, which we take to be the SM augmented by Weinberg's ∆L = 2 dimension-five operator [11], where u αβ is a 3 × 3 matrix, L = (ν L e L ) T is the left-handed SU (2) lepton doublet, H is the Higgs doublet, α, β ∈ e, µ, τ , and i, j, m, n are SU (2) indices. This operator induces a Majorana mass matrix for neutrinos, of the form m αβ = −u αβ (v 2 /Λ LNV ), where v = ( √ 2G F ) −1/2 246 GeV is the Higgs vacuum expectation value: for Λ LNV v this is the well known "seesaw" relation. Neglecting QED and weak neutral-current effects, the low-energy effective Lagrangian at scale µ Λ χ ∼ 1 GeV is given by The second term in (2) represents the Fermi charged-current weak interaction. The last two terms encode LNV through the neutrino Majorana mass, given by m ββ = i U 2 ei m i in terms of mass eigenstates and elements of the neutrino mixing matrix, and a dimension-nine ∆L = 2 operator generated at the electroweak threshold: O L =ē L e c Lū L γ µ d LūL γ µ d L , with e c L = Cē T L . Since , the effect of the latter term on the 0νββ amplitude is suppressed by (k F /M W ) 2 (where k F ∼ O(100) MeV is the typical Fermi momentum of nucleons in a nucleus) compared to light-neutrino exchange and can be safely neglected at this stage.
The interactions of Eq. (2) induce ∆L = 2 transitions (such as π − π − → e − e − , nn → ppe − e − , 76 Ge → 76 Se e − e − , 136 Xe → 136 Ba e − e − , ....) through the non-local effective action obtained by contracting the neutrino fields in the two weak vertices, where is the scalar massless propagator. Computing matrix elements of S ∆L=2 eff in hadronic and nuclear states is a notoriously difficult task. The multi-scale nature of the problem can be seen more explicitly by going to the Fourier representation 1 The amplitude (5) receives contributions from neutrino virtualities ranging from the weak scale all the way down to the IR scale of nuclear bound states. Roughly speaking one can identify three regions, whose contributions can be conveniently described in terms of appropriate effective theories: This contribution is controlled by the quark-level short-distance behavior of the correlator (6). An Operator Product Expansion analysis shows that integrating out hard neutrinos and gluons generates a local term in the effective action proportional to O L , with Wilson coefficient This short-distance component is currently missing in all calculations of 0νββ, which start from the nucleon-level realization of the weak currents in the correlator (6). Within such approaches, the new effect can be estimated by considering the hadronic realization of O L , sensitive to pion-range and short-range nuclear effects, that has been studied in the context of TeV-sources of LNV [20][21][22][23]. In what follows we adopt a chiral EFT approach and the effect of hard modes will be encoded in local counterterms of the low-energy effective chiral Lagrangian, transforming as O L under the chiral group. 1 To obtain (5) we have approximatedēL(x)γ µ γ ν e c L (y) ēL(x)γ µ γ ν e c L (x) = g µνē L(x)e c L (x), which amounts to neglecting the difference in electron momenta, a safe assumption given that |p1 − p2|/kF 1. 2 Within Lattice QCD, OL captures O(a 2 ) discretization effects in the calculation of the amplitude (5). OL would appear in the Symanzik's action [24,25] with a pre-factor scaling as O(αsa 2 ) near the continuum limit. Similar contributions relevant to the case of two-neutrino double beta decay (2νββ) have been discussed in Ref. [26].
(ii) A soft and potential region with Here the appropriate hadronic degrees of freedom are pions and nucleons, described by chiral EFT. In analogy with the strong and electroweak interactions in the SM, integrating out pion degrees of freedom and neutrinos with soft (k 0 ∼ |k| ∼ m π ∼ k F ) and potential (k 0 ∼ k 2 F /m N , |k| ∼ k F ) scaling of their 4-momenta generates nucleon level ∆L = ∆I = 2 potentials that mediate 0νββ between nuclear states.
(iii) Ultrasoft or "radiation" region, with neutrino momenta scaling as k 0 ∼ |k| k F . Here the effective theory contains as explicit degrees of freedom nucleons interacting via appropriate potentials (see (ii) above), electrons, and essentially massless neutrinos, whose ultrasoft modes cannot be integrated out (similarly to gauge fields in NRQED and NRQCD [27][28][29]). These modes do not resolve the nuclear constituents and this part of the amplitude is sensitive to nuclear excited states and transitions among them induced by the electroweak currents.
Contributions to 0νββ from regions (ii) and (iii) are included in all existing calculations albeit within certain approximations and not fully in the spirit of EFT. In particular, we have identified corrections that can not be parameterized through the single nucleon form factors. We next discuss the 0νββ amplitude in the context of chiral EFT, in which the contributions from region (i) are captured by local counterterms, the contributions from region (ii) can be explicitly evaluated and lead to appropriate potentials, and the contributions from region (iii) can be displayed in terms of non-perturbative nuclear matrix elements of the weak charged current and bound state energies.
Chiral symmetry and its spontaneous and explicit breaking strongly constrain the form of the interactions among nucleons and pions. In the limit of vanishing quark masses, the χPT Lagrangian is obtained by constructing all chiral-invariant interactions between nucleons and pions. Pion interactions are derivative, allowing for an expansion in p/Λ χ , where p is the typical momentum scale in a process and Λ χ ∼ m N ∼ 1 GeV is the intrinsic mass scale of QCD. One can order interactions according to the chiral index ∆ = d+n/2−2, where d counts the number of derivatives and n counts the number of nucleon fields [30,34]. Chiral symmetry is explicitly broken by the quark masses and charges, and, in our case, by electroweak and ∆L = 2 operators. However, the explicit breaking is small and can be systematically included in the power counting by considering m q ∼ m 2 π ∼ p 2 . In presence of lepton fields we generalize the definition of chiral index to ∆ = d + n/2 − 2 + n e , where n e denotes the number of charged leptons in the interaction vertex. With this definition, the lowest order 0νββ transition operators have chiral index ∆ = 0. For nuclear physics applications, one has p ∼ k F ∼ m π and the expansion parameter is χ = m π /Λ χ . For 0νββ there is an additional infrared scale. The reaction Q value, the electron energies E 1,2 , and the energy differences E n − E i of the bound nuclear states have typical size O(5 − 10) MeV, to which we assign the scaling k 2 F /m N ∼ k F χ . Our building blocks are the pion field u = exp (iπ · τ /(2F 0 )) (where F 0 is the pion decay constant in the chiral limit, and F π = 92.2 MeV) and the nucleon doublet N = (p n) T , transforming as u → LuK † (π) = K(π)uR † and N → K(π)N under the SU (2) L × SU (2) R chiral group [37,38].
The effective Lagrangian of Eq. (2) maps onto the following operators with zero chiral index, where , and g 0 A is the LO contribution to the nucleon axial coupling which is measured to be g A = 1.2723(23) [39]. Tree-level diagrams involving the above interactions and Majorana neutrino exchange generate ∆L = 2 amplitudes such as π − π − → ee and nn → ppee, scaling as O(G 2 F m ββ ). At the one-loop level UV divergences appear which require the introduction of ∆L = ∆I = 2 local operators with chiral index ∆ = 2. We find three independent structures with the correct transformation properties: Here L µ = uu µ u † , the dots stand for terms involving more than two pions, and three a priori unknown O(1) low-energy constants (LECs) appear: g ππ ν , g πN ν , and g N N ν . In the mesonic and single-nucleon sector of the theory, all momenta and energies are typically ∼ p, and the perturbative expansion of the χPT Lagrangian and power counting of loops [30] implies that the scattering amplitudes can also be expanded in p/Λ χ . For systems with two or more nucleons the energy scale p 2 /2m N becomes relevant and the corresponding amplitudes do not have a homogeneous scaling in p. Therefore, the perturbative expansion of interactions does not guarantee a perturbative expansion of the amplitudes [34,35]. Indeed, the so-called "reducible" diagrams (in which the intermediate state consists purely of propagating nucleons) are enhanced by factors of m N /p with respect to the χPT power counting and need to be resummed. On the other hand, loop diagrams whose intermediate states contain interacting nucleons and pions -"irreducible"-follow the χPT power counting [34,35]. Reducible diagrams are then obtained by patching together irreducible diagrams with intermediate states consisting of A free-nucleon propagators. This is equivalent to solving the Schrödinger equation with a potential V defined by the sum of irreducible diagrams. For external perturbations, such as electroweak currents and ∆L = 2 interactions, one can similarly identify irreducible contributions, that admit an expansion in p/Λ χ [40].
While the scaling of irreducible loop diagrams is unambiguous, the power counting for fournucleon operators has been the object of much debate in the literature [40][41][42]. In the Weinberg power counting [34,35], the scaling is determined by naive dimensional analysis, and the lowest order four-nucleon operators in the strong and ∆L = 2 sectors are given, respectively, by (9). While phenomenologically successful [43], the Weinberg power counting is not fully consistent. Inconsistencies appear in some channels, such as the 1 S 0 channel, where the cutoff dependence of the solution of the Lippman-Schwinger equation cannot be absorbed by the counterterms that appear at lowest order [41,42]. Various solutions to this problem have been proposed, including treating pion exchange in pertubation FIG. 1. Diagram contributing to the leading-order neutrino potential. Double and single lines denote, respectively, nucleons and lepton fields. The black square denotes an insertion of the neutrino Majorana mass, while the gray circle denotes the SM weak charged-current interaction. theory (perturbative pion or "KSW" scheme [41,44,45]), expanding the nuclear forces around the chiral limit [46], or, for processes at low enough energy, integrating out pions and working in the pionless EFT [40,47]. In Section III A, we will discuss 0νββ in the Weinberg power counting, and we will extend the treatment to the pionless EFT in Section III B.
Finally, matching the chiral EFT framework to many-body quantum mechanics one obtains the following nuclear Hamiltonian appropriate for calculating 0νββ amplitudes: (10) The first term (H 0 ) encodes the strong interaction. In the Weinberg counting, the leading-order strong potential is given by one pion exchange plus the contact terms C S,T [34,35], which in momentum space reads (q is the momentum conjugate to Here and in the following, we replace the LO couplings and decay constants by their physical values, g 0 A → g A , F 0 → F π , etc., which can be consistently done to the order we are working in the chiral expansion. The second term in (10) is the usual charged-current weak interaction. From now on, we set g V = 1, neglecting small isospin-breaking corrections. Note that light Majorana neutrinos and electrons with ultrasoft momenta are active degrees of freedom in the low-energy theory.
The third term in (10) directly mediates ∆L = 2 amplitudes, and we discuss it next.
A. The ∆L = ∆I = 2 potential The potential V ν encodes physics from hard scales (the counterterms of Eq. (9)) as well as soft scales, obtained by integrating out pions and Majorana neutrinos with soft and potential scaling of their 4-momenta. In practice V ν is given by the sum of "irreducible" diagrams mediating nn → ppee in chiral EFT. As discussed above, V ν admits a chiral expansion: The LO neutrino potential is obtained by tree-level neutrino exchange, which involves the single-nucleon currents (see Fig. 1). In momentum space one finds [15] Loop diagrams contributing to an effective ππe − e − vertex (upper panel), and to an effective npπe − e − vertex (lower panel). Pions are denoted by dashed lines, the remaining notation is as in Fig. 1.
The diagrams give rise to corrections to the ∆L = 2 potential when the pions are connected to external nucleon lines.
Analogously to the strong-interaction case, the neutrino potential V ν depends only on the momentum scale q ∼ k F and not on infrared scales such as the excitation energies of the intermediate odd-odd nucleus in 0νββ, often approximated by their averageĒ − 1/2(E i + E f ). Note that the commonly used neutrino potential [16,19] reduces to V ν,0 whenĒ − 1/2(E i + E f ) is set to zero. At N 2 LO in the Weinberg power counting several new contributions arise. These consist of (a) corrections to single-nucleon currents, which are often included in the literature via momentumdependent form factors; (b) genuine N 2 LO two-body effects, such as loop corrections to Fig. 1, which induce the short-range neutrino potential V ν,2 that has never been considered in the literature. Note that two-nucleon effects in the weak currents [48,49], which induce three-nucleon potentials in Eq. (12), start contributing to 0νββ at N 3 LO, once one takes into account that where we have introduced the tensor operator The functions h F , h GT , and h T are expressed in terms of the isovector vector, axial, induced pseudoscalar, and magnetic nucleon form factors as [19] h In the literature, the dipole parameterization of the vector and axial form factors is often used 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 where κ 1 = 3.7 is the nucleon isovector anomalous magnetic moment. Expanding Eqs. (16) and (17) for small |q|, one recovers the LO and, for g A (q 2 ), the N 2 LO χPT expressions of the nucleon form factors. In the case of g V , g P and g M , the N 2 LO χPT results, given for example in Ref. [50], deviate from Eqs. (16) and (17). However, any parameterization that satisfactorily describes the observed nucleon form factors can be used in the neutrino potential (14). The potential V ν,2 is induced by one-loop diagrams with a virtual neutrino and pions contributing to nn → ppee, built out of the leading interactions of Eqs. (8). They can be separated into three classes, involving the ππ → ee (Fig. 2, upper panel), n → pπ + ee (Fig. 2, lower panel), and nn → ppee (Fig. 3) effective vertices. Note that for diagrams such as Fig. 3(b) or 3(d) we include only the two-nucleon irreducible component. We regulate the loops dimensionally, with scale µ, and subtract the divergences according to the MS scheme. The UV divergences are absorbed by the counterterms of Eq. (9), which cancel the µ dependence of the loops and also provide finite contributions.
V ν can be thought of as the matching coefficient between the chiral EFT and the low-energy nuclear EFT described by Eq. (10), containing non-local potentials and ultrasoft neutrino modes. The matching is achieved by subtracting the low-energy theory diagrams depicted in Fig. 4, involving ultrasoft neutrino exchange and insertions of the LO strong and ∆L = 2 potentials, from the chiral EFT diagrams of Figs. 2 and 3. Since the two EFTs have the same IR behavior, the IR divergences, stemming from diagrams (m) of Fig. 2 and (a), (b) of Fig. 3, cancel in the matching. We have checked this by regulating the IR divergences with a neutrino mass. Moreover, ultrasoft neutrino loops in Fig. 4 contain UV divergences, which we deal with in dimensional regularization and MS subtraction, with renormalization scale µ us . Thus the matching leads to a term in the potential V (a.b) ν,2 that depends logarithmically on µ us . As we show in Section III C below, the dependence on µ us cancels once one includes the ultrasoft contribution to the 0νββ amplitude.
Finally, the counterterm potential reads The µ dependence of g ππ ν , g πN ν , and g N N ν cancels the µ dependence of L π in Eqs. (19) and (20). We will discuss strategies to estimate the finite parts of the LECs in Section III D below.

B. The ∆L = ∆I = 2 potential in the pionless EFT
The previous discussion assumed the Weinberg power counting, that, while phenomenologically successful [43], is not formally consistent [41,42]. Few-body systems and processes characterized by scales p m π can be studied in pionless EFT, a low-energy EFT in which pion degrees of freedom are integrated out [40,47]. For physical pion masses, pionless EFT converges very well for the A = 2, 3 systems, and works satisfactorily well for up to A = 6 [52]. While the application of this EFT to nuclei with A > 6 needs to be studied in more detail, it is interesting to extend the framework developed in the previous section to pionless EFT, especially in the light of a possible matching to lattice calculations of 0νββ matrix elements performed at heavy pion masses. A similar matching between lattice and pionless EFT for strong interaction and electroweak processes has been carried out in Refs. [26,[53][54][55][56][57][58]. While the lattice QCD calculations relevant for 2νββ were performed at a single lattice spacing of a ∼ 0.145 fm and at m π ∼ 806 MeV [26,58], they represent the first step for the field and are very promising. We are optimistic that in the near future lattice calculations of electroweak processes and 0νββ in the two nucleon system will reach control over all lattice systematics, as recently achieved for the nucleon axial coupling g A [59][60][61][62].
Pionless EFT describes physics at the scale p smaller than the cutoff of the theory Λ / π ∼ m π . For power counting purposes, we introduce the scale ℵ ∼ p Λ / π . The leading-order Lagrangian is given by Eq. (8c), and the fine tuning of the S-wave nucleon-nucleon scattering lengths is accounted for by assigning the coefficients C S,T the scaling Using dimensional regularization with power divergence subtraction (PDS) [45], at the scale µ the couplings C S,T can be expressed in terms of the spin-singlet 1 S 0 and spin-triplet 3 S 1 scattering lengths a s and a t according to where a −1 s ∼ −8.3 MeV, a −1 t = 36 MeV. Higher-order operators involve additional derivatives and are related to additional parameters (effective range, shape parameter, . . . ) of the effective-range expansion. Note that in pionless EFT the three-body nucleon force is a leading-order effect [40].
The leading ∆L = 2 potential in the pionless EFT has the form The first term comes from long-range neutrino exchange, as in Eq. (13), with the difference that the contributions of the induced pseudoscalar form factor are subleading. In addition, the scaling of the nucleon-nucleon coupling g N N ν , introduced in Eq. (9), is modified. This coupling connects two S-waves and thus is enhanced in the pionless theory [40], scaling as This scaling can be understood by studying the scattering amplitude for two neutrons to turn into two protons with the emission of two zero-momentum electrons. At leading order in the pionless EFT, the scattering amplitude in the 1 S 0 channel receives contributions from the diagrams in Fig. 5, where the two neutrons and two protons in the initial and final states are dressed by bubble diagrams with insertions of the leading order contact interaction C s . These contributions to the amplitude have the schematic form where T is the leading-order, strong-interaction scattering amplitude in the 1 S 0 channel, which is scale independent, and the dots in Eq. (31) denote additional scale-independent contributions. I 2 is the dimensionless two-loop integral that appears in the first diagram of Fig. 5. The loop is logarithmically divergent in d = 4, giving, in the PDS scheme, where E is the energy of the two neutrons in the initial state, and P the center-of-mass momentum. This is the same UV divergence that appears in Coulomb corrections to proton-proton scattering [63]. The amplitude (31) can be made independent of the renormalization scale µ by rescaling whereg N N ν = O(1) and satisfies d d log µg Eq. (33), together with (28), confirms the power-counting expectation of Eq. (30).
Beyond leading order in the pionless EFT there appears a single four-nucleon operator contributing to V ν,2 at N 2 LO, which is conveniently expressed in terms of the 1 S 0 projectors P √ 2, and κ given in Eq. (9). The LEC g N N ν, 2 scales as Operators connecting two neutrons and two protons in the P -waves are not enhanced by ℵ −2 , and appear at even higher order. Additional corrections to the potential V ν,0 arise from loop diagrams (3n) and (3o). These diagrams are scaleless and vanish in dimensional regularization. If the infrared divergence is regulated by a neutrino mass m ν , the m ν dependence is canceled by the diagrams in Fig. 4, and one obtains The µ dependence in (37) is reabsorbed by a sub-leading term in g N N ν . Eq. (37) shows that it is always possible to choose µ us so that the correction to the potential vanishes, and the effect of diagrams (3n) and (3o) is all encoded in the ultrasoft contribution. Note that the loop (37) and the ultrasoft amplitude are suppressed by p/(4πm N ) ∼ (p/Λ / π ) × 1/(4π) 2 with respect to the LO, and are thus smaller than corrections from Eq. (35), scaling as p 2 /Λ 2 / π . The relevance of the ultrasoft region can be also understood diagrammatically. Indeed, while diagram (3o) is suppressed with respect to the LO, diagrams with an arbitrary number of insertions of C s and C t between the emission and absorption of the neutrino are not suppressed with respect to (3o). These diagrams need to be resummed, and correspond to building up the intermediate states. They are captured by the ultrasoft contribution discussed in Section III C below.

C. The 0νββ amplitude
Starting from the nuclear Hamiltonian of Eq. (10), one calculates the full 0νββ amplitude as the sum of two contributions 3 where we defined T lept = 4G 2 F V 2 ud m ββūL (p e1 )Cū T L (p e2 ). The first term represents a single insertion of the ∆L = ∆I = 2 potential (third term in Eq. (10)). On the other hand T usoft arises from double insertions of the weak interaction (second term in Eq. (10)), which involves the exchange of ultrasoft Majorana neutrinos, with four-momenta scaling as k 0 ∼ |k| k F . The diagram contributing to T usoft is given in Fig. 6 and its expression is 3 The amplitude T f i is related to the S-matrix element by   6. The ultrasoft contribution to 0νββ amplitude. The thick shaded lines represent nuclear bound states. The remaining notation is as in Fig. 1.
is the lowest-order nuclear weak current and |n represent a complete set of nuclear states (eigenstates of H 0 ) with threemomentum ±k + (1/2)(p i + p f ) (the ± refer to the first and second term in Eq. (39), respectively). The quantum numbers of J µ imply that, for given 0 + even-even initial and final states, |n spans the set of eigenstates of the intermediate odd-odd nucleus. Since we are in the ultrasoft regime, we expand f |J µ |n n|J µ |i in k and keep only the k = 0 term, noting that finite momentum terms would produce upon integration additional positive powers of the IR scale E 1,2 + E n − E i , and therefore additional suppression.
Evaluating the loop integral in dimensional regularization with MS subtraction, we find The UV divergence and the associated logarithmic dependence on µ us are reabsorbed by the term in the potential proportional toṼ AA . To verify this, using the completeness relation for the eigenstates of H 0 we write the term proportional to log µ us in (40) as a double commutator [28] and evaluate it using the lowest order chiral potential in H 0 , finding withṼ (a,b) AA given in (21). The µ us -independence of the total amplitude implied by Eqs. (41) is a non-trivial consistency check for our calculation and allows us to pick a convenient scale, such as µ us = m π , which eliminates the contribution ofṼ AA . Moreover, the cancellation implies that T usoft has the same chiral scaling as (V   (a,b) ν,2 ) f i , and is thus two orders down compared to the leading contribution (V   (a,b) ν,0 ) f i . This suppression can also be seen by directly comparing the scaling of T usoft and (V ν,0 ) f i . In fact, (V ν,0 ) f i ∼ 1/(4πR A ) ∼ k F /(4π) 4 , which leads to T usoft /T 0 ∼ n (E 1,2 + E n − E i )/(4πk F ) × f |J µ |n n|J µ |i . Note that the dimensionless transition matrix elements f |J µ |n n|J µ |i also control the 2νββ decay amplitude, through a different E n -dependent weighted sum [19]. The overlap matrix elements quickly die out for E n − E i > 10 MeV, as borne out in several explicit calculations using different many-body methods [64][65][66]. Therefore, recalling the scaling E n − E i ∼ k 2 F /m N we recover T usoft /T 0 ∼ 2 χ . In summary, in the chiral EFT framework one expects the following hierarchy of contributions to the 0νββ amplitude of Eq. (38): • The leading contribution is given by T 0 = −T lept × (V ν,0 ) f i with V ν,0 given in Eqs. (12) and (13). This leading term is not sensitive to the intermediate states of the odd-odd nucleus. V ν,0 corresponds to the standard neutrino potential [19] • A commonly included (but incomplete) N 2 LO contribution is obtained by inserting momentum dependent form factors in (13), as shown in Eq. (14) and the subsequent discussion.
• The new N 2 LO contribution is given by given in Eqs. (12), (18) and T usoft (µ us ) from Eq. (40). With the choice of renormalization scale µ us = m π ,Ṽ AA drops out of the calculation. Note that T usoft requires the same nuclear structure input needed in 2νββ calculations, namely f |J µ |n n|J µ |i and the excited energy levels of the intermediate nucleus (E n ).
We suggest that many-body calculations be organized according to this hierarchy, with the aim of (i) comparing results of various methods order by order in chiral EFT and (ii) checking to what degree the chiral counting is respected in large nuclei.
Finally, note that in evaluating (V ν,2 ) f i in chiral EFT and pionless EFT (and (V ν,0 ) f i in pionless EFT), one encounters a priori unknown counterterms, which can be estimated in naive dimensional analysis. In the next section we discuss how to go beyond this rough estimate.

D. Estimating the Low Energy Constants
Chiral EFT: Interestingly, ππ and πN interactions similar to those in Eq. (9) are encountered when considering electromagnetic corrections to meson-meson and meson-nucleon interactions [67][68][69][70][71]. In the electromagnetic case, these operators arise from two insertions of the electromagnetic interaction, which involves the exchange of hard photons. In the case considered here, the operators are generated by the exchange of hard neutrinos. However, the neutrino propagator and weak vertices combine to give, up to a factor, a massless gauge boson propagator in Feynman gauge (see Eqs. (3) and (5)). This formal analogy can be exploited to relate the LECs needed for 0νββ (two insertions of the τ + weak current) to the LECs associated with the ∆I = 2 component of the product of two electromagnetic currents, that belongs to the 5 L × 1 R irreducible representation of the chiral SU(2) group. Based on this observation, we have identified the operators of Refs. [67][68][69][70][71] that correspond to g ππ ν and g πN ν . Explicitly, the relation between our couplings renormalized in the MS scheme and those of e.g. Ref. [69] (which are in the modified MS scheme commonly employed in χPT [72]), is given by Our results for the anomalous dimensions of these couplings are in agreement with Ref. [67], At present the LEC g πN ν remains undetermined, while several estimates exist for g ππ ν [71,73]. For example, Ref. [71] finds in Feynman gauge κ r 3 (µ = m ρ ) = 2.7 · 10 −3 , which corresponds to g ππ ν (µ = m ρ ) = −7.6. We expect this estimate to be accurate at the 30-50% level, as it relies on a large-N C inspired resonance saturation of the correlators. Finally, electromagnetic counterterms in the two nucleon sector have been classified in Ref. [74], but as far as we know no estimate of the finite parts exist, which would give us a handle on g N N ν . A first-principles evaluation of g ππ ν , g πN ν , g N N ν based on Lattice QCD is also possible. g ππ ν and g πN ν can be determined by computing the S-matrix elements for the processes π − π − → ee and n → pπ + ee on the lattice and matching to the corresponding chiral EFT expressions. On the lattice side, one needs to compute matrix elements of the non-local effective action in Eq. (3) between appropriate external states. As discussed above, the calculation is formally very similar (modulo the Lorentz and isospin structure of the currents) to the one required to compute virtual photon corrections to hadronic processes. Techniques being developed in that context [75,76] might prove useful for 0νββ. On the EFT side, one needs to compute full S-matrix elements, not potentials. As an illustration, and because the ππ matrix element would probably be the first to be tested on the lattice, we report the N 2 LO S-matrix result for π − (q)π − (−q) → ee, with q 2 = m 2 π : The S-matrix element for n → pπ + ee cannot be readily extracted from our matching calculation, because we used off-shell "potential" pions in the external legs of Fig. 2 (bottom panel). Similarly, g N N ν can be determined by matching the lattice calculation of nn → ppee to the chiral EFT one, with a few caveats: (i) the chiral EFT S-matrix elements requires a non-perturbative calculation, which we have not performed. (ii) Matching at N 2 LO requires subtracting the ultrasoft contribution for the specific channel nn → ppee. In principle, all the ingredients (E n and pp|J µ |pn np|J µ |nn ) to evaluate T usoft can be computed in Lattice QCD, and a first step in this direction has been made in Ref. [26,58] in the context of the pionless EFT. (iii) Finally, one needs to subtract the contributions from g ππ ν and g πN ν , or alternatively extract these from the m π dependence of the lattice nn → ppee amplitude.
Pionless EFT: To determine the LO (N 2 LO) couplings g N N ν (g N N ν,2 ), one would have to match a Lattice QCD calculation of the ∆L = 2 nn → ppee scattering amplitude to a full LO (N 2 LO) calculation of the same amplitude in the pionless EFT (or the analogue amplitude for the 0νββ decay of the bound nn state, relevant for heavy pion lattices [58]). While obtaining the LO (N 2 LO) nn → ppee amplitude in pionless EFT is beyond the scope of this work, we note that part of the LO amplitude is given in Eq. (31). We also note that in performing the matching up to N 2 LO one can ignore the contribution from the ultrasoft amplitude, which in pionless EFT contributes beyond N 2 LO (see discussion in Sect. III B). Should one need to evaluate T usoft , the input quantities E n and pp|J µ |pn np|J µ |nn can be computed on the lattice [26,58].

IV. CONCLUSIONS
We have presented the first comprehensive effective theory analysis of 0νββ induced by light Majorana-neutrino exchange, describing the physics from the scale Λ LNV all the way down to the nuclear energy scale. The full 0νββ amplitude receives contributions from hard, soft, potential, and ultrasoft neutrino virtualities. Starting from the quark-level description, we have performed the matching to chiral EFT. In this context, contributions from hard modes are captured by local counterterms, while the contributions from soft and potential modes can be explicitly evaluated and lead to appropriate nuclear potentials -insensitive to properties of intermediate nuclear statesfor which we have derived LO and N 2 LO expressions (see Eqs. (12), (13), (18)). We have identified new contributions that can not be captured by parameterizations of single nucleon form factors and by power-counting arguments are as large as terms usually included in the 0νββ amplitude. The contributions from ultrasoft modes appear at N 2 LO and can be displayed in terms of nuclear matrix elements of the weak current and excitation energies of the intermediate odd-odd nuclei, that also control the 2νββ amplitude (see Eq. (40)).
In Section III D we discuss strategies to determine the low-energy constants (LECs) that appear in the potentials. We have worked out a connection to the electromagnetic LECs encoding the effect of hard virtual photons in hadronic processes, that can be obtained from model estimates, lattice QCD, and, at least in principle, from data. We have also discussed a strategy to match directly to ∆I z = 2 hadronic amplitudes that could be calculated in Lattice QCD.
While the bulk of our discussion is based on the Weinberg version of chiral EFT, in Section III B we also present the potential in pionless EFT to LO and N 2 LO. We plan to study the consistency of the Weinberg power counting for 0νββ decay in future work.
In Section III C we have discussed the hierarchy of chiral EFT contributions to the "master formula" for the 0νββ amplitude, Eq. (38), describing their relation (when applicable) to the standard treatment of 0νββ matrix elements in the literature. We advocate that many-body calculations with existing methods [77][78][79][80][81][82][83][84][85][86][87][88], as well as with methods under development [89], should be organized according to the EFT power counting scheme, isolating LO, N 2 LO, and ultrasoft contributions. It will be particularly interesting to quantify the impact of the new N 2 LO potential and ultimately assess the validity of the chiral framework.

ACKNOWLEDGEMENTS
The idea for this work was formulated at a meeting of the Topical Collaboration on "Nuclear Theory for Double-Beta Decay and Fundamental Symmetries", supported the the DOE Office of Science. VC and EM acknowledge support by the LDRD program at Los Alamos National Laboratory. WD acknowledges support by the Dutch Organization for Scientific Research (NWO) through a RUBICON grant. The work of AWL was supported in part by the Lawrence Berkeley National Laboratory (LBNL) LDRD program, the US DOE under contract DE-AC02-05CH11231 under which the Regents of the University of California manage and operate LBNL, the Office of Science Advanced Scientific Computing Research, Scientific Discovery through Advanced Computing (SciDAC) program under Award Number KB0301052 and by the DOE Early Career Research Program. VC acknowledges partial support and hospitality of the Mainz Institute for Theoretical Physics (MITP) during the completion of this work. We thank the Institute for Nuclear Theory at the University of Washington for its hospitality and the Department of Energy for partial support during the completion of this work. We thank Martin Hoferichter for suggesting to us the connection between the electromagnetic low-energy constants and the ones needed for 0νββ, and Luigi Coraggio for informative exchanges. We thank Javier Menéndez for preliminary studies of the N 2 LO potential in 136 Xe. We thank Bira van Kolck for several interesting discussions and for comments on the manuscript. We acknowledge discussions at various stages of this work with Joe Carlson, Zohreh Davoudi, Jordy de Vries, Jon Engel, Michael Graesser, Saori Pastore, and Martin Savage.