Analysis of Light Neutrino Exchange and Short-Range Mechanisms in $0\nu\beta\beta$ Decay

Neutrinoless double beta decay ($0\nu\beta\beta$) is a crucial test for lepton number violation. Observation of this process would have fundamental implications for neutrino physics, theories beyond the Standard Model and cosmology. Focussing on so called short-range operators of $0\nu\beta\beta$ and their potential interplay with the standard light Majorana neutrino exchange, we present the first complete calculation of the relevant nuclear matrix elements, performed within the interacting boson model (IBM-2). Furthermore, we calculate the relevant phase space factors using exact Dirac electron wavefunctions, taking into account the finite nuclear size and screening by the electron cloud. The obtained numerical results are presented together with up-to-date limits on the standard mass mechanism and effective $0\nu\beta\beta$ short-range operators in the IBM-2 framework. Finally, we interpret the limits in the particle physics scenarios incorporating heavy sterile neutrinos, Left-Right symmetry and R-parity violating supersymmetry.


I. INTRODUCTION
The nature of neutrinos and especially the origin of their masses are a crucial open question. While the Standard Model (SM) successfully explains the masses of the charged fermions it must be extended to incorporate neutrino masses. It would either require the presence of sterile neutrino states or effective lepton number violating (LNV) interactions.
The first scenario allows the generation of Dirac neutrino masses analogous to those of the charged fermions. While certainly feasible, given the stringent upper limits m ν O(0.1) eV on the absolute neutrino masses from Tritium decay [1,2] and cosmological observations [3], tiny Higgs Yukawa couplings are required. Also, total lepton number L will no longer be an accidental symmetry. Unless L symmetry is imposed, the sterile neutrinos would acquire an LNV Majorana mass. The most popular example for such a scenario is the Seesaw mechanism where the sterile neutrinos have such a large Majorana mass M ≈ 10 14 GeV naturally leading to light neutrino masses m ν ≈ 0.1 eV [4][5][6][7][8].
High-scale seesaw mechanisms, or more generally scenarios where L is broken at very high scales, are not the only way to generate light Majorana neutrino masses; other possibilities include LNV at low scales in secluded sectors, at higher loop order and when allowing higherdimensional effective interactions. If L-breaking occurs close to the electroweak (EW) scale, higher-dimensional LNV operators can be important. From a phenomenological point of view, searching for processes that violate total L thus play a crucial role in neutrino and Beyond-the-SM (BSM) physics. We here focus on the search for 0νββ decay as the most sensitive approach to probe Majorana neutrino masses. Currently, the most stringent limit on the 0νββ decay half life T 1/2 is set in the Germanium isotope 76 32 Ge [9], T 1/2 ( 76 Ge) ≡ T 1/2 76 32 Ge → 76 34 Se + e − e − > 1.8 × 10 26 yr.
However, Majorana neutrino masses are not the only contribution from BSM physics to 0νββ decay. We can generally consider the 0νββ decay rate by expressing high scale new physics contributions in terms of effective low-energy operators [10][11][12][13]. This only assumes that there are no exotic particles beyond the SM below the 0νββ energy scale of m F ≈ 100 MeV.
In this paper, we concentrate on so called short-range operators and their interplay with the standard light Majorana neutrino mass mechanism. As context, we provide here a brief overview of the possible mechanisms for 0νββ decay which can be categorized in two main classes: (i) Long-range transitions via exchange of a light neutrino. This includes the so-called standard mass mechanism in Fig. 1 (a), for which the 0νββ decay rate can be estimated Here, G F is the SM Fermi coupling and the phase space scales as Q 5 ββ with the kinetic energy release Q ββ = O(1 MeV) for typical double beta decays. Specifically, the mass mechanism of 0νββ decay is sensitive to the effective neutrino mass m ββ = i U 2 ei m ν i , summing over the light neutrino masses m ν i weighted by the square of the charged-current leptonic mixing matrix elements U ei . The inverse 0νββ decay half life in a given isotope is then conventionally expressed as with the phase space factor (PSF) G ν and the nuclear matrix element (NME) M ν . The normalization with respect to the electron mass m e yields a small dimensionless parameter | ν | = |m ββ |/m e . The current experimental bound in Eq. (1) sets a limit |m ββ | 79 − 180 meV at 90% confidence level (CL) and an unquenched axial coupling g A = 1.27 [9], with the uncertainty mainly due to the NMEs in different nuclear structure models. Future experiments will probe |m ββ | ≈ 20 meV [15], corresponding to the lowest possible value for an inverse hierarchy of the light neutrinos.
In BSM scenarios with exotic neutrino interactions, a neutrino mass insertion is not necessarily required, cf. Fig. 1 (b). In such cases, the decay rate is estimated as Γ 0νββ ββ ∼ (10 5 GeV/Λ O 7 ) 6 (10 26 yr) −1 , with the SM Higgs vacuum expectation value (VEV) v = 246 GeV and the scale Λ O 7 of the exotic dimension-7 operator. Such long-range mechanisms have received considerable attention, see e.g. [16][17][18][19][20], as the suppression at dim-7 is still fairly low and 0νββ decay is sensitive to high scales. We note, though, that due to the neutrino helicity-flip intrinsic in the exotic operator, typical mechanisms are indirectly suppressed by the light neutrino masses. It is generally difficult to have a dim-7 operator where the exotic long-range contribution dominates over the standard mass mechanism [21].
(ii) Short-range contribution where all mediating particles are heavier than m F ≈ 100 MeV, cf. Fig. 1 (c), represented as contact interactions with six external fermions.
These are the main focus of our analysis and they are generated by dim-9 and higher odd-dimensional operators. For a dim-9 operator, the decay rate can be estimated as Γ 0νββ SR ∼ Λ −10 O 9 m 6 F Q 5 ββ ∼ (5 TeV/Λ O 9 ) 10 (10 26 yr) −1 , with the operator scale Λ O 9 . The inverse 0νββ decay half life triggered by such a mechanism is expressed similarly to Eq. (2) as T −1 1/2 = | I | 2 G I |M I | 2 , with the PSF G I and NME M I , both depending on the Lorentz structure of the effective operator. The coupling constant I parametrizes the particle physics dynamics, i.e. the masses of the heavy states integrated out and their couplings.
A detailed analytic derivation of the relevant NMEs for short-range operators was provided in our previous paper [22], where we included additional NMEs that become important when the latest values of the nucleon form factors are taken into account. Moreover, we calculated PSFs using the exact radial wave functions and we presented the single electron energy and angular correlation distributions for the exotic short-range 0νββ decay mechanisms. In the present paper, we numerically evaluate all relevant NMEs within the IBM-2 framework. This will allow us to set upper limits on the effective couplings I where we will highlight the exchange of heavy sterile neutrinos as an important example. Within the same framework, we also provide updated NMEs for the standard light neutrino exchange and we analyse its interference with short-range mechanisms. The NMEs for the 0νββ transitions are generally difficult to calculate and the limits derived are affected for any contribution.
The paper is organized as follows. We summarize the effective short-range Lagrangian at the quark level in Sec. II together with examples of underlying particle physics scenarios.
The calculation of the 0νββ NMEs in the IBM-2 NME framework is outlined in Sec. III and that of the PSFs in Sec. IV. We then present our numerical results in Sec. V where we provide up-to-date limits on the standard mass mechanism and effective short-range 0νββ operators. Sec. VI concludes our discussion with a summary and an outlook.
In general, new physics where lepton number is broken at a a high scale will induce SM effective operators of dimension-5, 7, 9 and higher [33,34]. After EW symmetry breaking, this will give rise to long-and short-range contributions to 0νββ decay as outlined in the introduction, cf. Fig. 1. In this work we focus on short-range contributions and their potential interplay with the standard mass mechanism.

