Phase Space Factors for Double-Beta Decays

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νβ⁻ β⁻)

1b) Two neutrino double-positron decay

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

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

2) Decays occurring with LNV by two units

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

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

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

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

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]:

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:

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:

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].

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]:

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]:

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:

for β⁻ decay and

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:

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:

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

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

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

where

with

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:

where the charge density is

Ψ_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.

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:

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:

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:

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:

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:

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:

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.

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:

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

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

Also,

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

are changed within averaged values

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

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]:

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

from

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

and

with the products

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

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:

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].

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

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]:

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

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

where

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

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:

The PSF expression for 2νECβ decay mode is

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:

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:

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.

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.

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