New Phase Space Calculations for β-Decay Half-Lives

We revisit the computation of the phase space factors (PSF) involved in the positron decay and EC processes for a large number of nuclei of experimental interest. To obtain the electron/positron wave functions, we develop a code for solving accurately the Dirac equation with a nuclear potential derived from a realistic proton density distribution in the nucleus. The finite nuclear size (FNS) and screening effects are included through recipes which differ from those used in previous calculations. Comparing our results with previous calculations, performed with the sameQ-values, we find a close agreement for positron decays, while, for the EC process, there are relevant differences. For the EC process, we also find that the screening effect has a notable influence on the computed PSF values especially for light nuclei. Further, we recomputed the same PSF values but using the most recent Q-values reported in literature. In several cases, the new Q-values differ significantly from the older ones, leading to large differences in the PSF values as compared with previous results. Our new PSF values can contribute to more reliable calculations of the beta-decay rates, used in the study of nuclei far from the stability line and stellar evolution.


Introduction
The phase space factors for beta decay and electron capture were calculated since a long time [1][2][3] and were considered to be evaluated with sufficient accuracy. However, in those works, the distortion of the electron wave functions (w.f.) by the Coulomb field of the nucleus was taken into account through Fermi functions which were expressed in terms of approximate radial solutions of the Dirac equation at the nuclear surface. Also, other corrections were introduced in the calculations in approximate ways. Thus, the screening effect on spectrum was included by various recipes, for example, by replacing ( ) potential with a momentum dependent screening (for low energy positrons) [3] and by modifying the electron radial w.f. [4,5]. Also, the finite size of the nucleus (FNS) was taken into account by adding to the Fermi functions obtained in the "point-nucleus" approximation corrections that depend on -particle energy and nuclear charge [6,7]. Also, for the nuclear radius, older formula has been used [3,8]. For the EC process, the electron bound-state radial w.f. were also obtained as approximate solution of the Dirac equation evaluated at the nuclear surface. They were improved by including exchange and overlap corrections, which were obtained within a relativistic HF approach.
In this work, we revisit the computation of the PSF involved in the positron decay and electron capture (EC) processes for light and heavy nuclei of experimental interest. The Dirac equation is solved numerically with a Coulomb potential derived from a realistic proton distribution in the nucleus which includes the FNS correction. The numerical procedure follows the power series method described in [9] and is similar to that described in [10,11]. The screening effect was introduced by using a screened Coulomb potential, obtained by multiplying the Coulomb potential by function ( ), solution of the Thomas-Fermi equation obtained by the Majorana method [12]. The accuracy imposed in our numerical algorithms used to solve the Dirac equation always exceeds the convergence criteria given in those references. 2 Advances in High Energy Physics Also, a more efficient procedure to identify the electron bound states without ambiguity was developed.
In order to make a comparison between the actual PSF values found in literature and ours, the same PSF are also computed with the approach described in [2,3] and using the same -values. For positron decays, our results are in close agreement with the other previous results, while, for the EC process, we found significant differences. For these processes, we also find that the screening effect has a notable influence on the computed PSF values for light nuclei. Further, we recomputed the same PSF values using updated -values, reported recently in literature [13,14], which, for several light nuclei, differ significantly from the older ones. As an example we cite the maximum -particle energy (referred to as 0 throughout this paper) stated in Table 2 of [15]. These 0 -values differ considerably from those given in [13,14]. One reason for this big difference could be that Wilkinson and Macefield, in order to compare their calculation with those performed earlier by Towner and Hardy [16], restricted their phase space to only pure Fermi transitions. In other words, the Gamow-Teller window was not accessed in phase space calculation of [15]. Thus, in this paper, we propose new PSF values computed with a more accurate method and using updated -values, for a large number of nuclei of experimental interest. Our calculations can be useful for more reliable computation of the beta-decay rates of nuclei far from the stability line, as well as for better understanding of the stellar evolution.
Our work is further motivated by similar calculations done for the double beta-decay (DBD) process. The PSF for DBD were also considered for a long time to be computed with enough accuracy and were used as such for predicting DBD lifetimes. However, recently, they were recalculated with improved methods, especially for positron and EC decay modes [17,18], and several differences were found as compared to previous calculations where approximate electron/positron w.f. were used.
The paper is organized as follows. In Section 2, we present briefly the two approaches used to compute our PSF values. Our results are reported in Section 3. Here, we compare them with experimental data and previous results and discuss the differences. Finally, we summarize the main points and present our conclusions in Section 4.

Formalism
Following essentially the formalism from [3], we give here the necessary equations which we use to calculate the PSF.

Phase Space
Factors for + Transitions. The probability per unit time that a nucleus with atomic mass and charge decays for an allowed -branch is given by where is the weak interaction coupling constant, is the momentum of -particle, = √ 2 + 1 is the total energy of -particle, and 0 is the maximum -particle energy. 0 = − 1, in + decay ( is the mass difference between initial and final states of neutral atoms). Equation (1) is written in natural units (ℏ = = = 1) so that the unit of momentum is , the unit of energy is 2 , and the unit of time is ℏ/ 2 . Shape factors 0 ( , ) for allowed transitions which appear in (1) are defined as where 0,1 are the nuclear matrix elements and the Fermi functions 1 ( , ). Thus, for calculating the + decay rates, one needs to calculate the nuclear matrix elements and the PSF that can be defined as For the allowed decays, the Fermi functions are expressed as where −1 ( , ) and 1 ( , ) are the large and the small radial components of the positron radial wave functions evaluated at nuclear radius which can be obtained by solving the Dirac equation: where is the central potential for the positron and = ( − )(2 + 1) is the relativistic quantum number. We note that (5) is also written in natural units.
An important step in the PSF calculation for + decay is the method of obtaining the positron continuum radial functions. For this, we develop a new method (code) of solving the Dirac equation, which is adapted from the method used previously for the computation of PSF for DBD process [18,20].
We solved (5) in a nuclear potential ( ) derived from a realistic proton density distribution in the nucleus. This is done by solving the Schrodinger equation with a Woods-Saxon potential. In this case, where the charge density is where Ψ is the proton (Woods-Saxon) w.f. of the spherical single particle state and V is its occupation amplitude. Factor (2 + 1) reflects the spin degeneracy. The screening effect is taken into account by multiplying the expression of ( ) with function ( ), which is the solution of the Thomas-Fermi equation: 2 / 2 = 3/2 /√ , with = / , ≈ 0.8853 0 −1/3 and 0 = Bohr radius. It is calculated within the Majorana method [12]. The boundary conditions are (0) = 1 and (∞) = 0. As mentioned above, the screening effect is taken into account by a method developed in [12]. The possible ways in which the screening function modifies the Coulomb potential depend on the specific mechanism and its boundary conditions.
For the case of + -decay process, the potential used to obtain the electron w.f. is to take into account the fact that decay releases a final negative ion with charge −1. ( , ) is positive. In our approach, we considered the solution of the Thomas-Fermi equation as a universal function, giving an effective screening.
Here, product ℏ = 1, for atomic units. The asymptotic potential between an positron and an ionized atom is + = −1. In this case, the charge number = 0 − 1 corresponds to the daughter nucleus, 0 being the charge number of the parent nucleus. Asymptotically, ( ) tends to zero.
In this case, the radial solutions of the Dirac equations should be normalized in order to have the following asymptotic behavior: where is the speed of the light, / are the electron mass/energy, = /ℏ is the electron wave number, is the phase shift, and is the Coulomb interaction energy between the electron and the daughter nucleus.
On the other side, we also calculated the PSF for positron decays with the method described in [3]. and functions were calculated by solving the Dirac equation for a point-nucleus unscreened Coulomb potential, for which the equation has analytical solutions. The finite nuclear size and screening effects were introduced as corrections, after the recipe described in [3]. The finite-size correction was introduced by means of an empirical deviation that depended 4 Advances in High Energy Physics on atomic mass and energy [6,7]. The screening correction was given by the following replacement [4,5]: where = + 0 , = √( ) 2 − 1 and 0 was taken as a -dependent screening potential. For further details of this formalism, we refer to [3]. No electromagnetic corrections were undertaken in this calculation of PSF.

