*3.1. The Anatomy of the* 0*νββ NMEs*

The nuclear matrix elements needed in Equations (1)–(4) describe the transition from an initial nucleus <sup>|</sup>*i* <sup>=</sup> <sup>|</sup>0<sup>+</sup> *<sup>i</sup>* to a final nucleus <sup>|</sup> *<sup>f</sup>* <sup>=</sup> <sup>|</sup>0<sup>+</sup> *<sup>f</sup>* , and the matrix elements can also be presented as a sum over intermediate nuclear states <sup>|</sup>*κ* <sup>=</sup> <sup>|</sup>*J<sup>π</sup> <sup>κ</sup>* with certain angular momentum *Jκ*, parity *π*, and energy *Eκ*:

$$M\_{\mathfrak{A}}^{0\upsilon} = \sum\_{\kappa} \sum\_{1234} \langle 13 | \mathcal{O}\_{\mathfrak{A}} | 24 \rangle \langle f | \mathfrak{E}\_3^{\dagger} \mathfrak{E}\_4 | \kappa \rangle \langle \kappa | \mathfrak{E}\_1^{\dagger} \mathfrak{E}\_2 | i \rangle,\tag{5}$$

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


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* |*c*ˆ † <sup>3</sup>*c*ˆ4|*κ* and *κ*|*c*ˆ † <sup>1</sup>*c*ˆ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<sup>+</sup> 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 <sup>O</sup>*<sup>α</sup>* <sup>→</sup> <sup>O</sup>˜ *<sup>α</sup>* ≡ 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:

$$
\sum\_{\kappa} \langle f | \mathfrak{c}\_3^{\dagger} \mathfrak{c}\_4 | \kappa \rangle \langle \kappa | \mathfrak{c}\_1^{\dagger} \mathfrak{c}\_2 | i \rangle = \langle f | \mathfrak{c}\_3^{\dagger} \mathfrak{c}\_4 \mathfrak{c}\_1^{\dagger} \mathfrak{c}\_2 | i \rangle. \tag{6}
$$

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.
