Computation of products of phase space factors and nuclear matrix elements for double beta decay

The double-beta decay (DBD) 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 released leptons, there are several possible DBD modes that can be classified in two categories. One category is where two anti-neutrinos or two neutrinos are emitted in the final states along with the two electrons ($ 2\nu\beta^-\beta^- $) or two positrons ($ 2\nu\beta^+\beta^+ $). The double-positron decays can also be accompanied by one or two electron capture processes ($ 2\nu\beta^+ EC $, $ 2\nu ECEC $). These decay modes occur with lepton number conservation (LNC) and are allowed within the SM. In the other category, the decay processes are similar with the one described above, however no anti-neutrinos or neutrinos are emitted in the final states. They are generically called neutrinoless DBD processes ($ 0\nu\beta\beta $), such that we may have $ 0\nu\beta^-\beta^- $, $ 0\nu\beta^+\beta^+ $, $ 0\nu\beta^+ EC $, and $ 0\nu ECEC $ decays in this category. All these processes violate LNC, hence they are not allowed within the original framework of the SM, however they can appear in theories that are more general than the SM. The discovery of any $ 0\nu\beta\beta $ decay mode would first demonstrate the lepton number violation by two units, but would also provide us with valuable information on other beyond SM processes. From the $ 2\nu\beta\beta $ decay study, information about nuclear structure can be obtained, different nuclear methods can be tested, and the violation of the Lorentz symmetry can be investigated in the neutrino sector. Meanwhile, from the $ 0\nu\beta\beta $ decay study, the neutrino character can be decided (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 obtain information about neutrino mass hierarchy, the existence of heavy neutrinos, right-handed components in the weak interaction currents, etc. Therefore, the DBD study is a highly important and timely topic.

The first step in the theoretical study of the DBD process is to derive half-life expressions and calculate the quantities therein, for each possible decay mode and for different transitions and mechanisms that may contribute to the $ 0\nu\beta\beta $ decay mode. With a good approximation, the DBD half-life formulas can be written in factorized form, as follows [4], [5]:

where $ G^{(2,0)\nu} $ denotes the phase space factor (PSF), $ M^{(2,0)\nu} $ is the NME for the $ (2,0)\nu $ decay modes, and $ \langle\eta_{l}\rangle $ is a parameter related to the specific mechanism l that can contribute to the $ 0\nu\beta\beta $ decay. We note that the half-life expressions above are written such that the product of the nuclear (NME) and atomic part (PSF) is expressed in $ [\rm yr^{-1}] $. Moreover, we note that the axial-vector constant to the forth power is separated from other components. Such form of the half-life expressions allows facilitated use of the theoretical results for interpreting DBD data and the 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 the axial-vector constant, which is not precisely known. As previously shown, to estimate/predict DBD lifetimes and derive beyond SM parameters, a precise, reliable computation of both the PSF and NME is mandatory. The largest uncertainties in the DBD calculations originate from the NME. They are calculated with different nuclear methods, the most currently employed method being pnQRPA [3, 6–10], Shell Model [11–14], IBA2 [15–17], PHFB [18], and GCM with EDF[19]. They differ mainly by the choice of model spaces and the type of correlations taken into account in the calculation. Each of these methods has its own advantages and drawbacks, and errors in the NME computation associated with each of method have been extensively debated in the literature over time [3, 6–19]. The differences in the NME values computed by different methods may stem from different sources such as 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 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 of the strong dependence of the half-life on this constant. We note that errors originating from the different choice of values of these parameters can significantly increase uncertainty in the half-life computation, hence appropriate attention should also be paid to this source or errors.

In contrast, the PSF have long been considered to be computed with sufficient accuracy [3, 20–25]. However, newer calculations [4–5, 26] using more rigorous methods, i.e., employing exact electron Dirac wave functions (w.f.) and improving the method by which the finite nuclear size (FSN), electron screening effects, and a more realistic form of the Coulomb potential is considered, revealed notable differences in the PSF values as compared with older results, especially for heavier nuclei, positron emission and $ {\rm EC} $ decay modes, and transitions to excited states.

Errors in the PSF computation originate from i) the method of calculation of the electron w.f., namely - the non-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 obtaining the electron radial functions and the integration of 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, $ \langle E_N \rangle $, etc. Moreover, when these quantities are calculated separately, different groups have sometimes used different values for these parameters. Furthermore, 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 theoretical predictions and the interpretation of experimental data.

