Comparison between Variational Monte Carlo and Shell Model Calculations of Neutrinoless Double Beta Decay Matrix Elements in Light Nuclei

Benchmark comparisons between many-body methods are performed to assess the ingredients necessary for an accurate calculation of neutrinoless double beta decay matrix elements. Shell model and variational Monte Carlo (VMC) calculations are carried out for $A=10$ and $12$ nuclei. Different variational wavefunctions are used to evaluate the uncertainties in the {\it ab initio} calculations, finding fairly small differences between the VMC double beta decay matrix elements. For shell model calculations, the role of model space truncation, radial wavefunction choices, and short-range correlation are investigated and all found to be important. Based on the detailed comparisons between the VMC and shell model approaches, we conclude that accurate descriptions of neutrinoless double beta decay matrix elements require a proper treatment of both long-range and short-range correlations.


I. INTRODUCTION
Neutrinoless double beta decay (0νββ) is a process in which two neutrons in a nucleus decay into two protons, with the emission of two electrons and no neutrinos, thus violating lepton number (L) by two units.The observation of 0νββ would imply that neutrinos are Majorana particles [1], shed light on the mechanism of neutrino mass generation, and give insight into leptogenesis scenarios for the generation of the matter-antimatter asymmetry in the universe [2].For these reasons, 0νββ is the subject of intense experimental research programs [3][4][5][6][7][8][9][10][11][12][13][14].Current experimental limits are already very stringent, constraining the half-lives of 76 Ge, 130 Te and 136 Xe to be larger than 8 • 10 25 yr [13], 1.5 • 10 25 yr [14] and 1.1 • 10 26 yr [7], respectively.The next-generation ton-scale experiments will improve these limits by one or two orders of magnitude.
0νββ rates depend not only on nuclear properties, but also on unknown fundamental lepton number violating (LNV) parameters, such as the Majorana masses of light neutrinos.Extracting the values of these parameters from experiments requires the evaluation of nuclear matrix elements (NMEs) of 0νββ transition operators.Isotopes of experimental interest for 0νββ searches, e.g. 48Ca, 76 Ge, 82 Se, 124 Sn, 128 Te, 130 Te and 136 Xe, are medium and heavy open-shell nuclei with very complex nuclear structure.In addition, for these nuclei, one is forced by current computational limitations to utilize approximate methods to solve the nuclear many-body problem, and to work in a truncated model spaces where correlations and many-body terms in both the nuclear interactions and currents may be insufficient or neglected.Consequently, different theoretical models can give systematically different results.For example, the nuclear matrix elements of 48 Ca have been calculated using the shell model [15][16][17][18], energy density functionals [19], the quasiparticle random-phase approximation (QRPA) [20], and the interacting boson model (IBM) [21], leading to results that differ by a factor of two or three.Similar variations are observed for other 0νββ candidates [22].
Recently, a set of ab initio variational Monte Carlo (VMC) calculations of 0νββ decay matrix elements in A = 6−12 nuclei has been reported by some of the present authors [23].Within the ab initio VMC framework, the many-body problem is solved for nuclei up to A = 12, with a nuclear Hamiltonian consisting of two-and three-body interactions, namely the Argonne-v18 (AV18) and Illinois-7 (IL7), respectively, and associated electroweak many-body currents.While the 0νββ transitions studied in Ref. [23] are not directly relevant from an experimental point of view because the masses studied are too low, they provide a useful reference for more approximated nuclear methods.The ab initio framework used in the VMC approach accurately explains, qualitatively and quantitatively, the observed electroweak properties of light nuclei [24][25][26] over a broad range of momentum transfers, so that VMC results provide an important benchmark to test other many-body methods that can be extended to the heavy nuclei of experimental interest.
The goal of the present work is to benchmark shell model calculations of the relevant 0νββ NMEs in light nuclei to the aforementioned VMC results.We wish to examine the model dependence and uncertainties involved in shell model approaches and to identify the degree of sophistication that needs to be included in such a calculational approach to 0νββ.In general, shell model calculations involve a number of choices, including the size of the model space, the effective nucleon-nucleon interaction, the radial wave functions, and the short-range correlation (SRC) functions used.In the current work, we concentrate on nuclei of mass 10 and 12, and we examine the role played by these different nuclear structure inputs in determining the predicted 0νββ NMEs.
The paper is organized as follows.In Section II, we present the two-body transition operators that mediate 0νββ, and the shell model framework that is used to evaluate matrix elements of these operators.In Section III, we describe two sets of shell model calculations, and present a detailed analysis of uncertainties arising from different choices of the basis, radial functions, and SRC functions.Finally, in Section IV we discuss our results and present our conclusions.

