Theoretical Advances in Beta and Double-Beta Decay

Abstract

Weak interaction processes continue to be hot topics in fundamental physics research. In this paper, we briefly review some recent advances in the theoretical study of beta and double-beta decays that include both the nuclear and atomic part of these processes. On the nuclear side, we present a statistical approach for the computation of the nuclear matrix elements (NME) for neutrinoless double-beta ($0\nu \beta \beta$). A range of NME values, the most probable value for NME, and the associated theoretical uncertainty are given. Correlations with other related observables are shown as well. On the atomic side, we first briefly review the methods used to obtain the electrons’ wave functions. Further, we use them for the computation of some relevant kinematic quantities such as Fermi functions, electron spectra, and angular correlation between the emitted electrons. Then, we present applications of these calculations to the experimental data analysis related to the search of the Lorentz invariance violation in two-neutrino double-beta ($2\nu \beta \beta$) decay and description of the decay rates and decay rate ratios for allowed and unique forbidden electron capture (EC) processes.

1. Introduction

The study of weak interaction transitions plays an essential role in physics research because they can provide critical information on various hot topics in atomic, nuclear, particle physics, and astrophysics. For example, the study of beta transitions is essential in describing the reactor processes [1] and the properties of the nuclei far from stability [2], as well as in understanding the evolution of the stars [3,4]. Also, the study of double-beta decay (DBD) offers a broad range of very current topics related to the nuclear structure, validity of some symmetry laws, yet-unknown properties of neutrinos, and constraining physics scenarios beyond the Standard Model (SM) [5,6,7,8,9,10].

Beta and double-beta decay rates can be written as a product of an atomic and a nuclear part. The atomic part, namely, the phase space factors (PSF), encompasses the perturbation of the movement of the emitted leptons (electrons, protons) due to their interaction with the Coulomb field and electron cloud of the daughter nucleus, while for the nuclear part, the nuclear matrix elements (NME) is related to the nuclear structure of the nuclei involved in the decay. For the ($0\nu \beta \beta$) decay channel, a third term appears as follows:

$$\begin{array}{c}\hfill {\left(T_{1/2}^{2\nu }\right)}^{-1}=G^{2\nu }(Q_{\beta \beta },Z)g_A^4\mid m_e{c}^2{M}^{2\nu }{\mid }^2\\ \hfill {\left(T_{1/2}^{0\nu }\right)}^{-1}={\mathrm{\Sigma }}_k{G}^{0\nu }(Q_{\beta \beta },Z)g_A^4\mid M_k^{0\nu }{\mid }^2{\left(<{\eta }_k>\right)}^2\end{array}$$

(1)

where $G^{2\nu ,0\nu }$ and $M^{2\nu ,0\nu }$ are PSF and NME for the $2\nu$ and $0\nu$ decay modes, $g_A$ is the axial vector coupling constant, and <${\eta }_k$> is a beyond SM parameter associated with different possible mechanisms that can trigger the $0\nu \beta \beta$ decay.

For beta transitions, the probability per unit time that a nucleus with atomic mass A and charge Z decays for a nth-forbidden $\beta$-branch is given by [11]:

$$\begin{array}{c}\hfill {\lambda }_n=g^2/2{\pi }^3{\int }_1^{W_0}pW{(W_0-W)}^2{S}_n(Z,W)dW,\end{array}$$

(2)

where g is the weak interaction coupling constant, p is the momentum of $\beta$-particle, W = $\sqrt{p^2+1}$ is the total energy of $\beta$-particle, and $W_0$ is the maximum $\beta$-particle energy. $W_0$ = $Q-1$, in $\beta$ decay (Q is the mass difference between initial and final states of neutral atoms). Equation (2) is written in natural units ( $\hslash =m=c=1$ ) so that the unit of momentum is $mc$, the unit of energy is $m{c}^2$, and the unit of time is ℏ /$m{c}^2$. The shape factors $S_n(Z,W)$, which appear in Equation (2), are defined as:

$$S_n(Z,W)=\sum _{k_e,k_v=1}^{\infty }\left[{\lambda }_k{M}_L^2+{\lambda }_k^{\prime }{m}_L^2-2{\lambda }_k^{\prime \prime }{M}_L{m}_L\right]$$

(3)

where $M_L$ and $m_L$ contain the nuclear matrix elements and ${\lambda }_k$, ${\lambda }^{\prime }$, and ${\lambda }^{\prime \prime }$ are bilinear combinations of radial components of the electron wave function.

As seen from the above formulas, in order to reliably predict beta and double-beta decay half-lives, interpret the measured electron spectra, extract information about neutrino properties, or constrain different BSM scenarios for $0\nu \beta \beta$ decay occurrence, both the nuclear and atomic parts require accurate computation. The two terms, PSF and NME, are calculated separately with atomic and nuclear methods. NME calculation is a hot topic in the description of the weak interaction processes that lasts for a long time and that is still not satisfactorily solved. For DBD, a second-order weak interaction process, this problem is crucial since NME appears in the half-life formulas at a power of two and any uncertainty in the computation is amplified in terms of half-lives prediction. For beta decay, the NME calculation is important as well for quantitative predictions of the decay rates, and it is not an easy task since the nuclei undergoing beta transitions belong to different nuclear mass regions and have various nuclear structures.