In this study, we propose a new approach of calculating the NME and PSF entering the DBD half-lives to calculate their product directly, in an unique formula, instead of calculating them separately. This is indeed natural, since to predict half-lives and obtain information beyond SM physics from the DBD study, we need to precisely know the product $ {\rm NME} \times {\rm PSF} $. The computation of the product as a whole has some advantages. By calculating its values in units of $ [\rm yr^{-1}] $, the prediction and interpretation of experimental DBD data is facilitated, thereby removing any confusion related to the units, whose components are reported when these are calculated separately. Moreover, the formula of the product has the unique dependence on a certain parameter, which assumes a single value. Thus, the computation of the atomic and nuclear part of the DBD half-life obtains coherence, which has not been paid attention to so far. Finally, we note that the separation of the $ g_{A}^4 $ factor in the half-life expressions also provides 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.

We define the products as follows:

therefore, the half-life expressions become:

where $ g_{A,{\rm eff}} $ is the effective value of the $ g_{A} $ constant, which can be different for different nuclei and decay modes, as it can depend on nuclear medium and many-body effects. Hence, providing the products $ P^{(2,0)\nu} $ in $ [\rm yr^{-1}] $, one can use them easily for predicting half-lives and/or constraining beyond SM parameters. The detailed expressions of these products read:

where $ G $ is the Fermi constant, $ \theta_C $ is the Cabbibo angle, $ Q_{\beta\beta} $ denotes the Q-value for the DBD, $ m_e $ depicts the electron mass, and $ \epsilon_{1,2} $ and $ \omega_{1,2} $ are the electron and neutrino energies, respectively. Moreover,

where $ \langle E_N \rangle $ is an average energy of the states $ E_I $ in the odd-odd intermediate nucleus that contribute to the decay. $ \langle K_N \rangle $ and $ \langle L_N \rangle $ are quantities that depend on the electron and neutrino energies, as well as on the energies $ \langle E_N \rangle $ and $ E_I $ [21]. $ f_{11}^{(0)} $ depicts the combinations of the radial electron functions $ g_k $ and $ f_k $, which are solutions of the Dirac equations [26]. Finally, $ M^{(2,0)\nu} $ is the NME for the $ 2\nu $ and $ 0\nu $ decay modes.

To compute the products $ P^{2\nu} $ and $ P^{0\nu} $, we construct numerical codes by taking advantage of our previous codes to separately compute the NME and PSF quantities [13, 26, 27]. The expression of the product $ P^{(2,0)\nu} $ contains factors outside the integrals that stem from the simplification and separate multiplication of the nuclear and kinetic parts. Further, their kinetic part (PSF) and the nuclear part (NME) have common input parameters $ R_A $, $ \langle E_N \rangle $, and $ g_A $.

First, we refer to the $ P^{2\nu} $ computation. The kinetic part is computed following the main lines of the approach developed in our previous work from Refs. [5-26]. Here, we shortly review the main ingredients of the code and computation. We first use a subroutine where the electron wave functions are obtained as radial solutions ($ g_k $ and $ f_k $) with appropriate asymptotic behavior of the Dirac equations with a Coulomb-type potential, including the finite nuclear size and electron screening effects. The Coulomb-type potential is deduced from a realistic proton density in the daughter nucleus. To obtain the single particle densities inside the daughter nucleus, we solve the Schrödinger equation for a spherical Woods-Saxon potential with a spin orbit and Coulomb terms [5]–[26]. Subsequently, the PSF part of the code is completed by performing the integrals over the electron phase factors constructed with the Dirac radial functions. The code exhibits an improved numerical accuracy for finding the electron w. f. and a better interpolation procedure for integrating the PSF final expressions, particularly at low electron energies. For the NME part, we employ a code similar to that from Ref. [13] for computing the double Gamow-Teller transitions, assuming the following effective nucleon-nucleon interactions: GXPF1A [28] for $ ^{48}{\rm{Ca}} $, JUN45 [29] for $ ^{76}{\rm{Ge}} $, and $ ^{82}{\rm{Se}} $ and gcn50:82 [30] for $ ^{130}{\rm{Te}} $ and $ ^{136}{\rm{Xe}} $.