Phase Space Factors for Electron Capture (EC).
Electron capture is always an alternate decay mode for radioactive isotopes that do not have sufficient energy to decay by positron emission. This is a process which competes with positron decay. In order for electron capture leading to a vacancy in, say, K-shell, the atomic mass difference between initial and final states, , must be greater than the binding energy of K-shell electron in the daughter atom, . The energy carried off by the neutrino is then given by If the energy requirement > is satisfied, electron capture from K-shell is more probable than that from any other shell because of the greater density at the nucleus of Kshell electrons. The total K-shell capture rate can be expressed as where where 2 is a constant (with dimensions of −1 ), 0,1 are specific combinations of nuclear matrix elements, is the large component of the bound-state radial w.f. of the captured K-shell electron (evaluated at the nuclear surface ), is the neutrino energy in units of 2 , and is the "exchange" correction factor for the K-shell. In analogy with (12), L-shell total capture rate will be where denotes a particular L-subshell. The contribution of 1 pertaining to 2 1/2 orbital is the most important, so we keep in our calculations only the contribution of this subshell to the calculated our PSF. The expressions for 0 1 can be obtained from (13) by the replacement of , by 1 , 1 . Electron capture from M-, N-, and higher shells may be defined in a similar fashion, but they have negligible contributions in comparison with K-and L-shells.
Hence, for an allowed transition, the PSF expression of electron capture within the approximation stated above can be written as For / 1 quantities, we used the expression where EC is value of + decay in 2 units, are the binding energies of 1 1/2 and 2 1/2 electron orbitals of the parent nucleus, and are their radial densities on the nuclear surface.
≈ 1 represent the values of the exchange correction. These are due to an imperfect overlap of the initial and final atomic states caused by the one unit charge difference [21]. In our method, we consider these exchange corrections to be unity, for the nuclei considered, the estimated error in doing that being under 1%. Relation 0 = EC − 1 holds. / 1 are the electron bound states, solutions of the Dirac equation (5), and correspond to eigenvalues ( is the radial quantum number). Quantum number is related to total angular momentum = | | − 1/2. These w.f. are normalized such that For simplicity, we consider solutions of Dirac equations , and , that are divided by radial distance . An asymptotic solution is obtained by means of the WKB approximation and by considering that potential is negligibly small: In our calculations, we use number nodes = 0 and = 1, for orbitals 1 1/2 and 2 1/2 , respectively, being −1. Numerically, 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 = 0.
The radial density of the bound-state electron w.f. on the nuclear surface is where = 1.2 1/3 0 is given in fm, 0 being the Bohr radius. For 1 1/2 and 2 1/2 electron orbitals, we use 2 = 2 0,−1 and 2 1 = 2 1,−1 , respectively. For the EC processes, the potential used to obtain the electron w.f. reads Advances in High Energy Physics 5 and charge number = 0 corresponds to the parent nucleus. ( , ) is negative.
The numerical solutions of the Dirac equation were obtained within the power series method of [9], by using similar numerical algorithm as that of [10,11]. The method is able to provide numerical solutions of the Dirac equation for central fields. We provide a grid with values of the potential for different radial distances. The radial w.f. is expanded in an infinite power series that depends on the radial increment and the potential values. The w.f. is calculated step by step in the mesh points. The increment and the number of terms in the series expansion determine the accuracy of the solutions. In our calculations, the increment interval is 10 −4 fm and at least 100 terms are taken into account in the series expansion. These values exceed the convergence criteria of [10]. To renormalize the numerical solutions, we made use of the fact that, at very large distances, the behavior of the w.f. must approach that of the Coulomb function. Therefore, the amplitudes and the phase shifts can be extracted by comparing the numerical solution and the analytical ones. For discrete states, the asymptotic behavior of the w.f. gives a guess for the inverse solutions. The eigenvalue is obtained when the direct solutions and the inverse ones match each other. We constructed an adequate procedure to find the bound states of the electron up to an accuracy of 0.3 keV, or lower, by searching solutions up to 130 keV binding energies. In this range of energies, all the possible bound-state energies are found. We calculated the solutions starting outward from = 0 and inward from a very large value of radius . The bound states should be obtained when both solutions are equal in an intermediate point, for the two components of the wave function. We found these energies by interpolation. We selected radial wave functions , and , that have the same number of nodes = 0 or 1.
For the PSF computation, all integrals in (5) were performed accurately with Gauss-Legendre quadrature in 32 points. We calculated up to 49 values of the radial functions in Q value energy interval, that were interpolated with spline functions.
We also calculated the PSF for EC process using (15) but employing essentially the formalism adopted by [2]. Here, we used the electron radial density (and density ratios) as given in Table 2 of [2]. Exchange corrections were taken as unity. Binding energies were also taken from the same reference.

