Sterile neutrinos and neutrinoless double beta decay in effective field theory

We investigate neutrinoless double beta decay ($0\nu\beta\beta$) in the presence of sterile neutrinos with Majorana mass terms. These gauge-singlet fields are allowed to interact with Standard-Model (SM) fields via renormalizable Yukawa couplings as well as higher-dimensional gauge-invariant operators up to dimension seven in the Standard Model Effective Field Theory extended with sterile neutrinos. At the GeV scale, we use Chiral effective field theory involving sterile neutrinos to connect the operators at the level of quarks and gluons to hadronic interactions involving pions and nucleons. This allows us to derive an expression for $0\nu\beta\beta$ rates for various isotopes in terms of phase-space factors, hadronic low-energy constants, nuclear matrix elements, the neutrino masses, and the Wilson coefficients of higher-dimensional operators. The needed hadronic low-energy constants and nuclear matrix elements depend on the neutrino masses, for which we obtain interpolation formulae grounded in QCD and chiral perturbation theory that improve existing formulae that are only valid in a small regime of neutrino masses. The resulting framework can be used directly to assess the impact of $0\nu\beta\beta$ experiments on scenarios with light sterile neutrinos and should prove useful in global analyses of sterile-neutrino searches. We perform several phenomenological studies of $0\nu\beta\beta$ in the presence of sterile neutrinos with and without higher-dimensional operators. We find that non-standard interactions involving sterile neutrinos have a dramatic impact on $0\nu\beta\beta$ phenomenology, and next-generation experiments can probe such interactions up to scales of $\mathcal O(100)$ TeV.

the SMEFT operators are evolved to the electroweak scale where heavy SM particles (top, Higgs, W, Z) are integrated out of the EFT. The resulting operators are evolved to the GeV scale, after which they are matched to LNV hadronic operators in chiral perturbation theory (χPT). The χPT Lagrangian is used to calculate the nn → pp ee transition operator, which, once inserted into many-body nuclear wave functions, gives rise to 0νββ of atomic nuclei. The final result is a so-called 'Master formula' that relates SMEFT operators to 0νββ decay rates in a systematic expansion in (v/Λ) α (Λ χ /v) β (m π /Λ χ ) γ , where m π is the pion mass and Λ χ 1 GeV the chiralsymmetry-breaking scale, and α, β, γ exponents that depend on the original LNV source. The formula expresses 0νββ rates in terms of a set of phase-space factors, nuclear matrix elements, hadronic low-energy constants, QCD evolution factors, and the original LNV Wilson coefficients. With this formula, any BSM model for which the SMEFT framework is applicable (Λ v and no light BSM degrees of freedom) can be directly connected to 0νββ rates.
In this work, we extend the above-sketched framework to an important class of BSM scenarios: models with additional light sterile neutrinos. Such models have been considered in light of low-scale leptogenesis [55][56][57][58], the possibility of sterile neutrinos as a dark matter candidate [55][56][57][58][59][60], and to account for anomalies in neutrino-oscillation experiments [61]. More generally, the presence of neutrino masses hints towards the existence of sterile neutrinos but not towards a specific mass scale. As such, it is important to extend the framework developed in Refs. [50,51] to include the option of light sterile neutrinos and allow for non-standard interactions that could originate at scales not too far from the EW scale, Λ ∼ 1-100 TeV. 0νββ in presence of light sterile neutrinos is not a new topic and has been investigated extensively in the literature, see e.g. Refs. [62][63][64][65][66][67]. Here we wish to go beyond these studies in several directions. First of all, we perform a systematic study in the framework of the sterile-neutrino-extended SMEFT [68][69][70]. We extend the SM not only with sterile neutrinos and the usual renormalizable interactions with SM fields, but we also include higher-dimensional operators arising from integrating out non-neutrino states that are assumed to be heavy compared to the electroweak scale. This is relevant to describe a vast class of models, from left-right symmetric models [71][72][73], where sterile neutrinos interact with heavy SU (2) R gauge bosons, to leptoquark models [74,75] and Grand Unified Theories [76].
To describe physics at the EW and lower scales, the heavy mediators can be integrated out, leading to the appearance of effective operators such as right-handed Fermi-like interactions. We extend the SM with a full set of dimension-six and -seven gauge-invariant operators including the light gauge-singlet sterile neutrinos, where light means a mass of order of the electroweak scale or below. Depending on the mass m ν R we proceed in different ways. For Λ χ < m ν R ≤ v we can integrate out sterile neutrinos before matching to hadronic operators. After integrating out ν R we obtain effective dim-3, -6, -7, and -9 operators that have already been studied in Refs. [50,51]. The resulting 0νββ rates can then be readily read from the Master formula in those works.
The situation becomes more complicated for m ν R Λ χ . In this case, ν R remains a propagating degree of freedom at hadronic scales and needs to be considered explicitly in calculations of nn → pp ee transition operators. In this paper, we derive the 0νββ transition operators induced by light sterile neutrino in the framework of chiral EFT. In particular: induced by the exchange of virtual sterile neutrinos, which in several cases contribute to the transition operator at leading order. We organize these interactions in the chiral EFT power counting, and study the neutrino mass dependence of the associated low-energy constants (LECs) • We derive the transition operators in a consistent power counting, which guarantees that LNV scattering amplitudes are properly renormalized. We further identify the set of nuclear matrix elements (NMEs) required to calculate the 0νββ half-life.
• The NMEs and the LECs in the chiral Lagrangian depend on the mass of the sterile neutrinos. We derive effective interpolation formulae grounded in QCD and χPT which allow us to smoothly interpolate between the m ν R Λ χ and m ν R Λ χ regimes. These formulae can be systematically improved by calculating pion, nucleon, and two-nucleon LNV matrix elements with nonperturbative methods for different neutrino masses.
• We address sources of theoretical uncertainties on the 0νββ half-lives. In addition to the large uncertainty on the NMEs, the often-neglected uncertainties of the LECs, which originate when matching the EFT at the quark/gluon level to Chiral EFT, is very significant. We estimate this uncertainty by conservatively varying the unknown LECs.
Our main result is an extension of the master formula obtained in Refs. [50,51] that now includes the contributions from light sterile neutrinos. As we provide a direct matching to the UV scale, the applied framework can be matched to any BSM scenario involving sterile neutrinos and be readily connected to other probes of sterile neutrinos such as LHC searches, oscillation experiments, and meson decays. In particular, our results can be used in models where sterile neutrinos play a role as dark matter or in producing the universal matter/antimatter asymmetry via leptogenesis. To illustrate the use of the framework, we end by studying several simple scenarios involving sterile neutrinos and the associated 0νββ phenomenology. The organization of the paper is as follows. In Sect. 2 we introduce the sterile-neutrinoextended SMEFT framework and discuss its evolution to the GeV scale. We discuss the matching to the low-energy EFT where heavy SM fields are integrated out, and the effects of integrating out sterile neutrinos with masses between the GeV and electroweak scale. In Sect. 3 we match the operators at the quark level to the hadronic level using χPT, the low-energy EFT of QCD. We discuss the hadronic input required to describe LNV processes at low-energies. In Sect. 4 we derive the resulting nn → pp ee transition operators by considering soft-and hard-neutrino exchange between nucleons. In Sect. 5 we present our formulae for 0νββ decay rates as a function of phase-space factors, hadronic low-energy constants, nuclear matrix elements, the neutrino mass eigenvalues, and the Wilson coefficients of higher-dimensional operators. The nuclear matrix elements and their neutrino-mass dependence are discussed in Sect. 6. The neutrino-mass dependence of hadronic low-energy constants and so-called subamplitudes are studied in Sect. 7. In Sect. 8 we illustrate some applications of the developed framework by considering several scenarios involving light sterile neutrinos. We summarize and conclude in Sect. 9. Several appendices are devoted to technical issues.

The Lagrangian in the Standard Model Effective Field Theory
We consider a Lagrangian at the scale of BSM physics Λ v that consists of the SM Lagrangian supplemented by a right-handed gauge-singlet neutrino and higher-dimensional operators. In this work, we consider operators up to dimension seven. To be precise, when discussing gaugeinvariant operators in the SM-EFT, we follow Ref. [51] and denote their dimensions by dim-n with n = 5, 6, 7. After EWSB, the EFT operators are only SU (3) c × U (1) em invariant and we refer to them, without the overline, as dim-n operators where n = 3, 6, 7, 9. dim-9 LNV operators play an important role in the phenomenology of 0νββ but the relevant operators do not involve neutrinos. Their contribution has been studied in detail in Ref. [51] and is not affected by the inclusion of light sterile neutrinos. The Lagrangian we consider is then +L (5) ν L + L (5) ν R + L (6) ν L + L (6) ν R + L (7) ν L + L (7) in terms of the lepton doublet L = (ν L , e L ) T , whileH = iτ 2 H * with H the Higgs doublet where v = 246 GeV is the Higgs vacuum expectation value (vev), h(x) is the Higgs field, and U (x) is a SU (2) matrix encoding the Goldstone modes. ν R is a column vector of n right-handed sterile neutrinos. Y ν is a 3 × n matrix of Yukawa couplings andM R a general symmetric complex n × n mass matrix. Without loss of generality we will work in the basis where the charged leptons e i L,R and quarks u i L,R and d i R are mass eigenstates (i = 1, 2, 3). This implies d i L = V ij d j, mass L , where V is the CKM matrix. The relation between the mass and weak eigenstates for the neutrinos will be discussed below. We define Ψ c = CΨ T for a field Ψ in terms of the charge conjugation matrix C = −C −1 = −C T = −C † . We use the definition for chiral fields Ψ c L,R = (Ψ L,R ) c = CΨ L,R T = P R,L Ψ c , with P R,L = (1 ± γ 5 )/2.
We now turn to the higher-dimensional operators. In general they contain all generations of quarks, but for 0νββ the most important operators are those involving the first generation of quarks. We therefore focus on operators with just u and d quarks, which implies that the Wilson coefficients will carry indices in lepton flavor only. We make one further truncation of the set of effective operators by focusing on interactions containing just one neutrino field. The only exception are operators that contribute to neutrino masses after EWSB and thus contain two neutrino fields.
The dim-5 operators obeying the above criteria are written as which after EWSB contribute to Majorana mass terms for active and sterile neutrinos. For n ≥ 2 there appears a dim-5 transition dipole operator but it does not play an important role in 0νββ. The number of operators grows when going to higher dimensions, but the operators that match at tree level to 0νββ operators is not that large. In Tables 1 and 2 we list the operators in L (6) ν L and L (7) ν L , which involve active left-handed neutrinos and were constructed in Ref. [77] and [78,79], respectively. The operators appearing in L (6) ν R and L (7) ν R involve sterile neutrinos and were first constructed in Ref. [70], they are given in Tables 3 and 4. We use the convention of Ref. [51] for the covariant derivative.  Table 1: LNC dim-6 operators [77] involving active neutrinos that affect 0νββ at tree level. Table 2: LNV dim-7 operators [78] involving active neutrinos that affect 0νββ at tree level.  Table 3: LNC dim-6 operators [70] involving a sterile neutrino that affect 0νββ at tree level. Table 4: LNV dim-7 operators [70] involving a sterile neutrino that affect 0νββ at tree level.

Evolution to the electroweak scale
To evolve the higher-dimensional operators from µ = Λ to the electroweak scale, µ m W , we briefly discuss the required renormalization group equations (RGEs). Although for several classes of operators the complete set of RGEs have been derived, in particular for the dim-6 [80][81][82] and dim-7 [79] operators without right-handed neutrinos, here we only consider the RGEs due to oneloop QCD effects. Most of the operators in Tables 1-4 do not undergo QCD renormalization or consist of a quark bilinear and evolve either like a scalar or tensor operator. For these operators the RGEs are rather simple LedQ , C LLQuH } , C T ∈ {C (6) uW , C (6) dW , C LeQu3 } , where C F = (N 2 c − 1)/(2N c ) for N c = 3, the number of colors. In addition, there are several cases for which only combinations of couplings follow a simple RGE S,T = ∓ LνQd , C where C (i) S,T follow the same RGEs as C S,T in Eq. (4). For some of the operators involving ν R the linear combinations that run like a scalar or a tensor current are more involved and lead to the following RGEs dνQLD , C dQνeH /y e , QuνL , C dLνuH /y d , C QeνuH /y e , C QuνeH /y e , where QLνuD and C dνQLD couplings induce additional operators that only contribute to neutral currents and are therefore not shown.

Matching at the electroweak scale
After EWSB where the Higgs field takes its vacuum expectation value and after integrating out SM particles with masses of the order of the electroweak scale, the SMEFT operators can be matched to a new EFT. Operators in this EFT are only invariant under SU (3) C × U (1) em gauge symmetries and only involve light quarks, charged leptons, neutrinos, photons, and gluons. For purposes of 0νββ, operators containing first-generation quarks and no photons or gluons are the most interesting and we focus on this subset of operators.
Below the electroweak scale, the Lagrangian in Eq. (1) can be matched to the following effective Lagrangian ∆L=0 + L ∆L=2 + L ∆L=0 + L ∆L=2 , where L SM now refers to interactions of dim-4 and lower of light SM fields, and M D = vY ν / √ 2. The relevant higher-dimensional operators are given by where Apart from dim-6 and -7 operators, several dim-9 operators can be induced as well. Although only a small subset is induced at the electroweak scale, almost all can be populated if a righthanded neutrino with Λ χ < m ν < m W is integrated out. We therefore list the complete set 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. [83,84] and can be written as where τ ± = (τ 1 ± iτ 2 )/2 with τ i the Pauli matrices and α, β are color indices.
where the second column of operators is related to the first column by a parity transformation. As mentioned, we only included operators with first-generation quarks and no photons. In principle, there appear dipole-type operators containing F µν and operators with heavier quarks. We have kept all generations of leptons for now. To derive the matching contributions, we applied the equations of motion of the various fields Here M e = diag(m e , m µ , m τ ) which appears in the Lagrangian as L Me = −ē L M e e R + h.c., while m u,d are the masses of the up and down quarks. Before giving the explicit matching conditions, it is convenient to first rotate to the mass basis of the neutrino fields.

