*3.2. 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\_{\mathfrak{A}}^{\mathbb{O}\nu} = \sum\_{j\_{\mathcal{V}} j\_{\mathcal{V}} j\_{n} j\_{n'} I \pi} \text{TBTD} \left( j\_{\mathcal{V}} j\_{p'}, j\_{n} j\_{n'}; J\_{\pi} \right)$$

$$\left\langle j\_{p} j\_{p'}; J^{\pi} T \mid \pi\_{-1} \pi\_{-2} O\_{12}^{\mathfrak{a}} \mid j\_{n} j\_{n'}; J^{\pi} T \right\rangle\_{\prime} \tag{7}$$

where *O<sup>α</sup>* <sup>12</sup> 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 *jp jp* , *jn jn* ; *J<sup>π</sup>* 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 Δ*Tz* = −2 to Δ*Tz* = 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<sup>α</sup> ββ*, 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 48Ca, thus becoming an isobar analog state of the daughter nucleus that has isospin *T*<sup>&</sup>lt; = *T*<sup>&</sup>gt; − 2, e.g., in 48Ti. Denoting by <sup>|</sup> <sup>0</sup><sup>+</sup> *<sup>i</sup>*<sup>&</sup>lt; *T*<sup>&</sup>gt; > the transformed state, one can now calculate the many-body matrix elements of the transformed 0*νββ* operator,

$$M\_{\mathfrak{a}}^{0v}(T\_z = T\_{<}) = \left\langle 0\_f^+ \, T\_{<} \, \big|\, H\_{\mathfrak{f}\mathfrak{f}}^{\mathfrak{a}} \mid 0\_{\mathfrak{i} < \mathfrak{i}}^+ \, T\_{\ge} \right\rangle. \tag{8}$$

Choosing <sup>|</sup> <sup>0</sup><sup>+</sup> *<sup>i</sup>*<sup>&</sup>lt; *T*<sup>&</sup>gt; > as a starting Lanczos vector and performing one Lanczos iteration with *H<sup>α</sup> ββ*, one obtains

$$H\_{\beta\beta}^a \mid 0\_{i<}^+ \ T\_{\succeq}> = a\_1 \mid 0\_{i<}^+ \ T\_{\succeq}> + b\_1 \mid L\_1 >,\tag{9}$$

where | *L*<sup>1</sup> > is the new Lanczos vector. Then, one can calculate the matrix elements in Equation (8):

$$M\_a^{0\upsilon}(T\_z = T\_{<}) = b\_1 \langle 0\_f^+ | T\_{<} \mid L\_1 \rangle. \tag{10}$$

The transition matrix elements in Equation (7) can then be recovered using again the Wigner–Eckart theorem,

$$M\_{\mathfrak{a}}^{0\upsilon} = M\_{\mathfrak{a}}^{0\upsilon}(T\_{\overline{z}} = T\_{\leq}) \times \mathbb{C}\_{T\_{>}}^{T\_{>}} \stackrel{2}{\underset{T\_{\leq}}{T\_{\leq}}} \stackrel{T\_{\leq}}{T\_{\leq}} / \mathbb{C}\_{T\_{\leq}}^{T\_{>}} \stackrel{2}{\underset{T\_{\leq}}{T\_{\leq}}} \tag{11}$$

where *CT*<sup>1</sup> *<sup>T</sup>*<sup>2</sup> *<sup>T</sup> Tz*<sup>1</sup> *Tz*<sup>2</sup> *Tz* 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 resummation approach (see, e.g., [96] and references therein) also provide effective operators, which breaks the symmetry of the two-body matrix elements,

$$\left\langle j\_p j\_{p'} ; J^\pi T \mid \pi\_{-1} \pi\_{-2} O\_{12}^\kappa \mid j\_n j\_{n'} ; J^\pi T \right\rangle\_+$$

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<sup>+</sup> 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<sup>α</sup> ββ* 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 48Ca in the *p f* model space (0 *f*7/2, 1*p*3/2, 0 *f*5/2, 1*p*1/2), I use GXPF1A effective Hamiltonian [97] with *E* <sup>=</sup> 0.5 MeV; for 76Ge and 82Se in the *jj*44 model space (<sup>0</sup> *<sup>f</sup>*5/2, 1*p*3/2, 1*p*1/2, 0*g*9/2), I choose JUN45 [98] with *E* <sup>=</sup> 3.4 MeV; and for 124Sn, 130Te, and 136Xe in the *jj*<sup>55</sup> model space (0*g*7/2, 1*d*5/2, 1*d*3/2, 1*s*1/2, 0*h*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 48Ca I use KB3G [100] with *E* <sup>=</sup> 2.5 MeV, for 76Ge I use 82Se GCN.28-50 with *E* <sup>=</sup> 10 MeV, and for 130Te and 136Xe 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.

## **4. 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]:

$$M\_{2\nu} = \sum\_{k} \frac{\langle 0\_f^+ || \sigma r^- || 1\_k^+ \rangle \langle 1\_k^+ || \sigma r^- || 0\_i^+ \rangle}{(E\_k + E\_0)}.\tag{12}$$

Here, *Ek* is the excitation energy of the 1<sup>+</sup> *<sup>k</sup>* state of the intermediate odd-odd nucleus and *E*<sup>0</sup> = <sup>1</sup> <sup>2</sup>*Qββ* <sup>+</sup> <sup>Δ</sup>*M*. *<sup>Q</sup>ββ* is the Q-value corresponding to the *ββ* decay to the final 0<sup>+</sup> *f* state of the grand-daughter nucleus, and Δ*M* is the mass difference between the parent (e.g., 48Ca) and the intermediate nucleus (e.g., 48Sc). The most common case is the decay to the 0<sup>+</sup> <sup>1</sup> g.s. of the grand-daughter, but decays to the first excited 0<sup>+</sup> <sup>2</sup> state were also observed [80].

The 2*νββ* decay half-life is given by

$$\left[T\_{1/2}^{2\nu}\right] = G\_{2\nu} \cdot g\_A^4 \cdot \left(m\_c c^2 \cdot M\_{2\nu}\right)^2\tag{13}$$

In Ref. [51], I fully diagonalized 250 1<sup>+</sup> states of the intermediate nucleus 48Sc in the *p f* valence space to calculate the 2*νββ* NME for 48Ca. 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 48Ca NME calculations when taking into account up to 2*h*¯ *ω* excitations *sd*-*p f* valence space are larger than 1 billion (716 million for 48Sc). These large dimensions also require a method more efficient than direct diagonalization. The pioneering work on 48Ca [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*± <sup>1</sup> either the initial or final states in the decay (only 0<sup>+</sup> to 0<sup>+</sup> transitions are considered here), on which is applied the Gamow–Teller operator,

$$|\sigma \tau^- 0\_i^+ >= c\_-|dw\_- > \equiv c\_-|L\_1^- > \,, \tag{14}$$

$$|\sigma \pi^+ 0\_f^+ \succ = c\_+|dw\_+ \succ \equiv c\_+|L\_1^+ \succ \ . \tag{15}$$

The results are the "door-way" states |*dw*<sup>±</sup> > 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:

$$M\_{2\nu}(0^+) \approx 3\mathfrak{c} + \mathfrak{c} - \sum\_{\mathfrak{m}} \mathfrak{g}\_{\mathfrak{m}}^- < d\mathfrak{w} + \left| L\_{\mathfrak{m}}^- \right. \\ \left. \right| \equiv M\_{2\nu}^{\mathrm{GT}-} \,, \tag{16}$$

$$M\_{2\nu}(0^+) \approx 3\mathcal{c} + \mathcal{c}\_- \sum\_m \mathcal{g}\_m^+ < L\_m^+|dw\_-> \equiv M\_{2\nu}^{\text{GT}+}.\tag{17}$$

Here, the sum is over the Lanczos vectors *Lm*. One can show that the *g*± *<sup>m</sup>* factors can be calculated with the following formula after *N* Lanczos iterations:

$$g\_m^{\pm} = \sum\_{k=1}^{N} \frac{V\_{1k}^{\pm} V\_{mk}^{\pm}}{E\_L^N(1\_k^+) - E\_{\mathbb{g}.s.} + E\_0} \,. \tag{18}$$

Here, *Vmk* are the eigenvectors of the N-order Lanczos matrix corresponding to eigenvalue *E<sup>N</sup> <sup>L</sup>* (1<sup>+</sup> *<sup>k</sup>* ). 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*± *<sup>m</sup>* 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<sup>+</sup> 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<sup>+</sup> 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–51,70] indicates that for the selected effective Hamiltonians, only quenching factors between 0.6 and 0.74 are needed.

#### **5. 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.

**Funding:** This research received support from the US Department of Energy grant DE-SC0022538 "Nuclear Astrophysics and Fundamental Symmetries". The author acknowledges a Central Michigan University Faculty Research and Creative Endeavors (FRCE) type A grant for travel support.

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