Double Beta Decay: A Shell Model Approach

: Studies of weak interaction in nuclei are important tools for testing different aspects of the fundamental symmetries of the Standard Model. Neutrinoless double beta decay offers an unique venue of investigating the possibility that neutrinos are Majorana fermions and that the lepton number conservation law is violated. Here, I use a shell model approach to calculate the nuclear matrix elements needed to extract the lepton-number-violating parameters of a few nuclei of experimental interest from the latest experimental lower limits of neutrinoless double beta decay half-lives. The analysis presented here could reveal valuable information regarding the dominant neutrinoless double beta decay mechanism if experimental half-life data become available for different isotopes. A complementary shell model analysis of the two-neutrino double beta decay nuclear matrix elements and half-lives is also presented.


Introduction
The recent experimental discovery of neutrino oscillations [1,2] proved that neutrinos have mass, and this discovery was awarded a Nobel prize in 2015 [3,4]. Neutrino oscillation experiments can only provide information about the squared mass differences, while other properties of neutrinos, such as their mass hierarchy, their absolute masses, or their fermionic signatures, Dirac or Majorana, remain to be determined. However, this new information coming from the neutrino oscillations experiments has led to new interest in neutrino physics and in particular in their nature as Dirac or Majorana fermions that may be unraveled by neutrinoless double beta decay investigations.
Neutrinoless double beta decay (0νββ) is one of the best experimental approaches for identifying processes that violate the lepton number conservation, thus signaling beyond the Standard Model (BSM) physics. If neutrinoless double beta transitions occur, then the lepton number conservation is violated by two units, and the black-box theorems [5][6][7][8] indicate that the light left-handed neutrinos are Majorana fermions. As a consequence, the BSM extension of the Standard Model Lagrangian would be significantly different from that where neutrinos are Dirac fermions. Theoretical investigations of 0νββ decay combine lepton number violation (LNV) amplitudes with leptonic phase-space factors (PSFs) and nuclear matrix elements (NMEs). The NMEs are computed using a large variety of nuclear structure methods and specific models. Among the LNV models considered, the left-right symmetric model [9][10][11][12][13] is among the most popular, and its predictions are currently investigated at the Large Hadron Collider [14]. In some recent papers [15][16][17], I have investigated observables that could identify the contributions of different left-right symmetric model mechanisms to the 0νββ decay rate, such as the angular distribution and the energy distribution of the two outgoing electrons that could be measured. A more general approach is effective field theory (EFT), which considers an expansion of the BSM Lagrangian consistent with the Standard Model symmetries and including LNV and neutrino mass mechanisms. This approach has the advantage of being independent of specific models, and it can be used to describe in a unified manner BSM-sensitive observables, including those related to 0νββ decay. One can then use the existing data/limits from different experiments to evaluate the energy scales up to which the effective field operators are not broken and limits for effective low-energy couplings.
The theoretical analysis of the 0νββ decay process has many steps, including the nuclear structure calculation of the NMEs. However, in the first step, the weak interaction of quarks and leptons described by the BSM EFT Lagrangian is considered in the lowest order (see the diagram in Figure 1). In the next step, the hadronization process to nucleons and exchanging pions is considered, leading to the diagram in Figure 2. Furthermore, the nucleons are treated in the impulse approximation leading to free space 0νββ transition operators, and the nucleon dynamics inside the nuclei are treated using nonperturbative nuclear wave functions, which are later used to obtain the nuclear matrix elements needed to calculate the 0νββ observables, such as half-lives and two-electron angular and energy distributions [15]. A modern approach that can be used to make the transition from quarks and gluons to nucleons and pions is based on the chiral effective field theory of pions and nucleons [18,19]. This approach introduces a number of effective low-energy couplings, which in principle can be calculated from the underlying theory of strong interaction using lattice QCD techniques [18] or may be extracted within some approximation from the known experimental data [19]. These couplings may have new complex phases, and they could include effective contributions from the exchange of heavier mesons. The lattice QCD approach is in progress (see, e.g., Ref. [20]), but it has proven to be difficult for extracting some of the necessary weak nucleon coupllowestngs, even the known g A [20].  Figure 1, the nucleon-level diagrams of 0νββ decay process: (a) the typical 0νββ decay process nucleon-level diagram presents the generic description of the process and (b) the light left-handed neutrino exchange diagram shows the light left-handed neutrino exchange. Here, ". . ." stands for other higher-order effective field theory (EFT) diagrams (see Figure 2 of Ref. [17]).
Here, as in Ref. [17], I use the formalism of Refs. [21][22][23][24] that provides a general EFT approach to the BSM Lagrangian. It also provides a somewhat older hadronization scheme, which is needed to obtain the neutrinoless double beta decay transition operators. To extract new limits for the effective Majorana mass and for the low-energy EFT couplings from the current experiment for the isotopes listed in Table 1 below, I use the assumption that only one single coupling in the BSM Lagrangian may dominate the 0νββ amplitude.
In the analysis, about 20 nuclear matrix elements and nine phase-space factors are needed. I use the existing neutrinoless double beta decay data to extract the limits for the BSM EFT couplings and limits of validity for the energy scale of the BSM Lagrangian. In addition, the calculated ratio of half-lives for different isotopes could be useful in guiding the experimental effort, in estimating their scales and costs, in fine-tuning the experimental searches for the 0νββ transition mechanism, and also in providing a better view and comparison of the status of various experimental efforts. Our analysis suggests that the experimental confirmation of 0νββ decay rates for several isotopes could possibly help in identifying the dominant mechanism responsible for the transition.  [25,26] and T 0ν 1/2 limits (in years), and the calculated PSFs G 2ν [27] and G 01 (G 02 -G 09 can be found elsewhere [17] [32] 1100 [33] One important step in describing the 0νββ decay observables is obtaining the appropriate NMEs. The nuclear structure methods used for NME calculations are the interacting shell model [34][35][36][37][38][39][40][41][42][43][44][45][46][47][48][49][50][51][52], proton-neutron quasi random phase approximation (pnQRPA) [21][22][23][24][53][54][55][56][57], interacting moson model [58][59][60][61], projected Hartree-Fock-Bogoliubov [62], energy density functional [63], and relativistic energy density functional method [64]. The NMEs calculated with different methods and by different groups show sometimes large variations by a factor of 3-5 [65,66]. Most references only provide NMEs for the light left-handed Majorana neutrino exchange mechanism, but some provide results for the right/left heavy neutrino exchange and some more exotic mechanisms. Ref. [50] provides tables and plots that compare results for the light left-handed neutrino exchange and for the heavy righthanded neutrino exchange, while Ref. [17] provides tables with all NMEs necessary for the EFT approach. I calculate the NMEs using shell model techniques [36,[41][42][43][44][45][46][47][48][49][50][51] and a preferred set of effective Hamiltonians that were tested for a wide set of nuclei. The shell model calculations of NMEs use a relatively small single-particle model space, but they are better suited and more reliable for 0νββ decay calculations because they take into account all the correlations around the Fermi surface, respect all nuclear many-body problem symmetries, and can take into account the effects of the missing single particle space via many-body perturbation theory (the effects were shown to be small [67]). In addition, it was shown [68,69] that the QRPA approaches using the same model spaces and effective Hamiltonian as in the shell model produce NMEs within 25% of the shell model results. Furthermore, I test the shell model methods and the effective Hamiltonians by comparing the calculations of spectroscopic observables for the nuclei involved in the transition to the experimental data, as presented in Refs. [41,50,70]. I do not consider any quenching for the bare 0νββ operator in these calculations. Such a choice is different from that for the simple Gamow-Teller operator used in the single beta and two-neutrino double beta decay (2νββ), where a quenching factor of about 0.7 is necessary [69]. For the PSFs, I use an effective theory based on the formalism of Ref. [71], but fine-tuned as to take into account the effects of a Coulomb-field-distorting finite-size proton distribution in the daughter nuclei. Table 1 provides information relevant for the main nuclei that can be calculated using shell model techniques (see Equations (1) and (13) below for a precise definition of the PSFs used).
In this paper, I mostly review the shell model techniques needed to accomplish the plan outline above. The numerical results and their analysis are available in different papers that are appropriately cited below. Although most material described below reviews results already published, some new results can be found at the end of Section 3.2 and in Section 4. The paper is organized as follows: Section 2 analyzes the contributions of several BSM mechanisms to neutrinoless double beta decay, and it presents the framework of effective field theory for neutrinoless double beta decay; Section 3 presents an analysis of the 0νββ nuclear matrix elements in the shell model approach; Section 4 presents an analysis of the 2νββ nuclear matrix elements in the shell model approach; Section 5 is dedicated to conclusions.

Neutrinoless Double Beta Decay And Neutrino Physics
The main mechanism considered to be responsible for neutrinoless double beta decay is the mass mechanism that assumes that the neutrinos are Majorana fermions and relies on the assumption that the light left-handed neutrinos have mass. However, the possibility that right-handed currents could contribute to neutrinoless double beta decay (0νββ) has been already considered for some time [71,72]. Recently, 0νββ studies [13,73] have adopted the left-right symmetric model [11,74] for the inclusion of right-handed currents at quark level. In addition, the R-parity-violating ( R p ) supersymmetric (SUSY) model can also contribute to the neutrinoless double beta decay process [75][76][77].

LNV Models Contributing to 0νββ
In the framework that includes the left-right symmetric model and R-parity-violating SUSY model, after hadronization, the 0νββ half-life can be written as a sum of products of PSFs, BSM LNV parameters, and their corresponding NMEs [15]: Here, G 01 is a phase-space factor that can be calculated with good precision for most cases [27,28,78,79], g A is the axial vector coupling constant, η 0ν = m ββ /m e , effective Majorana neutrino mass (see Equation (3)), and m e is the electron mass. η L N R and η R N R are the heavy neutrino parameters with left-handed and right-handed currents, respectively [13,36], ηq and η λ are R p SUSY LNV parameters [80], and η λ and η η are parameters for the so-called "λ-" and "η-mechanisms", respectively [13]. M 0ν and M 0N are the light and the heavy neutrino exchange NMEs, Mq and M λ are the R p SUSY NMEs, and X λ and X η denote the combinations of NMEs and other PSFs (G 02 -G 09 ) corresponding to the the λ-mechanism involving right-handed leptonic and right-handed hadronic currents and the η-mechanism with right-handed leptonic and left-handed hadronic currents, respectively [15]. Assuming a seesaw type I dominance [81], the term η L N R is considered negligible if the heavy mass eigenstates are larger than 1 GeV [52], and I ignore it here. For consistency with the literature, the remaining term η R N R is labeled as η 0N . Here, I exclusively describe transitions from the spin/parity J π = 0 + ground state (g.s.) of the parent nucleus to the final J π = 0 + ground state of the daughter nucleus. There is also the possibility of 0νββ decay to the excited states of the daughter, such as the first J π = 2 + . This alternative is rarely considered in the literature, mainly because besides a significant reduction in the effective Q-values for most isotopes, thus reducing the corresponding phase space factors, it has also been known for some time that based on a general analysis the NMEs for this transition are suppressed for the mass mechanism [72]. In addition, the initial numerical estimates of the NMEs corresponding to the η η and η λ in Equation (1) showed that they were also suppressed [82]. Recently, it was found that more up-to-date QRPA calculations of these right currents' contributions could lead to a significant increase in the matrix elements for the lambda mechanism that might compete with the transition to the J π = 0 + ground state, at least for case of 136 Xe [83,84]. These new findings are clearly interesting, and I plan to investigate them using shell model techniques similar to the ones described below and report them in future publications. Table 1 presents the Q ββ values, the most recent experimental half-life limits, and the nine PSFs for the 0νββ transitions to the ground states of the daughter nucleus for five isotopes considered in this investigation. The PSFs were calculated using a new effective method described in detail in Ref. [27]. G 01 values were calculated with a screening factor (s f ) of 94.5, while for G 02 -G 09 I used s f = 92.0, which was shown to provide results close to those of the more accurate approach described in Ref. [85].
As indicated in Equation (1), the main observable related to 0νββ decay is the half-life of the process. It is unlikely that this unique observable, even if measured for several isotopes, could provide enough information to identify different mechanisms that may contribute to this process. In Ref. [15], I investigated other observables that could be used to disentangle contributions from different mechanisms, such as the two-electron angular and energy distributions, in addition to the half-life data from several isotopes. I considered the case where one mechanism dominates, i.e., there is one single term in the decay amplitude of Equation (1). Table 2 of Ref. [17] shows the shell model values of the NMEs that enter Equation (1). Details regarding the definitions of specific NMEs can be found in Refs. [17,49]. All NMEs were calculated using the interacting shell model (ISM) approach [36,[43][44][45][46]49,52] (see also Ref. [49] for a review) and included short-range correlation effects based on the CD-Bonn parametrization [41], finite-size effects [80], and, when appropriate, optimal closure energies [70] (see Section 3 for more details). Table 2 of Ref. [17] also presents the upper limits for the corresponding LNV parameters extracted from the lower limits of the half-lives under the assumption of one-mechanism dominance. However, less general analyses are available based on QRPA [71,80,[85][86][87], NMEs, and other interactive shell model NMEs [34][35][36][37].
If only the main diagram in Figure 2b is considered, the associated mechanism is known as the light neutrino exchange mechanism and the half-life of Equation (1) becomes with the effective neutrino mass given by following sum over the light mass eigenstates: where U ei are the complex matrix elements of the first row in the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) neutrino mixing matrix. This quantity is very often used in the literature as an example of how one could potentially extract additional information about neutrino physics parameters, such as neutrino mass ordering and the mass of the lowest mass eigenstate, from the experimental value of T 0ν 1/2 [88].

EFT Approach to 0νββ Decay
As mentioned in the introduction, a more general approach could be constructed based on the effective field theory extension of the Standard Model. Such an EFT analysis is preferable because it does not rely on specific models and the parameters could be constrained by the existing 0νββ data and by data from the Large Hadron Collider and other experiments. In addition, the models considered in Equation (1) always lead to a subset of terms in the low-energy (∼200 MeV) effective field theory Lagrangian. EFT considers all terms in the BSM Lagrangian allowed by the symmetries, some of them corresponding to the model terms incorporated in Equation (1), but the couplings might have a wider meaning. Other terms in the EFT Lagrangian are new, not directly identifiable with those originating from specific models.
At the quark level, Figure 1 shows the generic 0νββ Feynman diagrams contributing to the 0νββ process. I consider contributions coming from the light left-handed Majorana neutrino (Figure 1b) and a long-range part coming from the low-energy fourfermion charged-current interaction (see Ref. [17] for details). After hadronization (see Figure 2), the extra terms in the Lagrangian require the knowledge of about 20 individual NMEs [22][23][24]75,80,89]. One can write the half-life in a factorized compact form: Here, the E i contain the neutrino physics parameters, E 1 = η 0ν represent the exchange of light left-handed neutrinos, , ε 4 , ε 5 , η 1π , η 2π } denote the short-range LNV parameters at the quark level (see Ref. [17] for definitions of notations and details). The contributions of pion-exchange diagrams are also included in the so-called "higher-order term in nucleon currents" [80]. However, they are constrained by partial conservation of axial current (PCAC) and are only included in the light neutrino exchange contribution in Figure 2a. This contribution changes the associated NMEs by only 20%, and one concludes that it does not represent a serious double counting issue.
In Equation (4), M 2 i and M ij are combinations of NMEs and integrated PSFs [27] denoted with G 01 -G 09 (see Ref. [17] for definitions and details). In some cases, the interference terms E i E j M ij are small [90] and can be neglected, but not all of them [91]. In Ref. [15], I analyzed a subset of terms contributing to the half-life formula, with Equation (1) originating from the left-right symmetric model. In that restrictive case, I showed that one can disentangle different contributions to the 0νββ decay process using two-electron angular and energy distributions as well as the half-lives of two selected isotopes.