In this paper, we give a brief review of some recent advances in the computation of NME for DBD and kinematic quantities such as Fermi functions, electron spectra, and angular correlation between electrons relevant for data analysis. In Section 2, we present an analysis and new estimations of the NME for $0\nu \beta \beta$ decay. We comment on the present status of the NME computation and briefly present a statistical model aiming to study the stability of NME calculated with the interacting Shell model (ISM) method, at small, random variations of some model effective Hamiltonians. Also, we predict a range of the NME values and give the most likely value of it with a predicted theoretical error. In Section 3, we resume different approaches for obtaining the wave function (w.f.) of the emitted electrons, which are the basic ingredients in the computation of other relevant kinematic quantities for beta and double-beta transitions. In Section 4, we resume the theoretical framework developed for investigating the Lorentz invariance violation (LIV) in $2\nu \beta \beta$ decay, while in Section 5 we present a study of the electron capture (EC) processes using the Dirac–Hartree–Fock–Slater (DHFS) method. Finally, we end up with conclusions to this review.

2. Nuclear Matrix Elements

The computation of the NME for DBD is a long-standing problem that still needs to be satisfactorily solved. Until recently, their computation was performed with nuclear methods that use phenomenologically effective NN potentials. The most well known are: the interaction shell model (ISM) [12,13,14,15,16,17], proton–neutron quasi-random phase approximation (pnQRPA) methods [7,18,19,20,21,22], interacting boson model (IBM) [23], Energy Density Functional [24], PHFB [25], Coupled-Cluster method (CC) [26], in-medium generator coordinate method (IM-GCM) [27], and valence-space in-medium similarity renormalization group method (VS-IMSRG) [28]. The current situation is that there are still significant differences between NME values calculated with different methods [29]). These methods have their strengths or weaknesses, depending on the different approximations on which they are built. These are mainly related to the dimension of the model spaces for the single-particle basis and the number and type of the NN correlations involved in the calculation. For example, pnQRPA and IBM can use large model spaces for the single-particle states but only a few types of correlations between nucleons, while ISM can use only a restricted model space (only one shell) that includes all types of correlations and preserves the symmetries. As mentioned above, the common feature of these methods is that they all employ phenomenological effective NN potentials in the NME computation. Here, an important aspect is that while pnQRPA derives the NN effective potential from infinite matter using the Brueckner approach and adapts their form by multiplying with specific parameters, ISM uses effective NN potentials built for specific model spaces for different nuclear regions and checks their reliability by comparing the theoretical estimations with the data for other observables in the same nuclear mass region. However, as mentioned above, these methods still give differences between NME values, sometimes calculated with the same method (pnQRPA or ISM), which can amount to factors of three or more, depending on the specific isotope and values of the input parameters considered in the calculation. Because of these uncertainties, there is a need for more fundamental approaches for the NME computation for DBD. In this regard, there is currently a coherent program for advancing such calculations [30], with the main idea being to obtain NME values with minimal dependence on the computational model and provide a quantified theoretical uncertainty. It is hoped that this will result in progress in developing methods based on particle and nuclear effective field theories (EFT) [31], lattice quantum chromodynamics (QCD) [32], and ab-initio nuclear-many-body techniques [33]. However, the current NME values obtained with these methods are still far from good predictions. The available results are only for the NME of ${}^{48}\mathit{Ca}$ [26,34], ${}^{76}\mathit{Ge}$ [35], and ${}^{82}\mathit{Se}$ [35,36]. In this context, the traditional methods still predict $2\nu \beta \beta$ NME values close to the experimental measurements. One missing thing would be finding a reliable range of NME values with a predicted uncertainty to provide experimenters with a better plan to calibrate their set-ups and interpret the results.

In this context, we propose a statistical analysis of $0\nu \beta \beta$ NME and applied it to ${}^{48}\mathit{Ca}$ [37] and ${}^{136}\mathrm{Xe}$ [38] isotopes. We only consider the standard light LH neutrino exchange mass mechanism, which is most likely to contribute to the $0\nu \beta \beta$ decay process. Details of these models can be found in the cited references, and here we give only the relevant ingredients. We employed three independent effective Hamiltonians for each isotope, suitable for each nuclear mass region, namely, FPD6, GXPF1A, and KB3G for ${}^{48}\mathit{Ca}$, and SVD, GCN5082, and jj55t for ${}^{136}\mathit{Xe}$, and we investigated the effect of small, random variations of the shell model effective Hamiltonians on $0\nu \beta \beta$ NME. These effective Hamiltonians are described by a small number of single-particle energies and a finite number of two-body matrix elements. The wave functions produced by these Hamiltonian can also be used to describe and predict other observables, such as the electromagnetic and Gamow–Teller (GT) transition probabilities, nucleon occupation probabilities, and spectroscopic factors, using relative simple changes of the transition operators in terms of effective charges and quenching factors. These effective charges and quenching factors are calibrated to the existing data. For $0\nu \beta \beta$ NME, such calibrations are not yet possible due to the lack of data. However, different existing effective Hamiltonians for nuclei involved in a given $0\nu \beta \beta$ decay produce smaller ranges of the NME values.

The main goals of our study were: (a) for each starting effective Hamiltonian find correlations between $0\nu \beta \beta$ NME and the other observables that are accessible experimentally; (b) find theoretical ranges for each observable; (c) establish the shape of different distributions for each observable and starting Hamiltonians; (d) use this information to find weights of contributions from different starting Hamiltonians to the “optimal” distribution of the $0\nu \beta \beta$ NME; and (e) find an “optimal” value of the $0\nu \beta \beta$ NME and its predicted probable range (theoretical error). We refer further to the relevant results obtained for ${}^{136}\mathit{Xe}$, which is one of the most experimentally studied isotope.

For example, in

Table 1. Correlation matrix for 9 observables of ${}^{136}\mathrm{Xe}$ using the SVD effective Hamiltonian and comparison between experimental data and theoretical calculations for the same observables [38].

	${\mathit{M}}_{0\mathit{\nu }}$	${\mathit{M}}_{2\mathit{\nu }}$	D${\mathit{E}}_{2^{+}}$	D${\mathit{E}}_{6^{+}}$	D${\mathit{E}}_{4^{+}}$	D$\mathbf{GT}$	P$\mathbf{GT}$	P$\mathit{B}\left(\mathit{E}2\right)$	D$\mathit{B}\left(\mathit{E}2\right)$
