Impact of the leading-order short-range nuclear matrix element on the neutrinoless double-beta decay of medium-mass and heavy nuclei

We evaluate the leading-order short-range nuclear matrix element for the neutrinoless double-beta ($0\nu\beta\beta$) decay of the nuclei most relevant for experiments, including $^{76}$Ge, $^{100}$Mo, $^{130}$Te and $^{136}$Xe. In our calculations, performed with the nuclear shell model and proton-neutron quasiparticle random-phase approximation (pnQRPA) methods, we estimate the coupling of this term by the contact charge-independence-breaking coupling of various nuclear Hamiltonians. Our results suggest a significant impact of the short-range matrix element, which is about $15\%-50\%$ and $30\%-80\%$ of the standard $0\nu\beta\beta$-decay long-range matrix element for the shell model and pnQRPA, respectively. Combining the full matrix elements with the results from current $0\nu\beta\beta$-decay experiments we find that, if both matrix elements carry the same sign, these searches move notably toward probing the inverted mass ordering of neutrino masses.


Introduction
Neutrinoless double-beta (0νββ) decay is a hypothetical nuclear process where two neutrons turn into two protons and only two electrons are emitted. Since two leptons are created, an observation of this decay would point to an event beyond the standard model (SM) of particle physics. Further, 0νββ decay can only occur if the neutrino is a Majorana particle-its own antiparticle-unlike any other fundamental particle known. The observation of 0νββ decay would therefore shed light on beyond-SM physics such as the matter-antimatter asymmetry of the universe [1,2], and the nature of the neutrino [3][4][5][6][7].
Furthermore, recently Refs. [47,48] introduced a previously unacknowledged short-range matrix element which appears at leading order in lightneutrino-exchange 0νββ decay. This brings an additional, potentially significant uncertainty to 0νββ-decay NMEs, especially because the value of the hadronic coupling associated with the new term is not known. First quantum Monte Carlo studies in very light 12 Be, estimating the new coupling from the charge-independence-breaking (CIB) coupling of nuclear Hamiltonians, indicate that this term could amount to as much as 80% of the standard long-range NME [48]. The new term was also studied in 48 Ca with coupled cluster theory [31]. More recently, Refs. [49,50] provide synthetic data to fix the coupling of the short-range term in ab initio calculations, a procedure leading to a 43(7)% enhancement of the 48 Ca IMSRG NME [51].
In this Letter, we extend these studies and explore for the first time the short-range 0νββ-decay NME in a wide range of ββ emitters, including all nuclei used in the most advanced experiments [9][10][11][12][13][14] and proposals for next generation searches [52][53][54]. We perform many-body calculations of medium-mass and heavy nuclei with nucleon number A = 48 − 136 with the pnQRPA and large-scale NSM frameworks commonly used to obtain long-range NMEs, and we take CIB couplings to estimate the size of the short-range NMEs. Finally, we analyze the impact of this new term on current 0νββ-decay searches combining in a consistent manner the likelihood functions of the most constraining experiments with the full NMEs for different nuclei.

