The interpolating formula for the $0\nu\beta\beta$-decay half-life in the case of light and heavy neutrino mass mechanisms

We revisit the"interpolating formula"proposed in our previous publication. It allows one to calculate the neutrinoless double beta decay ($0\nu\beta\beta$-decay) half-life for arbitrary neutrino mass without involvement of the complicated results for nuclear matrix elements (NME) obtained within specific nuclear structure models. The formula derives from the finding that the value of a properly normalized ratio of the NMEs for the light and heavy neutrino mass mechanisms weakly depends on isotope. From this fact it follows, in particular, that the light and heavy neutrino mass mechanisms can hardly be distinguished in a model independent way searching for $0\nu\beta\beta$-decay of different nuclei. Here we show that this formula holds for all the known nuclear structure approaches. We give a mathematical justification of our results examining analytical properties of the NMEs. We also consider several simplified benchmark scenarios within left-right symmetric models and analyze the conditions for the dominance of the light or heavy neutrino mass mechanisms in $0\nu\beta\beta$-decay.


I. INTRODUCTION
Neutrinoless double beta (0νββ) decay is a Lepton Number Violating (LNV) process changing lepton number in two units ∆L = 2. It is forbidden in the Standard Model (SM), where lepton number is conserving. Basically there are two sources of LNV: Majorana neutrino mass and LNV vertices. The latter may emerge from numerous high-scale models giving rise to the corresponding mechanisms of 0νββ-decay. Once this process is observed, the question of distinguishing between the dominant mechanisms will arise. Certainly, this task is highly non-trivial. One may hope that measurements of 0νββ half-life with different isotopes would facilitate its solution due to variability of Nuclear Matrix Elements (NME) of particular mechanisms from one isotope to the other. In the present paper we show that at least the light and heavy Majorana neutrino mass mechanisms are indistinguishable in this way without additional hypothesis. This fact becomes especially comprehensible in terms of the so called "interpolating formula" (IntF) [8] merging the light and heavy neutrino mass ranges in the NMEs and allowing a transparent physical interpretation of the above fact. The IntF is a simple analytical formula representing with an accuracy of 30% or better the NME as a function of the Majorana neutrino mass. This accuracy is sufficient for practical purposes taking into account the limited accuracy of the available nuclear structure approaches to the NME calculations. In what follows we will show that that the IntF is valid for all these nuclear structure approaches with the above-indicated accuracy and elucidate some of its other useful properties. On the particle physics side we adopt a generic scenario with Majorana neutrinos of arbitrary value masses and consider their contribution to 0νββ-decay via mass mechanism mediated by both leftand right-handed weak currents. Then for the sake of concreteness we consider the neutrino mass mechanism within the left-right symmetric models (LRSM) [3,4] and extent our analysis towards some more particular scenarios.