A. Effective Lagrangian
The general effective short-range interaction Lagrangian can be written in terms of five different Lorentz-invariant classes of fermion current products [11], where the sum is over all unique combinations of the quark and electron currents involved, Here, the 4-component Dirac spinor operators representing the up-quark, down-quark and electron are denoted by u, d and e, respectively. Quark SU (3) C colour indices are denoted by a, and each quark current forms a colour singlet in our parametrization. As the lepton current must violate lepton number by two units, the charge conjugate electron field e c appears there. Note that the chirality assignment in j R,L is flipped, i.e. the index L is associated with 1 + γ 5 . This is due to the appearance of the charge-conjugated electron field and, for example, the operatorē(1 + γ 5 )e c describes the creation of two left-handed electrons. Furthermore, the usual definition σ µν = i 2 [γ µ , γ ν ] is used. The normalization of the Lagrangian by the factor G 2 F cos 2 θ C /(2m p ) with the Fermi constant G F , the SM Cabibbo angle θ C and the proton mass m p is conventional and results in dimensionless couplings χ i . In principle, each unique current combination will be associated with a separate coupling, . Note that in Ref. [11], the Lagrangian is defined without the factor cos 2 θ C .
We chose to include it as the resulting PSFs can be defined in the same way as that for standard light neutrino exchange, cf. Sec. IV C.
Not all possible combinations of chiralities have to be considered in the Lagrangian Eq. (3), as redundancies and cancellations occur. First, the identity implies that terms corresponding to RLL

B. Example New Physics Scenarios with Short-Range Contributions
To illustrate the generation of different short-range contributions, we consider three well know scenarios beyond the SM.

