Phase Space Factors for Double-Beta Decays

REVIEW published: 15 February 2019 doi: 10.3389/fphy.2019.00012 Phase Space Factors for Double-Beta Decays Sabin Stoica 1,2* and Mihail Mirea 1,2 1 International Centre for Advanced Training and Research in Physics, Magurele, Romania, 2 Horia Hulubei National Institute of Physics and Nuclear Engineering, Magurele, Romania Edited by: Filipe Rafael Joaquim, Universidade de Lisboa, Portugal Reviewed by: Hiroyasu Ejiri, Osaka University, Japan Mario Andrés Acero O, University of Atlántico, Colombia *Correspondence: Sabin Stoica [email protected] Specialty section: This article was submitted to High-Energy and Astroparticle Physics, a section of the journal Frontiers in Physics Received: 12 November 2018 Accepted: 21 January 2019 Published: 15 February 2019 Citation: Stoica S and Mirea M (2019) Phase Space Factors for Double-Beta Decays. Front. Phys. 7:12. doi: 10.3389/fphy.2019.00012 Frontiers in Physics | www.frontiersin.org 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 halflives, 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 nonrelativistic 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. Keywords: double-beta decay, phase space factors, nuclear matrix elements, neutrinos, Dirac equation 1 February 2019 | Volume 7 | Article 12 Stoica and Mirea Phase Space Factors for DBD 1. INTRODUCTION 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 halflife 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]: 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/antineutrinos. 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 doublebeta 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νβ − β − ) (A, Z) → (A, Z + 2) + 2e− + 2ν̄ (1) −1  2ν = Gn2ν (E0 , Z)gA4 | me c2 M2ν |2 T1/2 2  −1  0ν ν> = Gn0ν (E0 , Z)gA4 | M0ν |2 <m T1/2 me 1b) Two neutrino double-positron decay + (A, Z) → (A, Z − 2) + 2e + 2ν (2) where Gn2ν(0ν) and M2ν(0ν) are the phase space factors and matrix elements for the 2ν (0ν) decay modes, me and < mν > are the electron and effective neutrino masses, gA 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:  2 gV 0ν M 0ν = MGT − · MF0ν − MT0ν (10) gA 1c) Two neutrino electron capture (EC) positron emitting decay (2νECβ + ) e− + (A, Z) → (A, Z − 2) + e+ + 2ν (3) 1d) Two neutrino double electron capture decay (2νECEC) 2e− + (A, Z) → (A, Z − 2) + 2ν (4) 2) Decays occurring with LNV by two units 2a) Neutrinoless double-electron decay (0νβ − β − ) (A, Z) → (A, Z + 2) + 2e− (5) 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 gA , while the NME depend on this constant through the (gV /gA ) 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 halflives. 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 2b) Neutrinoless double-positron decay (0νβ + β + ) (A, Z) → (A, Z − 2) + 2e+ 2c) Neutrinoless decay (0νECβ + ) electron capture positron (6) emitting e− + (A, Z) → (A, Z − 2) + e+ (7) 2d) Neutrinoless double electron capture decay (0νECEC) 2e− + (A, Z) → (A, Z − 2) (8) 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. Frontiers in Physics | www.frontiersin.org (9) 2 February 2019 | Volume 7 | Article 12 Stoica and Mirea Phase Space Factors for DBD charge Z to a plane wave, evaluated at the origin [2]. They are expressed in an analytical form as follows: 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 gA 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. F NR (Z, ǫ) = (11) 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]. F R (Z, ǫ) = 2(1 + γ )(2pR0 )2(γ −1) eπ η | Ŵ(γ + iη) 2 | Ŵ(2γ + 1) (12) with γ = [1 − (αZ)2 ]1/2 and R0 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 L0 that can be taken from Behrens and Janecke [34]. For electrons L0 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]: F NR (Z, ǫ) Ŵ(2γ + 1)2 1 1 Ŵ(1 + iη) 2 × = × ×| | R 2(γ −1) F (Z, ǫ) L0 (2pR0 ) 4γ Ŵ(γ + iη) (13) 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 nonrelativistic and relativistic Fermi factors can be quite large. For example, in the case of 130 Te→130 Xe decay the use of nonrelativistic 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 238 U →238 Pu. The same expression for the Fermi factors as Equation (8) but in the limit (αZ)2 ≪1 are used by Suhonen and Civitarese [5] when they show how one calculates the PSF for DBD. 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 Coulombtype 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 Frontiers in Physics | www.frontiersin.org 2πη = eπ η | Ŵ(1 + iη) |2 1 − e−2π η 3 February 2019 | Volume 7 | Article 12 Stoica and Mirea Phase Space Factors for DBD (Coulomb) potential: 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 61 (peff R0 )2 , where peff is the “effective” momentum of the leptons at the nuclear surface (peff = 3Zα/2R0 for electrons and peff = 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. + (r) = 9ǫκµ  µ gκ (ǫ, r)χκ µ ifκ (ǫ, r)χ−κ  (15)  (16) for β − decay and − 9ǫκµ =  −µ ifκ (ǫ, r)χ−κ −µ −gκ (ǫ, r)χκ for β + decay where κ = (l − j)(2j + 1) is the relativistic quantum µ number and χκ 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: dgκ (ǫ, r) κ +1 ǫ − V + me c 2 =− gκ (ǫ, r) + fκ (ǫ, r) dr r ch̄ ǫ − V − me c 2 κ −1 dfκ (ǫ, r) =− gκ (ǫ, r) + fκ (ǫ, r) dr ch̄ r 3. RECENT APPROACHES (17) 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: 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 Zd , V(r) = −Zd (α h̄c)/r. The FNS effect is taken into account with a Coulomb potential given by an uniform charge distribution in a sphere of radius RA = r0 A1/3 with r0 = 1.2 fm [11, 37]:  q   e c2 sin(kr − l π − η ln(2kr) + δ ) h̄e−iδk  ǫ+m gk (ǫ, r) k  2ǫ 2 q ∼ fk (ǫ, r) r→∞ ǫ−me c2 pr cos(kr − l π2 − η ln(2kr) + δk ) 2ǫ (18) 4.1. Electron Radial Functions Here c is the p light velocity, me /ǫ are the electron mass/energy, k = p/h̄ = ǫ 2 − (me c2 )2 /(ch̄) is the electron wave number, η = Zeff e2 /h̄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 effectivepcharge number is Zeff = ±Z for β ∓ decays. The energy is ǫ = (me c2 )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 δC = arg[ζ (ǫ +m0p c2 )−i(κ + π λ)kh̄c]−(λ−l−1) 2 +arg Ŵ(λ−iη)−Sζ ,κ π, with λ = κ 2 − ζ 2 , ζ = ±Zeff α, 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  b µ  gn,κ (r)χκ 9ǫbn κµ (r) = (19) b (r)χ µ ifn,κ −κ For free states we use relativistic scattering electron/positron wave functions, solutions of the Dirac equation in a central are solutions of the Dirac equation (17) and correspond to the eigenvalues ǫn (n is the radial quantum number). The quantum V(Z, r) = ( − Zαrh̄c ,   r ≥ RA , 2 A) , r < RA , −Z(α h̄c) 3−(r/R 2R (14) 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. Frontiers in Physics | www.frontiersin.org 4 February 2019 | Volume 7 | Article 12 Stoica and Mirea Phase Space Factors for DBD number κ is related to the total angular momentum jκ =| b (r) and κ | −1/2. For simplicity, the quantities gn,κ (r) = rgn,κ b fn,κ (r) = rfn,κ (r) are used in the following. These wave functions are normalized such that Z ∞ 2 2 [gn,κ (r) + fn,κ (r)]dr = 1. (20) 0 An asymptotic solution is obtained by means of the WKB approximation and by considering that the potential V is negligible small: fn,κ ch̄ = gn,κ ǫ + me c 2  ′ gn,κ κ + gn,κ r  (21) where ′ gn,κ 1 = − µ′ µ−1 − µ gn,κ 2 (22) with µ=  ǫ + me c 2 h̄2 c2 (V − ǫ + me c2 ) + κ2 r2 1/2 . FIGURE 1 | (A) Realistic proton density ρe for 150 Nd represented in cylindrical coordinates ρ and z. (B) Profile of the realistic proton density ρe for 150 Nd (thick line) compared with that given with the constant density approximation (dot-dashed line). (23) In our calculations we use n=0 and n=1 number of nodes, for the orbitals 1s1/2 and 2s1/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: V(Z, r) = α h̄c Z ρe (rE′ ) E′ dr | Er − rE′ | (24) where the charge density is ρe (Er) = X (2ji + 1)vi2 | 9i (Er) |2 (25) i 9i is the proton (Woods-Saxon) wave function of the spherical single particle state i and vi is its occupation amplitude. The factor (2ji + 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 150 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. Frontiers in Physics | www.frontiersin.org 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 150 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 5 February 2019 | Volume 7 | Article 12 Stoica and Mirea Phase Space Factors for DBD function φ(r), which is the solution of the Thomas Fermi equation: √ d2 φ/dx2 = φ 3/2 / x, (26) with x = r/b, b ≈ 0.8853a0 Z−1/3 and a0 = 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 ThomasFermi 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 φβ − β − (∞) = 2/Zd , Zd = Z0 − 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: rVβ − β − (Z, r) = (rV(Z, r) + 2) × φ(r) − 2 FIGURE 3 | The effective screening function φeff as function of the dimensionless parameter x calculated for the 136 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. 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 = Z0 − 1. In the case of the ECEC process, the potential used to obtain the electron wave function is: VECEC (Z, r) = V(Z, r) × φ(r) (27) and Z = Z0 , the final system being neutral. In Figure 3, the effective screening φeff given by the ratios Vβ − β − (x)/V(x), Vβ + β + (x)/V(x), Vβ + (x)/V(x), VEC (x)/V(x) and VECEC (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. 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 α h̄c=1. The asymptotic potential between an electron and a double ionized atom is rVβ − β − = −2. In this case, the charge number Z = Z0 +2 corresponds to the daughter nucleus, Z0 being the charge number of the parent nucleus. In the case of the β + β + process, the potential used to obtain the positron wave functions is: rVβ + β + (Z, r) = (rV(Z, r) + 2) × φ(r) − 2 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: (28) dW2ν = (a(0) + a(1) cos θ12 )w2ν dω1 dω2 dǫ1 dǫ2 d cos θ12 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 = Z0 − 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: rVEC (Z, r) = (rV(Z, r) + 1) × φ(r) − 1 Frontiers in Physics | www.frontiersin.org (32) where θ12 is the angle between the electrons. Both parameters a0 and a1 represents summation over intermediate nuclear states of products between electron phase factors and nuclear matrix elements summed over intermediate states N a(0) (29) and the charge number Z = Z0 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: rVβ + (Z, r) = (rV(Z, r) + 1) × φ(r) − 1 (31) a(1)   2 2 1 X 1 (0)  X MN (KN + LN ) + MN (KN − LN )  = f11 4 3 N N (33)   2 2 1 (1)  X 1 X = f11 MN (KN + LN ) − MN (KN − LN )  4 9 N (30) N (34) 6 February 2019 | Volume 7 | Article 12 Stoica and Mirea Phase Space Factors for DBD (i) Here, f11 represent the electron phase factors and the matrix elements between initial (I) and final (F) states are denoted + + + + + MN = h0+ F || τ σ || 1N ih1N || τ σ || 0I i (1) F2ν = (44) that correspond to the two parameters a(0) and a(1) of Equation (32). An energy dependent angular correlation α is defined as (35) Also, (gA G cos θC )4 2 2 w2ν = ω1 ω2 (p1 c)(p2 c)ǫ1 ǫ2 δ(ǫ1 +ǫ2 +ω1 +ω2 −W0 ) 64π 7 h̄ (36) with W0 = Qββ + 2m0 c2 . The previous expression depends on the electron energies ǫ1 , ǫ2 , electron momenta p1 , p2 and on the neutrino energies ω1 , ω2 . Usual values of the parameters can be considered as G = 1.16637 × 10−5 GeV−2 (the Fermi coupling constant), cos θC =0.9737 (θC standing for the CabbidoKobayashi-Maskawa mixing angle between d and s quarks) and gA ≈1 (the axial vector form factor). By considering that the excitations energies are of the order of giant Gamow-Teller resonance à = 1.12A1/2 MeV, the closure approximation is adopted. That is, the energy of the intermediate state EN is replaced with an average value hEN i. Therefore the energy denominators 1 1 + ǫ 1 + ω 1 + EN − EI ǫ 2 + ω 2 + EN − E I 1 1 + = ǫ 1 + ω 2 + EN − EI ǫ 2 + ω 1 + EN − E I KN = (37) LN (38) 2 R R Q +2me c2 −ǫ1 R Qββ +2me c2 −ǫ1 −ǫ2 2Ã2 Qββ +me c dǫ1 m ββ dǫ2 m c2 dω1 2 9 ln 2 me c2 ec e (1) 2 2 ×f11 w2ν (2(hKN i + hLN i ) + 5hKN ihLN i) α= (1) /dǫ1 dF2ν (45) (0) dF2ν /dǫ1 from (0) dF d2 W2ν = |M2ν |2 2ν (1 + α cos θ12 ) d cos θ12 dǫ1 dǫ1 (46) In this work, we report normalized phase factors, G(i) 2ν = (i) F2ν /[gA4 (me c2 )2 ] in units yr−1 . We mention that by this renormalization the quantities G(i) do not depend on the axialvector constant gA . Also, as one can see from the formulas, the PSF (G(i) (2ν,0ν) depend very little to the average energy < EN >, 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 (0) f11 =| f −1−1 |2 + | f11 |2 + | f −1 1 |2 + | f1 −1 |2 are changed within averaged values (47) and 1 W0 + hEN i − EI hKN i 21 2 W0 + E N − E I (39) 1 W0 + hEN i − EI LN ≈ hLN i 21 2 W0 + E N − E I (40) KN ≈ (1) ∗ + f −1 1 f1 −1∗ = −2ℜ f −1−1 f11 f11 with the products f −1−1 = g−1 (ǫ1 )g−1 (ǫ2 ) ; f11 = f1 (ǫ1 )f1 (ǫ2 ), f −11 = g−1 (ǫ1 )f1 (ǫ2 ) ; f1−1 = f1 (ǫ1 )g1 (ǫ2 ) In this way, a factorization of the Equation (32) can be realized. By considering à ≈ 21 W0 + hEN i − EI and the nuclear matrix element on the form M2ν = X h0+ || τ + σ || 1+ ih1+ || τ + σ || 0+ i F N N I N 1 2 W0 + E N − EI (0) F2ν = g−1 (ǫ) = g−1 (ǫ, RA ) ; f1 (ǫ) = f1 (ǫ, RA ) (41) (49) (50) (51) where RA = 1.2A1/3 fm is the nuclear radius. For the 2νββ decay mode and transitions to ground states (g.s.), the PSF expression reads: ββ G2ν (0+ → 0+ ) = (42) Z Qββ +me c2 2Ã2 dǫ1 3 ln 2gA4 (me c2 )2 me c2 Z Qββ +2mce2 −ǫ1 −ǫ2 Z Qββ +2me c2 −ǫ1 dω1 dǫ2 × me c 2 (0) ×f11 w2ν (hKN i2 2 R R Q +2me c2 −ǫ1 R Qββ +2me c2 −ǫ1 −ǫ2 2Ã2 Qββ +me c dǫ1 m ββ dǫ2 m c2 dω1 2 3 ln 2 me c2 ec e (0) ×f11 w2ν (hKN i2 + hLN i2 + hKN ihLN i) 0 + hLN i2 + hKN ihLN i) (52) where Qββ = M(A, Z0 ) − M(A, Z0 − 2) − 4me c2 is the kinetic energy released in the process. hKN i, hLN i are expressions that depend on the electron/positron (ǫ1,2 ) and (43) Frontiers in Physics | www.frontiersin.org (48) The values of the f and g functions are approximated with the solutions at the nuclear surface (the method I from [11]). 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]: dW2ν (0) (1) = |M2ν |2 F2ν + |M2ν |2 F2ν cos θ12 d cos θ12
7 February 2019 | Volume 7 | Article 12 Stoica and Mirea Phase Space Factors for DBD In Equations (52, 56, 57), it is convenient to redefine the PSF by a renormalization that eliminates the constant gA and correlates (by dividing by 4R2A ) the dimension of G0ν with the NME which are dimensionless. These are reflected in the lifetimes formulas (9). Thus, the PSF are also reported in [yr−1 ]. A similar expression is employed in the PSF calculation for the transitions ββ by Qββ (0+ ). The formula to excited 0+ 1 1 states, but replacing Q 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 EN in the intermediate odd-odd nucleus by a suitable average energy. 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]. 1 1 + ǫ1 + ω1 + hEN i − EI ǫ2 + ω2 + hEN i − EI 1 1 hLN i = + ǫ1 + ω2 + hEN i − EI ǫ2 + ω1 + hEN i − EI hKN i = (53) (54) Here, the difference in energy in the denominator can be obtained from the approximation Ã2 = [W0 /2 + hEN i − EI ]2 , where the empirical systematics à = 1.12A1/2 (in MeV) [2] gives approximately the energy of the giant Gamow-Teller resonance in the intermediate nucleus [11]. The quantity W0 is related to the Qββ value of the process and w2ν = gA4 (G cos θC )4 2 2 ω1 ω2 (p1 c)(p2 c)ǫ1 ǫ2 , 64π 7 h̄ 4.3.2. The ECβ + Case For the ECβ + decays the energy released in the process is QECβ = M(A, Z0 ) − M(A, Z0 − 2) − 2mce2 . If the numerical solutions of the Dirac equation are obtained in Bohr units a0 , the probability that an electron is found on the surface of a nucleus of radius RA can be defined as:  3  2 1 h̄c a0 2 2 2 Bn,κ = [gn,κ (RA ) + fn,κ (RA )] (62) 4π(me c2 )3 a0 RA (55) where ω1 and ω2 = Qββ − ǫ1 − ǫ2 − ω1 + 2me c2 are the neutrino energies. The PSF are finally renormalized to the electron rest energy and are reported in [yr−1 ]. The PSF for the 2νββ decay mode and transitions to excited 0+ 1 states is calculated with a formula similar to (52), but replacing ββ − E (0+ ), which is the kinetic energy Qββ by Q(0+ x 1 1) = Q + released in this transition. Ex (0+ 1 ) is the energy of the excited 01 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]: ββ G2ν (0+ → 2+ 1) = 2Ã6 Z The PSF expression for 2νECβ decay mode is × 2 Qββ (2+ 1 )+me c (56) + where Qββ (2+ 1 ) = Q − Ex (21 ). For the 0νββ decay and transitions to g.s. the PSF reads: ββ G0ν (0+ → 0+ ) = 2 4gA4 R2A ln 2 Z Qββ +me c2 me c 2 (0) f11 w0ν dǫ1 ECβ + (57) G0ν where w0ν = 0 (63) = 1 2 (G cos θ )4 2 (h̄c )(me c2 )5 4R2A ln 2 4π 3 X 2 B2i,−1 [g−1 (ǫp,i ) + f12 (ǫp,i )]pp cǫp × (64) i=0,1 gA4 (G cos θC )4 (me c2 )2 (h̄c2 )(p1 c)(p2 c)ǫ1 ǫ2 16π 5 where ǫp,i denotes the maximal value of the positron associated to the state i. (58) (0) For transitions on the 2+ 1 state, instead of f11 in Equation (57), the next expression is integrated f1+ + f1− = 3( me c 2 where ǫn,κ are the binding energies of the electron while pp and ǫp are the linear momentum and the energy of the positron. Here, the expressions for hKN i and hLN i are similar to those c from Equations (53)–(54), but where ǫ1 is replaced by ǫi,−1 = 2 me c −ǫi,−1 , the energy of the captured electron and ǫ2 is replaced by ǫp , the energy of the emitted positron. For the 0νECβ decay process the PSF expression is: 0 (0) ×f11 w2ν (hKN i − hLN i)2 = 2 ×[g−1 (ǫp ) + f12 (ǫp )] ×(hKN i2 + hLN i2 + hKN ihLN i)ω12 ω22 pp cǫp dω1 dǫp dǫ1 ln 2gA4 (me c2 )6 me c2 Z Qββ (2+ )+2me c2 −ǫ1 −ǫ2 Z Qββ (2+ )+2me c2 −ǫ1 1 1 dω1 dǫ2 × me c 2 X 2A2 (G cos θ )4 B2i,−1 (me c2 ) 5 3 ln 2 16π h̄ i=0,1 Z QECβ +ǫi,−1 +me c2 Z QECβ +ǫi,−1 −ǫp +me c2 ECβ + G2ν 4.3.3. The 2νECEC Case The PSF expression is defined as: h̄c )2 (|f −2−1 |2 + |f −1−2 |2 + |f21 |2 + |f12 |2 ). m e c 2 RA (59) f −2−1 = g−2 (ǫ1 )g−1 (ǫ2 ) ; f21 = f2 (ǫ1 )f1 (ǫ2 ), f −1−2 = g−1 (ǫ1 )g−2 (ǫ2 ) ; f12 = f1 (ǫ1 )f2 (ǫ2 ) Frontiers in Physics | www.frontiersin.org GECEC = 2ν (60) X 2Ã2 (G cos θ )4 B2i,−1 B2j,−1 (me c2 )4 × 3 3 ln 2 16π h̄ i,j=0,1 Z QECEC +ǫi,−1 +ǫj,−1 (hKN i2 + hLN i2 × 0 + hKN ihLN i)ω12 ω22 dω1 (61) 8 (65) February 2019 | Volume 7 | Article 12 Stoica and Mirea Phase Space Factors for DBD where QECEC = M(A, Z0 ) − M(A, Z0 − 2) is the energy released in the process. The expressions for hKN i and hLN i are similar to those from Equations (53)–(54), but where ǫ(1,2) are replaced by c ǫ(i,j)−1 = me c2 − ǫ(i,j)−1 , the energies of the captured electrons. 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 me c2 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. 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−4 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 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 TABLE 1 | PSF for 2νβ − β − decays to final g.s. Nucleus β−β− β−β− Qg.s. G2ν (g.s.) (10−21 yr−1 ) (MeV) [39] [11] [3, 35, 36] [5] [46] 48 Ca 4.267 15,536 15,550 16,200 16,200 14,805 76 Ge 2.039 46.47 48.17 53.8 52.6 45.1 82 Se 2.996 1,573 1,596 1,830 1,740 1,503 96 Zr 3.349 6,744 6,816 7,280 6,420 100 Mo 3.034 3,231 3,308 110 Pd 2.017 132.5 137.7 116 Cd 2.813 2,688 2,764 128 Te 0.8665 0.2149 0.2688 0.35 0.344 130 Te 2.528 1,442 1,529 1,970 1,940 136 Xe 2.458 1,332 1,433 2,030 1,980 1,337 150 Nd 3.371 35,397 36,430 48,700 48,500 34,675 238 U 1.144 98.51 14.57 Frontiers in Physics | www.frontiersin.org 9 3,860 3,600 3,106 127.8 2,990 2,588 1,427 February 2019 | Volume 7 | Article 12 Stoica and Mirea Phase Space Factors for DBD TABLE 2 | PSF for 0νβ − β − decays to final g.s. β−β− Nucleus β−β− Qg.s. G0ν (MeV) 48 Ca 4.267 76 Ge 2.039 82 Se 2.996 [39] [11] 24.65 [3, 35, 36] 24.81 2.372 96 Zr 3.349 20.48 20.58 3.034 15.84 15.92 [47] 26.0 2.62 10.16 100 Mo [5] 26.1 2.363 10.14 (g.s.) (10−15 yr−1 ) [46] 24.83 2.55 24.55 2.37 2.28 10.18 9.96 11.4 11.1 23.1 20.62 20.45 18.7 45.6 15.95 15.74 4.83 4.66 18.9 16.73 16.57 110 Pd 2.017 116 Cd 2.813 128 Te 0.8665 130 Te 2.528 14.24 14.22 19.4 16.7 14.25 14.1 136 Xe 2.458 14.54 14.58 19.4 17.7 14.62 14.49 150 Nd 3.371 61.94 63.03 85.9 78.4 63.16 66.0 238 U 1.144 32.53 33.61 4.915 4.815 16.62 16.70 0.5783 0.5878 0.748 0.671 TABLE 3 | PSF for β − β − decays to final excited 0+ 1 states. Nucleus Q β−β− 0+ 1 β−β− G2ν (MeV) [39] β−β− ) (10−21 yr−1 ) (0+ 1 [11] G0ν [5] [46] ) (10−15 yr−1 ) (0+ 1 [39] [11] 48 Ca 1.270 0.3518 0.3627 0.376 0.343 0.3041 0.2989 76 Ge 0.9171 0.06129 0.06978 0.0769 0.064 0.1932 0.1776 0.9440 82 Se 1.508 4.170 4.80 4.194 96 Zr 2.201 169.4 175.4 190 163.38 4.594 4.566 100 Mo 1.904 57.08 60.55 101 56.21 3.168 3.162 110 Pd 0.5472 116 Cd 1.056 3.3 ×10−3 4.8 ×10−3 0.7590 0.8737 4.3 ×10−3 0.89 0.8 0.1223 0.08844 0.7585 0.7163 130 Te 0.7335 0.05460 0.07566 18.6 0.069 0.3651 0.3086 136 Xe 0.8790 0.2823 0.3622 0.485 0.31 0.6746 0.6127 4850 4064 150 Nd 2.631 4116 4329 238 U 0.2032 1.5 ×10−4 4.6 ×10−4 26.96 27.27 0.8229 0.7534 TABLE 4 | PSF for β − β − decays to final excited 2+ 1 states. Nucleus Q β−β− β−β− 2+ 1 G2ν (MeV) [39] β−β− ) (10−21 yr−1 ) (2+ 1 [37] [5] [46] G0ν ) (10−15 yr−1 ) (2+ 1 [39] [37] 48 Ca 3.284 4074 4410 4400 3816 57.09 76 Ge 1.480 0.384 0.48 0.49 0.40 1.66 82 Se 2.219 69.6 90.6 85 71.16 12.13 96 Zr 2.571 745.5 850 730.8 33.87 100 Mo 2.494 569.0 690 585 32.1 110 Pd 1.359 0.46 0.46 2.41 116 Cd 1.520 1.88 128 Te 0.4255 130 Te 1.990 6.8 ×10−7 136 Xe 1.640 7.68 150 Nd 3.037 30308 238 U 1.099 2.66 Frontiers in Physics | www.frontiersin.org 79.6 2.3 1.36×10−6 116 1.3 ×10−6 120 15 45600 49000 2.11 1.84 13.8 4.28 0.049 81.09 18.34 9.03 8.31 31964 60.4 223 0.067 22.8 301 26.3 10 February 2019 | Volume 7 | Article 12 Stoica and Mirea Phase Space Factors for DBD TABLE 5 | PSF for β + β + decay mode. Nucleus β+β+ + + Qβ β G2ν β+β+ (10−29 yr−1 ) G0ν [12] (10−20 yr−1 ) (MeV) [39] [12] [35, 36] [39] [35, 36] 78 Kr 0.8023 9,159 9,770 13,600 243.2 250 293 96 Ru 0.6706 942.3 1,040 1,080 80.98 84.5 90.7 106 Cd 0.7314 1,794 2,000 1,970 91.75 92.6 102 124 Xe 0.8203 4,261 4,850 4,770 107.8 114 123 130 Ba 0.5748 91.54 110 47.9 23.82 25.7 136 Ce 0.3345 0.205 0.267 0.559 2.13 21 2.42 3.55 TABLE 6 | PSF for ECβ + decay mode. Nucleus QECβ + (MeV) ECβ + G2ν ECβ + (10−24 yr−1 ) ǫ0,−1 ǫ1,−1 (keV) (keV) [39] [12] [35, 36] G0ν (10−18 yr−1 ) [39] [12] 78 Kr 1.824 17.7 3.1 338 385 464 6.90 6.34 96 Ru 1.693 26.2 4.9 372 407 454 11.30 9.62 [35, 36] 7.11 10.8 106 Cd 1.753 31.1 5.9 741 702 779 15.39 13 124 Xe 1.842 39.4 7.8 1,235 1,530 1,720 17.10 19.7 22.9 130 Ba 1.597 42.4 8.5 740 580 549 22.89 17.6 19.8 136 Ce 1.357 45.6 9.2 236 190 253 21.45 15.3 18.7 1.62 ×10−6 1.16 ×10−6 0.376 0.0887 1.16 1.62 1.21 2.21 ×10−8 3.81 ×10−9 0.546 0.0507 0.9558 0.230 0.782 0.729 3.07 1.94 3.53 1.92 50 Cr 0.1469 8.3 1.2 58 Ni 0.9043 11.1 1.2 64 Zn 0.07269 12.6 2.0 74 Se 0.1872 16.0 2.7 84 Sr 0.7677 19.7 3.5 92 Mo 0.6298 23.9 4.4 0.245 0.206 102 Pd 0.1499 28.6 5.4 112 Sn 0.8978 33.7 6.5 1.40 ×10−5 1.62 ×10−6 120 Te 0.7082 36.5 7.1 1.406 0.730 7.61 3.92 144 Sm 0.7604 52.3 10.7 3.56 2.49 10.38 8.11 1.037 2.2 ×10−5 6.31 1.11 1.09 ×10−5 4.95 2.101 0.287 8.41 5.20 156 Dy 0.9840 59.6 12.4 48.1 25.3 162 Er 0.8250 63.4 13.3 15.4 6.4 168 Yb 0.3873 67.4 14.0 0.07689 71.6 15.2 9.8 ×10−3 16.48 174 Hf 9.1 ×10−3 10.81 0.0272 184 Os 0.4289 80.4 17.3 0.321 0.0299 26.21 7.04 190 Pt 0.3625 85.2 18.2 0.00588 27.61 5.57 7.5 ×10−6 0.127 1 ×10−9 26.01 15.2 25.03 12.9 1.30 15.8 4.23 of results, ranging in the interval (10–40)%. The difference for the G2ν in the case of 238 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 G0ν and G2ν . 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 the electron wave functions at the nuclear surface and without inclusion of the screening effect. The PSF values (G2ν , G0ν ) for β − β − decays are presented + in Tables 1–4 for transitions to g.s and excited 0+ 1 and 21 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 128 Te the difference betwee the G2ν values amounts to 22%. In the case of the 238 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 (76 Ge, 110 Pd, 116 Cd, 130 Te, 136 Xe) with notable differences between the two sets Frontiers in Physics | www.frontiersin.org 6.69 14.7 11 February 2019 | Volume 7 | Article 12 Stoica and Mirea Phase Space Factors for DBD TABLE 7 | PSF for 2νECEC decay mode. Nucleus QECEC ǫ0,−1 ǫ1,−1 (MeV) (keV) (keV) GECEC (10−24 yr−1 ) 2ν [39] [12] [35, 36] 78 Kr 2.846 17.7 3.1 410 660 774 96 Ru 2.715 26.2 4.9 1,450 2,400 2,740 106 Cd 2.775 31.1 5.9 4,299 5,410 6,220 124 Xe 2.864 39.4 7.8 15,096 17,200 20,200 130 Ba 2.619 42.4 8.5 14,773 15,000 16,300 136 Ce 2.379 45.6 9.2 12,223 12,500 15,800 50 Cr 1.169 8.3 1.2 0.238 0.422 58 Ni 1.926 11.1 1.2 9.90 15.3 64 Zn 1.095 12.6 2.0 1.03 1.41 74 Se 1.209 16.0 2.7 3.41 5.65 84 Sr 1.790 19.7 3.5 64.62 93.6 92 Mo 1.652 23.9 4.4 82.32 208 102 Pd 1.172 28.6 5.4 42.09 46 112 Sn 1.920 33.7 6.5 869.7 1,150 120 Te 1.730 36.5 7.1 840.3 888 144 Sm 1.782 52.3 10.7 6,436 5,150 156 Dy 2.006 59.6 12.4 22,078 17,600 162 Er 1.847 63.4 13.3 20,085 15,000 168 Yb 1.409 67.4 14.0 7,872 4,710 174 Hf 1.099 71.6 15.2 3,432 1,580 184 Os 1.451 80.4 17.3 24,222 12,900 190 Pt 1.384 85.2 18.2 28,153 12,900 2.9 ×10−4 36 Ar 0.4326 5.0 1.2 40 Ca 0.1935 5.9 1.2 54 Fe 0.6798 9.7 1.4 108 Cd 0.2718 31.1 126 Xe 0.9195 132 Ba 1.02 ×10−5 1.25 ×10−5 5.9 0.0682 0.0207 39.4 7.8 60.59 46.1 0.8439 42.3 7.7 61.98 39.1 138 Ce 0.6930 45.6 9.2 34.47 18.4 152 Gd 0.05570 55.9 11.6 1.12 ×10−2 158 Dy 0.2829 59.6 12.4 164 Er 0.02506 63.4 13.3 180 W 0.1433 75.8 16.0 196 Hg 0.8206 89.9 19.5 0.03021 3.191 8.3 ×10−3 18,100 0.0469 0.183 1.4781 0.00156 3587 821 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 QECβ 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. 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 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. The PSF values for the β + β + decays are displayed in Table 5. Differences of about 11–30% between G2ν values of Neacsu and Horoi [11] and Mirea et al. [39] are retrieved, in majority of cases. The agreement concerning the G0ν values is very good with one exception: in the case of the 136 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 Frontiers in Physics | www.frontiersin.org 17 12 February 2019 | Volume 7 | Article 12 Stoica and Mirea Phase Space Factors for DBD (0) (0) FIGURE 4 | The normalized energy distributions dG2ν /dǫ1 , dG2ν /d(ǫ1 + ǫ2 ) and the angular correlation parameter α as function of the electron energies ǫ1 and ǫ2 for 82 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 | Same as Figure 4 for the 136 Xe parent nucleus. FIGURE 6 | Same as Figure 4 for the 150 Nd parent nucleus. 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. are enumerated. The third sector includes the nuclei which can undergo only 2νECEC decays. The differences between the G2ν 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 Frontiers in Physics | www.frontiersin.org 13 February 2019 | Volume 7 | Article 12 Stoica and Mirea Phase Space Factors for DBD 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 136 Xe, 82 Se, and 150 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. 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 me c2 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−4 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 1s1/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. 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 136 Xe, 82 Se and 150 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. 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 (0) dG(0) 2ν /dǫ1 , dG2ν /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 136 Xe nucleus. Due to the screening, the values of the differential 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, CNCSUEFISCDI, project number PN-III-P4-ID-PCE-2016-0078. REFERENCES 2. Haxton WC, Stephenson GJ. Double beta decay. Prog Part Nuclear Phys. (1984) 12:409–79. 3. Doi M, Kotani T, Takasugi E. Double beta Decay and Majorana neutrino. Prog Theor Phys Suppl. (1985) 83:1–175. 1. Primakoff H, Rosen SP. Double beta decay. Rep Prog Phys. (1959) 22:121–66. Frontiers in Physics | www.frontiersin.org 14 February 2019 | Volume 7 | Article 12 Stoica and Mirea Phase Space Factors for DBD 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 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. 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 33. Haxton WC, Stephenson GJ, Strottman D. Lepton-number conservation and the double-β decay of 130 Te. Phys Rev. (1982) 25:2360–9. 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. 36. Doi M, Kotani T. Neutrinoless modes of double beta decay. Prog Theor Phys. (1993) 89:139–59. 37. Doi M, Kotani T, Nishiura H, Takasugi E. Double beta decay. Prog Theor Phys. (1983) 69:602–35. 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. 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 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. 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 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 43. Esposito S. Majorana solution of the Thomas-Fermi equation. Am J Phys. (2002) 70:852–6. doi: 10.1119/1.1484144 44. Doi M, Kotani T, Nisiura H, Okuda K, Takasuki E. Neutrino mass, the righthanded interaction and the double beta decay. I: formalism. Prog Theor Phys. (1981) 66:1739–64. 45. Buhring W. Computational improvements in phase shift calculations of elastic electron scattering. Zeitschrift fur Physik (1965) 187:180–96. 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 47. Stefanik D, Dvornicky R, Simkovic F, Vogel P. The light neutrino exchange mechanism of the 0νββ-decay with left-and righthanded leptonic and currents revisited. Phys Rev. (2015) 92:055502. doi: 10.1103/PhysRevC.92.055502 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 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 50. Arnold R, Augier C, Baker JD , Barabash AS, Basharina-Freshville A, Blondel S, et al. Search for neutrinoless double-beta decay of 100 Mo with the NEMO-3 detector. Phys Rev D (2014) 89:111101. doi: 10.1103/PhysRevD.89. 111101 4. Tomoda T. Double beta decay. Rep Prog Phys. (1991) 54:53–126. 5. Suhonen J, Civitarese O. Weak-interaction and nuclear-structure aspects of double beta decay. Phys Rep. (1998) 300:123–214. 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 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 8. Rodejohann W. Neutrino-less double beta decays and particle physics. Int J Modern Phys. (2011) 20:1833–930. doi: 10.1142/S0218301311020186 9. Vogel P. Nuclear structure and double beta decay. J Phys G. (2012) 39:124002. doi: 10.1088/0954-3899/39/12/124002 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 11. Kotila J, Iachello F. Phase-space factors for double-β decay. Phys Rev. (2012) 85:034316. 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 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 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 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 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 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 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 19. Kortelainen M, Suhonen J. Improved short-range correlations and 0νββ nuclear matrix elements of 76 Ge and 82 Se. Phys Rev. (2007) 75:051303(R). doi: 10.1103/PhysRevC.75.051303 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 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. 22. Retamosa J, Caurier E, Nowacki F. Neutrinoless double beta decay of 48 Ca. Phys Rev. (1995) 51:371–8. 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 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 25. Horoi M, Stoica S, Brown BA. Shell-model calculations of two-neutrino double-β decay rates of 48 Ca with the GXPF1A interaction. Phys Rev. (2007) 75:034303. doi: 10.1103/PhysRevC.75.034303 26. Horoi M, Stoica S. Shell-model analysis of the neutrinoless double-β decay of 48 Ca. Phys Rev. (2010) 81:024321. doi: 10.1103/PHYSREVC.81.024321 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 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 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 Frontiers in Physics | www.frontiersin.org 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. Copyright © 2019 Stoica and Mirea. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms. 15 February 2019 | Volume 7 | Article 12