II. NEUTRINO MASS MECHANISM OF 0νββ-DECAY
We start with a generic Majorana neutrino mass mechanism of 0νββ-decay induced by the low-energy effective Lagrangian with the left/right-handed hadronic J L/R and leptonic j L/R currents. As usual, G β = G F cos θ C , where G F and θ C are Fermi constant and Cabbibo angle, respectively. The dimensionless parameter λ depends on the underlying high-scale model. In the particular case of the Left Right Symmetric (LRS) models, based on the SU (2) L ⊗ SU (2) R ⊗ U (1) B−L gauge group [3,4], the Lagrangian (1) appears at low energies after integrating out W ± L,R massive gauge bosons. In this model where M W L and M W R (M L < M R ) are masses of W L and W R gauge bosons, respectively. The current constrain on the mass of W R is M W R ≥ 2.9 TeV [15] sets the limit The upper limit λ = 7.7 × 10 −4 we use everywhere in the present paper as a reference value for this parameter.
Since we focus on the mass mechanism we discarded in Eq. (1) the j L,R J R,L terms irrelevant in this case (for a review see, for instance [1]). In Eq. (1) the explicit form of the left-and right-handed hadronic currents J † L,R in nuclei can be found, e.g., in Ref. [7]. The leptonic currents are given by The ν eL and ν eR are the weak eigenstate electron neutrinos, which are expressed as superpositions of the light and heavy Majorana mass eigenstate neutrinos ν j and N k as where the unitary matrix is the generalizations of the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrix, which diagonalizes the general (3 + n) × (3 + n) neutrino mass matrix in the basis (ν eL , ν µL , ν τ L , N C 1R , ..., N C nR ). Here M L,R and M D are Majorana and Dirac mass terms, respectively. After diagonalization one should end up with 3 light ν i (i=1,2, 3) and n heavy N k (k=1,...,n) Majorana neutrino mass eigenstates with the masses m i and M k , respectively. In the LRS models n = 3. The smallness of m i can be guaranteed by the seesaw-I condition M R M D . As is well known this leads to very heavy states N k with masses M k 1 TeV being beyond the experimental reach. In the scenarios with n > 3 the inverse seesaw mechanism can be implemented. In this case among N k , accompanying the light ν i states, there can appear moderately heavy or even light Majorana states. Actually, their masses can be of arbitrary value. This is the case of our particular interest, for which we designed the above-mentioned interpolating formula.
Assuming the dominance of the mass mechanism we write down the 0νββ-decay half-life The proton mass is denoted by m p and g A is the unquenched value of axial-vector coupling constant (g A = 1.269). The phase-space factor G 0ν is tabulated for various 0νββ-decaying nuclei in Ref. [9]. The NMEs M 0ν as functions of neutrino mass m ν (m ν = m i or M k ) are given by [8] M 0ν Here, R and m e are the nuclear radius and the mass of electron, respectively. We use as usual R = r 0 A 1/3 with r 0 = 1.2 fm. Initial and final nuclear ground states with energies E I and E F are denoted by |0 + I and |0 + F , respectively. The summation runs over intermediate nuclear states |n with energies E n . The weak one-body nuclear charged current J L,R [7,8] depends on the effective value of axial-vector coupling constant g eff A of the nucleon, which is renormalized to a smaller, the so called, "quenched" value, g eff A [10].   (13), (14) for a given isotope and their average value p 2 a used in Eq. (17) with the variance σ (in parentheses) calculated within different nuclear structure approaches: interacting shell model (ISM) (Strasbourg-Madrid (StMa) [16] and Central Michigan University (CMU) [17] groups), interacting boson model (IBM) [18], quasiparticle random phase approximation (QRPA) (Tuebingen-Bratislava-Caltech (TBC) [19,20] and Jyväskyla (Jy) [21] groups), projected Hartree-Fock Bogoliubov approach (PHFB) [22], and covariant density functional theory (CDFT) [23]. The Argonne, CD-Bonn and UCOM two-nucleon short-range correlations are taken into account. The non-quenched value of weak axial-vector coupling gA is assumed. with Here the NMEs M 0ν ν , M 0ν N are derived from the NME M 0ν in Eq. (8) in the following way In the case of a neutrino spectrum with mass states N k of an arbitrary mass value one has to apply Eq. (8) for the NME calculations, resulting in a complicated function of the neutrino mass. This is a real hassle for use in practice. Fortunately, there is a very good approximate analytical representation for Eq. (8) proposed in Ref. [8] (and references therein) and having a remarkably simple form This is what we call "interpolating formula" (IntF) since it interpolates two limiting cases (11), (12) and is valid to a good accuracy for an arbitrary value of m ν . Eq. (13) contains the parameter with the dimension of (mass) 2 . The form of Eq. (13) suggests the interpretation of p 2 as the mean square momentum of the virtual neutrino propagating between two β-decaying nucleons. Therefore, it is expected to be of the order of p 2 F ∼ (200 MeV) 2 . The current values of the matrix elements M 0ν ν and M 0ν N calculated within different nuclear structure approaches can be found in Tables 6 and 7 of Ref. [10]. The value of corresponding parameter p 2 is given for various isotopes together with its averaged value p 2 a with variance σ in Table  I. The unquenched value of axial-vector coupling constant is assumed: g eff A = g A = 1.25 − 1. 27. We see that the value of p 2 depends noticeable on the chosen nuclear structure method and considered choice of twonucleon short-range correlation function. The values of p 2 a are displayed for different nuclear structure approaches and types of two-nucleon short-range correlations in Fig. 1. The largest value of the parameter p 2 a 200 MeV is found for the QRPA with isospin restoration and CD-Bonn two-nucleon short range correlations. Surprisingly, within all the considered nuclear structure approaches the variance σ is very small being of the order of 3-10 %, i.e. the value of p 2 is practically the same for all isotopes of experimental interest and can be replaced with averaged value p 2 a . In Appendix we discuss this finding from the view point of the analytical properties of the NME in Eq. (8) as a function in the complex plane of m ν . The above conclusion is also supported by the statistical treatment of M 0ν ν and M 0ν N NMEs performed in Ref. [11].
Using the parameter p 2 a in the "interpolating formula" (13) we can write to a good accuracy the 0νββdecay half-life for the Majorana neutrino exchange mechanism as where and for arbitrary mass M k . The sum runs over j = 1, 2, 3 and k = 1, ..., n. The values of parameter C νN are given for various isotopes in Table II. The "interpolating formula" in Eq. (15) reproduces the "exact" QRPA result with rather good accuracy except for the transition region where its deviation, as seen from Fig. 8, amounts 20% -25%. The parameter η νN is a general LNV parameter for the light and heavy neutrino mass mechanisms, which is practically independent of the isotope under consideration.

