Phase space factors for double-\beta decay

Phase space factors for double- decay


A complete and improved calculation of phase space factors (PSF) for and decay is presented. The calculation makes use of exact Dirac wave functions with finite nuclear size and electron screening and includes life-times, single and summed electron spectra, and angular electron correlations.

23.40.Hc, 23.40.Bw, 14.60.Pq, 14.60.St

I Introduction

Double- decay is a process in which a nucleus decays to a nucleus by emitting two electrons (or positrons) and, usually, other light particles


Double- decay can be classified in various modes according to the various types of particles emitted in the decay.

Figure 1: Double- decay mechanism for (a) two-neutrino, (b) neutrinoless and (c) neutrinoless decay with Majoron emission.

For , the process , Fig. 1a,


is allowed by the standard model and expected to occur with calculable probability. In recent years, the process , Fig. 1b,


has become of great interest, due to the discovery of neutrino oscillation SUP01 (); SNO2002 (); KAM2003 (). The process is of utmost importance for obtaining the neutrino mass since its decay probability is proportional to the square of the average neutrino mass . A third process has been also considered, , Fig. 1c,


in which a massless Nambu-Goldstone boson, called a Majoron, is emitted. However, most of the interest in this mode has disappeared in recent years and hence it will not be considered here. For decay, the corresponding modes , , are


In this case, there are also the competing modes in which either one or two electrons are captured from the electron cloud, , , , .

For processes allowed by the standard model (, , ) the half-life can be, to a good approximation, factorized in the form


where is a phase space factor and the nuclear matrix element. For processes not allowed by the standard model the half-life can be factorized as


where is a phase space factor, the nuclear matrix element and contains physics beyond the standard model through the masses and mixing matrix elements of neutrino species. For both processes, two crucial ingredients are the phase space factors and the nuclear matrix elements. Recently, we have initiated a program for the evaluation of both quantities. For the nuclear matrix elements we have developed an approach based on the microscopic interacting boson model (IBM-2) and presented some results in barea (). Additional preliminary results have been presented in barea11 (); iac11 () and will be discussed in a forthcoming publication barea12 (). In this article, we concentrate on phase space factors.

A general theory of phase space factors in DBD was developed years ago by Doi et al. doi81 (); doi83a () following previous work of Primakoff and Rosen primakoff () and Konopinski konopinski (). It was reformulated by Tomoda tom91 () whose work we follow here. Tomoda also presented results in a selected number of nuclei. These results were obtained by approximating the electron wave functions at the nuclear radius and without inclusion of electron screening. In this article we take advantage of some recent developments in the numerical evaluation of Dirac wave functions and in the solution of the Thomas-Fermi equation to calculate more accurate phase space factors for double- decay in all nuclei of interest. Our results are of particular interest in heavy nuclei, large, where relativistic and screening corrections play a major role. Studies similar to ours were done for single- decay in the 1970’s single (). In this article we report results for , which at the moment is the most promising decay mode. In a subsequent publication, we will present results for , , , which is very recently attracting some attention positron ().

Ii Electron wave functions

The key ingredients for the evaluation of phase space factors in single- and double- decay are the (scattering) electron wave functions. (For EC the bound wave functions.) The general theory of relativistic electrons can be found e.g., in the book of Rose rose (). We use, for decay, positive energy Dirac central field wave functions,


where are spherical spinors and and are radial functions, with energy , depending on the relativistic quantum number defined by . Given an atomic potential the functions and satisfy the radial Dirac equations:


The electron scattering wave function, denoted here by , where is the projection of the spin, can then be expanded in terms of spherical waves as




The large and small components and , respectively, with of the radial wave functions are normalized so that they asymptotically oscillate with




is the electron wave number, is the Sommerfeld parameter and is the phase shift. (For the neutrino wave functions appearing in the decay mode the limit is taken, in which case the wave functions become the spherical Bessel functions.)

The radial wave functions are evaluated by means of the subroutine package RADIAL sal95 (), which implements a robust solution method that avoids the accumulation of truncation errors. This is done by solving the radial equations by using a piecewise exact power series expansion of the radial functions, which then are summed up to the prescribed accuracy so that truncation errors can be completely avoided. The input in the package is the potential . This potential is primarily the Coulomb potential of the daughter nucleus with charge , . As in the case of single- decay single () we include nuclear size corrections and screening.

The nuclear size corrections are taken into account by an uniform charge distribution in a sphere of radius with fm, i.e.


The introduction of finite nuclear size has also the advantage that the singularity at the origin in the solution of the Dirac equation is removed. (Other charge distributions, for example a Woods-Saxon distribution, can be used if needed.)