$M_{0\nu }$	1	0.8	0.78	0.74	0.72	0.38	0.25	−0.18	−0.42
$M_{2\nu }$	0.8	1	0.47	0.43	0.42	0.52	0.44	−0.049	−0.24
D$E_{2^{+}}$	0.78	0.47	1	0.95	0.96	0.53	0.34	−0.54	−0.79
D$E_{6^{+}}$	0.74	0.43	0.95	1	0.99	0.5	0.32	−0.56	−0.84
D$E_{4^{+}}$	0.72	0.42	0.96	0.99	1	0.55	0.35	−0.55	−0.88
D$GT$	0.38	0.52	0.53	0.5	0.55	1	0.76	−0.23	−0.58
P$GT$	0.25	0.44	0.34	0.32	0.35	0.76	1	−0.25	−0.43
P$B\left(E2\right)$	−0.18	−0.049	−0.54	−0.56	−0.55	−0.23	−0.25	1	0.65
D$B\left(E2\right)$	−0.42	−0.24	−0.79	−0.84	−0.88	−0.58	−0.43	0.65	1
		[39]	[40]	[40]	[40]	[41]	[41]	[42]	[42]
Data	N/A	0.018	0.819	2.207	1.867	0.012	0.15	0.286	0.413
SVD	1.763	0.025	0.662	2.157	1.389	0.001	0.163	0.154	0.342
GCN	2.645	0.069	0.842	2.196	1.873	0.009	0.545	0.121	0.194
jj55t	2.314	0.06	0.917	2.502	2.113	0.004	0.512	0.096	0.158

, we present the correlations that we found between several relevant observables. The notations are obvious, and the letters “P” and “D” added to the name of the observables denote “parent” and “daughter” nuclei. The calculations were carried out with the ISM in the jj55 model space consisting of the 0g7/2, 1d5/2, 1d3/2, 2s1/2, and 0h11/2 orbitals that assume ${}^{100}\mathit{Sn}$ as a core, covering the sector of the nuclear chart between (N, Z) = 50–82. The correlation coefficients between observables shown in Table 1 were obtained using the Pearson method and the approach described in [37]. A set of 1000 Hamiltonians was generated from each starting Hamiltonian, after applying random perturbations on the two-body matrix elements of maximum amplitude $\pm 10\%$. This size of the maximum amplitude was chosen so as not to produce changes in the single particle energies on which the initial Hamiltonians were fitted. The methodology of calculating the $0\nu \beta \beta$ NME ($M^{0\nu }$) within ISM was extensively described elsewhere [13], and here we do not repeat the details of the calculation. First of all, as expected, $M^{0\nu }$ and $M^{2\nu }$ are strongly correlated, so an accurate computation of $M^{2\nu }$ may help in the computation of $M^{0\nu }$. Comparing the theoretical results with the experimental values in Table 1, one can see that the SVD starting Hamiltonian produces $M^{2\nu }$ NMEs that are more accurate, thus needing the least amount of quenching when compared to those of GCN5082 or jj55t. For the PGT and DGT, one observes that the SVD results are closest to the experimental data for the parent nucleus, overestimating the result by much less than GCN5082 and jj55t. However, for the daughter’s GT, GCN5082 was best, with SVD underestimating the result the most. For PB(E2)↑ and DB(E2)↑, SVD shows values closest to the experiment. The excitation energies are better described by GCN5082, in large part because the GCN5082 starting Hamiltonian was fine-tuned with data for more nuclei and energy levels than SVD and jj55t. Interestingly, the correlations between the $0\nu \beta \beta$ NME and the strengths of the parent and daughter GT transitions to the first $1^{+}$ state in ${}^{136}\mathit{Cs}$ are significantly reduced, while the correlation with the $2\nu \beta \beta$ NME is very strong. One explanation for this phenomenon could be related to the fact that the product of the GT matrix elements describing transitions to the first $1^{+}$ state in ${}^{136}\mathit{Cs}$ does not significantly contribute to the total sum of all excited $1^{+}$ states in the intermediate nucleus.

Figure 1 (taken from [38]) shows the probability distribution functions (PDF) for the three starting effective Hamiltonians and their weighted sum. The distribution weights were obtained using the Bayesian Model Averaging method, described in Section IV of [38]. Of the three Hamiltonians, only the SVD Hamiltonian accurately reproduces the experimental $2\nu \beta \beta$ decay NME with the same quenching factor typical for shell model GT calculations (q∼0.7). Thus, its contribution to the $0\nu \beta \beta$ decay NME is expected to have a larger weight in the probability distribution functions than the other Hamiltonians. Based on our statistical analysis summarized in this Figure, one can infer that with 90% confidence, the NME lies between 1.55 and 2.65, with a mean value of about 1.99 and a standard deviation of 0.37. We note that correlations between $0\nu \beta \beta$ NMEs and other observables were also studied in other works but in the context of finding useful experimental data for extracting information about $0\nu \beta \beta$ NME values. In [43], the $0\nu \beta \beta$ NMEs and double GT (DGT) transitions were calculated for ${}^{48}\mathit{Ca}$ with shell model using GXPF1B and KB3G effective Hamiltonians and then for heavier isotopes with the EDF and pnQRPA methods. In [44], the relation between double gamma ($\gamma \gamma$) decay NMEs and $0\nu \beta \beta$ NMEs was investigated. The calculations were performed within the shell model approach for a wide range of nuclei, using different effective Hamiltonians for the configuration spaces specific to different nuclear regions. In both of these works, the authors found a good correlation between the $0\nu \beta \beta$ NME values and the other observables (DGT and $\gamma \gamma$), revealing valuable information tools for $0\nu \beta \beta$ NMEs.