Light and Heavy Neutrinos
As discussed in the introduction, the exchange of light active Majorana neutrinos is the most prominent mechanism for 0νββ decay. As a long-range contribution, it is not represented in the Lagrangian Eq. (3) but arises from the SM charged current The sum is over the three SM neutrino mass eigenstates ν i , constructed as the Majorana spinors ν i = ν i,L + ν c i,L from the SM active left-handed neutrinos ν i,L and their chargeconjugates. This gives rise to the mass mechanism of 0νββ decay sensitive to the effective Majorana neutrino mass The 0νββ decay half life in a given isotope is then conventionally expressed as in Eq. (2).
One of the most attractive extensions of the SM involves adding fermionic states ν i,S (i = The sterile states participate in the leptonic charged current due to mixing with the active neutrinos, where V eN i are the elements of the active-sterile mixing matrix. If the sterile neutrinos are much lighter than the nuclear physics scale p F ≈ 100 MeV, their contributions to 0νββ decay will be completely analogous to that of the active neutrinos and they can be included in Eq. (2) by replacing Note that the U ei , and hence m ββ , as well as the V eN i are in general complex numbers and cancellations can occur. In fact, if the Majorana states N i are solely responsible for the light neutrinos masses in a Seesaw scenario, the active and sterile contributions cancel to zero.
If instead the sterile states are much heavier than the nuclear physics scale, m N i 100 MeV, they can be integrated out, resulting in a contribution of the type J µ L J L,µ j L and the associated coupling LLL 3 is matched with the underlying physics parameters as Note that the above considerations apply for sterile neutrinos that are Majorana fermions.
This includes quasi-Dirac states that can be described by pairs of Majorana neutrinos (N 1 , N 2 ) with a small mass splitting |m N 1 − m N 2 | m N 1,2 and a relative CP phase of π/2, In the limit of Dirac sterile neutrinos with m N 1 = m N 2 , the contributions to 0νββ decay cancel.

Left-Right Symmetry
The minimal Left-Right symmetric model (LRSM) is based on the extended gauge sym- [35][36][37]. It has a rich neutrino and 0νββ decay phenomenology as it naturally contains right-handed Majorana neutrinos N i (i = 1, 2, 3) that are charged under the SU (2) R part of the gauge group, forming a doublet together with the right-handed leptons. This gives rise to right-handed charged currents, mediated by a right-handed W R boson with the gauge coupling strength g R of the SU (2) R group. The angle θ R C and the mixing matrix U R are the right-handed equivalents of the Cabibbo angle and the PMNS matrix. The LRSM gauge group is understood to be spontaneously broken to the SM at a high scale giving masses to the right-handed W R boson and neutrinos N i . In turn, the active SM neutrino acquire masses via mixing with the heavy neutrinos (Seesaw type I) as well as via the VEV of an electroweak triplet Higgs scalar present in the model (Seesaw type II).
Hence, the standard light neutrino and the sterile heavy neutrino contribution described above are generally present. In addition, the equivalent diagram with a heavy neutrino and two W R bosons contributes, giving rise to the short-range operator associated with RRR 3 matched to the underlying physics parameters as where g is the SM SU (2) L gauge coupling strength. Note that the contribution is not suppressed by the small light-heavy neutrino mixing but instead by the expectedly high W R mass m W R . The right-handed mixing matrix U R is approximately unitary with elements of order one, although cancellations due to complex phases can occur.
The SM W and the W R boson are also expected to mix with an angle as large as sin θ W LR g R m 2 W /(gm 2 W R ). This permits the right-handed lepton current to couple with a left-handed quark current mediated by the SM W giving rise to the contributions With the W mixing taking the generic value sin θ W LR ≈ g R m 2 W /(gm 2 W R ) ≈ f LR , all three effective couplings are of the same order. As mentioned, the LRSM also has the standard contribution from m ββ and the sterile neutrino contribution LLL 3 in Eq. (11). In addition, the LRSM in principle also gives rise to the remaining contributions of type 3 , namely = RLL 3 and RRL 3 but these are suppressed by both the light-heavy neutrino mixing and the high W R mass. Furthermore, the LRSM gives rise to additional long-range contributions that are not directly suppressed by the light neutrino masses.
Finally, the LRSM has contributions from the electroweak triplet scalars ∆ L,R that acquire analogous to Eqs. (13) and (14). Here, the heavy neutrino masses m N i appear because the couplings of the triplet Higgs to the gauge boson and electrons are proportional to v R and the heavy neutrino Yukawa coupling, m N ∼ y N v R . Likewise, there are contributions from the left-handed ∆ ++ L but they are additionally suppressed by the light neutrino masses (instead of m N i ) and thus negligible.

R-Parity Violating Supersymmetry
As the final example of an ultraviolet-complete theory, we consider the minimal supersymmetric Standard Model (MSSM) with R-parity violation [38,39]. Without explicitly imposing invariance under the discrete R symmetry where each field carries the multiplicative quantum number R = (−1) 3B+L+2S , with the baryon number B, total lepton number L and spin S, the MSSM allows for the R-parity breaking terms in the superpotential. Here, the indices i, j, k denote flavour generations of the superfields L,Ē, Q,D andŪ , associated with the SM weak lepton doublet L, the lepton singlet e c , the quark doublet Q and the quark singlets d c , u c . Short-range contributions to 0νββ are induced by the second term in Eq. (16), namely that associated with λ 111 for the first lepton and quark generations [40]. are generally functions of all supersymmetric particle masses and couplings involved. We here follow the assumptions of gluino dominance [41] where the diagrams involving gluinos and squarks contribute, Here we also assume degeneracy of squark masses mq = mũ L = md R in line with Ref. [41]. In addition, mg is the gluino mass and α s = 0.127 is the strong fine structure constant at m W .
Note that the gluino dominance assumption is based on the relevant NME values and limits on supersymmetry particle masses from other sources and may thus not be appropriate in light of new results. We nevertheless adopt it for simplicity and to compare with Ref. [41].

III. DETERMINATION OF NUCLEAR MATRIX ELEMENTS
The NMEs for short-range mechanisms have been analytically derived in [22]. We follow the approach therein and summarize the basic formalism using nucleon form factors.

A. Nucleon Form Factors
The nucleon matrix elements of the colour-singlet quark currents in Eq. (3) have the structure [42] p|ū where τ + denotes the isospin-raising operator which converts a neutron into a proton, and the tensor J µν in Eq. (21) is defined as The above matrix elements generally depend on the neutron and proton momenta p n = p N and p p = p N , respectively. The nucleon form factors are then functions of the momentum transfer q = p p − p n . The most general parametrization of the vector current in Eq. (20) would include also induced scalar and axial-tensor terms -these can be, however, safely neglected, since they vanish in the isospin-symmetric limit and they are not enhanced by any other effects [43].
The momentum dependence in Eqs. (19) - (21) is encoded in the nucleon form factors F X (q 2 ) with X = S, P , V, W, A, P, T 1 , T 2 , T 3 , usually parametrized in the so-called dipole Here, the so called charge g X represents the value of the form factor at zero momentum transfer, g X ≡ F X (0), and the scale m X determines the shape of the form factor. We apply this parametrization to all form factors except for the pseudoscalar form factors F P (q 2 ) and F P (q 2 ) which are enhanced by the pion resonance.
The form factors with their corresponding parametrizations and charges are given by [42].
The shape parameters are m V = 0.84 GeV, m A = 1.09 GeV [46] and the pion mass is m π = 0.138 GeV. The form factors F V (q 2 ), F W (q 2 ) and F A (q 2 ) can be determined experimentally and the parametrizations shown above provide a good description in the range 0 ≤ |q| ≤ 200 MeV of interest in 0νββ decay. On the other hand, as it is not possible to directly obtain the induced pseudoscalar form factor from experiment, we use the parametrization suggested in Ref. [45], which is based on the partially conserved axialvector current (PCAC) hypothesis. The corresponding value of the free g P charge agrees with the recent chiral perturbation theory analysis [47], which yields the value g P = 233.
The value is also consistent with measurements of muon capture. With the muon mass m µ = 0.105 GeV, the resulting value of F P (−0.88m 2 µ ) = 8.0 agrees well with the measured value of F P (−0.88m 2 µ ) = 8.06 ± 0.55 [48]. The scalar and pseudoscalar charges, g S and g P , come from recent lattice QCD calculations [44]. As there is not much information on the q 2 -dependence of the corresponding form factors, we use the dipole parametrization, which, in the Breit frame, is the Fourier transform of the matter distribution. In the case of the pseudoscalar form factor we also include the monopole factor 1/(1 + q 2 /m 2 π ) used in chiral perturbation theory. As for the tensor form factors, only F T 1 enters our calculations. The value of the corresponding charge g T 1 quoted by Ref. [44] reads 0.987 ± 0.055. We emphasize that the charges in Eqs. (23) - (29) are applicable at the free nucleon level. When calculating the 0νββ decay NMEs we will use an effective axial-vector charge g A = 1.0 and, consequently, an induced pseudoscalar charge g P (g A = 1.0) = 182 to approximately account for quenching in the nuclear medium.

B. Nuclear Matrix Elements
The five different types of quark current products appearing in Eq. (3) are mapped to the nucleon matrix elements according to Eqs. (19) - (21). By virtue of a non-relativistic expansion and the closure approximation, the resulting product of nucleon matrix elements is then mapped to the nuclear matrix element between the final and initial 0 + nuclear states involved in the 0νββ decay. This procedure is described in Ref. [22] and we here summarize the definition of NMEs involved. One should note that in the following expressions the relative sign between GT and T terms is different than in our previous papers [23][24][25]49] and other available literature taking into account tensor terms using the formulation in [45].
The confusion about the relative sign arises from Eqs. (13) and (22) in [45], where in Eq. (13) a minus sign is used in front of the tensor term, while in Eq. (22) the plus sign is used. The tensor term contributes very little to the standard long range mechanism but in case of short range mechanisms it has a notable effect. Thus we have checked the derivation and concluded that the following signs should be used.
The NMEs for the five short-range operators will generally depend on the chiralities of the two quark currents involved. For the first three operators associated with χ 1 , χ 2 and For the operators associated with χ 4 and χ 5 , the two quark currents involved have different Lorentz structures and thus all four possible combinations of chiralities have to be considered in principle: RR, LL, RL and LR. Again, it turns out that the NMEs only distinguish between the case where the quark chiralities are the same (RR, LL → upper sign) or different (RL, LR → lower sign), In the above expressions, we have explicitly factored the form factor charges g X = F X (0).
The q-dependence arising from the product of the reduced form factors F X (q 2 )/g X is still to be included in the various matrix elements appearing in Eqs. (30)- (34). The individual Fermi (M F ), Gamow-Teller (M GT ) and tensor (M T ) NMEs along with the associated reduced form factor productsh(q 2 ) are given in Table I. The numerical values of these NME will be given in Sec. III C but we would like to note that the so called recoil NMEsM AP GT andM AP T , and the NMEs explicitly depending on the temporal momentum transfer q 0 , M q 0 P P GT , M q 0 P P T are difficult to evaluate exactly. We instead assume that the sum of nucleon spatial momenta [16][17][18], approximately applicable in an average sense considering that the NME is calculated summing over all nucleons in the nucleus. Similarly, we take the average value q 0 ∼ q 2 /m p ≈ 10 MeV [18] for the temporal component of the momentum transfer. This allows to reduce the corresponding NMEs as indicated in Table I.
In addition to the product of the reduced nucleon form factors, the NMEs listed in Table I also contain the so called neutrino potential describing the q dependence of the underlying particle physics mediator of 0νββ decay. Here we follow the formulation of [45] and [24] where the two-body transition operator is constructed in momentum space as the product of the neutrino potential v(q) times the product of the reduced form factorsh(q 2 ). In the case of the short-range mechanisms we consider here, the neutrino potential is especially simple; as point-like operators, they are described by a Dirac delta function in configuration space, δ(r a − r b ), hence in momentum space it is a q-independent constant. Following the usual normalization the short-range neutrino potential is [24,45] v( We also consider the standard light neutrino exchange mechanism with the NME Note that this is fully analogous to M 3 in Eq. (32) in the case where the quark currents have the same chirality, but the crucial difference is that the NMEs in Eq. (36) are calculated with the appropriate neutrino potential in momentum space [24], Here, the neutrino mass has been neglected in comparison with the neutrino momentum q ∼ 100 MeV, andÃ is the closure energy, taken from Ref. [50] or estimated by the systematics, parametrizing the non-perturbative nature of QCD. Thus, the LECs play a role similar to that of the nuclear form factors arising in hadronization and their reliable determination, e.g. using lattice QCD input, is necessary to calculate the 0νββ decay rate in the chiral EFT framework. The benefit of this approach is that one can avoid the factorization of the nucleon currents, which is a necessary approximation in the hadronization procedure.

C. Determination of NMEs in the IBM-2
In order to evaluate the NMEs we make use of the microscopic interacting boson model (IBM-2) [56,57] which has the advantage that it can be used to all nuclei of interest.
The interacting boson model has been one of the most successful models in reproducing collective features of the low-lying levels of medium as well as heavy nuclei, and is one of the few models that can be used consistently to all nuclei of interest. We have already studied different mechanisms systematically using the microscopic interacting boson model (IBM-2) [23-25, 49, 58-60] and this study adds the short-range non-standard mechanisms of double beta decay to the list.
The method of evaluation is discussed in detail in [23,25]. Here we briefly mention the logic of the method, which is a mapping of the fermion operator H onto a boson space and its evaluation with bosonic wave functions. The mapping [61] can be done to leading order The values of parameters used in the current calculations are given in Appendix A.
The single-particle and single-hole energies and strengths of interaction were evaluated and discussed in detail in Ref. [62] where the occupancies of the single particle levels were calculated in order to satisfy a twofold goal: to asses the goodness of the single particle energies and check the reliability of the used wave functions. Both tests are particularly important in the case of nuclei involved in double beta decay, as they affect the evaluation of the NMEs and then their reliability [63]. The energies of the single particle levels constitute a very important input for the calculation of the occupancies in the method used in Ref. [62]. In principle those energies can be considered as input parameters that can be fitted to reproduce the experimental occupancies. Instead of fitting, the single particle energies were extracted from experimental data on nuclei with a particle more or one particle less than a shell closure.
These single particle energy sets were then used to calculate the occupancies of several nuclei of interest in double beta decay. Finally, the results were compared with other theoretical calculations and experimental occupancies, when available, and good correspondence was obtained. As part of the calculation single particle energies for several major shells were updated to values given in Appendix B.
Finally, an additional improvement is the introduction of short-range correlations in the nuclear structure calculation. These are of crucial importance for short-range non-standard mechanisms and they can be taken into account by multiplying the potential v(r) in coordinate space by a correlation function f (r) squared. The most commonly used correlation function is the Jastrow function, with a = 1.1 fm −2 , b = 0.68 fm −2 and c = 1 for the phenomenological Miller-Spencer parametrization [64], and a = 1.59 fm −2 , b = 1.45 fm −2 and c = 0.92 for the Argonne parametrization [65]. Since our formulation is in momentum space, we take short-range correlations into account by using the Fourier-Bessel transform of f J (r). By specifically separating the value of g A we allow for the possibility of a quenching of the axial-vector coupling. Even though quenching of g A goes beyond the topic of this study, we would like to remind that it is well known from single beta decay and electron capture  that g A is renormalized in models of nuclei. Quenching of g A in 2νββ-decay, consistent with single-beta decay, has also been observed [24,25] (for a review see [66]). However, the question of whether or not g A in 0νββ decay is renormalized as much as in 2νββ is of much debate. This problem is currently being addressed both experimentally, by employing single and double charge exchange reactions [67,68], and theoretically, by using effective field theories to estimate the effect of non-nucleonic degrees of freedom [69]. Quenching of g A arises from the omission of non-nucleonic degrees of freedom and from the limited model space in which the calculations are done. The former effect is not expected to be present in 0νββ decay since the average neutrino momentum is ∼ 100 MeV, while in 2νββ decay is of the order of 1 − 2 MeV. The latter effect instead appears both in 0νββ and 2νββ decays.
This consideration suggests to use an effective value of g eff A = 1.0, in between the free value g A = 1.269 and the value observed in 2νββ decay, g A ∼ 0.6. We henceforth use this value.

Comparison with Earlier Results
From the NMEs in Tables II and III one can calculate the NMEs for the standard mass mechanism, M ν and heavy neutrino exchange M ν h = M LL 3 to compare with earlier calculations. To this end, it is convenient to introduce the quantities  [25] using the quenched value g A = 1.0 and the convention that M ν < 0. The "old" F , GT and T NMEs of Table I in [25] are combined in the NMEs M old ν andM old ν using a negative and positive sign of the tensor NME relative to that of the GT NME, respectively. and write M ν as and similarly for M ν h .
The values of the NMEs in the present work are compared with those in [25] in Table IV for light neutrino exchange and in Table V   improved single particle energies is sizeable in 76 Ge, 82 Se, 96 Zr, 150 Nd and small otherwise.
The main difference between the calculation reported in [25] and the present one is in the sign of the tensor matrix element in Eq. (41). The present derivation gives a sign of the M T term relative to that of M F which is opposite to the one employed in Ref. [25]. This correction has little effect on the standard mass mechanism, for which M T is small, but has considerable effect on the short-range mechanisms. Additionally, one can see that the In Figs. 2 and 3 we, respectively, compare the compound NMEs M ν and M ν h in the different calculations: present work (red circles), Ref. [25] (filled blue squares) and Ref. [25] with the correct sign for the tensor term (empty blue squares). This allows disentangling the effect of the new single particle energies from that induced by the sign of the tensor NME.
As already mentioned, for light neutrino exchange (Fig. 2), the sign of the tensor term has relatively little impact, whereas the single particle energies lead to a sizeable increase of M ν in lighter isotopes. On the other hand, Fig. 3 demonstrates the strong effect of the sign of the tensor term in essentially all isotopes.
In comparison with calculations other than IBM-2 we note that our improved results for the standard mass mechanism are very similar to those in QRPA in all isotopes [26], but still differ from those of the Shell Model [70]. For the short-range mechanisms the obtained numbers are again similar to the QRPA in case in which both neutron and proton are particle-like (p-p) or hole-like (h-h), while differ in the case neutrons are hole-like and protons are particle-like or vice versa (p-h and h-p). We also note that, although not discussed here,  Table V) in reasonable agreement with the QRPA result, except for the overall sign in Eq. (41), which, as indicated above, is opposite to that of QRPA.