IV. LIGHT VS HEAVY NEUTRINO MASS MECHANISMS
From the conclusion of the previous section and Eq. (17) it follows that contrary to the previous expectations in the literature (see for instance Refs. [12,13]) the dominance of light or heavy neutrino mechanisms of 0νββ-decay cannot be recognized just by observation of this process with different isotopes. An additional theoretical or experimental input about neutrino masses and mixing is needed to shed light on the particular role of each of these mechanisms.
Let us give a couple of examples of model inputs allowing us to distinguish two above-mentioned mechanisms.
For a scenario with three SM singlet neutrinos ν e,µ,τ R the 6 × 6 mixing matrix U in Eq. (5) is completely parameterized with 15 angles, 10 Dirac and 5 Majorana CP violating phases. Let consider some viable structures of this mixing matrix.

A. Uncoupled light and heavy neutrino sectors
In the particular case of  there is no mixing between heavy and light neutrino sectors. Then we have  [16] and Central Michigan University (CMU) [17] groups), interacting boson model (IBM) [18], quasiparticle random phase approximation (QRPA) (Tuebingen-Bratislava-Caltech (TBC) [19,20] and Jyväskyla (Jy) [21] groups), projected Hartree-Fock Bogoliubov approach (PHFB) [22] and covariant density functional theory (CDFT) [23] are considered. The Argonne, CD-Bonn and UCOM twonucleon short-range correlations are taken into account. The non-quenched value of the weak axial-vector coupling gA is assumed.  In this scenario U 0 can be identified with the PMNS mixing matrix U . Thus we assume U 0 = U . The mixing matrix V 0 for the heavy neutrinos is unknown, but it is similar to U 0 in the light neutrino sector, then V 0 = U is frequently assumed. For sake of simplicity we consider two different cases for the heavy neutrino masses: In the case of the constant products ζ p = m i M i we have for the LNV parameter in Eq. (19): Thus, in this scenario the presence of heavy neutrinos leads to a vertical shift of the standard plot in Fig. 2 by a constant factor κ. As a result, the 0νββ-decay experimental upper bound on m ββ is significantly less stringent, if ζ p λ p 2 a 24 MeV 2 . In our estimation we used the upper-bound-value in (3), i. e. λ = 7.7 · 10 −4 , and p 2 a = 175 MeV calculated within the QRPA by assuming Argonne potential and g A = 1.27 (see Table II).
In the case of the constant ratios ζ r = m i /M i in Eq. (22) the effective Majorana neutrino mass M ββ is shown in Fig. 3. Contribution of M ββ becomes comparable to m ββ as soon as ζ r = 10 −17 , which corresponds to M i ∼ 10 16 eV = 10 4 TeV. λ = 7.7 · 10 −4 is assumed again. Notice the reversed behavior of the mass hierarchies: NH no longer exhibits a region unbounded from below, while IH does.