The contribution of screening to the phase space factors was extensively investigated in single- decay wil70 (); buh65 (). The screening potential is of order and thus gives a contribution of order relative to the pure Coulomb potential . We take it into account by using the Thomas-Fermi approximation. The Thomas-Fermi function , solution of the Thomas-Fermi equation


with and


where is the Bohr radius, is obtained by solving Eq. (15) for a point charge with boundary conditions


This takes into account the fact that the final atom is a positive ion with charge . With the introduction of this function, the potential including screening becomes


This can be rewritten in terms of an effective charge where now depends on . In order to solve Eq. (15), we use the Majorana method described in esp02 () which is valid both for a neutral atom and a positive ion. The method requires only one quadrature and is thus amenable to a simple solution. It is particularly useful here, since we want to evaluate screening corrections in several nuclei. The Thomas-Fermi electron density is approximate, especially at the origin. However, the screening correction is only of order relative to the Coulomb potential and the error on this small correction is therefore negligible. (A better method would be to do an atomic Hartree-Fock calculation and then fit the result to the expansion


where as in Eq. (15). However, it has been shown in single- decay that this method gives results comparable to the Thomas-Fermi approximation buh65 (), except in very light nuclei, , which we do not discuss here.) We also do not consider radiative corrections to the phase space factors which are of order and thus negligible to the order of approximation we consider in this article.

In order to show the improvement in our calculation as compared with the approximate solution used in the literature we show in Fig. 2 a comparison of the radial wave functions for Nd decay, , at MeV.

Figure 2: Electron radial wave functions , (left panel) and , (right panel) for , MeV and fm (vertical line). The notations WF1, WF2, and WF3 correspond to leading finite size Coulomb, exact finite size Coulomb and exact finite size Coulomb with electron screening, respectively.

Iii Phase space factors in double- decay

iii.1 Two neutrino double- decay

The decay, Fig. 1a, is a second order process in the effective weak interaction. It can be calculated in a way analogous to single- decay. Neglecting the neutrino mass, considering only S-wave states and noting that with four leptons in the final state we can have angular momentum , and, , we see that both and decays can occur. We denote by the -value of the decay, by the excitation energy in the intermediate nucleus, and by the excitation energy with respect to the average of the initial and final ground states,


The situation is illustrated in Fig. 3.

Figure 3: Notation used in this article. The example is for Nd decay.


The differential rate for -decay is given by (primakoff (); konopinski (); doi81 (); doi83a (); tom91 (); haxton ())


where and are the electron energies, and the neutrino energies, the angle between the two emitted electrons, and


The quantities and are a sum of the contributions of all the intermediate states and depend on the energy of the intermediate state in the odd-odd nucleus and on the nuclear matrix elements . Introducing the short-hand notation


where is a suitably chosen excitation energy in the odd-odd nucleus, one can write tom91 (), to a good approximation,


where are the nuclear matrix elements and and are products of radial wave functions. Since Eq. (24) is an approximation to the exact expression, which is, however, of crucial importance for the separation of the decay probability into a phase space factor and a nuclear matrix element we have investigated the dependence of and on the energy . Since appears both in the denominator of Eq. (24) through and and in the numerator through , the dependence on cancels almost completely, as already remarked years ago by Tomoda tom91 (), and as it is shown by explicit calculation in the following paragraphs.

The functions and are defined as




The functions and are obtained from the electron wave functions. We have used several ways to obtain and following an approach similar to that used in single- decay. We write


In approximation (I) we use the weighing function in which case


that is the electron wave functions are evaluated at the nuclear radius . This is the simplest approximation and is commonly used in single- decay. We adopt it in this article. In approximation (II) we use the weighing function for and for (an uniform distribution of radius ). This is not a good approximation, since the inner states cannot decay due to Pauli blocking and the decay occurs at the surface of the nucleus. Nevertheless, it is sometimes used. It essentially amounts to an evaluation of and at a radius , as one can show by explicitly evaluating


The third and most accurate approximation (III) is that in which the weighing function is the square of the wave function, , of the nucleon undergoing the decay,


By using harmonic oscillator wave functions and assuming that only one orbital is involved, the integrals in Eq. (30) can be easily evaluated. The approximation (III) essentially amounts to an evaluation of and at a radius . For harmonic oscillator wave functions




one has


This approximation has the disadvantage that it must be done separately for each nucleus. Since in this paper we are seeking greater generality and do not wish to make a commitment to definite nucleonic orbitals, we make use of approximation (I). However, our computer program is written in such way as to allow the possibility of using Eq. (30) instead of Eq. (28). Also in Sect. IV we study in a specific case, Pd, where the transition is between and orbitals, the error we make by using Eq. (28) instead of Eq. (30).