A. 0νββ operators
We consider 0νββ transitions induced by a Majorana mass term for light, left-handed neutrinos where m ββ = U 2 ei m i , m i are the masses of the neutrino mass eigenstates, and U ei are elements of the Pontecorvo-Maki-Nakagawa-Sato (PMNS) matrix.C = iγ 2 γ 0 denotes the charge conjugation matrix.The operator in Eq. ( 1) arises from the Weinberg operator [27] after electroweak symmetry breaking, and the SU (2) L invariance of the Standard Model implies that m i ∼ v 2 /Λ, where v is the Higgs vacuum expectation value and Λ is the high-energy scale at which LNV arises.Eq. ( 1) is the first term in an expansion in v/Λ, and, in general, 0νββ will receive contributions from operators of dimension higher than five [28][29][30][31][32][33][34].We limit ourselves to study nuclear matrix elements of the 0νββ transition operator induced by m ββ , since its structure is rich enough to capture several features of non-standard mechanisms.
The 0νββ transition operator, or "neutrino potential", induced by m ββ has been derived in several papers [35][36][37].In general, several operators M β α contribute, with α ∈ {F, GT, T }, and β ∈ {ν, AA, AP, P P, M M } labeling the various components of the Fermi(F), Gamow-Teller(GT) and tensor(T) matrix elements of the vector (ν), axialaxial(AA), axial-pseudoscalar (AP), pseudoscalar-pseudoscalar (PP), and magnetic-magnetic (MM) operators.The leading operators are, where R A = 1.2A 1/3 fm is the nuclear radius.The form of the neutrino potentials V β α (r 12 ) in coordinate space is given in Ref. [23], where the corresponding potential in momentum space is also listed.The functions V β α in Eq. ( 2) have different radial dependencies.In particular, V ν F and V AA GT , which give the largest contributions to the NMEs, have a Coulombic 1/r dependence, up to corrections from the axial and vector form factors.The AP and P P components are induced by the induced pseudoscalar form factor, which is dominated by the pion pole, and thus have pion range.The M M components are generated by neutrinos coupling to the weak magnetic form factor, and have short-range.Finally, it was pointed out in Ref. [38] that, in chiral effective field theory (χEFT), the 0νββ transition operator should be supplemented by a leading-order short-range component, in addition to M M , with unknown coefficient.

