Computation of products of phase space factors and nuclear matrix elements for the Double Beta Decay

The nuclear matrix elements (NME) and phase space factors (PSF) entering the half-life formulas of the double-beta decay (DBD) process are two key quantities whose accurate computation still represents a challenge. In this paper we propose a new approach of calculating them, namely to compute directly their product as an unique formula. This procedure allows a more coherent treatment of the nuclear approximations and input parameters appearing in both quantities and avoids possible confusion in interpreting the DBD data due to different individual expressions adopted for PSF and NME (and consequently their reporting in different units) by different authors. Our calculations are performed for both two neutrino ($2\nu\beta\beta$) and neutrinoless ($0\nu\beta\beta$) decay modes, and for five nuclei of most experimental interest. Further, using the most recent experimental limits for $0\nu\beta\beta$ decay half-lives, we provide new constraints on the light mass neutrino parameter. Finally, by separating in the half-lives formulas the factor representing the axial-vector constant to the forth, we advance suggestions on how to reduce the errors introduced in calculation by the uncertain value of this constant by exploiting the DBD data from different isotopes and/or decay modes.


I. INTRODUCTION
The double-beta decay is a rare nuclear process intensively studied due to its potential to test nuclear structure methods and investigate beyond Standard Model (SM) physics [1]- [3]. According to the number and type of the released leptons there are several possible DBD modes, that can be classified in two categories. One category is that where two anti-neutrinos or two neutrinos are emitted in the final states besides the two electrons (2νβ − β − ) or two positrons (2νβ + β + ). The double-positron decays can also be accompanied by one or two electron capture processes (2νβ + EC, 2νECEC). These decay modes occur with lepton number conservation (LNC) and are allowed within the SM. In the other category enter decay processes similar with the above ones, but where no anti-neutrinos or neutrinos are emitted in the final states. They are generically called neutrinoless DBD processes (0νββ), so we may have 0νβ − β − , 0νβ + β + , 0νβ + EC and 0νECEC decays in this category. All these processes violate LNC, hence they are not allowed within the original framework of the SM but can appear in theories more general than the SM. The discovery of any 0νββ decay mode would firstly demonstrate the lepton number violation by two units, but would also provide us with valuable information on other beyond SM processes. From the 2νββ decay study one can get information about nuclear structure, test different nuclear methods and investigate the violation of Lorentz symmetry in the neutrino sector, while from the 0νββ decay study one can * sabin.stoica@cifra.infim.ro decide about the neutrino character (is it a Dirac or a Majorana particle?), one can constrain beyond SM parameters associated with different mechanisms that may contribute to this decay mode and one can get information about neutrino mass hierarchy, existence of heavy neutrinos, of right-handed components in the weak interaction currents, etc. That is why, the DBD study is a very important and timely topic.
The first step in theoretical study of the DBD process is to derive half-lives expressions and calculate the quantities therein, for each possible decay mode and for different transitions and mechanisms that may contribute to the 0νββ decay mode. With good approximation, the DBD half-lives formulas can be written in factorized forms, as follows [4], [5]: where G (2,0)ν are the PSF, M (2,0)ν are the NME, for the (2, 0)ν decay modes and η l is a parameter related to the specific mechanism l that can contribute to the 0νββ decay. We note that the half-lives expressions from above are written such that the product of the nuclear (NME) and atomic part (PSF) is expressed in [yr −1 ]. Also, we note that the axial-vector constant to the forth power is separated from the other components. Such a form of the half-lives expressions allows an easy using of the theoretical results for interpreting the DBD data and possibility to make connections between data from different decay modes and experiments in an attempt to find solutions to reduce the errors in computation related to the value of axial-vector constant which is not precisely known. As seen, for estimating/predicting DBD lifetimes and deriving beyond SM parameters, a precise, reliable computation of both the PSF and NME is mandatory. The largest uncertainties in the DBD calculations come from the NME. They are calculated with different nuclear methods, the most currently employed being pn-QRPA [3], [6]- [10], Shell Model [11]- [14], IBA2 [15]- [16], PHFB [18], GCM with EDF [19]. They differ each other mainly by the choice of the model spaces and type of correlations taken into account in calculation. Each of these methods has its own advantages and drawbacks, and errors in the NME computation associated with each of them have been extensively debated in the literature over time [3], [6]- [19]. The differences in the NME values computed with different methods may come from different sources such us i) the choice of the model space of single-particle orbitals and type of the nucleon-nucleon correlations included in calculation which are specific to different nuclear methods, ii) the nuclear structure approximations associated with the short range correlations (SRC), finite nucleon size (FNS), higher order terms in the nucleon currents (HOC), inclusion of deformation, etc., or iii) the use of input parameters whose values are not precisely known, like nuclear radius, the average energy of the virtual states in the intermediate odd-odd nuclei or the value of the axial-vector constant, g A , etc. Particularly important is the value of g A (which can be 1.0 = quark value; 1.273 = free nucleon value; or other quenched value (0.4-0.9) because the dependence of the half-lives on this constant is strong. We note that errors coming from the different choice of values of these parameters can increase significantly the uncertainty in the half-lives computation, hence appropriate attention should also be paid to this source or errors. On the other hand, the PSF have been considered during long time to be computed with enough accuracy [3], [20]- [25]. However, newer calculations [4]- [5], [26] performed with more rigorous methods, i.e. by using exact electron Dirac wave functions (w.f.) and improving the way of taking into account the finite nuclear size (FSN), electron screening effects and more realistic form of the Coulomb potential, revealed notable differences of the PSF values as compared with older results, especially for heavier nuclei, for positron emitting and EC decay modes and for transitions to excited states.
The errors in the PSF computation can come from i) the method of calculation of the electron w.f., namelynon-relativistic approach [20]; -relativistic approach with approximate electron w.f. [3]; -relativistic approach with exact electron w.f. [4]- [5], [26]; ii) numerical accuracy both in the resolution of the Dirac equations for getting the electron radial functions and in the integration of the PSF expressions, for different decay modes.
We also note that some input parameters appear both in the NME and PSF expressions, such as the axial-vector constant g A , the nuclear radius R A (R A = r 0 A 1/3 ), the value of the average energy of the virtual states in the intermediate odd-odd nucleus, used in the closure approximation, E N , etc. Also, when these quantities are calculated separately, different groups have used sometimes different values for these parameters. Moreover the NME and PSF have been reported in different units depending on which factors were included in their expressions, and this led sometimes to some confusion/difficulty in the theoretical predictions and interpretation of the experimental data.
In this paper we propose a new approach of calculating the NME and PSF entering the DBD half-lives, namely to calculate directly their product, in an unique formula, instead of calculating them separately. This is actually natural, since for predicting half-lives and getting information about beyond SM physics from the DBD study, we need to know precisely the product N M E × P SF . The computation of the product as a whole has some advantages. Calculating its values in units of [yr −1 ] one facilitates its using in predicting and interpreting the DBD experimental data, by removing any confusion related to the units in which its components are reported when they are calculated separately. Also, the formula of the product has an unique dependence of a certain parameter, for which one takes a single value. Thus, the computation of the atomic and nuclear part of the DBD half-lives gets coherence, of which has not been paid attention to so far. Finally, we note that the separation of the g 4 A factor in the half-life expressions can also have advantages. For example, by combining experimental data and information from different DBD isotopes and/or decay modes and transitions, one can reduce the uncertainty of the calculation related to this parameter.