The values for the products $ P^{2\nu} $ are presented in the third column of Table 1 for the five nuclei of experimental interest. With the values of $ g^{2\nu}_{A, {\rm eff}} $, 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\nu}_{A,{\rm exp}} $ values taken from Refs. [34–35], 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 percentages between the $ g_{A,{\rm eff}} $ values obtained within our calculations to reproduce the experimental DBD half-lives and those obtained by fitting the B(GT) experimental data, estimated in percentages, $ \epsilon = \left(g^{2\nu}-g^{2\nu}_{A,{\rm eff}}\right)/g^{2\nu} $. As observed, the two sets of values are close to each other, the smallest differences being in the case of $ ^{76}{\rm{Ge}} $ and $ ^{48}{\rm{Ca}} $ nuclei.

nucleus	$T_{1/2}^{2\nu}/\rm yr$	$P^{2\nu}/\rm yr^{-1}$	$g^{2\nu}_{A,{\rm eff}}/g^{2\nu}_{\rm exp}$	$\epsilon $(%)
48Ca	$6.40\times 10^{19}$ [31]	$123.81\times10^{-21}$	0.65/0.71[34]	8.45
76Ge	$1.92 \times 10^{21}$ [32]	$5.16 \times 10^{-21}$	0.56/0.60[35]	6.60
82Se	$0.92 \times 10^{20}$ [32]	$186.62 \times 10^{-21}$	0.49/0.60[35]	18.33
130Te	$8.20 \times 10^{20}$ [33]	$25.26 \times 10^{-21}$	0.47/0.57[35]	17.33
136Xe	$2.16 \times 10^{21}$ [31]	$20.30 \times 10^{-21}$	0.39/0.45[35]	13.33

Table 1.  Results for $2\nu\beta\beta $ decay mode.

Subsequently, we calculated the $ P^{0\nu} $ products in the case of the light Majorana neutrino exchange mechanism, with $ \langle\eta_l\rangle = \langle m_{\nu}\rangle/{m_e} $ and the light neutrino parameter defined as:

where $ U_{ek} $ denotes the first row elements of the Pontecorvo-Maki-Nakagawa-Sakata neutrino mixing matrix, and $ m_k $ depicts the light neutrino masses [36]. The expression of the nuclear matrix elements can be written as a sum of Gamow-Teller ($ GT $), Fermi ($ F $), and tensor ($ T $) components [9, 27]:

where $ M^{0 \nu}_{GT} $, $ M^{0 \nu}_F $, and $ M^{0 \nu}_T $ are these components. These are defined as follows:

$ O^\alpha_{mn} $ are transition operators ($ \alpha = GT,F,T $) and the summation covers all nucleon states. Correspondingly, the two-body transition operators $ O^{\alpha}_{12} $ can be expressed in factorized form as [27]:

where $ N_\alpha $ is a numerical factor including the coupling constants, and $ S_\alpha $, $ R_\alpha $, and $ C_\alpha $ 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^{\alpha}_{12} $ operators contain three components, namely the spin, center-of-mass, and the radial component, and the expectation values of the first two components are easily managed. The radial part is the most difficult to calculate, as it contains neutrino potentials written in different approximations, where the expectation values are double integrals over these potentials. Moreover, short-range correlations and finite nucleon size corrections are introduced in this part of the computation. The neutrino potentials depend weakly on intermediate states, and they 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.2{\rm fm} $, $ j_{0,2}(qr) $ are the spherical Bessel functions, and the integrals are over the neutrino exchange momentum $ {\it q} $. In our calculations, we use the closure approximation. $ \langle E_N\rangle $, as mentioned above, represents the average energy of the virtual states in the intermediate odd-odd nucleus included in the description of the decay. Furthermore, we note that the factor $ 2R $a is canceled by a similar factor from the denominator of the PSF expression, such that $ P^{0\nu} $ 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].