Rotation to the neutrino mass basis
After EWSB the mass terms can be written as where N = (ν L , ν c R ) T and M ν is a N × N symmetric matrix (since M L and M R are symmetric matrices), with N = 3 + n. The mass matrix can be diagonalized by a single N × N unitary matrix, U , where U is the PMNS matrix containing N (N + 1)/2 phases and N (N − 1)/2 rotation angles and the m 1 , . . . , m N are real and positive. The kinetic and mass terms of the neutrinos can be written as in terms of the Majorana mass eigenstates ν = N m + N c m = ν c . The rotation to the mass basis is given by where P and P s are 3 × N and n × N projector matrices The above rotations lead to the following form for the SM charged and neutral currents, Both currents involve the combination P U , which is an 3 × N non-unitary matrix, implying that the neutral current is no longer necessarily diagonal or universal. In general, the matrix P U contains (N − n)(N + n + 1)/2 phases. In the absence of higher-dimensional operators N − n = 3 of these phases can be absorbed by the charged-lepton fields, leading to (N − n)(N + n − 1)/2 = 3(n + 1) phases and an equal number of angles [85]. In the presence of higher-dimensional operators the same re-phasings of the electron fields can still be performed, but will result in redefinitions of the Wilson coefficients of these operators. After rotating to the neutrino mass basis the operators in L ∆L=0 and L ∆L=2 , and L ∆L=0 and L ∆L=2 can be written in rather compact form. We combine the dim-6 operators into while for the dim-7 operators we obtain The dim-9 operators contain no neutrino fields and are unaffected. The dim-6 and -7 operators are now mixtures of LNC and LNV terms, as the ν fields do not have a definite lepton number. The Wilson coefficients of the dim-6 operators are given by VL P U * +c VL P s U * , C VR P U * +c VRL =C VR P s U + c VR P U , C SRL =C T P U * +c and those of the dim-7 operators become C VLL = c VL P U +C VL P s U , C VRL = c VR P U +C VR P s U , C VL P U * +c (7) VL P s U * , C VR P U * +c VR P s U * , C TL1 =C (7) TL P s U , C TL2 = c TL P U , C TR1 =C (7) TR P s U , C TR2 = c TR P U .
The operators involving ν c L,R and ν R,L contribute to the same terms in the mass basis, the only difference results from the flavor indices that are summed over, i.e. whether P or P s appears. This notation will help simplify the calculation of the nn → pp ee transition operators.
Before discussing the matching conditions, we note that although the dim-7 operators are in principle independent, this is no longer true in the approximation that the charged leptons carry zero momenta. In this case, derivatives on the charged-lepton fields can be dropped which allows one to neglect C TR1 . In addition, the derivatives in the vector-like operators can now be moved onto the quark bilinears, which, after using the equations of motion, gives rise to interactions that have the same form as the scalar dim-6 operators. As a result, the contributions of the dim-7 vector operators can be captured by the following shifts of the dim-6 scalar operators, VRL , C VLR , C VLL . (26) While working within this approximation, we will often employ the above shifts to obtain the contributions from the dim-7 vector operators, instead of writing them out explicitly.

Matching contributions to the neutrino mass terms
Finally we explicitly give the matching conditions for the various effective interactions. We only consider tree-level relations. Some one-loop matching results can be found in Ref. [86]. The mass terms are given by νH , such that the Majorana mass of active neutrinos gets dim-5 and dim-7 contributions. Additional dim-7 contributions are induced at the loop level and discussed in Ref. [50]. The Majorana mass of sterile neutrinos gets a direct dim-3 contribution and higher-dimensional corrections. The Dirac mass gets a direct dim-4 contribution and a dim-6 correction. In principle, one expects the lowest-dimensional contribution to each mass term to dominate the mass, but power-counting estimates can be violated by small dimensionless numbers such as Yukawa couplings.

Matching conditions for operators involving active neutrinos
We now turn to the dim-6 operators involving active neutrinos ν L . The Wilson coefficients of LNC operators are given by The first contribution to c V L is the SM contribution. The remaining contributions to c V L and the other couplings are from BSM interactions and can be probed in β-decay experiments [69,87]. The matching conditions for LNC dim-7 operators are such that the right-handed dim-7 coupling is not generated. The tensor operators are generated by coupling to leptons through the SM charged current and are therefore diagonal in lepton flavor space. The analogous conditions for LNV interactions involving ν L are given by LLQdH 2,ji + C LLQdH 2,ij + C LLQdH 2,ji + C for dim-6 operators and LHD 1,ji + C LHD 2,ji + 8C LLduD 1,ij + C for dim-7 operators. The expressions for dim-6 and -7 LNV operators were obtained earlier in Ref. [50], with the exception of the dim-6 contribution proportional to M T D and M L . Note that the contribution proportional to M L C (6) eW scales as Λ −3 as M L = O(v 2 /Λ) as given in Eq. (27), so that all terms scale as Λ −3 .

Matching conditions for operators involving sterile neutrinos
In analogous fashion we obtain the Wilson coefficient of the operators involving sterile neutrinos ν R . For the LNC dim-6 operators we find and for the LNC dim-7 operatorsc Analogous to Eq. (29) the right-handed coupling is not induced, while dim-7 tensor couplings are not generated for sterile neutrinos either. In the case of the active neutrinos, such tensor couplings arise from the product of a dim-6 operator involving quark fields and the SM weak interaction, the latter of which only couples to active neutrinos. The dim-6 LNV conditions are given bȳ All contributions scale as Λ −3 except for the first contribution toC (6) VL which scales as Λ −2 and is proportional to a LNC dim-6 coefficient, while the LNV source is the Majorana mass of the sterile neutrino.
Finally the dim-7 LNV Wilson coefficients become

Matching conditions for dim-9 operators without neutrinos
The matching conditions for the dim-9 operators can be taken from Ref. [50] 1 v 3 C Other dim-9 operators are induced from dim-9 contributions as discussed in Ref. [51].
For the scalar dim-9 couplings we have [51,88,89] d d ln µ C The RGEs do not depend on the lepton chirality, and we therefore omitted the subscripts L, R in Eq. (38). The equations for the C 1,2,3 coefficients are equivalent to those in Eq. (38), while the RGEs for the vector operators are given by

Integrating out sterile neutrinos with
In case one or more neutrinos have masses in the range Λ χ < m ν ≤ v, we should integrate them out before matching onto chiral perturbation theory. We can do so by writing the Lagrangian involving the heavy neutrinos as where i is a neutrino mass eigenstate index that runs over the heavy neutrinos, i.e. Λ χ < m ν i ≤ v for i ∈ {N − n H + 1, N }, with n H the number of heavy neutrinos. Furthermore, J i incorporates the interactions of the i-th neutrino that are present in L (6,7) 1 . When integrating out the heavy neutrinos, combinations of the interactions in L (6,7) will give rise to dimension-nine operators. These can be derived by making use of the equations of motion, wherem −1 ν is a diagonal N × N mass matrix for the heavy neutrinos,m −1 , and we neglected the kinetic term of the heavy neutrinos, which produces terms that are suppressed by q/m ν i . Making the same approximation for the interactions in J allows us to drop the dim-7 terms. Appendix D discusses corrections to this approximation. We obtain whereJ i are the interactions that arise from the hermitian conjugate in L (6,7) , i.e. terms involving e L,R instead ofē L,R . This leads to the following effective Lagrangian in which we can replace J → J in case we are only interested in the terms that have L = −2 (not L = 0, 2). This procedure leads to the following matching conditions for the scalar operators 3 R = 16vC 3 R = 0 , C Analogous matching contributions arise for the C i L operators, which can be obtained from the above by replacing C (9) i R → C (9) i L , C i R → C (9) i L , C 4,5 R → C ARR , C The matching conditions for the vector operators are given by Furthermore, integrating out a heavy neutrino in principle induces four-quark two-lepton operators with an additional derivative. We discuss such terms in Appendix D.

Chiral perturbation theory with (sterile) neutrinos
Below the GeV scale, a description in terms of quarks and gluons as degrees of freedom breaks down. We therefore match to an effective description in terms of pions and nucleons. To keep the connection to QCD and the higher-dimensional operators we apply the framework of chiral perturbation theory (χPT) [90][91][92][93]. χPT is the low-energy EFT of QCD and the χPT Lagrangian consists of all interaction among the effective low-energy degrees of freedom consistent with the chiral and space-time symmetry properties of the underlying microscopic theory. We apply two-flavored χPT in which pions appear as pseudo-Goldstone bosons of the approximate chiral symmetry of QCD. Up to small chiral-symmetry-breaking corrections, pionic interactions involve space-time derivatives. This feature allows for a perturbative expansion in χ = p/Λ χ where p is the momentum scale of a process. For p ∼ m π only a finite number of interactions need to be considered. Each interaction is proportional to a coupling constant, often called a lowenergy constant (LEC), whose value cannot be obtained from symmetry considerations alone. The LECs can be fitted to experimental data, calculated using nonperturbative QCD methods such as lattice QCD, or estimated based on the power-counting scheme with naive dimensional analysis (NDA) [94]. The application of χPT to neutrinoless double beta decay was developed in Refs. [50,51,83,84,95].
The extension of χPT to systems with more than one nucleon, as required for our purposes, is often called chiral EFT (χEFT) [96] and has a more complicated power counting. The nuclear scale p 2 /2m N becomes relevant in diagrams in which the intermediate state consists purely of propagating nucleons, which are enhanced with respect to the χPT counting. This leads to the need to resum certain classes of diagrams to all orders, which manifests in the appearance of bound states: atomic nuclei. The need to resum certain nuclear interactions also has important consequences for external currents [97]. Currents that are sandwiched between nuclear interactions that must be resummed can appear at lower order in the power counting than expected based on NDA. For example, the exchange of a light Majorana neutrino between two nucleons leads to a nn → pp ee transition operator whose matrix element between nuclear wave functions diverges [98,99]. This implies that a counterterm must be present, in the form a short-range nnppee operator, to absorb the associated divergence. In this work, we determine the scaling of nucleon-nucleon currents by explicitly enforcing that the nn → pp ee amplitude is renormalized.
The discussion of 0νββ mediated by light neutrinos requires the consideration of two more scales: the energy of the outgoing electrons and the neutrino mass. The electron energy is determined by the Q-value of the 0νββ reactions. All isotopes of experimental interested have Q-values around Q ∼ O(5 MeV) which is small compared to the typical momentum exchange between nucleons q ∼ k F ∼ m π where k F is the Fermi momentum of a nucleus. nn → pp ee transition operators that explicitly depend on the lepton momenta therefore give rise to suppressed 0νββ amplitudes. To explicitly consider this suppression in the power-counting scheme we assign the counting rule Q ∼ m e ∼ m π 2 χ [51]. The electron energy is of similar size as the excitation energy of the nuclear intermediate states, which are related to violation of the so called "closure approximation". In chiral EFT, corrections to closure are associated to the propagation of ultrasoft (usoft) neutrinos, with q 0 ∼ |q| ∼ Q, and, in the standard mechanism 2 , are suppressed by Q/(4πk F ) ∼ 3 χ [95] 3 . The mass of sterile neutrinos is a varying parameter, which can go from m i Q, similar to the standard mechanism, all the way to m i Λ χ . For almost massless neutrinos, m i Q, the 0νββ amplitude receives contributions from "potential" neutrinos, with (q 0 , q) ∼ (0, k F ), from hard neutrinos, with q 0 ∼ q ∼ Λ χ , and from the usoft regime discussed above. The usoft region is suppressed, unless the LO potential contribution cancels, as happens when the active neutrinos have no Majorana mass, M L = 0, all sterile neutrinos are light, and higher-dimensional operators are turned off. In this case, the potential region is suppressed by m 2 i /k 2 F [64,65], while the usoft region is comparatively less suppressed, by m 2 i /Q 2 , and the two become similar. In this paper, we do not include the contributions from the ultrasoft region, which are phenomenologically important only in the narrow region m i Q, and very small M L . We will address the intricacies of the usoft region in a forthcoming study. In several cases, the hard region gives contributions that are comparable to those from potential neutrinos. We will consider the corresponding chiral Lagrangian in Sect. 3.4.
If we increase the neutrino mass to m i ∼ m π , the usoft region disappears. In this case, soft neutrinos and pions with q 0 ∼ q ∼ m π are explicit degrees of freedom in the theory, but they correct the 0νββ transition operator only at loop level, implying a suppression by factors of χ [95]. Instead, in the region m i ∼ Λ χ , such loop corrections become large, ∼ m i /Λ χ , making this the most complicated region to describe rigorously. Finally, if m i Λ χ , the sterile neutrinos can be integrated out in perturbation theory, as discussed in Sect. 2.3.1.
Within the framework of χEFT, extended with the additional scale considerations mentioned above, the nn → pp ee transition operators have been derived for several sources of LNV. For example, Refs. [95,98,99] calculated the transition operator in the standard mechanism up to next-to-next-to-leading order in the χ expansion. Refs. [50,51] calculated the first non-vanishing contribution for, respectively, dim-7 and dim-9 LNV operators. In this work we extend these calculations to contributions from sterile neutrinos.