Compound NMEs for Short-range Mechanisms
Concluding our calculation of NMEs, in Tab Table IV.
We note that the NMEs M 1 and M 5 are generally enhanced due to the large pseudo-scalar charges g P and g P , though in M 5 this is often compensated by the suppressed component NMEs. The enhancement is especially strong in isotopes with the same sign for M P P GT and M P P T which arises from particle-particle versus particle-hole configurations of the nucleons. This, along with the large PSFs discussed below, makes 100 Mo an ideal isotope to probe the corresponding mechanisms from a theoretical point of view.

A. Leptonic Matrix Elements
Besides the NMEs, the calculation of the 0νββ decay rate requires the calculation of the so called leptonic phase space factors. Here, we follow the numerical approach of Ref. [72].
Because of Pauli-blocking of the inner states, the nucleons are expected to decay largely at the surface of the nucleus, which means that the electron wave function can be approximated by its value at the nuclear radius r = R A .
Since we are interested only in 0 + → 0 + 0νββ transitions, nucleon operators of a certain parity must be combined with partial leptonic wave functions of the same parity. Specifically, the parity-even operators will be accompanied by S 1/2 − S 1/2 and P 1/2 − P 1/2 electron wave functions, while parity-odd ones will go together with the S 1/2 − P 1/2 combination of wave functions. In this study we restrict ourselves only to the S 1/2 − S 1/2 approximation, which allows us to drop the parity-odd nucleon operators from our calculation. The leptonic squared matrix elements for S 1/2 − S 1/2 wave functions, summed over the electron spins s 1 and s 2 , then read [22] s 1 ,s 2 where the scalar product between the asymptotic electron momentum vectors is parametrized asp 1 ·p 2 = cos θ with the opening angle 0 ≤ θ ≤ π. The term (1 − P e 1 e 2 ) indicates that the matrix element is anti-symmetrized over the electrons. The result for Eq. (42) when both currents are left-handed is the same as the one shown when both currents are right-handed.
Since we are interested only in 0 + → 0 + transitions in the S 1/2 − S 1/2 approximation, we omit phase space factors corresponding to µ = j or ν = j. Further, in Eqs. (42) -(44) we have used the quantities f Here, the definitions in terms of electron wave functions g −1 (E) and f 1 (E) evaluated at the nuclear surface apply, . When compared to Refs. [11] and [18] our results agree but we also introduce additional factors f (0,1) 11− which appear as a result of the interference between the left-and right-handed scalar electron currents. In fact, these terms are not independent of the others as they can be expressed as f (0,1) 66± . In determining the squared leptonic matrix elements, we numerically calculate the electron wave functions according to [72], taking into account the finite nuclear size and electron cloud screening corrections.