The vector and axial vectors form factors are specified as $ \Lambda_V = 850 \;{\rm{MeV}} $ and $ \Lambda_A = 1086 \;{\rm{MeV}} $ [2].

For computing the radial matrix elements $ \langle nl|H_{\alpha}|n'l'\rangle $, we use the harmonic oscillator $ {\rm HO} $ wave functions $ \psi_{nl}(lr) $ and $ \psi_{n^\prime l^\prime}(r) $ corrected by a factor $ [1 + f(r)] $, which takes into account the SRC effects induced by the nuclear interaction [27]:

For the correlation function, we take the functional form

For the $ a $, $ b $, and $ c $ constants we use the parametrization employed in Ref. [37].

Taking into account the HOC and FNS effects, the radial matrix elements of the neutrino potentials become:

We note that in the case of $ P^{0\nu} $ products, the axial-vector constant also enters the expressions of the neutrino potentials, in addition to the factor $ g_{A,{\rm eff}}^4 $, such that the half-life expression for the $ 0\nu\beta\beta $ decay, Eq. (5), contains a "double" dependence on this constant. Naturally, for coherence, the same values of the $ g_{A,{\rm eff}} $ constant should be assumed in both places, i.e., both in the $ {\it P^{0\nu}} $ and in the half-life computation. We note that these values may differ from the values of this constant used in the $ 2\nu\beta\beta $ decay mode. Because presently we do not know the correct value of $ g_{A,{\rm eff}} $ for the $ 0\nu\beta\beta $ decay mode, we calculate the $ P^{0\nu} $ products for the free nucleon value (1.273). As input parameters, the $ P^{0\nu} $ values can be easily computed for other effective values of the g_A,eff constant. The obtained values of the products $ P^{0\nu} $, in $ [\rm yr^{-1}] $ units, are presented in the third column of Table 2.

nucleus	$T_{1/2}^{0\nu} /\rm yr$	$P_\nu^{0\nu} /\rm yr^{-1}$	$\langle m_\nu\rangle$/eV
48Ca	$> 2.0 \times 10^{22}$ [31]	$7.30 \times 10^{-15}$	<26.49
76Ge	$> 8.0 \times 10^{25}$ [38]	$9.95 \times 10^{-15}$	<0.29
82Se	$> 3.6\times10^{23}$ [39]	$34.45 \times 10^{-15}$	<2.87
130Te	$> 4.0\times10^{24}$ [40]	$71.45 \times 10^{-15}$	<0.59
136Xe	$> 1.8\times10^{25}$ [41]	$71.01 \times 10^{-15}$	<0.28

Table 2.  Results for the $0\nu\beta\beta$ decay mode.

The values of the products $ P^{(2,0)\nu} $ from Tables 1 and 2, obtained with the approach described above, are very close to the computed values, obtained by multiplying the separately calculated values of NME and PSF. This is understandable, as in their calculation by the two methods, we employed the same values of input parameters and the same nuclear approximations and parametrizations. The small differences stem from the numerical precision of numerical codes we used. We emphasize, however, that the importance of our current approach is the elimination of the incoherence by using NME and PSF values calculated separately with different values for common nuclear parameters, which can introduce significant errors in the evaluation of $ {\rm NME} \times {\rm PSF} $ product as a whole. The errors in the evaluation of these products can indeed be significant if different values of $ g_A $ are assumed in the computation of $ {\rm NME} $ and $ {\rm PSF} $, and if these values are not the same as the value used in the $ g_A^4 $ 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. [27] for $ ^{48}{\rm{Ca}} $, $ ^{76}{\rm{Ge}} $, and $ ^{82}{\rm{Se}} $ nuclei, and found to be within 10%–14% (without the factor $ g_A^4 $). 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. Thus, relevant errors may occur in calculating the products $ {\rm NME} \times {\rm PSF} $ when the $ {\rm NME} $ and $ {\rm PSF} $ values are taken from separate calculations reported in the literature.

