Phase Space Factors for Double-Beta Decays

Double-beta decay is presently a very studied process both theoretically and experimentally due to its potential to provide valuable information about important, but still unknown issues related to the neutrino properties and conservation of some symmetries. In the theoretical study of the double-beta decay two key quantities entering the half-life formulas are important, namely the phase space factors embedding the influence of the Coulomb field of the daughter nucleus on the emitted electrons/positrons, and the nuclear matrix elements embedding the nuclear structure effects of the nuclei participating in the decay. Accurate calculation of both of them are needed for good predictions of the double-beta decay half-lives and transitions still unmeasured, and for constraining various beyond Standard Model parameters associated with mechanisms that may contribute to the neutrinoless double-beta decay modes. During time much attention has been paid to the nuclear matrix elements that were considered to bring the largest uncertainties in the computation of the double-beta decay half-lives, while the phase space factors were considered until the recent past to be computed with enough precision. However, newer computation of the phase space factors performed with more precise methods revealed relevant deviations from their values reported previously, especially for heavier nuclei and for positron emitting and electron capture decay modes. In this paper we review the progress made in the computation of the phase space factors for double beta decay. We begin with the non-relativistic approaches, continue with the relativistic approaches which use approximate electron/positron wave functions, and end up with recent, more precise, computations of the phase space factors where exact electron wave functions are obtained from the resolution of a Dirac equation in a Coulomb-type potential and with inclusion of finite nuclear size and screening effects. We report an up-dated and complete list of the phase space factors (PSF) for the following DBD modes: β⁻β⁻, β⁺β⁺, ECβ⁺, and ECEC and for transitions to final ground and first excited 2⁺ and 0⁺ states of the daughter nuclei. We also make a comparison between different values of the phase space factors found in literature and discuss the differences between these results.

1. Introduction

Nuclear double beta decay (DBD) is the process with the longest lifetime measured so far, by which an even-even nucleus decays naturally into another even-even nucleus with the same atomic mass A but with the electric charge changed by two units. Within the Standard Model (SM) it can occur through several decay modes with lepton number conservation, namely with the emission of two electrons/positrons and two neutrinos/anti-neutrinos. However, theories beyond SM predict that this process may also occur without emission of neutrinos/anti-neutrinos, namely with lepton number violation (LNV) by two units, and these decay modes are generically called neutrinoless double-beta decays. The DBD modes can be classified according to the number and type of the released leptons, and can be divided in two categories:

1) Decays occurring with LNC

1a) Two neutrino double-electron decay (2νβ⁻ β⁻)

$\begin{array}{ll}(A,Z)\to (A,Z+2)+2e^{-}+2\overline{\nu }& (1)\end{array}$

1b) Two neutrino double-positron decay

$\begin{array}{ll}(A,Z)\to (A,Z-2)+2e^{+}+2\nu & (2)\end{array}$

1c) Two neutrino electron capture (EC) positron emitting decay (2νECβ⁺)

$\begin{array}{ll}{e}^{-}+(A,Z)\to (A,Z-2)+e^{+}+2\nu & (3)\end{array}$

1d) Two neutrino double electron capture decay (2νECEC)

$\begin{array}{ll}2e^{-}+(A,Z)\to (A,Z-2)+2\nu & (4)\end{array}$

2) Decays occurring with LNV by two units

2a) Neutrinoless double-electron decay (0νβ⁻ β⁻)

$\begin{array}{ll}(A,Z)\to (A,Z+2)+2e^{-}& (5)\end{array}$

2b) Neutrinoless double-positron decay (0νβ⁺ β⁺)

$\begin{array}{ll}(A,Z)\to (A,Z-2)+2e^{+}& (6)\end{array}$

2c) Neutrinoless electron capture positron emitting decay (0νECβ⁺)

$\begin{array}{ll}{e}^{-}+(A,Z)\to (A,Z-2)+e^{+}& (7)\end{array}$

2d) Neutrinoless double electron capture decay (0νECEC)

$\begin{array}{ll}2e^{-}+(A,Z)\to (A,Z-2)& (8)\end{array}$

The neutrinoless DBD presents particularly a great interest because it can provide us with information about some conservation laws (lepton number, Charge-Parity, Lorentz symmetry) and about fundamental properties of neutrino as their character (Dirac or Majorana), absolute mass and mass hierarchy, number of flavors (are there other neutrino species than the three ones known until present?), etc.

Complete information about the progress in the study of DBD can be obtained from several excellent reviews written during time given. See for examples the Primakoff and Rosen [1], Haxton and Stephenson [2], Doi et al. [3], Tomoda [4], Suhonen and Civitarese [5], Avignone et al. [6], Ejiri [7], Rodejohann [8], Vogel [9], Vergados et al. [10], which, in turn, contain a comprehensive list of other useful references in domain. The theoretical study of DBD implies firstly the derivation of the half-lives for the above mentioned decay modes and transitions to ground (0⁺) and excited 2⁺ and 0⁺ states of the daughter nucleus. In a good approximation the half-life expressions can be written as products of the phase space factors (PSF) embedding the influence of the Coulomb field of the nucleus on the emitted electrons/positrons and the nuclear matrix elements (NME) embedding the nuclear structure effects of the nuclei participating in the decay. In the case of neutrinoless DBD an additional factor related to the specific mechanism that can contribute to the occurrence of these decays enters as well [11–13]:

$\begin{array}{l}{\left(T_{1/2}^{2\nu }\right)}^{-1}=G_{2\nu }^n(E_0,Z)g_A^4\mid m_e{c}^2{M}_{2\nu }{\mid }^2\\ {\left(T_{1/2}^{0\nu }\right)}^{-1}=G_{0\nu }^n(E_0,Z)g_A^4\mid M_{0\nu }{\mid }^2{\left(\frac{<m_{\nu }>}{m_e}\right)}^2& (9)\end{array}$

where $G_{2\nu (0\nu )}^n$ and M_2ν(0ν) are the phase space factors and matrix elements for the 2ν (0ν) decay modes, m_e and < m _ν > are the electron and effective neutrino masses, g_A is the axial vector coupling constant and n denotes one the DBD modes (1–2). For the 2νββ decay mode the nuclear matrix elements are only on the Gamow-Teller type (GT) (are given by the GT operator), while for the 0νββ they are a sum of GT, Fermi (F), and tensor contributions:

$\begin{array}{ll}{M}^{0\nu }=M_{GT}^{0\nu }-{\left(\frac{g_V}{g_A}\right)}^2·M_F^{0\nu }-M_T^{0\nu }& (10)\end{array}$

The half-life expressions from above are written such as the “nuclear” parts (NME) are dimensionless quantities while the “atomic” parts (PSF) are given in [yr⁻1] units. Also, it worth to mention that the PSF are quantities which do not depend on g_A, while the NME depend on this constant through the (g_V/g_A) factor. Much attention has been paid over time to the NME calculation, since it was considered to bring the largest uncertainties in the computation of the DBD half-lives. There are several nuclear methods used more often to calculate the NME, including the quasiparticle random phase approximation (QRPA), the interacting shell model, the interacting boson model, the generator coordinate method, and the projected Hartree-Fock Bogolibov model, etc. These uncertainties come mainly from the choice of the model spaces and type of correlations taken into account in calculation. The results and the discrepancies between the NME obtained with these methods and possible sources of uncertainties in calculations have been extensively discussed in the literature. For the interested reader we suggest a number of references that include the most important results [14–31]. The PSF were considered until the recent past to be computed with enough precision. However, newer computation of these PSF performed with more precise methods revealed relevant deviations from the earlier calculations, especially for heavier nuclei and for positron emitting and EC decay modes. Accurate PSF calculations are very needed now in the study of DBD processes where precision measurements of electron spectra and correlation functions are under way, since they may bring useful information, as for example about possible violation of the Lorentz symmetry or about the effective value of the g_A constant.

The purpose of this work is to give a review on the PSF calculations for all decay modes and transitions to final ground and first 0⁺ and 2⁺ excited states and to give a complete list of the calculated PSF values (however, we except the decay mode 0νECEC that can not occur to the order of approximation that is presently considered in literature). We begin with the presentation of the earliest calculations using non-relativistic approaches and simple expression of the Fermi factors that allows the analytical computation of the PSF. Then, we continue with relativistic approaches, first with those that use approximate analytical expressions of the Fermi factors, and then with more advanced methods where the Coulomb distortion of the electron wave functions are obtained by solving numerically the Dirac equation with inclusion of the finite nuclear size (FNS) and electron screening effects. In the main part of the review we present our recent PSF calculations where the following ingredients are included in addition to other (recent) similar calculations: (i) the use of a Coulomb potential derived from a realistic proton density distribution in nucleus; (ii) development of new routines with improved numerical precision both for solving the Dirac equations and integrating the PSF expressions, and (iii) inclusion in calculation the most recent Q^ββ values in Wang et al. [32]. Also, we present the first calculations of the PSF values for transitions to final excited $2_1^{+}$ states computed with exact Dirac electron functions. Then, comparing different PSF calculations found in literature, performed with relativistic approaches, we find with few exceptions an agreement within 10% between the PSF values calculated with most recent methods and those performed with approximate electron functions for β⁻β⁻ and β⁺β⁺ decays of lighter nuclei. For heavier nuclei and particularly for EC decay modes we got more and larger discrepancies, whose possible sources are discussed as well.