B. Differential Decay Distributions
The NMEs presented in the previous section and the squared leptonic matrix elements shown in Eqs. (42) - (44) can now be combined to calculate the rate of 0 + → 0 + 0νββ decay.
The fully differential rate is expressed as [16][17][18] and where E 2 , p 1 = E 2 1 − m 2 e and p 2 = E 2 2 − m 2 e are understood to be functions of E 1 due to overall energy conservation, E 2 = Q ββ + 2m e − E 1 . Here, Q ββ is the so called double beta decay Q value of the given isotope, i.e. the kinetic energy release of the electrons.
The coefficients a(E 1 ) and b(E 1 ) in Eq. (48) are, respectively, given by a(E 1 ) = f and the energy-dependent angular correlation is introduced as α(E 1 ) = b(E 1 )/a(E 1 ). The latter has the property −1 < α(E 1 ) < +1 and as it appears in front of the cos θ term, it describes the likelihood for the electrons to be emitted back-to-back (α(E 1 ) −1), collinearly (α(E 1 ) +1) or isotropically (α(E 1 ) ≈ 0). Defining and their ratio K = B/A, the angular distribution reads With the given information we determine the single electron distribution dΓ/dE 1 and the angular correlation α(E 1 ) for the three relevant phase space factors that occur for shortrange operators: f also apply for the standard mass mechanism, calculated in the closure approximation.
The resulting single energy distribution and angular correlation were already presented in Ref. [22] for several isotopes, but in Fig. 4 (left) we illustrate the normalized single energy distributions for 76 Ge as functions of the kinetic energy E kin 1 = E 1 −m e of one of the electrons, i.e. the range is from zero up to Q ββ value. As can be seen, the term f (0) 11− produces an energy distribution virtually indistinguishable from that of f (0) 16 . All mechanisms produce a hill-like shaped energy distribution and observing the single energy spectrum is not expected to help distinguish between the standard mass mechanism (corresponding to f (0) 11+ ) and any of the short-range mechanisms. The angular correlation α(E kin 1 ), shown in Fig. 4 (right), can distinguish between different mechanisms, namely short-range mechanisms of type i = 4, 5 produce electrons that are emitted collinearly whereas for i = 1, 2, 3 and the standard mass mechanism, they are dominantly back-to-back. As mentioned, the factors f vanish. There is therefore no change of the angular correlation due to interference and the angular correlation is an incoherent sum over contributions.

