Neutrino Self-Interactions and Double Beta Decay

Neutrino Self-Interactions ($\nu$SI) beyond the Standard Model are an attractive possibility to soften cosmological constraints on neutrino properties and also to explain the tension in late and early time measurements of the Hubble expansion rate. The required strength of $\nu$SI to explain the $4\sigma$ Hubble tension is in terms of a point-like effective four-fermion coupling that can be as high as $10^9\, G_F$, where $G_F$ is the Fermi constant. In this work, we show that such strong $\nu$SI can cause significant effects in two-neutrino double beta decay, leading to an observable enhancement of decay rates and to spectrum distortions. We analyze self-interactions via an effective operator as well as when mediated by a light scalar. Data from observed two-neutrino double beta decay is used to constrain $\nu$SI, which rules out the regime around $10^9\, G_F$.


I. INTRODUCTION
The discrepancy between Cosmic Microwave Background (CMB) and local measurements of the Hubble constant, known as the Hubble tension, has grown to about 4σ [1][2][3][4][5]. If indeed a physical fact, it would imply that nonstandard particle physics or cosmology is required. Introducing a neutrino self-interaction (νSI), i.e. a fourneutrino contact interaction, to inhibit neutrino freestreaming in the early Universe can resolve the Hubble tension. The required strength of νSI needs to be much larger than the Fermi effective interactions predicted in the Standard Model (SM) [6][7][8][9][10]. Writing the interaction as G S (νν) (νν), 1 there are two regimes for the coupling G S : a strongly interacting regime with G S = 3.83 +1. 22 −0.54 × 10 9 G F and a moderately interacting regime with 1.3 × 10 6 < G S /G F < 1.1 × 10 8 [9].
In this work, we propose to search for νSI in double beta decay experiments 2 . These experiments search for the lepton number violating, and thus SM-forbidden, neutrinoless double beta decay (0νββ) [30,31]. The standard diagram of this process is the exchange of a massive Majorana neutrino, see Fig. 1 (left). As part of this effort, the SM-allowed two-neutrino double beta (2νββ) decay is measured with increasing precision and may itself be used 1 Here we adopt the Weyl spinor notation, with ν being a twocomponent spinor and the combination (νν) is a scalar product. 2 For future prospects in beta decay experiments, see Ref. [29]. to probe physics beyond the Standard Model [32]. In the presence of νSI, two neutrinos can be emitted via the effective νSI operator, see Fig. 1 (right). The final state of this νSI-induced double beta (2ν SI ββ) decay is identical to that of 2νββ decay. We will here discuss the total decay rate of 2ν SI ββ decay and the energy as well as angular distributions of the emitted electrons. Such a study is warranted because a dimensional analysis estimate reveals that the total decay rates of 2νββ and 2ν SI ββ are Γ 2ν ∼ G 4 F (0.1p F ) −2 Q 11 and Γ νSI ∼ G 2 S G 4 F p 2 F Q 11 , respectively. Here, the Fermi momentum p F ≈ 100 MeV represents the nuclear scale and Q ≈ (1 − 4) MeV is the isotope-dependent kinetic energy release (Q-value) in double beta decay. The factor 0.1 in Γ 2ν takes into account that without a virtual neutrino line only states up to 10 MeV are excited in 2νββ decay. Assuming that 2ν SI ββ rates of order Γ νSI > ∼ Γ 2ν can be seen experimentally, it is expected that couplings G S > ∼ 10 p −2 F ≈ 10 8 G F can be probed.
In previous studies of double beta decay, indirect constraints on νSI mediated by light scalars were obtained [33][34][35][36][37][38][39][40]. It was assumed that the scalar is emitted in the decay, hence is lighter than the Q-value of double beta decay. For a scalar particle φ that couples with strength g φ to two electron neutrinos, one finds from searches for so-called Majoron emitting double beta decays that g φ < ∼ 10 −4 − 10 −5 [35,36]. Taking m φ = 1 MeV, this bound on the Yukawa coupling corresponds to G S < ∼ (10−10 3 ) G F .
If νSI are not mediated by light scalars or the scalar mass is larger than the Q-value, this bound does not apply. In this case, the effect of νSI operators on double beta decay becomes more important, which we will investigate here.
In the next section we will study 2ν SI ββ decay in the effective operator language. In Sec. III we will generate the operator with an s-channel mediator whose mass is larger than the Q-value, and show how the distributions are affected. We conclude in Sec. IV, and various technical details are delegated to the Appendix.