2. Early Approaches of the Phase Space Factors

The phase space factors represent the distortion of the electron plane waves in the Coulomb field of the nucleus. Thus, to compute PSF for DBD modes we need first to obtain the wave functions of the electron(s)/positron(s) emitted or electron(s) captured in the decay, which are distorted by the Coulomb-type potential of the nucleus. In the first derivations of the DBD half-life formulas the electron wave functions were corrected by multiplying them with the so-called Coulomb functions (or Fermi factors) F(Z, ϵ) obtained in a non-relativistic approximation as the square of the ratio of the solution of a Schrödinger equation in a Coulomb potential given by a point charge Z to a plane wave, evaluated at the origin [2]. They are expressed in an analytical form as follows:

$\begin{array}{ll}{F}^{NR}(Z,ϵ)=\frac{2\pi \eta }{1-e^{-2\pi \eta }}=e^{\pi \eta }\mid \Gamma (1+i\eta ){\mid }^2& (11)\end{array}$

where η = ±αZβ, β = p/ϵ, with p = |p| and ϵ the electron momentum and energy. The +(−) signs are taken for outgoing electrons (positrons).

In their calculations of the DBD half-lives, Primakoff and Rosen [1] used an even simpler expression of these factors obtained from Equation (7) by taking β = 1 in the exponent of the above expression. This permits an analytic evaluation of the phase space integrals resulting in suggestive formulas for DBD half-lives, namely proportional to seventh through eleventh powers of the Q^ββ values. This approximation is an acceptable one, especially for the DBD of nuclei with low Z with emission of two electrons (β⁻β⁻) as it was shown Haxton and Stephenson [2]. Then, going to improve the accuracy in estimating the PSF, in the simplest relativistic approach the Fermi factors are given by the square of the ratio of the solution of a Dirac equation in a Coulomb potential given by a point charge Z to a plane wave, evaluated at the nuclear surface [33].

$\begin{array}{ll}{F}^R(Z,ϵ)=2(1+\gamma ){(2p{R}_0)}^{2(\gamma -1)}{e}^{\pi \eta }\mid \frac{\Gamma (\gamma +i\eta )}{\Gamma (2\gamma +1)}{\mid }^2& (12)\end{array}$

with γ = [1 − (αZ)²]1/2 and R₀ the nuclear radius. Further, a better estimation should include the screening and finite nuclear size effects. This can be done by multiplying the Fermi factors of Equation (8) by a factor L₀ that can be taken from Behrens and Janecke [34]. For electrons L₀ is considered close to unity for the energies and nuclei of interest, so it is not taken into account in these DBD calculations. For positrons screening enhancements at small β can lead to values substantially larger than unity [34]. To see how large can be the relativistic correction we take the ratio of the two Fermi factors [2]:

$\begin{array}{ll}\frac{F^{NR}(Z,ϵ)}{F^R(Z,ϵ)}=\frac{1}{L_0}\times \frac{1}{{(2p{R}_0)}^{2(\gamma -1)}}\times \frac{\Gamma {(2\gamma +1)}^2}{4\gamma }\times \mid \frac{\Gamma (1+i\eta )}{\Gamma (\gamma +i\eta )}{\mid }^2& (13)\end{array}$

First, one observes that this expression is invariant to the change η → −η, so one expects relativistic corrections to be similar for electron and positron emitting modes. As it is discussed in Haxton [33], in the limit of small Z this ratio is approaching to unit, thus the non-relativistic correction to the electron/positron wave functions used for computing PSF for light nuclei can be a good approximation. However, when we are dealing with DBD nuclei, the difference between PSF values calculated with non-relativistic and relativistic Fermi factors can be quite large. For example, in the case of ¹³⁰Te → ¹³⁰Xe decay the use of non-relativistic Fermi factors instead of relativistic ones leads to underestimation of the DBD half-life by a factor 5 [33], which increases to a factor 200 for the decay ²³⁸U → ²³⁸Pu. The same expression for the Fermi factors as Equation (8) but in the limit (αZ)²≪1 are used by Suhonen and Civitarese [5] when they show how one calculates the PSF for DBD.

Suhonen and Civitarese have also used these expressions in their DBD calculations presented in their DBD review [5]. Also, Further, Tomoda used another approach in Tomoda [4], obtaining the Fermi factors necessary for PSF computation as combinations of radial solutions of the Dirac equation in spherical coordinates in a Coulomb potential. The electron radial wave functions are expanded in powers of r and the leading constant terms are retained. In this approximation they estimate the error to be $\frac{1}{6}{(p_{eff}{R}_0)}^2$, where p_eff is the “effective” momentum of the leptons at the nuclear surface (p_eff = 3Zα/2R₀ for electrons and p_eff = k ~ 1MeV for neutrinos). The radial functions are determined by boundary conditions, namely they have to be finite at r = 0 and become a plane wave plus an incoming spherical wave at infinity. Doi et al. also consider analytical (approximate) Fermi factors in their papers [3, 35–37], with similar expressions as those from Equation (8), but also, without inclusion of the screening effects.

3. Recent Approaches

Kotila and Iachello recalculated the PSF by using exact solutions of the relativistic Dirac equation for Coulomb fields modified by the screening effect [11, 12]. The radial electron/positron wave function were calculated numerically by using the subroutines package RADIAL [38]. The input in the package is the potential V which is primarily the Coulomb potential of the daughter nucleus with charge Z_d, V(r) = −Z_d(αℏc)/r. The FNS effect is taken into account with a Coulomb potential given by an uniform charge distribution in a sphere of radius $R_A=r_0{A}^{1/3}$ with r₀ = 1.2 fm [11, 37]:

$\begin{array}{ll}V(Z,r)=\{\begin{array}{ll}-\frac{Z\alpha \hslash c}{r},\hfill & r\ge R_A,\hfill \\ -Z(\alpha \hslash c)\left(\frac{3-{(r/R_A)}^2}{2R}\right),\hfill & r<R_A,\hfill \end{array}& (14)\end{array}$

The screening effect was taken into account by multiplying the above expression of V(r) with a function ϕ(r), solution of the Thomas Fermi equation. The procedure will be described later. This more rigorous relativistic treatment, with the inclusion of FNS and screening effects was also adopted by us in Stoica and Mirea [13],Mirea [39], but with the following improvements: (i) the use of a Coulomb-type potential, obtained from a realistic distribution of the protons in the nucleus, (ii) building up new, more precise and under control, routines; (iii) using more recent Q_ββ values. In the following, we review our works from Stoica and Mirea [13] and Mirea [39].

4. Theoretical Treatment

The first step in the calculation of the PSF is the determination of the wave function of the electron in the Coulomb field of the nucleus.

4.1. Electron Radial Functions

For free states we use relativistic scattering electron/positron wave functions, solutions of the Dirac equation in a central (Coulomb) potential:

$\begin{array}{ll}{\Psi }_{ϵ\kappa \mu }^{+}(r)=\left(\begin{array}{l}{g}_{\kappa }(ϵ,r){\chi }_{\kappa }^{\mu }\hfill \\ i{f}_{\kappa }(ϵ,r){\chi }_{-\kappa }^{\mu }\hfill \end{array}\right) & (15)\end{array}$

for β⁻ decay and

$\begin{array}{ll}{\Psi }_{ϵ\kappa \mu }^{-}=\left(\begin{array}{l}i{f}_{\kappa }(ϵ,r){\chi }_{-\kappa }^{-\mu }\hfill \\ -g_{\kappa }(ϵ,r){\chi }_{\kappa }^{-\mu }\hfill \end{array}\right)& (16)\end{array}$

for β⁺ decay where κ = (l − j)(2j+1) is the relativistic quantum number and ${\chi }_{\kappa }^{\mu }$ are spherical spinors. The quantities g_κ(ϵ, r) and f_κ(ϵ, r) are the large and small components of the radial wave functions which satisfy the radial equations:

$\begin{array}{l}\frac{d{g}_{\kappa }(ϵ,r)}{dr}=-\frac{\kappa +1}{r}{g}_{\kappa }(ϵ,r)+\frac{ϵ-V+m_e{c}^2}{c\hslash }{f}_{\kappa }(ϵ,r)\\ \frac{d{f}_{\kappa }(ϵ,r)}{dr}=-\frac{ϵ-V-m_e{c}^2}{c\hslash }{g}_{\kappa }(ϵ,r)+\frac{\kappa -1}{r}{f}_{\kappa }(ϵ,r)& (17)\end{array}$

where V can be negative/positive (for β⁻/β⁺). For the continuum spectrum, these functions are normalized so that they have the asymptotic behavior of the Coulomb wave functions:

$\begin{array}{ll}{\left(\begin{array}{l}{g}_k(ϵ,r)\hfill \\ f_k(ϵ,r)\hfill \end{array}\right)}_{r\to \infty }~\frac{\hslash e^{-i{\delta }_k}}{pr}\left(\begin{array}{l}\sqrt{\frac{ϵ+m_e{c}^2}{2ϵ}}\sin (kr-l\frac{\pi }{2}-\eta \ln (2kr)+{\delta }_k)\hfill \\ \sqrt{\frac{ϵ-m_e{c}^2}{2ϵ}}\cos (kr-l\frac{\pi }{2}-\eta \ln (2kr)+{\delta }_k)\hfill \end{array}\right)& (18)\end{array}$

Here c is the light velocity, m_e/ϵ are the electron mass/energy, $k=p/\hslash =\sqrt{{ϵ}^2-{(m_e{c}^2)}^2}/(c\hslash )$ is the electron wave number, $\eta =Z_{eff}{e}^2/\hslash v$, is the asymptotic Sommerfeld parameter, δ_κ is the global phase shift and V is the Coulomb interaction energy between the electron and the daughter nucleus. For pure Coulomb fields (no screening), the effective charge number is Z_eff = ±Z for β^∓ decays. The energy is $ϵ=\sqrt{{(m_e{c}^2)}^2+{(pc)}^2}$.

The phase shifts are obtained by matching the inner numerical solution to the analytic function. The global phase shift δ = δ_C + δ_k [37] is obtained by renormalizing the numerical solution of the Dirac equation in the asymptotic region. The global phase shift depends on the electron energy. It should be taken into consideration when the electron phase factors are calculated [40]. It includes the Coulomb phase shift ${\delta }_C=arg\left[\zeta (ϵ+m_0{c}^2)-i(\kappa +\lambda )k\hslash c\right]-(\lambda -l-1)\frac{\pi }{2}+arg\Gamma (\lambda -i\eta )-S_{\zeta ,\kappa }\pi$, with $\lambda =\sqrt{{\kappa }^2-{\zeta }^2}$, ζ = ±Z_effα, and the shift due to distorsion of the Coulomb potential δ_k [38, 41]. The phase shift due to the distorsion of the Coulomb potential is calculated numerically during the renormalization of the numerical solutions to a sinusoidal form given by the Equation (18). The normalization procedure is described in Salvat [38] and Salvat and Mayol [42].

The bound states wave functions for the electron

$\begin{array}{ll}{\Psi }_{{ϵ}_n\kappa \mu }^b(r)=\left(\begin{array}{l}{g}_{n,\kappa }^b(r){\chi }_{\kappa }^{\mu }\hfill \\ i{f}_{n,\kappa }^b(r){\chi }_{-\kappa }^{\mu }\hfill \end{array}\right)& (19)\end{array}$

are solutions of the Dirac equation (17) and correspond to the eigenvalues ϵ_n (n is the radial quantum number). The quantum number κ is related to the total angular momentum j_κ = ∣κ∣−1/2. For simplicity, the quantities $g_{n,\kappa }(r)=r{g}_{n,\kappa }^b(r)$ and $f_{n,\kappa }(r)=r{f}_{n,\kappa }^b(r)$ are used in the following. These wave functions are normalized such that

$\begin{array}{ll}{\displaystyle {\int }_0^{\infty }[g_{n,\kappa }^2(r)+f_{n,\kappa }^2(r)]dr=1.}& (20)\end{array}$

An asymptotic solution is obtained by means of the WKB approximation and by considering that the potential V is negligible small:

$\begin{array}{ll}\frac{f_{n,\kappa }}{g_{n,\kappa }}=\frac{c\hslash }{ϵ+m_e{c}^2}\left(\frac{{g^{\prime }}_{n,\kappa }}{g_{n,\kappa }}+\frac{\kappa }{r}\right)& (21)\end{array}$

where

$\begin{array}{ll}\frac{{g^{\prime }}_{n,\kappa }}{g_{n,\kappa }}=-\frac{1}{2}{\mu }^{\prime }{\mu }^{-1}-\mu & (22)\end{array}$

with

$\begin{array}{ll}\mu ={\left[\frac{ϵ+m_e{c}^2}{{\hslash }^2{c}^2}(V-ϵ+m_e{c}^2)+\frac{{\kappa }^2}{r^2}\right]}^{1/2}.& (23)\end{array}$

In our calculations we use n=0 and n=1 number of nodes, for the orbitals 1s_1/2 and 2s_1/2, κ being -1. The eigenvalues of the discrete spectrum are obtained by matching two numerical solutions of the Dirac equation: the inverse solution that starts from the asymptotic conditions and the direct one that starts at r=0.

4.2. The Coulomb Potential

The influence of the nuclear structure was taken into consideration in Stoica and Mirea [13], Mirea [39]. by using a potential V(r) derived from a realistic proton density distribution in the nucleus. This proton density in the nucleus is obtained by solving the Schrödinger equation for a Woods-Saxon potential. In this case the potential becomes:

$\begin{array}{ll}V(Z,r)=\alpha \hslash c{\displaystyle \int \frac{{\rho }_e(\overrightarrow{r^{\prime }})}{\mid \overrightarrow{r}-\overrightarrow{r^{\prime }}\mid }}d\overrightarrow{r^{\prime }}& (24)\end{array}$

where the charge density is

$\begin{array}{ll}{\rho }_e(\overrightarrow{r})={\displaystyle \sum _i(2j_i+1)v_i^2}{\left|{\Psi }_i(\overrightarrow{r})\right|}^2& (25)\end{array}$

Ψ_i is the proton (Woods-Saxon) wave function of the spherical single particle state i and v_i is its occupation amplitude. The factor (2j_i + 1) reflects the spin degeneracy.

As an example, the difference between the behavior of the constant charge density ρ_e and the realistic charge density is displayed in Figure 1 for the daughter nucleus ¹⁵⁰Nd. We computed the Coulomb potential with formula (24). In this case, the differences given by the charge densities are translated in a shift of 0.5 MeV energy in the potential at r = 0, as it can be observed in Figure 2. This difference in energy vanishes when r increases, but is able to affect the values of the wave functions.

FIGURE 1

Figure 1. (A) Realistic proton density ρ_e for ¹⁵⁰Nd represented in cylindrical coordinates ρ and z. (B) Profile of the realistic proton density ρ_e for ¹⁵⁰Nd (thick line) compared with that given with the constant density approximation (dot-dashed line).

FIGURE 2

Figure 2. The Coulomb potential V obtained by using the realistic charge density is plotted with a thick full line as function of the radius r. The Coulomb potential obtained with the finite nuclear size approximation is displayed with a thin line. The calculations were performed for the ¹⁵⁰Nd.

The Coulomb field felt by the β particle is damped by the electron cloud. A simple way to take into account the screening effect is to multiply the potential V(r) with an universal screening function ϕ(r), which is the solution of the Thomas Fermi equation:

$\begin{array}{ll}{d}^2\varphi /d{x}^2={\varphi }^{3/2}/\sqrt{x},& (26)\end{array}$

with x = r/b, $b\approx 0.8853a_0{Z}^{-1/3}$ and a₀ = Bohr radius. A solution of this equation having the boundaries conditions ϕ(0) = 1 and ϕ(∞) = 0 was calculated within the Majorana method in Esposito [43]. A solution of the Thomas-Fermi equation was also used Kotila and Iachello [11] to take into account the screening effect. The Equation (26) does not depend on physical constants or on the charge of the parent nucleus, being an universal equation. The condition ϕ(0) = 1 means that the Coulomb potential is not perturbed by the electron cloud in the nuclear region. Asymptotically, ϕ(∞) = 0 means that the Coulomb field has completely disappeared. Therefore, the modality in which the universal screening function modifies the Coulomb potential should be adapted to the boundary conditions related to our specific mechanisms. For example, in the case of the β⁻β⁻ decay, the asymptotic boundary condition should be ${\varphi }_{{\beta }^{-}{\beta }^{-}}(\infty )=2/Z_d$, Z_d = Z₀−2 being the charge number of the daughter. That is, the Coulomb field is suppressed by the screening such that an interaction between β⁻-particles and an positive ion of charge 2⁺ is obtained asymptotically. Therefore, we use different methods to introduce the screening accordingly to the universal function ϕ(x).

In the case of the β⁻β⁻ process, the potential used to obtain the electron wave functions is:

$\begin{array}{ll}r{V}_{{\beta }^{-}{\beta }^{-}}(Z,r)=(rV(Z,r)+2)\times \varphi (r)-2& (27)\end{array}$

to take into account the fact that DBD releases a final positive ion with charge +2. V(Z, r) is negative. In the relations related to screening we use the atomic units, that is, the product αℏc=1. The asymptotic potential between an electron and a double ionized atom is $r{V}_{{\beta }^{-}{\beta }^{-}}=-2$. In this case, the charge number Z = Z₀+2 corresponds to the daughter nucleus, Z₀ being the charge number of the parent nucleus. In the case of the β⁺β⁺ process, the potential used to obtain the positron wave functions is:

$\begin{array}{ll}r{V}_{{\beta }^{+}{\beta }^{+}}(Z,r)=(rV(Z,r)+2)\times \varphi (r)-2& (28)\end{array}$

where the final configuration is characterized by an ion with charge -2. V(Z, r) is positive. In this case, the daughter nucleus has the charge number Z = Z₀ − 2. In both approaches, at r = 0 the potential is unscreened because ϕ(0) = 1. Asymptotically ϕ(r) tends to 0 and we are left with the charge number of the final system. In the ECβ⁺ reaction, a simultaneous combination of two processes is inferred. In the case of the EC process, the potential used to obtain the electron wave functions reads:

$\begin{array}{ll}r{V}_{EC}(Z,r)=(rV(Z,r)+1)\times \varphi (r)-1& (29)\end{array}$

and the charge number Z = Z₀ corresponds to the parent nucleus. V(Z, r) is negative. In the case of the β⁺ process, the potential used to obtain the positron wave functions reads:

$\begin{array}{ll}r{V}_{{\beta }^{+}}(Z,r)=(rV(Z,r)+1)\times \varphi (r)-1& (30)\end{array}$

to take into account that in the final configuration we have an ion with charge −1. V(Z, r) is positive. In this case the daughter nucleus has the charge number Z = Z₀ − 1. In the case of the ECEC process, the potential used to obtain the electron wave function is:

$\begin{array}{ll}{V}_{ECEC}(Z,r)=V(Z,r)\times \varphi (r)& (31)\end{array}$

and Z = Z₀, the final system being neutral. In Figure 3, the effective screening ϕ_eff given by the ratios $V_{{\beta }^{-}{\beta }^{-}}(x)/V(x)$, $V_{{\beta }^{+}{\beta }^{+}}(x)/V(x)$, $V_{{\beta }^{+}}(x)/V(x)$, V_EC(x)/V(x) and V_ECEC(x)/V(x). The behavior of our effective screenings agree to that obtained in Kotila and Iachello [11] where the Thomas-Fermi equation was solved for different boundaries.

FIGURE 3

Figure 3. The effective screening function ϕ_eff as function of the dimensionless parameter x calculated for the ¹³⁶Xe parent nucleus. The thin full line corresponds to the β⁻β⁻ decay, the thick dashed line is for the EC process, the thick full line is for the ECEC process, the dot-dashed thick line is for the β⁺ process and the dotted thick line is for the β⁺β⁺ decay.

4.3. Phase Space Factors Formulas

As given in Primakoff and Rosen [1], Doi et al. [44], Haxton et al. [33], and Tomoda [4], the differential ββ decay rate for 0⁺ → 0⁺ transitions reads:

$\begin{array}{ll}d{W}_{2\nu }=(a^{(0)}+a^{(1)}\cos {\theta }_{12})w_{2\nu }d{\omega }_{1}d{\omega }_{2}d{ϵ}_{1}d{ϵ}_{2}d\cos {\theta }_{12}& (32)\end{array}$

where θ₁₂ is the angle between the electrons. Both parameters a₀ and a₁ represents summation over intermediate nuclear states of products between electron phase factors and nuclear matrix elements summed over intermediate states N

$\begin{array}{ll}{a}^{(0)}=\frac{1}{4}{f}_{11}^{(0)}\left({\left|\text{}{\displaystyle \sum _N{M}_N(K_N+L_N)}\right|}^2+\frac{1}{3}{\left|\text{}{\displaystyle \sum _N{M}_N(K_N-L_N)}\right|}^2\text{}\right)& (33)\end{array}$

$\begin{array}{ll}{a}^{(1)}=\frac{1}{4}{f}_{11}^{(1)}\left(\text{}{\left|{\displaystyle \sum _N{M}_N}(K_N+L_N)\right|}^2-\frac{1}{9}{\left|{\displaystyle \sum _N{M}_N}(K_N-L_N)\right|}^2\text{}\right)& (34)\end{array}$

Here, $f_{11}^{(i)}$ represent the electron phase factors and the matrix elements between initial (I) and final (F) states are denoted

$\begin{array}{ll}{M}_N=\langle 0_F^{+}\Vert {\tau }^{+}\sigma \Vert 1_N^{+}\rangle \langle 1_N^{+}\Vert {\tau }^{+}\sigma \Vert 0_I^{+}\rangle & (35)\end{array}$

Also,

$\begin{array}{ll}{w}_{2\nu }=\frac{{(g_{A}G\cos {\theta }_C)}^4}{64{\pi }^7\hslash }{\omega }_1^2{\omega }_2^2(p_{1}c)(p_{2}c){ϵ}_1{ϵ}_2\delta ({ϵ}_1+{ϵ}_2+{\omega }_1+{\omega }_2-W_0)& (36)\end{array}$

with $W_0=Q_{\beta \beta }+2m_0{c}^2$. The previous expression depends on the electron energies ϵ₁, ϵ₂, electron momenta p₁, p₂ and on the neutrino energies ω₁, ω₂. Usual values of the parameters can be considered as G = 1.16637 × 10⁻⁵ GeV⁻² (the Fermi coupling constant), cos θ_C=0.9737 (θ_C standing for the Cabbido-Kobayashi-Maskawa mixing angle between d and s quarks) and g_A ≈1 (the axial vector form factor).

By considering that the excitations energies are of the order of giant Gamow-Teller resonance à = 1.12A^1/2 MeV, the closure approximation is adopted. That is, the energy of the intermediate state E_N is replaced with an average value 〈E_N〉. Therefore the energy denominators

$\begin{array}{ll}{K}_N=\frac{1}{{ϵ}_1+{\omega }_1+E_N-E_I}+\frac{1}{{ϵ}_2+{\omega }_2+E_N-E_I}& (37)\end{array}$

$\begin{array}{ll}{L}_N=\frac{1}{{ϵ}_1+{\omega }_2+E_N-E_I}+\frac{1}{{ϵ}_2+{\omega }_1+E_N-E_I}& (38)\end{array}$

are changed within averaged values

$\begin{array}{ll}{K}_N\approx \langle K_N\rangle \frac{\frac{1}{2}{W}_0+\langle E_N\rangle -E_I}{\frac{1}{2}{W}_0+E_N-E_I}& (39)\end{array}$

$\begin{array}{ll}{L}_N\approx \langle L_N\rangle \frac{\frac{1}{2}{W}_0+\langle E_N\rangle -E_I}{\frac{1}{2}{W}_0+E_N-E_I}& (40)\end{array}$

In this way, a factorization of the Equation (32) can be realized. By considering $Ã\approx \frac{1}{2}{W}_0+\langle E_N\rangle -E_I$ and the nuclear matrix element on the form

$\begin{array}{ll}{M}_{2\nu }={\displaystyle \sum _N\frac{\langle 0_F^{+}\Vert {\tau }^{+}\sigma \Vert 1_N^{+}\rangle \langle 1_N^{+}\Vert {\tau }^{+}\sigma \Vert 0_I^{+}\rangle }{\frac{1}{2}{W}_0+E_N-E_I}}& (41)\end{array}$

the half-live obtained with Equation (32) can be written as a product between phase space factors an nuclear matrix elements. By integrating Equation (32), two phase space factors are obtained [4, 11]:

$\begin{array}{ll}\frac{d{W}_{2\nu }}{d\cos {\theta }_{12}}=|M_{2\nu }{|}^2{F}_{2\nu }^{(0)}+|M_{2\nu }{|}^2{F}_{2\nu }^{(1)}\cos {\theta }_{12}& (42)\end{array}$

$\begin{array}{l}{F}_{2\nu }^{(0)}=\frac{2{\tilde{A}}^2}{3\ln 2}\text{}{\displaystyle {\int }_{m_e{c}^2}^{Q_{\beta \beta }+m_e{c}^2}\text{}}\text{}d{ϵ}_1\text{}{\displaystyle {\int }_{m_e{c}^2}^{Q_{\beta \beta }+2m_e{c}^2-{ϵ}_1}\text{}}\text{}d{ϵ}_2\text{}{\displaystyle {\int }_{m_e{c}^2}^{Q_{\beta \beta }+2m_e{c}^2-{ϵ}_1-{ϵ}_2}\text{}}d{\omega }_1\\ \times f_{11}^{(0)}{w}_{2\nu }\left({\langle K_N\rangle }^2+{\langle L_N\rangle }^2+\langle K_N\rangle \langle L_N\rangle \right)& (43)\end{array}$

$\begin{array}{l}{F}_{2\nu }^{(1)}=\frac{2{\tilde{A}}^2}{9\ln 2}\text{}{\displaystyle {\int }_{m_e{c}^2}^{Q_{\beta \beta }+m_e{c}^2}\text{}}d{ϵ}_1\text{}{\displaystyle {\int }_{m_e{c}^2}^{Q_{\beta \beta }+2m_e{c}^2-{ϵ}_1}\text{}}\text{}d{ϵ}_2\text{}{\displaystyle {\int }_{m_e{c}^2}^{Q_{\beta \beta }+2m_e{c}^2-{ϵ}_1-{ϵ}_2}\text{}}\text{}d{\omega }_1\\ \times f_{11}^{(1)}{w}_{2\nu }\left(2\left({\langle K_N\rangle }^2+{\langle L_N\rangle }^2\right)+5\langle K_N\rangle \langle L_N\rangle \right)& (44)\end{array}$

that correspond to the two parameters a⁽⁰⁾ and a⁽¹⁾ of Equation (32). An energy dependent angular correlation α is defined as

$\begin{array}{ll}\alpha =\frac{d{F}_{2\nu }^{(1)}/d{ϵ}_1}{d{F}_{2\nu }^{(0)}/d{ϵ}_1}& (45)\end{array}$

from

$\begin{array}{ll}\frac{d^2{W}_{2\nu }}{d\cos {\theta }_{12}d{ϵ}_1}=|M_{2\nu }{|}^2\frac{d{F}_{2\nu }^{(0)}}{d{ϵ}_1}(1+\alpha \cos {\theta }_{12})& (46)\end{array}$

In this work, we report normalized phase factors, $G_{2\nu }^{(i)}=F_{2\nu }^{(i)}/[g_A^4{(m_e{c}^2)}^2]$ in units yr⁻¹. We mention that by this renormalization the quantities G⁽ⁱ⁾ do not depend on the axial-vector constant g_A. Also, as one can see from the formulas, the PSF ($G_{(2\nu ,0\nu )}^{(i)}$ depend very little to the average energy < E_N >, because this constant enters both to the nominator and denominator in all the PSF expressions.

4.3.1. Double Electron and Double Positron Decay Modes

Here, the electron phase factors are

$\begin{array}{ll}{f}_{11}^{(0)}=\mid f^{-1-1}{\mid }^2+\mid f_{11}{\mid }^2+\mid f^{-1}{ }_1{\mid }^2+\mid f_1{ }^{-1}{\mid }^2& (47)\end{array}$

and

$\begin{array}{ll}{f}_{11}^{(1)}=-2\Re \left(f^{-1-1}{f}_{11}^{*}+f^{-1}{ }_1{f}_1{ }^{-1*}\right)& (48)\end{array}$

with the products

$\begin{array}{ll}{f}^{-1-1}=g_{-1}({ϵ}_1)g_{-1}({ϵ}_2) ; f_{11}=f_1({ϵ}_1)f_1({ϵ}_2),& (49)\end{array}$

$\begin{array}{ll}{f}_{ 1}^{-1}=g_{-1}({ϵ}_1)f_1({ϵ}_2) ; f_1^{ -1}=f_1({ϵ}_1)g_1({ϵ}_2)& (50)\end{array}$

The values of the f and g functions are approximated with the solutions at the nuclear surface (the method I from [11]).

$\begin{array}{ll}{g}_{-1}(ϵ)=g_{-1}(ϵ,R_A) ; f_1(ϵ)=f_1(ϵ,R_A)& (51)\end{array}$

where $R_A=1.2A^{1/3}$ fm is the nuclear radius. For the 2νββ decay mode and transitions to ground states (g.s.), the PSF expression reads:

$\begin{array}{l}{G}_{2\nu }^{\beta \beta }(0^{+}\to 0^{+})=\frac{2{\tilde{A}}^2}{3\ln 2g_A^4{(m_e{c}^2)}^2}{\displaystyle {\int }_{m_e{c}^2}^{Q^{\beta \beta }+m_e{c}^2}d}{ϵ}_1\\ \times {\displaystyle {\int }_{m_e{c}^2}^{Q^{\beta \beta }+2m_e{c}^2-{ϵ}_1}d}{ϵ}_2{\displaystyle {\int }_0^{Q^{\beta \beta }+2m{c}_e^2-{ϵ}_1-{ϵ}_2}d}{\omega }_1\\ \times f_{11}^{(0)}{w}_{2\nu }\left({\langle K_N\rangle }^2+{\langle L_N\rangle }^2+\langle K_N\rangle \langle L_N\rangle \right)& (52)\end{array}$

where $Q^{\beta \beta }=M(A,Z_0)-M(A,Z_0-2)-4m_e{c}^2$ is the kinetic energy released in the process. 〈K_N〉, 〈L_N〉 are expressions that depend on the electron/positron (ϵ_1,2) and neutrino (ω_1,2) energies, and on the g.s. energy of the parent nucleus and on the excited states energy of the intermediate nucleus [2, 3, 5, 11, 35, 36].

$\begin{array}{ll}\langle K_N\rangle =\frac{1}{{ϵ}_1+{\omega }_1+\langle E_N\rangle -E_I}+\frac{1}{{ϵ}_2+{\omega }_2+\langle E_N\rangle -E_I}& (53)\end{array}$

$\begin{array}{ll}\langle L_N\rangle =\frac{1}{{ϵ}_1+{\omega }_2+\langle E_N\rangle -E_I}+\frac{1}{{ϵ}_2+{\omega }_1+\langle E_N\rangle -E_I}& (54)\end{array}$

Here, the difference in energy in the denominator can be obtained from the approximation ${Ã}^2={\left[W_0/2+\langle E_N\rangle -E_I\right]}^2$, where the empirical systematics à = 1.12A^1/2 (in MeV) [2] gives approximately the energy of the giant Gamow-Teller resonance in the intermediate nucleus [11]. The quantity W₀ is related to the Q^ββ value of the process and

$\begin{array}{ll}{w}_{2\nu }=\frac{g_A^4{(G\cos {\theta }_C)}^4}{64{\pi }^7\hslash }{\omega }_1^2{\omega }_2^2(p_{1}c)(p_{2}c){ϵ}_1{ϵ}_2,& (55)\end{array}$

where ω₁ and ${\omega }_2=Q^{\beta \beta }-{ϵ}_1-{ϵ}_2-{\omega }_1+2m_e{c}^2$ are the neutrino energies. The PSF are finally renormalized to the electron rest energy and are reported in [yr⁻¹].

The PSF for the 2νββ decay mode and transitions to excited $0_1^{+}$ states is calculated with a formula similar to (52), but replacing Q^ββ by $Q(0_1^{+})=Q^{\beta \beta }-E_x(0_1^{+})$, which is the kinetic energy released in this transition. $E_x(0_1^{+})$ is the energy of the excited $0_1^{+}$ state of the daughter nucleus x.

For the 2νββ decay mode and transitions to excited $2_1^{+}$ state, the PSF formula reads [5, 11–13, 39]:

$\begin{array}{l}{G}_{2\nu }^{\beta \beta }(0^{+}\to 2_1^{+})=\frac{2{\tilde{A}}^6}{\ln 2g_A^4{(m_e{c}^2)}^6}{\displaystyle {\int }_{m_e{c}^2}^{Q^{\beta \beta }(2_1^{+})+m_e{c}^2}d}{ϵ}_1\\ \times {\displaystyle {\int }_{m_e{c}^2}^{Q^{\beta \beta }(2_1^{+})+2m_e{c}^2-{ϵ}_1}d}{ϵ}_2{\displaystyle {\int }_0^{Q^{\beta \beta }(2_1^{+})+2m_e{c}^2-{ϵ}_1-{ϵ}_2}d}{\omega }_1\\ \times f_{11}^{(0)}{w}_{2\nu }{\left(\langle K_N\rangle -\langle L_N\rangle \right)}^2& (56)\end{array}$

where $Q^{\beta \beta }(2_1^{+})=Q-E_x(2_1^{+})$.

For the 0νββ decay and transitions to g.s. the PSF reads:

$\begin{array}{ll}{G}_{0\nu }^{\beta \beta }(0^{+}\to 0^{+})=\frac{2}{4g_A^4{R}_A^2\ln 2}{\displaystyle {\int }_{m_e{c}^2}^{Q^{\beta \beta }+m_e{c}^2}{f}_{11}^{(0)}}{w}_{0\nu }d{ϵ}_1& (57)\end{array}$

where

$\begin{array}{ll}{w}_{0\nu }=\frac{g_A^4{(G\cos {\theta }_C)}^4}{16{\pi }^5}{(m_e{c}^2)}^2(\hslash c^2)(p_{1}c)(p_{2}c){ϵ}_1{ϵ}_2& (58)\end{array}$

For transitions on the 2${}_1^{+}$ state, instead of $f_{11}^{(0)}$ in Equation (57), the next expression is integrated

$\begin{array}{ll}{f}_{1+}+f_{1-}=3{(\frac{\hslash c}{m_e{c}^2{R}_A})}^2(|f^{-2-1}{|}^2+|f^{-1-2}{|}^2+|f_{21}{|}^2+|f_{12}{|}^2).& (59)\end{array}$

$\begin{array}{ll}{f}^{-2-1}=g_{-2}({ϵ}_1)g_{-1}({ϵ}_2) ; f_{21}=f_2({ϵ}_1)f_1({ϵ}_2),& (60)\end{array}$

$\begin{array}{ll}{f}^{-1-2}=g_{-1}({ϵ}_1)g_{-2}({ϵ}_2) ; f_{12}=f_1({ϵ}_1)f_2({ϵ}_2)& (61)\end{array}$

In Equations (52, 56, 57), it is convenient to redefine the PSF by a renormalization that eliminates the constant g_A and correlates (by dividing by $4R_A^2$) the dimension of G_0ν with the NME which are dimensionless. These are reflected in the lifetimes formulas (9). Thus, the PSF are also reported in [yr⁻¹]. A similar expression is employed in the PSF calculation for the transitions to excited $0_1^{+}$ states, but replacing Q^ββ by $Q^{\beta \beta }(0_1^{+})$. The formula used for the PSF computation for 2νβ⁺β⁺ decay mode is similar to that used for 2νβ⁻β⁻ decay, but ϵ_1,2 are now the positron energies. Also, we use the same approximations as described above, to evaluate the radial positron wave functions (g and f) at the nuclear surface and replace the excitation energy E_N in the intermediate odd-odd nucleus by a suitable average energy.

4.3.2. The ECβ⁺ Case

For the ECβ⁺ decays the energy released in the process is $Q^{EC\beta }=M(A,Z_0)-M(A,Z_0-2)-2m{c}_e^2$. If the numerical solutions of the Dirac equation are obtained in Bohr units a₀, the probability that an electron is found on the surface of a nucleus of radius R_A can be defined as:

$\begin{array}{ll}{B}_{n,\kappa }^2=\frac{1}{4\pi {(m_e{c}^2)}^3}{\left(\frac{\hslash c}{a_0}\right)}^3{\left(\frac{a_0}{R_A}\right)}^2[g_{n,\kappa }^2(R_A)+f_{n,\kappa }^2(R_A)]& (62)\end{array}$

The PSF expression for 2νECβ decay mode is

$\begin{array}{l}{G}_{2\nu }^{EC{\beta }^{+}}=\frac{2A^2}{3\ln 2}\frac{{(G\cos \theta )}^4}{16{\pi }^5\hslash }(m_e{c}^2){\displaystyle \sum _{i=0,1}{B}_{i,-1}^2}\\ \times {\displaystyle {\int }_{m_e{c}^2}^{Q^{EC\beta }+{ϵ}_{i,-1}+m_e{c}^2}{\displaystyle {\int }_0^{Q^{EC\beta }+{ϵ}_{i,-1}-{ϵ}_p+m_e{c}^2}}}\\ \times [g_{-1}^2({ϵ}_p)+f_1^2({ϵ}_p)]\\ \times \left({\langle K_N\rangle }^2+{\langle L_N\rangle }^2+\langle K_N\rangle \langle L_N\rangle \right){\omega }_1^2{\omega }_2^2{p}_{p}c{ϵ}_{p}d{\omega }_{1}d{ϵ}_p& (63)\end{array}$

where ϵ_{n, κ} are the binding energies of the electron while p_p and ϵ_p are the linear momentum and the energy of the positron. Here, the expressions for 〈K_N〉 and 〈L_N〉 are similar to those from Equations (53)–(54), but where ϵ₁ is replaced by ${ϵ}_{i,-1}^c=m_e{c}^2-{ϵ}_{i,-1}$, the energy of the captured electron and ϵ₂ is replaced by ϵ_p, the energy of the emitted positron. For the 0νECβ decay process the PSF expression is:

$\begin{array}{l}{G}_{0\nu }^{EC{\beta }^{+}}=\frac{1}{4R_A^2}\frac{2}{\ln 2}\frac{{(G\cos \theta )}^4}{4{\pi }^3}\left(\hslash c^2\right){\left(m_e{c}^2\right)}^5\\ \times {\displaystyle \sum _{i=0,1}{B}_{i,-1}^2}[g_{-1}^2({ϵ}_{p,i})+f_1^2({ϵ}_{p,i})]p_{p}c{ϵ}_p& (64)\end{array}$

where ϵ_{p, i} denotes the maximal value of the positron associated to the state i.

4.3.3. The 2νECEC Case

The PSF expression is defined as:

$\begin{array}{l}{G}_{2\nu }^{ECEC}=\frac{2{\tilde{A}}^2}{3\ln 2}\frac{{(G\cos \theta )}^4}{16{\pi }^3\hslash }{\left(m_e{c}^2\right)}^4\times {\displaystyle \sum _{i,j=0,1}{B}_{i,-1}^2}{B}_{j,-1}^2\\ \times {\displaystyle {\int }_0^{Q^{ECEC}+{ϵ}_{i,-1}+{ϵ}_{j,-1}}{(\langle K_N\rangle }^2}+{\langle L_N\rangle }^2\\ +\langle K_N\rangle \langle L_N\rangle ){\omega }_1^2{\omega }_2^{2}d{\omega }_1& (65)\end{array}$

where $Q^{ECEC}=M(A,Z_0)-M(A,Z_0-2)$ is the energy released in the process. The expressions for 〈K_N〉 and 〈L_N〉 are similar to those from Equations (53)–(54), but where ϵ_{(1, 2)} are replaced by ${ϵ}_{(i,j)-1}^c=m_e{c}^2-{ϵ}_{(i,j)-1}$, the energies of the captured electrons.

5. Numerical Details

The numerical solutions of the Dirac equation are based on the Buhring method [45] which considers a power series expansion of the radial wave functions. Essentially, the numerical algorithm of our code follows the lines given in Salvat and Mayol [42] and Salvat et al. [38]. In such a calculation, the uncertainties are introduced by the interpolation procedure which determine potential values and by the round-off errors. A very dense grid of potential values was taken into account in order to increase the accuracy. The finite size and the screening effects modify the phase shifts. These phase shifts are obtained through a normalization of the free wave functions as explained in Salvat et al. [38, 41].

For bound and free states wave functions, the potential energy corrected with the screening function is calculated in different mesh points defined by an increment x along the distance r. Further, the potential is approximated with a spline cubic function. As mentioned, the radial wave function is expanded as an infinite power series. Successive numerical values of the radial wave function depend on the increment and the coefficients of the spline function by mean of recurrence relations. At r = 0, the solution is determined from the regularity condition. The wave functions are determined step by step in the mesh points. If the round-off errors are considered negligible, the numerical precision depends mainly on the number of terms in the series expansion and on the value of the increment. An increment interval of 10⁻⁴ fm and at least 100 terms in the series expansion were considered in our numerical code, exceeding the convergence criteria given Salvat and Mayol [42]. At very large distances, the behavior of the free wave functions determined numerically must resemble to that of the asymptotic Coulomb function, characterized by a value of the Sommerfeld parameter modified by the screening effect. By comparing the behavior of the numerical and the analytical asymptotic solutions, a renormalization to unity and a determination of the phase shift is obtained. For discrete states, the potential provides the boundaries of the solutions. The initial conditions of the reverse solution of the wave function is determined by the WKB approximation. The direct and the reverse wave functions are calculated numerically for different electron energies. The eigenvalue is obtained when the direct solutions and the inverse ones match together. In order to find the bound states of the electron we spanned a range 0.3 MeV smaller than $m_e{c}^2$ in steps of 0.0002 MeV. The bound solutions of interest are found in this interval.

Concerning the PSF evaluations, it should be noted that all integrals were performed accurately with Gauss-Legendre quadrature in 32 points. Up to 49 values of the radial wave functions on the surface of the nucleus were determined in the Q^ββ-value interval. The electron phase spaces are obtained by spline interpolations of the values calculated for the electron wave functions. Recently reported Q^ββ in Wang et al. [32] are used.

6. Results and Discussions

6.1. Phase Space Factors

The PSF values reported in the literature within different models are presented in Tables 1–7. As we mentioned, the PSF values reported by KI in Kotila and Iachello [11], Kotila and Iachello [12], Stoica and Mirea [13] are obtained with a similar approach as ours, namely, the use of exact electron wave functions obtained by solving numerically the Dirac equation and with inclusion of the finite nuclear size and electron screening effects. The results reported in Doi et al. [37], Doi et al. [3], Doi and Kotani [35], and Suhonen and Civitarese [5] are obtained by approximating the electron wave functions at the nuclear surface and without inclusion of the screening effect.

TABLE 1

Table 1. PSF for 2νβ⁻β⁻ decays to final g.s.

TABLE 2

Table 2. PSF for 0νβ⁻β⁻ decays to final g.s.

TABLE 3

Table 3. PSF for β⁻β⁻ decays to final excited 0${}_1^{+}$ states.

TABLE 4

Table 4. PSF for β⁻β⁻ decays to final excited 2${}_1^{+}$ states.

TABLE 5

Table 5. PSF for β⁺β⁺ decay mode.

TABLE 6

Table 6. PSF for ECβ⁺ decay mode.

TABLE 7

Table 7. PSF for 2νECEC decay mode.

The PSF values for the β⁺β⁺ decays are displayed in Table 5. Differences of about 11–30% between G_2ν values of Neacsu and Horoi [11] and Mirea et al. [39] are retrieved, in majority of cases. The agreement concerning the G_0ν values is very good with one exception: in the case of the ¹³⁶Ce nucleus, a difference of 13% is found. The differences betweenthe recent PSF values calculated exactly and those reported in Doi and Kotani [35] are larger for most of the results.

The results for ECβ⁺ decay mode are shown in Table 6 and they are grouped in two sectors: the first sector includes the six nuclei that can also undergo a β⁺β⁺ decay and the second one includes the other nuclei which can undergo only ECβ⁺ decays. One can see that the differences between KI results and ours are relevant in many cases. For the nuclei from the first sector (the most interesting experimentally cases, having the largest Q^ECβ values), we got several PSF values which differ by more than 20% compared to those reported by KI and other authors. For the “pure” ECβ⁺ decays, we got even more significant differences in several cases as compared with KI results, while other (older) results reported are a few.

The PSF values (G_2ν, G_0ν) for β⁻β⁻ decays are presented in Tables 1–4 for transitions to g.s and excited $0_1^{+}$ and $2_1^{+}$ states, respectively. For the transitions to g.s. (Tables 1, 2) the results reported in Kotila and Iachello [11, 12] are in very good agreement. Apart two exceptions, the differences between the two set of results do not exceed 5%. In the case of nucleus ¹²⁸Te the difference betwee the G_2ν values amounts to 22%. In the case of the ²³⁸U nucleus the difference is of about 7 times. We double checked our calculations for this nucleus and did not found errors. For the transitions to excited $0_1^{+}$ states (Table 3), there are several cases, especially for heavier nuclei (⁷⁶Ge, ¹¹⁰Pd, ¹¹⁶Cd, ¹³⁰Te, ¹³⁶Xe) with notable differences between the two sets of results, ranging in the interval (10–40)%. The difference for the G_2ν in the case of ²³⁸U is the largest one. The comparison between the recent results and those previously reported in Doi et al. [37], Doi et al.[3], and Suhonen and Civitarese [5] reveal also larger discrepancies in several cases. For the transitions to excited $2_1^{+}$ states (Table 4) recent results are reported in Neacsu and Horoi [46], by using an effective screening method, and with approximate wave functions. Comparing our PSF values for $2_1^{+}$ final states with those reported in Doi et al. [37] and Suhonen and Civitarese [5], one observes that the differences range in a 30% interval, both for G_0ν and G_2ν. A much better agreement is obtained between the data of Neacsu and Horoi [46] and Mirea et al. [39]. It is worth to mention that, although they are not yet measured, the β⁻β⁻ decays to excited $2_1^{+}$ states are of interest to probe alternative mechanisms for occurrence of the 0ν-decay mode, to the (most common) mechanism of exchange of light left-handed neutrinos between two nucleons inside the nucleus.

For the 2νECEC decay mode, the PSF values are given in Table 7. Three groups of nuclei are considered. The first sector includes six nuclei that can also undergo a β⁺β⁺ decay. In the second sector, nuclei that can also undergo ECβ⁺ decays are enumerated. The third sector includes the nuclei which can undergo only 2νECEC decays. The differences between the G_2ν values of the recent calculations reported in Kotila and Iachello [11] and Mirea et al. [39] concerning the nuclei from the first sector are within (30-55)%. In the remaining two sectors, the discrepancies are even larger. These differences between the PSF values that are obtained with similar approaches on the one side, and with other approximations on the other side, might be explained by the use of less rigorous approaches, i.e., of approximate electron wave functions, non-inclusion of the screening effect, (possible) less efficient routines at that time, etc.

The differences between results obtained with exact solutions of the Dirac equation [11, 39] may arise from the additional ingredients used in the latter calculation. For example, we got relevant differences between their results and ours especially for decays where the Q^ββ values are small. For small interval for energy integration, the numerical precision of the PSF expressions depend on values of the free electron wave functions for energies close to $m_e{c}^2$ where the solutions of the Dirac equation diverges. These situations should be treated carefully. A large number of mesh points in the numerical integration over Q^ββ values is needed, especially at energies close to the electron mass, in order to get stable results. The calculations performed in Mirea et al. [39] are realized with up to 49 interpolation points. The authors of Mirea et al. [39] observed that after 40 interpolation points the results show no variation. Around the zero kinetic energy of the released electron/positron, an interval energy of about 10⁻⁴ MeV was used, because in this region a sudden variation of the wave function amplitude is produced. The use of a Coulomb potential constructed from the density of protons can also contribute to differences of a few percents, especially for EC decays. The use of different reported Q^ββ values may also contribute to the observed differences up to 3% as remarked in Mirea et al. [39], were the tables of Wang et al. [32] were used. In the cases involving bound electron states, important is the numerical accuracy in the identification of electron/positron bound state 1s_1/2 that intervene in the Equation (62). A different procedure was built in Mirea et al. [39] to get the correct electron/positron bound states. The mentioned possible sources of uncertainties can act cumulatively to determine differences between the results obtained by solving numerically the Dirac equation. In our opinion, differences larger that 10% between the values of different PSF computations are important for precise DBD calculations/estimations, hence the evaluation of the PSF for DBD is still a challenge.

6.2. Energy and Angular Distributions

From the DBD decay rates formulas for the energy and angular distributions of the two electrons/positrons released in the decay can also be derived. The normalized energy distributions $d{G}_{2\nu }^{(0)}/d{ϵ}_1$, $d{G}_{2\nu }^{(0)}/d({ϵ}_1+{ϵ}_2)$ and the angular correlation parameter α are calculated for a series of nuclei of interest and are plotted with a thick full line in Figures 4–6 for β⁻β⁻ decay. The range of energies span the Q^ββ-value for each investigated nucleus. We also performed calculations without considering the screening effect. The distributions without screening were plotted with a thin full line. Moreover, a calculation of the angular correlation parameter α is realized by considering that the global phase shift is zero and it is plotted with a dashed thin line. The effect of the screening becomes more pronounced for heavier nuclei, as for ¹³⁶Xe nucleus. Due to the screening, the values of the differential distributions decrease with few percents. The screening is also responsible for a small attenuation of the angular distribution α at energies close to the Q^ββ-value. A similar behavior was also remarked in Kotila and Iachello [11]. If the global phase shift is set to zero, the absolute values of the angular correlation α become larger. Our results can be compared with single-electron, summed energy spectra and angular correlation spectra provided in Kotila and Iachello [11] for the parent nuclei ¹³⁶Xe, ⁸²Se, and ¹⁵⁰Nd. The behavior of the differential decay rates exhibited by the results given Kotila and Iachello [11] agree with our data, but differ in magnitude. In general, the differential decay rates reported here are larger than those obtained in Kotila and Iachello [11] with about 10%. Also, the α parameter calculated by disregarding the global phase shift in this work seems to agree better with the results published in Kotila and Iachello [11]. New calculations concerning the energy distributions were recently reported in Simkovic et al. [48], where modifications of the shape of the energy distributions were evidenced.

FIGURE 4

Figure 4. The normalized energy distributions $d{G}_{2\nu }^{(0)}/d{ϵ}_1$, $d{G}_{2\nu }^{(0)}/d({ϵ}_1+{ϵ}_2)$ and the angular correlation parameter α as function of the electron energies ϵ₁ and ϵ₂ for ⁸²Se parent nucleus β⁻β⁻ decay are plotted with a full thick line. The results obtained by neglecting the screening effect are plotted with a full thin line. A thin dashed line is used for calculations in which the global phase shift of the Dirac wave functions is set to zero.

FIGURE 5

Figure 5. Same as Figure 4 for the ¹³⁶Xe parent nucleus.

FIGURE 6

Figure 6. Same as Figure 4 for the ¹⁵⁰Nd parent nucleus.

7. Conclusions

A review of calculated PSF involved in the β⁻β⁻, β⁺β⁺, ECβ⁺ and 2νECEC DBD modes is presented. Different models, from non-relativistic approaches with approximate wave functions calculations to relativistic approaches with exact electron wave functions obtained by solving numerically the Dirac equation, and with inclusion of finite nuclear size and electron screening effects are presented. For each decay mode, we present tables of PSF values for the majority of the isotopes decaying double-beta. We also compare our results with other similar ones from literature and the differences are analyzed and discussed. Finally, using our method, we also calculate the energy and angular distributions for three isotopes, namely ¹³⁶Xe, ⁸²Se and ¹⁵⁰Nd. Accurate values of the PSF are required in the DBD investigations, in order to improve the lifetimes predictions, or the constraint the neutrino parameters [49], or for experimentalists [50], to plan their set-ups.

Author Contributions

SS and MM elaborated a new code to calculate the values of phase spaces factors involved in double-beta decay, starting with the resolution of the Dirac equation. The values reported as Stoica and Mirea [13] and Mirea et al. [39] are calculated by the authors.

Funding

The authors would like to acknowledge the support of a grant of the Romanian Ministry of Research and Innovation, CNCS-UEFISCDI, project number PN-III-P4-ID-PCE-2016-0078.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

References

1. Primakoff H, Rosen SP. Double beta decay. Rep Prog Phys. (1959) 22:121–66.

Google Scholar

2. Haxton WC, Stephenson GJ. Double beta decay. Prog Part Nuclear Phys. (1984) 12:409–79.

Google Scholar

3. Doi M, Kotani T, Takasugi E. Double beta Decay and Majorana neutrino. Prog Theor Phys Suppl. (1985) 83:1–175.

Google Scholar

4. Tomoda T. Double beta decay. Rep Prog Phys. (1991) 54:53–126.

Google Scholar

5. Suhonen J, Civitarese O. Weak-interaction and nuclear-structure aspects of double beta decay. Phys Rep. (1998) 300:123–214.

Google Scholar

6. Avignone FT, Elliott SR, Engel J. Double beta decay, Majorana neutrinos, and neutrino mass. Rev Modern Phys. (2008) 80:481–516. doi: 10.1103/RevModPhys.80.481

CrossRef Full Text | Google Scholar

7. Ejiri H. Double beta decays and neutrino nuclear responses. Prog Part Nuclear Phys. (2010) 64:249–57. doi: 10.1016/j.ppnp.2009.12.021

CrossRef Full Text | Google Scholar

8. Rodejohann W. Neutrino-less double beta decays and particle physics. Int J Modern Phys. (2011) 20:1833–930. doi: 10.1142/S0218301311020186

CrossRef Full Text | Google Scholar

9. Vogel P. Nuclear structure and double beta decay. J Phys G. (2012) 39:124002. doi: 10.1088/0954-3899/39/12/124002

CrossRef Full Text | Google Scholar

10. Vergados J, Ejiri H, Simkovic F. Theory of neutrinoless double-beta decay. Rep Prog Phys. (2012) 75:106301. doi: 10.1088/0034-4885/75/10/106301

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Kotila J, Iachello F. Phase-space factors for double-β decay. Phys Rev. (2012) 85:034316.

Google Scholar

12. Kotila J, Iachello F. Phase space factors for β⁺β+ decay and competing modes of double-β decay. Phys Rev. (2013) 87:024313. doi: 10.1103/PhysRevC.87.024313

CrossRef Full Text | Google Scholar

13. Stoica S, Mirea M. New calculations for phase space factors involved in double-β decay. Phys Rev. (2013) 88:037303. doi: 10.1103/PhysRevC.88.037303

CrossRef Full Text | Google Scholar

14. V A Rodin AF, Simkovic F, Vogel P. Uncertainty in the 0νββ decay nuclear matrix elements. Phys Rev. (2003) 68:044302. doi: 10.1103/PhysRevC.68.044302

CrossRef Full Text | Google Scholar

15. Simkovic F, Faessler A, Rodin VA, Vogel P, Engel J. Anatomy of the 0νββ nuclear matrix elements. Phys Rev. (2008) 77:045503. doi: 10.1103/PhysRevC.77.045503

CrossRef Full Text

16. Simkovic F, Faessler A, Muther H, Rodin VA, Stauf M. 0νββ-decay nuclear matrix elements with self-consistent short-range correlations. Phys Rev. (2009) 79:055501. doi: 10.1103/PhysRevC.79.055501

CrossRef Full Text

17. Fang DL, Faessler A, Rodin V, Simkovic F. Neutrinoless double-β decay of deformed nuclei within quasiparticle randem-phase approximation with a realistic interaction. Phys Rev. (2011) 83:034320. doi: 10.1103/PhysRevC.83.034320

CrossRef Full Text | Google Scholar

18. Kortelainen M, Civitarese O, Suhonen J, Toivanen J. Short-range correlations and neutrinoless double beta decay. Phys Lett. (2007) 647:128–32. doi: 10.1016/j.physletb.2007.01.054

CrossRef Full Text | Google Scholar

19. Kortelainen M, Suhonen J. Improved short-range correlations and 0νββ nuclear matrix elements of ⁷⁶Ge and ⁸²Se. Phys Rev. (2007) 75:051303(R). doi: 10.1103/PhysRevC.75.051303

CrossRef Full Text | Google Scholar

20. Stoica S, Klapdor-Kleingrothaus HV. Critical view on double-beta decay matrix elements with Quasi Random Phase Approximation-based methods. Nucl Phys A (2001) 694:269–94. doi: 10.1016/S0375-9474(01)00988-5

CrossRef Full Text | Google Scholar

21. Caurier E, Zuker AP, Poves A, Martinez-Pinedo G. Full pf shell model study of A=48 nuclei. Phys Rev. (1994) 50:225–36.

PubMed Abstract | Google Scholar

22. Retamosa J, Caurier E, Nowacki F. Neutrinoless double beta decay of ⁴⁸Ca. Phys Rev. (1995) 51:371–8.

PubMed Abstract | Google Scholar

23. Caurier E, Menendez J, Nowacki F, Poves A. Influence of pairing on the nuclear matrix elements of the neutrinoless ββ decays. Phys Rev Lett. (2008) 100:052503. doi: 10.1103/PhysRevLett.100.052503

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Menendeza J, Poves A, Caurier E, Nowacki F. Disassembling the nuclear matrix elements of the neutrinoless ββ decay. Nuclear Phys. (2009) 818:139–51. doi: 10.1016/j.nuclphysa.2008.12.005

CrossRef Full Text | Google Scholar

25. Horoi M, Stoica S, Brown BA. Shell-model calculations of two-neutrino double-β decay rates of ⁴⁸Ca with the GXPF1A interaction. Phys Rev. (2007) 75:034303. doi: 10.1103/PhysRevC.75.034303

CrossRef Full Text | Google Scholar

26. Horoi M, Stoica S. Shell-model analysis of the neutrinoless double-β decay of ⁴⁸Ca. Phys Rev. (2010) 81:024321. doi: 10.1103/PHYSREVC.81.024321

CrossRef Full Text | Google Scholar

27. Barea J, Iachello F. Neutrinoless double-β decay in the microscopic interacting boson model. Phys Rev. (2009) 79:044301. doi: 10.1103/PhysRevC.79.044301

CrossRef Full Text | Google Scholar

28. Barea J, Kotila K, Iachello F. Limits on neutrino masses from neutrinoless double-β decay. Phys Rev Lett. (2012) 109:042501. doi: 10.1103/PhysRevLett.109.042501

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Rodriguez TR, Martinez-Pinedo G. Energy density functional study of nuclear matrix elements for neutrinoless ββ decay. Phys Rev Lett. (2010) 105:252503. doi: 10.1103/PhysRevLett.105.252503

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Rath PK, Chandra R, Chaturvedi K, Raina PK, Hirsch JG. Uncertainties in nuclear transition matrix elements for neutrinoless ββ decay within the projected-Hartree-Fock-Bogoliubov model. Phys Rev. (2010) 82:064310. doi: 10.1103/PhysRevC.82.064310

CrossRef Full Text | Google Scholar

31. Neacsu A, Stoica S. Shell model calculations for neutrinoless double beta decay through the exchange of heavy neutrinos. Roman Rep Phys. (2014) 66:376–81.

Google Scholar

32. Wang M, Audi G, Wapstra AH, Kondev FG, MacCornick M, Xu X, et al. The AME2012 atomic mass evaluation. Chinese Phys C. (2012) 36:1603–2014. doi: 10.1088/1674-1137/36/12/003

CrossRef Full Text | Google Scholar

33. Haxton WC, Stephenson GJ, Strottman D. Lepton-number conservation and the double-β decay of ¹³⁰Te. Phys Rev. (1982) 25:2360–9.

Google Scholar

34. Behrens H, Janecke J. Numerical Tables for Beta Decay and Electron Capture. Berlin: Springer (1969).

35. Doi M, Kotani T. Neutrino emitting of double beta decay. Prog Theor Phys. (1992) 87:1207–31.

Google Scholar

36. Doi M, Kotani T. Neutrinoless modes of double beta decay. Prog Theor Phys. (1993) 89:139–59.

Google Scholar

37. Doi M, Kotani T, Nishiura H, Takasugi E. Double beta decay. Prog Theor Phys. (1983) 69:602–35.

Google Scholar

38. Salvat F, Fernandez-Varea JM, Williamson W. Accurate numerical solution of the radial Schrodinger and Dirac wave equations. Comp Phys Commun. (1995) 90:151–68.

Google Scholar

39. Mirea M, Pahomi T, Stoica S. Values of the phase space factors involved in double beta decay. Roman Rep Phys. (2015) 67:872–89. doi: 10.1063/1.4934913

CrossRef Full Text | Google Scholar

40. Doi M, Kotani T, Nishiura H, Takasugi E. The energy spectrum and the angular correlation in the ββ decay. Prog Theor Phys. (1983) 70:1353–74.

Google Scholar

41. Salvat F, Jablonski A, Powell CJ. ELSEPA-Dirac partial-wave calculation of elastic scattering of electrons and positrons by atoms, positive ions and molecules. Comp Phys Commun. (2005) 165:157–90. doi: 10.1016/j.cpc.2004.09.006

CrossRef Full Text | Google Scholar

42. Salvat F, Mayol R. Accurate numerical solution of the Schrodinger and Dirac wave equations for central fields. Comp Phys Commun. (1991) 62:65–79. doi: 10.1016/0010-4655(91)90122-2

CrossRef Full Text | Google Scholar

43. Esposito S. Majorana solution of the Thomas-Fermi equation. Am J Phys. (2002) 70:852–6. doi: 10.1119/1.1484144

CrossRef Full Text | Google Scholar

44. Doi M, Kotani T, Nisiura H, Okuda K, Takasuki E. Neutrino mass, the right-handed interaction and the double beta decay. I: formalism. Prog Theor Phys. (1981) 66:1739–64.

Google Scholar

45. Buhring W. Computational improvements in phase shift calculations of elastic electron scattering. Zeitschrift fur Physik (1965) 187:180–96.

Google Scholar

46. Neacsu A, Horoi M. An effective method to accurately calculate the phase space factors for β⁻β⁻ decay. Adv High Energy Phys. (2016) 2016:7486712. doi: 10.1155/2016/7486712

CrossRef Full Text | Google Scholar

47. Stefanik D, Dvornicky R, Simkovic F, Vogel P. The light neutrino exchange mechanism of the 0νββ-decay with left-and right-handed leptonic and currents revisited. Phys Rev. (2015) 92:055502. doi: 10.1103/PhysRevC.92.055502

CrossRef Full Text | Google Scholar

48. Simkovic F, Dvornicky R, Stefanik D, Faessler A. Improved description of the 2νββ-decay and a possibility to determine the effective axia-vector coupling constant. Phys Rev. (2018) 97:034315. doi: 10.1103/PhysRevC.97.034315

CrossRef Full Text | Google Scholar

49. Barabash AS. Average and recommended half-life values for two-neutrino beta decay. Nucl Phys. (2015) 935:52–64. doi: 10.1016/j.nuclphysa.2015.01.001

CrossRef Full Text | Google Scholar

50. Arnold R, Augier C, Baker JD, Barabash AS, Basharina-Freshville A, Blondel S, et al. Search for neutrinoless double-beta decay of ¹⁰⁰Mo with the NEMO-3 detector. Phys Rev D (2014) 89:111101. doi: 10.1103/PhysRevD.89.111101

CrossRef Full Text | Google Scholar