Subsequently, we revise the limits of the light neutrino mass parameter $ \langle m_{\nu}\rangle $ using our calculations and the most recent experimental limits reported for the $ 0\nu\beta\beta $ decay half-lives. These results are presented in the last column of Table 2. One observes that presently, the most stringent constraints on this parameter stem from the nuclei $ ^{76}{\rm{Ge}} $ and $ ^{136}{\rm {Xe}} $, because of both the experimental results (the lowest limits measured at present for the $ 0\nu\beta\beta $ decay half-lives [38], [41]) and the 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 remains unknown, for accurate half-life predictions and constrains of beyond SM parameters, the goal is to reduce the errors associated with this constant. One suggestion is to apply information from different decay modes and/or from DBD experiments on different nuclei. For example, for a particular nucleus, the ratio of the $ 2\nu $ and $ 0\nu $ half-life expressions reads:

As seen from the above formula, any information that we obtain regarding the relative magnitude of the $ g_{A, {\rm eff}} $ values for the $ 2\nu $ and $ 0\nu $ 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. Further, referring to two different nuclei, denoted with $ m $ and $ n $, the ratio of their half-lives reads:

For the $ 2\nu\beta\beta $ decay mode:

For the $ 0\nu\beta\beta $ decay mode:

Therefore, the $ g^{2\nu}_{A,{\rm eff}} $ for one particular nucleus can be deduced if one knows with (more) precision the value of this constant for another nucleus, using the experimental half-lives and the calculated $ P^{2\nu} $ for both nuclei. For example, one can take advantage of the possible experimental determination of this parameter for some particular isotopes, as recently proposed in Ref. [42]. Similar considerations, i.e., the exploitation of data from several experiments, are valid for predicting $ 2\nu\beta\beta $ decay half-life for a nucleus that was not yet measured, if accurate data for another nucleus and good estimations of the $ g^{2\nu}_{A,{\rm eff}} $ value from other (non-DBD) experimental data are available. For such predictions, information on DBD half-lives not yet measured obtained from empirical formulas, as proposed in Ref. [43], is valuable.

Similarly, for the $ 0\nu\beta\beta $ decay, more information about the effective value of $ g^{0\nu}_A $ can be deduced for a particular nucleus if we this value is known for another nucleus. For example, we might know $ g^{0\nu}_{A,{\rm eff}} $ with more precision in the case of nuclei where the 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}{\rm Mo} $ case), and where $ g^{2\nu}_{A,{\rm eff}} $ and $ g^{0\nu}_{A,{\rm eff}} $ might have close values.

We proposed a new approach for calculating the NME and PSF for DBD, by directly computing their product. The product as a whole can be computed more consistently, with a unique dependence on some parameters that were previously entered separately in the NME and PSF expressions, which thus assumed single values. The values of the product are given in the same units as $ T_{1/2}^{-1} $ (i.e. [$ \rm yr^{-1}] $), removing any possible confusion by using the theoretical calculations for interpretation of DBD data. The new codes for calculating $ {\rm NME} \times {\rm PSF} $ product include improved routines used in our previous studies for the separate computation of these two quantities. We provide the values of these products for the $ 2\nu $ and $ 0\nu $ DBD modes for the five nuclei of most experimental interest. Subsequently, using our calculations and the newest half-life values for the $ 0\nu\beta\beta $ decays reported in literature, we revise the upper limits for the light mass neutrino parameter. In the half-life formulas, we separate the strong dependence on the axial-vector constant, i.e., the factor $ g_A^4 $, which brings a large uncertainty in the calculation and suggest some ways to reduce/avoid the errors related to the uncertain value of this constant. This can be done by employing ratios of $ g^{(2,0)\nu}_{A, {\rm eff}} $ and $ P^{(2,0)\nu} $ (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 of experimentally determining this constant for some particular isotopes. We hope that our work proves to be a step forward for more consistent DBD calculations, which will aid in predicting and interpreting experimental data.