C. Total Decay Rate
Finally, we can integrate over the whole electron phase space to determine the total decay rate Γ and the decay half life T 1/2 , To facilitate calculation of the total rate under the presence of one or more mechanisms, we define the integrated PSFs [72] G (0,1) ij = 2C ln 2 With the above PSFs, the inverse 0νββ decay half-life can be written

A. Bounds on the Effective Neutrino Mass
With ν = m ββ /m e and the other short-range I set to zero, Eq. (57) simplifies to the well know formula for light neutrino exchange,   Using the updated NME values for the light neutrino exchange mechanism shown in Tab. IV (last column) we can set new limits on the effective 0νββ mass |m ββ |. For isotopes with existing experimental bounds on the 0νββ decay rate, the resulting limits at 90% CL are summarized in Table VIII. As mentioned, the axial coupling is set to g A = 1.0. Generally, the limits have improved compared to the previous analysis [25]. This is a consequence of the better experimental limits for 76 Ge, 82 Se, 130 Te and 136 Xe as well as of the updated single particle energies in the NMEs for 76 Ge, 82 Se, 96 Zr and 150 Nd.
In Fig. 5, we compare the existing limit and future sensitivities in a plot correlating the 0νββ mass |m ββ | with the sum of neutrino masses Σm ν = m ν 1 + m ν 2 + m ν 3 for the standard picture of three active neutrinos. The shaded regions indicate, as usual, the allowed parameter space for normally (NO) and inversely (IO) ordered neutrino spectra by varying over the Majorana CP phases, where we take the best fit values of the oscillation angles and mass-squared differences as given in [81]. Using our NMEs, the currently best limit is set by the KamLAND-Zen collaboration T 1/2 ( 136 Xe) > 1.  |m ββ | < 114 meV at 90% CL for g A = 1.0. The recent final result from GERDA with T 1/2 ( 76 Ge) > 1.8 × 10 26 yr [9] corresponds to an essentially equal limit of |m ββ | < 118 meV at 90% CL. In Fig. 5 we also illustrate the corresponding limit assuming no quenching with g A = 1.27, giving |m ββ | < 76 meV. In addition to the current limit we also show two examples of prospective sensitivities T 1/2 ( 100 Mo) = 5 × 10 26 yr expected at AMoRE-II [82] and T 1/2 ( 76 Ge) = 10 28 yr for LEGEND-1000 [83]. The latter will probe the full IO regime and a large chunk of the NO regime.
Neutrino masses are also probed by the cosmological effect of the relic neutrino background on the cosmic microwave background and the structure of the universe. Current observations are compatible with no effect arising from neutrino masses setting stringent limits on Σm ν down to Σm ν < 150 meV at 90% CL [84]. The limit generally depends on the neutrino ordering due to different priors in the statistical analysis and it is affected by the choice of the astrophysical data. It can also be weakened if an underlying cosmological model other than the standard minimal ΛCDM is used. In Fig. 5 we show the most conservative limits arising from a choice of cosmological models surveyed in Ref. [84]. Namely, Σm ν < 280 meV (NO) and Σm ν < 290 meV (IO) at 95% CL arise in the ΛCDM with nonzero neutrino masses and a free scaling of the so-called weak lensing amplitude A lens (ΛCDM + Σm ν + A lens ). These limits correspond to |m ββ | < 89 meV (NO) and |m ββ | < 101 meV.