II. PRODUCTS OF PHASE SPACE FACTORS AND NUCLEAR MATRIX ELEMENTS
We define the products as follows: so, the half-lives expressions become: where g A,ef f is the effective value of the g A constant that can be different for different nuclei and decay modes, because it can depend on nuclear medium and many-body effects. Hence, providing the products P (2,0)ν in [yr −1 ] one can use them easily for predicting half-lives and/or constraining beyond SM parameters. The detailed ex-pressions of these products read: i where G is the Fermi constant, θ C the Cabbibo angle, Q ββ the Q-value for the DBD, m e the electron mass, and ǫ 1,2 and ω 1,2 are the electron and neutrino energies, respectively. Also: where E N is an average energy of the states E I in the odd-odd intermediate nucleus that contribute to the decay. K N and L N are quantities that depend on the electron and neutrino energies, as well as on the energies E N and E I [21]. f 11 are combinations of the radial electron functions g k and f k , solutions of the Dirac equations [26]. Finally, M (2,0)ν are the NME for 2ν and 0ν decay modes.
For computing the products P 2ν and P 0ν we build up numerical codes taking advantage of our previous codes for computing separately the NME and PSF quantities [13], [26], [36]. The expressions of the products P (2,0)ν contain factors outside the integrals stemming from the multiplication and simplification of factors that multiply separately the nuclear and kinetic parts. Also, their kinetic part (phase space factors) and the nuclear part (NME) have common input parameters as R A , E N and g A .
Firstly, we refer to the P 2ν computation. The kinetic part is computed following the main lines of the approach developed in our previous works from refs. [5]- [26] and here we shortly review the main ingredients of the code and computation. We first use a subroutine where the electron wave functions are got as radial solutions (g k and f k ) with appropriate asymptotic behavior of the Dirac equations with a Coulomb-type potential, and including the finite nuclear size and electron screening effects. The Coulomb-type potential is obtained from a realistic proton density in the daughter nucleus. For getting the single particle densities inside the daughter nucleus, we solve the Schrodinger equation for a spherical Woods-Saxon potential with spin orbit and Coulomb terms [5]-  [26]. Then, the PSF part of the code is completed by performing the integrals over the electron phase factors build up with the Dirac radial functions. The code has an improved numerical accuracy for finding the electron w. f. and a better interpolation procedure for integrating the PSF final expressions, especially at low electron energies. For the NME part we use a code similar to that from ref. [13] for computing the double Gamow-Teller transitions, using the following effective nucleon-nucleon interactions: GXPF1A [27] for 48 Ca, JUN45 [28] for 76 Ge and 82 Se and gcn50:82 [29] for 130 T e and 136 Xe. The values for the products P 2ν are presented in the third column of the Table 1 for five nuclei of experimental interest. With the values of g 2ν A,ef f written in the forth column of the table as first entries, we reproduced the most recent measured half-lives found in literature, which are displayed in the second column. In the forth column we also show the g 2ν A,exp values taken from Refs.
[33]- [34], which were obtained by comparing the theoretical B(GT) strengths with the experimental ones extracted from charge-exchange reactions. In the last column we present the difference in percentage between the g A,ef f values obtained within our calculations to reproduce the experimental DBD half-lives and those obtained by fitting the B(GT) experimental data, estimated in percentage, ǫ = g 2ν − g 2ν A,ef f /g 2ν . As seen, the two sets of values are close to each other, the smallest differences being in the case of 76 Ge and 48 Ca nuclei.
Then, we calculated the P 0ν products in the case of the light Majorana neutrino exchange mechanism, with η l = m ν /m e and the light neutrino parameter defined as: where U ek are the first row elements of the Pontecorvo-Maki-Nakagawa-Sakata neutrino mixing matrix and m k are the light neutrino masses [35]. The expression of the nuclear matrix elements can be written as a sum of Gamow-Teller (GT ), Fermi (F ) and tensor (T ) components [9], [36]: where M 0ν GT , M 0ν F and M 0ν T are these components. These are defined as follows: O α mn are transition operators (α = GT, F, T ) and the summation is over all the nucleon states. Correspondingly, the two-body transition operators O α 12 can be expressed in a factorized form as [36]: where N α is a numerical factor including the coupling constants, and S α , R α and C α are operators acting on the spin, relative and center-of-mass wave functions of the two-particle states, respectively. Thus, the calculation of the matrix elements of these operators can be decomposed into products of reduced matrix elements within the two subspaces [14]. The expressions of the two-body transition operators are: The O α 12 operators contain three components, namely, the spin, center-of-mass and radial ones, and the expectation values of the first two components can be easily managed. The radial part is the most difficult to be calculated because it contains neutrino potentials written in different approximations, and the expectation values are double integrals over them. Also, short-range correlations and finite nucleon size corrections are introduced in this part of computation. The neutrino potentials depend weakly on the intermediate states, and are defined by integrals over momentum carried by the virtual neutrino exchanged between the two nucleons [9]. They include Fermi (F), Gamow-Teller (GT) and tensor (T) components: where R = r 0 A 1/3 fm, with r 0 = 1.2f m, j 0,2 (qr) are the spherical Bessel functions and the integrals are over the neutrino exchange momentum q. In our calculations we use the closure approximation and E N , as mentioned above, represents the average energy of the virtual states in the intermediate odd-odd nucleus included in the description of the decay. Also, we note that the factor 2R is canceled by the similar one from the denominator of the PSF expression, so the P 0ν does not depend on the nuclear radius. The expressions of neutrino potentials h F,GT,T can be found in many references (see for example [9]). These expressions include FNS effects taken into account through vector and axial-vector form factors, G V and G A [9].
(16) We take the following values for the vector and axial vectors form factors: Λ V = 850M eV and Λ A = 1086M eV [2].
For computing the radial matrix elements nl|H α |n ′ l ′ we use the harmonic oscillator HO wave functions ψ nl (lr) and ψ n ′ l ′ (r) corrected by a factor [1 + f (r)], which takes into account the SRC effects induced by the nuclear interaction [36]: For the correlation function we take the functional form For the a, b and c constants we use the parametrization used in ref. [37].
Including HOC and FNS effects the radial matrix elements of the neutrino potentials become: We note that in the case of P 0ν products, the axial-vector constant enters also the expressions of the neutrino potentials, in addition to the factor g 4 A,ef f , so the half-life expression for the 0νββ decay, Eq. (5), contains a "double" dependence on this constant. Of course, for coherence, the same values of the g A,ef f constant should be taken in both places, i.e. both in the P 0 ν and in the half-life computation. We note that these values may differ from the values of this constant used in the 2νββ decay mode. Because we do not know until now what is the correct value of g A.ef f for 0νββ decay mode, we calculated the P 0ν products for the free nucleon value  Table 2.
At this point it is worth to mention that the values of the products P (2,0)nu from Tables 1 and 2, obtained with the approach described here, are very close to the values that we obtained when we computed these products by multiplying the values of NME and PSF calculated separately, as one should be. This is understandable because in their calculation by the two methods we used the same values of the input parameters and the same nuclear approximations and parametrizations. The small differences come from the numerical precision of the numerical codes, we used. We emphasize, however, that the importance of our current approach is that it can eliminate the incoherence of using NME and PSF values calculated separately with different values for common nuclear parameters, which can introduce significant errors in the evaluation of N M E × P SF product as a whole. The errors in the evaluation of these products can indeed be significant if one takes different values of g A in the computation of N M E and P SF and if these values are not the same with the value used in the g 4 A factor. For example, the errors introduced in the NME computation by the use of a quenched (1.0) or an unquenched (1.27) value of g A were analyzed in Ref. [36] for 48 Ca, 76 Ge and 82 Se nuclei, and found to be within 10-14% (without the factor g 4 A ). The use of different values for the other (common) parameters involved in calculations as the nuclear radius, < E N >, etc. can bring additional uncertainties of the same order. The errors can be amplified by the use of different values of these parameters in the PSF computation. So, there may be relevant errors in calculating the products N M E × P SF when the N M E and P SF values are taken from separate calculations reported in literature.
Then, we revise the limits of the light neutrino mass parameter m ν using our calculations and the most recent experimental limits reported for the 0νββ decay half-lives. These results are presented in the last column of the Table 2. One observes that presently the most stringent constraints on this parameter come from the nuclei 76 Ge and 136 Xe, due to both the experimental results (the lowest limits measured at present for the 0νββ decay half-lives [38], [41]) and to accurate theoretical calculations. An important issue in this case remains the use of a correct value for the g A constant. As far as this value is still unknown, for accurate half-lives predictions and constrains of beyond SM parameters, the goal is to reduce the errors associated with this constant. One suggestion is to use information from different decay modes and/or from DBD experiments on different nuclei. For example, for a particular nucleus the ratio of the 2ν and 0ν half-lives expressions reads: As seen from the above formula, any information that we can get about the relative magnitude of the g A,ef f values for the 2ν and 0ν decays in the same nucleus, can be exploited to improve the constraints on the neutrino mass parameter, when improved calculations of P (2,0) are available. Also, referring to two different nuclei, denoted with m and n, the ratio of their half-lives reads: For 2νββ decay mode: For 0νββ decay mode we have: As seen, one can deduce g 2ν A,ef f for one particular nucleus if one knows with (more) precision the value of this constant for another nucleus, using the experimental halflives and the calculated P 2ν for both nuclei. For example, one can take advantage of the possible experimental determination of this parameter for some particular isotopes, as it was recently proposed in Ref. [42]. Similar considerations, i.e. the exploitation of data from several experiments are valid for predicting 2νββ decay half-lives for a nucleus that was not yet measured, if we have accurate data for another nucleus and good estimations of the g 2ν A,ef f value from other (non-DBD) experimental data. For such predictions information on DBD half-lives not yet measured obtained from empirical formulas, as is proposed in Ref. [43], are valuable. Similarly, for the 0νββ decay one can deduce more information about the effective value of g 0ν A for a particular nucleus if we know this value for another. For example, we might know g 0ν A,ef f with more precision in the case of nuclei where single state dominance (SSD) approximation is valid, where the half-life can be computed with reasonable precision by taking into account only one state in the intermediate odd-odd nucleus (for example 100 Mo case), and where g 2ν A,ef f and g 0ν A,ef f might have close values.

III. CONCLUSIONS
We proposed a new approach of calculating the NME and PSF for DBD, by computing directly their product. The product as a whole can be computed more consistently, having an unique dependence of some parameters entering previously separately the NME and PSF expressions and taking thus single values for them. The values of the product are given in the same units as T −1 1/2 (i.e. [yr −1 ]) removing any possible confusion in using the theoretical calculations for interpreting the DBD data. The new codes of calculating the N M E × P SF products include improved routines used in our previous papers for the separately computation of these two quantities. We provide values of these products for 2ν and 0ν DBD modes for five nuclei of most experimental interest. Then, using our calculations and the newest half-life val-ues for the 0νββ decays reported in literature, we revise the upper limits for the light mass neutrino parameter. In the half-lives formulas we separate the strong dependence on the axial-vector constant, i.e. the factor g 4 A , that bring a large uncertainty in calculation and suggest some ways to reduce/avoid the errors related to the uncertain value of this constant. This could be done by using ratios of g (2,0)ν A,ef f and P (2,0)ν (instead of their individual values) and exploiting data on the same nucleus but for different decay modes and/or DBD data from experiments on different nuclei, including the possibility that this constant to be determined experimentally for some particular isotopes. We hope our work will be a forward step for more consistent DBD calculations and which will be easier used in predicting and interpreting the experimental data.