Phase space factors for double-\beta decay

# Phase space factors for double-β decay

## Abstract

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.

###### pacs:
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

 (A,Z)→(A,Z±2)+2e∓+anything. (1)

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

For , the process , Fig. 1a,

 (A,Z)→(A,Z+2)+2e−+2¯ν (2)

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

 (A,Z)→(A,Z+2)+2e− (3)

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,

 (A,Z)→(A,Z+2)+2e−+M0 (4)

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

 (A,Z)→(A,Z−2)+2e++2ν(A,Z)→(A,Z−2)+2e+. (5)

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

 [τ2ν1/2]−1=G2ν|M2ν|2, (6)

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

 [τ0ν1/2]−1=G0ν|M0ν|2|f(mi,Uei)|2, (7)

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,

 ψϵκμ(r)=(gκ(ϵ,r)χμκifκ(ϵ,r)χμ−κ,), (8)

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:

 dgκ(ϵ,r)dr=κrgκ(ϵ,r)+ϵ−V+mec2cℏfκ(ϵ,r),dfκ(ϵ,r)dr=−ϵ−V−mec2cℏgκ(ϵ,r)+κrfκ(ϵ,r). (9)

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

 es(ϵ,r)=eS1/2s(ϵ,r)+eP1/2s(ϵ,r)+eP3/2s(ϵ,r)+... (10)

where

 (11)

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

 (gκ(ϵ,r)fκ(ϵ,r))∼e−iδκℏpr⎛⎜ ⎜⎝√ϵ+mec22ϵsin(kr−lπ2−ηln(2kr)+δκ)√ϵ−mec22ϵcos(kr−lπ2−ηln(2kr)+δκ)⎞⎟ ⎟⎠, (12)

where

 k≡pℏ=√ϵ2+(mec2)2cℏ (13)

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.

 V(r)=⎡⎢ ⎢⎣−Zd(αℏc)r,r≥R−Zd(αℏc)(3−(r/R)22R),r

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

 d2φdx2=φ3/2√x (15)

with and

 b=12(3π4)2/3ℏ2mee2Z−1/3d≃0.8853a0Z−1/3d, (16)

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

 φ(0)=1,φ(∞)=2Zd. (17)

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

 V(r)≡φ(r)×⎡⎢ ⎢⎣−Zd(αℏc)r,r≥R−Zd(αℏc)(3−(r/R)22R),r

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

 V(r)=(−Zd(αℏc)/r)∑iaiexp(−bix), (19)

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.

## 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,

 ~A=12W0+EN−EI=12(Qββ+2mec2)+EN−EI. (20)

The situation is illustrated in Fig. 3.

#### 0+→0+12νββ-decay

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

 dW2ν=(a(0)+a(1)cosθ12)w2νdω1dϵ1dϵ2d(cosθ12) (21)

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

 w2ν=g4A(GcosθC)464π7ℏω21ω22(p1c)(p2c)ϵ1ϵ2. (22)

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

 (23)

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

 a(0)=14f(0)11|M2ν|2~A2[(⟨KN⟩+⟨LN⟩)2+13(⟨KN⟩−⟨LN⟩)2],a(1)=14f(1)11|M2ν|2~A2[(⟨KN⟩+⟨LN⟩)2−19(⟨KN⟩−⟨LN⟩)2], (24)

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

 f(0)11=|f−1−1|2+|f11|2+|f−11|2+|f1−1|2,f(1)11=−2Re[f−1−1f∗11+f−11f1−1∗]. (25)

with

 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)g−1(ϵ2). (26)

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

 g−1(ϵ)=∫∞0w(r)g−1(ϵ,r)r2dr,f1(ϵ)=∫∞0w(r)f1(ϵ,r)r2dr. (27)

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

 g−1(ϵ)=g−1(ϵ,R)f1(ϵ)=f1(ϵ,R),(I) (28)

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

 g−1(ϵ)=3R3∫R0g−1(ϵ,r)r2drf1(ϵ)=3R3∫R0f1(ϵ,r)r2dr.(II) (29)

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,

 g−1(ϵ)=∫∞0|Rnl(r)|2g−1(ϵ,r)r2drf1(ϵ)=∫∞0|Rnl(r)|2f1(ϵ,r)r2dr.(III) (30)

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

 Rnl(r)=√2n!b3Γ(n+l+3/2)(rb)le−r2/2b2Ll+1/2n(r2/b2) (31)

with

 b2=ℏMω≃1.0A1/3fm2, (32)

one has

 ⟨r2⟩nl=b2(2n+l+32). (33)

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

 F(0)2ν=2~A23ln2∫Qββ+mec2mec2∫Qββ+mec2−ϵ1mec2∫Qββ−ϵ1−ϵ20f(0)11×(⟨KN⟩2+⟨LN⟩2+⟨KN⟩⟨LN⟩)w2νdω1dϵ2dϵ1, (34)
 F(1)2ν=2~A29ln2∫Qββ+mec2mec2∫Qββ+mec2−ϵ1mec2∫Qββ−ϵ1−ϵ20f(1)11×[2(⟨KN⟩2+⟨LN⟩2)+5⟨KN⟩⟨LN⟩]w2νdω1dϵ2dϵ1, (35)

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

 G(i)2ν=F(i)2νg4A(mec2)2. (36)

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

 [τ2ν1/2]−1=G(0)2νg4A∣∣mec2M2ν∣∣2. (37)

(ii) The differential decay rate

 dW2νdϵ1=N2νdG(0)2νdϵ1, (38)

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

 dW2νd(ϵ1+ϵ2−2mec2)=N2νdG(0)2νd(ϵ1+ϵ2−2mec2). (39)

(iv) The angular correlation between the two electrons

 α(ϵ1)=dG(1)2ν/dϵ1dG(0)2ν/dϵ1. (40)

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

 M2ν≃(gVgA)2MF2ν~AF−MGT2ν~AGT, (41)

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.

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

 MGT2ν=⟨0+F||τ+→σ||1+1⟩⟨1+1||τ+→σ||0+I⟩12(Qββ+2mec2)+E1+1−EI. (42)

Finally, using our program, one can evaluate the sum

 ∑NG(i)2ν,N⟨0+F||τ+→σ||1+N⟩⟨1+N||τ+→σ||0+I⟩12(Qββ+2mec2)+EN−EI (43)

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.

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 ().

#### 0+→0+22νββ-decay

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

 Qββ−Ex(0+2)=Qββ(0+2) (44)

The results of this calculation are shown in Table 2.

#### 0+→2+12νββ-decay

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

 F(0)0+→2+12ν=2~A6ln2∫Qββ(2+1)+mec2mec2∫Qββ(2+1)+mec2−ϵ1mec2∫Qββ(2+1)−ϵ1−ϵ20f(0)11×(⟨KN⟩−⟨LN⟩)2w2νdω1dϵ2dϵ1, (45)

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

 [τ2ν1/2(0+→2+)]−1=F(0)0+→2+2ν∣∣M(2+)2ν∣∣2. (46)

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

 M(2+)2ν≃−MGT(2+)2ν~A3 (47)

where

 MGT(2+)2ν=⟨2+F||∑nn′τnτn′[→σn⊗→σn′](2)||0+I⟩. (48)

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.

#### 0+→0+10νββ-decay

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

 dW0ν=(a(0)+a(1)cosθ12)w0νdϵ1d(cosθ12) (49)

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

 w0ν=g4A(GcosθC)416π5(mec2)2(ℏc2)(p1c)(p2c)ϵ1ϵ2 (50)

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 ()

 (51)

, 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

 F(i)0ν=2ln2∫Qββ+mec2mec2f(i)11w0νdϵ1, (52)

where is determined as , and defining the quantities

 G(i)0ν=F(i)0νg4A(4R2) (53)

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

 (54)

(ii) the single electron spectrum

 dW0νdϵ1=N0νdG(0)<