B. Variational Monte Carlo Method
In this work, VMC calculations used the same wave functions as constructed in Ref. [23].Here, we summarize the calculations scheme adopted in that reference.
The evaluation of the matrix elements defined in Eq. ( 2) is carried out using VMC computational algorithms [24].The VMC wave function Ψ(J π ; T, T z )-where J π and T are the spin-parity and isospin of the state-is constructed from products of two-and three-body correlation operators acting on an antisymmetric single-particle state of the appropriate quantum numbers.The correlation operators are designed to reflect the influence of the two-and threebody nuclear interactions at short distances, while appropriate boundary conditions are imposed at long range [39,40].
The Ψ(J π ; T, T z ) has embedded variational parameters that are adjusted to minimize the expectation value which is evaluated by Metropolis Monte Carlo integration [41].In the equation above, E 0 is the exact lowest eigenvalue of the nuclear Hamiltonian H for the specified quantum numbers.The many-body Hamiltonian is given by where K i is the non-relativistic kinetic energy of nucleon i and v ij and V ijk are, respectively, the Argonne v 18 (AV18) [42] two-body potential and the Illinois-7 (IL7) [43] three-nucleon interaction.The AV18+IL7 model reproduces the experimental binding energies, charge radii, electroweak transitions and responses of A = 3-12 systems in numerically exact calculations based on Green's function Monte Carlo (GFMC) methods [24][25][26]44].
A good variational wave function, that serves as the starting point of GFMC calculations, can be constructed with The Jastrow wave function Ψ J is fully antisymmetric, translationally invariant, has the (J π ; T, T z ) quantum numbers of the state of interest, and includes a product over pairs of a central correlation function f (r ij ) that is small at short distances, peaks around 1 fm, and decays exponentially at long range [45].The U ij and Ũijk are the two-and three-body correlation operators, and S is a symmetrization operator.The two-body correlation operators [24,45] can be schematically written as where are the main static operators that appear in the two-nucleon potential and the f p are functions of the interparticle distance r ij generated by the solution of a set of coupled differential equations containing the bare two-nucleon potential with asymptotically-confined boundary conditions [24].
The results presented below for A = 10 nuclei use the VMC wave functions that serve as starting trial functions for the GFMC calculations in Ref. [24].For A = 12 systems, we use both shell-model-like wave functions (denoted in the figures by "VMC-1") and clusterized wave functions (denoted by "VMC-2") which were originally used in the calculations reported in Ref. [46].The difference in these two types of variational wave functions is described below.
For A > 4 nuclei, the Jastrow wave function Ψ J of Eq. ( 5) is written as a sum over different LS-space symmetry combinations 2S+1 L J [n] (where [n] denotes the spatial symmetry Young diagram), each of which is an antisymmetric sum over partitions of A nucleons into four s-shell nucleons and A − 4 p-shell nucleons.For the shell-model-like wave function, each partition has a product over all pairs of a set of three central correlation functions, f xy (r ij ), where xy = ss indicates both nucleons are in the s-shell, xy = pp indicates both are in the p-shell, and xy = sp that one is in each shell.Each of the f xy has a different radial dependence with variational parameters to be optimized, but all p-shell particles are treated equally.Thus for a 12-body nucleus, there is a product of 6 f ss , 32 f sp , and 28 f pp functions.In addition, the parametrization of the f xy are allowed to be different for different 2S+1 L J [n] components.For the 12 C ground state we use 1 S 0 [444] and 3 P 0 [4431] spatial symmetry combinations, with coefficients determined in a 2 × 2 energy diagonalization.For 12 Be we use 1 S 0 [4422] and 3 P 0 [4332] combinations.
For the cluster-type wave functions, we allow for the formation of clusters in the p-shell by further partitioning the A − 4 nucleons into subgroups and having multiple types of f pp and f sp .For the 12 C 1 S 0 [444] ground state there are two f pp correlation functions, used according to whether both nucleons are in the same or different [4] clusters; effectively this builds a triple-alpha structure which is a major part of the 12 C ground state.For the 12 Be 1 S 0 [4422] ground state there are three subgroups in the p-shell, effectively one alpha-like [4] cluster and two dineutron-like [2] clusters, with three f pp and two f sp correlation functions.These cluster-type wave functions are more sophisticated in construction, but the increase in the number of variational parameters to optimize is a burden, and only these highest spatial symmetry states have been used so far.
At present, the shell-model-like VMC-1 wave function for 12 C with two spatial symmetry components gives a slightly better energy, while getting a charge radius closer to experiment than the cluster-like VMC-2 wave function with one spatial symmetry.The VMC-1 and VMC-2 wave functions for 12 Be give comparable charge radii close to experiment, but the VMC-2 allows the neutron distribution to be more diffuse with a slightly better energy.
In what follows, we calculate matrix elements defined in Eq. ( 2), and their associated transition distributions in r-space, C α,β (r), and q-space, Cα,β (q) defined as where ρ α β (r) is the transition density associated with the transition operator O α β (r).

C. NMEs within the shell model framework
Within the shell model, the matrix elements of the M α β operators between many-body wave functions | i > and | f > can be calculated as a sum of products of two-body transition densities (TBTDs) between many-body states and antisymmetrized two-body matrix elements TBMEs between two-particle states, as Here the labels k stands for the set of spherical quantum numbers (n, l, j) describing the single nucleons making up the two-particle wave functions, and TBTD are matrix elements between the many-body wave functions given by, where A + is a two-particle creation operator of rank J 0 , as

