A new leading contribution to neutrinoless double-beta decay

Within the framework of chiral effective field theory we discuss the leading contributions to the neutrinoless double-beta decay transition operator induced by light Majorana neutrinos. Based on renormalization arguments in both dimensional regularization with minimal subtraction and a coordinate-space cutoff scheme, we show the need to introduce a leading-order short-range operator, missing in all current calculations. We discuss strategies to determine the finite part of the short-range coupling by matching to lattice QCD or by relating it via chiral symmetry to isospin-breaking observables in the two-nucleon sector. Finally, we speculate on the impact of this new contribution on nuclear matrix elements of relevance to experiment.

Within the framework of chiral effective field theory we discuss the leading contributions to the neutrinoless double-beta decay transition operator induced by light Majorana neutrinos. Based on renormalization arguments in both dimensional regularization with minimal subtraction and a coordinate-space cutoff scheme, we show the need to introduce a leading-order short-range operator, missing in all current calculations. We discuss strategies to determine the finite part of the shortrange coupling by matching to lattice QCD or by relating it via chiral symmetry to isospin-breaking observables in the two-nucleon sector. Finally, we speculate on the impact of this new contribution on nuclear matrix elements of relevance to experiment.

Introduction:
Neutrinoless double-beta decay (0νββ) is the most sensitive laboratory probe of lepton number violation (LNV). In 0νββ L is violated by two units when two neutrons in a nucleus turn into two protons, with the emission of two electrons and no neutrinos. The observation of 0νββ would demonstrate that neutrinos are Majorana fermions [1], shed light on the mechanism of neutrino mass generation [2][3][4], and give insight into leptogenesis scenarios for the generation of the matter-antimatter asymmetry in the universe [5].
0νββ is actively being searched for in a number of even-even nuclei for which single-β decay is energetically forbidden. Current experimental limits [6][7][8][9][10][11][12][13][14][15] on the half-lives are at the level of T 1/2 > 5.3 × 10 25 y for 76 Ge [12] and T 1/2 > 1.07 × 10 26 y for 136 Xe [10], with next-generation ton-scale experiments aiming at improvements in sensitivity by two orders of magnitude. 0νββ can be generated by a variety of dynamical LNV mechanisms, which in an effective field theory (EFT) approach to new physics are parametrized by ∆L = 2 operators of odd dimension greater than four [16][17][18][19][20][21][22]. If the mass scale associated with LNV is much higher than the electroweak scale, the only low-energy manifestation of this new physics is a Majorana mass for light neutrinos, encoded in a single gauge-invariant dimension-five operator [16], which induces 0νββ through light Majorananeutrino exchange [23,24]. To interpret positive or null 0νββ results in this minimal scenario it is crucial to have good control over the relevant hadronic and nuclear matrix elements. Current knowledge of these is not satisfactory [25], as various many-body approaches lead to estimates that differ by a factor of two to three and most calculations are not based on a modern EFT analysis. In Ref. [26] a first step was presented towards the analysis of 0νββ induced by a light Majorana neutrino in the chiral EFT framework [27][28][29], which provides a systematic expansion of hadronic amplitudes . The 0νββ transition operators were derived up to next-to-next-to-leading order (N 2 LO) in Weinberg's power-counting scheme [30,31].
In this letter we demonstrate that Weinberg's scheme for 0νββ assumed in Ref. [26] breaks down and any consistent power counting requires a leading-order (LO) short-range ∆L = 2 operator, whose effect is missing in all current calculations. Our argument is based on renormalization. Using two different schemes (dimensional regularization with minimal subtraction and a coordinate-space cutoff) we show that once the strong nucleon-nucleon scattering amplitude is made finite and independent of the ultraviolet regulator, an additional ∆L = 2 contact operator with coupling g NN ν has to be introduced to make the nn → ppee amplitude finite and regulator-independent. The finite part of g NN ν , which encodes hard-neutrino exchange, can be determined by (i) matching the chiral EFT nn → ppee amplitude to future lattice QCD calculations; (ii) relating it via chiral symmetry to electromagnetic low-energy constants (LECs) that control isospin-breaking in the two-nucleon sector. A combination of couplings involving g NN ν can be fit to nucleon-nucleon charge-independence-breaking (CIB) observables, confirming the LO scaling of this coupling. Based on this, we speculate on the impact of g NN ν on nuclear matrix elements of relevance to experiments.
The need for an LO short-range ∆L = 2 interaction: We consider a scenario in which LNV at low energy is dominated by the electron-neutrino Majorana mass where C = iγ 2 γ 0 denotes the charge conjugation matrix. The nuclear effective Hamiltonian can be written as in terms of the Fermi constant G F and the V ud element of the CKM matrix [32,33]. The neutrino potential V ν can be obtained from two-nucleon irreducible diagrams mediating nn → ppee to a given order in p/Λ χ . Within Weinberg's power counting the only LO contribution [26] comes from the exchange of potential neutrinos, with q 0 |q|, where g A 1.27 is the nucleon axial coupling, m π the pion mass, and q the momentum transfer. N 2 LO terms arise from corrections to the single nucleon weak currents, irreducible one-loop diagrams, and contact interactions mediating ππ → ee, n → pπ + ee, and nn → ppee. In particular, the short-range potential includes a two-nucleon term [26] where the LEC g NN ν is O((4πF π ) −2 ) in Weinberg's counting and F π = 92.2 MeV is the pion decay constant. However, it is known that Weinberg's power counting leads to inconsistent results in nucleon-nucleon scattering [34][35][36][37] and nuclear processes mediated by external currents [38], due to a conflict between naive dimensional analysis and nonperturbative renormalization. We therefore investigate the scaling of g NN ν by studying the amplitude A(nn → ppee) ≡ A ∆L=2 with strong interactions, H strong , included nonperturbatively.
We work at LO in chiral EFT, and focus on the scattering of two neutrons to two protons in the 1 S 0 wave, where H strong has short-range and Yukawa components, [31,34,35]. We have checked that transitions involving higher partial waves such as 3 P 0,1 → 3 P 0,1 are correctly renormalized and do not require enhanced ∆L = 2 counterterms.
The contributions to A ∆L=2 from the exchange of a light neutrino (A (ν) ∆L=2 ) are shown in Fig. 1. The blue ellipse denotes the iteration of the Yukawa potential V π (q). The diagrams in the second and third rows include an infinite number of bubbles, dressed with iterations of V π . Without loss of generality for our arguments, we use the kinematics n(p) n(−p) → p(p ) p(−p ) e(p e1 = 0) e(p e2 = 0), with |p| = 1 MeV and correspondingly ∆L=2 can be expressed in terms of the Yukawa "in" and "out" wavefunctions χ ± p (r) and the propagators In the counterterm amplitude (fourth line) the black square represents g NN ν . The · · · in the second to fourth lines denote diagrams with arbitrary numbers of bubble insertions. [34,37]. Observing that the bubble diagrams in Fig. 1 are related to G + E (0, 0), while the triangles dressed by Yukawas are related to χ + p (0) and χ − p (0) * = χ + p (0) [34], the LO amplitude reads where A A , A B , and A C denote the first diagram in the first, second, and third rows of Fig. 1, respectively (without the wavefunctions at 0, in the case of A B and A C ). A B is similar to A B and not shown in Fig. 1.
To study the renormalization of the ∆L = 2 amplitude, we now discuss the divergence structure of A [34]. We note that: (i) All diagrams in A A are finite. The tree level is finite and each V π iteration improves the convergence by bringing in a factor of d 3 k/(k 2 ) 2 , where one k 2 comes from the pion propagator and the other from the twonucleon propagator.
(ii) All the diagrams in A B andĀ B are finite. The first loop goes as d 3 k/(k 2 ) 2 , while V π insertions further improve the convergence.
(iii) The first two-loop diagram in A C has a logarithmic divergence, which stems from an insertion of the most singular component of the neutrino potential, namelỹ The two-loop diagram with insertion of V ν,0 −Ṽ ν and higher-loop diagrams are convergent. We focus on A C and write A C = A (div) C + δA C . In dimensional regularization, The counterterm amplitude, shown in the fourth line of Fig. 1, reads and we can renormalize A ∆L=2 by replacing A C → A C + 2g NN ν /C 2 in Eq. (6). In the MS scheme, after defining the dimensionless coupling This coupling obeys the renormalization-group equation A similar enhancement also occurs in four-nucleon couplings induced by higher-dimensional LNV operators. Treating V π as a subleading correction [35,39] is equivalent to working to LO in pionless EFT, and does not affect our conclusions about the importance of g NN ν [26]. Details on how to obtain δA C will be provided in future work [40].
A ∆L=2 in a cutoff scheme: The need for an LO counterterm can be demonstrated also in a coordinatespace scheme that makes no direct reference to Feynman diagrams. In this approach we regulate the short-range part of V 0 with a smeared δ-function,  and obtain ψ − p (r) and ψ + p (r) by solving the Schrödinger equation. We determineC(R S ) by requiring that the 1 S 0 scattering length be reproduced (C ≈ −0.4/F 2 π at R S = 0.8 fm). We find that 1/C(R S ) has linear (1/R S ) and logarithmic divergences [35] and that the 1 S 0 phase shifts at nonzero momentum are indeed R S -independent.
We then compute where V ν,0 (r) is obtained by Fourier-transforming the 1 S 0 projection of Eq. (3). In Fig. 2 we plot A (ν) ∆L=2 as a function of R S . The plot displays a logarithmic dependence on R S (analogous to the log µ dependence in Eq. (10)) as well as milder power-like behavior. Therefore, to obtain a physical, regulator-independent amplitude one needs to include an LO counterterm, given in R S (r). The corresponding amplitude, is also regulator-dependent. As expected from Eq. (9), we find its leading divergent behavior to be well reproduced by 1/C(R S ) 2 . We can then make ∆L=2 finite for R S → 0 and R S -independent by choosingg NN ν (R S ) = −(a/2)(1 + 2g 2 A ) log R S + b + cR S + · · · , with the coefficient of the logarithm quite close to the MS expectation a = 1.
Relating g NN ν to electromagnetic isospin violation: The finite part of g NN ν can be obtained by matching the chiral EFT amplitude to a lattice QCD calculation performed at the same kinematic point, as it is done in the strong-interacting sector [41]. First lattice results related to double-beta decay are starting to appear [42,43].
We now discuss a complementary estimate based on the fact that the short-range operators and associated LECs arising in 0νββ and electromagnetic processes are closely related [26]. In the electromagnetic case, the short-range hadronic operators arise from amplitudes in the underlying theory involving two insertions of the electromagnetic current with exchange of hard virtual photons [44,45]. In the ∆L = 2 case, up to a proportionality factor, the same operators are generated by the insertion of two weak currents with exchange of hard neutrinos. This comes about because the neutrino propagator and weak vertices combine to give a massless gauge-boson propagator in Feynman gauge, multiplied by 8G 2 F V 2 ud m ββēL e c L [26]. The LECs needed for 0νββ are therefore related to the LECs associated with the isospin I = 2 component of the product of two electromagnetic currents, which belongs to the 5 L × 1 R irreducible representation of chiral SU (2) L × SU (2) R .
Only two independent four-nucleon operators that transform as I = 2 objects exist: where Q L = u † Q L u, Q R = uQ R u † , u = exp(iτ · π/(2F π )), and Q L,R are "spurions" transforming under the chiral group as Q L → LQ L L † , Q R → RQ R R † . In the electromagnetic case Q L = Q R = τ 3 /2, while in 0νββ Q L = τ + , Q R = 0. In our conventions O 1 enters the ∆L = 2 Lagrangian with coefficient 2G 2 F V 2 ud m ββ g NN ν . Defining the electromagnetic LECs multiplying O 1,2 as e 2 C 1,2 /4, chiral symmetry dictates g NN ν = C 1 . In the electromagnetic case, O 1 and O 2 only differ at the multipion level, and an isospin-breaking two-nucleon observable, such as the I = 2 combination of scattering lengths a CIB = (a nn + a pp )/2 − a np , only constrains the sum C 1 + C 2 . Extracting this combination from data provides a rough estimate of g NN ν under the assumption C 1 ∼ C 2 . As in the ∆L = 2 case, we introduce the dimensionless couplingsC i ≡ [4π/(m NC )] 2 C i and compute the scattering lengths a pp,nn,np including the leading sources of isospin breaking -the Coulomb potential and pion mass splitting -andC 1 +C 2 . Similarly to the ∆L = 2 case, we find thatC 1 +C 2 needs to be promoted to LO and obeys the RGE while, of course,C 1 has the same RGE asg NN ν . By fitting to a CIB using a np = −23.7 fm, a nn = −18.9 fm, and a pp = −7.8 fm, we find (C 1 +C 2 )/2 ≈ 2.5 at µ = m π in the MS scheme. Using instead the R S scheme, we find (C 1 +C 2 )/2 ≈ 2.0 at R S = 0.5 fm. 1 This estimate, based on data and chiral symmetry, again confirms that g NN ν ∼ O(F −2 π ). Numerical impact: To roughly estimate the impact of the contact term, we assume for concreteness C 1 = C 2 and hence g NN ν = (C 1 + C 2 )/2 at someR S orμ −1 in the range 0.002 − 0.8 fm. The total two-nucleon amplitude ∆L=2 then becomes independent of the regulator, as illustrated in Fig. 2, where the widths of the horizontal bands reflect the ambiguity in the choice of the pointR S orμ where C 1 = C 2 is assumed. (They do not account for the uncontrolled error of the assumption itself.) The relative size of the two components depends on R S , with A ∆L=2 ∼ 30% at R S ∼ 0.1 fm, decreasing to ∼ 10% at R S ∼ 0.6 fm. More insight can be obtained by plotting the matrix-element densities ρ ν and ρ NN defined as Figure 3 (top panel) shows that ρ NN (r) is concentrated at smaller distances than ρ ν (r), and its contribution to the amplitude is thus partially diluted.
We have performed a similar analysis for A = 6, 12 nuclei, using Variational Monte Carlo nuclear wavefunctions [51] based on the AV18 two-nucleon [50] and IL7 three-nucleon [52] interactions. The mismatch between the short-range behaviors of existing strong-interaction potentials and our 0νββ interaction introduces additional model dependence, which we mitigate by: (i) Considering an alternative extraction of (C 1 + C 2 )/2 from the phaseshift analysis of Refs. [47,48] 2 , which employs the same regulator (13) with R S 0.6 − 0.8 fm, approximately the range of AV18's short-range part. (ii) Simply replacing our V ν,CT (r) with AV18's short-range CIB potential.
For ∆I = 0 transitions such as the 6 He → 6 Be shown in Fig. 3 (middle panel), we find A ∆L=2 ∼ 10%, similarly to the nn → ppee case. In realistic 0νββ transitions, however, the total nuclear isospin changes by two units, ∆I = 2. This implies the presence of a node in ρ ν (r) due to the orthogonality of the initial and final spatial wavefunctions. The resulting partial cancellation between the regions with r < ∼ 2 fm and r > ∼ 2 fm [51] leads to a relative enhancement of the short-range contribution, as illustrated in Fig. 3 (bottom panel) for 12 Be → 12 C. Numerically we find A ∆L=2 ∼ 25% (our fit), 1 Our result is consistent with analyses based on chiral [46][47][48][49] and phenomenological potentials such as AV18 [50], which also find that, except at very low energies, long-and short-range components of the CIB interaction induce effects of similar size. 2 C 1 + C 2 is related to the CIB coefficient C IT 0 of Refs. [47,48] by (C 1 + C 2 )/2 = −6C IT 0 /e 2 .  3. ρν (r) and ρNN (r) for the nn → ppee process (top), and for nuclear transitions with A = 6 (middle) and A = 12 (bottom). In the middle and bottom panels the green (ρNN ) and red (ρ (χ) NN ) bands correspond to g N N ν = (C1 + C2)/2 extracted from our analysis and from Refs. [47,48], respectively. ∼ 55% (fit from Refs. [47,48]), and ∼ 60% (AV18 representation of the short-range CIB potential). Because the node in the density is a robust feature of ∆I = 2 transition [53,54], we expect the effects in 12 Be → 12 C and experimentally relevant transitions to be of comparable size.
Conclusion: The above arguments suggest that the new short-range ∆L = 2 potential identified in this work can significantly impact 0νββ phenomenology and its implications for Majorana neutrino masses. We hope this will stimulate work towards a more controlled determination of g NN ν from lattice QCD and an assessment of the impact of the short-range potential in nuclei of experimental interest.
Acknowledgments: VC, WD, MG, EM, and SP acknowledge support by the U.S. Department of Energy, Office of Science, and by the LDRD program at Los Alamos National Laboratory. VC and EM acknowledge partial support from the DOE topical collaboration on "Nuclear Theory for Double-Beta Decay and Fundamental Symmetries". JdV acknowledges support by the Dutch Organization for Scientific Research (NWO) through a VENI grant. The work of UvK was supported in part by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics, under award No. DE-FG02-04ER41338, and by the European Union Research and Innovation program Horizon 2020 under grant agreement No. 654002. We acknowledge stimulating discussions with Will Detmold and Martin Savage, which triggered this research. We thank Joe Carlson, Jon Engel, Amy Nicholson, Maria Piarulli, Petr Vogel, Andre Walker-Loud, and Bob Wiringa for discussions at various stages of this work.