Chiral building blocks
The construction of the χPT Lagrangian is well documented [91,93,100], and the application to 0νββ is spelled out in Ref. [50]. Here we just repeat the main steps. It is convenient to write the QCD Lagrangian supplemented by the operators in Eqs. (22) and (23) as with q = (u d) T a doublet of quark fields, and M = diag(m u , m d ) is a diagonal matrix of the real quark masses. The external sources s, p, l µ , r µ , t µν L , and t µν R can be read from Eqs. (22) and (23). Neglecting the SM electromagnetic and neutral weak interaction, and focusing on terms that create electrons instead of positrons, we identify The quark-level Lagrangian is formally invariant under local SU (2) L × SU (2) R transformations, q L → Lq L and q R → Rq R with L and R general SU (2) matrices, provided that the spurions transform as The chiral Lagrangian that describes the exchange of potential neutrinos is then constructed by building the most general interactions that are invariant under these transformations.

The pion sector
In χPT pions are described by in terms of the Pauli matrices τ , the pion triplet π, and F 0 is the decay constant in the chiral limit. We use F π = 92.2 MeV for the physical pion decay constant, and, since we work at lowest order in χPT, we will use F π = F 0 . Under SU (2) L × SU (2) R transformations the pion field transforms as U → LU R † . It is convenient to define a covariant derivative that transforms in the same way where l µ and r µ are the external source terms given above. Quark masses explicitly break chiral symmetry and their effects are included by the spurion χ that transforms as χ → LχR † , and explicitly where B is a LEC, often called the quark condensate, related to the pion mass via m 2 π = B(m u + m d ). The LO chiral Lagrangian consists of the Lorentz-and chiral-invariant terms with the lowest number of derivatives By expanding the U field, we can immediately read off the interactions between pions, neutrinos, and electrons, that are induced by effective operators that contribute to l µ , r µ , s, and p. Contributions from the tensor sources require two additional derivatives and only appear at higher order. For those sources, interactions in the pion-nucleon sector are more relevant. Interactions with more derivatives or insertions of χ also appear at higher order, but will not be necessary for our purposes.

The pion-nucleon sector
We work with non-relativistic heavy-baryon nucleon fields denoted by N = (p, n) T characterized by the nucleon velocity v µ = (1, 0) and spin S µ = (0, σ/2). Under chiral symmetry the nucleon field transforms as N → KN with K an SU (2) matrix, belonging to the diagonal subgroup of The same matrix appears in the transformation of u = √ U → LuK † = KuR † . A nucleon covariant derivative can be defined as It is useful to introduce two more objects with convenient symmetry properties that transform as X → KXK † , with X ∈ {χ ± , u µ }. Operators relevant for 0νββ with the lowest number of derivatives are given by whereχ + = χ + − Tr(χ + )/2 and c 5 and g T are two LECs. c 5 is connected to the strong protonneutron mass splitting (m n − m p ) str and to the scalar charge g S via [101]. g T is nowadays known from lattice QCD calculations. The numerical values of all LECs are given in Table 5.
For certain LNV sources we also require the NLO corrections. Particularly important are the contributions from the nucleon isovector magnetic moment g M and the tensor form factor g T L (2) The m N in the definitions in Eq. (57) is conventional, and do not indicate that the LECs g M and g T are determined by reparameterization invariance. At the same order in the Lagrangian, Figure 1: Tree-level contributions to n → pe − ν, in the presence of non-standard vector, axial, scalar, pseudoscalar, and tensor currents. Nucleons and pions are denoted by double and dashed lines, the electron by a single line with an arrow, while ν, which is a Majorana mass eigenstate, by a single line (with no definite particle flow). The insertion of the non-standard current is denoted by a square, while strong-interaction vertices by a circle. In the case of the vector, scalar, and tensor currents, only the first topology appears because parity forbids the couplings of the current to a single pion. Both diagrams contribute to the axial current at LO, while the pseudoscalar current is dominated by the pion pole in the second diagram.
there arise recoil corrections to the axial, vector and tensor form factors Notice however that these terms contribute to the neutrino potentials only at N 2 LO, and we disregard them in what follows. Before turning towards the nucleon-nucleon sector, we first discuss the single neutron β-decay transition operator, which plays an important role in the descriptions of 0νββ induced by sterile neutrinos.

The neutron β-decay transition operator
The chiral Lagrangians of the pion and pion-nucleon sector can be used to derive the β-decay amplitude of a single neutron. This amplitude provides a building block towards deriving the 0νββ transition operators. Not all contributions to 0νββ can be captured in this way, since LNV interactions such as ππee, πN ee, andN NN N ee operators, contribute to the 0νββ transition operator without the exchange of a neutrino. The long-distance contributions from sterile neutrinos with masses below Λ χ , however, can be captured by combining two neutron β-decay transition operators that are derived here. At tree-level in χPT, there are two types of diagrams that contribute to n → p+e+ν (where ν denotes a Majorana mass eigenstate) depicted in Fig. 1. In the region m i Λ χ , loop corrections appear at next-to-next-to-leading order and will be neglected. They have been considered in Ref. [95] for the case of the exchange of a light Majorana neutrino with SM couplings. The amplitude can be written in compact form where N (p) denotes a spinor of a non-relativistic nucleon field with three-momentum p. The sources are given in Eq. (48), and include both LNC and LNV terms. Up to NLO in the chiral expansion we obtain where p and p stand for the momenta of the incoming neutron and outgoing proton. J µν T L and J µν T R start to differ only at higher orders in the chiral expansion. We define q µ = (q 0 , q) = p µ − p µ , and ε µναβ is the totally antisymmetric tensor, with ε 0123 = +1. We have written the currents in terms of form factors that depend on the momentum transfer q 2 . Up to the order we work, most form factors are constants with the important exception of g P (q 2 ). Explicitly, we obtain and the values of the (combinations of) LECs g V , g A , g M , g S , g T , and g T are given in Table 5. All form factors are O(1), except for g P (q 2 ) that is enhanced by m N /m π . We stress that the form factors g P (q 2 ), g M (q 2 ), and g T (q 2 ) are associated with an inverted power of m N in the contributions to the hadronic currents in Eq. (59).

Chiral Lagrangian induced by dimension-nine operators
The chiral Lagrangian induced by the dimension-nine operators in Eqs. (13) and (14) was discussed in Refs. [51,84]. These operators induce LNV couplings of two pions, two nucleons and one pion, or four nucleons to two electrons. The ππ interactions only obtain significant contributions from the scalar operators. Neglecting terms with more than two pions, which are only relevant at loop level or in multi-nucleon operators, the pionic Lagrangian can be written as As we will see in the next subsection, these ππ interactions, as well as the πN and N N interactions discussed below, are not only induced by dim-9 operators, but also receive contributions from the exchange of hard neutrinos. In anticipation of these additional contributions, we will write the above couplings (as well as those to be introduced below) as where we use c α β to denote the dim-9 contributions, while the remaining terms will be discussed in Sect. 3.4. The contributions from the dim-9 operators can then be written as where the couplings with right-handed electron fields are obtained by the replacement L → R. The LECs, g ππ i , were defined in [51] and their sizes can be estimated using NDA The LECs in Eq. (61) were computed in Ref. [102], and are found to be in agreement with these expectations. We report the values of the LECs in Table 5.
Pion-nucleon couplings are induced by both scalar and vector operators. For scalar operators, the πN couplings are subleading, with the exception of the operator O 1 . For vector operators, they contribute to the LO 0νββ transition operator. Expanding in pion fields, the Lagrangian has the form where C 9 . The LECs g πN 1,V and g πN V were defined in Ref. [51], and they are O(1). Finally, both scalar and vector operators induce nucleon-nucleon interactions. Following the definitions of Ref. [51] and again expanding in pion fields, we have The scaling of the nucleon-nucleon couplings follows the NDA expectation for g N N 1,6,7 , while g N N 2,3,4,5 need to be enhanced with respect to NDA in order to renormalize the nn → pp ee amplitude [51]. Explicitly, we have Currently, only NDA estimates are available for the πN and N N LECs.

Chiral Lagrangian from the exchange of hard neutrinos
In addition to the long-range contributions originating from the exchange of potential neutrinos between nucleons, mediated by the currents in Eq. (58), the 0νββ half-lives receive corrections from short-range operators, induced by the insertions of two currents connected by the exchange of hard, virtual neutrinos. The origin of these contributions can be understood by considering the effective action induced by two insertions of the interactions in Eqs. (22) and (23) In terms of Eq. (68), the long-distance potential, derived in Sects. 3.1.1 and 3.1.2, arises from the region |x − y| 1/k F where factorizing the two interactions is a good approximation. These long-distance contributions do not necessarily capture the region where |x − y| 1/Λ χ . In fact, as we will argue below, NDA and renormalization imply that this region contributes at leading order in several cases. In order to correctly describe 0νββ, the constructed chiral Lagrangian should be able to reproduce the amplitudes that result from inserting S eff between initial and final states. In cases where the |x − y| 1/Λ χ region is important, this implies that additional short-distance interactions, of the same form as those induced by the dim-9 operators, are needed at LO in the chiral Lagrangian.

Double insertions involving the vector couplings C (6) VLR,VRR
Before discussing these contributions in generality, let us consider the amplitude h f e 1 e 2 |S eff |h i , where h i,f are hadronic states, for the example of the insertion of two vector operators. Since we are interested in amplitudes without initial-or final-state neutrinos, the neutrino fields in L (6,7) will be contracted among each other. Using this fact, and neglecting electron momenta, the Dirac algebra for the leptonic part can be performed, leading to where the dots stand for terms proportional to other Wilson coefficients, as well as terms that arise from the / q term in the propagator. In this example, we will focus on the terms ∝ m i . As mentioned above, Eq. (69) will induce operators of the same form as those induced by the dimension-nine operators. In particular, the C where we introduced , and we have explicitly given it a dependence on m i . With this scaling, g ππ LR (m i ) contributes at LO to 0νββ, meaning that the |x − y| 1/Λ χ region in Eq. (68) significantly contributes.
Very similar short-distance LECs are generated by the insertions of two electromagnetic currents, where hard virtual photons are exchanged instead of neutrinos. As explained in Refs. [95,98,99], this analogy can be made precise in the limit m i → 0, which allows for a relation between g ππ LR (0) and the pion mass splitting, explicitly confirming the NDA expectations. The C (6) VLR,VRR 2 terms in principle give rise to pionic operators involving derivatives, which however induce subleading corrections to the longdistance neutrino potentials. None of the terms in Eq. (69) induce πN couplings at leading order and we neglect them here.
Additional interactions appear in the nucleon-nucleon sector. All combinations of couplings in Eq. (68) give rise to short-distance nucleon-nucleon couplings, which are expected at N 2 LO by NDA. However, as discussed in Refs. [51,98,99], in the case of the standard mechanism and several dim-9 operators they must appear at LO to guarantee that nn → pp ee amplitudes are properly renormalized and regulator independent. The chiral Lagrangian is given by where the C (6) VLR,VRR 2 terms are related by parity and therefore come with the same LEC. In addition, we omitted traces that vanish for the form of Q L,R relevant for 0νββ, but in principle could be non-zero for other isospin components. From NDA, one finds g NN i (m i ) ∼ Λ −2 χ which implies the short-range operators contribute at N 2 LO. To absorb divergences in the scattering amplitudes, however, the scaling needs to be modified into g NN i ∼ F −2 π , so that the g NN LR,ν operators in L N N contribute at LO. The coupling g NN ν was already encountered in Refs. [98,99], since it also appears in the standard mechanism.
As was the case of the ππ interactions, the LECs that appear in the N N sector can be related to LECs that appear due to the insertion of two electromagnetic currents. In this case, it is the sum of g NN ν and g NN LR that is related to electromagnetic LECs, which affect isospin-breaking observables in nucleon-nucleon scattering. As a result, this combination of couplings can be obtained from the charge-independence breaking combination of scattering lengths, ∼ a nn + a pp − 2a np , as detailed in Ref. [99]. Within pionful chiral EFT, at m i = 0 and in the MS scheme, this leads tõ C i , while the electromagnetic couplings C 1,2 were defined in Ref. [99]. Furthermore, C(µ) = O(F −2 π ) is the nucleon-nucleon contact interaction that appears at LO in the 1 S 0 channel within chiral EFT. This coupling can be obtained by fitting to the isospin conserving nucleon-nucleon scattering lengths, and, within pionful EFT and using the MS scheme, one has [99], The above equations implyg NN . This example explicitly confirms the arguments below Eq. (72).

