Short-Range Neutrinoless Double Beta Decay Mechanisms

Neutrinoless double beta decay can significantly help to shed light on the issue of non-zero neutrino mass, as observation of this lepton number violating process would imply neutrinos are Majorana particles. However, the underlying interaction does not have to be as simple as the standard neutrino mass mechanism. The entire variety of neutrinoless double beta decay mechanisms can be approached effectively. In this work we focus on a theoretical description of short-range effective contributions to neutrinoless double beta decay, which are equivalent to 9-dimensional effective operators incorporating the appropriate field content. We give a detailed derivation of the nuclear matrix elements and phase space factors corresponding to individual terms of the effective Lagrangian. Using these, we provide general formulae for the neutrinoless double beta decay half-life and angular correlation of the outgoing electrons.


I. INTRODUCTION
While the Standard Model (SM) gauge group SU (3) C ×SU (2) L ×U (1) Y perfectly explains the interactions we observe, its breaking also provides masses to the charged fermions via the Higgs mechanism. The discovery of the Higgs at the Large Hadron Collider (LHC) [1,2] allows us to probe and test this mass mechanism in the SM. Yet, neutrinos continue to evade our understanding as only left-handed neutrinos exist in the SM and they therefore cannot acquire a so-called Dirac mass like the other SM fermions. Neutrino oscillation experiments [3] have unambiguously shown though that at least two of the three known neutrino species have finite masses and while they are not sensitive to the absolute neutrino masses, they point to mass scales of order 10 −2 eV to 5 × 10 −2 eV. In addition, cosmological observations set an upper limit on the sum of neutrino masses Σm ν 0.15 eV [4], assuming the standard cosmological model and with the exact value depending on the observational data considered.
Neutrinos could be of Dirac type as the other SM fermions, but this requires a new right-handed neutrino ν R and tiny Yukawa couplings 10 −12 , which is rather unnatural.
Because the right-handed neutrinos would be completely sterile with respect to the SM gauge interactions, it is on the other hand theoretically indicated that they acquire a Majorana mass M of the form Mν L Cν T L . It is generically expected to be of the order of a large new physics scale Λ NP ≈ M associated with the breaking of lepton number L symmetry. Via the Yukawa couplings between left and right-handed neutrinos, it will induce an effective dimension-5 operator, Λ −1 NP (LLHH) [5], where L and H represent the SU (2) L left-handed lepton and the Higgs doublets, respectively. After electroweak (EW) symmetry breaking, a small effective Majorana mass m ν ∼ m 2 EW /Λ NP is generated for the active neutrinos. This corresponds to the famous seesaw mechanism [6][7][8][9][10], with a scale Λ NP naturally of the order 10 14 GeV to explain the light neutrino masses m ν ≈ 0.1 eV.
While the most prominent scenario, the high-scale seesaw mechanism is not the only possibility to generate light neutrino masses; there are numerous other ways by incorporating lepton number violation at low scales in secluded sectors, at higher loop order and when allowing higher-dimensional effective interactions beyond the Weinberg operator. If the L breaking occurs closer to the EW scale, higher-dimensional L-breaking operators will be important for phenomenology, and specifically they will potentially induce neutrinoless double beta 0νββ decay. contribution, (c) 9-dim operator leading to short-range contribution. Adapted from [14].
The search for 0νββ decay is the most sensitive approach to probe Majorana neutrino masses. The experimentally most stringent lower limit on the decay half life T 1/2 is derived using the Xenon isotope 136 54 Xe, physics which can induce it. As hinted on above, other mechanisms of 0νββ decay where the LNV originates from LNV masses and couplings of new particles appearing in various possible extensions of the SM. The same couplings and states will also induce light neutrino masses due to the Schechter-Valle black box argument [11], but the resulting contribution will not be necessarily dominant. Instead, we consider the 0νββ decay rate by expressing high scale new physics contributions in terms of effective low-energy operators [12,13]. As a basis of our subsequent discussion, we provide a brief overview of the possible effective contact interactions at the Fermi scale m F ≈ 100 MeV at which 0νββ decay occurs. These are likewise triggered by effective SM invariant operators violating ∆L = 2 of dimension 5, 7, 9, 11, etc.. Fig. 1 schematically shows the contribution of such operators. They can in general be categorized in two main classes: (i) Long-range transitions proceeding through the exchange of a light neutrino. This includes the so-called standard neutrino mass mechanism via Majorana neutrinos in Fig. 1 (a) but also via exotic interactions incorporating rightchiral neutrinos Fig. 1 (b). In the SM with only left-handed neutrinos, these operators violate ∆L = 2 and they incorporate a helicity flip through inclusion of a Higgs field. (ii) Short-range transitions with no mediating particle lighter than ≈ 100 MeV. As contact interactions with six external fermions, they are of dimension 9 and higher odd dimensions.
The most prominent scenario where such an operator is generated is through the inclusion of heavy sterile neutrinos [15]. In the above classification we do not include the case where additional light states are either mediating the decay or are emitted in it (e.g. Majorons).
Probing exotic ∆L = 2 transitions is crucial for our understanding of how light neutrinos acquire their tiny masses. If exotic contributions were to be observed in upcoming experiments, it would indicate that the origin of light neutrino masses is around the TeV scale.
Likewise, the non-observation of 0νββ puts strong constraints on neutrino mass mechanisms close to the EW scale. It is not just neutrino physics that can be probed, though. Operators violating ∆L = 2, or the underlying physics responsible for them, can also erase an asymmetry between the number of leptons and anti-leptons throughout the thermal history of the early universe. Together with the sphaleron transitions in the SM violating the sum of total baryon and lepton number (B + L), this will also erase an asymmetry between baryons and anti-baryons. The rate of this washout can be related to the half life of 0νββ decay for a given operator. The observation of non-standard 0νββ decay mechanisms can thus generally falsify baryogenesis mechanisms operating at scales above the EW scale [14,16].
A similar argument applies to other process probing lepton number violation around the TeV scale, such as searches for same sign di-leptons at the LHC [17,18].
Before discussing the exotic short-range contributions of our interest, we remind the reader that the mass mechanism of 0νββ decay is sensitive to the effective Majorana neutrino where the sum is over all light neutrinos with masses m ν i , weighted by the square of the charged-current leptonic mixing matrix U . This quantity is equal to the (ee) entry of the Majorana neutrino mass matrix. The inverse 0νββ decay half life in a given isotope can then be expressed by where G ν is the phase space factor (PSF) and M ν the corresponding nuclear matrix element (NME) of the process. The normalization with respect to the electron mass m e yields a small dimensionless parameter ν = m ν /m e . The current experimental results lead to a limit m ν 0.06 − 0.17 eV [19], with an uncertainty due to the different NMEs in various nuclear structure models. Future experiments will probe m ν ≈ 0.02 eV, corresponding to the lowest value for an inverse hierarchy of the light neutrino states. A popular modification is through the inclusion of light sterile neutrinos with masses in the range from eV to MeV [20][21][22][23], in which case the half life is still given by Eq. (3), but with different masses m ν i and couplings U ei [24][25][26].
The NMEs of the nuclear 0νββ transitions are notoriously difficult to calculate and limits derived from 0νββ decay are affected for any contribution. Detailed treatments using different nuclear structure model approaches can for example be found in [27][28][29][30][31][32][33][34][35][36]. Despite tremendous efforts to improve the nuclear theory calculation, the latest matrix elements obtained using various approaches differ in many cases by factors of ∼(2−3). Experimentally, the most stringent bounds on 0νββ decay are currently from 76 Ge [37] and 136 Xe [19]. The results presented below are using the recent results in 76 Ge of T Ge 1/2 ≥ 5.3 × 10 25 y and in 136 Xe of T Xe 1/2 ≥ 1.07 × 10 26 y at 90% confidence level (CL). Planned future experiments searching for 0νββ decay are expected to reach sensitivities of the order of T 1/2 ≈ 10 27 y.
For example, the recent comparative analysis [38] quotes a discovery sensitivity at 3σ of T Xe 1/2 = 4.1 × 10 27 y for the planned nEXO experiment [39]. For more details on the effective 0νββ interaction, see for example the review [40] and references therein. General up-to-date reviews of 0νββ decay and associated physics can be found in [41], while a more specific recent review on 0νββ NMEs is available in [42].
Besides the light and heavy neutrino exchange, exotic long-range mechanisms have received most attention so far [43][44][45][46][47]. This is reasonable as the underlying SM operators already occur at dimension 7. It is important to note though that due to the helicity-flip involved, such operators are typically suppressed by the smallness of the light neutrino masses [48]. In our work, we instead focus on short-range mechanisms that often do not suffer such a strong suppression at a similar level. In both cases, the half life triggered by a single mechanism may be generically expressed similarly to Eq. (3), where G I is the nuclear PSF and M I the NME, both generally depending on the Lorentz structure of the effective operator in question. The coupling constant I parametrizes the underlying particle physics dynamics, e.g. the couplings to and the masses of the heavy states integrated out. NMEs and PSFs for dimension-6 operators were given in [47]. In this paper we present a detailed derivation of NMEs and PSFs for dimension-9 effective operators.
The paper is arranged as follows. After presenting the general effective Lagrangian at the quark level in Section II, we outline the calculation of the 0νββ differential decay rate in Section III. We then give the derivation of the NMEs in Section IV. Section V details the calculation of the leptonic PSFs. The results of these calculations are then combined in Section VI to give explicit expressions for the decay rate and angular correlations. Limits on the effective couplings I are also derived therein, assuming one contribution is different from zero at a time. Section VII contains some concluding remarks.