II. νSI-INDUCED DOUBLE BETA DECAY
In the standard 0νββ mechanism, two neutrinos produced in double beta decay annihilate due to a Majorana mass term, leaving only electrons in the leptonic final states, as shown in Fig. 1 (left). Under the presence of νSI operators, where α, β = e, µ, τ are flavor indices, Fig. 1 (right) implies that the two electron-antineutrinos (ν e ) generated by neutron decay can take part in the νSI interaction and convert in the 2ν SI ββ process to ν α ν β or ν α ν β . Note that both the lepton number conserving (LNC) and violating (LNV) interactions in Eq. (1) can lead to 2ν SI ββ.
Assuming that the momenta of leptonic final states are negligible compared to the momenta of the neutrino propagators (the typical values of the former and the latter are of order Q = O(1) MeV and p F = O(100) MeV, respectively), it can be shown that the two processes in Fig. 1 share the same nuclear matrix elements (NMEs), see Appendix A. Consequently, we can compute the decay rate of 2ν SI ββ using the NME of 0νββ: Here m e denotes the electron mass and R = 1.2A 1/3 fm is the radius of the nucleus with nucleon number A. The structure of the 0νββ NME M 0ν is explained in Appendix A 2. The quantity G νSI is the 2ν SI ββ phase space factor, which is derived in Appendix A 1. It reads where p 1 and p 2 are the momenta of the two electrons. Neglecting the final state lepton momenta in the calculation of the 2νββ and 2ν SI ββ NMEs, the phase space factors are related as G νSI = G 2ν /(4π) 2 . The Q-value is given in Tab. I for various isotopes and F 2 (p 1 , p 2 ) stands for the Fermi function correction caused by the Coulomb potential of the nucleus. Finally, T 12 = E 1 + E 2 − 2m e is the total kinetic energy of both electrons, implicitly depending on p 1 and p 2 , and neutrino masses in the final state have been neglected. The constant c νSI appearing in the above equation reads where θ C denotes the Cabibbo angle. Note that the electron mass m e and nuclear radius R are included in Eq. (2) so that the normalization of the NME and phase space factor conforms with that adopted in the literature. Using Eqs. (2) and (3), it is straightforward to compute Γ νSI . It should be noted, however, that the electron spectrum of 2ν SI ββ is very similar to that of 2νββ decay. We will comment below on the potential differences arising in the case of a light s-channel scalar mediator inducing the νSI.
If the energy and angular resolution of detectors cannot distinguish the electron spectrum of 2ν SI ββ decay from that of 2νββ decay, then only the change of the total decay rate can be probed. The total rate Γ 2ν of 2νββ decay has been measured precisely for many isotopes. For example, the 2νββ rate of 136 Xe has been measured to a 3% level [48]. Nonetheless, there remains a considerable uncertainty in the theoretical prediction of the 2νββ decay rate arising from the NMEs. Writing the theoretical prediction for the total decay rate approximately as the sensitivity to G S will largely depend on the uncertainty of the NME ratio |M 0ν | 2 /|M 2ν | 2 . While some of the nuclear uncertainties are expected to drop out from this, unresolved issues such as the quenching of the effective nuclear axial coupling g A in 0νββ decay likely provide a major error source. Note that in the above equation, we neglect the effect of interference between the 2νββ and 2ν SI ββ diagrams. If two electron antineutrinos are being emitted in 2ν SI ββ decay, such an interference will generally take place. We proceed by constraining the new interaction requiring that the 2ν SI ββ rate is less than the one for 2νββ, where Γ ex 2ν stands for the experimentally measured value of 2νββ. This roughly corresponds to an assumed uncertainty in the NME ratio |M 0ν | 2 /|M 2ν | 2 within a factor of two. If this uncertainty can be reduced in future theory NME determinations, the sensitivity on G S will improve accordingly. As mentioned, the uncertainty depends on g A . In our calculations we implicitly assume the unquenched value g A = 1.269 as used in the calculation of the NMEs.
Taking the best-fit value of G S = 3.83 × 10 9 G F of the strongly interacting regime and the Interacting Boson Model (IBM-2) NMEs [42], we compute the decay rate Γ νSI and compare it with Γ ex 2ν in Tab. I. By requiring Γ νSI /Γ ex 2ν < 1, we obtain the corresponding constraints on Table I. Estimate of 2νSIββ decay rates for several isotopes. Here, Q is the corresponding Q-value, T 2ν 1/2 represents the experimental 2νββ decay half-lives adopted from Ref. [41] that can be translated to the experimental 2νββ decay rates using Γ ex 2ν = log 2/T 2ν 1/2 and ΓνSI denotes the theoretical prediction for the 2νSIββ decay rates computed from Eq. (2), assuming GS = 3.83 × 10 9 GF . Bounds on GS obtained according to Eq. (6) are presented in the last row. Nuclear matrix element values from IBM-2 [42]  νSI favored by cosmological data Figure 2. Upper limit on the νSI coupling GS from 2νββ decay data for several isotopes and three different NME calculations as indicated. The blue band corresponds to the strongly interacting regime GS = 3.83 +1.22 −0.54 × 10 9 GF favored by cosmological data, which here is excluded by observations of 2νββ of various isotopes. G S , which is presented in Fig. 2. Here, we also use NME values computed in the Interacting Shell Model (Shell) [49] and Quasi-particle Random Phase Approximation (QRPA) model [50]; the corresponding limits on G S are shown in Fig. 2, indicating the uncertainty arising from nuclear theory uncertainties. As one can see, the strongly interacting regime for G S favored by the cosmological data causes Γ νSI /Γ ex 2ν > 1 for all the isotopes listed in Tab. I. For some isotopes, Γ νSI can be even one or two orders of magnitude higher than Γ ex 2ν . Even including the theoretical NME uncertainties, most isotopes can fully exclude the cosmologically favored strongly interacting regime band, given the premise that two ν e are involved in the νSI.