The general case
We now discuss the general chiral Lagrangian induced by hard neutrino exchange. This involves other combinations of Wilson coefficients, as well as the terms induced by the ∝ / q term in the neutrino propagator, both can be constructed along similar lines. The induced interactions will have the form of the ππ, πN , and N N Lagrangians of Sect. 3.3 such that all of these effects can be captured by the couplings defined in that section. For the non-derivative pion couplings of Eq. (61) we have Contributions from dim-7 operators are captured by c (ν)ππ 7 i L , which are discussed in Appendix B, and all LECs scale as g ππ . The contributions that are explicitly proportional to m i arise from choosing the ∼ m i part of the neutrino propagator when performing the lepton contractions, as in Eq. (69). The remaining terms arise from the / q part of the propagator, but only contribute at the dim-7 level. The right-handed coupling c νππ i R can be obtained from c νππ i L by interchanging the L, R labels on the Wilson coefficients, L ↔ R, while leaving those on the LECs unchanged, and dropping c νππ 7 i L . The derivative ππ couplings are given by where c ππ i R can again be obtained from c ππ i L with the interchange L ↔ R and c νππ 7 i L is given in Appendix B. The LECs related to these derivative couplings scale as O(1), one of which was already encountered in Ref. [51] where it was called g ππ T = m N Λχ g ππ T,VLL (0). It should be noted that the terms proportional to g ππ S,VLL , g ππ S,VRL , and g ππ S1,S2 are generally suppressed by F 2 π /Λ 2 χ compared to the long-distance amplitudes in the limit m i → 0. In this limit these pieces only significantly contribute if the pseudo-scalar and axial couplings, that induce the long-distance contribution, are suppressed compared to the scalar and vector couplings.
The pion-nucleon couplings of (65) can be written as, where g πN α = O(1), the right-handed coupling c πN i R is given by c πN i L with L ↔ R, and the dimension-7 contributions are again relegated to Appendix B. Several of the above LECs are connected to those of Ref. [51], for which we have g πN VLL, Finally, the contributions to the nucleon-nucleon couplings in Eq. (66) are given by, The right-handed couplings c NN,νNN  1). However, apart from the terms proportional to g NN VLL,VLR , g NN S,VLL , g NN S,VRL , g NN T,VLL , and g NN T,VRL , one has to enhance the scaling of all NN LECs by Λ 2 χ /F 2 π in order to obtain renormalized amplitudes. We report the RGEs for the enhanced LECs in Appendix C. Finally, two of the above LECs are related to those discussed in Ref. [51], namely g NN VLL,

Summary
The LECs needed to construct the neutrino potential at LO, and their current determinations, are summarized in Table 5. The LECs that enter the neutron β decay operators discussed in Sect. 3.2 are well determined, either from experiment, as in the case of g A,M , which appear in SM currents, or from lattice QCD, in the case of g S , g T and B. The one exception is g T , which contributes to the tensor current at recoil order and is not very important in β decays. The evaluation of this LEC could be pursued with the same methods discussed in Refs. [103][104][105]. The ππ couplings induced by dim-9 operators have been computed in Lattice QCD [102], with uncertainty better Table 5: The low-energy constants relevant for the dim-3, dim-6, dim-7, and dim-9 operators. The headings show the type of long-distance (n → peν, π → eν) or short-distance processes the LECs induce, while the labels O (9) and O (6,7) ⊗ O (6,7) indicate whether the corresponding LECs are induced by dim-9 operators or by the insertion of two dim-6(-7) interactions. Whenever known, we quote the values of the LECs at µ = 2 GeV in the MS scheme.
than 10%. The ππ, πN , and N N couplings induced by dim-6 and dim-7 operators are functions of the neutrino mass. In the case of g ππ LR , both the small-and large-m i behavior are known, allowing us to obtain a reliable interpolation formula, as we will discuss in Sect. 7. In several other cases, only the large m i behavior is known. The calculation of these couplings as a function of m i could use techniques similar to the Hubbard-Stratanovich transformation proposed in Ref. [106] for the g ππ 1,...,5 couplings, with the difference that the scalar particle σ introduced in Ref. [106] is kept light. The determination of the pion-nucleon and nucleon-nucleon couplings, induced by dim-6, -7 and -9 operators, is much more uncertain. At the moment, only the combination g NN ν (0) + g NN LR (0) is known, in a variety of renormalization schemes [51,99], via its relation to charge-independence breaking in nucleon-nucleon scattering. All other couplings require dedicated Lattice QCD calculations of LNV nucleon-nucleon scattering amplitudes. In the literature, the LECs in Table 5 are often estimated using uncontrolled assumptions such as "factorization" of the product of two weak currents. While this might be unavoidable at the moment, we will show that varying the LECs in a range suggested by their NDA scaling introduces uncertainties in the 0νββ half-lives that are as large as those in the nuclear matrix elements, and should not be neglected. Figure 2: Tree-level contributions to the 0νββ transition operator arising from the exchange of a light neutrino. The notation for nucleons, electrons, and neutrinos is as in Fig. 1. The squares denote the nucleon vector, axial, scalar, pseudoscalar, and tensor form factors, which, at LO in chiral EFT, include one or both diagrams in Fig. 1. The currents acting on the two nucleons can be different, which we denoted by hatching one of the two squares. LNV arises from the mass of the neutrinos, which in general are Majorana eigenstates, or from the couplings of the neutrinos to the nucleons, which receive LNV contributions at dim-7.

The nn → pp ee transition operator including sterile neutrinos
We now turn to the main part of this work: the derivation of the nn → pp ee transition operator. This transition operator will be inserted between nuclear wave functions and is sometimes called the "neutrino potential". The transition operator is not necessarily due to the exchange of a neutrino as other mechanisms exist, for instance via the contact ππee, npπee and nnppee interactions discussed in Sect. 3.3 and 3.4. Such mechanisms have been discussed in detail in Ref. [51] and the derivation of the potential in the presence of sterile neutrinos amounts to generalizing the couplings C ( )ππ , C πN and C N N as in Eq. (62), to include the contributions of hard-neutrino exchange. We therefore focus here on the neutrino potential arising from the exchange of a light neutrino, with mass below the chiral-breaking scale Λ χ . In general the induced neutrino potential arises from the four diagrams in Fig. 2, where any combination of hadronic currents can be used. The top (bottom) incoming and outgoing nucleons have momenta p 1 (p 1 ) and p 2 (p 2 ), respectively, and we define q 1,2 = p 1,2 − p 1,2 . The top (bottom) electron has outgoing four-momenta k 1 (k 2 ). In diagrams (a) and (b) the neutrino then carries momentum q 11 = q 1 −k 1 = −q 2 +k 2 . In diagrams (c) and (d) the neutrino carries momentum q 12 = q 1 −k 2 = −q 2 + k 1 . In most cases, we can neglect the electron momenta in the neutrino propagators and hadronic currents. In those cases, q 11 = q 12 = q 1 = −q 2 ≡ q. Finally, we define the notation J x (i), where x = {V, A, S, P, T R , T L } and i = {1, 2}, that implies that the expression in Eq. (59) should be evaluated for nucleon i using the momenta p i , p i , and q i . We begin by studying the so-called standard mechanism of 0νββ which is the exchange of a light Majorana neutrino. We review how to derive the well-known form of the neutrino potential appearing in this scenario and how it is affected by the presence of additional sterile neutrinos that interact via left-handed currents. This warm-up calculation provides a useful guide towards obtaining the neutrino potential arising from other interactions. The calculation of the remaining terms is in principle straightforward, however, as it is rather lengthy we have checked our results by use of the Mathematica package FeynCalc [108,109].

The standard mechanism with sterile neutrinos
We start by considering the nn → pp ee transition operator arising from potential neutrinos that interact via the C (6) VLL term in Eq. (22). This term includes the SM weak interaction as can be seen from Eqs. (24) and (28). We use the same vertex for the top and bottom nucleon propagators such that diagrams (a) and (b) add coherently and the resulting factor 2 is cancelled by the 1/2! from using the same vertex twice. Diagrams (a) and (b) then sum to where u(k i ) denotes an electron spinor with momenta k i , and we introduced n L = N − n H , such that the sum runs over all neutrino eigenstates with masses below Λ χ . The potential is related to the amplitude via A = −V . This expression can be simplified into where the dots denote corrections proportional to the lepton momenta or nucleon energy which are suppressed by additional powers of χ . Similarly, the remaining two diagrams sum to where the overall sign difference is from exchanging the two electrons. Summing all diagrams then gives The product of hadronic currents can be explicitly calculated from Eq. (59) and contains parity-even and parity-odd components. As the remaining part of V ν is an even function of q, only the parity-even parts contribute to the 0 + → 0 + transitions of experimental interest. The relevant hadronic currents are therefore where we have defined the tensor operator It is useful to split the Fermi (F), Gamow-Teller (GT), and Tensor (T) operators into their separate contributions arising from vector, axial, pseudoscalar, and magnetic currents, as the corresponding nuclear matrix elements are reported in the literature. We define the combinations For the F, GT, and T functions, we have and We then obtain for the neutrino potential This expression reduces to the familiar expression for the neutrino potential for the case of 3 light Majorana neutrinos. We set n = 0 and we turn off all higher-dimensional operators except for the active Majorana mass (which is formally a dim-5 operator). In this limiting case, C where we used Eq. (17) and used m i q. The neutrino potential becomes proportional to the Majorana mass of the active neutrinos and agrees with the usual result for the standard mechanism