II. EFFECTIVE PARTICLE PHYSICS LAGRANGIAN
The contributions to 0νββ decay can be parametrized by effective operators of dimension 6 and 9 [12,13], corresponding to short-range and long-range interactions, respectively. The general Lagrangian of 0νββ decay consists of long-range and short-range parts, corresponding to point-like vertices at the Fermi scale ≈ 100 MeV.
In this work, we concentrate on the short-range contributions for which the general effective interaction Lagrangian schematically reads [13] where the sum and the place holders • indicate that the currents involved can have different chiralities and there is a separate effective coupling • i for each such combination. Specifically, the hadronic and leptonic currents in Eq. (5) can be with σ µν = i 2 [γ µ , γ ν ]. The fields u, d and e are 4-component Dirac spinor operators representing the up-quark, down-quark and electron, respectively. The field e c = Ce denotes the charge conjugate, corresponding to the fact that all lepton currents violate the electron lepton number by two units. While the currents involved in Eq. (5) can have different chiralities as denoted in Eq. (6), the results will not depend on many of the specific choices.
As convention, the Lagrangian Eq. (5) is normalized by the factor G 2 F /(2m p ) with the Fermi constant G F and the proton mass m p . As a result, the effective coupling constants i are dimensionless.
In the Lagrangian Eq. (5) one does not have to consider all possible combinations of the chiralities of the currents, as some of them are redundant or vanish. In order to prevent any confusion about the basis of low-energy dimension-9 operators we are considering, we spell these explicitly out in Table I Similarly, the Lagrangian Eq. The total number of these operators agrees with the result we obtained as a consistency check from a calculation using the Hilbert series method [49,50] and with the results in [51].
Despite comments in the latter reference we do not see the need of inclusion of operators containing quark bilinears transforming as colour octets. These can be shown to be related by Fierz transformation to the herein presented operators with tensor Lorentz structure.
For example, for operator O LLL 2 one can find the following Fierz identity If we further apply the "coloured Fierz identity" based on the well-known group-theoretical formula onto the first term on the right-hand side of Eq. (9), we get the following transformation of Hence, we see that if the operator containing colour octets (the first term on the righthand side of the above equation) is neglected, for the fact it does not contribute to 0νββ decay as argued in [52], then the operator O LLL , and thus redundant.
In a similar way all the operators with tensor quark bilinears can be traded for operators containing colour octets. If these are neglected, then only the operators consisting of (both left-and right-handed) scalar and vector colour-singlet quark bilinears are left.