B. Seesaw-mixed light and heavy neutrino sectors
Assuming for simplicity the flavor universal mixing between the active and sterile neutrino sectors the seesaw mixing matrix U takes the form Here, ζ = m D m LN V , where m D is the typical scale of the charged leptons masses and m LN V is the LNV scale of the order of the Majorana masses M i of the heavy neutrinos. As in the previous scenario U 0 can be identified with the PMNS U matrix. Thus we assume U = U 0 . For V 0 , analogue of U 0 in the heavy neutrino sector, we find from the unitarity conditions and It is assumed that a small violation of the unitarity of U 0 and V 0 matrices is beyond the current accuracy of phenomenological determination of elements of the PMNS matrix. The matrix V 0 takes the form We note that each element of the first row is multiplied by the same phase factor e −iα1 . Analogously, the second raw is multiplied by e −iα2 . Therefore, the Majorana phases α 1,2 do not affect the heavy neutrino LNV parameter M R ββ in this case. On the contrary, the Dirac phase δ, which does not affect the light neutrino LNV parameter m ββ will impact the value of M R ββ . The seesaw structure of (24) implies m i m 2 D /m LN V and M i m LN V . For a product of light and heavy neutrino masses let assume m i M i m 2 D . If the LNV scale is significantly larger than p 2 a we find with We note that for m D 5 MeV the coefficient λ p 2 a /m 2 D enetring M R ββ in Eq. (29) is close to unity and it might be that contributions from the light and heavy neutrinos to η νN are comparable. However, M R ββ is not proportional to m ββ as off-diagonal elements of matrices U 0 and (U † ) are different. Therefore, a detailed analysis is needed to establish an useful constraint on the Yukawa potential associated with neutrinos. In  Within the seesaw structure one can also assume For ζ 2 = 10 −17 and λ = 7.7 · 10 −4 the effective mass M R ββ in Eq. (30) is plotted in Fig. 5. We see again that for a chosen set of parameters the value of M R ββ can be comparable with m ββ (see Fig. 2).
In Table III we show upper bounds on η νN in Eq. (28) derived from the current 0νββ-decay experiments. As seen, the most stringent bound comes out from 136 Xe 0νββ-decay experiment Ref. [31]. For this bound we analyzed the separate contributions of the light and heavy neutrinos to 0νββ-decay. Fig. 6 displays the corresponding results in the plane of the parameters m D and m 0 (m i m 2 D /M i is assumed) for the cases of the normal (left panel) and inverted hierarchy (right panel) of neutrino masses. We see that in the considered scenario for normal (inverted) hierarchy the values m D ≤ 1.5 MeV (m D ≤ 2.9 MeV) are already excluded by the existing experimental data on 0νββ-decay. We also see that in the case of normal (inverted) neutrino mass hierarchy the heavy neutrino exchange mechanism cannot dominate over the light one in the region m 0 ≥ 0.08 eV (m 0 ≥ 0.065 eV). The contraint from the 0νββ-decay experiment imply that the limit on the mass of lightest heavy neutrino is M 3 > 38 TeV and M 2 > 171 TeV in the cases of normal and inverted hierarchy, respectively. Fig. 7 shows results of an analysis similar to the above-discussed one, but for ζ 2 = m i /M i scenario. In this case the heavy neutrino mechanism cannot dominate in practically the same domain of m 0 as previously. It is concluded that in the case of normal (inverted) hierarchy ζ ≤ 1.75×10 −8 (ζ ≤ 1.65×10 −8 ). We note that within the considered seesaw scenario within the LRSM the effective Majorana neutrino mass m ββ can not be identified with the first element of the Dirac-Majorana mass (see Appendix B) (M L ) ee , which contains additional term ζ 2 M 1 in magnitude comparable with m 1 . The corresponding term in m ββ has been neglected as it is suppressed by properties of neutrino propagator for large neutrino mass. Due to the same reason M ββ can not be identified with (M R ) ee .