The general neutrino transition operator with sterile neutrinos
The neutrino potentials arising from the other interactions in Eqs. (22) and (23) can be obtained in analogous fashion to the calculation in the previous subsection. We give here the results for all combinations of interactions that lead to a non-vanishing potential when the electron mass and momenta are neglected. The limited cases in which the first non-vanishing contribution to the transition operator involves lepton momenta are discussed in Ref. [50].
The neutrino potential can be divided into three separate leptonic structures We separate the structures into three parts. The part with superscript (6) denotes contributions from the dim-6 operators in Eq. (22) and is given by VLL C Instead, V L,R,M arise from the dim-7 operators in Eq. (23). As these terms are parametrically suppressed by a power of m π /v or Λ χ /v, we relegate the explicit expressions to Appendix A.
Finally, the short-distance part of the potentials, V 5 The neutrinoless double beta decay master formula including sterile neutrinos Armed with the neutrino potentials we define the nn → pp ee amplitude by where V (q 2 ) are the neutrino potentials from the previous section, and the sum extends over all the nucleons in the nucleus. r = r n − r m is the distance between nucleons m and n and |r| = r, and the potentials are inserted between the 0 + initial-and final-state nuclei of experimental interest. The leptonic part of the neutrino potentials can be taken outside of the nuclear wave functions and we define where R A = 1.2 A 1/3 fm is the nuclear radius in terms of the atomic number A. This factor is introduced to align the definitions of the NMEs to those in the literature. The subamplitudes A L,R,M (m i ) and A L,R,M (m i ) depend on nuclear and hadronic matrix elements, the neutrino masses, and the Wilson coefficients of the higher-dimensional operators. They are discussed in detail below. We have explicitly separated contributions from light neutrinos (with masses m i < Λ χ ) and heavy neutrinos (m i > Λ χ ). If corrections due to electron masses and momenta are kept, additional terms appear, see e.g. Ref. [50,51]. [ Table 6: Phase space factors in units of 10 −14 yr −1 obtained in Ref. [110]. The last row shows the Q value of 0νββ for various isotopes, where With the definitions of the amplitudes in Eqs. (95) and (96), we express the inverted half-life for 0 + → 0 + transitions as For a derivation of this formula we refer to Refs. [50,112]. Here G 0j are electronic phase-space factors given in Table 6 that have been calculated in the literature [110,113,114]. The subamplitudes depend on a product of Wilson coefficients and hadronic and nuclear matrix elements. To keep the expressions somewhat compact, we list here only the contributions from the standard mechanism and dim-6 interactions. Contributions from dim-7 interactions are given in Appendix A. The amplitude A L , which includes the standard mechanism, is given by where M i (m i ) are combinations of LECs and NMEs defined below. Most of the terms above describe the long-distance contributions, while A (ν) L is due to the exchange of hard neutrinos and given by where the combinations of couplings, C α i,β , are defined in Sect. 3.4. The subamplitude for righthanded electrons, which does not appear for the standard mechanism, has a very similar structure. At the dim-6 level, this amplitude can be obtained by exchanging L ↔ R where A (ν) R (m i ) can be obtained by replacing C α i L → C α i R in Eq. (99). Once dim-7 operators are included there appear differences between A L (m i ) and A R (m i ) due to the dim-7 tensor operators. The explicit formulae are given in Appendix A.
The "magnetic" subamplitude A M is given by 4 Here A

(ν)
M (m i ) again describes the contributions from hard neutrinos, Finally, we have the subamplitudes related to heavy neutrinos. Since they are induced by the same ππ, πN , and N N interactions as those arising from hard-neutrino exchange, the resulting amplitudes are very similar to those mentioned above. In particular, one can obtain the dim-9 amplitudes using the following replacements, The combinations of couplings c β α are defined in Sect. 3.3. In the above expressions we have defined the combinations of NMEs and LECs which arise from the insertions of the same currents on the nucleon lines, and the combinations which appear when two different currents interfere. We explicitly denoted the dependence on the neutrino mass in Eqs. (104) and (105).

Nuclear matrix elements
While the expressions for the subamplitudes given in the previous section look complex, they actually only depend on a relatively small set of structures. It is useful to introduce the Fouriertransformed functions in r-space We define the nuclear matrix elements (NMEs) from these functions via where the tensor in coordinate space is defined as The set of NMEs have been calculated with various nuclear many-body methods and for different isotopes in the limit m i → 0. We use calculations in the quasi-particle random phase approximation (QRPA) [115], the Shell Model [116], and the interacting boson model [117,118] and their  Table 7: Comparison of NMEs computed in the quasi-particle random phase approximation [115], shell model [116], and interacting boson model [117,118] for several nuclei of experimental interest. All NMEs are evaluated at m i = 0. The NMEs are defined in Eq. (107).
values are given in Table 7. We focus on these particular calculations because, in those works, the results were presented in terms of the different components (i.e. AA, AP , P P , M M ) of the F , GT , and T long-and short-distance matrix elements. At leading order in Chiral EFT, not all NMEs are independent. The momentum dependence of g V (q 2 ) and g A (q 2 ) is a higher-order effect in χPT. Neglecting this dependence gives leading-order relations such as This and other relations were obtained in Ref. [50].

Interpolation formulae
To calculate 0νββ decay rates when sterile neutrinos are present, we require an understanding of the m i dependence of the NMEs. For certain linear combinations of the NMEs given above, this mass dependence has been explicitly calculated [62,66], but the results are often not split up as in Table 7. Furthermore, in certain cases we also require the derivative of the NMEs with respect Eq. (113) (solid) using the quasi-particle random phase approximation [115] (red) and the Shell Model [116] (blue).
to the neutrino masses. Inspired by Ref. [66], we therefore construct an interpolation formula for the m i dependence of the NMEs. For most NMEs, we know the behavior in the small and large neutrino-mass limits. For instance, for M F (m i ) we have Note that when we give an NME without an (m i ) dependence, it is implied that it is an NME given in Table 7 corresponding to m i = 0. We stress that the meaning of an NME becomes ambiguous for m i Λ χ . For example, in the standard mechanism, the contributions are proportional to the NME M V + M A . However, for large neutrino masses, m i Λ χ , the neutrino mass eigenstate should be integrated our before matching onto χEFT, which would contribute via the dim-9 operators in Eqs. (44)- (46), and its effects are no longer captured by M V + M A alone. This implies that the correct large-m i limit is not necessarily equivalent to naively taking the m i → ∞ limit in the NMEs. We discuss this issue in detail in the next section.
To nevertheless capture the m i dependence of the NMEs in the m i Λ χ region, we can construct a simple Padé approximation of order (0, 1) that interpolates between the limits in Eq. (110) We can do the same for M AA GT,T int (m i ), M AP GT,T int (m i ), and M P P GT,T int (m i ). The formulae for the tensor NMEs are less reliable due to smallness and large model dependence of the tensor NMEs. The functional form in Eq. (111) was used in Ref. [66] and was shown to agree well with the explicit m i dependence calculated in the interacting boson model.
In case of M AA GT (m i ) we can use additional information to further constrain the interpolation formula With additional constraints, we can construct an order (1, 2) Padé approximation For an NME, M , with the large m i behavior of Eq. (110), the coefficients are given by where In Fig. 3 we plot M AA GT int (m i ) and M AA GT int2 (m i ) for neutrino masses between 1 MeV and 1 GeV for 136 Xe based on two nuclear many body methods. The two interpolation formulae agree within 25% over the whole range of neutrino masses, and the associated spread is smaller than the spread between different many-body methods. We therefore use Eq. (111) for the NMEs where we do not have sufficient information to construct the more accurate approximation in Eq. (113).
Armed with these interpolation formulae it is straightforward to obtain the m i dependence of the remaining NMEs in Table 7. For the magnetic GT NME we use while for the short-distance NMEs we obtain 6.1.1 O(m 2 i ) corrections in the small neutrino mass limit From the functional form of Eqs. (111) and (113) it is obvious that the NMEs quickly saturate for small neutrino masses, m i m π , and become constant. However, in certain interesting cases it is important to understand how fast the functions become constant. For instance, as observed in Ref. [62], in scenarios with light, m i m π , sterile neutrinos and no additional higher-dimensional operators, i.e. M L = 0, the leading nn → pp ee transition operator is proportional to and vanishes. In this case, the m 2 i /q 2 correction is necessary to get a non-vanishing result. This correction can in principle be estimated from expanding the interpolation formulae in the small  Table 8: Comparison of the derivative of 136 Xe NMEs with respect to m 2 i in the quasi-particle random phase approximation [115] and shell model [116]. The NMEs are defined in Eq. (118).
Using the interpolation formula in Eq. (111), we can directly calculate M ab {F,GT,T } and we give the results in Table 8 for 136 Xe. We stress that the results for the derivatives for the NMEs in the small m i regime are associated with significant uncertainties even beyond those from the dependence on the nuclear many-body method. By using the interpolation formulae in Eq. (113) for M AA GT (m i ) instead of Eq. (111) leads to O(100%) corrections in M AA GT . More importantly, for neutrino potentials scaling as 1/q 4 , contributions from ultrasoft neutrinos can be as important as those from potential neutrinos that are considered here. We leave these corrections to future work, but stress that our results for M should be taken as order-of-magnitude estimates.

Neutrino mass dependence of subamplitudes
The master formula in Eq. (96), combined with the results presented in Ref. [51], describes all possible contributions to 0νββ from sterile and active neutrinos, capturing both the regime of heavy sterile neutrinos, m i Λ χ , and light sterile neutrinos m i Λ χ . In the former regime, the heavy neutrino is integrated out at the quark level, while in the latter regime it has to be kept as a degree of freedom in chiral EFT. In the region m i ∼ Λ χ , however, both approaches are questionable. Within chiral EFT, diagrams arising at the n-loop level give corrections ∼ (m i /Λ χ ) 2n and the loop expansion breaks down. Instead, when integrating out the heavy neutrino at the quark level, operators involving additional derivatives, ∼ ∂ m i n × O (9) , cannot be neglected as they induce corrections ∼ (Λ χ / m i ) n . This implies the m i ∼ Λ χ regime is beyond the reach of chiral EFT and of perturbative QCD methods, and it is therefore hard to treat rigorously. In this section we discuss the m i dependence of the amplitudes in more detail and employ what is known of the amplitudes in the two regimes to suggest approximate interpolation formulae that link the low-and high-mass regions.
Before discussing the interpolation it is useful to consider the two regimes in an example involving one neutrino mass eigenstate with mass m i which couples to left-handed electrons, and to right-and left-handed quark vector currents. The low-energy amplitudes depend on the neutrino masses in two ways, explicitly through the light neutrino propagator, and implicitly via the mass dependence of the low-energy constants in Eqs. (70) and (72). The resulting 0νββ amplitude, valid in the region m i Λ χ , can be written as where we include the contributions from the hard neutrino exchange amplitudes, A 1L and C (9) 4L . Using Eq. (103), we find . Although the amplitudes in Eqs. (119) and (120) look rather different from one another, one can see that they take similar forms in the large-m i limit. To naively take this limit for the long-range contributions, as discussed in Sect. 6.1, we use and neglect the magnetic contributions which lead to short-range derivative operators, subleading in the power counting. Using the above expressions, we obtain naive estimates of the LECs g ππ 1 , g πN 1 , and g NN 1 in Eq. (120). By matching terms that depend on the same NMEs in Eqs. (119) and (120), which is equivalent to matching the ππ → ee and nn → pp ee amplitudes, we obtain These equations were obtained by setting A L = A L in the regime m i Λ χ where the left-hand side should be reliable, while the right-hand side receives large contributions from loop diagrams that were neglected. This implies that Eq. (122) can only give an order-of-magnitude estimate. Nevertheless, neglecting the g NN ν (m i ) contribution, we see that these estimates are consistent with the NDA expectation and coincide with the "factorization" approximation. The value of g ππ 1 extracted from the lattice, g ππ 1 (µ = 2 GeV) = 0.36, is about 40% smaller than Eq. (122). While these estimates are not very accurate, they at least give the right scaling. This is not so clear for the LECs in the third line of Eq. (120). For instance, we can obtain where the first two terms on the right-hand side are O(1), whereas which is not clear from the right-hand side. Similarly, it is not obvious that the g NN ν term in Eq. (122) scales the same as the left-hand side which is O(1). In all of these cases, the comparison of the naive limit of Eq. (119) with Eq. (120) suggests that the hard-neutrino LECs should have a non-trivial m i dependence. As we will argue below for the case of g ππ LR , it is indeed the m i -dependence of these LECs which ensures that the matching relations are restored.

A dispersive representation
In the case of the g ππ LR we can investigate its m i dependence by exploiting the isospin relation to the pion mass splitting. Modifying the dispersive representation derived in Ref. [119] to account for a massive neutrino, we find 5 where Π LR is the vacuum matrix element of the time-ordered product of a left-handed and righthanded current, see e.g. Ref. [119]. The correlator Π LR is exactly zero in perturbation theory and in the chiral limit, making it an order parameter for spontaneous symmetry breaking. As such, the correlator falls off rapidly, Q 4 Π LR (Q 2 ) → 0, as Q 2 → ∞ and the integral in Eq. (125) is finite. This behavior leads to the Weinberg sum rules, which are discussed for example in Ref. [119]. In the large-N c limit, the correlator can be modeled by an infinite sum of narrow axial and vector resonance contributions, subjected to the Weinberg sum rules In this parametrization the integral in Eq. (125) can be done explicitly and, after imposing the Weinberg sum rules, we obtain Considering a two-resonance model with one axial and one vector resonance, the Weinberg sum rules allow us to solve for f V,A in terms of the resonance masses and F π . Using m V = m ρ = 770 MeV and m A = m a 1 = 1.24 GeV, and taking the massless neutrino limit, we obtain 5 Here we used g ππ LR (mi) as an 'effective' LEC that captures both the hard-neutrino exchange contributions as well as the loop diagrams ∼ (mi/Λχ) n which become non-negligible for mi ∼ Λχ. Explicitly, we have A(π(q)π(−q) → eLeL)|q=0 = 8C Lu c L , such that g ππ LR (mi)| eff = g ππ LR (0) in the limit mi → 0.
Here and in what follows we use the notation g ππ LR (mi)| eff → g ππ LR (mi).
which is in reasonable agreement with the determination from the pion mass splitting, g ππ LR (0) 0.8 F 2 π , see Eq. (71). Considering the large-m i limit one instead obtains Using m A ∼ m V ∼ O(Λ χ ) and (4π) 2 = Λ 2 χ /f 2 π , for for large neutrino masses the LEC scales as and the left-and right-hand sides of Eq. (124) are of the same size. We can be more precise by taking into account additional constraints on g ππ LR . Firstly, we can consider the asymptotic limit of Π LR for Q 2 → ∞. In the chiral limit, the correlator can be obtained through the operator product expansion and it is dominated by the matrix elements of dimension-six operators [119,120] lim Using this, we can rewrite Eq. (125) for m i Λ χ as in the first term, m i Q, and we can drop Q 2 in the denominator, while in the second term we use the asymptotic expression in Eq. (132) for the correlator. This gives The dependence on Λ and µ cancels after applying the sum rule in Eq. (132). Secondly, at large values of m i , the expression for g ππ LR has to reproduce the amplitude obtained from integrating out the neutrino at the quark level as done in Sect. 2.3.1. In this case, matching the ππ → ee amplitudes gives, see Eq. (124), Together with Eq. (134) this gives an expression for g ππ 4 (m i ), which can be written in terms of g ππ 4 (µ) using the RGE of this LEC, see Eq. (38), and its perturbative solution By combining Eqs. (134)-(137) we can derive an expression for the LEC g ππ 4 (µ), Note that the log m 2 i dependence has dropped out in the above expression, as one would expect. By using the two-resonance approximation and the corresponding Weinberg sum rules, the above already allows for an estimate of the LEC, for µ = 2 GeV. This is in reasonable agreement with the direct lattice-QCD calculation g ππ 4 (µ = 2 GeV) = −1.9 GeV 2 [102].
The final condition that can be imposed on g ππ LR is the known behavior in the limit m i → 0, where we have, g ππ LR (m i = 0) 0.8F 2 π , see Eq. (71). All in all, we then have five constraints, the two Weinberg sum rules Eq. (127), knowledge from the large Q 2 behavior of the correlator in Eq. (132), and the high-and low-mass limits of g ππ LR , Eqs. (135) and (71). Clearly, not all of these constraints can be satisfied in the two-resonance approximation, even though the predictions for the high-and low-mass limits, Eqs. (139) and (129) are not very far off. To incorporate all constraints we can take into account five resonances, namely the ρ, ρ(1450), ρ(1700), a 1 (1260), and a 1 (1640) [121], and we use the physical masses while fixing their couplings through the five constraints. We show the resulting g ππ LR in Fig. 4 in green, while a three-resonance approximation, for which we do not include ρ(1700) and a 1 (1640) and do not impose Eqs. (71) and (129), is shown in blue. As is clear from the figure, the three-and five-resonance case approximations agree very well.

A naive interpolation formula
It will prove useful to construct a simpler interpolation formula that can be applied to LECs for which we have less information or where the resonance model is not applicable. We follow a similar strategy as in the start of this section and impose where A gα→ḡα involves 'effective' LECs,ḡ α (m i ), that capture both the hard-neutrino exchange contributions as well as the loop corrections that become large in the m i ∼ Λ χ regime. In the limit of zero neutrino mass, we haveḡ α (0) = g α (0). In this way the interpolated amplitude has the correct limiting behavior in both the low-and high-mass regions.
The condition Eq. (140) can be used to obtain expressions for these 'effective' hard-neutrino LECs at a scale m i Λ χ in terms of LECs arising from dim-nine operators and possible longdistance contributions. To do so, we employ the large-m i limits of the NMEs, see Eq. (121), and demand that the contributions proportional to each NME in Eq. (140) match (this is equivalent to demanding that Eq. (140) not only holds for the nn → pp ee amplitude, but also for the ππ → ee and n → pπ ee subamplitudes). We then construct an interpolation formula where m 0 2 GeV is a scale at which the procedure of integrating out the heavy neutrino becomes reliable and we use Eq. (140) to obtain expressions forḡ α (m 0,i ). The 'effective' LECs in the large-m i region scale as m −2 i , due to the m 2 i in the denominator of the neutrino propagator. This scaling ensure that g α (m i ) naive reduces toḡ α (m i ) for m i → ∞, while it reproduces g α (0) in the opposite, m i → 0, limit. Using g α (m i ) naive in the interpolated amplitude, A A int (m i ), then allows us to obtain amplitudes for any value of m i .
Applying this procedure to the case of g ππ LR , we again obtainḡ ππ (140), which, in combination with Eq. (141), leads to where on the right-hand-side we use the observed value for g ππ LR (0) = 0.8F 2 π (see Eq. (71)). At m i = m 0 = 2 GeV this function approximates Eq. (135) and it has the correct logarithmic dependence on m i in the perturbative QCD regime due to the θ(m i − m 0 )g ππ 4 (m i ) term. This naive interpolation formulae is shown in red in Fig. 4, where it is compared to the results from the three-and five-resonance models. The formulae agree over the whole range of m i within 20%. Based on this success, we will use similar naive interpolation formulae for the other LECs.
It remains to understand the relations in Eqs. (122) and (123) and the m i dependence of g NN ν (m i ) and g NN LR (m i ). These equations imply that g NN LR is enhanced with respect to g NN ν in the m i ≥ Λ χ region. This can be understood from the different RGEs of these couplings. For m i Λ χ , the RGE for g NN LR receives contributions from both light neutrino exchange and the ππ → e − e − coupling, while g NN ν only from the former where While the first term is independent of m i , the second term in the RGE of g NN LR scales as Λ 2 χ /m 2 i for large m i as can be seen explicitly from Eq. (142). It is this behavior which ensures that the left-and right-hand side of Eq. (123) match.
We can now construct interpolation formulae for g NN ν (m i ) and g NN LR (m i ) similar to Eq. (142), by using the matching conditions at large neutrino masses to obtain where we use the following RGEs to express g NN 1,4 (m i ) in terms of g NN 1,4 (µ 2 GeV) Armed with the interpolation formulae we are now ready to calculate the 0νββ amplitude in Eq. (119), starting with the combination of NMEs and LECsM V ±M A . The combination M V +M A , relevant for scenarios without higher-dimensional operators, depends on NMEs such as M F (m i ) for which we use interpolation formulae of the form in Eq. (111). In addition, there is a dependence onḡ NN ν (m i ) for which we use Eq. (145). Unfortunately, the latter formula depends on two LECs, g NN ν (0) and g NN 1 (m 0 ), that have not been determined with nonperturbative methods. For a discussion of the required lattice QCD calculation, as well as recent steps towards such a determination, we refer to Refs. [122][123][124][125][126]. We use two reasonable choices for the LECs to assess the associated uncertainty. We use the NDA estimates where the latter choice is guided by the factorization approximation in Eq. (122). The choice for g NN ν (0) is dictated by the NDA expectation g NN ν (0) ∼ g NN LR (0), and by the extraction of g NN ν (0) + g NN LR (0) from isospin breaking in nucleon-nucleon scattering [99]. As discussed in Ref. [99], the value of (2F π ) 2 (g NN ν (0) + g NN LR (0)) varies between 0.2 and ∼ 2.5 in various high-quality chiral interactions, depending on the form and value of the ultraviolet regulators.
For the left-right combinationM V −M A we require, in addition to the usual NMEs, the 'effective' LECsḡ ππ LR (m i ) andḡ NN LR (m i ) for which we use the interpolations in Eqs. (142) and (144). The former is reasonably well understood, see Fig. 4, while the latter is as uncertain as g NN ν (0). In this case, we include this uncertainty by assigning a 50% error on the contribution from g ππ LR (m i ), as both effects are expected to appear at the same order. The dependence of M V ±M A on m i is depicted in the top panels of Fig. 5. The blue and red bands correspond, respectively, toM V ±M A for 136 Xe NMEs obtained with QRPA (solid) or the Shell Model (dashed). The bands are obtained by varying the LECs as discussed above and it is clear that these LECs, and not the NMEs, that provide the dominant uncertainty. At small neutrino masses,M V +M A andM V −M A are of similar size with the former being a bit larger, but this behavior changes drastically once m i increases. Around a few hundred MeV,M V −M A becomes larger due to m i dependence of g ππ LR (m i ). In the bottom panels of Fig. 5 we plot the combinations (m i /4m e )(M V ±M A ) that appear in the subamplitude A L . The amplitudes peak in the hundreds of MeV to GeV range, but the exact location depends on the underlying operators. The uncertainties are sizable, at the order-of-magnitude level, and dominated by uncertainties in the LECs. The amplitudes show a non-trivial behavior on m i which is in part due to the m i dependence of the 'effective' hardneutrino LECs. The contributions from these LECs are not included in interpolating formulae in the literature, which are purely based on the m i -dependence of the NMEs [66]. This leads to significant differences for the case ofM V −M A where g ππ LR dominates the m i Λ χ region. Similar interpolation formulae can be derived for the other LECs introduced in Sect. 3.4.2 by matching via Eq. (140) and employing the interpolation formula in Eq. (141). This procedure allows us to smoothly interpolate between the m i Λ χ and m i Λ χ limits. In Appendix E we discuss several cases that we require in Sect. 8.3.

Phenomenology
We now turn to applications of the master formula in Eq. (97) by investigating several scenarios involving sterile neutrinos. We emphasize that our purpose is not to find phenomenologically viable models of neutrino masses, but mainly to illustrate the use of the framework developed in this work. The search for sterile neutrinos is a very rich field with searches in a wide range of experiments, see e.g. Refs. [127][128][129][130][131][132], of which 0νββ is only a small, but crucial, part. The framework presented here can be used directly in future global analyses of sterile neutrinos.
In what follows we study several relatively simple scenarios. We start by considering minimal scenarios in which we extend the SM by one or two sterile neutrinos that are gauge-singlets and do not interact via higher-dimensional interactions. In these so-called 3 + 1 and 3 + 2 models, 0νββ arises solely from the Majorana masses of the sterile neutrinos. We begin by studying whether 0νββ can be measured in these minimal models and discuss the m i dependence of the resulting decay rates. After considering these cases, we turn on several higher-dimensional operators that are induced in BSM scenarios involving leptoquarks and determine the impact of such interactions on the 0νββ predictions.
The current experimental bounds on the half-lives of various isotopes are summarized in Table 9, where the expected future sensitivities are also shown. In our numerical analyses, we use the limit on T 0ν 1/2 ( 136 Xe) obtained by KamLAND-Zen, which is the strongest one at present,  Table 9: Current and future experimental limits on T 0ν 1/2 at 90% C.L. and take into account the following future prospects The prospects for the LEGEND experiment are of high interest as well, with an expected sensitivity of T 0ν 1/2 76 Ge ∼ 10 28 yr. Before we begin analyzing the scenarios with sterile neutrinos mentioned above, we first briefly discuss the scenario of 3 light Majorana eigenstates, that is, the standard mechanism. This case corresponds to only turning on the LNV operators that generate Majorana masses for the lefthanded neutrinos, M L in Eq. (7). We use the standard parametrization and write in terms of the rotation matrices W ab (θ ab , δ ab ) ij = δ ij + (δ ia δ jb e iδ ab − δ ib δ ja e −iδ ab )s ab + (δ ia δ ja + δ ib δ jb )(c ab − 1) and R ab (θ ab ) = W ab (θ ab , 0), so that in terms of the sines (cosines) of the neutrino mixing angles, s ij (c ij ), the Dirac phase δ 13 , and the Majorana phases, λ 1,2 . We set the mixing angles to their central values [121] (see Table 10). The relevant subamplitude is A L (m i ) that depends on the combinationM V (m i ) +M A (m i ), but since all mass eigenstates are at the eV scale or below, we actually only requireM V (0) +M A (0). This combination depends on the unknown LEC g NN ν (0) associated with the exchange of hard neutrinos. This contribution is usually not considered in the literature, but as demonstrated in Refs. [98,99] appears at the same order as the exchange of potential neutrinos. To calculate decay rates we marginalize over the Majorana phases and vary g NN ν (0) between ±(2F π ) −2 as discussed below Eq. (145).
The resulting 136 Xe half life is depicted in the left panel of Fig. 6 for the normal hierarchy (NH) in blue and inverted hierarchy (IH) in red as a function of the lightest neutrino mass, m. We used two sets of NMEs obtained with the QRPA (solid) and Shell Model (dashed). Around m = 5 · 10 −3 eV, the usual 'funnel' appears for the NH due to a possible cancellation in m ββ . For smaller m, the uncertainty on the half life is roughly two orders of magnitude for both the NH and IH. This uncertainty arises roughly in equal parts from the uncertainties in the LECs and Majorana phases, and in smaller amount from the change in NMEs between the QRPA and Shell Model. This can be seen more clearly by comparing to the right panel of Fig. 6 where we have set g NN ν (0) = 0 and thus ignored contributions from hard-neutrino exchange as usually done in the literature. The bands in both the NH and IH are significantly smaller highlighting the importance of pinning down the value of g NN ν (0) with nonperturbative methods.

3+1 model
The simplest scenario we investigate is the 3 + 1 model where the SM is extended by one gaugesinglet neutrino. That is n = 1 and N = 4. This model leads to two massless neutrinos and is thus ruled out by the combined atmospheric, ∆m 2 ATM 2.4 · 10 −3 eV 2 and solar, ∆m 2 SOL 7.5 · 10 −5 eV 2 , squared mass differences. Nevertheless, due to its simplicity, the scenario provides a useful toy model. In the flavor basis, we write the neutrino mass matrix as and we set M * D, The neutrino mass matrix is then described by just two parameters M * D and M R . Note that in the absence of higher-dimensional operators, (M ν ) ij = 0 for i, j = 1, 2, 3 is required by gauge invariance. This simple setup predicts two massless neutrinos, m 1 = m 2 = 0, and two massive neutrinos described by where we assumed m 3 < m 4 , while the inverted relations are It is straightforward to diagonalize the mass matrix to obtain the PMNS matrix, which can be parametrized as [133] where W ij and R ij are defined below Eq. (150) and α D,R = Arg M D,R . In this fairly simple case nearly all the mixing angles can be expressed in terms of the neutrino masses where t ij = s ij /c ij , while s 12 is unconstrained due to the two vanishing neutrino masses. One useful result of the above is As mentioned, this model of neutrino masses is too restrictive to reproduce the oscillation data. To nevertheless investigate the effect a sterile neutrino would have on 0νββ in this scenario we approximate ∆m 2 ATM 0, and set m 3 = ∆m 2 31 0.05 eV. For m 4 < Λ χ the 0νββ half life is then simply given by (Eq. (97)) T 0ν where (Eqs. (98), (119)) and (Eqs. (24) and (28)) C All other Wilson coefficients vanish. The combinations of NMEs and LECsM V,A (m i ) are defined in Eq. (104), with the inclusion of the short-range pieces in Eq. (99), and their dependence on the neutrino masses is discussed in Sects. 6.1 and 7. For m 4 > Λ χ , instead, the half life becomes T 0ν as the fourth mass eigenstate is integrated out on the quark level. In this case A L (m 4 ) is defined in Eqs. (103) and the neutrino mass dependence enters via The interpolation formulae described in Sect. 7.2 ensure that the m 4 ≤ Λ χ and m 4 ≥ Λ χ limits smoothly match. The result for the 0νββ 136 Xe rate is depicted in the left panel of Fig. 7 as a function of m 4 . It can be divided into three regions. For m 4 m π the life time increases as m 4 4 for decreasing neutrino masses. For m 4 ≥ Λ χ the lifetime becomes independent of the neutrino mass. The intermediate region is more complicated and shown in more detail in the right panel of Fig. 7. The behavior in the three regions can be understood from the neutrino mass dependence of M V (m i ) +M A (m i ) shown in Fig. 4. For small neutrino masses m 4 m π ∼ k F , the NMEs become almost independent on the neutrino mass. The dominant contribution, however, is proportional to as can be seen explicitly from Eq. (157), see also Eq. (117) and Ref. [62]. The first non-vanishing contributions in this regime are suppressed by m 2 i /m 2 π as discussed in Sect. 6.1.1 and the half life is proportional to in the regime m 4 m 3 . For large neutrino masses, m 4 > Λ χ , we need to compare the two terms in Eq. (161) which are both non-zero. The first term involves the sum over light neutrinos, but only the m 3 contribution is nonzero, and depends on U e3 −1/ √ 3 which is constant in the 3 + 1 model up to m 3 /m 4 corrections, see Eq. (157). Similarly, we can see that m 4 U 2 e4 is roughly independent of m 4 from Eq. (157), but in this case the amplitude scales as A  From Fig. 7 it is clear that for small neutrino masses, the 3 + 1 toy model predicts extremely slow 0νββ rates, orders of magnitude away from present or projected sensitivities. For m 4 ≥ 1 GeV, the half lives range from roughly 1.7 · 10 27 yr to 1.5 · 10 28 yr depending on the choice of NMEs and LECs. The uncertainty in the half life is at the order-of-magnitude level and mainly due to our poor knowledge of the short-distance LEC g NN ν (m i ). Although the uncertainty in the small-m 4 region looks small this is probably unrealistic. In this regime, the amplitudes depend on the derivative of the NMEs with respect to m i , which we estimated by expanding the interpolation formula. However, as discussed in Sect. 6.1.1, this expression might not be accurate in this region and a realistic uncertainty would likely be at the order-of-magnitude level as well. Since the decay rates are immeasurably small, this uncertainty is not too relevant.

3+2 model
We now consider a minimal 3 + 2 model, where we extend the SM with two sterile neutrinos (n = 2). This model is more realistic as it can readily accommodate the measured neutrino mass splittings, mixing angles, and CP phase. We closely follow the analysis of Ref. [134] to conveniently parametrize the 5 × 5 mixing matrix in terms of neutrino masses and PMNS parameters. The original 5 × 5 mass matrix is given by This mass matrix leads to one massless and four massive neutrinos and can be diagonalized by a 5 × 5 unitary matrix U that consists of physical parameters [134] where, in the normal hierarchy, we have Here, m l and M h are 2 × 2 mass matrices and U PMNS and R are 3 × 3 and 2 × 2 matrices The above form of U assumes m 1 = 0 making it directly applicable to the case of the NH, while in the case of the IH we instead have m 3 = 0. To account for this change we can replace m 2,3 → m 1,2 in the above, leading to a solution of U T M ν U = diag(0, m 1 , m 2 , m 4 , m 5 ), after which the mass matrix can be brought into its usual ordering by rearranging the columns of U . Although this procedure leads to a perfectly adequate parametrization of U , it does not lead to the familiar identification of the mixing angles s ij with the solar and reactor angles. To ensure that the usual U P M N S appears as the upper left-hand block in our parametrization of U (in the limit m 4,5 → ∞) we simply follow the steps of the derivation in Ref. [134], starting from Eq. (166) with m 3 = 0 instead of m 1 = 0. This leads to a form of U which can again be written as in Eq. (166), but now   3) of Ref. [134]. The rest of the derivation then closely follows that of Ref. [134], leading to a similar form of U as in the NH case, but with additional factors of R 3 . We do not wish to perform a fully general analysis of the parameters in the mixing matrix. Instead, we fix all parameters except for m 4 and m 5 . We work in the normal hierarchy and set m 2 = ∆m 2 SOL 8.58 · 10 −3 eV and m 3 = ∆m 2 ATM 0.05 eV. We pick the best-fit values for the PMNS mixing angles and Dirac phase [121]. For simplicity we set the Majorana phase to zero, α = 0. While the choice of θ 45 does not affect the unitary matrix drastically, U e2 and U e3 can deviate from the experimental values if γ 45 O(1). Taking into account this restriction, we pick moderate values for θ 45 and γ 45 . All choices of parameters are given in Table 10.
We show the lifetime T 0ν 1/2 ( 136 Xe) in the case of the NH as a function of m 4 for four different values of m 5 in the left-panel of Fig. 8. We use the QRPA NMEs of Ref. [115] and a specific value of the short-distance LECs to not clutter the plots too much. We set as discussed in Sect. 7.2. We can observe a few things. For small m 5 < m π and m 4 m 5 the half life becomes independent of m 4 and scales as m 4 5 , similar to the behavior of the 3 + 1 scenario for small m 4 (the left part of Fig. 7). This is the 'cancellation regime', where the NMEs and LECs become m i independent and For m 4 m π and m 5 ≥ Λ χ (the right part of the green and orange lines), the scenario becomes similar to a standard seesaw scenario with 3 light Majorana neutrinos and 'decoupled' heavy states, see the horizontal gray line. The lifetime becomes roughly 10 29 yr, the predicted lifetime in the NH for a massless lightest neutrino and vanishing Majorana phases. Shorter half lives are possible if one neutrino is heavy while the other is at or below the pion mass. The right sides of the blue and red curves correspond to such a scenario with a half life of roughly 2 · 10 28 yr, not too far from the projected nEXO limits. Finally, for almost degenerate m 4 and m 5 there is a cancellation leading to a peak in the half life. In the right panel of Fig. 8, we show the same results for three choices of m 5 , now including the uncertainty from the short-distance LECs and we show results for two choices of NMEs. The uncertainties are at the order-of-magnitude level. The case of the IH is shown in Fig. 9, which shows a behavior that is very similar to that of the NH scenario. In particular, at large values of m 5 the half life becomes almost independent of m 4,5 , while, in the small mass region, the scaling ∼ m 4 5 reappears. We again see a cancellation when m 4 ∼ m 5 , and the half life approximates that of a seesaw scenario with 3 light neutrinos in the large-m 4,5 region. In contrast, the absolute value of the half life does differ between the scenarios and is roughly an order of magnitude smaller in the IH than in the NH. In fact, nEXO will be sensitive to minimal 3 + 2 scenarios in the IH for which at least one mass eigenstate has mass larger than roughly 500 MeV.
We show the dependence on the Majorana phase α and the mixing angle θ 45 in the NH in the left-and right-panel of Fig. 10 for three choices of m 4 and m 5 . The blue lines correspond to decoupled mass eigenstates with m 4, 5 1 GeV. In this case, the Majorana phase is still relevant and the half life varies by roughly a factor 5 depending on α, but only has a mild θ 45 dependence. The red line corresponds to m 4 = 1 eV while m 5 = 500 GeV is decoupled. In this case, the dependence on α is significant and for α = π the half life increases by roughly six orders of magnitude compared to α = 0, while the dependence on θ 45 is far less severe.

