An effective method to accurately calculate the phase space factors for $\beta^- \beta^-$ decay

Accurate calculations of the electron phase space factors are necessary for reliable predictions of double-beta decay rates, and for the analysis of the associated electron angular and energy distributions. We present an effective method to calculate these phase space factors that takes into account the distorted Coulomb field of the daughter nucleus, yet allows one to easily calculate the phase space factors with good accuracy relative to the most exact methods available in the recent literature.


I. INTRODUCTION
Double-beta decay (ββ) processes are of considerable importance for the study of neutrinos. They change the charge Z of a nucleus by two units, releasing two electrons, while the mass A remains unchanged. The ββ decay with two associated electron antineutrinos in the final state conserves the lepton number, and is permitted within the Standard Model (SM). This process, called two-neutrino double-beta decay (2νββ), has been experimentally observed for several isotopes with transitions to both ground states and to excited states of the daughter nuclei [1]. Should the lepton number conservation be violated, then theories beyond the standard model (BSM) predict that the ββ decay transition could occur without antineutrinos in the final state, called neutrinoless double-beta (0νββ), and this implies that the neutrino is a Majorana fermion [2]. The 0νββ transitions have not yet been confirmed experimentally, but there are many recent experimental and theoretical efforts dedicated to their discovery. Recent reviews on this matter are in Refs. [3][4][5]. There are several mechanisms that could contribute to the 0νββ decay rate, of which the simplest and most studied one involves the exchange of light Majorana neutrinos in the presence of left-handed weak interaction. Other, more complex, mechanisms include contributions from right-handed currents [6,7], and mechanisms involving super-symmetry [5,8].
The phase space factors (PSF) that enter the ββ lifetimes expressions were considered for a long time as being accurately calculated (see e.g. Refs. [9,10]). Recent reevaluations of the PSF, using methods that take into account the proton distributions distorting the Coulomb field of the daughter nucleus [11][12][13][14], have shown considerable differences in some cases when compared to the previous results [9,10]. A very recent paper [14] presents four of the different methods commonly used to calculate the PSF, and compares their results for the case of 0νββ ground state (g.s.) transitions. Table 1 and Fig. 2 of Ref. [14] show that the Coulomb distortion of the electron wavefunction by inclusion of the finite nuclear size and electron screening effects can produce differences of up to 100%, compared to the point-charge formalism of Ref. [9] (see for example the 0νββ PSF G 08 of 150 Nd in Ref. [14]). However, taking into account the charge distributions in the daughter nuclei and solving numerically the Dirac equation with finite nuclear size is very slow and plagued by convergence issues. This makes these complex methods unattractive for the calculations of electron angular and energy distributions, such as those presented in Ref. [15,16].
In this paper, we propose an effective method for treating the distortion of the Coulomb field in the daughter nucleus. This method uses the well known formalism of Ref. [9], but provides accurate results that are in good agreement with those of Refs. [11][12][13][14]. This method could be particularly useful when performing complex investigations involving PSF to test BSM physics due to different possible underlying mechanisms contributing to the 0νββ process. These investigations often involve calculations of electron distributions [15,16], where components of the PSF enter the equations, and it is not possible to only use the tabulated values of Refs. [11][12][13][14].
The paper is organized as follows: Section II shows the formalism for 0νββ transitions to ground states, and for 2νββ transitions to ground and excited states. In Section III we present our effective method for the treatment of the distorted Coulomb field in the daughter nucleus. Section IV is dedicated to the results, and Section V shows our conclusions. The Appendixes summarize the pointcharge formalism from Refs. [9,10] that we adjusted to calculate the 0νββ and 2νββ PSF.