Neutrinoless Double Beta Decay Nuclear Matrix Elements
From previous Sections, one can conclude that the analysis of main experimental data regarding the 0νββ, the half-lives of multiple isotopes, and the two-electron angular and energy distributions [15] require a set of nuclear matrix elements. In this Section, I describe different techniques for calculating NMEs, starting with the direct summation on the states in the intermediate nucleus (Z − 1, N − 1), where Z denotes the atomic number and N the number of neutrons in a nucleus, and continuing with the often used closure approximation. An alternative method that performs a summation on the intermediate states in the (Z − 2, N) or (Z, N − 2) nuclei is described in Ref. [46].

The Anatomy of the 0νββ NMEs
The nuclear matrix elements needed in Equations (1)-(4) describe the transition from an initial nucleus |i = |0 + i to a final nucleus | f = |0 + f , and the matrix elements can also be presented as a sum over intermediate nuclear states |κ = |J π κ with certain angular momentum J κ , parity π, and energy E κ : where operators O α -with α denoting Gamow-Teller (GT), Fermi (F), tensor (T), etc. operators-contain neutrino potentials, spin and isospin operators, and explicit dependence on the intermediate state energy E κ . The most common of the operators can be found in Refs. [17,43], and they include vector and axial nucleon form-factors that take into account nucleon size effects. The calculation details for two-body matrix elements, 13|O α |24 , are discussed in Appendix D of Ref. [43]. Let us note that the two-body wave functions in the matrix elements (5) are not antisymmetrized, as one would expect for nuclear two-body matrix elements. The wave functions should be understood as |24 = |2 · |4 and |13 = |1 · |3 , where 1, 2, 3, and 4 represent single-nucleon quantum numbers, e.g., 1 = {τ 1z , n 1 , l 1 , j 1 , µ 1 } Calculations using a summation on intermediate states is very time-consuming, due to the need for obtaining a large number of intermediate states κ and the associated one-body transition densities f |ĉ † 3ĉ 4 |κ and κ|ĉ † 1ĉ 2 |i in Equation (5), which can only be conducted efficiently in J-scheme codes such as NuShellX code [92]. The results and analyses for most of the nuclei in Table 1 can be found in Refs. [16,43,45,48,70].
Although time-consuming, this method has the advantage of being applicable for a large class of effective nuclear Hamiltonian and transition operators. For example, it can be used for isospin-breaking nuclear Hamiltonians and with transition operators that are treating asymmetrically the initial neutron single particle (s.p.) states vs. the final proton s.p. states, such as the in-medium similarity renormalization group and realistic shell model methods. This method is always applicable for transitions to the 2 + states in the daughter nucleus, even in cases when the transition operator is not a rotational scalar anymore [83,84].
If one replaces the energies of the intermediate states in the form-factors by an average constant value, one obtains the closure approximation. The operators O α →Õ α ≡ O α ( E ) become energy-independent and the sum over the intermediate states in the nuclear matrix element (5) can be taken explicitly using the completeness relation: The advantage of this approximation is significant because it eliminates the need for calculating a very large number of states in the intermediate nucleus, which could be computationally challenging, especially for heavy systems. One needs only to calculate the two-body transition densities (see Section 3.2) between the initial and final nuclear states. This approximation is very good due to the fact that the values of q that dominate the matrix elements are of the order of 100-200 MeV, while the relevant excitation energies are only of the order of 10 MeV. The obvious difficulty related to this approach is that I have to find a reasonable value for this average energy, E , which can effectively represent the contribution of all the intermediate states. This average energy needs to account also for the symmetric part of the two-body matrix elements 13|O α |24 in Equation (7) below. Indeed, the two-body wave functions |13 and |24 are not antisymmetric; by replacing the energies of the intermediate states with a constant, only the antisymmetric parts of these matrix elements are taken into account.
Most reported calculations are using closure approximation with some closure energies taken from Ref. [93]. By comparing the closure and the summation method results for different isotopes in different model spaces, I find [48,70] the optimal closure energies for a given model space and effective Hamiltonian (see end of Section 3.2 for examples). The optimal closure energies for a given model space and effective Hamiltonian can then be found by performing a calculation for a (fictitious) 0νββ NME of lower complexity.