A leptoquark scenario
In this section we illustrate the use of the EFT framework by performing the matching in the case of an explicit model of BSM physics. Our main goal is to illustrate the framework and we consider a simple SM extension. We extend the SM with right-handed neutrinos that interact with leptoquarks (LQs). Leptoquarks are hypothetical particles that convert quarks to leptons and vice versa. All possible representations of LQs are summarized in Ref. [75], and among them is a scalar LQ that transforms as an SU (3) c triplet, an SU (2) L doublet, and carries nonzero hypercharge:R (3, 2, 1/6). The LQ interaction with sterile neutrinos is given by where a, b and i, j are flavor and SU (2) indices, respectively. In addition to these interactions, we include the right-handed neutrino Majorana mass terms and Yukawa interactions as in Eq. (1). Integrating outR, one dim-6 operator in Table 3 is generated where with m LQ being the mass of theR LQ. The above operator, as well as O LνQd generated through the RGEs of Eq. (4), induces the dim-6 scalar and tensor operators in Eq. (9) below the electroweak scale. We focus on operators involving the first-generation quarks and charged leptons where the subscripts e and a denote the charged-lepton and neutrino flavor, and where a runs from 1 to n.
Of course, we also have to include the SM weak interactions and together we obtain the following matching to the operators in Eq. (24) With these nonzero Wilson coefficients the 0νββ decay rate can be directly read from the master formula in Eq. (97). We use the procedure outlined in Sect. 7.2 to obtain matching relations for the relevant LECs and we explicitly give the resulting interpolation formulae in App. E.