III. ENERGY AND ANGULAR DISTRIBUTIONS
We now consider possible distortions of the electron energy and angular distributions arising from the νSI-induced contribution. For an exact contact interaction of four neutrinos and neglecting final state lepton momenta, one can show that the electron spectra of 2ν SI ββ decay are the same as that of 2νββ decay, see Appendix A 1. However, considering that the νSI operator may be generated by light mediators, the corresponding energy dependence of G S , can cause observable spectral distortions, as we shall discuss below. We should mention here that the spectral distortions depend on the underlying model for νSI, which currently still lacks comprehensive exploration. There are various possibilities to generate the νSI effective operator, as shown in Fig. 3, where both tree and one-loop level diagrams are illustrated. At tree level, the νSI operator can be opened via either an s-channel (diagram I) or a t-channel (diagram II) scalar mediator. For vector mediators, most of the discussions below apply as well 3 . For the sand t-channel diagrams, G S has the following energy dependence where m φ is the mediator mass and s ≡ p 2 , t ≡ −q 2 with p and q being the momenta of the mediators in the tree diagrams. In the context of 2ν SI ββ, they are of order t ∼ p 2 F and s < ∼ Q 2 , respectively. The values of G S at zero momentum transfer are denoted as G 0 S = g 2 φ /m 2 φ , with the coupling g φ between φ and the neutrinos. Note that in Eq. (7) we omit the small effect of the φ decay width.
At the one-loop level, the νSI operators can be generated e.g. by the box diagram in Fig. 3. The corresponding energy dependence of G S is much more complicated than in the tree-level case. In general, it depends on both s = p 2 and t = −q 2 . However, for most loop diagrams, there are no simple analytical expressions similar to Eqs. (7) and (8). For the box diagram, we can obtain a simple result assuming all the particles running in the loop have the same mass m φ and that s t m 2 φ . With these assumptions, following the calculation in Ref. [52], we get Compared to Eq. (8), where the expansion in t yields , the t/m 2 φ term in Eq. (9) has a different coefficient but the same sign. In addition to the box diagram, other one-loop diagrams are also possible, as illustrated by diagrams IV and V in Fig. 3. Such diagrams can be roughly regarded as treelevel diagrams with energy-dependent couplings or mediator masses, which may cause more elusive effects in probing νSI in experiments of different energy scales. Here we only mention these possibilities and refrain from further discussions 4 .
Among the aforementioned possibilities, only the schannel case in Eq. (7) can be analyzed without involving novel nuclear physics calculations. Other t-dependent scenarios necessarily involve integrals over q that are different from the one in 0νββ decay, which calls for a dedicated study in the future. Here we proceed only with the s-channel case, specifically for m φ > ∼ Q. While the t-channel may also contribute in this regime, its effect is expected to be considerably smaller due to the large t suppression in the propagator. For G S in Eq. (7), we have derived the 2ν SI ββ differential decay rate in Appendix A 1 yielding the dependence dΓ νSI dp 1 dp 2 dcos θ 12 4 We note that νSI may also lead to significant corrections to the neutrino self-energy, which is not fully identical to the neutrino mass in 0νββ [53]. The effect is quite model-dependent, and can be studied if a complete model of νSI has been constructed. Here, cos θ 12 = p 1 · p 2 /(p 1 p 2 ) with the angle 0 ≤ θ 12 ≤ π between the two emitted electrons and β i = p i /E i are the electron velocities. The effect of the s-channel mediating scalar is captured by the function where ξ = 2 arcsin((Q − T 12 )/m φ ). It is a function of the total electron kinetic energy T 12 , and as we will see it can cause distortions of both the energy and angular distributions of the electrons. In the limit m φ → ∞, the effective operator is recovered and the dependence approaches I s (T 12 ) ∝ (Q − T 12 ) 5 yielding a phase space factor equivalent to 2νββ decay. We first consider the energy distribution. All modern double beta decay experiments measure the differential decay rate dΓ νSI /dT with respect to the total electron kinetic energy. This rate is computed by integrating over cos θ 12 , p 1 and p 2 with the total kinetic energy T = T 12 fixed at a given value. As noted, in the limit m φ → ∞, dΓ νSI /dT will have the same energy distribution as that of 2νββ decay.
In Fig. 4, we show the electron energy distributions of 2ν SI ββ and 2νββ decay for the isotope 100 Mo with an s-channel mediator mass m φ = Q + 0.1m e , slightly above the kinematic threshold. For comparison, we also show a vertical line corresponding to 0νββ decay, and the distribution for Majoron emission (0νββφ) taken from Ref. [39]. As can be seen in Fig. 4, the energy spectrum of 2ν SI ββ decay is shifted towards lower energies when compared to the 2νββ spectrum. The shift can be understood qualitatively. With increasing T the energy taken away by neutrinos is smaller, leading to a smaller value of s and hence a smaller value of the schannel G S . To determine the experimental sensitivity to such distortion, we have performed a simple χ 2 -fit to the NEMO-3 data [54] as detailed in Appendix C. We find that for Γ νSI = r 2 Γ 2ν with r = 16%, the χ 2 -value is changed by ∆χ 2 = 9 (3σ), which implies that if the spectral distortion is taken into account, the bound on G S can be approximately improved by one order of magnitude. We emphasize that this applies for the specific mediator mass m φ = Q + 0.1m e and the sensitivity will decrease for larger masses.
In addition to the energy distribution, the angular distribution can also be measured in dedicated experiments such as NEMO-3 [54] and SuperNEMO [55]. From Eq. (10), the angular distribution of 2ν SI ββ decay is obtained by integrating out p 1 and p 2 . The result takes the form where the angular correlation k νSI θ is computed in Appendix A 1. For 2νββ and 0νββ decay, the electron angular distributions take the same form as Eq. (12), but with different angular coefficients k θ [55,56]. We refer to those as k 0ν θ and k 2ν θ , respectively. Their expressions are given in Appendix A 1 as well. For 100 Mo, the numerical values are k νSI θ = 0.58 (using again m φ = Q + 0.1m e ), k 2ν θ = 0.65 and k 0ν θ = 0.88. With these values, we show in Fig. 5 the angular distributions of electrons for the three processes. Again, we have performed a χ 2 -fit to the NEMO-3 data [54] and the result is r < 29% at 3σ confidence level. This indicates that the angular distribution is less sensitive than the energy distribution to distortions from 2ν SI ββ. This is in fact interesting, as among the running and future experiments only one (Su-perNEMO) has sensitivity on the angular distribution.