Results and Discussion
We perform PSF computations for + decay and the EC process with the method described in the previous section that we call TW (this work), for a large number of nuclei of experimental interest.
For + decays, we found previous PSF results computed with approximate methods [15,16], for sixteen nuclei of astrophysical interest. In Table 1, we display the PSF values for these nuclei calculated with our new method (TW) and, for comparison, the values taken from [15,16]. Also, we present the PSF values computed by us using the recipe described in [3]. All calculations were done with 0 value indicated in   [15]. One can see that the agreement between TW results and the other results is in general under 1%, except the last two (heavy) nuclei where the differences reach ∼3%.
In Table 2, we display our computed PSF with the new method for few heavy nuclei, for which we did not find Table 3: Calculated phase space factors EC for electron capture (assuming exchange corrections to be equal to 1). The value of maximum -decay energy is taken from [15] for pure Fermi transitions. The electron densities, their ratios, and binding energies are also provided for orbitals 1 1/2 and 2 1/2 , including those given in [2]. Binding energies are given in units of keV. previous results. For comparison, we computed the same PSF values with the recipe adopted from [3]. 0 -values were taken from [19] for both sets of calculations. We found a rather good agreement between the two sets of results, with differences within, generally, a few percent. There was one exception, 105 Ag, where the difference was large (∼a factor 10). This is a case where 0 -value is very small (0.325 MeV), and this might make our numerical routine inaccurate at such small values. However, this discrepancy may not be so significant, as long as the calculated PSF value is small enough to have little contribution to the corresponding beta-decay rates In Table 3, we preset our results for EC for the same set of nuclei. -values for positron decay were taken from [15] for nuclei marked with ⋆. For the rest of nuclei,values were taken from [19]. Together with the PSF values for EC, the electron densities, , 1 , their ratios, and binding energies for orbitals 1 1/2 and 2 1/2 are also given in Table 3. We compare the results performed with the new method (TW) with those calculated using the recipe of [2]. For these transitions, the differences between the two sets of results are significantly larger than for the positron decays, ranging from a few percent to about a mammoth 35%. We attribute these differences in the calculated PSF values mainly to electron densities, , whose values, calculated with the "old" and "new" methods, differ significantly from each other. We also checked the influence of the screening effect on the PSF values. We found that while, for the positron decays, this effect is very small, for the EC transitions, there are some differences between the "screened" and "unscreened" PSF values. Figure 1 shows this effect on the electron density, , and on the final PSF values. For small values of , the results without screening give PSF values that are 10-15% larger than those listed in Table 3. For heavier nuclei, these differences are only up to 2-3%. The screening effect in PSF calculation is more important for light nuclei and leads to a decrease in the PSF values up to 15%. Finally, in Table 4, we present PSF values for EC transitions recomputed with updated -values taken from [13,14]. We propose to use these new computed values of PSF for calculation of -decay rates.

Summary and Conclusion
In summary, we constructed a new code for computing PSF values for positron decays and EC processes. In our approach, we get positron-free and electron bound w.f. by solving a Dirac equation with a Coulomb-type potential, obtained from a realistic distribution of protons in the daughter nuclei. The FNS and screening effects are addressed as well by our new recipe. Using the same -values, we compare our results with previous calculations where electron/positron w.f. were obtained in an approximate way. For positron decays, the agreement with older results is quite good, while, for EC processes, the differences between "new" and "old" PSF values are as big as 35%. We further found that the screening effect is important for EC processes, especially for light nuclei, having an impact up to 10-15% on the calculated PSF values. Finally, using our new method, we recomputed the PSF for all nuclei using updated -values. We hope that these computed PSF values will prove useful in more accurate estimations of the beta-decay rates. We are currently working on the impact of newly computed PSF values on -decay half-lives and hope to report our findings in the near future.
Finally, we mention that our numerical formalism for getting exact electron/positron w.f. for free and bound states is also suited for the treatment of other charged leptons, such as muons, for example, in any central field. So, it can be also used in applications such as the calculation of muon conversion capture rates. In this respect, there are works where the exact muon w.f. are obtained by solving the Dirac equation in a Coulomb-like potential, within the context of genetic algorithms and neural network techniques [22][23][24]. Within this method, one can also take into account deviations from a pure central Coulomb field by using experimental finite-size charge densities for the attracting nucleons, a procedure which differs from that described in this work.