III. CALCULATIONS AND DISCUSSIONS
We focus on the 0νββ NMEs for the ∆T = 0 10 Be→ 10 C and ∆T = 2 12 Be→ 12 C transitions.The T = 0 transition is between isobaric analog states, a situation that is never realized in 0νββ experiments, but is nonetheless very useful for comparisons between models.In particular, as discussed below, the radial contributions to the A=10 0νββ matrix elements involve no nodes, in contrast to the A=12 case.As in Eq. ( 8), the radial contribution to the matrix elements is defined as C(r), where ∞ 0 C(r)dr = M β α .Our shell model calculations are carried out using the code Oxbash [47].The PSDMWK shell model Hamiltonian [48,49] is used for the p-plus sd-shell (psd) model space, with up to 4p4h excitations.Lawson's prescription is adopted for treating center of mass (C.M.) excitations [50].The same NMEs are also calculated by VMC, for which the Argonne v 18 two-nucleon potential plus Illinois-7 three-nucleon interaction are used.As discussed in Sec.II, two different sets of variational wave functions are used in the VMC calculations, namely the shell-model-like (VMC-1) and cluster-model-like (VMC-2).

A. The wave function normalization
We begin our study with a comparison of the wave function normalization densities  Talmi-Brody-Moshinsky brackets [51,52].However, H.O. wave functions exhibit the wrong asymptotic behavior at large interparticle distance, r, which affects the calculation of 0νββ NMEs.Thus, we examine the effect of using more realistic Woods-Saxon (WS) radial wave functions.In the case of WS wave functions, we take the p-shell neutron and proton asymptotic wave function behavior to be determined by the corresponding experimental separation energies, while the unbound sd-shell particles are assumed to be bound by 0.05 MeV.For the calculations of two-body matrix elements, we expand the WS wave functions into a H.O. basis, and then apply the Talmi-Moshinsky method.
As seen in both panels of Fig. 1, the extended shell model space results in a higher first (and lower second) peak in C(r).Enlarging the model space introduces more correlations, making the nuclei more bound, and the wave functions more concentrated at short distances.The more realistic WS wave functions provide a better description at large r.We note that the normalization of C(r) is unity for the A=10 T = 0 isobaric analog transition, and zero for ∆T = 2 A=12 transition.The position of the nodes in the ∆T = 2 C(r) function has an important effect on the predictions for the 0νββ NMEs.
As seen in panel (b) of Fig. 1, the first node appears at about 2 fm in both the shell model and VMC calculations.However, the smaller p-shell model space has an extra node at around 5 fm, in contrast to the predictions of the larger model spaces.Neither VMC calculations predict this second node.Thus, we conclude that the size of the model space used is crucial for shell model estimates of 0νββ NMEs.

B. The radial distribution of NMEs
We next consider the matrix elements of V ν F and V AA GT , which have a simple Coulombic 1/r dependence and give the largest contributions to the NMEs.Detailed comparisons between the shell model and VMC calculations are shown in Fig. 2. Apart from the differences in magnitude and sign, the distribution C(r) for V ν F and V AA GT have similar shapes.Because of the Coulombic 1/r dependence, the first peak of the distributions is much larger than the second peak, and the functions die off rapidly at large distance r.As with the normalization functions, the shell model distributions from the larger psd model spaces show higher peaks than those from the smaller p-shell model spaces, when the same radial wave functions are used.In general, the predicted shape of the functions C(r) is poor in the case of the small pure p-shell calculations, and we again conclude that small shell model spaces are unlikely to provide accurate 0νββ predictions.When a larger model space is used in combination with the more realistic WS wave functions, the shell model predictions are in reasonable agreement with the VMC results, except at very small r < 1 fm.
) r ( f m ) nuclei.F-ν and GT-AA matrix element contributions are labeled as "F" and "GT ", respectively.The plots show the radial distributions C(r) determining the 0νββ matrix elements, where The shell model results are for the lowest order 0hω p-shell calculations and psd-shell calculations with up to four particle excitations, and for two choices of radial wave functions, H.O. and WS, without short-range correlations.VMC results are the full calculations and include statistic error bars.VMC results with shell-model-like wave functions are labeled as "VMC-1", and those with cluster-like wave functions are labeled as "VMC-2".The radius factor RA of Eq. ( 2) is not included in the curves.

C. The short-range correlation functions
Short-range correlations (SRC) arising from the repulsive hard core of the nucleon-nucleon interaction are absent from the radial wave functions used in shell model calculations.A standard approximate method for correcting this is to introduce a Jastrow-like correlation function, f (r ≡| r 1 − r 2 |), of the form, The parameters a, b, c have been determined using different assumptions.The Miller-Spencer (M.S.) SRC, modeled in 1976 [53], is widely used in the literature.Newer sets of SRC have been obtained from 1 S 0 (n = 0) correlated two-body wave function derived by the coupled-cluster method (CCM) [54].More recently, through the study of paired density distributions in 16 O and 40 Ca by Cluster Variational Monte Carlo (CVMC) [55,56], isospin dependent correlation functions have been suggested [57].Yet another set of correlation functions was obtained in Ref. [58] from nuclear matter variational calculations, and we also study the effects of using these functions.In Figs. 3 and 4, two-body correlation functions from Ref. [58] with and without three-body and higher order correlations are labeled with "Fab+abc" and "Fab", respectively.In Fig. 3, we show a comparison between the different Jastrow correlation functions which we used in calculating the shell model NMEs.This allows us to assess the sensitivity of our shell model results with respect to variations in 0 .0 0 .5 1 .0 1 .5 0 .0 0 . 2 0 .4 0 .6 0 .8 The square of the Jastrow-like correlation functions (Eq.( 13)).The different correlation functions are described in the text.The points correspond to the ratio of the VMC to shell-model (with psd model space and WS radial wave functions) normalization functions for A=10 and 12.
the SRCs.The points in Fig. 3 denote the ratio of the normalization density C(r), defined in Eq. ( 12), computed in the VMC and in the shell-model for A=10 and 12.This ratio might be regarded as a rough estimate of short-range dynamics present in VMC versus shell model calculations.
As seen in the figure, most of short-range correlations overshoot unit at intermediate distance, and exhibit a correlation "hole" at short distance.Such overshoot is needed to preserve the wave function normalizations, as has been studied in Ref. [59].The Miller-Spencer SRC function peaks at about 1.5 fm, while more modern functions, CCM SRC and CVMC SRC, peak at 1.0 fm.The CVMC SRC exhibits a more pronounced peak and hole relative to the CCM SRC (Av18).As discussed in the next section, the difference in NME predictions obtained using the different SRC functions reflects a general uncertainty in the shell model calculations.

D. Results
In Fig. 4, shell model matrix elements using the different SRC functions are shown.In these calculations, psd shells with WS wave functions are used.As expected, the degree of suppression at small r of the function C(r) is dependent on the SRC used.The resulting NMEs are given in Table .I. The VMC calculations with the two different type of wave functions give very similar results.However, within a shell model framework, several important ingredients must be included for realistic predictions.These include the size of the model space, the choice of radial wave functions and the inclusion of realistic SRC functions.
For A=10, the extension of shell model space results in larger (in absolute value) NMEs.For A=12, the extension of the model space results in smaller matrix elements, because of the cancellations from contributions above and below the node in C(r).The use of more realistic radial wave functions, WS versus HO, also reduces the predicted NMEs.The inclusion of a SRC function is important to further reduce the shell model predictions and to obtain realistic functions C(r).In Ref. [60], it was reported that for 48 Ca and 76 Ge, a 30 to 40 % reduction of NMEs arose from the inclusion of the M.S. SRC.In the present study of light nuclei, the M.S. SRC reduces the NMEs by about 20 %.Among more modern SRC functions, respecting the preservation of wave-function normalization [59], SRC fitted from CVMC also gives similar reduction of NMEs.