V. CONCLUSIONS
In summary, we have shown that the ratio of nuclear matrix elements for the light and heavy neutrino mass mechanisms exhibits practically no dependence on isotope for all favoured nuclear structure methods. This quantity, when properly scaled, can be identified with squared average neutrino momentum p 2 of the interpolating formula including light and heavy neutrino exchange mechanisms. The universality of the averaged value of p 2 for a set of isotopes allows determination of a new LNV parameter η νN , which is a coherent sum of squared LNV parameters m ββ and M R ββ characterizing the light and heavy neutrino exchange mechanisms, respectively. Thus, the observation of 0νββ-decay on two and more nuclear isotopes will allow one to deduce information about the size of η νN , but not about the relative contribution of light or heavy neutrino-exchange mechanism to the decay rate. An additional theoretical or experimental input about neutrino masses and mixing is needed to shed light on the particular role of each of these mechanisms. As an example we considered a simplified see-saw type 6 × 6 neutrino mixing matrix (24), which implies that the 3 × 3 mixing matrix of heavy neutrinos is the hermitian conjugate of the 3 × 3 PMNS mixing matrix of light neutrinos. Assuming several viable seesaw relations among the light m i and heavy M i neutrino masses (i=1,2 and 3) useful constraints on the parameters, in particular Dirac neutrino mass m D , entering these relations have been  Here we give some comments on the possible improvement of our interpolating formula in Eq. (13), which we call the "monopole" approximation. Numerically the latter is already a very good approximation to the "exact" NMEs given by Eq. (8) and calculated in the framework of any specific nuclear structure approach. However in certain cases one may need an approximate formula having not only a good numerical precision, but also the analytical properties in the complex plane of m ν the same or maximally close to the "exact" NME defined in expression (8).
Obviously the monopole approximation (13) has two imaginary poles in the complex plane of m ν , while they are absent in the exact expression (8). Below we describe a class of approximations with the analytic properties of the "exact" NME (8).
Let us rewrite Eq. (8) in the form and ϕ(x, y) = 1 m p m e 4πR g 2 The function ϕ(x, y) describes a distribution of currents inside the nucleus. In Eq. (A1) the neutrino mass enters the denominator of the integrand. Analytic properties of functions defined in terms of a contour integral are fixed by the Landau rules [33,34].
The singular points of the first kind are associated with singular behavior of the integrand at the end points of the integration contour. In the case of Eq. (A1), these singularities could occur provided that χ(p) ≡ E p (E p + ∆) = 0 for p = 0 or ∞. This equation can be fulfilled for p = 0 only to give m = 0 and m = ±∆. The points m = ±∆ are located on the different sheets of the Riemann surface of M 0ν LL,RR (m ν ). It is clear that model dependent features of the nuclear structure entering ϕ(x, y) do not affect the end point singularities.
Singular points of the second kind are associated with the pinch singularities of the integrand. To find them, the equations χ(p)/ϕ(p) = 0 and (χ(p)/ϕ(p)) = 0 are to be solved, which localise high-order poles of the integrand in the complex p-plane. These singularities depend on ϕ(x, y) and thereby on the nuclear structure model.
Analytic properties of M 0ν LL,RR (m ν ) as a function of ∆ are particularly simple. Changing the variable in Eq. (A4) to p = m sinh θ, we arrive at the dispersion integral is an analytic function in the complex ξ-plane with the cut (1, +∞) corresponding to the cut (−∆, 0) in m ν . Provided ϕ(p) is an analytic function for |p| < ∞ and the integral (A4) converges, M 0ν LL,RR (m) turns out to be an analytic function in the complex m ν -plane with the cut (−∆, 0). On the second sheet of the Riemann surface one finds a branch point m = +∆.
As we discussed before the monopole parametrization (13) is numerically very accurate. This parametrization corresponds to an approximation of the spectral function with the delta function: φ m (p) ∼ δ(p 2 − p 2 ). Then for the formula with the correct analytical properties, which we are going to construct here, we chose the spectral function in a form close to the φ m (p) to guarantee its numerical accuracy comparable with the monopole parameterization. We may choose φ(p) = exp − ρ 2 p 2 2 sinh(ρ 2 pp 0 ) ρ 2 pp 0 .
with the free parameters, ρ and p 0 ∼ p 2 1/2 , which can be fixed by normalization to the exact values at zero and at infinity. The function φ(p) for p = p 0 is close to the maximum, the value of 1/ρ determines the width of momentum distribution. This spectral function is analytic for |p| < ∞ and it generates model independent end-point singularities only. The corresponding interpolating formula appears to be an analytic function in the complex m ν -plane with the cut (−∆, 0). The cut position is model-independent. The discontinuity depends on φ(p) and is model dependent. A particularly strong effect on the behavior of analytic functions in a fixed domain comes from nearest singularities. Taking into account that ∆ ∼ 10 MeV, an improved description of the neutrino mass dependence can be expected around zero neutrino mass in the circle with a radius of a few tens of MeV. This scale is smaller than the characteristic momentum transfer p 0 ∼ 200 MeV. Reasonable accuracy is also expected for large m ν domain, provided the spectral function (A5) approximates closely the monopole spectral function found to be successful phenomenologically.
The ratio between the interpolating formula of Eq. (13) and the exact calculation for 76 Ge is shown on Fig. 8. The result is compared to the interpolating formula of the spectral function (A5) with ρ = 5 fm and p 0 = 0.84/fm. For low neutrino masses up to about 40 MeV the analytic interpolation formula approximates the exact result with a better accuracy. For higher masses the nuclear structure at about 200 MeV becomes important, which could reflect a contribution of the model dependent pinch singularities which we do not consider.

Appendix B: Dirac-Majorana neutrino mass term within seesaw in LRSM
In this Appendix the Dirac-Majorana neutrino mass term associated with the see-saw mass mechanism within the LRSM and particular case of neutrino mixing given in Eq. (24) is presented. We have