B. Bounds on Effective Short-Range Mechanisms
We can likewise assume that only a single short-range contribution is present by setting all other coefficients to zero and specifically assuming that the standard light neutrino contribution is negligible. Equation (57) then reduces to with the PSF G I and NME M I depending on the type of contribution. From the current nonobservation of 0νββ decay we can then set upper limits on the effective I couplings. These are also shown in Table VIII  matrix elements also plays an important role.
The limits in Table VIII on the effective couplings apply at the QCD scale Λ QCD ≈ 1 GeV.
As described in [22] following [85,86]    The effective short-range operator couplings can be interpreted in terms of effective New Physics operator scales Λ I where we simply match using the couplings c I defined at the electroweak scale. In Fig. 6

C. Interference between Light Neutrino Exchange and Short-Range Mechanisms
So far we have only considered one mechanism (operator) to be present at a given time, either the light neutrino exchange or one of the short-range operators. We now discuss the effect of two or more mechanisms operating at the same time. A large number of combinations are of course possible but at least the standard light neutrino contribution is expected to be present at some level in any case. This is because any New Physics scenario that generates a ∆L = 2 short-range operator is also expected to generate Majorana neutrino masses at a level to explain neutrino oscillations. Therefore, it is reasonable to look into the interference of one of the non-standard short-range mechanisms with the standard light neutrino exchange. We here discuss a few illustrative scenarios.
We first consider the interference with the operator associated with LLL

3
. , Eq. (57) simplifies to Because light neutrino exchange and the operator associated with LLL Here, A = G (0) 16 |M ν ||M RR 5 |/m e are positive coefficients (we consider the NMEs to be real with M ν , M RR 5 having opposite signs, cf. Tab. VI), and α, β are the complex phases of m ββ , RR 5 , respectively. We again consider that the relative phase between m ββ and RR 5 is α − β = 0, π in which case Eq. (62) is a quadratic function in |m ββ | and RR 5 and for a given value of T −1 1/2 represent an ellipse. This is shown in Fig. 8 (left) where the tilting is determined by the size of the PSF G  Table VI. Fig. 8 (right) shows the equivalent plot for the effective operator scale Λ RR 5 . As can be seen in Table VII, the PSFs G 16 applicable to all contributions of type 4,5 are generally quite sizeable resulting in a comparatively strong interference. On the other hand, the PSF G (0) 11− regulates the interference with operators of type R 1,2,3 with right-handed lepton currents, see Eq. (57), which is suppressed by the electron mass compared to the beta decay Q ββ value.
The above constraints on the effective neutrino mass and short-range operator couplings can be interpreted in terms of the New Physics scenarios introduced in Sec. II B.

Light and Heavy Sterile Neutrinos
In the sterile neutrino case discussed in Sec. II B 1, we consider the simplified scenario where a single sterile neutrino of mass m N with mixing V eN to the electron neutrino contributes to 0νββ decay. The limiting cases where the sterile neutrino is much lighter and heavier than 100 MeV were discussed in Sec. II B 1. Currently, the most stringent limit in with the average momentum transfer q 2 = m p m e |M LL 3 /M ν |. In Fig. 9, we show the current limit and future sensitivity in the (m N , |V eN | 2 ) parameter space. The region above the 0νββ bottom-most contours indicated are ruled out by the corresponding observation, assuming that the sterile neutrino is of a Majorana nature. We compare the 0νββ decay constraints with other searches for sterile neutrinos which are being pursued in neutrino oscillations, single beta decays, meson decays, at colliders and in electroweak precision measurements. The most recent searches are summarized in Ref. [91]. The shaded area is excluded by current data and the dashed lines give examples of sensitivities in future searches.
This includes the Tritium decay experiment KATRIN [92], searches for long-lived particles (LLP, the shape is mainly determined by the planned DUNE [93], SHiP [94] and FCC-ee collider [95]) and high energy colliders FCC-hh [96], ILC [97] and CLIC [98]. As can be seen, future 0νββ decay searches at a level of T 1/2 ≈ 10 28 yr will be able to probe mixing strengths expected for light neutrino neutrino mass generation via the Seesaw mechanism, We stress that this strong sensitivity of 0νββ decay searches applies to purely Majorana neutrinos, which are difficult to reconcile with the lightness of active neutrinos for m N 1 GeV. For sterile neutrinos with such masses it is more natural that they form Quasi-Dirac states where LNV is suppressed by a small mass splitting. In Fig. 9, we also show the sensitivity towards such a scenario where two Majorana neutrinos with a relative mass splitting of ∆m N /m N = 10 −4 form a Quasi-Dirac pair, partially cancelling their contributions to 0νββ decay. While the sensitivity is strongly reduced, 0νββ decay searches are still competitive at this level for m N ≈ 1 MeV and m N ≈ 100 GeV.