Neutrinoless Double-Beta Decay
The 0νββ-decay half-life can be written as [5,47] [t 0ν where G 0ν is a phase-space factor for the finalstate leptons [55], and M 0ν L and M 0ν S are the longand short-range NMEs, with unknown relative sign. The effective Majorana mass m ββ = i U ei m i (normalized to the electron mass m e ) characterizes the lepton-number violation and depends on the neutrino masses m i and mixing matrix U .
The matrix element M 0ν L denotes the standard light-neutrino-exchange matrix element, which can be written in the familiar way [5] M 0ν with Gamow-Teller, Fermi and tensor contributions M 0ν GT , M 0ν F and M 0ν T , and vector coupling g V = 1.0. The calculation of the matrix elements involves, in addition to the initial and final states 0 + i and 0 + f , a sum over intermediate states, carried out explicitly in the pnQRPA [20]. Alternatively, for our NSM results we use the closure approximation to perform this sum analytically, so that the dominant Gamow-Teller term reads with sum over the spin σ and isospin τ − operators of all A nucleons, momentum transfer q, a Bessel function j 0 , and R = 1.2A 1/3 fm. We use as average energy for the intermediate states E = 0. The matrix element also depends on a neutrino potential, with h GT (0) = g 2 A and additional q-dependent subleading terms, regularized with a dipole as in previous pnQRPA [20,56] and NSM [15] studies. The Fermi and tensor parts follow similar expressions to Eq. (3) [5]. Finally, we correct our manybody states with two-nucleon short-range correlations (SRCs) following the so-called CD-Bonn and Argonne parametrizations [57].
The short-range matrix element connects directly the initial and final nuclei [47] where we choose to regularize the contact term with a Gaussian in the neutrino potential: with Λ the scale of the regulator. The coupling g NN ν , not part of the SRCs, can only be fixed by fitting to lepton-number-violating data-currently unavailable-or synthetic data [49,50]-only accessible to ab initio calculations. Here we follow Ref. [48] and estimate its value by considering the CIB term of different nuclear Hamiltonians, restricted to cases with a Gaussian regulator. This strategy carries some uncertainty since two low-energy constants are needed to fix g NN ν , while only one of them can be extracted from the CIB term-the other one is assumed to have the same value. Nonetheless, the empirical CIB was well reproduced in Refs. [49,50] with the same strategy used to obtain their synthetic lepton-numberviolating data. Table 1 shows the g NN ν − Λ pairs considered in our work.

Many-Body NME Calculations
We perform NSM calculations with the coupled code NATHAN [60]. We use the KB3G [61] interaction in the pf -shell (0f 7/2 , 1p 3/2 , 0f 5/2 and 1p 1/2 orbitals) for A = 48, the RG.5-45 [62] interaction with 1p 3/2 , 0f 5/2 , 1p 1/2 and 0g 9/2 orbitals for A = 76, 82 and the GCN5082 [62] interaction with 0g 7/2 , 1d 5/2 , 1d 5/2 , 2s 1/2 and 0h 11/2 orbitals for A = 124 − 136. In all cases our valence space is common to protons and neutrons. Overall, with the NSM we study seven 0νββ decays: 48 Ca, 76 Ge, 82 Se, 124 Sn, 128,130 Te and 136 Xe. We consider all nuclear configurations in the full valence space except in 124 Te (final nucleus of the 124 Sn ββ decay) which is limited to seniority v ≤ 5 states (up to five broken zero-angular-momentum pairs) instead On the other hand, we use the spherical pn-QRPA method as in Refs. [21,63]. The large nocore single-particle bases consist of 18 48 Ca because the pnQRPA does not describe doubly-magic nuclei reliably. We take the singleparticle energies from a Coulomb-corrected Woods-Saxon potential optimized for nuclei close to the β-stability line [64], but in the vicinity of the Fermi surface we slightly modify them to better reproduce the low-lying spectra of neighboring odd-mass nuclei. The quasiparticle spectra, needed in the pn-QRPA diagonalization, follow the solution of the BCS equations for protons and neutrons. We use the two-body interaction derived from the Bonn-A potential [65], fine-tuning the proton and neutron pairing parameters to reproduce the phenomenological pairing gaps. The residual Hamiltonian for the pnQRPA calculation contains two adjustable scaling factors: the particle-hole g ph and particleparticle g pp parameters. We fix g ph to reproduce the centroid of the Gamow-Teller giant resonance, and g pp to the two-neutrino ββ-decay half-life according to the partial isospin-symmetry restoration scheme introduced in Ref. [66].

Results and Discussion
We calculate the 0νββ-decay short-and longrange NMEs for ten heavy nuclei, listed in Table 2. For both NSM and pnQRPA, in all transitions the standard matrix element M 0ν L is larger than the new term M 0ν S . Nonetheless, Table 2 shows that in both many-body frameworks the contribution of the short-range matrix element is significant: in the pnQRPA the ratios of the short-over longrange NMEs typically range between 30% − 80%; in the NSM the ratios are slightly more moderate, between 15% − 50%. Within a given method, the relative size of M 0ν S is in general rather stable. Our results indicate that the short-range contribution can considerably impact the expected rates of current and future 0νββ-decay experiments. There- values. The NME ranges in Table 2 are much wider in the case of the short-range term than for the standard matrix element, as the difference between the lower and upper M 0ν S values can be up to a factor of three for both methods. This partially reflects the variety of couplings g NN ν and regulator scales Λ in Table 1: the smallest short-range values are always given by g NN In order to study the short-and long-range NMEs in more detail, Fig. 1 shows their radial and momentum distributions, denoted by C L/S (r) and C L/S (q), for 76 Ge. The distributions satisfy where r = |r n − r m | is the distance between the two decaying nucleons. Figure 1 shows the short-range distributions calculated with the two extreme neutrino potentials: g NN  both NSM and pnQRPA. As expected, the radial distribution of M 0ν S involves shorter internucleon distances than the one of the long-range M 0ν L , and its momentum distribution reaches larger momentum transfers. Apart from the consistently smaller NME values obtained with the NSM, the overall behavior of the matrix-element distributions is quite similar in both frameworks. Figure 1 also shows differences between manybody methods. In the pnQRPA, the radial distribution of the long-range NME gets a sizeable cancellation from distances r 2.5 fm, which is much milder in the NSM. This rather well-known feature of the pnQRPA [20,56,57] partly explains the relatively larger size of M 0ν S with respect to M 0ν L , since there are no cancellations in the short-range NME radial distribution. Furthermore, Fig. 1 also highlights that the short-range pnQRPA NME extends to longer distances than the NSM one, whereas the positive contribution to the pnQRPA long-range NME is concentrated at shorter distances. This behavior also leads to larger pnQRPA M 0ν S /M 0ν L ratios. In momentum space, the pnQRPA long-range NME distribution reaches larger momentum transfers, while the NSM one does not vanish at q = 0 because of our closure energy E = 0.
Two transitions stand out with the largest relative short-range M 0ν S values: 100 Mo for the pn-QRPA, and 48 Ca for the NSM. Figure 2 shows the radial short-and long-range NME distributions for these two cases. Apart from the different scales, the 48 Ca radial distribution resembles the 76 Ge pn-QRPA long-range one in Fig. 1: there is a sizeable cancellation in C L (r) at distances r ≈ (2 − 5) fm, not observed in any other NSM decay. Such cancellation never occurs for the short-range C S (r), which 0.  [68,69], compared to the exclusion (blue) bands which combine data [70] from 0νββ-decay experiments [9][10][11]14] and pnQRPA or NSM NMEs. The upper, middle and lower exclusion bands correspond to NME ranges for M 0ν L −M 0ν S , M 0ν L and M 0ν L +M 0ν S , accordingly. The cyan bands correspond to a "quenching" scenario, see text for details. explains the larger relative contribution of M 0ν S for this nucleus. The relative size of our 48 Ca shortrange NME is similar to the ab initio result from Ref. [51], however obtained with a different coupling and regulator scheme. Figure 2 also shows a more marked cancellation in the pnQRPA longrange 100 Mo NME than in 76 Ge. This exceptionally large cancellation, not present in any other nucleus, is explained by a notable negative contribution at low momenta which reduces the value of M 0ν L . This behavior is driven by the 1 + multipole which dominates at low-q values, as observed in previous pn-QRPA works [20,56]. A similar feature appears in light nuclei studied with quantum Monte Carlo [67] and the NSM.
The M 0ν L matrix elements in Table 2 assume g A = 1.27. Related NSM and pnQRPA β and two-neutrino ββ decay rates obtained this way are known to be overestimated, calling for corrections usually known as "g A quenching". While the implications to 0νββ-decay NMEs are not clear [5], they would only affect the long-range NME, leading to a larger relative impact of the short-range term. We consider such "quenching" scenario to provide more conservative estimates. On the one hand, we calculate pnQRPA (g eff A /g A ) 2 M 0ν L values with g pp fitted to reproduce two-neutrino ββ-decay half-lifes with g eff A = 1. The results are reduced by about 20%. In a similar spirit, we multiply the NSM M 0ν GT terms by 1/g 2 A 0.6, which reduces M 0ν L by 30% or so. For 48 Ca, this brings the NSM in good agreement with ab initio theory [30,32,51]. With these rough estimates, the short-range NME contribution increases by about 25% in the pnQRPA, reaching about 40%−100% in 76 Ge, 130 Te and 136 Xe. For the NSM, the impact of the short-range term would be enhanced by about 50%, typically up to 25%−70%.
Finally, we explore the impact of our M 0ν S results on the current reach of the experimental 0νββdecay program in terms of bounds on m ββ . In order to obtain stronger limits, we follow Ref. [70] to combine our 76 Ge, 130 Te and 136 Xe NMEs with the parameterized likelihood functions of the 0νββ-decay rate, extracted from the CUORE [9], GERDA (Phase II) [10], EXO-200 [11] and KamLAND-Zen [14] experiments. To combine the likelihood functions, we convert decay rates into effective Majorana masses according to Eq. (1) with our NMEs and the phase-space factors of Ref. [55]. This way, we obtain 90% confidence level (CL) upper bounds on m ββ from the 90% CL upper bound of the combined likelihood function (related to the Bayesian rather than the frequentist limit [70]).
We consider three different scenarios to derive bounds on m ββ : a baseline using the standard M 0ν L ; an optimistic scenario assuming common signs for the short-and long-range NMEs (with M 0ν L + M 0ν S in Table 2); and a pessimistic one where the shortrange part cancels the standard matrix element (with M 0ν L − M 0ν S ). The first consistent determination of the short-range matrix element in 48 Ca supports the optimistic scenario [51]. We take a matrix-element uncertainty given by the extreme values of M 0ν L and M 0ν L ±M 0ν S obtained from the set of 12 calculations corresponding to the six g NN ν − Λ pairs in Table 1 and two SRCs. For 76 Ge, 130 Te and 136 Xe, the extreme values are always given by the same g NN ν − SRC combinations. Figure 3 compares the constraints on m ββ in the three scenarios with the bands corresponding to the normal and inverted neutrino-mass orderings, obtained with neutrino-oscillation data [68] as described in Ref. [69]. The widths of the blue bands in Fig. 3 correspond to the ranges of M 0ν L (middle), M 0ν L − M 0ν S (top) and M 0ν L + M 0ν S (bottom) in Table 2, and are much larger once the new short-range term is included. The pnQRPA bands are dominated by the likelihood function of GERDA partly due to the large 76 Ge NME, while the next-constraining experiments are KamLAND-Zen, EXO-200 and CUORE, in that order. In the NSM the hierarchy is similar. In the scenario that both matrix elements carry the same sign, our results indicate that the reach of current experiments approaches the inverted mass-ordering region notably. This feature is more marked when using the pnQRPA NMEs, and may be more moderate if our results obtained with g A = 1.27 somewhat underestimate the decay half-lives. A more conservative "quenching" scenario is shown by the cyan bands in Fig. 3. On the other hand, if the signs of the two matrix elements are opposite, experiments would still be far from exploring m ββ values corresponding to the inverted mass ordering. In this case, the pnQRPA and NSM NMEs would be similar within uncertainties.

Summary
We have calculated for the first time the shortrange NME which contributes at leading order to the 0νββ decay of medium-mass and heavy nuclei including 76 Ge, 100 Mo, 130 Te and 136 Xe. Since the value of the coupling g NN ν of this short-range term is not known, we estimate it by a set of CIB couplings of different Hamiltonians, together with the corresponding regulators. We find that the new short-range NME values are a significant fraction of the standard long-range ones: typically between 30% − 80% in the pnQRPA and 15% − 50% in the NSM. These ranges are driven by the different couplings g NN ν considered, and are rather stable among all nuclei. The only exceptions are 100 Mo and 48 Ca where the ratios are notably larger due to cancellations in the standard long-range NME. Since these cancellations are typically larger in the pnQRPA than in the NSM, for the former the relative impact of the short-range matrix element is also larger.
The new short-range term can also affect the interpretation of present and future 0νββ-decay searches. To this end, we derive constraints on m ββ using our pnQRPA and NSM NMEs to combine the likelihood functions of the most constraining experiments. We observe that if the long-and shortrange NMEs carry the same sign, as suggested by a recent determination [51], the m ββ values constrained by these searches clearly approach the inverted neutrino-mass region.