IV. CONCLUSION AND DISCUSSION
The search for 0νββ decay constitutes one of the most important aspects to determine the nature and properties of neutrinos. As we have demonstrated in Figs. 2, 4, and 5, in the presence of νSI involving two ν e , there can be significant effects not only on the total rates of 2νββ decay, but also on the spectrum shapes. If only the total rates are considered, we find that 136 Xe currently has the best sensitivity to νSI. The observed 2νββ rate implies G S < (0.32 − 0.43) × 10 9 G F , which is significantly lower than the cosmologically favoured value G S = 3.83 +1. 22 −0.54 × 10 9 G F in the strongly interacting regime. However, one should note that this bound does not apply if only ν µ and ν τ participate in νSI.
Including spectral distortions could further improve the sensitivities. This is of interest when the particle that mediates the self-interactions has a mass that is larger than the available Q-value of the double beta decay. The distortions are caused by the energy dependence of the effective coupling G S and hence are affected by the underlying models for νSI. In this work, we only consider an s-channel mediating scalar, which allows us to evade nuclear physics calculations and to quantitatively show spectral distortions of the energy and angular distributions. For other possibilities containing a t-channel dependence, very different spectral distortions could appear, which will be addressed when a more dedicated study involving nuclear physics calculations is carried out.
In our calculations, we neglected the interference between the SM 2νββ and exotic 2ν SI ββ contributions, which in principle would be present if two electron antineutrinos are emitted in the νSI-mediated process. We estimate that the interference contributes in Appendix B. When the theoretical determination of the 2νββ rate becomes more precise, it will be important to include this interference term, but currently it does not improve the sensitivity to G S .
In summary, our work shows that strong νSI favored by the cosmological data might have an impact on 2νββ decay experiments. Precision measurements of 2νββ decay spectra combined with more theoretical effort in computing NMEs have the potential of probing hidden interactions of neutrinos. Furthermore, it demonstrates the importance of having access to energy and angular distributions of electrons in double beta decay experiments.