The 0νββ NME in Closure Approximation
In the closure approximation, the 0νββ NME can be reduced to a sum of the products of two-body transition densities (TBTD), defined by the right-hand side of Equation (6), and antisymmetrized two-body matrix elements, M 0ν α = ∑ j p j p j n j n J π TBTD j p j p , j n j n ; J π j p j p ; J π T | τ −1 τ −2 O α 12 | j n j n ; J π T , where O α 12 are the two-body operators corresponding to different transitions (here denoted by α = F, GT, T, Fq, GTq, . . . [17]) contributing to some of the diagrams of the 0νββ process in Figure 2 (see Ref. [17] for details). One should not confuse the isospin T in the two-body matrix elements with the tensor operator notation for α.
Having the two-body matrix elements ready, one can calculate the NME in Equation (7) if two-body transition densities TBTD j p j p , j n j n ; J π are known. Most of the shell model codes do not provide two-body transition densities. One alternative approach is to take advantage of the isospin symmetry that most of the effective interactions have, which creates wave functions with good isospin. The approach described below works also when the proton and neutron are in different shells. If the above conditions are satisfied, one can transform the two-body matrix elements of a change in isospin ∆T = 2 operator using the Wigner-Eckart theorem, from a change in isopspin projection ∆T z = −2 to ∆T z = 0, which can be further used to describe transitions between states in the same nucleus. Then the transformed matrix elements preserve spherical symmetry and they can be used as a piece of a Hamiltonian, H α ββ , which violates isospin symmetry, but it is a scalar with respect to rotational group. One can then lower by two units the isospin projection of the g.s. of the parent nucleus that has the higher isospin T > , e.g., that of 48 Ca, thus becoming an isobar analog state of the daughter nucleus that has isospin T < = T > − 2, e.g., in 48 Ti. Denoting by | 0 + i< T > > the transformed state, one can now calculate the many-body matrix elements of the transformed 0νββ operator, Choosing | 0 + i< T > > as a starting Lanczos vector and performing one Lanczos iteration with H α ββ , one obtains where | L 1 > is the new Lanczos vector. Then, one can calculate the matrix elements in Equation (8): The transition matrix elements in Equation (7) can then be recovered using again the Wigner-Eckart theorem, where C T 1 T 2 T T z1 T z2 T z are isospin Clebsch-Gordan coupling coefficients. This procedure can be implemented in most nuclear shell model codes. The transformation of the g.s. of the parent to an analog state of the daughter can be performed very quickly, and one Lanczos iteration represents a small load as compared with the calculation of the g.s. of the daughter. The additional calculations described in Equations (9)-(11) require smaller resources than those necessary to calculate the TBTDs.
The form of the NME described in Equation (7) assumes that the underlying manybody Hamiltonian and the resulting wave functions have good isospin symmetry. That might not be the case when the Coulomb interaction is included or/and ab initio approaches to obtain the effective shell model Hamiltonian, such as the in-medium similarity renormalization group [94] or realistic shell model [95], are used. In that case, one could project the parent and daughter wave functions on good isospin components and extend the above procedure for each pair of isospin components considering the appropriate jump in isospin (which might be a difference of 2). In practice, the contributions from the main isospin components described in the above procedure dominate.
There are also some limitations to this method. For example, the in-medium similarity renormalization group and realistic shell model methods, as well as the G-matrix-like re-summation approach (see, e.g., [96] and references therein) also provide effective operators, which breaks the symmetry of the two-body matrix elements, j p j p ; J π T | τ −1 τ −2 O α 12 | j n j n ; J π T between the nn and pp two-body states of the initial and final nucleus. In that case, one could consider the average of the corresponding two-body matrix elements [96]. This method is not always applicable for transitions to the 2 + states in the daughter nucleus because if one uses the contributions from the right-handed currents, such as those of the λ and η mechanisms (see Section 2.1), then the transition operator is not a rotational scalar anymore [83,84]. In those cases, the H α ββ cannot be defined and the TBTDs are needed. It would be interesting to compare the half-lives of two isotopes to identify the dominant mechanisms contributing to 0νββ decays. It is often the case that even within the shell model approach, using two different effective Hamiltonians leads to significantly different results. Given that the conclusion of two-isotope analysis is sensitive to the accuracy of NMEs, it is important to consider at least two sets of effective Hamiltonians. In addition, for consistency, I use the optimal closure energies [43,45,48], with E corresponding to each Hamiltonian and model space. One set of NMEs is obtained using the Hamiltonians preferred by our CMU (Central Michigan University) group: for 48 Ca in the p f model space (0 f 7/2 , 1p 3/2 , 0 f 5/2 , 1p 1/2 ), I use GXPF1A effective Hamiltonian [97] with E = 0.5 MeV; for 76 Ge and 82 Se in the jj44 model space (0 f 5/2 , 1p 3/2 , 1p 1/2 , 0g 9/2 ), I choose JUN45 [98] with E = 3.4 MeV; and for 124 Sn, 130 Te, and 136 Xe in the jj55 model space (0g 7/2 , 1d 5/2 , 1d 3/2 , 1s 1/2 , 0h 11/2 ), I use SVD effective Hamiltonian [99] with E = 3.5 MeV. The second set of NMEs I calculate using the Hamiltonians preferred by the Strasbourg-Madrid group: in this case, for 48 Ca I use KB3G [100] with E = 2.5 MeV, for 76 Ge I use 82 Se GCN.28-50 with E = 10 MeV, and for 130 Te and 136 Xe I use GCN.50-82 with E = 12 MeV [101] (see Section 3.1 for the definition of the optimal closure energies E ). The numerical analysis is given in Ref. [17], where I find that using the ratio of experimental half-lives one could identify if a selected few mechanisms may be dominant.