III. NEUTRINOLESS DOUBLE BETA DECAY RATE
The differential rate of 0νββ decay can be written as where |R| 2 is the full matrix element of the 0νββ decay process summed over the spin Here, we neglect the recoil energy of the final nucleus which is of the order Because of the overall rotational invariance and energy conservation, the differential rate can be expressed in terms of the energy m e < E 1 < Q ββ + m e of one of the electrons and the angle cos θ =p 1 ·p 2 , 1 0 ≤ θ ≤ π, between the two electrons. The energy of the other electron is given by E 2 = Q ββ + 2m e − E 1 and the magnitudes of the electron 3-momenta are The full matrix element of the process can be formally expressed as Here, O + F e p 1 s 1 e p 2 s 2 denotes the final state composed of the 0 + daughter nuclear state and the two emitted electrons, and O + I the initial nuclear state. It is understood that the wave function of the two electrons in Eq. (13) is anti-symmetrized; the same holds for the wave function of the O + I and O + F state in terms of the nucleons. In Eq. (13), we allow the most general form of the quark level Lagrangian Eq. (5) for L SR , which we symbolically express as where the summation over K and Ξ collectively denotes the different electron-quark-quark current combinations jJJ (including different chiralities) and the Lorentz contractions, respectively.
The evaluation of the matrix element R is rather complicated since, in general, the leptonic part is nested with the hadronic part. For the long-range case a detailed calculation was given by Doi et al. [43,44] and Tomoda [53]. Since the hadronic part is the product of two currents, a sum over a set of intermediate states |N must be performed. This is a daunting task since for 0νββ decay all states up to an energy E ≈ 100 MeV contribute. It is therefore customary to treat the above summation in the closure approximation, i.e. sum over a complete set of states, This approximation is very well justified in the case of short-range operators as the intermediate transition occurs at very high energies |q| ≈ 100 MeV, corresponding to the inter-nucleon distance, compared to the nuclear transition itself at Q ββ ≈ 1 MeV.
Another problem is that the leptonic and hadronic parts are entangled. In order to disentangle them, an approximation is made wherein the electron wave functions are evaluated at the surface of the nucleus [43,54]. This approximation can be improved by using simplified nucleon wave functions and calculating the weighted average electron position [54], but the approximation we employ here does not introduce a sizeable error.
In this approach, the overall matrix element R is factorized into (i) the product of the leptonic matrix element that will be integrated over the two-electron phase space yielding the so-called phase space factor (PSF) that depends only on the leptonic current and the electron wave function at the surface of the nucleus 2 , and (ii) the nuclear matrix element (NME).
For the latter, we first reduce each nucleon current J K Ξ to its non-relativistic form by means of a Foldy-Wouthuysen (FW) transformation, then take the product of the two currents and evaluate the matrix elements of the corresponding two-body operator in the nuclear many-body wavefunctions. In the FW transformation we take terms to order |q|/m p . This is a good approximation, since the momentum transfer in the process is of order |q| ∼ 100 MeV, and therefore |q|/m p ∼ 0.1. In certain cases of enhanced form factors higher order terms in the products of hadronic currents are also taken into account.
Altogether, this yields the full matrix element The leptonic matrix elements will be evaluated using the appropriate electron wavefunctions in Section V. The nucleon matrix elements are evaluated as discussed above, including appropriate q 2 -dependent form factors, in Section IV A-B, and the nuclear matrix elements in Section IV C.
Putting together the PSFs and NMEs, one can write the fully differential rate for 0 + → 0 + 0νββ decay as [43][44][45] 3 with where E 2 , p 1 and p 2 are expressed as functions of E 1 .
Following the notation of [43][44][45]  The total decay rate Γ and the decay half life T 1/2 are then given by From Eq. (17), one can calculate the single electron energy distribution, and the energy-dependent angular correlation Introducing the integrated quantities and their ratio K = B/A, one obtains the angular distribution Although both the single electron energy distribution and the angular correlation require dedicated experimental setups, we calculate them nonetheless since they contain important information on the underlying mechanism.