Leptonic Part
The leptonic part of 0νββ matrix element reads where m ee = U 2 ei m νi is the usual effective mass with the neutrino masses m νi and the charged-current leptonic mixing matrix U . The above expression is calculated using two massless neutrino propagators and one mass insertion. If the whole neutrino line in 0νββ decay was considered as a single propagator, then with the P L projectors one would obtain which is approximately equivalent for |q 2 | m 2 νi . By comparing the two diagrams in Fig. 1 in the main text, the leptonic matrix element of 2ν SI ββ decay can be written in an analogous way. Instead of a mass insertion one employs the νSI vertex, which additionally gives two extra external neutrino legs: Here, ψ ν3 and ψ ν4 denote the external lines of neutrinos in the 2ν SI ββ diagram and p = p ν3 + p ν4 is the sum of final state neutrino momenta. Note that we are here assuming a lepton number violating νSI interaction; the result for the conserving case is the same. Assuming that the momenta of the final state leptons are negligible compared to the momentum q of the neutrino propagators, s = p 2 q 2 , one can immediately relate the amplitude of 2ν SI ββ decay to that of 0νββ decay, The above leptonic matrix elements are to be contracted with their nuclear counterparts. The structure of the latter is fully identical between the 2ν SI ββ and 0νββ cases.
With the same q-dependence giving rise to the same neutrino potential, the resulting NMEs are the same, In calculating the leptonic phase space factor, we will take the S 1/2 approximation for the outgoing electrons and neutrinos, where p and E e,ν denote the electron and neutrino 3momenta and energies, and u e,ν stands for the usual Dirac spinor. For the electron wave function we include the Fermi function F 0 (Z f , p e ), taking into account the interaction of the emitted electron with the final nucleus of charge Z f = Z + 2. It can be approximated for the purposes of our numerical calculations as [56] F 0 (Z f , p e ) = 4 (2pR) 2(γ0−1) with γ 0 = 1 − (Z f α) 2 and y = αZ f E e /p, where α denotes the fine-structure constant, R ≈ 1.2A 1/3 fm is the nuclear radius (A denotes the atomic number of the decaying isotope) and Γ(x) is the Gamma function.
For a 0 + → 0 + nuclear transition and the S 1/2 approximation of the wave functions of the emitted electrons with momenta p 1 and p 2 we therefore have for the 0νββ matrix element Here, M 0ν denotes the nuclear part of the full matrix element and F 2 (p 1 , p 2 ) = F 0 (Z f , p 1 )F 0 (Z f , p 2 ). Consequently, combining Eq. (A4) with the above leads to the 2ν SI ββ matrix element In general, if the two outgoing neutrinos are replaced by any two massless fermions with a scalar product connected to the G S vertex, one would always get the product (p 3 · p 4 ) in Eq. (A9).
With the above we can express the decay widths as where the factor m 2 e /(4R 2 ) is included to make the NME M 0ν dimensionless and the phase space G νSI have units of yr −1 , and to conform to the usual conventions employed in 0νββ decay calculations. Furthermore, G 0ν and G νSI are the phase space factors, which we can derive starting with the following integrals where E I and E F denote the energies of the initial and final nuclei, respectively, with the Q-value defined as Q = E I − E F − 2m e . Next, we transform the phase space integral from Cartesian to polar coordinates, where in the second row we use p i to denote |p i | for simplicity, and c ij to denote the cosine of the angle between p i and p j . With the above replacement and a similar one for i = 3 and 4, the integrals in Eqs. (A12) and (A13) become and Since the two final state neutrinos in 2ν SI ββ decay are not visible, we need to integrate over their kinematic parameters p 3 , p 4 and c 43 . We include the s-channel dependence of G S given in Eq. (7) in the main text and evaluate the following part of the phase space integral, Assuming the two outgoing neutrinos are massless, we have , and thus Integrating over c 43 , E 3 and E 4 sequentially gives with Eq. (A16) then simplifies to × p 2 1 p 2 2 dp 1 dp 2 dc 21 (2π) 4 I s (T 12 ).
Note also that here I s (T 12 ), with T 12 = E 1 + E 2 − 2m e , is an implicit function of p 1 and p 2 .
In the limit of large mass m φ of the assumed scalar mediator we have Hence, collecting all the prefactors the resulting phase space factor reads The expression in Eq. (A23) is almost the same as the phase space factor of standard 2νββ decay when neglecting the final state lepton momenta in the corresponding nuclear matrix element, For 0νββ decay, we can analogously write where c 0ν = G 4 F cos 4 θ C m 2 e /(16π 5 ). Note that employing the above definitions of the phase space factors and comparing Eq. (A10) with Eq. (A11) we get the ratio relating the total decay widths of 2ν SI ββ and 0νββ decay. Using Eq. (A21), we also obtain the differential decay rate dΓ νSI dp 1 dp 2 dc 21 = c νSI |m 2 If only the total kinetic energy of the electrons (T = E 1 + E 2 − 2m e ) is measured, as in most double beta decay experiments, one must calculate the corresponding differential rate dΓ νSI /dT by integrating Eq. (A28) over dp 1 , dp 2 and dc 21 while keeping T at a given value. Here the integral over dc 21 can be done analytically, while the dp 1 dp 2 part has to be evaluated numerically.
Likewise, to derive the angular distribution dΓ/dc 12 , we integrate Eq. (A28) over p 1 and p 2 to obtain the general form where Γ is the total decay rate and k θ the angular correlation for the mode in question (2ν SI ββ, 0νββ, 2νββ).