Two-Neutrino Double Beta Decay
Two-neutrino double beta decay (2νββ) is an associated process allowed by the Standard Model, which was observed experimentally for about 10 isotopes, including most in Table 1. Here, I describe an improved spectra-function technique for calculating associated NMEs in very large model spaces in which a direct summation on intermediate states is not practical. For the 2νββ mode, the relevant NMEs are of Gamow-Teller type, and have the following expression for decays to states in the grand-daughter that have the angular momentum J = 0 [93]: Here, E k is the excitation energy of the 1 + k state of the intermediate odd-odd nucleus and E 0 = 1 2 Q ββ + ∆M. Q ββ is the Q-value corresponding to the ββ decay to the final 0 + f state of the grand-daughter nucleus, and ∆M is the mass difference between the parent (e.g., 48 Ca) and the intermediate nucleus (e.g., 48 Sc). The most common case is the decay to the 0 + 1 g.s. of the grand-daughter, but decays to the first excited 0 + 2 state were also observed [80].
The 2νββ decay half-life is given by In Ref. [51], I fully diagonalized 250 1 + states of the intermediate nucleus 48 Sc in the p f valence space to calculate the 2νββ NME for 48 Ca. This method of the direct diagonalization of a large number of states can be used for somewhat heavier nuclei using the J-scheme shell model code NuShellX [92], but for large-dimension cases one needs an alternative method. In particular, the m-scheme dimensions needed for the 48 Ca NME calculations when taking into account up to 2hω excitations sd-p f valence space are larger than 1 billion (716 million for 48 Sc). These large dimensions also require a method more efficient than direct diagonalization. The pioneering work on 48 Ca [102] used a strengthfunction approach that converges after a small number of Lanczos iterations, but it requires large-scale shell model diagonalization when one wants to check the convergence. Ref. [103] proposed an alternative method which converges very quickly, but it did not provide full recipes for all its ingredients, and it was never used in practical calculations. Here, I propose a simple numerical scheme to calculate all coefficients of the expansion proposed in Ref. [103]. Following Ref. [103], I choose as a starting Lanczos vector L ± 1 either the initial or final states in the decay (only 0 + to 0 + transitions are considered here), on which is applied the Gamow-Teller operator, |στ The results are the "door-way" states |dw ± > multiplied by the constants c ± , which represent the square roots of the respective Gamow-Teller sum rule. Ref. [103] showed that the matrix element in Equation (12) could be calculated using one of the following two equations: Here, the sum is over the Lanczos vectors L m . One can show that the g ± m factors can be calculated with the following formula after N Lanczos iterations: Here, V mk are the eigenvectors of the N-order Lanczos matrix corresponding to eigenvalue E N L (1 + k ). The advantage of using Equations (14)- (18) is that in order to check the convergence at each iteration one only needs the Lanczos vectors, which have to be stored anyway, and not the eigenvectors of the many-body Hamiltonian. The g ± m can be calculated very quickly, and only the last overlap in the sum of Equation (16) or Equation (17) needs to be calculated at each iteration. This algorithm can provide a gain in efficiency by a factor of about two as compared with the strength-function approach of Ref. [102].
Another advantage of this method is that it can be used with both M-scheme and Jscheme shell model codes, while a direct summation in Equation (12) on the 1 + states in the intermediate nucleus can only be performed using J-scheme codes. The method described here requires about 20 Lanczos iterations for convergence. I estimate (see, e.g., [51]) that good convergence for the direct summation in Equation (12) requires about 300-500 1 + that usually can be achieved with about 5000-10,000 iterations. Given the input/output burden associated with so many iterations, I estimate computational speed improvement by a factor of about 1000 in the present method as compared with the direct summation method.
It is known that a good comparison of the shell model results with experimental data requires a multiplicative quenching factor for the Gamow-Teller operator. This numerical analysis when compared with the experimental data [45,[49][50][51]70] indicates that for the selected effective Hamiltonians, only quenching factors between 0.6 and 0.74 are needed.