We mention that the present analysis can help ab-initio studies, such as [26,34,35], to better identify correlations and further reduce the uncertainties of the $0\nu \beta \beta$ NME.

3. Electron Wave Functions

To calculate the kinematic part of the beta and double-beta transitions, one needs to obtain reliable w.f. for the emitted electrons, which includes their interaction with the Coulomb-type field given by the protons and electronic cloud of the daughter nucleus. They are the essential ingredients to compute Fermi functions, PSF, electron spectra, decay rates, and all other needed kinematic quantities (experimentally measurable). Until now, several methods have been used to calculate the electron w.f. for continuous and bound states to describe the beta and double-beta transitions. Here, we cover some of them.

Approach A. In a non-relativistic treatment used in early calculations of the DBD decay rate [45,46], the electron w.f were taken as solutions of the Schrödinger equation in a Coulomb potential given by a point-like nucleus, and the Fermi function was built as a scattering solution for a point charge $Z_f$ to a plane wave, evaluated at the origin:

$$F^{NR}(Z_f,ϵ)=\frac{2\pi \eta }{1-e^{-2\pi \eta }}$$

(4)

where $\eta =\pm \alpha Z_fϵ/p$, “+” (“−”) are for electrons (positrons), $\alpha \simeq 1/137$ is the fine structure constant, $ϵ$ is the energy of the outgoing lepton, and $p=\left|\overrightarrow{p}\right|$ is its momentum.

This approximation fails badly for heavy nuclei, but it has the advantage that the DBD rate can be integrated analytically, resulting in a polynomial in powers of Q-values.

Approach B. In a relativistic treatment, the electron w.f. were obtained as solutions of the Dirac equation with a Coulomb potential given by a point charge, and the Fermi function can be written as [7]:

$$\begin{array}{c}\hfill F_0(Z_f,ϵ)=4{\left(2pR\right)}^{2(\gamma -1)}{e}^{\pi y}\frac{{\left|\Gamma \left({\gamma }_k+iy\right)\right|}^2}{{\left[\Gamma (2\gamma +1)\right]}^2},\end{array}$$

(5)

with

$$\gamma =\sqrt{1-{\left(\alpha Z_f\right)}^2},y=\pm \alpha Z_fϵ/p$$

(6)

with the same sign convention as above. R is the cut-off radius in the evaluation of the Dirac equation, which is taken to match the radius of the daughter nucleus (i.e., $R=1.2A^{1/3}$ fm). In another version of this relativistic treatment, the electron w.f. and Fermi functions were obtained by solving a Dirac equation with Coulomb potential given by a spherical finite-size charge (nucleus). The Fermi function can be also obtained in an analytical form by retaining the first term in a series expansion of the electron w.f. in spherical functions [5].

$$F_{k-1}(Z,\epsilon )={\left[\frac{\Gamma (2k+1)}{\Gamma \left(k\right)\Gamma \left(2{\gamma }_k+1\right)}\right]}^2{\left(2pR\right)}^{2\left({\gamma }_k-k\right)}{\left|\Gamma \left({\gamma }_k+iy\right)\right|}^2{e}^{\pi y}$$

(7)

These analytical forms of the Fermi functions were extensively used for computing DBD half-lives in many works [5,7]. We also note that the screening effect, which proved to be important, was not included in both of these methods.

Approach C. The Fermi functions are built from the exact electron w.f. [47], which are obtained by numerically solving the Dirac equation in a spherical Coulomb potential given by a finite-size charged nucleus and including the electron screening effect. The details of the calculations can be found in refs. [47,48,49].

$$\begin{gathered} \begin{array}{c}\hfill \left(\frac{d}{dr}+\frac{\kappa +1}{r}\right)g_{\kappa }(ϵ,r)=(ϵ+V\left(r\right)+1)f_{\kappa }(ϵ,r)\end{array} \\ \begin{array}{c}\hfill \left(\frac{d}{dr}+\frac{\kappa -1}{r}\right)f_{\kappa }(ϵ,r)=-(ϵ+V\left(r\right)-1)g_{\kappa }(ϵ,r)\end{array} \end{gathered}$$

(8)

where $R_A=1.2A^{1/3}$ fm is the nuclear radius. To compute the Fermi function and PSF, we have to obtain the electron factors $f=f_{11}^{\left(0\right)}$ [6]:

$$f_{11}^{\left(0\right)}=\mid f^{-1-1}{\mid }^2+\mid f_{11}{\mid }^2+\mid f_1^{-1}{\mid }^2+\mid f_1^{-1}{\mid }^2$$

(9)

with

$$\begin{array}{c}\hfill f^{-1-1}=g_{-1}\left({ϵ}_1\right)g_{-1}\left({ϵ}_2\right);f_{11}=f_1\left({ϵ}_1\right)f_1\left({ϵ}_2\right),\end{array}$$

(10)

$$\begin{array}{c}\hfill f_1^{-1}=g_{-1}\left({ϵ}_1\right)f_1\left({ϵ}_2\right);f_1^{-1}=f_1\left({ϵ}_1\right)g_1\left({ϵ}_2\right)\end{array}$$

(11)

which are related to the electron radials w.f., f, and g:

$$\begin{array}{c}\hfill g_{-1}\left(ϵ\right)=g_{-1}(ϵ,R_A);f_1\left(ϵ\right)=f_1(ϵ,R_A)\end{array}$$

(12)

for which g and f are the radial solutions of the above Dirac equation.