II. BRIEF FORMALISM OF THE ββ DECAY
For the 0νββ decay, one usually writes the inverse halflife as products of electron PSF, nuclear matrix elements (NME) that depends on the nuclear structure of the parent and that of the daughter nuclei, and lepton number violation (LNV) parameters of the BSM mechanisms taken into account. Considering the existence of righthanded currents, one would find several additional con-tributions to the decay rate [7,9]. The most studied mechanism is that of the light left-handed neutrino exchange, but other mechanisms could be of importance [5]. One popular model that includes contributions of right-handed currents is the left-right symmetric model [17,18]. This model assumes the existence of heavy particles that are not included in the Standard Model (SM). Within this framework, the 0νββ half-life expression is given by where g A is the Axial-Vector coupling strength, η ν , η L NR , η R NR , η λ , and η η are neutrino physics parameters defined in Ref. [19], M 0ν and M 0N are the light and heavy neutrino-exchange nuclear matrix elements [5,20], and X λ and X η represent combinations of NME and phase space factors. G 0ν 01 is a phase space factor [10] that can be calculated with relatively good precision in most cases [11,12,14]. Other possible contributions, such as those of R-parity violating SUSY particle exchange [5,20], etc, are neglected here. With some simplifying notations the half-life expression [9] (here we omit the contribution from the η L NR and η R NR terms, which share the same PSF as η 2 ν term, G 2ν 01 , and have the same energy and angular distribution as the η ν term) is written as where φ 1 and φ 2 are the relative CP-violating phases [19], and M 0ν GT is the Gamow-Teller part of the light left-handed neutrino-exchange NME. Different processes give rise to several contributions: C 1 comes from the lefthanded leptonic and currents, C 4 from the right-handed leptonic and right-handed hadronic currents, and C 5 from the right-handed leptonic and left-handed hadronic currents. The interference between these terms is represented by the C 2 , C 3 and C 6 contributions. Neglecting the very small tensor contributions in the mass mechanism, the C 1−6 components are defined as products of PSF and NME [9]: with The fractions of NME are defined [9] as and P indicating other NME. All these nine NME were calculated by several methods, including the interacting shell model (ISM) [16,21,22] and quasiparticle random phase approximation (QRPA) [23]. The light-neutrino mass mechanism M 0ν GT and M 0ν F NME have been extensively studied with many nuclear structure methods, such as interacting boson model (IBM-2) [24][25][26][27], interacting shell model (ISM) [20,21,[28][29][30][31][32][33][34][35][36][37], quasiparticle random phase approximation (QRPA) [38][39][40][41][42], projected hartree fock bogoliubov (PHFB) [43], energy density functional (EDF) [44], and the relativistic energy density functional (REDF) [45] method. The NME calculated with different methods and by different groups still present large differences, and that has been a topic of many debates in the literature (see e.g. [46,47]). Expressions for the G 0ν 01 − G 0ν 09 PSF are given in A. For the 2νββ process, the half-life for the transition to a state of angular momentum J (J = 0 or 2) of the daughter nucleus is given to a good approximation by [48] T 2ν where G 2ν (J) is a phase space factor [9, 10, 13] described in B, m e is the electron mass, and M 2ν (J) is the 2νββ NME, which can be calculated as [10,48] M 2ν Here the k-sum is taken over the 1 + k states with excitation energies E k in the intermediate odd-odd nucleus.
is the Q-value for the transition to the state of angular momentum J in the daughter nucleus, and ∆M is the difference in mass between the intermediate nucleus and the decaying nucleus.

III. DESCRIPTION OF THE EFFECTIVE "SCREENING" METHOD
Our approach is based of the formalism from Ref. [9,10] where the nuclear charge is considered point-like, but we replicate the effects of a finite size proton distribution distorting Coulomb field in the daughter nucleus by modifying the charge of the final nucleus (Z f ).
We multiply the Z f with a parameter, called "screening factor" (S f ) in what follows, to obtain an effective "screened charge" (Z s = S f 100 Z f ). For large enough energies, the tail of the Coulomb field plays a less significant role when compared to its part close to the nucleus, and the effect resembles a charge screening. The PSF calculated with Z s for each nucleus are compared to those of Refs. [11][12][13][14] (called "data" below), which were obtained with methods that consider Dirac electron wave functions calculated with finite nuclear size and atomic electron screening. Refs. [11,14] take into account the finite nuclear size by an uniform charge distribution of radius R, while Refs. [12,13] consider a more realistic Woods-Saxon proton distribution inside the nucleus. It was shown [11] that the atomic electron screening effect is small, of the order of 0.1%. The relative deviations between our results and the data, expressed in percentages (∆ = 100 PSF -data data ), are denoted with ∆G 0ν 01−09 for the 9 PSF of 0νββ transitions to ground states, ∆G 2ν g.s. for PSF of 2νββ transitions to ground states, ∆G 2ν , and those for which p-wave electrons contribute G 0ν 02 − G 0ν 09 . Here, we treat them separately naming them s-PSF and p-PSF, respectively. We consider the largest deviation (∆ max ) between the PSF of a certain class and the corresponding data, and we search for the value of S f that minimizes it. Our goal is to maintain ∆ max ≤ 10%. This value of the maximum deviation is considerably lower than the uncertainties of the NME contributing to the decay rate, Eq. (2). We find that controlling the maximum deviation provides more stable and predictable results than minimizing a χ 2 distribution.
In our analysis data is selected as follows: G 0ν 01−09 PSF are chosen from Table III of Ref. [14]. Other recent results for G 0ν 01 [11][12][13] are within a few percent of these values. For take into account the finite size effects of the charge distribution. The G 2ν g.s. data is taken from Table 1 of the very recent Ref. [13], and it is in very good agreement with the results of Ref. [11]. For the 2νββ transitions to the first excited 0 + states we take the data from Table 2 of Ref. [13]. There are four cases ( 110 Pd, 124 Sn, 130 Te, and 136 Xe) of PSF in Ref. [13] that are in significant disagreement with those of Ref. [11]. We do not take them into account in our analysis. The data for 2νββ transitions to the first excited 2 + states is taken from Table 3 of Ref. [13]. In this case, there are three PSF values ( 116 Cd, 124 Sn, and 136 Xe) that seems to deviate significantly from the model results. These 2νββ PSF were not confirmed by other groups, and they were often readjusted [49]. We do not include them in the analysis, but we compare them with our prediction in Table III. The 2νββ data of 124 Sn attributed to Ref. [13] in Tables  I-III was provided to us as private communications by the authors of Ref. [13].

IV. RESULTS AND DISCUSSIONS
For the analysis of the s−PSF we consider the G 0ν 01 of Ref. [14] and G 2ν g.s. , G 2ν of Ref. [13]. We find that the smallest maximum deviations from the data can be obtained using an optimal "screening factor" S f = 94.5. Fig. 1 shows how the maximum deviation reaches a minimum when one gets close to the optimal "screening factor". Table I  and ∆G 2ν g.s. for transitions to ground states. The adjusted charge of the daughter nucleus is also presented, together with the Q ββ values. We find very good agreement for these s−PSF and the data, with deviations smaller than 5%. Should one consider the point-charge formalism [9], the largest deviation goes up to 40% for the case of G 0ν 01 of 150 Nd (see, e.g., Table 1 columns A and D of Ref. [14]).  formalism deviations exceed 38%. The PSF marked with the "*" symbol correspond to the four nuclei not included in our analysis. Our results for these cases can be considered as predictions for these cases that are not yet validated by other methods. The last two columns show the results of Ref. [11] and the corresponding deviations, for comparison. Ref. [11] provides no value for 82 Se. The case of the 124 Sn is more complicated because of the different values used in the literature for the energy of the first excited 0 + state (see Ref. [50] for details). We include here with * the phase space factor corresponding to the Q ββ used in Ref. [11], and with (*) the phase space factor corresponding to the Q ββ considered in Ref. [51] (see discussion in Ref. [50]). The G 2ν 2 + 1 PSF and their deviations are displayed in Table III. We find the largest deviation ∆G 2ν 2 + 1 = 6.3% for 48 Ca. Neglecting finite nuclear size effects, one would get a deviation of 47% for 150 Nd. Similar to the previous table, the three results excluded from the analysis are presented for comparison and marked with the "*" symbol.
When calculating the p−PSF, G 0ν 02 − G 0ν 09 we find a different optimal "screening factor", S f = 92, corresponding to a larger maximum deviation, ∆ max = 18.1%. Fig.  2 presents the evolution of ∆ max close to the "optimum screening factor" for p−PSF. We attribute this larger deviations to the different kinematic factors of the nine 0νββ PSF (see Eqs. (A5)).
To further minimize the deviations, we obtain eight optimal "screening factors" corresponding to G 0ν 02 − G 0ν 09 , as seen in Fig. 3. The best results are presented in Table  IV, where we show the optimal "screening factor", the maximum deviations, and the PSF values. The last line presents the optimal "screening factor" that minimizes the deviations of all the p−PSF. Alongside G 0ν 02 − G 0ν 09 , we display the G 0ν 01 obtained with the optimal "screening factor" common for all s−PSF.

V. CONCLUSIONS
In this paper we present an effective method to calculate the phase space factors of the β − β − transitions, which can provide results close to those of methods that consider the finite size of the proton charge and the atomic electron screening. It modifies the point-charge formalism of Refs. [9,10], by considering a constant multiplicative screening factor for the charge of the daughter nucleus. The main advantage of our method consists in its simplicity given its accuracy, and its potential to be extended to calculations of the energy and angular electron distributions needed for the analysis of the contributions of the right-handed currents to the 0νββ decay.
Our method works well for PSF of 0νββ and 2νββ transitions to ground stated, and also for 2νββ transitions to the first excited 0 + and 2 + states. For PSF where only s−wave electrons contribute, an effective "screening factor", S f = 94.5, was obtained. Using this S f value one finds a maximum deviation of 6.3% between our results and other results in recent literature [11][12][13][14]. In the  . 3: The behaviour of maximum deviations from the reference values of the 9 PSF involved in the 0νββ transitions to ground states. Subfigures (a) to (i) correspond to G 0ν 01 to G 0ν 09 . The horizontal axis represents the "screening factor" values, and the vertical axis represents the maximum deviations expressed in percentages.
case of the PSF where p−wave electrons contribute, we obtained another optimal "screening factor", S f = 92. This corresponds to a maximum deviation of 18.1% between our results and those of Ref. [14]. We attribute this large deviation to the kinematic factors of Eqs. (A5). The deviations are greatly reduced, to less than 10%, when considering individual "screening factors" for each specific PSF (G 0ν 02 −G 0ν 09 ). It is remarkable that in the case of the G 08 , the original point-charge formalism [9] PSF deviates by over 100% for 150 Nd, while it is significantly reduced in our model. Similar spectacular reductions are found for other p − P SF . We also provide predictions for the PSF of some isotopes, which can be also used as guidance in cases of disagreement between the more precise methods.
In addition, using S f = 92 one gets the largest max- imum deviation of 18.1 for all neutrinoless double-beta decay PSF, G 0ν 01 − G 0ν 09 . This information is relevant for the calculation of the two-electron energy and angular distributions [16].
We conclude that this method is well suited for fast and accurate calculations of the ββ decay PSF, with uncertainties much lower than those of the associated NME. One could envision to further reduce these PSF uncertainties by considering a mass-dependent screening factor. A Mathematica notebook that can be used to obtain all these phase space factors can be downloaded from [52].  Table I. The last line shows the optimal "screening factor" S f for all 8 p−PSF (G 0ν 02 − G 0ν 09 ). For G 01 , the s−PSF optimal "screening factor", S f = 94.5 of Table I, was used. Shown in the last column are the maximum deviations between our calculations with the indicated parameters and the results from Ref. [14]. 48