All quantities of interest are obtained by integration of Eq. (21). In the approximation described above, all quantities are separated into a phase space factor (independent of nuclear matrix elements) and the nuclear matrix elements. The two phase space factors are


where is determined as . It has become customary to normalize these to the electron mass . Also since the axial vector coupling constant is renormalized in nuclei it is convenient to separate it from the phase space factors and define quantities


These quantities are then in units of y. From these, we obtain:
(i) The half-life


(ii) The differential decay rate


where .
(iii) The summed energy spectrum of the two electrons


(iv) The angular correlation between the two electrons


We can evaluate the phase space factors for any value . The dependence of on is shown in Fig. 4b for the specific case of Pd decay. We see that depends mildly on () except very close to threshold , where the dependence is . A similar situation occurs for . We have done a calculation of and in the list of nuclei shown in Table 1 with from Ref. haxton () or estimated by the systematics MeV, which approximately represents the energy of the giant Gamow-Teller resonance in the intermediate odd-odd nucleus. The obtained values are also shown in Fig. 5 where they are compared with previous calculations boe92 (). These values of are those estimated in the closure approximation and should be combined with the closure matrix elements


where and . Here is the closure energy for states in the odd-odd intermediate nucleus and it can be approximately taken as the energy of the isobaric analogue state.

Figure 4: Panel a) Skeleton of the Pd decay scheme. The ground state of the intermediate Ag nucleus is leading to the lowest possible value for to be  MeV.
Panel b) Behaviour of the phase phase factor as a function of . The value obtained using single state dominance hypothesis, MeV, is denoted by a red circle and the value obtained using MeV MeV is denoted by a blue square.

In recent years, it has been suggested, that in some nuclei, the lowest intermediate state dominates the decay. This is called the single state dominance hypothesis (SSD) aba84 (); gri92 (); civ98 (); sim01 (); dom05 (). This situation is likely to occur in Zr, Mo, Pd, and Cd, where protons occupy mostly the level and neutrons mostly the level, and in Te, where protons occupy mostly the level and neutrons mostly the level, which are spin-orbit partners of each other. In the SSD model the energy is that of the single state . We have done a calculation of and for the nuclei mentioned above in the SSD case. This is also shown in Table 1 in columns 3 and 5. In this case, and should be combined with the matrix elements


Finally, using our program, one can evaluate the sum


if the individual GT matrix elements are known from a calculation, and a similar sum for Fermi matrix elements. In this case, there is no separation between phase space factors and nuclear matrix elements.

Nucleus y y y y (MeV) (MeV) (MeV)
Ca 15550. -11930. 4.27226(404) 7.717
Ge 48.17 -26.97 2.039061(7) 9.411
Se 1596. -1075. 2.99512(201) 10.08
Zr 6816. 7825. -4831. -5477. 3.35037(289) 10.97 2.203
Mo 3308. 4134. -2263. -2762. 3.03440(17) 11.20 1.685
Pd 137.7 146.9 -79.56 -84.45. 2.01785(64) 11.75 1.893
Cd 2764. 3176. -1857. -2108. 2.81350(13) 12.06 1.875
Sn 553.0 -342.7 2.28697(153) 12.47
Te 0.2688 0.2727 -0.1047 -0.1061 0.86587(131) 12.53 1.685
Te 1529. -993.9 2.52697(23) 13.27
Xe 1433. -927.2 2.45783(37) 13.06
Nd 324.8 -195.5 1.92875(192) 13.63
Nd 36430. -26860. 3.37138(20) 13.72
Sm 9.591 -4.816 1.21503(125) 13.90
Gd 193.8 -114.2 1.72969(126) 14.17
Pt 15.36 -8.499 1.04717(311) 15.76
Th 11.31 -6.779 0.84215(246) 17.06
U 14.57 -9.543 1.14498(125) 17.28
Table 1: Phase space factors and obtained using screened exact finite size Coulomb wave functions. The -values are taken from experiment ( Ref. mou10 (), Ref. rah08 (), Ref. fin11 (), Ref. rah11 (), Ref. sci09 (), Ref. red07 (), Ref. kol10 ()) when available, or from tables of recommended values. is taken from Ref. haxton () or estimated by the systematics, MeV, where A without tilde denotes the mass number. Phase space factors and correspond to values obtained using the SSD model, in which case the used is listed in the last column.
Figure 5: Phase space factors in units y. The label ”approximate” refers to the results obtained by the use of approximate electron wave functions. The figure is in semilogarithmic scale.