A 3+1 scenario with leptoquark interactions
We first consider the 3+1 scenario, but now with the additional LQ interactions described above. For simplicity we set the couplings y LR 11 y RL * 1e = 1 as deviations can be absorbed in m LQ . For the contributions from the standard mechanism we set the unknown LECs as in Eq. (175). The LQ interactions induce a large number of contributions to the subamplitude A L as given in Eq. (199). Unfortunately most of these contributions are associated to LECs we do not control. To simplify the analysis somewhat we only include contributions from purely pionic operators, where lattice results are available that can help constrain the LECs. The missing contributions from pionnucleon and nucleon-nucleon interactions can be added once more information about the LECs is obtained. The missing contributions appear at the same order as the pionic contributions we do include, and thus correspond to a significant uncertainty. For the pionic operators we require the interpolation formulae discussed in Eq. (206) of App. E. These depend on g ππ 1,2,3,4,5 that are known (see Table 5) and the unknown couplings g ππ S1 (0), g ππ TT (0), g ππ S,VLL (0), g ππ T,VLL (0). In this section we use the NDA estimates g ππ S1 (0) = −g ππ where the signs were chosen in such a way that the LECs do not change sign when varying the neutrino mass. We plot the resulting 0νββ half-life of 136 Xe in the left panel of Fig. 11 for m LQ = 10 TeV as a function of m 4 . The red dashed line denotes the limit of m LQ → ∞ corresponding the 3 + 1 scenario discussed above. We have set the Dirac and Majorana phases to zero and used QRPA NMEs. The theoretical uncertainty due to NMEs and LECs is not shown, but is similar to the 3+1 and 3+2 scenarios and thus roughly one-to-two orders of magnitude on the half life. The plot shows that for a light fourth neutrino, m 4 < 100 GeV, the LQ interactions completely dominate over the standard mechanism. In fact, the current experiments rule out 10 TeV interactions for 100 keV < m 4 < 100 GeV. In the LQ scenario, 0νββ experiments are most constraining at m 4 170 MeV, where we find m LQ > 56 TeV. Future experiments could push this towards m LQ > 150 TeV.
In the right panel of Fig. 11 we show the half life as a function of m LQ for three specific values of m 4 . We set m 4 = 10 eV, with the eV-scale being motivated by anomalies in neutrino experiments [135][136][137][138] (but see Ref. [139] as well), m 4 = 1 keV, as motivated by models of sterile neutrino DM [140,141], and m 4 = 5 GeV as motivated by studies of low-scale leptogenesis [58]. For these cases, the limits on m LQ are respectively m LQ > 1 TeV, m LQ > 3 TeV, and m LQ > 21 TeV. While it is difficult for present and future 0νββ experiments to probe scenarios with just light sterile neutrinos, these results illustrate that prospects are much better in scenarios where the sterile neutrinos have non-standard interactions.

A 3+2 scenario with leptoquark interactions
We now repeat the analysis in the 3 + 2 model. In Sect. 8.2 we concluded that detection of 0νββ is not possible with the next generation of experiments in the pure 3 + 2 model for the NH (see Fig. 8). It is interesting to investigate what the scale of BSM physics should be to make a detection possible. We take the LQ interactions as an example of BSM physics in order to determine the relevant mass of m LQ for different choices of m 4 and m 5 . To reduce the number of parameters we set y LR 11 y RL * 1e = y LR 12 y RL * 1e = 1 and use the same phases, angles and LECs as in Sect. 8.2.
In the top-left panel of Fig. 12, we set m LQ = 10 TeV and m 5 = 10 eV (blue), m 5 = 1 keV (red), m 5 = 1 MeV (green) and vary m 4 . Current experiments already rule out such scenarios Our main goal is not to exhaustively investigate the parameter space of this LQ model. Instead, we wish to highlight that future 0νββ experiments are sensitive to a class of BSM models with light sterile neutrinos if they have weak non-standard interactions with SM fields. How weak these interactions can be depends on the exact value of the neutrino masses, but scales of 1-100 TeV can be accessed. We illustrate this in Fig. 13. In scenarios with m 4 and m 5 around 1 eV, the 136 Xe half life from just the standard mechanism is roughly 10 60 yr, but nEXO can make a detection for m LQ ≤ 1 TeV. In this particular case, LHC experiments already strongly constrain such LQ masses [142,143], but for heavier neutrinos much higher LQ mass scales can be probed. For instance, for m 4 at the keV scale and m 5 at the GeV scale, next-generation experiments are sensitive to leptoquark masses up to 100 TeV.
This discussion highlights the important role that 0νββ experiments can play in understanding the nature and interactions of sterile neutrinos. While a positive 0νββ signal alone cannot unambiguously be interpreted as evidence for sterile neutrinos, 0νββ experiments, in conjunction with direct observations in laboratory experiments or astrophysical probes, can put severe constraints on the possible interactions of these particles that are very competitive with other low-energy probes and high-energy collider experiments.