A. From Quarks to Nucleons
We are interested in matrix elements induced by the quark bilinears appearing in the Lagrangian (5), i.e. by the left-and right-handed scalar, vector and tensor quark currents.
Considering the nucleon isodoublet N = p n , the nucleon matrix elements of these color singlet quark currents have according to [55] the following structure Here we define and τ + is isospin-raising operator, transforming a neutron into a proton. The matrix elements are in general functions of the neutron and proton momenta p n = p N and p p = p N , respectively, and the momentum transfer entering the form factors is defined as q = p p − p n .
In Eq. (26) we omit the induced scalar and axial-tensor terms as the corresponding currents can be safely neglected because they vanish in the isospin-symmetric limit [56]. Moreover, they are suppressed by a factor 1/m p and are not enhanced by a pion resonance.
We parametrize these except F P S (q 2 ) and F P (q 2 ) in the so-called dipole form, where the coupling constants g X give the value of the form factor at zero momentum transfer, For example, the vector form factor can be experimentally determined from the electromagnetic form factor and from the conserved vector current (CVC) hypothesis, This parametrization provides a good description of F V (q 2 ) in the range 0 ≤ |q| ≤ 200 MeV of interest in 0νββ decay. A better parametrization, important for large q 2 1 GeV 2 , is given in [57], but it is of no interest for the purposes of the present paper.
The induced form factor F W (q 2 ) can also be determined from experiment, since it is related to the Pauli form factor F 2 (q 2 ) [57] and to the isovector anomalous magnetic moment of the nucleon, The axial vector form factor can also be parametrized in dipole form and is obtained from experiment, The value of g A is determined in neutron decay [58] and m A is obtained from neutrino scattering [59].
The induced form factor F P (q 2 ) cannot be directly obtained from experiment. We use the parametrization suggested in [15], based on the partially conserved axial-vector current (PCAC) hypothesis, with the pion mass m π = 0.138 GeV.
From Eq. (33) we have g P ≡ F P (0) = 231. This formula is consistent with a recent analysis in chiral perturbation theory [60], which gives g P = 233 and with recent measurements in muon capture, which give F P (q 2 ) |q 2 | 2mp at |q| = 0.88m µ , where m µ = 0.105 GeV is the muon mass. The calculated value is 8.0, while the measured value is 8.06 ± 0.55 [61].
A considerable amount of attention has been devoted recently to the form factors F S (q 2 ) and F P S (q 2 ), in particular to the values at zero momentum transfer, g S = F S (0) and g P S = F P S (0). Quoted values are g S = 1.02 ± 0.11 and g P S = 349 ± 9 [62]. Not much is known about the q 2 -dependence; for the scalar form factor, which, in the Breit frame, is the Fourier transform of the matter distribution, a reasonable parametrization is in dipole form with g S = 1 and m S = m V = 0.84 GeV, The value of the pseudo-scalar form factor F P S (q 2 ) diverges at q 2 = 0 in the chiral limit and results of lattice calculations depend on the extrapolation procedure. We take The question of whether or not the value of g P S is enhanced as in Ref. [62] is beyond the scope of this paper. The parametrization (35) reduces to the simple monopole form 1/ (1 + q 2 /m 2 π ) used in chiral perturbation theory, but it includes the finite size of the nucleon.
No experimental information is available for the tensor form factors. Ref. [62] quotes a value of 0.987 ± 0.055 for F T 1 (0) ≡ g T 1 . An old calculation [55] estimates, from the MIT bag with g T 1 = 1 (and g T 2 = −3.30, g T 3 = 1.34 estimated from [55]). The two form factors F T 2 (q 2 ) and F T 3 (q 2 ) do not enter the results of this paper, but are quoted here for completeness.

