Present Status of Nuclear Shell-Model Calculations of Neutrinoless Double-Beta Decay Matrix Elements

Neutrinoless double beta decay searches are currently among the major foci of experimental physics. The observation of such a decay will have important implications in our understanding of the intrinsic nature of neutrinos and shed light on the limitations of the Standard Model. The rate of this process depends on both the unknown neutrino effective mass and the nuclear matrix element associated with the given neutrinoless double-beta decay transition. The latter can only be provided by theoretical calculations, hence the need of accurate theoretical predictions of the nuclear matrix element for the success of the experimental programs. This need drives the theoretical nuclear physics community to provide the most reliable calculations of the nuclear matrix elements. Among the various computational models adopted to solve the many-body nuclear problem, the shell model is widely considered as the basic framework of the microscopic description of the nucleus. Here, we review the most recent and advanced shell-model calculations of the nuclear matrix elements considering the light-neutrino-exchange channel for nuclei of experimental interest. We report the sensitivity of the theoretical calculations with respect to variations in the model spaces and the shell-model nuclear Hamiltonians.


Introduction
Neutrinoless double beta (0νββ) decay is a process in which two neutrons inside the nucleus transform into two protons with the emission of two electrons and no neutrinos. With only two electrons in the final state, lepton number conservation needs to be violated by two units for the process to occur. The observation of 0νββ decay would imply that neutrinos are Majorana particles-that is, they are their own antiparticles [1]-and give insight into leptogenesis scenarios for the generation of the observed matter-antimatter asymmetry in the universe [2]. In fact, 0νββ decay is the most promising laboratory probe of lepton number violation and it is, in fact, the subject of intense experimental activities [3][4][5][6][7][8][9][10][11][12][13][14]. The current reported 0νββ half-lives of 76 Ge, 130 Te and 136 Xe are larger than 8 · 10 25 yr [13], 1.5 · 10 25 yr [14] and 1.1 · 10 26 yr [7], respectively, with next generation ton scale experiments expecting a two orders of magnitude improvement in the half life sensitivity.
Since 0νββ decay measurements use atomic nuclei as a laboratory to test the extent of the Standard Model, nuclear theory plays a crucial role in correctly interpreting the experimental data and disentangling arXiv:2011.14734v1 [nucl-th] 30 Nov 2020 nuclear physics effects from unknown lepton number violating mechanisms. The half-life of the 0νββ decays is given by where G 0ν is a phase-space (or kinematic) factor [15,16], M 0ν is the nuclear matrix element (NME), and f (m i , U ei ) is a function of the neutrino masses m i and their mixing matrix elements U ei that accounts for beyond-standard-model physics. Within the light-neutrino exchange mechanisms, f (m i , U ei ) has the following expression: where g A is the axial coupling constant, m e is the electron mass, and m ν = ∑ i (U ei ) 2 m i is the effective neutrino mass. It is then clear that access to unknown neutrino properties is granted only if M 0ν is calculated with great accuracy. Currently the calculated matrix elements for nuclei of experimental interest are characterized by large uncertainties. For nuclear systems with medium masses and beyond, the many-body nuclear problem cannot be solved exactly with the available computational resources. For these systems, one is inevitably forced to truncate the model spaces and reduce or neglect the effects of many-body correlations and electroweak currents. As a result, different computational methods provide calculated M 0ν swhich differ by a factor of two [17]. On the above grounds, it is clear that reliable calculations of M 0ν 's are a prime goal of nuclear many-body investigations.
The ab initio framework for nuclei allows to retain the complexity of many-body correlations and currents. Within this approach nuclei are described as systems made of nucleons interacting via twoand three-body forces. Interactions with external probes, such as electrons, neutrinos, and photons are described using many-body current operators. One-and two-body current operators describe external probes interacting with individual nucleons and pairs of correlated nucleons, respectively. This scheme has been implemented successfully to study light to intermediate mass nuclei within several many-body computational approaches. Due to their prohibitive computational cost, ab initio methods have been used to study 0νββ transitions in light nuclei instead. While transitions in light nuclei do not have a direct experimental interest, these studies provide us with an important benchmark to test other many-body methods that can be used to calculate transition matrix elements for heavy-mass nuclei of experimental interest. Further, they allow us to size the importance of the different lepton number violating mechanisms leading to 0νββ decay processes, and to quantify the effect of the various approximations used in the many-body methods for medium to large nuclear systems. Studies along this line have been carried out, for example, in Refs. [18][19][20]. Only very recently, the ab initio community is venturing calculations of 0νββ decay matrix element of experimental relevance, as reported, for example, in Ref. [21], where the authors calculate the M 0ν of 48 Ca-the lightest system where the Q-value is compatible with the decay-combining the in-medium similarity renormalization group with the generator-coordinate method [22].
Besides the exceptions mentioned above, the nuclear physics community has been primarily focused on employing approximated many-body methods to access heavy open-shell nuclei of experimental interest. These approximated methods generally invoke a truncation of the full Hilbert space of configurations. To account for missing dynamics and degrees of freedom in the nuclear wave functions, the nuclear Hamiltonian is then replaced by an effective or renormalized Hamiltonian, i.e., H eff . This procedure is carried out, in general, by fitting parameters inherent the given nuclear model to spectroscopic properties of the nuclei under investigation. Nuclear models adopted to study 0νββ decay of nuclei of experimental interest are: the interacting boson model [23][24][25], the quasiparticle random-phase approximation (QRPA) [26][27][28], the energy density functional methods [29], the covariant density-functional theory [30][31][32][33][34][35][36], and the shell model (SM) [37][38][39][40][41]. These models agree within a factor of two (see, for example, Fig. 5 of Ref. [17] and references therein) when calculating 0νββ decay matrix elements of A = 48 − 150 nuclei. The difference is mostly to be ascribed to the different renormalization procedures adopted by the different models.
In addition to renormalize the nuclear Hamiltonian, in this scheme one has to renormalize the free constants that appear in the definitions of the decay operators-e.g., proton and neutron electric charges, spin and orbital gyromagnetic factors, etc. For example, the axial coupling constant g f ree A = 1.2723 [42] needs to be quenched by a factor of q [43], because all the aforementioned models usually overestimate Gamow-Teller (GT) rates when compared to the experimental data [44]. The choice of q depends on the nuclear structure model, the dimensions of the reduced Hilbert space, and the mass of the nuclei under investigation [45]. The common procedure to handle the quenching of g A is to fit GT related data (e.g., single-β decay strengths, two-neutrino double-β decay rates, etc.), and some authors argue that the value of q required to reach agreement between theoretical and experimental values should be also employed to calculate M 0ν (see for instance Refs. [45,46]). In passing, it is worth mentioning that within the ab initio framework one can utilize the free nucleonic charges, magnetic moments, and axial coupling constant without having to resort to quenching, provided that corrections from two-body currents and two-body correlations are accounted for [47][48][49][50][51][52][53][54]. In this work, we review the most recent and advanced SM results of M 0ν for nuclei currently candidates for the detection of the 0νββ decay in many laboratories around the world. We focus on the sensitivity of the calculations with respect to variations in the model spaces and the shell-model nuclear Hamiltonians, as well as to the variations in the "short-range correlations" which reveal the role of SM correlated wave functions.
The paper is organized as follows: in Section 2 we outline the basics theory of the nuclear SM and short-range correlations, and provide the analytical expressions of the nuclear matrix elements, for both neutrinoless and two-neutrino double beta decay. The latter are reported to assess the validity of the adopted nuclear wave functions. In fact, a comparison with experimental data is clearly possible for two-neutrino double beta decays. Section 3 is devoted to the results of the latest SM calculations for 48 Ca→ 48 Ti, 76 Ge→ 76 Se, 82 Se→ 82 Kr, 130 Te→ 130 Xe, and 136 Xe→ 136 Ba 0νββ and 2νββ decays. Comparisons between experimental and calculated M 2ν s are reported at the end of Section 3.2 with a discussion on the g A quenching. Our conclusions are given in Section 4.

The nuclear shell model
The nuclear shell model allows for a microscopic description of the structure of the nucleus [55,56], and it is the root of most current ab initio approaches (No-Core Shell Model, Coupled Cluster Method, In-Medium Similarity Renormalization Group). It is based on the ansatz that each nucleon inside the nucleus moves independently in a spherically symmetric mean field generated by all other constituents. The mean field is usually described by a Woods-Saxon or a harmonic oscillator (HO) potential supplemented by a strong spin-orbit term.
This basic version of the shell model successfully explains the appearance of protons and/or neutrons "magic numbers"-characterizing nuclei bounded more tightly with respect to their neighbors-along with several nuclear properties [57], including angular momenta and parity for ground-states of odd-mass nuclei. Within this framework, nucleons arrange themselves in well defined and separated energy levels, i.e., the "shells". It is worth emphasizing that shell-model wave functions do not include correlations induced by the strong short-range two-nucleon interaction. We will come back to this point later, when we discuss the "short-range correlations" (SRC).
The SM can be further improved, especially its description of low-energy nuclear structure, introducing the "interacting shell model" (ISM) picture. In the ISM, the complex nuclear many-body problem is reduced to a simplified one where only few valence nucleons interact in the reduced model space spanned by a single major shell above an inert core. The valence nucleons interact via a two-body "residual interaction", that is the part of the interaction which is not already accounted for in the central potential. The inclusion of the residual interaction removes the degeneracy of the states belonging to the same configuration and produces a mixing of different configurations.
The SM Hamiltonian consists of one-and two-body components and associated parameters, namely the single-particle (SP) energies and the two-body matrix elements (TBMEs) of the residual interaction. These parameters account for the degrees of freedom that are not explicitly included in the truncated Hilbert space of configurations. As a matter of fact, SP energies and TBMEs should be determined to include, in an effective way, the excitations of both the core and the valence nucleons into the shells above the model space.
The construction of the effective SM Hamiltonian, H eff , can be carried out into two distinct ways. In one approach one starts from realistic two-and three-nucleon forces (see Ref. [58] and references therein for a review on realistic two-and three-nucleon potentials) and derive the effective Hamiltonian from them. The H eff will then have eigenvalues that belong to the set of eigenvalues of the full nuclear Hamiltonian, defined in the whole Hilbert space. The alternative approach is phenomenological. In this case, the SM Hamiltonian one-and two-body components are adjusted to reproduce a selected set of experimental data. For the fitting procedure one could i) use adjustable parameters entering the analytical expressions of the residual interaction, ii) directly consider the Hamiltonian matrix elements as free parameters (see, e.g., Refs. [59,60]), or iii) fine tune the TBMEs of a realistic H eff to reproduce the experimental results. The phenomenological approach has been widely utilized since its formulation in the fifties, and it successfully reproduces a huge amount of data and describes some of the most fundamental properties of the structure of atomic nuclei [61].
The SM provides suitable and well tested nuclear wave functions for the initial and final states entering the calculation of M 0ν associated with 0νββ decays. SM results based on both the realistic and phenomenological H eff are reported in Section 3.

Short-range correlations
Short-range correlations (SRCs) are required to account for physics that is missing in all models that expand nuclear wave functions in terms of a truncated non-correlated SP basis. In particular, two-body decay operators-such as those entering 0νββ decays-acting on an unperturbed (uncorrelated) wave function yield results that are intrinsically different from those obtained acting the real (correlated) nuclear wave function [62,63]. Due to the highly repulsive nature of the short-range two-nucleon interaction, and in order to carry out nuclear structure calculations, one is forced to perform a consistent regularization of the two-nucleon potential, V NN , and of any two-body transition operator [64].
In nuclear structure calculations based on realistic potentials, one has to deal with non-zero values of the non-correlated wave function, Φ, at short distances. This can be appreciated in Fig. 1. However, because of the repulsive nature of the V NN interaction at small inter-particle distances (or equivalently, the repulsive V NN behavior at high-momentum) the correlated wave function, Ψ, has to approach zero as the inter-nucleon distance diminishes, and as fast as the core repulsion increases, see  To remedy to this shortcoming, one has to renormalize the short-range (high-momentum) components of the V NN potential whenever a perturbative approach to the many-body problem is pursued. The most common way to soften the matrix elements of the 0νββ decay operator is to include SRC given by Jastrow type functions [65,66]. Recently, SRC have been modeled using the Unitary Correlation Operator Method (UCOM) [38,63]. This approach prevents the overlap between wave functions of a pair of nucleons [67].
Another approach has been proposed by some of the present authors in Refs. [68,69], where the renormalization of the the 0νββ two-body decay operator is carried out consistently with the V low-k procedure [70] adopted to renormalize the repulsive high-momentum components of the V NN potential. In particular, the renormalization of V NN occurs through a unitary transformation, Ω, which decouples the full momentum space of the two-nucleon Hamiltonian, H NN , into two subspaces; the first one is associated with the relative-momentum configurations below a cutoff Λ and is specified by a projector operator P, the second one is defined in terms of its complement Q = 1 − P [68]. Being a unitary transformation, Ω preserves the physics of the original potential for the two-nucleon system, e.g., the calculated values of all nucleon-nucleon observables are the same as those reproduced by solving the Schrödinger equation for two nucleons interacting via V NN .
The two-body 0νββ operator, Θ, is calculated in the momentum space and renormalized using Ω. This ensure a consistency with the V NN potential, whose high-momentum (short range) components are dumped by the introduction of the cutoff Λ. The Θ vertices appearing in the perturbative expansion of thê Θ box are substituted with the renormalized Θ low-k operator. The latter is defined as Θ low-k ≡ PΩΘΩ −1 P for relative momenta k < Λ, and is set to zero for k > Λ. The magnitude of the overall effect of this renormalization procedure is comparable to using the SRC modeled by the Unitary Correlation Operator Method [38], that is a lighter softening of M 0ν with respect to the one provided by Jastrow type SRC [68].

The 0νββ-decay operator for the light-neutrino exchange
We now turn our attention to the vertices of the bare 0νββ operator, Θ, for the light-neutrino-exchange channel [17].
We recall that the formal expression of M 0ν α -where α stands for Fermi (F), Gamow-Teller (GT), or tensor (T) decay channels-is written in terms of the one-body transition-density matrix elements between the daughter and parent nuclei (grand-daughter and daughter nuclei) k|a † p a n |i ( f |a † p a n |k ). Here, p and n denote proton and neutron states, and i, k, f refer to the parent, daughter, and grand-daughter nuclei, respectively, while M 0ν α reads [71,72]: where the tilde denotes a time-conjugated state,ã jm = (−1) j+m a j−m . The operators Θ α given by [17]: where H α are the neutrino potentials defined as: In the equation above, R = 1.2A 1/3 fm, j n α (qr) is the spherical Bessel function, n α = 0 for Fermi and Gamow-Teller components, while n α = 2 for the tensor component. For the sake of clarity, the explicit expressions [17] of neutrino form functions, h α (q), for light-neutrino exchange are reported below: For the vector, g V (q 2 ), axial-vector, g A (q 2 ), and weak-magnetism, g M (q 2 ), form factors we use the dipole approximation: where , and the cutoff parameters Λ V = 850 MeV and The total nuclear matrix element M 0ν can be then written as The expression in Eq. (2) can be easily calculated within the QRPA computational approach, while all other models-including most of the SMs-have to resort to the closure approximation. This approximation is based on the observation that the relative momentum q of the neutrino, appearing in the propagator of Eq. (6), is of the order of 100-200 MeV [17], while the excitation energies of the nuclei involved in the transition are of the order of 10 MeV [71]. It is then customary to replace the energies of the intermediate states, E k , appearing in Eq. (6), by an average value E k − (E i + E f )/2 → E . This allow us to simplify both Eqs. (2) and (6). In particular, M 0ν α can be re-written in terms of the two-body transition-density matrix elements f |a † p a n a † p a n |i as M 0ν α = ∑ j n j n j p j p f |a † p a n a † p a n |i and the neutrino potentials become Most SM calculations adopt the closure approximation to define the Θ operators given in Eqs.
(3)-(5), and take the average energies E from the evaluations of Refs. [73,74]. It is important to point out that the authors of Ref. [71] performed SM calculations of M 0ν for 48 Ca both within and beyond the closure approximation, and found that in the second case the results are ∼ 10% larger. In most cases, short-range correlations are included when computing the radial matrix elements of the neutrino potentials ψ nl (r)|H α |ψ n l (r) . In particular, the HO wave functions ψ nl (r) and ψ n l (r) are corrected by a factor [1 + f (r)], which takes into account the short range correlations induced by the nuclear interaction The functional form of the correlation function f (r) is usually written using a Jastrow-like parametrization as where a, b and c are parameters whose values depend on the renormalization procedure adopted to renormalize the non-correlated HO wave functions, (e.g., Jastrow or UCOM schemes, see Section 2.2 for details). In Table 1 we report the values of the a, b and c constants commonly employed in SM calculations. In addition to the values proposed by Miller and Spencer [65], we show those based on the modern nucleon-nucleon interactions CD-Bonn and AV18 and derived in Ref. [26].

The 2νββ-decay operator
As already pointed out in the Introduction, because of the impossibility to compare the theoretical values of M 0ν with the experiment, one has to find another way to check the reliability of the computed results. A viable route that is often considered in literature is the calculation of the standard or ordinary two-neutrinos double beta decay transitions where one observes the emission of two electrons and two antineutrinos. Two-neutrino double beta decays are simply the occurrence of two single beta decay transitions inside a nucleus. They differ from 0νββ decays in the characteristic value of momentum transfer, which is negligible in ordinary decays and of the order of hundreds of MeVs in 0νββ decay. Here, we list the expressions of the GT and Fermi components of the two-neutrinos double beta decay matrix elements M 2ν , namely In the equation above, E n is the excitation energy of the J π = 0 + n , 1 + n intermediate state, and E 0 = It should be pointed out that the Fermi component is zero in Hamiltonians that conserve the isospin, and most of the SM effective Hamiltonians do. It would play a marginal role only when isospin violation mechanisms are introduced, for example, to account for the effects of the Coulomb force acting between the valence protons [73,75]. In practice, in most calculations, the Fermi component is neglected altogether.
An efficient way to calculate M 2ν is to resort to the Lanczos strength-function method [61], which allows to include the intermediate states required to obtain a given accuracy for the calculated values.
The theoretical values are then compared with the experimental counterparts, that are extracted from the observed half life T 2ν One can base the evaluation of the M 2ν GT on the closure approximation, commonly adopted to study 0νββ-decay NMEs [73]. Within this approximation, one can avoid to explicitly calculate the intermediate J π = 1 + n states. The drawback is that, in using the closure on the intermediate states, the two one-body transition operators become a two-body operator.
This approximation is more adapt to evaluate M 0ν where the neutrino's momentum is about one order of magnitude greater than the average excitation energy of the intermediate states. This allows to safely neglect intermediate-state-dependent energies from the energy denominator appearing in the neutrino potential (see discussion in Section 2.3). Conversely, the closure approximation has turned out to be unsatisfactory when used to calculate M 2ν , and that is because the momentum transfer in 2νββ process are much smaller. Theoretical calculations of M 2ν are discussed in the next session.

Shell-model results
In this Section, we report SM results for M 0ν based on both the phenomenological and realistic H eff 's. All the calculations are based on the light-neutrino-exchange hypothesis, and the values of all the input parameters are the same as reported in Section 2.3. The only exception is the g A parameter, whose adopted value is equal to 1.254 in some reported calculations. It is worth pointing out, however, that in Ref. [76], where it can be found a detailed analysis of the sensitivity of the M 0ν results on the values of the input parameters, it has been shown that the effects of such a tiny difference in g A are negligible. We will focus our attention on the 48 Ca, 76 Ge, 82 Se, 128 Te, 130 Te and 136 Xe emitters. These results have been obtained performing a complete diagonalizations of H eff . The latter has been defined in different valence spaces tailored for the specific decay under investigation. All the calculations based on phenomenological interactions are performed starting from Brueckner G-matrix elements "fine tuned" to reproduce some specific set of spectroscopic data.

Results from phenomenological H eff 's
We test different phenomenological H eff 's. All these interactions have been derived modifying the matrix elements of a G-matrix so as to reproduce a chosen set of spectroscopic properties of some nuclei belonging to the mass region of interest. With this procedure one can end up with results that provide similar descriptions of the nuclei under consideration, nevertheless the phenomenological TBMEs are quite different each other.
It is worth stressing out that the calculated M 0ν s , reported in this section, are obtained using free value of the axial coupling constant g A without any quenching factor.
The double-magic nucleus 48 Ca is the lightest emitter investigated in regular ββ decay searches. The SM calculation for M 0ν is obtained using the model space spanned by four neutron and proton single-particle orbitals 0f 7/2 , 1p 3/2 , 1p 1/2 , and 0f 5/2 . It is worth mentioning that the regular ββ decay of 48 Ca is a paradigm for shell-model calculations. Because within the p f model space all spin-orbit partners are present, the Ikeda sum rule is satisfied [77].
Several phenomenological SM effective interactions have been developed to describe p f -shell nuclei. Among these are the GXPF1 [78], GXPF1A [79], KB3 [80], KB3G [81], and FPD6 [82] interactions. In Table 2, we compare the most recent results for the M 0ν of 48 Ca obtained using the GXPF1A [71] and KB3G interactions [83]. For the medium-mass emitters 76 Ge and 82 Se, the calculations adopt the valence space with the four neutron and proton single-particle orbitals 0f 5/2 , 1p 3/2 , 1p 1/2 , 0g 9/2 outside doubly-magic 56 Ni, as for instance in Refs. [84] and [83], where the effective interactions GCN2850 [38] and JUN45 [85] have been employed. These results are given in Tables 3 and 4. Finally, in Tables 5 and 6 we report and compare the calculated M 0ν for 130 Te and 136 Xe. These are based on two different effective interactions, namely the SVD [86] and the GCN5082 [38] defined in the jj55 valence space spanned by the neutron and proton orbitals 0g 7/2 , 1d 5/2 , 1d 3/2 , 2s 1/2 , and 0h 11/2 . Results with the SVD and GCN5082 interactions are taken, respectively, from Refs. [40] and [83]. Table 5. Same as Table 2, but for 130 Te. Results are from Refs. [40,83]. From the results shown above, it can be inferred that the effect associated with using different SRCs does not exceed 10%, while different effective interactions can provide results differing up to 50%. The results reported in this Section are based on the closure approximation. As discussed above, Senkov and coworkers have shown in a series of papers [71,84,87] that in going beyond this approximation, the 0νββ decay M 0ν becomes about 10% larger.

SVD (a) SVD (b) GCN5082 (a) GCN5082 (b)
We recall that the results reported in Tables 2-6 are obtained without quenching the axial coupling constant. However, the calculations based on the empirical SM Hamiltonians so far considered need a quenching factor q different from 1 to reproduce the experimental values of the nuclear matrix elements of the corresponding 2νββ-decays M 2ν s. This can be appreciated in Table 7 where we list the M 2ν s calculated with the empirical effective Hamiltonians GXPF1A, KB3G, JUN45, GCN2850 and GCN5082, and compare them with the experimental data. In these calculations, that are performed employing the Lanczos strength-function method [61], it has been used the unquenched value of g A (or equivalently a quenching factor q = 1), and, as expected, the theory is systematically overpredicting the experimental data. Table 7. M 2ν s for 48 Ca, 76 Ge, 82 Se, 130 Te, and 136 Xe 2νββ-decay calculated with GXPF1A, KB3G, JUN45, GCN2850, and GCN5082 interactions and compared with data [88], there are no published results based on the SVD interaction. These calculations use a quenching factor q = 1. The values of M 2ν s are reported in MeV −1 .

Results from realistic H eff 's
In the realistic SM (RSM), H eff is constructed from realistic V NN potentials. This is achieved via a similarity transformation utilized to constrain both the SM Hamiltonian and the SM transition operators. More details on this procedure can be found in the papers by B. H. Brandow [91], T. T. S. Kuo and coworkers [92,93], and K. Suzuki and S. Y. Lee [94,95]. Perturbative and non-perturbative derivations of H eff have been most recently reviewed in Ref. [96] and Ref. [97], respectively. The derivation of effective SM decay operators carried out consistently with H eff is discussed in Refs. [98,99]. Fundamental contributions to the field have been made by I. S. Towner, who has extensively investigated the role of many-body correlations induced by the truncation of the Hilbert space, especially for spin-and spin-isospin-dependent one-body decay operators [43,100].
The first calculations of M 0ν starting from realistic V NN potentials and associated effective SM decay operators, were made by Kuo and coworkers in the middle of 1980's for 48 Ca [64]. In that reference, H eff and the associated transition operators were based on the Paris and the Reid potentials [101,102]. The short-range repulsive behavior was renormalized by calculating the corresponding Brueckner reaction matrices [103]. Many-body perturbation theory was then implemented to derive both the TBMEs of H eff and the effective 0νββ-decay operator. The effect of the SRC was embedded in the defect wave function [104], consistently with the renormalization procedure from the Paris and Reid potentials. Finally, the authors calculated the half lives of 48 Ca 0νββ-decay, for both light-and heavy-neutrino exchange, as a function of the neutrino effective mass.
In more recently, J. D. Holt and J. Engel calculated effective SM operators Θ eff 's from modern chiral effective field theory V NN potentials. In particular, they started from the chiral V NN potential by Entem and Machleidt [105] and cured its perturbative behavior using the V low-k procedure [70]. The Θ eff was expanded up to third order in perturbation theory and used to calculate M 0ν for 76 Ge, 82 Se [106], and 48 Ca [107]. The effects of SRC was included via an effective Jastrow function obtained from Brueckner-theory calculations [26]. For 76 Ge and 82 Se decays, the authors employed the empirical GCN2850 [38] and JUN45 [85] SM interactions, respectively, and for the 0νββ decay of 48 Ca they used the GXPF1A H eff [79]. The values obtained by Holt and Engel in the light-neutrino exchange channel are M 0ν =1.30 for 48 Ca , M 0ν =3.77 for 76 Ge, and M 0ν =3.62 for 82 Se [106,107].
In Ref. [106], the authors calculated also the 76 Ge 2νββ matrix element. The calculation used the closure approximation. However, as we discussed in Sessions 2.3 and 2.4, this approximation is not robust when applied to study 2νββ processes where the values of momentum transfer are small. In fact, the authors obtain a result for 76 Ge M 2ν GT that is about two times larger than the one calculated with the Lanczos strength-function method [41,108], and about 5 times larger than the experimental value [88].
RSM calculations based on the high-precision CD-Bonn NN potential [109] have been recently carried out in Ref. [69], where the repulsive high-momentum components have been integrated out through the V low-k technique with "hard cutoff" Λ = 2.6 fm −1 [70]. The M 0ν 's have been calculated within the SM using H eff , and effective decay operators Θ eff 's up to the third order in perturbation theory. Two-body matrix elements entering the 0νββ-decay operator have been renormalized consistently within the V low-k to account for short-range correlations (see Section 2.2 for more details), and Pauli-principle violations in the effective SM operator. It should be pointed out that calculations of systems with a number n of valence nucleons require the derivation of n-body effective operators, that take into account the evolution of the number of valence particle in the model space P. The correlation between P-space configurations and those belonging to Q space is affected by the filling of the model-space orbitals, and in a perturbative expansion of SM operators this is considered by way of n-body diagrams. This is called the "Pauli-blocking effect" and calculations in Ref. [69], where all SM parameters have been consistently derived from the realistic V NN potential without any empirical adjustments, take it into consideration by including the contribution of three-body correlation diagrams to derive Θ eff .
The results for 0νββ decay in the light-neutrino exchange channel of 48 Ca, 76 Ge, 82 Se, 130 Te, and 136 Xe are reported in Table 8. It is worth pointing out that this approach to SM calculations has been successfully tested on energy spectra, electromagnetic transition strengths, GT strength distributions, and nuclear matrix elements for the two-neutrino ββ decay [110,111], without resorting to effective proton/neutron charges and gyromagnetic factors, or quenching of g A . In Table 9 we report the the calculated values of M 2ν from Ref. [111] and compare them with experimental data [88]. Table 9. Same as in Table 7, but the calculated values [111] are obtained employing effective Hamiltonians and decay operators derived starting from the CD-Bonn realistic potential (see text for details) and compared with experiment [88]. 48  A few comments are now in order. As we already pointed out in the Introduction, SM calculations overestimate M 2ν and Gamow-Teller transition strengths. To remedy to this deficiency one introduces a quenching factor q that is multiplied to g A to reduce the values of the calculated matrix elements. This factor depends on i) the mass region of the nuclei involved in the decay process; and ii) the dimension of the model space used in the calculation. The quenching factor has on average the empirical value q ≈ 0.7 [45]).
The quenching factor accounts for missing correlations and missing many-body effects in the transition operators. The truncation of the full Hilbert space to the reduced SM space has the effect of excluding all correlations between the model-space configurations and the configurations belonging to either the doubly-closed core or the shells placed in energies above the SM space. In addition, SM calculations are based on the single-nucleon paradigm for the transition operators. However, two-body electroweak currents [47][48][49][50][51][52][53][54][112][113][114][115][116][117][118][119][120] are found to play a role in several electroweak observables. These involve the exchange of mesons and nucleonic excitations.
I. S. Towner extensively studied how to construct effective β-decay operators that account for the degrees of freedom that are not explicitly included in the model space (see Ref. [43]). This has been more recently investigated in Refs. [110,111]). The results that are reported in Table 9 demonstrate a satisfactory agreement with the data can be obtained without resorting to quenching factors q if one employs effective GT operators within the SM.
Moreover, in Ref. [69] the authors showed that, the renormalization procedure implemented to account for missing configurations plays a marginal role in the calculated M 0ν s, while it is relevant in 2νββ-decay induced by the one-body GT operator. This evidences that the mechanisms which underlies the renormalization of the one-body single-β and the two-body 0νββ decay operators are different.
In closing, we address the role of many-body electroweak currents in the redefinition of single-β-decay and 0νββ-decay operators. Within the shell model, one can use nuclear potentials derived within chiral perturbation theory [121,122], and include also the contributions of chiral two-body electroweak currents. Studies along these lines have been recently carried out in Refs. [54,123,124], where the authors find significant contributions from two-body axial currents in β-decay. Concurrently, a study reported in Ref. [125] argues that many-body electroweak currents should play a negligible role in standard GT transitions (namely, β-and 2νββ decays) due to the "chiral filter" mechanism [125]. The chiral filter mechanism may be no longer valid for 0νββ decay since the transferred momentum is ∼ 100 MeV, which will require further investigations to fully understand the role of many-body currents in SM calculations of M 0ν s.