Conclusions
In this paper, I provide an overview of the double beta decay process and describe in some detail the shell model approach for the calculation of the nuclear matrix elements necessary for the analysis of experimental data.
Analyzing the physics related to the neutrinoless double beta decay process, one observes that it would entail physics beyond the Standard Model, namely the lepton number violation, which may lead to the conclusion that neutrinos may be the only known Standard Model fermions that are of Majorana type. This information may be crucial for properly extending the Standard Model Lagrangian to describe the observed neutrino masses and other LNV processes.
I describe the 0νββ decay half-life using BSM mechanisms induced by new particles such as the left-right symmetric model or SUSY and also use a more general EFT approach that includes the most general LNV addition to the Standard Model Lagrangian. Both approaches lead to similar numbers of NMEs associated with either model-specific or EFT-linked LNV couplings.
The largest part of the paper is dedicated to the techniques for calculating the needed NMEs within the shell model approach. For 0νββ, I analyze the different scenarios under which the NMEs can be calculated in the closure approximation that is good within 10%. I also describe how to calculate the same NMEs beyond closure and identify optimal closure energies which can minimize the error of the less time-consuming closure approximation.
Two-neutrino double beta decay is an associated process allowed by the Standard Model, which was observed experimentally for about 10 isotopes. Here, I describe an improved spectra-function technique of calculating the associated NMEs in very large model spaces for which the dwerect summation on intermediate states is impractical.
Finally, although most of the paper reviews results already published, some new results regarding techniques for calculating 0νββ and 2νββ NMEs in extreme situations can be found at the end of Sections 3.2 and 4.

Conflicts of Interest:
The author declares no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: