{'fraction_non_alphanumeric': 0.03788765912439089, 'fraction_numerical': 0.016054092849318174, 'mean_word_length': 4.368890691543298, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Anatomy of nuclear matrix elements for neutrinoless double-beta decay

We show that, within the Quasiparticle Random Phase Approximation (QRPA) and the renormalized QRPA (RQRPA) based on the Bonn CD nucleon-nucleon interaction, the competition between the pairing and the neutron-proton particle-particle and particle-hole interactions causes contributions to the neutrinoless double-beta decay matrix element to nearly vanish at internucleon distances of more than 2 or 3 fermis. As a result, the matrix element is more sensitive to short-range/high-momentum physics than one naively expects. We analyze various ways of treating that physics and quantify the uncertainty it produces in the matrix elements, with three different treatments of short-range correlations.


I. INTRODUCTION
Neutrino oscillations are firmly established (see, e.g. [1,2,3,4,5]) and demonstrate that neutrinos have masses many orders of magnitude smaller than those of charged leptons. But since the masses are nonzero, neutrinoless double-beta (0νββ) decay experiments will likely tell us sooner or later whether neutrinos are Majorana or Dirac particles [6,7,8,9,10]. Moreover, the rate of the 0νββ decay, or limits on it, can tell us about the absolute neutrino-mass scale and to some extent about the neutrino mass hierarchy 1 . But to achieve these goals we need an accurate evaluation of the nuclear matrix elements that govern the decay.
In this paper, which builds on previous publications [13,14] (which we call I and II), we analyze some of the physics affecting the nuclear matrix element M 0ν -the competition between pairing and neutron-proton particle-particle correlations, the nonintuitive dependence of the decay amplitude on internucleon distance, and the treatment of short-range correlations and other high-momentum phenomena -that have not been sufficiently discussed before. As in our earlier papers (and most attempts to evaluate M 0ν ) we use the Quasiparticle Random Phase Approximation (QRPA) and its generalization, the Renormalized QRPA (RQRPA), with an interaction obtained from the G matrix associated with the realistic Bonn CD nucleon-nucleon interaction. That interaction, slightly renormalized, is used both as the like particle pairing and as the neutron-proton force. Where appropriate, we compare the results to those of the com-plementary Large-Scale Shell Model (LSSM).
The paper is organized as follows: In the next section, after briefly summarizing the relevant formalism, we show that the final value of M 0ν reflects two competing forces: the like particle pairing interaction that leads to the smearing of Fermi levels and the residual neutronproton interaction that, through ground state correlations, admixes "broken-pair" (higher-seniority) states. A partial cancellation between these interactions increases the sensitivity to their strengths. The same tendencies are present in the LSSM, as a recent paper shows [15]. (During the processing of this manuscript, a new paper on the LSSM appeared [16], emphasizing the competition again.) In section III we discuss the dependence of M 0ν on the distance between the two neutrons that are converted into two protons. We show that the competition mentioned above implies that only internucleon distances r ij < ∼ 2-3 fm contribute. That fact, not recognized before, explains the sensitivity of the decay rate to higher order terms in nucleon currents, nucleon form factors, and short-range nucleon-nucleon repulsion. We show that the surprising dependence on the internucleon distance occurs not only in the QRPA but also in an exactly solvable model [17] that contains many ingredients of real nuclear systems. Short-range correlations have recently inspired a lively discussion [19,20,21] and we devote Section IV to various ways of treating them. In Section V we present numerical results for nuclei of experimental interest, that include a comprehensive analysis, within the QRPA method and its generalization, of the total uncertainty of the 0νββ nuclear matrix elements, and compare with results of the LSSM. Section VI summarizes our findings. Finally, in Appendix A we present formulae for two ways of evaluating the matrix elements, one via the evaluation of unsymmetrized twobody matrix elements (the procedure usually used) and another one through the product of two one-body matrix elements. And in Appendix B we show how to calculate shell-model particle-hole decompositions so they can be compared with those calculated in the QRPA.

II. FORMALISM AND MULTIPOLE DECOMPOSITIONS
Throughout we assume that the 0νββ decay, if observed, is caused by the exchange of the Majorana neutrinos, the same particles observed to oscillate. The half-life of the decay is then where G 0ν (E 0 , Z) is a precisely calculable phase-space factor and M 0ν is the nuclear matrix element. The effective Majorana neutrino mass m ββ is related to the absolute mass scale and oscillation parameters through where U ei is the first row of the neutrino mixing matrix and the and α i are unknown Majorana phases. Any uncertainty in M 0ν makes the value of m ββ equally uncertain.
As stated above, we use the QRPA and RQRPA methods based on the G matrix derived from the realistic Bonn CD nucleon-nucleon force, i.e., the many body hamiltonian is We describe in detail in Section V below the input used to solve the corresponding well known equations of motion. In the QRPA (and RQRPA) M 0ν is written as a sum over the virtual intermediate states, labeled by their angular momentum and parity J π and indices k i and k f (explanations of the notation are in Appendix A, and II): and Tensor (T) contain neutrino potentials and spin and isospin operators, and RPA energies E ki,k f J π . The neutrino potentials, in turn, are integrals over the exchanged momentum q, . The results of LSSM (dark histogram) [23] and QRPA treatments (lighter histogram) are compared. Both calculations use the same single-particle spaces: (f 5/2 , p 3/2 , p 1/2 , g 9/2 ) for 82 Se and (g 7/2 , d 5/2 , d 3/2 , s 1/2 , h 11/2 ) for 130 Te. In the QRPA calculation the particle-particle interaction was adjusted to reproduce the experimental 2νββ-decay rate.
The functions f F,GT (qr 12 ) = j 0 (qr 12 ) and f T (qr 12 ) = j 2 (qr 12 ) are spherical Bessel functions (the sign of j 2 was given incorrectly in Ref. [14]). The functions h K (q 2 ) are defined in Appendix A and in II. The potentials depend explicitly, though rather weakly, on the energies of the virtual intermediate states, E k J π . The functionf (r 12 ) in Eq. 4 represents the effects of short range correlations. These will be discussed in detail in Section IV.
Two separate multipole decompositions are built into Eq. (4). One, already mentioned, is in terms the J π of the virtual states in the intermediate nucleus, the good quantum numbers of the QRPA and RQRPA. The other decomposition is based on the angular momenta and parities J π of the pairs of neutrons that are transformed into protons with the same J π (we drop the superscript π from now on for convenience). This latter representation is particularly revealing. In Fig. 1 we illustrate it both in the LSSM and QRPA, with the same single-single particle spaces in each. These two rather different approaches agree in a semiquantitative way, but the LSSM entries for J > 0 are systematically smaller in absolute value.
Ref. [15] makes the claim that QRPA results are too large because they omit configurations with seniority greater than 4, which are especially effective in canceling the pairing part of the matrix element. This statement is not correct. The QRPA does include configurations with higher seniority (4,8,12, etc.) and, as the Fig. 1 shows, the broken pair contributions to the matrix elements are as large or larger than in the LSSM. (Some of the difference might be due to differences in single-particle energies and occupation numbers, which are not identical in the two calculations even though the single-particle wave functions are.) The reason that the QRPA results presented in Section V are larger are a somewhat greater pairing contribution and contributions from negative parity multipoles that reinforce it. Most of the negative-parity contributions are absent from the shell model because of restrictions on the model space. These results suggest that the shell model is as likely to be missing important physics as is the QRPA. We return to this point in Section V.
In Fig. 2 we show the J decomposition for three nuclei in the QRPA, with single-particle spaces encompassing two major shells, a more natural span for this method. The cancellation between components with J = 0 pairs and with J = 0 pairs is always pronounced. The net M 0ν is considerably smaller than the pairing contribution and so depends rather sensitively on the pairing physics that determines the J = 0 part, as well as on the strength of the proton-neutron force that determines the J = 0 part.
From the structure of Eq. (4) and in particular from the form of the reduced matrix elements 0 it is obvious that only one of the two possible couplings between the neutron and proton operators in the two-body matrix element p(1), p ′ (2); J f (r 12 )O Kf (r 12 ) n(1), n ′ (2); J is realized. This means that this two-body matrix element should not be antisymmetrized. In the LSSM one typically uses the closure approximation, which represents M 0ν as the ground-state-to-ground-state transition matrix element of a two-body operator. M 0ν can then be rewritten purely in terms of the antisymmetrized twobody matrix elements. After antisymmetrization, however, it is not possible to recover the decomposition into the multipoles J π of the virtual intermediate states. In Appendix B we show how shell-model practitioners, by retaining unsymmetrized matrix elements, can decompose the matrix element into intermediate-state multipoles J π for comparison with QRPA calculations. We cannot, however, make the comparison here without more shell-model data than has been published.
When using the QRPA or RQRPA to evaluate M 0ν , one must fix several important parameters, the effects of which were discussed in detail in I and II. The strength g pp by which we renormalize the Bonn-CD G matrix in the neutron-proton particle-particle channel is particularly important. We argued in I and II that g pp should be chosen to reproduce the rate of two-neutrino ββ decay. This choice, among other things, essentially removes the dependence of M 0ν on the number of the single-particle states (or oscillator shells) in the calculations. The 2ν matrix element depends only on the 1 + multipole. In Fig. 3 we show that it is essentially this multipole that is responsible for the rapid variation of M 0ν with g pp . Fixing its contribution to a related observable (2ν decay) involving the same initial and final nuclear states appears to be an optimal procedure for determining g pp .

III. DEPENDENCE ON THE DISTANCE BETWEEN THE NUCLEONS INVOLVED IN
THE 0νββ TRANSITION.
The operators O K in Eq. (4) depend on the distance r 12 between the two neutrons that are transformed into protons. The corresponding neutrino potentials are the Fourier transforms over the neutrino momentum q as shown in Eq. (5). Obviously, the range of r 12 is restricted from above by r 12 ≤ 2R nucl . We show here, however, that in reality only much smaller values, r 12 < ∼ 2-3 fm, or equivalently larger values of q, are relevant. Thus a good description of the physics involving distances r 12 ∼ 1 fm, or q ∼ 200 MeV is important. That finding has not been recognized before, but perhaps it should be not so surprising that q ∼ p F ermi is the most relevant momentum transfer.
An example of the r 12 dependence of M 0ν is shown in Fig. 4 for three nuclei. The quantity C(r) is defined by evaluating M 0ν after multiplying H K (r ′ , E k J π ) by r 2 δ(r− r ′ ), so that C(r) is the contribution at r to M 0ν , with ∞ 0 C(r)dr = M 0ν . As the lower panel of the figure demonstrates, the cancellation between the J = 0 and J = 0 components is essentially complete for r 12 > ∼ 2-3 fm. Since the typical distance from a particular nucleon to its nearest neighbor is ∼ 1.7 fm (because R nucl = 1.2A 1/3 ) the nucleons participating in the 0νββ decay are mostly nearest neighbors. Short-range nucleon-nucleon repulsion, the finite nucleon size, represented by nucleon form factors, and components of the weak currents that are typically suppressed by q/M nucleon are therefore more important than one would naively expect.
Perhaps the most interesting thing about the figure is that the pairing and non-pairing parts of C(r) taken individually (as in the two panels of the figure) extend to significantly larger r. The cancellation between them, that we discussed earlier, is particularly effective beyond 2 or 3 fm, leaving essentially nothing there. Figure 5 shows that the shape of C(r), like the integrated matrix element, is essentially independent of the number of single-particle orbits included, as long as the truncation is not too severe (as it is with the dash-double-dot curve, for which important spin-orbit partners were omittedonly the 4 single particle states p 3/2 , p 1/2 , f 5/2 , g 9/2 were included) and the coupling constant g pp is chosen to re- produce the measured 2νββ lifetime. For other values of g pp the cancellation between the J = 0 and J = 0 contributions at r larger than 2 or 3 fm is not as complete as in Fig. 4. We return to this point shortly We show the r 12 dependence of the different parts of the M 0ν in Fig. 6. All individual contributions die out at r larger than 2 or 3 fm. The pseudoscalar-axial vector interference part has opposite sign from the other contributions, and essentially (and accidentally) cancels the contributions of the vector, weak magnetism and pure pseudoscalar pieces. The higher-order terms reduce the matrix element noticeably, and have to be included.
To gain some insight into the renormalization of the double-beta decay operator in the shell model, Ref. [17] employs a solvable model based on the algebra SO(5) × SO(5). The valence space contains two major shells (f pg 9/2 and sdg 7/2 ), split by an energy ǫ, with degenerate levels within each shell. The single-particle wave functions are taken from a harmonic oscillator withhω =9.2 MeV. The Hamiltonian (for this schematic model only, not in the rest of the paper) is where a, b = 1, 2 label the shells (lower and upper), ǫ is the energy difference between the shells,N 2 is the number operator for the upper shell, T a is total isospin operator for shell a, and Here π † α (ν † α ) creates a proton (neutron) in level α with angular momentum j α ,ĵ ≡ √ 2j + 1, and the square brackets indicate angular-momentum coupling. H contains only generators of SO(5)×SO(5), so its lowest lying eigenstates consist of configurations in which the nucleons are entirely bound in isovector S pairs of the type in Eq. (7). This model has no active spin, so it is only suitable for calculating Fermi (neutrinoless) double-beta decay. To simulate the effect of g pp on Gamow-Teller decay we change the Fermi matrix element by varying the strength of isovector neutron-proton pairing, in the same way that we change the Gamow-Teller matrix element in the realistic QRPA by varying the strength of isoscalar pairing. The advantage of this model is that we can solve it exactly rather than in the QRPA.
We can use the model to test the r 12 -dependence of the double-beta decay matrix element in an exact solution. (Analytic expressions for the necessary matrix elements are in Ref. [17].) The upper panel of Fig. 7 shows the dependence for several values of g pp , with g ph = 0 and ǫ = 10G. Just as in the realistic QRPA calculations, the contribution beyond r = 3 fm is very small for g pp around 1; it is too small to distinguish from zero in the figure beyond 5 fm for g pp = 0.8. However, for other values, as noted above, the large-r contributions can be substantial. The bottom panel divides the function into like-particle pairing and non-pairing parts for g pp = 0.8. The two cancel to high precision at large r. The suppression at long ranges we observe in the QRPA, then, appears to be fairly general. It happens even in a very simple model, solved exactly.(As noted above, a new preprint [18] appeared during the processing of this manuscript. In it the r 12 dependence of the M 0ν , as well as the dependence of the separated pairing and broken-pair contributions, was evaluated in the LSSM. That analysis, inspired by our work, yielded curves that are strikingly similar to those in Fig. 4.) The resulting M 0ν values are 5.32 when the effect of short range correlations is ignored, 5.01 when the UCOM transformation [24] is applied, 4.14 when thef (r12) from Fermi hypernetted-chain calculations [25] is used in Eq.(4), and 3.98 when the phenomenological Jastrowf (r12) is used [26].

IV. SHORT-RANGE CORRELATIONS AND OTHER HIGH-MOMENTUM PHENOMENA
Since only r 12 < ∼ 2 − 3 fm, (equivalently q >hc/(2 − 3 fm)), contributes to M 0ν , some otherwise negligible effects become important. These effects are not commonly included, or included only in rough approximation, in nuclear-structure calculations. For example, the dipole approximation for nucleon form factors and the corresponding parameters M V and M A come from electronand neutrino charged-current-scattering from on-shell nucleons. Nuclear structure deals with bound nucleons and virtual neutrinos that are far off-shell. Similarly, the induced pseudoscalar current, with its strength obtained from the Goldberger-Treiman relation, has been tested in muon capture on simple systems. Here we are using this current for off-shell virtual neutrinos. Short range nucleon-nucleon repulsion has been considered carefully when calculating nuclear binding energies, but here we need its effect on a transition operator connecting two different nuclear ground states. All these effects will introduce some uncertainty because their treatment is not well tested. Nevertheless, it is important to understand their size at least roughly.
To show the importance of high momenta explicitly, we display in Fig. 8 the q dependence C(q) -defined in complete analogy to C(r) -of M 0ν in 76 Ge, in a similar manner to which we exhibited the r 12 dependence earlier.
The cancellation between the J = 0 and J = 0 parts is particularly complete at lower values of q so that the resulting curve in the upper panel, although reduced in magnitude, is clearly shifted towards higher q.
The first high-momentum effect we examine is shortrange correlations. Fig. 9 displays the r 12 dependence of M 0ν for several methods of handling short-range physics.
For obvious reasons all methods reduce the magnitude of M 0ν . The Unitary Correlation Operator Method (UCOM) [24] leads to the smallest reduction, less than 5%. The phenomenological Jastrow-like functionf (r 12 ) in Eq.(4) (from Ref. [26]) reduces M 0ν by about 20%. We also display the results of using a microscopicallyderived Jastrow function [25]; its effect is similar to that of the phenomenological function. Since it is not clear which approach is best, we believe it prudent to treat the differences as a relatively modest uncertainty.
Nucleon form factors pose fewer problems because it turns out that once the short-range correlations effects are included, no matter how, the form factors are almost irrelevant as long as the cut-off masses M A,V are at least as large as the standard values (M A = 1.09 GeV and M V = 0.85 GeV). In Fig. 10 we show the dependence of M 0ν on the values of M A,V which for this purpose are set equal to each other, with three alternatives for treating short-range correlations. By 2 GeV the curves have essentially reached the infinite-mass limit. Since they are essentially flat past 1 GeV for both the UCOM and Jastrow-like prescriptions, including the form factors causes only minor changes in M 0ν . Only if the correlations are ignored altogether do the form factors make a significant difference.
Finally, there is little doubt that the higher order weak currents, induced pseudoscalar and weak magnetism, should be included in the calculation. Even though the Goldberger-Treiman relation has not been tested in twobody operators, the relation is sufficiently well established that we do not associate a sizable uncertainty with its use. V. NUMERICAL RESULTS Fig. 11 shows our calculated ranges for M ′0ν , defined as M ′0ν = (g A /1.25) 2 M 0ν to allow us to display the effects of uncertainties in g A . Such a definition allows us to use the same phase space factor G 0ν with g A = 1.25 when calculating the 0νββ-decay rate.
We consider the effects of short-range correlations an uncertainty, since we don't know the best way of treating them. The error bars in Fig. 11 represent the difference between the highest and lowest of 24 calculations -in either the QRPA or RQRPA, each with 3 different sets of single-particle states (usually 2, 3 and 4 oscillator shells), 2 values for the in-medium g A (1.0 and 1.25), and 2 treatments of short-range correlations (phenomeological Jastrow functions and the UCOM method) -and include the experimental uncertainty in the values of the 2ν lifetimes used to determine g pp . Thus the error bars displayed in Fig. 11 represent our estimate of the full uncertainty in the 0νββ matrix elements within the QRPA and RQRPA methods. Though the results are in reasonable agreement with those presented in Refs. [20,21], the uncertainties are different.
Reference [27] presents new LSSM results for M ′0ν , the values of which are shown also in Fig. 11. They are somewhat smaller than the QRPA values for 76 Ge and 82 Se and in a fair agreement for 116 Cd, 128 Te, 130 Te and 136 Xe. Various nuclear-structure effects may be responsible for the discrepancies. More complicated configurations in LSSM that are absent in QRPA typically reduce M 0ν and could make most of the difference. On the other hand, in the QRPA includes more single particle states than in LSSM. That has the tendency to increase M 0ν and could also be responsible for the discrepancy. As we said earlier, the notion that the QRPA omits highseniority states is not correct and shouldn't be used to argue that the LSSM results are more accurate.
Given the interest in the subject, we show the range of predicted half-lives corresponding to our full range of M ′0ν in Table I (for m ββ = 50 meV). As we argued above, this is a rather conservative range within the QRPA and its related frameworks. One should keep in mind, however, the discrepancy between the QRPA and LSSM results as well as systematic effects that might elude either or both calculations.

VI. CONCLUSIONS
The most important and novel result here is that the generic competition between J = 0 (pairing) and J = 0 (broken pair) multipoles leads to almost complete cancellation of the contribution to the matrix element from internucleon distances r > ∼ 2-3 fm. That explains why the effects that depend on smaller values of r, or equivalently larger momentum transfers q, become important. This competition also means that the final matrix elements have enhanced sensitivity to the strengths of these interactions. Despite the uncertainties associated with the short range effects, we conclude that a proper fitting of the QRPA and/or RQRPA parameters leads to a relatively narrow range for M ′0ν , with a smooth dependence ( 96 Zr being an exception) on atomic charge and mass.
We evaluate the values of the matrix elements for nuclei of experimental interest and display our best estimate of the corresponding spread. Part of that spread is associated with the difference in the size of the single particle space and whether QRPA or RQRPA is used, as discussed earlier in I and II. An interesting new conclusion is that short-range correlations, no matter how they are treated, essentially eliminate the effect of finite nucleon size on the matrix elements. But we still do not know the best way to treat the correlations, a fact that contributes about 20% to uncertainties presented above. The uncertainty in the effective value of g A contributes about 30%, with the rest due to choice of method and model space, and the experimental uncertainty in 2ν lifetimes. Here we outline the derivation of M 0ν , emphasizing induced (higher-order) currents. We assume light-neutrino exchange throughout, and the standard lepton×hadron weak charged-current Hamiltonian. The hadronic current, expressed in terms of nucleon fields Ψ, is where m p is the nucleon mass and q µ is the momentum transfer, i.e. the momentum of the virtual neutrino. Since in the 0νββ decay q 2 ≫ q 2 0 we take q 2 ≃ − q 2 . For the vector and axial vector form factors we adopt the usual dipole approximation g V ( MeV, and M A = 1086 MeV. We use the usual form for weak magnetism, g M ( q 2 ) = (µ p −µ n )g V ( q 2 ), and the Goldberger-Treiman relation, g P ( q 2 ) = 2m p g A ( q 2 )/( q 2 + m 2 π ), for the induced pseudoscalar term.
To derive the expression for the matrix element we follow the procedure outlined in Ref. [28], arriving after a few steps at an expression for the 0 + i → 0 + f ground state to ground state transition: We have made the (accurate) approximation that all electrons are emitted in the s 1/2 state, with energies equal to half the available energy (E i − E f )/2. The normalization factor 4π R/g 2 A , introduced for convenience, is compensated for by corresponding factors in the phase space integral.
The reduced matrix elements of the one-body operators c + pc n (the tilde denotes a time-reversed state) in Eq. (A8) depend on the BCS coefficients u i , v j and on the QRPA vectors X, Y [22]. The nuclear structure information resides in these quantities.
The final step, leading to Eq. (4) in the text, is to rewrite the product of two one-body matrix elements as an appropriately recoupled (with pairs of protons and neutrons coupled to angular momentum J ) unsymmetrized two-body matrix element. Without the complications of angular momentum, this step simply reads p|O(1)|n p ′ |O(2)|n ′ = p, p ′ |O ′ (1, 2)|n, n ′ . (A11) We then transform to relative and center-of-mass coordinates r 12 = r 1 − r 2 and R 12 = ( r 1 + r 2 )/2. Since the exchange potential depends only on r 12 = | r 12 | we end up with Eq. (4). The Fermi (F), Gamow-Teller (GT) and Tensor (T) operators in that equation are