IV. SUMMARY AND CONCLUSION
In the current study, we compare 0νββ NMEs predictions in light nuclei from shell model calculations with VMC calculations.In all cases the bare operators are used.The VMC calculations agree well with numerically exact GMFC   We study two very different double beta-decay transitions, 10 Be→ 10 C and 12 Be→ 12 C.The most significant difference between these two transitions is that nuclear structure causes the C(r) radial contributions to the 0νββ matrix elements in A=10 to be nodeless, whereas the same C(r) distributions for A=12 involve nodes and hence cancellations in the matrix elements.
We examine the role of various model-dependent choices that go into shell model calculations of 0νββ matrix elements, including the the size of the model space, radial wave functions, and the SRC functions.The largest impact on the shell model predictions comes from the many-body correlations that are introduced by increasing the size of the model space.In both A=10 and A=12, the transition densities C(r) determining the matrix elements of M β α from the larger shell model calculation come significantly closer to the predictions of the VMC calculations than do the predictions from the smaller shell model spaces.
The change in the absolute magnitude of the M F and M GT 0νββ matrix elements with increasing shell model space is different for A=10 and A=12.In both cases the magnitude of the peaks in the C(r) distribution increase.In the case of the nodeless A=10 C(r) function, this translates into an overall increase in the absolute magnitude of the 0νββ matrix elements.In contrast, the change in the magnitude of the peaks, coupled with the shift in the position of the node, causes the A=12 cancellations from contributions above and below the node of C(r) to reduce (increase) the absolute magnitude of the M GT (M F ) matrix element by a factor of 2 (1.4).
The use of more realistic WS radial wave functions, that take the separation energies of the transferred neutrons and protons into account, reduce the magnitudes of the NMEs because of the reduced overlap between the neutron and proton wave functions.By using the WS wave functions and increased model space, the long-range behavior of C(r) distribution is found to be in excellent agreement with the VMC predictions.
The SRC functions affect the contributions to the matrix elements at short distance (less than 1.5 fm).We have examined several different SRC functions.In general, the introduction of SRCs moves the shell model predictions closer to the VMC calculations, reducing the NMEs.The magnitude of this reduction is dependent on the SRC function used.
While the present calculations in light nuclei cannot be used to make definitive statements about the validity of shell model calculations in medium and heavy nuclei, they do suggest some trends.First, the use of H.O. radial wave functions will likely lead to an overestimate of matrix elements.Second, and perhaps more importantly, limited size model space calculations could affect the magnitude of the predicted 0νββ matrix elements, particularly for calculations constrained to a single shell.Third, the inclusion of a SRC function is needed.The best choice for this function requires further study.

displayed in Fig. 1 .
For the shell model calculations, we show the results for the p shell and the more extended psd model space.In addition, we show shell model results using different choices of radial wave functions.The SRC functions are not included at this stage in the shell model wave functions.Harmonic oscillator (H.O.) radial wave functions are commonly used in shell model calculations because H.O. basis have the advantage of a simple well-defined method of separating into center-of-mass and relative coordinates via

FIG. 1 :
FIG. 1: (Color online) The normalization densities C(r) = f | a<b δ(r − r ab )τ + a τ + b |i , for (a) A=10, and (b) A=12.The 0hω model space (p-shell only) and extended model space of psd shell model calculation results are shown.The different choices of radial wave functions, H.O., and WS, are also shown.VMC results with shell-model-like wave functions are labeled as "VMC-1", and those with cluster-like wave functions are labeled as "VMC-2".

FIG. 2 :
FIG. 2: (Color online) The radial distributions associated with F-ν and GT-AA operators, for (a)-(b) A=10, and (c)-(d) A=12nuclei.F-ν and GT-AA matrix element contributions are labeled as "F" and "GT ", respectively.The plots show the radial distributions C(r) determining the 0νββ matrix elements, where

FIG. 4 :
FIG.4: (Color online) The same as in Fig.2, but for the shell model calculations of NMEs involving SRC.

TABLE I :
Results of dimensionless matrix elements F-ν and GT-AA matrix elements, labeled as "F" and "GT ".VMC results with different variational wave functions are labeled as "VMC-1" and "VMC-2".VMC statistical errors are also listed.Shell model results with different radial wave functions (H.O. and WS), model space (p and psd) and SRC are given correspondingly.The radius factor RA from eqs. (2) is included., reproducing nuclear radii and electroweak distributions.In this sense, the VMC calculations act as a benchmark for comparison to other models. calculations