Nuclear matrix elements for 0 νββ decays with light or heavy Majorana-neutrino exchange

We compute the nuclear matrix elements (NMEs) corresponding to the neutrinoless double beta (0 νββ ) decays of nuclei which attract current experimental interest. We concentrate on ground-state-to-ground-state decay transitions mediated by light (l-NMEs) or heavy (h-NMEs) Majorana neutrinos. The computations are done in realistic single-particle model spaces using the proton-neutron quasiparticle random-phase approximation (pnQRPA) with two-nucleon interactions based on the Bonn one-boson-exchange G matrix. Both the l-NMEs and the h-NMEs include the appropriate short-range correlations, nucleon form factors, and higher-order nucleonic weak currents. In addition, both types of NMEs are corrected for the isospin symmetry by the recently proposed method in which the particle-particle proton-neutron interaction parameter ( g pp ) is decomposed into isoscalar ( g T = 0 pp ) and isovector ( g T = 1 pp ) parts. A detailed analysis of the l-NMEs and the h-NMEs is performed to benchmark our computer code and to compare with other recent calculations which produce h-NMEs that are in tension with each other.


I. INTRODUCTION
The neutrinoless double beta (0νββ) decay is a traditional way to search for physics beyond the standard model. It can access the absolute neutrino mass and the fundamental nature of the neutrino, i.e., whether it is a Dirac or a Majorana particle (see, e.g., [1][2][3]). In this way it is complementary to the neutrino oscillation experiments which can investigate the neutrino mass differences and neutrino mixing amplitudes to high precision already. The information to be extracted from the 0νββ experiments is subject to uncertainties arising from the uncertainties in the related nuclear matrix elements (NMEs). Hence, it is extremely important to be able to compute these NMEs as accurately as possible.
The half-lives of the standard-model variant of the 0νββ decay, the two-neutrino double beta (2νββ) decay, have a strong dependence on the NMEs which contain virtual transitions through the 1 + intermediate states of the decay. The 0νββ decays, on the other hand, proceed by virtual transitions through intermediate states of all multipoles J π , J being the total angular momentum and π being the parity of the intermediate state. Both the 2νββ and 0νββ NMEs have been computed by a number of different models [1] but lately the computations of the 0νββ NMEs by the following models have drawn the attention of the double beta community: the quasiparticle random-phase approximation (QRPA), in its proton-neutron version (pnQRPA) (see [4] and references therein) and its renormalized extensions [5,6], the interacting shell model (ISM) [7], the (proton-neutron) interacting boson model (IBA-2) [8], the Gogny-based energydensity functional approach (EDF) [9], and the projected Hartree-Fock-Bogoliubov mean-field scheme (PHFB) [10]. An extensive comparison of the double beta properties of the aforementioned models is performed in [11].
The pnQRPA model is a mature and established model to compute the NMEs of various double-beta-decay processes [12]. However, it suffers from particular shortcomings that are associated with its limited number of many-particle configurations (in comparison with the ISM), its sensitivity to the effective proton-neutron interaction in the particle-particle channel of two-nucleon interactions (the "g pp problem"), and the isospin breaking induced by the quasiparticle approximation lurking behind the pnQRPA treatment of many-body correlations. The limited number of configurations in the pnQRPA might not be that severe shortcoming for the computation of the double beta decays since the centroids of the strong transitions through the intermediate multipole states can be handled in sufficient accuracy by the pnQRPA, in particular since it can accommodate a large set of single-particle basis states including all the relevant spin-orbit partner states [13,14]. The "g pp problem" is usually handled by the inspection of the measured single-beta-decay rates [15,16] or 2νββ decay rates [17][18][19][20]. Lately, the problem of isospin breaking has been tackled in [21] where an isospin-restoration scheme for the pnQRPA has been proposed.
Referring to the above-mentioned shortcomings of pn-QRPA, in the present work we solve the "g pp problem" by fixing the value of g pp such that the measured 2νββ half-lives are reproduced by the computed 2νββ Gamow-Teller NMEs, in the spirit of [17][18][19][20][21]. At the same time we perform the isospin-symmetry restoration by the scheme proposed in [21]. Furthermore, we take into account the appropriate shortrange nucleon-nucleon correlations [22], and contributions arising from the induced currents and the finite nucleon size [23]. In [21] the NMEs corresponding to the exchange of light Majorana neutrinos (l-NMEs) were treated for the conservation of the isospin symmetry. Here we extend this work to the exchange of heavy Majorana neutrinos and treat the corresponding NMEs (h-NMEs) by the isospin-restoration recipe of [21]. The computation of the h-NMEs has been, and still seems to be, a problem. By this we refer to the very different and contradictory results obtained in [3], [23], [24], and [25]. In this work we will analyze and decompose the h-NMEs to enable clean comparison with other calculations with different methods.
This article is organized as follows: In Sec. II we give a brief introduction to the underlying formalism of the 2νββ and 0νββ decays. In Sec. III we discuss the determination of the model parameters, and display and discuss the obtained double-beta-decay results for the l-NMEs and h-NMEs. The final conclusions are drawn in Sec. IV.

II. COMPUTATIONAL FORMALISM
In this section the basic theoretical ingredients of the computations are reviewed. The double-beta-decay half-lives are defined in terms of the involved NMEs and the phase-space factors. The NMEs are given in terms of the single-particle matrix elements and the one-body transition densities. For the calculation of the phase-space factors we refer to the available literature at the appropriate places.
The overlap factor takes care of the matching of the corresponding states in the two sets of states based on the initial and final even-even reference nuclei. The amplitudes X and Y (X andȲ ) come from the pnQRPA calculation starting from the final (initial) nucleus of the double beta decay. In principle, the expression (2) should also contain the Fermi contribution, i.e., (2) should go over to where the Fermi NME (for pnQRPA) reads with Here the denominator and the overlap factor correspond to the expressions (7) and (8) with the 1 + intermediate states replaced by the 0 + intermediate states. However, we will force this contribution to be zero by the procedure described in Sec. III A. This is justified since in the case of isospin symmetry, obeyed by the nuclear forces to good extent, the Fermi contribution to the 2νββ NME (9) should vanish. Thus only the Gamow-Teller contribution is left and the final 2νββ NME will be the one of (2).

B. Neutrinoless double beta decays
Assuming that the light or heavy Majorana-neutrino exchange is the dominant mechanism over the other possible mechanisms, the 0νββ half-life for a ground-state-to-groundstate transition can be written as where G 0ν is a phase-space factor for the final-state leptons defined here without the axial-vector coupling constant g A . The quantity η x describes the physics beyond the standard model. For the light-neutrino exchange (sum over the mass eigenstates m j of light neutrinos) the η x factor denotes the effective mass and for the heavy-neutrino exchange (involving the mass eigenstates M j of heavy neutrinos) the η x factor stands for Here the amplitudes U l ej and U h ej are the components of the electron row of the neutrino-mixing matrix including both light and heavy neutrinos.
The nuclear matrix element M (0ν) in (13) is defined as where the double Fermi, Gamow-Teller, and tensor nuclear matrix elements are defined as where the tensor operator reads The summation over k in Eqs. (17), (18), and (19) runs over all the states of the intermediate odd-odd nucleus, r mn = |r m − r n | is the relative distance between the two decaying neutrons, labeled m and n, andr mn = (r m − r n )/r mn . The ground state of the initial even-even nucleus is denoted by 0 + i and the ground state of the final even-even nucleus is denoted by 0 + f , like in the two-neutrino case. The tensor matrix element M (0ν) T has been found to be small [18,19,27] but here we include it for completeness.
The neutrino potentials h K (r mn ,E k ), K = F,GT,T, are different for the light and heavy Majorana neutrinos and are given by the following two expressions: for the light neutrinos and for the heavy neutrinos. Here R A = 1.2A 1/3 fm is the nuclear radius, M i c 2 (M f c 2 ) is the ground-state mass energy of the initial (final) nucleus, E k is the energy of the nuclear state k of the intermediate nucleus and j λ is the spherical Bessel function with λ = 0 for the Fermi and Gamow-Teller NMEs and λ = 2 for the tensor NME. The factors h K (q 2 ) in (21) and (22) include the contributions arising from the induced currents and the finite nucleon size [23] and can be written as with The dipole form factors we define here without the vector and axial-vector coupling constants [they have already been included in (13) and (16)] aŝ where m π is the pion mass, μ p − μ n = 3.70 is the anomalous magnetic moment of the nucleon, and we take for the vector mass ( V c 2 ) 2 = 0.71 (GeV) 2 and for the axial mass ( A c 2 ) 2 = 1.19 (GeV) 2 .
The nuclear matrix elements can be written in the pnQRPA framework as where k 1 and k 2 label the different pnQRPA solutions for a given multipole J π . The operators O K inside the two-particle matrix element derive from (17), (18), and (19), and they can be written as where S T 12 is the tensor operator (20) and r = |r 1 − r 2 | is the distance between the participating nucleons. The energy E k is the average of the kth eigenvalues of the pnQRPA calculations based on the initial and final nuclei of the decay and the overlap factor in (38) is the one of (8). The factor f CD (r) takes into account the nucleon-nucleon short-range 024613-3 TABLE I. The key Hamiltonian parameters used in our present calculations. The first column gives the decaying nucleus and the second column gives the average of the proton and neutron pairing parameters in the initial and final nuclei. The third column lists the values of the isovector particle-particle parameter, and the last two columns the values of the isoscalar particle-particle parameters for the two representative values of the axial-vector coupling constant. correlations (SRC) [22,28] and in this work we use the CD-Bonn parametrization [29] of the SRC with

C. Transition densities in the pnQRPA
The pnQRPA states of the intermediate nucleus are written as where |QRPA is the QRPA vacuum. The operator a † p (a † n ) creates a proton (neutron) quasiparticle in the orbital p (n). The sum runs over all proton-neutron configurations in the chosen valence space. In the case 0 + f = 0 + gs the form (43) of the pnQRPA state leads to the transition densities where v (v) and u (ū) correspond to the BCS occupation and vacancy amplitudes of the initial (final) even-even nucleus.
The amplitudes X and Y (X andȲ ) come from the pnQRPA calculation starting from the initial (final) nucleus of the double beta decay.

III. RESULTS AND DISCUSSION
In this section we present and discuss the results of the calculations. The presentation of the material serves the purpose of a detailed comparison with other presently available and future calculations by different methods.

A. Determination of model parameters
In the present calculations we adopt the single-particle bases used in the 0νββ calculations of [18,19], i.e., the orbitals 1p − 0f − 2s − 1d − 0g − 0h 11/2 for the A = 76,82 systems, the orbitals 1p − 0f − 2s − 1d − 0g − 0h for the A = 96,100 systems, and the orbitals 1p − 0f − 2s − 1d − 0g − 2p − 1f − 0h for the A = 116,128,130,136 systems. The same orbitals are used for both neutrons and protons. In  The orbitals left out from the active single-particle space have BCS occupation probabilities which are essentially 1 (core states) or 0 (high-lying states). As long as the number of included single-particle states is reasonable, the adjustment of the pairing and pnQRPA parameters in each given singleparticle space guarantees that the produced 0νββ NMEs are stable and essentially of the same magnitude irrespective of the number of included single-particle states [17]. Hence, it is expected that the presently produced l-NMEs and h-NMEs are stable against increase in the single-particle model space.
The single-particle energies were obtained from the Coulomb-corrected Woods-Saxon potential with the parametrization of Ref. [30], optimized for nuclei close to the β stability line. This is a justified choice since the double-beta-decaying nuclei lie always close to the bottom of the valley of beta stability. Small adjustments of the proton and/or neutron single-particle energies in the vicinity of the respective Fermi surfaces were done in the way done in [18,19] and thus our present results are directly comparable with those of [18,19] in this respect. The adjustments were done to better reproduce the low-lying spectra of the neighboring odd-mass nuclei. This also improves the computed energy spectra of the   [21] in the third last column, the second last column showing the ratios between our present results and those of Ref. [21]. The last column shows the relative difference(|M (0ν) [21] −M (0ν) |/M (0ν) ) between the present results and the results of [21]. l-NMEs, present results l-NMEs [21] Nuclear transition We use the Bonn G matrix as the two-body interaction and fine tune it in the standard way (see, e.g., [31][32][33]): The pairing parameters of the BCS were adjusted by fitting the phenomenological pairing gaps, extracted from the nucleon separation energies. The residual Hamiltonian for the pnQRPA calculation contains two scaling parameters: one in the particle-hole channel (the g ph parameter) and one in the particle-particle channel (the g pp parameter). These two parameters carry over to the A and B matrices of the pnQRPA [31] such that the particle-hole contribution is proportional to the particle-hole matrix element, i.e., ∝ g ph pn −1 ; J π |V |p n −1 ; J π , where J π is the multipole of the states in the intermediate odd-odd nucleus, and the particleparticle contribution is proportional to the two-body matrix element, i.e., ∝ g pp pn; J π |V |p n ; J π .
Traditionally, the g ph parameter is fixed by fitting the centroid of the giant Gamow-Teller resonance (GTGR) in the 1 + channel of the calculations. This is also the procedure which we adopt in the present calculations. In the cases where data are available the computed beta strength, B(GT), of the GTGR is larger than the measured one. This is usually associatedwith the missing of strength in the charge-exchange reaction experiments probing the GTGR region. The missing of strength is reflected in the unsatisfactory fulfillment of the model-independent Ikeda 3(N − Z) sum rule by the data. Contrariwise, this sum rule is perfectly obeyed by the pnQRPA calculations. The 0νββ NMEs depend only weakly on the g ph parameter: an increase in g ph lifts the GTGR energy higher and reduces the Gamow-Teller strength at low energies. This, in turn, makes smaller g pp values reproduce the experimental 2νββ-decay half-life in a pnQRPA calculation.
The experimental location of the isobaric analog state (IAS) is not always well reproduced by the pnQRPA calculations. This owes to the fact that the present pnQRPA calculations are not self-consistent, i.e., the single-particle energies are not generated by the same Hamiltonian which is used in the pnQRPA calculations, as discussed in [34]. This might increase artificially the 0 + contributions in the 0νββ NMEs, but as seen in Sec. III B the 0 + contributions are quite small and thus under control.
The g pp parameter can be adjusted by fitting the measured 2νββ-decay half-life, compiled, e.g., in [35], and more recently in [36]. This procedure was followed in, e.g., [17][18][19][20][21]. Recently, an improved method was proposed in [21] where the NMEs corresponding to the exchange of light Majorana neutrinos (l-NMEs) were treated for the conservation of the isospin symmetry. There the particle-particle parts of the pnQRPA matrices, referring to the above discussion, were divided into isoscalar (T = 0) and isovector (T = 1) parts by the decomposition g pp pn; J π |V |p n ; J π → g T =1 pp pn; J π ; T = 1|V |p n ; J π ; T = 1 + g T =0 pp pn; J π ; T = 0|V |p n ; J π ; T = 0 . (46) We can now adjust the parameters g T =1 pp and g T =0 pp independently in the following way: The isovector parameter g T =1 pp can be adjusted such that the Fermi NME of (10) vanishes and thus the isospin symmetry is restored for the 2νββ decay. We then keep this adjusted value of g T =1 pp in the further calculations for the 0νββ decay. We can independently vary g T =0 pp to reproduce the measured 2νββ-decay half-life and again use this value in the calculation of the 0νββ NMEs. It should be noted here that for the A = 110 and A = 124 nuclear systems there is no 2νββ-decay data so we used beta-decay data for A = 110 and for A = 124 we took the 2νββ NMEs proposed in [21]. Like in [21], we found in the present study that the value of the isovector parameter is very close to the average of the pairing parameters of the BCS, fitted in the initial and final even-even nuclei for each 0νββ transition (see Table I of [21] for reference). We list the values of all these parameters for the present calculations in Table I.

B. Matrix elements for light-neutrino exchange
Let us start the presentation of our results by a discussion of the NMEs corresponding to the exchange of a light Majorana neutrino (the l-NMEs). The l-NMEs have been decomposed according to the terms of (23)- (25) in Table II The isoscalar parameter g T =0 pp , which was adjusted to reproduce the 2νββ decay data, forces the 1 + contribution to have a negative value in nuclei 96 Zr, 100 Mo, and 124 Sn. The interference with the rest of the contributions reduces the magnitude of the final l-NME for these decays. This effect was already noticed in [19] for the 96 Zr and 100 Mo cases. Here it has to be pointed out that although the 1 + contribution to the 0νββ l-NME is negative in the quoted three cases, the 2νββ NME is positive and it has not yet crossed zero to the negative side. This also means that the pnQRPA solutions are at a safe distance from their breaking point. Due to the neutrino potential present in (40) the sign of the total 1 + contribution can change when going from the 2νββ to the 0νββ decay. The leading multipole component is 1 − for a majority of decays exceptions being the decays of 76 Ge, 82 Se, and 136 Xe. For  show the decomposition of the presently calculated total NMEs (column 6) in terms of the Fermi, Gamow-Teller, and tensor contributions. The total NMEs are compared with those of [3] in the fourth last column, and with the results of [25] in the second last column. The third last and last columns show the ratios between our present matrix elements and those computed in Refs. [3] and [25].
h-NMEs, present results h-NMEs  [17] and [19]. We present our final l-NMEs (16) in Table III and compare them there with the results of Ref. [21] for two values of the axial-vector coupling constant, namely the effective one, g eff A = 1.00 and the bare one g b A = 1.26 (Ref. [21] uses g b A = 1.27 but the difference is insignificant). The agreement is quite good in most cases. However, there seems to be larger deviations for the nuclei 100 Mo and 124 Sn and perhaps also for 82 Se. In the case of 100 Mo the Tübingen-Caltec result is about 50% larger than ours. In the case of 124 Sn the magnitudes go in the opposite direction, the Tübingen-Caltec matrix element being only about 60% of our result. The differences between the two calculations vary quite much from one nucleus to the other. One reason may be the g pp sensitivity of the 1 + contribution: The details of the 1 + contribution may vary between the two calculations, in particular concerning the negative contribution which is very case dependent.
As can be seen from Table III the NMEs for g b A = 1.26 are much smaller than those for g eff A = 1.00. This stems from two sources: (a) The factor (g V /g A ) 2 in (16) diminishes the contribution from the Fermi matrix element (which is the same for both values of g A ), and (b) the magnitude of the Gamow-Teller NME diminishes due to the reduced 1 + contribution to it, caused by the fitting of the 2νββ half-life. The negative 1 + contribution in some of the 0νββ NMEs enhances the effect even further, like in the cases of 96 Zr, 100 Mo, and 124 Sn. The choice of the two bordering values g eff A = 1.00 and g b A = 1.26 is based on a kind of ad hoc convention: The bare value g b A = 1.26 is a natural upper limit and g eff A = 1.00 corresponds to the commonly adopted shell-model quenching [37,38]. The  TABLE VII. Effect of the isospin symmetry restoration on the l-NMEs for selected nuclei. In the calculations marked by "old" there is only one particle-particle parameter g pp = g T =0 pp , with values displayed in Table I. Here bare value of the axial-vector coupling was used.

Nuclear transition
Parameters issue of the effective value of g A is yet unsettled for the 0νββ decays. Some attempts have been made to address this issue by analyzing β − and 2νββ decays within the pnQRPA framework in [39][40][41][42][43] by using a much quenched effective value of g A . However, it is still unclear how this effective g A relates to the one used in the 0νββ calculations (see, e.g., the discussion in [41,43]).

C. Matrix elements for heavy-neutrino exchange
Here we present our results for the NMEs corresponding to the exchange of a heavy Majorana neutrino (the h-NMEs). These NMEs are decomposed according to (23)- (25) in Table IV. In Figs. 2(a)-2(h) we have plotted the QRPA multipole decomposition of the total h-NME M (0ν) for the decays of 76 Ge, 82 Se, 96 Zr, 100 Mo, 110 Pd, 116 Cd, 124 Sn, and 128 Te, and in Figs. 4(a) and 4(b) for decays of 130 Te and 136 Xe. All decompositions are again calculated with the bare value of the axial-vector coupling, g b A = 1.26. The multipole decompositions have much more irregular shape for the heavy neutrino exchange than for the light neutrino exchange. Multipole strength is distributed towards higher angular-momentum values. There seems to be no characteristic leading multipole, either. TABLE VIII. Effect of the isospin symmetry restoration on the h-NMEs for selected nuclei. In the calculations marked by "old" there is only one particle-particle parameter g pp = g T =0 pp , with values displayed in Table I We present our final h-NMEs (16) in Table V for g eff A = 1.00 and g b A = 1.26. We compare our NMEs with the results of Ref. [3] and [25] in the last four columns of the table (Ref. [3] uses g b A = 1.25 and Ref. [25] g b A = 1.269 but, as stated before, the difference is negligible). The primed matrix element M (0ν) of Ref. [3] is defined as As seen in the third last column of Table V our h-NMEs compare rather well with the Tübingen-Caltec matrix elements of Ref. [3]. There seems to be a significant difference between our present results and those computed by the Yale group [25]. The last column of Table V shows that an average constant of proportionality between our results and those of the Yale group is about 2.3.
In Table VI we investigate the evolution of the h-NMEs for nuclei 76 Ge, 82 Se, 100 Mo, and 130 Te when moving from the Miller-Spencer to the softer Argonne and CD-Bonn parametrizations of the short-range correlation function (42). We also compare this evolution with those given by Refs. [23], [3], and [25]. In the sixth and seventh column of Table VI the ratios are defined as the results of [3] divided by the results of [23]. Our M-S correlated h-NMES differ considerably from the old Tübingen calculations [23]. They also differ substantially from those computed by the Yale group [25], but the ratios listed in Table VI are  The effect of the isospin symmetry restoration procedure on the magnitudes of the light-exchange and heavy-exchange nuclear matrix elements is investigated in Tables VII and VIII. The Gamow-Teller and tensor parts are hardly affected at all in both cases. The effect on the Fermi matrix element, M (0ν) F , is ranging from small to negligible. For the h-NMEs there are no previous results to compare with but the results for the l-NMEs can be compared with those of [21]. The presently obtained small effect on M (0ν) F by the isospin restoration is at odds with the results of Ref. [21] where a significant reduction in the Fermi matrix element was observed when moving from the old to the new parametrization. There also the values of the old Fermi NMEs were considerably higher than in the present calculations. So, in fact, after the isospin restoration the Fermi l-NMEs of the present calculation and those of [21] are surprisingly close to each other.
The 0νββ transition operators O K in (38) depend on the relative distance r between the two decaying neutrons. In Figs. 5 and 6 we have plotted the radial dependence of the matrix elements M (0ν) F , M (0ν) GT , M (0ν) T , and the full NME M (0ν) for both the light and heavy neutrino exchange cases. The two figures show clearly the different nature of the two decay mechanisms. The light neutrino exchange, being a long range process, collects strength to the NME still beyond the nucleon separation of 8 fm. For the heavy neutrino exchange the neutrino potential is a contact interaction and the corresponding matrix element cuts off completely after about 4 fm. When doing computations with different parametrizations of the SRC it was noticed that the tensor part is not sensitive to the chosen parametrization. The reason for this behavior becomes clear from Figs. 5 and 6. In both decay modes the tensor part switches on only after the decaying nucleons are about 1 fm apart so the presence of the Jastrow function does not affect them much.

IV. CONCLUSIONS
In this work we have calculated the nuclear matrix elements of the neutrinoless double beta-minus decays mediated by the light (light NMEs) or heavy (heavy NMEs) Majorana neutrino. The matrix elements have been computed for ten key decays of immediate experimental interest. The calculations have been done by using realistic two-body interactions and single-particle bases, and including up-to-date nucleonnucleon short-range correlations, nucleon form factors and induced weak currents of nucleons. We also use a recently proposed method to improve on the isospin properties of the two-neutrino and neutrinoless double Fermi nuclear matrix elements. We decompose both the light and heavy NMEs in several different ways to benchmark our computer code and to enable clean comparison with the existing and future calculations done by using different model frameworks.
We found in the calculations that the light and heavy NMEs behave qualitatively and even quantitatively in a similar way as the ones calculated recently by the Tübingen-Caltec group [3,21]. The heavy NMEs differ considerably from the old Tübingen calculations [23] and the recent Yale calculations [25].
The uncertainties of the present calculations, stemming from the determination of the values of the model parameters, are well under control owing to the fact that experimental data were used at every step of the process. This ideology in determining the values of the pairing and pnQRPA parameters yields also the computed NMEs stable against variations in the size of the single-particle model space. The largest uncertainty pertains to the variation of the value of the axial-vector coupling constant g A . There is no known recipe how to determine the value of g A for the neutrinoless double beta decays of nuclei. The conservative range of variation g A = 1.00-1.26, adopted in this work, induces an interval of the l-NMEs and h-NMEs such that the endpoints of the interval can be inferred from Tables III and V, respectively. Due to the rather large width of the interval of the NME values, it becomes more and more important in the future to have a reliable estimate of the value of g A in the NME calculations.