B. Non-Relativistic Expansion
To obtain the nuclear matrix elements of interest, we have to calculate the non-relativistic expansion of the above nucleon matrix elements. This form is obtained by a Foldy-Wouthuysen transformation [63,64], which is an expansion in powers of the velocity v/c or equivalently in |p|/m p .
The resulting expressions are summarized in the following section where we use the spatial momentum difference q = p p − p n and momentum sum Q = p p + p n . Particular terms will be listed according to the order in |q|/m p , where we perform the expansion up to and including terms of order |q|/m p , except for terms incorporating F P (q 2 ) and F P S (q 2 ) which are enhanced as discussed above. In these cases we retain terms of order q 2 /m 2 p and even higher.

Scalar Bilinears
The non-relativistic expansion of the scalar and pseudo-scalar nucleon current corresponding to J S±P = (1 ± γ 5 ) can be written as Here, I denotes the 2×2 identity matrix and σ = (σ 1 , σ 2 , σ 2 ) T is the vector of Pauli matrices, both operating in the spin space of the nucleon.

Vector Bilinears
The vector currents corresponding to J µ V ±A = γ µ (1 ± γ 5 ) have four different components (vector, axial vector and induced pseudo-scalar and weak magnetism), which can be non-relativistically expanded as follows

Tensor Bilinears
The non-zero nuclear components corresponding to the tensor bilinears Terms containing the momentum sum Q are called recoil terms [45].
The nuclear currents can be obtained from Eqs. (37)-(40) by summing over all neutrons in the initial nucleus as where J a Ξ denotes any of the six nucleon currents. (For positron emission, τ + is replaced by τ − and the sum is over protons. The short-range 0νββ decay transition involves two such currents, each transforming one neutron into a proton. Having the five different terms in the effective Lagrangian Eq. (5) we thus need to evaluate five different products of nucleon currents in the non-relativistic expansion. Furthermore, we need to sum over all neutrons in the initial nucleus and the corresponding nuclear transition operators can be expressed as with the products of the relevant nucleon currents Π i Ξ , i = 1, . . . , 5 given in Appendix A. As before, Ξ, Ξ symbolically denote the Lorentz structure. In the case of the terms 1, 2 and 3, the nuclear transition operator has no free Lorentz indices whereas for terms 4 and 5, there remains one free index that is contracted with that of the electron current.

C. From Nucleons to the Nucleus
The final and most challenging step in the determination of the 0νββ NMEs concerns the calculation of the matrix elements at the nuclear level. This requires an understanding of nuclear structure and given the highly complex many-body problem, it is not solvable from first principle. Above, we have constructed the short-range nuclear 0νββ transition operators using the general nucleon operator with their q 2 and thus distance dependence parametrized by experimentally constrained form functions. We now define the nuclear matrix elements as where O + Using this notation the matrix elements for the five different short-range operators read In these expressions, we have kept all terms to order 1. In M 3 we have retained also the term proportional to (F V (q 2 ) + F W (q 2 )) 2 , which is a bit smaller, but may still represent an important contribution of the respective operator. We have also separated out from the form factors F X (q 2 ), the so-called charges, i.e. the values at q 2 = 0, F X (0) ≡ g X , the and treated separately the A, P and P S couplings which have a different q-dependence. If the axial-vector coupling A is present in the first or second power, the q-dependence reads h A (q 2 ) = 1 h AA (q 2 ) = 1 respectively.
Similarly, if a single power of pseudoscalar coupling in combination with the axial-vector coupling A or some other coupling X is present, then we havẽ h AP (q 2 ) = 1 h P (q 2 ) = 1 respectively, while in case of the second power of pseudoscalar coupling the q-dependence has the formh The Fermi, F , Gamow-Teller, GT , and tensor, T , matrix elements, appearing in Eqs.
(44)-(48), can be calculated in any nuclear structure model [28,30,65]. We follow in this article the formulation of [15] and [28] where the two-body transition operator H is constructed in momentum space as the product of the so-called neutrino potential, v(q), times the form factorsh(q 2 ). Since we consider short-range mechanisms with a δ-function in configuration space, the Fourier transform is a constant, and H the neutrino potential in momentum space is [15,28] v where we have used the standard normalization. Incidentally, for the long-range mechanism the neutrino potential is whereÃ is the closure energy. This formulation allows one therefore to calculate simultaneously all matrix elements, short-and long-range, by simply specifying the neutrino potential.
As a further aside, to do calculations in coordinate space, one simply takes the Fourier-Bessel transforms of the product of the neutrino potential, v, times the form factor,h, where λ = 0 for Fermi and Gamow-Teller contributions and λ = 2 for tensor contribution.
Finally, an additional improvement is the introduction of short-range correlations (SRC) in the nuclear structure calculation. These are of crucial importance for short-range nonstandard mechanisms. These 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 [66], and a = 1.59 fm −2 , b = 1.45 fm −2 and c = 0.92 for the Argonne parametrization [67]. Since the formulation described above is in momentum space, we take SRC's into account by using the Fourier-Bessel transform of f J (r). Introducing (where the placeholder • is used to note that the same redefinition is used for all the above defined types of q-dependencies) and the notation where H ab denotes any two-body operator, we can write the Fermi, M F , and Gamow-Teller, M GT , matrix elements appearing in Eqs. (44)- (48) as When the Gamow-Teller matrix element comes with one or two powers of the axial-vector coupling, we define In the third short-range operator, matrix elements M GT and M T appear, which are defined as Since q ∼ 100 MeV in 0νββ decay, these terms are suppressed by a factor of 0.01 relative to the standard terms M GT and M T . However, the enhancement of the corresponding form factor partly compensates for this fact; therefore, we include them. These matrix elements can be easily calculated since the neutrino potential is just a function of q 2 .
Similarly, the matrix elements M P GT , M P T and M AP GT , M AP T are given by respectively. Also these terms are smaller by a factor of 0.01 relative to the standard terms M GT and M T . Nonetheless, this suppression is compensated by the enhancement of the g P form factor.
Next, we have terms M P P GT , M P P T contributing to the first short-range operator, which These matrix elements are smaller by a factor of 0.01. However, if the pseudoscalar coupling g P S is larger by two orders of magnitude as claimed in [62], these terms become comparable to those with the Fermi and Gamow-Teller matrix elements M F and M GT , or even larger.
The matrix elements M P P GT and M P P T appearing in the third short-range operator can be written as Again, since q ∼ 100 MeV in 0νββ decay, these terms are smaller by a factor of 0.0001 relative to the standard terms M GT and M T . However, as before, this suppression is balanced by the enhancement of the form factor g P , which appears here in the second power. Also these terms can be easily calculated since the neutrino potential is just a function of q 4 .
TermsM AP GT andM AP T , also called recoil terms, are defined as 1 Although these terms are suppressed by a factor of 0.01, if the pseudoscalar coupling g P S is larger by two orders of magnitude, these terms become important. These matrix elements are difficult to calculate since the operator Q is not simply a function of q 2 . A good estimate can be however obtained by replacing (Q · q)/m 2 p by the expectation value in the state O + I , Q · q/m 2 p ∼ 0.01, in which case etc.
Finally, terms M q 0 P P GT , M q 0 P P T appearing in the fifth matrix element read Since q ∼ 100 MeV and q 0 ∼ 10 MeV in 0νββ decay, these terms are smaller by a factor of 0.0001 relative to the terms M GT and M T . However, this suppression is balanced by the enhancement of the form factors g P and g P S .
The above basic building blocks are written in terms of the Pauli spin operators σ, the nucleon momenta difference q and sum Q, and the direction unit vector between two nucleons,r ab = r ab /|r ab |. Furthermore, we use S ab = 3(σ a ·r ab )(σ b ·r ab ) − (σ a · σ b ).

V. LEPTONIC PHASE SPACE FACTORS
The leptonic phase-space factors describe the atomic part of the physics involved in 0νββ decay. They quantify the effect of the relativistic electrons emitted in the process.
where p is the asymptotic momentum of the electron at long distance and s denotes its spin projection. The S 1/2 and the P 1/2 waves on the right-hand side of the above expansion are respectively given by [45] e S 1/2 where g κ ( , r) and f κ ( , r) are the radial wavefunctions of the 'large' and 'small' components.
The electron energy at asymptotically large distance is E = p 2 + m 2 e and its spin state is described by the two-dimensional spinor χ s . As before, the hat denotes the unit vector along the respective direction and σ is the vector of Pauli matrices here operating on the electron spin space. The wavefunctions satisfy the asymptotic boundary condition where κ = ±(j + 1 2 ), l k = j ± 1 2 , y = αZE/p and ∆ c κ is a phase shift. Here, p = |p|, α is the fine structure constant, j is the angular momentum of the electron and Z stands for the atomic number of daughter nucleus. Inside the nucleus, the radial wavefunctions Eq. (81) can be expanded in r approximated by the leading terms as for S 1/2 and P 1/2 waves, respectively. Here A κ are normalization constants and R A is radius of the daughter nucleus. In limit Z F → 0 the radial wavefunctions acquire the form of spherical Bessel functions, while the normalization constants become A ±1 ≈ (E ∓ m e )/(2E). In our calculations, however, the above shown approximations are not employed, as we derive the phase-space factors using numerically calculated radial wavefunctions as described in [54].
Therein a numerical solution is performed using a piecewise exact power series expansion of the radial wave functions. On top of the Coulomb potential of the daughter nucleus (with charge Z F ), V (r) = −αZ F /r, nuclear size and electron cloud screening corrections are taken into account. As a result, the considered potential reads where ϕ(r) is the Thomas-Fermi function taking into account the electron screening. The non-trivial r-dependence of the above potential for r < R is a result of the finite nuclear size, when a uniform charge distribution in a sphere of radius R A = R 0 A 1/2 with R 0 = 1.2 fm is considered. In order to calculate the electron currents involved in 0νββ decay we, in principle, need to evaluate the wavefunctions at the position of the corresponding transition.
To be exact, this would require the wavefunction of the nucleon undergoing the respective decay, ideally from the nuclear structure model, or using a simplified harmonic oscillator wavefunctions [54].
Instead, we follow [54] and adopt the approximation of evaluating the electron wavefunction at the nuclear radius r = R A , This choice reflects the fact that nucleons largely decay at the surface of the nucleus due to Pauli-blocking of inner states.
For 0 + → 0 + transitions, even nucleon operators need to be combined with S 1/2 − S 1/2 and P 1/2 − P 1/2 electron wave functions, while odd operators need to be combined with S 1/2 − P 1/2 wave functions. The calculation of the leptonic squared matrix elements are outlined in Appendix B and the results for S 1/2 − S 1/2 wavefunctions are Here,p 1 ·p 2 = cos θ is the scalar product between the momentum vectors of the two electrons, yielding the opening angle 0 ≤ θ ≤ π. The quantities f We note that our results agree with those of Päs et al. [13] and Tomoda [45], except for the extra interference term f 11− in Eq. (89) between the left-and right-handed scalar electron currents, and the fact that these authors use the notation of Doi [43,44], while we have used that of Tomoda [45]. The phase space factors corresponding to µ = j or ν = j in Eqs. (87) and (88) are not shown, as their corresponding contributions to 0νββ decay do not trigger 0 + → 0 + transition (in the case of S 1/2 − S 1/2 approximation) we are interested in (although they are relevant when general 0 + → J + transitions are considered). All the above listed phase space factors are given in terms of the underlying energy-dependent wave functions of the two electrons

A. Decay Half-life and Angular Correlation
Combining the above results, the coefficients a(E 1 ) and b(E 1 ) in the fully differential rate Eq. (17) for 0 + → 0 + 0νββ decay are given by a(E 1 ) = 2f These matrix elements are given in dimensionless units, that is they are multiplied by the mass number dependent radius R A = R 0 A 1/3 of the nucleus where R 0 = 1.2 fm.
The PSFs are calculated using the formalism of [68]. Following [54], we define the phase space factors with g (0) 16 = 0. The factor 1/R 2 A has been introduced in Eq. (98) to conform with standard notation to compensate the corresponding factor in the NMEs as discussed above. The numerical values of the phase space factors ij , adopted from [68], are given in Tab. III. jk for selected nuclei, adopted from [68], and corresponding angular correlation factors K 11 and K 66 (K 16 = 0 is always identically zero). All PSFs are given in units of 10 −15 y −1 as indicated.
The 0νββ decay half-life is then given by From the phase space factors one can also calculate the integrated angular correlation factor as in the three different cases jk = 11, 66, 16. The numerical values of K 11 and K 66 are also listed in Tab. III (K 16 = 0 is always identically zero). In view of the opposite sign for 11 and 66 in Tab. III, an eventual measurement of the angular correlation will allow a discrimination of the two types of non-standard mechanisms.

B. Bounds on Couplings
In principle, a given underlying particle physics may give rise to several contributions, and/or mixing among the corresponding Wilson coefficients will induce contributions through radiative effects from the scale of new physics, through the electroweak scale and down to the scale of QCD. The above formulae for the decay rate take into account all In constructing these tables, we have used the values of M F , M GT and M T given in Table   II and values of other matrix elements M , M andM estimated by replacing (q/m p ) 2 by the average value 0.01.

VII. SUMMARY AND CONCLUSION
In this article we have developed a general formalism for short-range mechanisms contributing to neutrinoless double beta decay in an effective operator approach. Such contributions will arise when lepton number is broken at a new physics scale Λ NP much larger than the typical energy scale of 0νββ decay q ≈ 100 MeV. We have calculated the expected 0νββ half lives by making use of the phase space factors of [54,68] and nuclear matrix elements to leading order in q/m p of [29], where we especially elucidate the different con-  contributions to short-range 0νββ decay and to verify the strong limits we obtain on the effective new physics parameters ranging between I ≈ 10 −10 to 10 −7 , which correspond to new physics scales in the multi-TeV region. Short-range contributions scale as ∝ 1/Λ 5 NP , and thus an increase in sensitivity on by an order of magnitude will only improve a limit on Λ NP by a factor of ≈ 1.6. Nevertheless, accurate calculations of the limits and sensitivities are crucially important as they probe the phenomenologically interesting TeV scale. A robust map of potential sources of lepton number violation in this energy region will help us to get a better understanding of the mechanism of neutrino mass generation.
Note: While finalizing this manuscript, we noticed the preprint [72], which discusses short-range contributions as well, but using a complementary approach based on chiral effective field theory. where the term proportional to F 2 P S (q 2 ) can be re-coupled using the following relation with S ab = 3(σ a ·r ab )(σ b ·r ab ) − (σ a · σ b ).
Term 2: J µν J µν j For the second term of the short-range part of the Lagrangian we get the following approximation of the nuclear currents Term 3: J µ J µ j Approximating the nuclear currents for the third term we obtain where the term proportional to (F V (q 2 ) + F W (q 2 )) 2 can be re-coupled as follows (σ a × q)(σ b × q) = − 1 3 (σ a · σ b )q 2 − 1 6 q 2 − 1 3 (q ·r ab ) 2 S ab , and as before S ab = 3(σ a ·r ab )(σ b ·r ab ) − (σ a · σ b ).
[O(10) S − P ] ē 1 (1 ± γ 5 )e c 2 ≈ (ē p 1 s ) S 1/2 (1 ± γ 5 )(e C p 2 s ) S 1/2 = e S 1/2 p 1 s † γ 0 (1 ± γ 5 )iγ 2 e where all the matrices are considered in the standard Dirac representation In case of the interference term combining the left-handed electron scalar current with the right-handed one, the calculation is analogous to the procedure shown above, only the plus sign in the definition of f 11 would change to a minus sign. b. Term 4, 5 The electron current is j µ ≡ē 1 (x)γ µ γ 5 e C 2 (x).