Nuclear Part
As illustrated above, under the very good approximation that the momenta of the final state leptons can be neglected compared to the momentum flow of the internal neutrino propagators, the resulting NMEs of 2ν SI ββ and 0νββ decay will be identical. We here briefly summarize the method of calculation of the latter. For our numerical analysis we use the NME calculations in the IBM-2 [42], Shell Model [49] and QRPA [50] nuclear structure frameworks.
The key quantities entering the microscopic description of double beta decays are the nuclear matrix elements, values of which have to be obtained using demanding nuclear structure calculations. Let us identify Table II. Definitions of the double beta decay Fermi (MF ), Gamow-Teller (MGT ) and tensor (MT ) NMEs entering Eq. (A30) and the corresponding reduced form-factor producth(q 2 ). We use here the usual notation with τ + being the isospin-raising operator and |0 + I and |0 + F denoting the initial and final states of the nucleus, respectively. The q-dependent functions h•(q 2 ) = v(q 2 )h•(q 2 ) are enhanced by the neutrino potential of the standard light neutrino exchange in Eq. (A31). The subscript X collectively denotes the three possibilities X = V, W, T1 sharing the same form factor shape parameter mV . The Pauli matrices σ a,b stand for the spins of the individual nucleons a, b and the tensor NMEs include the definition S ab = 3(σa · q)(σ b · q) − (σa · σ b ).
now the elementary nuclear matrix elements necessary for computing the exotic neutrinoless double beta decay mechanism introduced in this text.
Under the assumption of negligible momenta of the outgoing neutrinos the nuclear matrix element |M νSI | entering Eq. (A11) can be taken to be approximately equal to the nuclear matrix element of the standard mass mechanism, M νSI ≈ M 0ν , thus allowing for writing the νSIββ decay rate as in Eq. (2). This approximation is reasonable, as the propagating momentum p ∼ 10 2 MeV, while the momenta of leptonic final states are of ∼ 1 MeV.
Following the standard literature [57,58] the NME M 0ν for the 0 + → 0 + transition can be written as