We note that in ref. [47], the Coulomb potential ($V\left(r\right)$) from the Dirac equation has the form given by a finite size charged nucleus, while in refs. [48,49] the Coulomb potential was built from a realistic distribution of the protons in the daughter nucleus. The difference between the Coulomb-type potentials can be seen in Figure 2 (taken from [48]). Notably, the differences in the PSF values for DBD with double-electron emission, computed with the two types of Coulomb potentials, amount at most to a few percent, while for other decay channels the differences may be relevant. We mention that in both of these methods, the screening effect was introduced as a renormalization of the Coulomb potential with a Thomas Fermi function [47,48,49]. To obtain the radial electron/positron w.f., the Dirac equation was solved using the subroutines package RADIAL [50].

There are differences when applying the approaches (A–C) to calculating the electron spectra. For example, in Figure 3 the summed energy electron spectra for the $2\nu \beta \beta$ decay of ${}^{136}\mathit{Xe}$ and its deviation due to LIV effects are presented, with methods A, B, and C as is specified in the caption. One can see a large difference between the non-relativistic and relativistic treatments and between the B and C approaches.

Approach D. In this approach, the electron w.f. are obtained with the DHFS method using a Coulomb potential given by a nucleus with a diffuse surface. More details can be found in [52]. This method also gives accurate calculations of the kinematic quantities that appear in the beta and double-beta decay studies and better applies to EC processes, as we will see in Section 5. The Dirac equation is solved iteratively: one starts with predefined w.f. and computes the potential, then the equation is solved to obtain the new w.f., and the potential is recalculated until a convergence criterion is reached. The DHFS method constructs the potential as a sum of three components: the nuclear potential generated by a nucleus with a realistic Fermi proton charge distribution, the electronic potential as a mean field generated by the charge distribution of the electronic cloud, and the exchange potential assuring the asymptotic condition at $r\to \infty$ as [50]:

$$V_{\mathrm{DHFS}}\left(r\right)=V_{\mathrm{nuc}}\left(r\right)+V_{\mathrm{el}}\left(r\right)+V_{\mathrm{ex}}\left(r\right),$$

(13)

In this approach, the screening effect is considered by including the mean-field treatment of the electronic potential component, $V_{\mathrm{el}}\left(r\right)$.

4. Lorentz Violation in the Neutrino Sector

Lorentz invariance violation is a current topic that joins the increasing effort to test the SM limits. The theoretical framework that introduces LIV contributions is the SM extension (SME) theory [53], including operators that break Lorentz invariance for all the particles in the SM. Direct observation of LIV implies an investigation of processes at Plank scale (∼10¹⁹GeV), which is not currently possible. However, LIV can manifest in processes at low energy scales by suppressed effects potentially observable with current or near future experimental techniques. In particular, the neutrino sector of SME provides the theoretical framework for a rich phenomenology for searching evidence of LIV [54]. The operators that couple to neutrinos in SME affect neutrino flavour oscillations, neutrino velocity, and the electron spectra of beta decay and double-beta decay. The first experimental searches of LIV in this sector were carried out in the neutrino oscillations experiments [55,56,57]. There are, however, distinct LIV signatures (so-called counter-shaded effects) associated with the oscillation-free operators of mass dimension three, which cannot be investigated in these experiments. An analysis of electron spectra in beta and double-beta decays offers the possibility to investigate the LIV effects related to the time-like (isotropic) component of this oscillation-free operator, whose size is controlled by the coefficient ${\stackrel{\circ }{a_{\mathrm{of}}}}^{\left(3\right)}$ In ref. [58], the LIV effects in $2\nu \beta \beta$ decay were calculated for the summed energy spectra of electrons employing a non-relativistic approximation for the electron radial w.f. Currently, the accuracy required in the DBD experiments far exceeds this approximation. Therefore, precise computation of the electron spectra is needed. So far several experiments have reported limits for the ${\stackrel{\circ }{a_{\mathrm{of}}}}^{\left(3\right)}$ parameter from the analysis of the summed energy spectra of electrons in $2\nu \beta \beta$ decays, using theoretical predictions of these spectra performed with approximate methods: EXO-200 [59] $-2.65\times {10}^{-2}<$ ${\stackrel{\circ }{a_{\mathrm{of}}}}^{\left(3\right)}<7.6\times {10}^{-3}$ MeV, CUPID-0 [60] ${\stackrel{\circ }{a_{\mathrm{of}}}}^{\left(3\right)}<4.1\times {10}^{-3}$ MeV, NEMO-3 [61] $(-4.2<$ ${\stackrel{\circ }{a_{\mathrm{of}}}}^{\left(3\right)}<3.5)\times {10}^{-4}$ MeV, GERDA [62] $(-2.7<$ ${\stackrel{\circ }{a_{\mathrm{of}}}}^{\left(3\right)}<6.2)\times {10}^{-3}$ MeV. In [51], we examined the LIV effects on summed energy spectra of electrons and quantities related to them using Fermi functions built with exact electron w.f. (Approach C is described in Section 3). First, we derived the complete formalism for LIV analyses in $2\nu \beta \beta$ decay. For the sake of consistency, here we give the final formulas, while the detailed derivation of them can be found in refs. [51,63,64].

The DBD decay rate can be written as a sum of the SM decay rate and its deviation due to LIV [58]:

$${\Gamma }_{\mathrm{SME}}={\Gamma }_{\mathrm{SM}}+\delta \Gamma$$

(14)

As was already mentioned earlier, the $2\nu \beta \beta$ decay rate can be written as a product of NME and PSF.

$${\Gamma }_{SM}=g_A^4{\left|m_e{c}^2{M}^{2\nu }\right|}^2{G}^{2\nu }$$

(15)

Since the nuclear part does not contain neutrinos, and the LIV effects appear through the change of the neutrino momentum, only the phase space factors $G^{2\nu }$ are affected, where the neutrino energy and momentum appear explicitly. The angular differential DBD rate within the SM can be expressed as [5,63]:

$$\frac{d{\Gamma }_{\mathrm{SM}}}{d\left(cos{\theta }_{12}\right)}=\frac{1}{2}\left({\Gamma }_{\mathrm{SM}}+{\Lambda }_{\mathrm{SM}}cos{\theta }_{12}\right)=\frac{1}{2}{\Gamma }_{\mathrm{SM}}\left(1+{\kappa }_{\mathrm{SM}}cos{\theta }_{12}\right).$$

(16)

A similar formula can be written within SME [63]:

$$\frac{d{\Gamma }_{\mathrm{SME}}}{d\left(cos{\theta }_{12}\right)}=C{G}_{\mathrm{SM}}\times \left[1+{\stackrel{\circ }{a_{\mathrm{of}}}}^{\left(3\right)}\frac{\delta G}{G_{\mathrm{SM}}}\right.\left.+\left({\kappa }_{\mathrm{SM}}+{\stackrel{\circ }{a_{\mathrm{of}}}}^{\left(3\right)}\frac{\delta H}{G_{\mathrm{SM}}}\right)cos{\theta }_{12}\right]$$

(17)

The spectrum part ( ${\Gamma }^{2\nu }$) and the angular correlation part (${\Lambda }^{2\nu }$) of the decay rate can be obtained by integrating Equations (17) and (18) over the lepton energies. C is a constant, G and H are the PSF associated with the two decay rate parts, ${\theta }_{12}$ is the angle between the two emitted electrons, $\kappa$ is the angular correlation coefficient defined as the ratio of $\Gamma$ and $\Lambda$ decay rates, and $a_{off}^{0\left(3\right)}$ is the coefficient that governs the size of the LIV effects. Besides the summed energy electron spectra, one can also define another quantity related to the summed energy electron spectra that can be compared with experimental measurements, ${\chi }^{+}\left(K\right)$:

$$\frac{d{\Gamma }_{\mathrm{SME}}}{dK}=C\frac{d{G}_{\mathrm{SM}}}{dK}\left(1+{\stackrel{\circ }{a_{\mathrm{of}}}}^{\left(3\right)}{\chi }^{(+)}\left(K\right)\right)$$

(18)

$${\chi }^{(+)}\left(K\right)=\frac{d\left(\delta G\right)}{dK}/\frac{d{G}_{\mathrm{SM}}}{dK}$$

(19)

The PSF expressions and their deviations due to LIV can be gathered as follows [63]:

$$\begin{array}{cc}\hfill \left\{\begin{array}{c}{G}_{\mathrm{SM}}\\ \delta G\end{array}\right\}=& \frac{{\tilde{A}}^2{G}_F^2{\left|V_{\mathrm{ud}}\right|}^2{m}_e^9}{96{\pi }^{7}ln2}\frac{1}{m_e^{11}}{\int }_{m_e}^{E_I-E_F-m_e}d{\epsilon }_1{\epsilon }_1{p}_1{\int }_{m_e}^{E_I-E_F-{\epsilon }_1}d{\epsilon }_2{\epsilon }_2{p}_2\hfill \\ \hfill & \times {\int }_0^{E_I-E_F-{\epsilon }_1-{\epsilon }_2}d{\omega }_1{\omega }_2^{2}a\left({\epsilon }_1,{\epsilon }_2\right)\left[{〈K_N〉}^2+{〈L_N〉}^2+〈K_N〉〈L_N〉\right]\left\{\begin{array}{c}{\omega }_1^2\\ 4{\stackrel{\circ }{a_{\mathrm{of}}}}^{\left(3\right)}{\omega }_1\end{array}\right\},\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill \left\{\begin{array}{c}{H}_{\mathrm{SM}}\\ \delta H\end{array}\right\}=& \frac{{\tilde{A}}^2{G}_F^2{\left|V_{\mathrm{ud}}\right|}^2{m}_e^9}{96{\pi }^{7}ln2}\frac{1}{m_e^{11}}{\int }_{m_e}^{E_I-E_F-m_e}d{\epsilon }_1{\epsilon }_1{p}_1{\int }_{m_e}^{E_I-E_F-{\epsilon }_1}d{\epsilon }_2{\epsilon }_2{p}_2\hfill \\ \hfill & \times {\int }_0^{E_I-E_F-{\epsilon }_1-{\epsilon }_2}d{\omega }_1{\omega }_2^{2}b\left({\epsilon }_1,{\epsilon }_2\right)\left[\frac{2}{3}{〈K_N〉}^2+\frac{2}{3}{〈L_N〉}^2+\frac{5}{3}〈K_N〉〈L_N〉\right]\left\{\begin{array}{c}{\omega }_1^2\\ 4{\stackrel{\circ }{a_{\mathrm{of}}}}^{\left(3\right)}{\omega }_1\end{array}\right\}\hfill \end{array}$$

(21)

Until now, the LIV investigations have been carried out in several DBD experiments and are limited to extracting the LIV coefficient from the analysis of the summed electron spectra. The deviation of the summed energy electron spectrum due to LIV effects was presented in Section 3, Figure 3 (taken from [51]), and it is displayed as a shift of the maximum of the spectrum to higher energy. This deviation was computed using improved electron w.f. (Method C) and up-dated Q-values [51]. In our next works [63,64], we computed the LIV deviations of the summed energy electron spectra for all DBD isotopes of experimental interest. Then, we also extended our study to the single electron spectra (Figure 4) and electron angular correlation (Figure 5b), thus providing additional possible LIV investigations by analyzing these spectra as well.