Conclusions
Sterile neutrinos are natural candidates to explain the origin of active neutrino masses via the type-I seesaw mechanism [40][41][42], which, however, does not definitively point towards a specific mass scale. While the smallness of neutrino masses can be realized with very massive ν R , m ν R ∼ 10 15 GeV, with O(1) Yukawa couplings to SM leptons, the sterile neutrino mass scale can be lowered to the TeV scale [42], inducing interesting signals at the Large Hadron Collider [144], or even to the GeV or sub-GeV region [57,58]. In such scenarios the sterile neutrinos could provide a dark matter candidate and induce the correct matter-antimatter asymmetry in the Universe [57][58][59][60]. While muon, meson, and nuclear β decays constrain the interactions of sterile neutrinos with SM particles to be more feeble than the weak interaction, in generic models ν R can be charged under new forces, mediated by bosons with masses much larger than the W and Z masses. These interactions, while weak, can leave traces in high-precision experiments, including searches for 0νββ.
In this paper, we studied the impact of light sterile neutrinos with mass smaller than the electroweak scale, m ν R < v, on 0νββ, in a systematic framework which relies on a tower of EFTs. The contribution of sterile neutrinos to 0νββ has been investigated extensively in the literature [62,63,[65][66][67], in the minimal scenario in which ν R has only renormalizable interactions, namely the Majorana mass term M R and a Yukawa interaction Y ν in Eq. (1). Here we extend these works in several, significant directions • We extend the SM by adding a (family) of sterile neutrinos ν R , which are singlets under the SM gauge group. At scales larger than the electroweak scale, we allow the ν R to interact with SM degrees of freedom via a dimension-four Yukawa interaction, and the most general set of gauge-invariant higher-dimensional operators up to dim-7, in the framework of the sterile-neutrino-extended SMEFT [68][69][70]. With the exclusion of the Majorana mass term in Eq. (1), we consider operators with at most one ν R , which are the most relevant for 0νββ. The dim-6 LNC operators with active and sterile neutrinos are listed in Tables  1 and 3, respectively, while the dim-7 LNV operators are given in Tables 2 and 4. The renormalization-group equations describing the QCD evolution for the coefficients of the dim-6 and dim-7 operators is given in Eqs. (4) and (6).
• We match the SMEFT operators onto SU (3) × U (1) em -invariant operators, by integrating out heavy SM degrees of freedom. For 2 GeV m ν R < v, the sterile neutrino is also integrated out in this step, and one is left with dim-9 operators in Eqs. (13) and (14), with coefficients given by Eqs. (44), (45), and (46). If ν R is lighter, it remains dynamical in chiral EFT. In this case, at the quark level, one finds the generalized LNC and LNV β-decay Lagrangians in Eqs. (8), (9), (10), and (11), which involve, vector, axial, scalar, pseudoscalar, and tensor currents as well as derivative operators. Rotating these interactions to the neutrino mass basis gives rise to the operators in Eqs. (22) and (23), whose coefficients can be expressed in terms of the gauge-invariant operators in Sect. 2.2.
• We systematically construct the chiral Lagrangian in the presence of light sterile neutrinos with non-standard interactions, starting from the quark-level operators in Eqs. (22) and (23). This Lagrangian includes terms with explicit neutrinos, that couple to pions and nucleons via the vector, axial, scalar, pseudoscalar, and tensor currents. These are constructed in Sects. 3.1 and 3.2, and the neutron β-decay transition operator is summarized in Eq. (58). In addition, the Chiral Lagrangian contains LNV operators that couple pions and nucleons to two leptons, which are induced by the exchange of virtual light neutrinos. In Sect. 3.4 we construct, for the first time, these operators in full generality and we show that the ππee and nnppee couplings contribute at LO for several higher-dimensional operators. We discuss the dependence of the LECs of these operators on the neutrino mass, and demonstrate that they are needed to guarantee that the 0νββ amplitude has a smooth dependence on the neutrino mass.
• We identify all the chiral EFT low-energy constants required for the derivation of the 0νββ operator at LO. They are summarized in Table 5. With the exception of the recoil-order LEC appearing in the tensor current, g T , all LECs in the neutron β-decay operator are well known, either from experiment or Lattice QCD. The ππ couplings induced by dim-9 operators are also well known [102], while ππ couplings induced by hard-neutrino exchange and πN and N N couplings induced by dim-9 operators and hard-neutrino exchange are at the moment mostly undetermined. The ignorance of these LECs causes a sizable hadronization uncertainty in the 0νββ half-lives, which, though usually neglected, is often as big as the error from nuclear matrix elements, see for example Figs. 5 and 6.
• We derive the 0νββ transition operator in Sect. 4 and the master formula for 0νββ in Sect. 5. In spite of the large number of operators, we find that the final form of the amplitude and the half-lives, Eqs. (96) and (97), are rather compact and involve structures that are very similar to those found for the exchange of active neutrinos in Refs. [50,51].
In particular, they involve the same nuclear matrix elements, with the difference that, for 1 MeV < m ν < Λ χ , the NMEs acquire a non-trivial dependence on the neutrino mass.
• We study the dependence of the 0νββ half-lives on the neutrino masses. For neutrino masses below Λ χ , the dependence arises in two ways, explicitly via the neutrino propagators in the neutrino potentials and implicitly via the LECs induced by hard-neutrino exchange. In contrast, if m ν R Λ χ , the mass dependence appears through the matching coefficients of dim-9 operators. We derive interpolation formulae, grounded in QCD and χPT, for both the nuclear matrix elements, in Sect. 6.1, and the LECs, in Sect. 7. These formulae allow us to smoothly interpolate between the m ν R Λ χ and the m ν R Λ χ regimes, as shown in Figs. 5 and 11, and to get the correct chiral EFT scaling of the amplitudes. The interpolation formulae for the LECs are partially phenomenological, due to the difficulties to treat the intermediate region m π m ν R Λ χ rigorously, but can be systematically improved by calculating pion, nucleon, and two-nucleon LNV matrix elements with nonperturbative methods for different neutrino masses.
• As a consequence of the systematic construction of the chiral Lagrangian and the study of the mass dependence of NMEs and LECs, we can address all sources of theoretical uncertainties on the 0νββ half-lives, induced by active and sterile neutrinos. An important finding of this work is that the hadronization uncertainty is very significant. We estimate this uncertainty by conservatively varying the unknown LECs in the range expected on the basis of naive dimensional analysis and internal theoretical consistency. We recommend to include the hadronization uncertainty in future analyses of 0νββ. As an example, in the left panel of Fig. 6 we show the usual inverted and normal hierarchy predictions for the 136 Xe half life as a function of the lightest neutrino mass with errors bands that include the hadronization uncertainty. This can be compared with the right panel where this hadronic uncertainty has been neglected.
To illustrate the use of the developed EFT framework, we studied several scenarios with light sterile neutrinos. Scenarios with two additional sterile neutrinos can reproduce all neutrino oscillation data. In the normal hierarchy, a minimal model where sterile neutrinos only interact via mixing, leads to 136 Xe half lives that for all choices of m 4 and m 5 are (slightly) above 10 28 yr, the prospected sensitivity of next-generation experiments. In the inverted hierarchy a detection will be possible if at least one of the neutrinos has a mass above O(500) MeV. These conclusions change drastically in the presence of higher-dimensional interactions. For instance, 0νββ experiments already exclude scalar interactions at a scale of 10 TeV if one of the neutrinos has a mass around the keV scale and the other at the GeV scale. Such mass ranges appear in scenarios where sterile neutrinos can account for both Dark Matter and the matter-antimatter asymmetry [59], although to account for both at least three sterile neutrinos are necessary. Depending on the exact neutrino masses, next-generation experiments are sensitive to scales up to O(100) TeV. The framework developed here, and the master formula in (97) in particular, can be used directly to assess the impact of 0νββ experiments on any BSM scenario with light sterile neutrinos and should prove useful when comparing 0νββ with other probes of sterile neutrinos.

A Long-distance potentials induced by dim-7 operators
Here we collect the potentials induced by dim-7 interactions. As the contributions of dim-7 vector operators can be obtained from the dim-6 potentials through a shift of the scalar interactions, see Eq. (26), we only explicitly list the terms due to C (7) TL ≡ C (7) TL1 + C (7) TL2 and C (7) TR ≡ C (7) TR2 . These are given by TR − C TL ei where h ab K, sd = h ab K q 2 /m 2 π and the dots indicate that the complete dimension-seven potentials involve contributions from the dim-7 vector operators.
The resulting left-handed amplitude can be written as, TR + C TL + C The amplitude for right-handed electrons does not obtain contributions from the dim-7 tensor operators, which implies A R (m i ) = 0 + . . . , while the "magnetic" amplitude A M is given by

B Hard-neutrino exchange contributions from dim-7 operators
In this appendix we collect the contributions from dim-7 tensor operators, C TL,TR , due to hardneutrino exchange, mentioned in Sect. 3.4. As noted in Sect. 2.2.1, all contributions from the dim-7 vector operators can again be obtained from a shift of the dim-6 scalar couplings, see Eq.
where all the LECs scale as g πN α = O(1).
Finally, the dim-7 contributions in the N N sector are given by, Of the above LECs that were not encountered in Sect. 3.4.2, g NN TL,TR , g NN TL and g NN TL,V , g NN TR,V are enhanced with respect to the NDA expectation and scale as O(1/F 2 π ) and O(Λ 2 χ /F 2 π ), respectively. The RGEs that signal the enhancement of these LECs are collected in Appendix C. The LECs g NN TL,T and g NN TR,T do not need to be enhanced and scale as O(1/Λ 2 χ ). As one would expect, the above equations lead to contributions that scale as v −5 (v −6 ) for one (two) insertions of dim-7 operators. Interestingly, the terms proportional to g ππ TL,TR and g ππ TL give rise to contributions that are enhanced by Λ 2 χ /F 2 π over the corresponding long-distance contributions. In addition, we note that C (7) TL,TR combined with dim-6 vector or scalar operators induces the same LECs that were already encountered in the insertions of dim-6 operators only (see Sect. 3.4.2), namely, g α T,VLL , g α T,VRL , g α T,SLL , and g α T,SRL . The hadronic parts of C TL,TR ×C TRR,TLL × C (6) i , and, although the former involve an additional derivative this can be supplied by a factor of / q from the neutrino propagator in the latter case. This leads to a relation between the amplitudes of C TL,TR × C (6) i , resulting from the ∝ m i term in the neutrino propagator, and those of C (6) TRR,TLL × C (6) i , which result from the / q term instead. After working out the leptonic Dirac algebra one finds that the S eff in several cases differ only by an overall constant ∼ m i /v. As these are relations between the effective actions they are also respected by the corresponding long-distance amplitudes.

C Non-perturbative renormalization of LECs
Here we collect the RGEs which give rise to the enhancement of LECs over their NDA expectations. These follow from the requirement that the LO nn → ppee amplitude is finite and µ independent, see Refs. [51,98] for details. The RGEs for the LECs related to two insertions of dim-6 operators are given by, D Derivative contributions to four-quark two-lepton operators Integrating out a heavy neutrino induces operators involving four-quark fields, two-lepton fields, and a derivative, in addition to the contributions discussed in Sect. 2.3.1. These terms can be generated through dim-6 operators, proportional tom −2 ν , or through dim-7 interactions, ∝m −1 ν . When such operators are matched onto Chiral Perturbation Theory, they will give rise to interactions that come with new, unknown, LECs. For the dimension-6 contributions one expects the derivative operators to be suppressed by m q /m ν compared to the operators of Sect. 2.3.1 6 . However, there are several Wilson coefficients for which the contributions in Sect. 2.3.1 scale as ∼ v 4 /Λ 4 , while the derivative couplings can be induced by interfering with the SM coupling, C (6) VLL . This implies that there is a regime, 1 Λ 2 mq v 2 mν , where the derivative interactions are important.
Here we will refrain from constructing the full set of derivative interactions that are induced in this way. The reason being that these interactions are only relevant in fairly specific scenarios and regimes of parameter space. In addition, their construction would lead to a sizable number of new operators all of which will come with unknown LECs, meaning that we cannot do better than an order-of-magnitude estimates in this regime. We therefore restrict ourselves to listing the contributions can be captured by the usual dim-9 interactions, i.e. the terms that, through the equations of motions, can be written as m q × O (9) .
The ∼ m q terms resulting from interference with the SM-like vector operators are given by δC (9) 1 L = −4vm d C These ∼ m q terms are expected to contribute at the same order as the original derivative couplings. Although NDA counts a derivative as a factor of ∼ Λ χ m q when matching onto ChiPT, the insertion of m q allows for a chirality flip which changes the chiral symmetry properties of the operators. In the above cases this allows one to trade derivative operators that do not induce ππ terms at LO, for interactions that do induce ππ terms proportional to m q . As a result, one expects the relative scaling of Λ χ ∼ (4π) 2 m q for the derivative and m q terms, respectively, allowing us to use the latter as an order-of-magnitude estimate for the total.
We can in principle derive similar ∼ m q terms for the interference of C VLR with C VLL . However, in this case a chirality flip does not lead to chiral enhancement (the induced operators are dim-9 vector operators which do not induce ππ terms at LO), implying that the m q terms cannot be used to estimate the derivative pieces in this case. If such interference terms are important we would have to construct the derivative terms to obtain an order-of-magnitude estimate.
We take a similar approach for the derivative operators involving dim-7 interactions. Namely, we only list the pieces ∼ m q that can be written in terms of the dim-9 operators discussed in Sect. 2.3.1. These pieces can be written as δC (9) 1 L = − C TL + m d C TL + m u C TL + m u C TL + m d C TL + m d C TR m −1 ν m d C TL + m u C while possible terms ∝m −2 ν do not appear for the combinations C VLL × C TL,TR and C (7) TL,TR 2 . We again stress that in all these cases one should in principle consider the additional operators involving derivatives, along with the unknown LECs that result from them, and the results above can only serve as an order-of-magnitude estimate.
wherec (6) SR ea = 4c and a runs from 1 to n. For simplicity we neglect the evolution to lower energies and simply give the matching to the operators in Eq. (24). We obtain The first three lines arise from the exchange of potential neutrinos while the last two lines arise from hard-neutrino exchange. The couplings are given by c νππ iL = −2g ππ S1 (m i ) C The resulting amplitude A