Left-Right Symmetry
In Fig. 10, we show the limits from 0νββ decay searches on the right-handed W R boson mass m W R and the heavy neutrino mass m N in the LRSM introduced in Sec. II B 2. Here, we consider a simplified scenario with one lepton generation, i.e. a single heavy neutrino N and U R e1 = 1. We also choose the so-called manifest left-right symmetric case with g R = g, cos θ R C = cos θ C and take m ∆ −− R = m W R for the mass of the doubly charged triplet Higgs.
The solid and dashed 0νββ curves give the lower limit on m W R where we additionally neglect the W boson mixing, sin θ W LR = 0, thus RRR 3 is the only contribution. The rise of the 0νββ curves to the right of m N ≈ 10 3 GeV in Fig. 10 (64)).
Both the LHC and SHiP limits were derived assuming negligible W boson mixing; those based on the lifetime of the heavy neutrino will be affected and need to be re-assessed.
We nevertheless also include the sensitivity of future 0νββ decay searches for sin θ W contribute. Future searches are then expected to be sensitive up to m W R ≈ 26 TeV.

R-Parity Violating Supersymmetry
Assuming gluino dominance, R-parity violating supersymmetry will induce the contributions in Eq. (18). Neglecting any other contributions, including those from light neutrinos, Eq. (57) simplifies to where the numerical factors in front of the NMEs take into account the effect of QCD run- Ref. [102] it is an effect of the QCD running and operator mixing.
If 0νββ decay is not observed in future experiments with a sensitivity approaching T 1/2 ( 100 Mo) = 10 27 yr, the limit will improve to This is mainly a result of the strong sensitivity to RR 1 especially in 100 Mo, see Sec. V B.
As mentioned, the derived limit is based on the assumption of gluino dominance. It will be important to re-evaluate the impact of 0νββ decay searches on the R-parity violating supersymmetry in light of the new results and the current constraints from direct searches for supersymmetric particles. We have presented a first complete numerical evaluation of the NMEs needed for the description of short-range non-standard mechanisms of 0νββ decay. The calculation is performed within the framework of IBM-2 with restoration of the isospin properties of the Fermi transition operator. We also use updated single particle energies extracted from experimental data on nuclei with one nucleon removed or added from shell closure. We include additional NMEs that become important when the latest values of the nucleon form factors are taken into account. However, the main difference between previous calculations and the present one is in the sign of the tensor NMEs. The present derivation gives a sign of the tensor term M T which is opposite to that in e.g. [25]. This change has little effect on the standard light neutrino mechanism, for which M T is small ≈ 1%, but it is sizeable for the short-range mechanisms.

VI. SUMMARY AND CONCLUSION
As noted, we have performed our calculation in the phenomenological framework of the interacting boson model, using nucleon currents in the impulse approximation including higher-order terms in the nucleon momentum transfer determined in [22]. We model pionmediated modes via enhanced pseudo-scalar nucleon form factors informed by PCAC and lattice QCD calculations. In our numerical results we consider a possible quenching of the axial-vector coupling by choosing g A = 1.0 compared to the unquenched value g A = 1.27.
We follow this classical approach in contrast to ab initio methods based on chiral EFT interactions [51]. Such formulations promise the determination of NMEs with controllable errors, e.g. may address part of the quenching problem [103]. Calculations of the standard light neutrino exchange NME M ν following this approach have become possible for the lightest double beta decay isotopes 48 Ca [104][105][106], 76 Ge and 82 Se [105], indicating noticeably smaller values than those from phenomenological models such as IBM-2, see [107] for a recent review. If confirmed, this will require an understanding for such a deviation as well further studies to apply ab initio methods to heavier nuclei. NMEs should ideally be verified experimentally by employing single and double charge exchange reactions [67,68]. Chiral EFT techniques have been used to reveal a potentially sizeable short-range contribution in standard light neutrino exchange [53] and to calculate exotic contributions [52,54].
In addition to the NMEs calculated in our approach we also present the full set of lep-tonic PSFs for all relevant isotopes, determined numerically including effects from the finite nuclear size and electron cloud screening corrections. This allows us to set updated limits on the effective couplings of all possible short-range operators contributing to 0νββ decay.
Considering one operator at a time, the current limits correspond to operator scales ranging between 3 TeV to 10 TeV, where the strongest sensitivity is achieved for operators enhanced by pion-mediated corrections, in agreement with previous analyses [54,[108][109][110], in our case arising from enhanced pseudo-scalar form factors. We further illustrate the interplay between different contributions by considering the interference between the standard light neutrino exchange with one short-range contribution I thus setting constraints on the com- A detailed description of the IBM-2 Hamiltonian is given in [57] and [112]. For most nuclei, the Hamiltonian parameters are taken from the literature [113][114][115][116][117][118][119][120][121][122][123][124][125][126]. The new calculations are done using the program NPBOS [112]. They include energies, B(E2) values, quadrupole moments, B(M1) values, magnetic moments, etc.. For the semi-magic nuclei 124 Sn and 136 Xe, we have obtained the parameters by a fit to the energy of the low lying states using the same procedure as in Ref. [121] for 116 Sn. A compilation of the used parameters is given in Table X. The reliability of single-particle and -hole energies as well as the interaction strengths in connection with IBM-2 wave functions was studied in [62] by comparing recently measured occupation probabilities of initial and final states of interest in double beta decay. The pair structure constants were generated as usual by diagonalizing the surface delta interaction (SDI) in the two identical particle states, pp, nn, where the strength of the (isovector) interaction, A 1 , is obtained by fitting the 2 + -0 + energy difference in nuclei with either two protons (proton holes) or two neutrons (neutron holes). The used single particle energies and A 1 values are given in Tables XI -XIV. 198 Hg [125]