Single electron spectra can also undergo LIV deviations, but they can only be investigated in experiments using detectors with tracking systems of individual electrons. Also, we suggested the analysis of the quantities ${\chi }^{+}$ (Equation (19)) in the vicinity of the Q-values, where larger LIV deviations of the SM spectra would be expected, as seen in Figure 5a. By comparing the electron and angular correlation spectra computed with and without LIV contributions, one can obtain information about the observability of the LIV effects with the current experimental statistics. Further, we proposed a new, alternative method to constrain the LIV parameter. This consists in comparing the measured value of the angular correlation coefficient k, which can be determined in the DBD experiments with electron tracking systems [65] by using the forward–backwards asymmetry, with its theoretical value. Equation (22) shows how the $k_{SME}$ coefficient can be determined by measuring the balance between the number of electrons emitted forward and backward.

$$A\equiv \frac{{\int }_{-1}^0\frac{d\Gamma }{dx}dx-{\int }_0^1\frac{d\Gamma }{dx}dx}{\Gamma }=\frac{N_{+}-N_{-}}{N_{+}+N_{-}}=\frac{1}{2}{k}_{\mathrm{SME}}$$

(22)

Finally, we mention that the search for LIV effects will have a greater chance of success in future DBD experiments as the statistics for the DBD events, especially in the vicinity of the Q-values, will improve. Also, there are more chances to be observed in DBD experiments with tracking detectors for individual electrons. We hope our works have strenghtened the theoretical support necessary for further research into LIV effects in DBD.

5. Decay Rates for EC Processes

Electron capture processes manifest across a broad spectrum of isotopes [66] and play a crucial role in practical applications such as radionuclide metrology [67] and nuclear medicine [68]. In these applications, the relevance is that most Auger electrons emitted during EC processes possess kinetic energies in the range of a few keV, enabling precise deposition within a confined area, thus allowing targeted irradiation of tumour sites. EC processes also play a determinant role in fundamental physics research like nuclear astrophysics [69] and Dark Matter and DBD experiments for accurate background characterization [70]. Notably, in the precise measurement of the $2\nu ECEC$ in ${}^{124}\mathit{Xe}$, a critical background contribution emerges from ${}^{125}\mathit{I}$ EC as its decay peak closely aligns with that of the $2\nu ECEC$ peak [71]. That is why there is a great interest in understanding these processes better, and in this regard, theoretical support is crucial.

In this context, we develop a program to calculate EC processes for a wide nuclear region, A = 7–204 [72]. Here, we give only the main ingredients of our work, highlighting the improvements and the obtained results. The decay rate ($\lambda$) for allowed and unique-forbidden transitions is defined as [73]:

$$\lambda =M_L^2\frac{(2L-2)!!}{(2L-1)!!}\sum _{n\kappa }\frac{n_{n\kappa }{p}_{n\kappa }^{2\left(\left|\kappa \right|-1\right)}{q}_{n\kappa }^{2\left(L-\left|\kappa \right|+1\right)}{\beta }_{n\kappa }^2{B}_{n\kappa }{S}_{n\kappa }}{\left(2\left|\kappa \right|-1\right)!\left[2\left(L-\left|\kappa \right|\right)+1\right]!},$$

(23)

where ${\beta }_{n\kappa }$ is the Coulomb amplitude, $q_{n\kappa }$ represents the energy of the emitted neutrino, $B_{n\kappa }$ denotes the overlap and exchange corrections, and $S_{n\kappa }$ is the shake-up and shake-off correction. The angular momentum of the electron capture transition is represented as L, and $M_L$ contains the nuclear matrix element. The total energy of the captured electron is given by $W_{n\kappa }=m_e-\left|E_{n\kappa }\right|$, with $E_{n\kappa }$ representing the binding energy and $m_e$ representing the electron rest-mass energy. The momentum of the bound electron is obtained as $p_{n\kappa }=\sqrt{m_e^2-W_{n\kappa }^2}$, and the momentum of the emitted neutrino is denoted $q_{n\kappa }$. The relative occupation number is obtained as $n_{n\kappa }=N_{n\kappa }/2\left|\kappa \right|$.

The decay energy Q is calculated as the atomic mass difference between the initial and final atom:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill Q& =q_{n\kappa }+R_{\gamma }+\left[B_{\mathrm{gs}}\left(Z^{\prime }\right)-B_{n\kappa }\left(Z^{\prime }\right)\right]\hfill \\ \hfill & =q_{n\kappa }+R_{\gamma }+R_{n\kappa },\hfill \end{array}\end{array}$$

(24)

where $R_{\gamma }$ and $R_{n\kappa }$ are the $\gamma$ and X-ray radiation energies released by the final nucleus/ atom. The total binding energy of the electronic cloud is represented as $B_{n\kappa }$.

More details can be found in [72]. The most important elements of the calculation are the electron w.f., the single-electron binding energies, and the total electron binding energies. We utilize the DHFS self-consistent framework for their computation. The framework is validated by comparing the single-electron energies and the total binding energies we obtained with the experimental ones, thus observing good agreement. By comparison, other models artificially induce the convergence of energies to specific values for each atomic number. A mention-worthy approach is the BETASHAPE (BS) program [74]. One significant improvement of our work is a refined energy balance in Equation (24). We consider the whole system (the initial and final atom), not only the nuclei and the captured electron. The binding energies of electrons that do not participate directly in the decay (spectator electrons) change between the initial and final state. This has an important effect on the energy balance, affecting the neutrino energy, which enters directly in the decay rate Equation (23). This approach yielded improved agreement with experimental values for EC fractions, particularly in low-energy transitions.

For our calculations, we computed the w.f. of all electrons in the initial atom in the ground state configuration and in the final atom, considering all possible configurations after an electron is captured. Then, we computed the shake-up and shake-off correction and the exchange and overlap correction using the w.f.s of the electrons.

In