Comparison between SM calculations
In Fig. 2 [83]. Green diamonds correspond to the calculations by Horoi and coworkers [40,71,84]. Results by Holt and Engel are indicated by the red triangles [106,107]. Black dots indicate results by some of the present authors reported in Ref. [69].
The scale on the y-axis is consistent with the one employed in Fig. 5 of Ref. [17]. We stress that our criteria rule out, for the sake of consistency, results from SM calculations where alternative approaches have been followed. For example, we recall that Horoi and coworkers have extensively performed calculations beyond the closure approximation [71,84,87]. In particular, as already mentioned several times, results for the M 0ν of 48 Ca calculated with and without the closure approximation differ by ∼ 10%. Likewise, we neglect the results of large-scale SM calculations, where model spaces larger than a single major shell are used. This is the case, for instance, of the work by Horoi and Brown [39] and by Iwata and coworkers [126]. In the former, it was shown how the enlargement of the standard g 9/2 dsh 9/2 model space by the inclusion of the spin-orbit partners of g 9/2 and h 9/2 orbitals leads to a 10-30% reduction of the calculated M 0ν for the 136 Xe emitter. In the latter, the authors reported the results for the M 0ν of 48 Ca based on large-scale shell-model calculations including two harmonic oscillator shells (sd and p f shells). They found that M 0ν is enhanced by about 30% with respect to p f -shell calculations when excitations up to 2hω are explicitly included.
The spread among different results is rather narrow, except from the 130 Te → 130 Xe 0νββ decay since the result in Ref. [40] are more than 30% larger than those in Refs. [69,83]. We observe that the results from Refs. [69,83] are close each other, and the M 0ν 's calculated by Holt and Engel [106,107] are consistently larger than those from other SM calculations. The computational methods reported in the figure use substantially different SM effective Hamiltonians, yet they all lead to equally satisfactory results for a large amount of structure data. This leads us to argue that the SM approach is reliable for the study of 0νββ decay.
Finally, comparing Fig. 2 with the compilation of results reported in Fig. 5 in Ref. [17], we confirm that also these more recent SM calculations provide values that are smaller than those obtained with other nuclear models such as the Interacting Boson Model, the QRPA, and Energy Density Functional methods. Since the advantage of the nuclear shell model with respect to other approaches is to include a larger number of nuclear correlations, one may argue that by enlarging the dimension of the Hilbert space of nuclear configurations it should be expected a reduction in magnitude of the predicted values of M 0ν 's, as it is found e.g., in Ref. [39].

Conclusions
In this paper, we have briefly reviewed the present status of SM calculations of 0νββ-decay nuclear matrix elements. More precisely, we have focused our attention to the 48 Ca→ 48 Ti, 76 Ge→ 76 Se, 82 Se→ 82 Kr, 130 Te→ 130 Xe, and 136 Xe→ 136 Ba decays. These transitions are relevant to the current and planned experimental programs. We have considered calculations performed with phenomenological SM Hamiltonians, as well as studies where SP energies, two-body matrix elements, and effective decay operators have been derived from realistic V NN potentials. For completeness, important aspects that characterize the calculation of M 0ν , such as the role of short-range correlations and the closure approximation, have been briefly discussed. These and other approximations, such as the size of the model space, may affect the results by 30%.
We showed how different SM calculations, notwithstanding the diversity of the effective Hamiltonians that are employed to calculate the nuclear wave functions, exhibit a rather narrow spread among the predictions of the nuclear matrix elements, making the SM a solid and reliable framework for M 0ν calculations.