We also have available upon request for all nuclei in Table 1 the single electron spectra, summed energy spectra and angular correlations between the two outgoing electrons. As examples we show the cases of Xe Ba decay, Fig. 6, of very recent interest to EXO experiment exo () and the case of Se Kr, Fig. 7, of interest to NEMO experiment nemo (). The use of our ”exact” calculation makes a considerable difference as shown in Fig. 8. For the SSD case there is a difference in the single electron spectra at small energies , as is shown in Fig. 9 for Pd, and previously emphasized in Refs. sim01 (); dom05 ().

Figure 6: Single electron spectra (left panel), summed energy spectra (middle panel) and angular correlations between two outgoing electrons (right panel) for the Xe Ba -decay. The scale in the left and middle panels should be multiplied by when comparing with experiment.
Figure 7: Same as Fig. 6 for the Se Kr -decay.
Figure 8: Same as Fig. 6 for the Nd Sm -decay. The figure also shows the difference between our ”exact” calculation and the previously used approximate calculation.
Figure 9: Single electron spectra for the Pd Cd -decay obtained using the two approximations discussed in the text, namely closure approximation and single state dominance hypothesis.


The decay to the excited state, (Fig. 3), is also of interest. The phase space factor for this decay can be calculated using the formulas of the previous subsection, with replaced by


The results of this calculation are shown in Table 2.

Nucleus y y y y (MeV) (MeV)
Ca 0.3627 -0.1505 2.99722(16) 1.27504(253)
Ge 0.06978 -0.02380 1.122283(7) 0.916757(167)
Zr 175.4 185.3 -103.8 -109.2 1.14813(7) 2.20224(296)
Mo 60.55 65.18 -33.54 35.89 1.13032(10) 1.90408(27)
Pd 0.004842 0.004864 -0.001371 -0.001377 1.47312(12) 0.54773(76)
Cd 0.8727 0.8878 -0.3642 -0.3701 1.756864(24) 1.056636(154)
Sn 0.01988 -0.006408 1.657283(22) 0.629687(1552)
Te 0.07566 -0.02705 1.79352(11) 0.73345(34)
Xe 0.3622 -0.1451 1.578990(23) 0.878840(393)
Nd 0.009911 -0.003339 1.42446(4) 0.50429(196)
Nd 4329. -2934. 0.740382(22) 2.630998(222)
Sm 0.01850 -0.006583 0.6806673(18) 0.5343627(12518)
Gd 0.006318 -0.002178 1.279941(23) 0.449749(1283)
Th 0.00004221 -0.00001944 0.69142(9) 0.15073(255)
U 0.0004635 -0.0002289 0.94146(8) 0.20352(133)
Table 2: Phase space factors and for decay to the first excited states, , obtained using screened exact finite size Coulomb wave functions. Phase space factors and correspond to values obtained using the SSD model.


The half-life for -decay is given by equations similar to those of sect. III.1.1molina (); doi81 (); haxton (). The lepton phase space factor is now


with , (Fig. 3), from which the life-time can be calculated


The nuclear matrix elements can be written, in the closure approximation, as




Since this decay contains the term , it is suppressed, due to cancellations, and it will not be considered further. Also, other models (SSD, no-closure) can be used, if needed.

iii.2 Neutrinoless double- decay

The theory of decay was first formulated by Furry furry () and further developed by Primakoff and Rosen primakoff (), Molina and Pascual molina (), Doi et al. doi81 (), and, Haxton and Stephenson haxton (). Here we follow mainly the formulation of Tomoda tom91 (). The phase space factors for decay are simpler than those of because of the absence of integration over the neutrino energies. Also, with two leptons in the final state and S-wave decay we can only form angular momentum , and therefore the decay to is forbidden.


The differential rate for the decay is given by doi81 (); tom91 ()


where and are the electron energies, the angle between the two emitted electrons, and


This decay is forbidden by the standard model and can occur only if the neutrino has mass and/or there are right-handed currents. In view of recent experiments on neutrino oscillations SUP01 (); SNO2002 (); KAM2003 () it appears that neutrinos have a mass and we therefore consider the phase space factors for this case. The quantities and in Eq. (49) can then be written as tom91 ()


, where is the nuclear matrix element and , are the quantities given in Eq. (25).

All quantities of interest are then given by integration of Eq. (49). Introducing


where is determined as , and defining the quantities


where , fm, is the nuclear radius, we can calculate:
(i) The half-life


(ii) the single electron spectrum