Table 2. Comparison between several theoretical models results and measured electron capture decay ratios for isotopes studied in [76]. The RD marked with ¹ are taken from [76].

Isotope	$\mathit{Q}-{\mathit{R}}_{\mathit{\gamma }}$ [76]	Quantity	BS [76]	KLI [76]	[72]	RD [75]
	(keV)			Frozen Orbitals
${}^7\mathrm{Be}$	861.89 (7)	${\lambda }_L/{\lambda }_K$	0.105 (8)	0.0509 (20)	0.11053 (3)	0.101 (13)
${}^{41}\mathrm{Ca}$	421.64 (14)	${\lambda }_L/{\lambda }_K$	0.09800 (40)	0.09078 (16)	0.1046 (2)	0.102 (10)
${}^{54}\mathrm{Mn}$	542.2 (10)	${\lambda }_L/{\lambda }_K$	0.11219 (31)	0.09590 (19)	0.1076 (6)	0.1066 (16)
		${\lambda }_K/\lambda$	0.88419 (34)	0.90005 (21)	0.8870 (5)	0.8896 (17)
${}^{55}\mathrm{Fe}$	231.21 (18)	${\lambda }_L/{\lambda }_K$	0.11629 (31)	0.10073 (20)	0.1121 (3)	0.1110 (15)
		${\lambda }_M/{\lambda }_K$	0.01824 (12)	0.014824 (45)	0.01909 (5)	0.01786 (29) 1
		${\lambda }_M/{\lambda }_L$	0.1568 (11)	0.14716 (49)	0.1704 (4)	0.1556 (26) 1
${}^{109}\mathrm{Cd}$	127.1 (18)	${\lambda }_K/\lambda$	0.8148 (14)	0.8164 (12)	0.810 (7)	0.812 (3)
		${\lambda }_{L+}/{\lambda }_K$	0.2274 (12)	0.2250 (12)	0.2344 (101)	0.2315 (8)
${}^{125}\mathrm{I}$	150.28 (6)	${\lambda }_K/\lambda$	0.79927 (41)	0.80376 (23)	0.7983 (18)	0.8011 (17)
${}^{138}\mathrm{La}$	312.6 (3)	${\lambda }_L/{\lambda }_K$	0.3913 (25)	0.4242 (49)	0.409 (7)	0.432 (6)
		${\lambda }_M/{\lambda }_K$	0.0965 (9)	0.1002 (11)	0.100 (2)	0.102 (3) 1
		${\lambda }_M/{\lambda }_L$	0.2465 (20)	0.2362 (24)	0.244 (1)	0.261 (9) 1

(taken from [72]), we present a comparison of the decay constant ratios for multiple shells. We employ the use of ratios in order to reduce the NME, thus obtaining measurable quantities independent of the nuclear model. We outline a few important cases compared to other experimental (RD) [75] and theoretic results. The experimental values marked with ¹ are taken from [76]. We also compare our results with other results from [76] computed using a self-interaction-corrected model, the Krieger–Li–Iafrate (KLI) [77]. In Table 2, we present five allowed transitions ${}^7\mathrm{Be}$, ${}^{54}\mathrm{Mn}$, ${}^{55}\mathrm{Fe}$, ${}^{109}\mathrm{Cd}$, and ${}^{125}\mathrm{I}$; a first unique forbidden transition ${}^{41}\mathrm{Ca}$; and a second unique transition ${}^{138}\mathrm{La}$.

The uncertainties in Table 2 were obtained using the pseudo-experiment technique. We took into account the experimental uncertainties of the Q-values and of the energy levels of the final nuclei in case of the ground state to excited state transitions, and we assumed a 3% error on the theoretical atomic relaxation energies. In the cases presented here, the largest source of uncertainty comes from the Q-values and energy levels of the final nuclei. In low Q-value cases, the uncertainty from the atomic relaxation energies becomes dominant.

We noted deviations from experimental values below 2% for the ${\lambda }_K/\lambda$ ratios across all models, increasing to 12% for captures from higher shells. The BS model and our model consistently provided the most accurate values. Additionally, our model’s validity was tested by comparing theoretical predictions to experimental values across a wide range of atomic numbers and transition energies, revealing excellent agreement (below two standard deviations) in most cases.

Lastly, we emphasize that refined energetics contribute to improved agreement between experimental and theoretical EC ratios, particularly in low-energy transitions. This characteristic has potential implications for determining the neutrino mass scale from electron capture processes. Our findings have significant implications for future studies in the EC field and related nuclear physics and astrophysics applications.

6. Conclusions

We briefly reviewed some recent advances in the theoretical study of beta and double-beta decays. A statistical approach for determining the NME for $0\nu \beta \beta$ is presented on the nuclear side. We studied the stability of the NME values against small, random variations of some effective Hamiltonians used in the computation and their correlation with other observables that are accessible experimentally. From this statistical study, we predicted a theoretical range of the NME values and determined an optimal NME value together with a predictable error. On the atomic side, we briefly reviewed the methods for obtaining the electron w.f. Then, we used them to compute relevant, measurable quantities for the data analysis in beta and double-beta decay experiments. First, we deduced the theoretical formulas for electron spectra deviations and angular correlation deviations in $2\nu \beta \beta$ decay due to LIV. We proposed the search of possible LIV effects in the summed electron energy and single electron spectra and angular correlation, as well as a new, alternative method to constrain the coefficient ${a_{\mathrm{of}}}^{\left(3\right)}$ that governs the size of the LIV effects. Then, we presented a formalism based on the DHFS approach to compute decay rates and decay ratios from atomic multi-shells, for allowed and unique forbidden EC processes. Our calculations display good agreement with experimental data. We hope that the results reviewed in this paper represent relevant advances that can improve data analyses in beta and double-beta decay experiments.