48
M P P GT + M P P T . (A30) In the above, g X = F X (0) denotes the form factor charge, i.e. the value of the nuclear form factor F X (q 2 ) at zero momentum transfer. The standard values used for the charges read: g V = 1, g A = 1.269, g W = 3.7 and g P = 231. Especially crucial is the axial vector coupling g A , quenching of which is the subject of ongoing discussions in the 0νββ community [59]. For the purpose of this work we consider the usually employed unquenched value for free nucleon, g A = 1.269. The q-dependence arising from the product of the reduced form factors F X (q 2 )/g X is included in the nuclear matrix elements appearing in Eq. (A30). 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 Tab. II.
In addition to the product of the reduced nucleon form factors, the NMEs listed in Tab. II also contain the socalled neutrino potential describing the q-dependence of the underlying particle physics mediator of 0νββ decay. In the standard formulation of [57] and [58] the two-body transition operator is constructed in momentum space as the product of neutrino potential, v(q), and the product of the reduced form factorsh(q 2 ). For the standard mass mechanism the neutrino potential reads v(q) = 2 π 1 q(q +Ã) , where neutrino mass has been neglected in comparison with the typical internal neutrino momentum q ∼ 100 MeV, andÃ is the closure energy, which can be adopted from Ref. [60] or estimated usingÃ = 1.12 A 1/2 MeV. The above expression therefore describes the long-range exchange of an essentially massless neutrino mediating the 0νββ decay.