Neutrinoless Double Beta Decay from Lattice QCD: The Short-Distance $\pi^-\rightarrow\pi^+ e^- e^-$ Amplitude

This work presents a determination of potential short-distance contributions to the unphysical $\pi^-\rightarrow\pi^+ e^- e^-$ decay through lattice QCD calculations. The hadronic contributions to the transition amplitude are described by the pion matrix elements of five Standard Model Effective Field Theory operators, which are computed on five ensembles of domain-wall fermions with $N_f = 2 + 1$ quark flavors with a range of heavier-than-physical values of the light quark masses. The matrix elements are extrapolated to the continuum, physical light-quark mass, and infinite volume limit using a functional form derived in chiral Effective Field Theory ($\chi\mathrm{EFT}$). This extrapolation also yields the relevant low-energy constants of $\chi\mathrm{EFT}$, which are necessary input for $\chi\mathrm{EFT}$ calculations of neutrinoless double beta decay of nuclei.


I. INTRODUCTION
Neutrinoless double beta (0νββ) decay, if observed, would unambiguously reveal the existence of physics beyond the Standard Model (BSM) [1]. In particular, it would imply that the difference between baryon number and lepton number (B − L) is not a fundamental symmetry of the universe [2], and would prove that the neutrino is a Majorana particle [3]. Moreover, observation of 0νββ decay would provide additional information about the matter-antimatter asymmetry in the universe [4], which may help to explain baryogenesis and further constrain the neutrino masses [5].
As such, experiments are underway worldwide to search for 0νββ decay, the most sensitive of which study 76 Ge and 136 Xe and constrain the half-lives of 0νββ decay in each isotope to be greater than 10 26 years [6][7][8][9][10]. Understanding the implication of these constraints for possible BSM physics scenarios requires input in the form of nuclear matrix elements (NMEs); which NMEs are relevant depends on the underlying mechanism of 0νββ decay. These mechanisms can be broadly divided into two categories: long-distance mechanisms, in which the decay is induced by a non-local interaction mediated by a light particle of mass much less than the hadronic scale [11,12]; and short-distance mechanisms, in which the decay is mediated by a heavy particle that can be integrated out in Effective Field Theory (EFT) to generate contact interactions [13,14]. In extensions of the Standard Model, long-distance mechanisms are typically assumed to be generated by the dimension-5 Weinberg operator, in which the mediating particle is generally a light Majorana neutrino (although other scenarios have been considered) [15][16][17][18], while short-distance mechanisms are described by operators of dimension greater than or equal to 9 [19]. The dominant mechanism of 0νββ decay will determine the scale Λ LNV at which leptonnumber violating physics is observed. In particular, if 0νββ decay is primarily described by a long-distance mechanism, then Λ LNV 1 TeV [20], while if 0νββ decay is primarily described by a short-distance mechanism, Λ LNV ∼ 1 TeV [21,22]. Both cases must be understood in order to draw conclusions about the underlying BSM physics from any experimental detection of 0νββ decay.
Calculations of long and short-distance 0νββ decay matrix elements have been performed with nuclear manybody methods [21,23]. These techniques are currently the only theoretical methods which can provide insight into 0νββ decay in nuclear isotopes which are experimentally relevant. The requisite NMEs for the long and short-distance 0νββ decay of 48 Ca, 76 Ge and 136 Xe have been computed, although large model dependence in the calculated NMEs remains a challenge for these techniques [6,7,[24][25][26]. To improve these calculations, connection to the Standard Model is required.
Lattice quantum chromodynamics (LQCD) is the only known method with which to compute NMEs directly from the underlying quark and gluon degrees of freedom. However, current LQCD calculations of nuclei suffer from a signal-to-noise problem [27,28] and a factorial increase in the number of quark contractions with atomic number [29], which make calculations of phenomologically relevant nuclei impractical in the absence of new algorithms and approaches. Instead of direct computation of large nuclei, recent work uses EFT [16,22,[30][31][32][33] to relate LQCD calculations of simpler processes such as the unphysical mesonic transition π − → π + e − e − and the twonucleon 0νββ decay n 0 n 0 → p + p + e − e − to nuclear 0νββ decay. Studies of the π − → π + e − e − transition in particular do not incur the technical challenges faced by LQCD calculations of nuclei. The long-distance pion matrix elements have been computed directly using LQCD with a domain-wall fermion action [34,35]. The associated short-distance pion matrix elements have been calculated from LQCD input with two approaches: relating the desired matrix elements to kaon-mixing matrix elements, assuming SU (3) chiral symmetry [36]; and computing the pion matrix elements directly using LQCD with a mixed action [37].

arXiv:2208.05322v2 [hep-lat] 16 Dec 2022
This work presents a direct LQCD computation of the π − → π + e − e − matrix elements of the leading shortdistance (dimension-9) operators, performed for m e = 0 and at threshold. This calculation uses domain-wall fermions, as their chiral symmetry properties yield matrix elements that have a simple renormalization structure. There is a mild tension between the results of the present calculation and the previous mixed-action LQCD calculation of the same matrix elements in Ref. [37], which may be due to the differences in the action used in each calculation. The ensembles used in this calculation are the same as those used in the first lattice computation of the long-distance π − → π + e − e − amplitude mediated by light Majorana neutrino exchange [34]. As such, both the long and short-distance contributions to π − → π + e − e − have now been computed in a consistent framework, allowing conclusions to be drawn regarding the relative importance of the two potential contributions, as discussed in Section IV.
The remainder of this paper is organized as follows. Section II details the EFT framework for the shortdistance π − → π + e − e − decay and the LQCD calculation of the hadronic part of the transition amplitude. Section III describes the procedure used to extrapolate the renormalized LQCD matrix elements to the physical point using a model based on chiral EFT (χEFT), and presents results for the extrapolated matrix elements and the extracted χEFT low-energy constants (LECs). Section IV summarizes the results and presents an outlook.

A. Short-distance operators
In the Standard Model EFT (SMEFT) framework, the Standard Model enters as the renormalizable sector of a non-renormalizable theory [38]. Potential short-distance contributions to π − → π + e − e − are induced by physics at the scale Λ LNV v, where v = 247 GeV is the electroweak scale set by the Higgs vacuum expectation value, and described in the SMEFT by operators with mass dimension greater than 4. At the quark level, any SMEFT operator that contributes to 0νββ decay must induce the process dd → uue − e − . Every such operator must therefore contain at least six fermion fields, and so have mass dimension d ≥ 9, with contributions to the π − → π + e − e − decay power-suppressed by a factor of Λ d−4 LNV . The dimension-9 lepton-number violating operators thus contribute to the decay at leading-order (LO) in inverse powers of Λ LNV .
There are fourteen SU (3) c × U (1) EM -invariant dimension-9 SMEFT operators which violate lepton number and may contribute to π − → π + e − e − ; they can be factorized into a 4-quark operator multiplying a leptonic operator. Of these operators, four have corresponding 4-quark operators that transform as Lorentz 4-vectors, and therefore match to the χEFT operator π(∂ µ π)eγ µ γ 5 e c + h.c., where the superscript c denotes charge conjugation and π and e represent the pion and electron fields.
Integration by parts shows that pionic matrix elements of this operator are proportional to one power of the electron mass and give subleading contributions to the decay π − → π + e − e − . Of the remaining ten operators, five have corresponding 4-quark operators with positive parity and contribute to π − → π + e − e − , while the five operators containing 4-quark operators of negative parity do not contribute. Consequently, at LO the decay is described with the Lagrangian where G F is the Fermi coupling constant, c k are dimensionless Wilson coefficients, and the operator basis with k ∈ {1, 2, 3, 1 , 2 } [31]. Here q L (x) and q R (x) are the left and right-handed components of the quark field isospin doublet, respectively, and is the isospin raising operator. The round and square brackets in Eq. (2) denote color contraction: for arbitrary Dirac matrices Γ 1 and Γ 2 , the operators (2) is named the BSM basis and is typically used in phenomenological calculations of 0νββ decay [33]. Although the π − → π + e − e − transition is unphysical, it has phenomenological importance as it can be related to the nuclear decays with χEFT [22]. In particular, the two-nucleon decay n 0 n 0 → p + p + e − e − is induced in χEFT by the diagrams in Fig. (1) and has LO contributions from the ππ and N N vertices [16,33]. 1 The associated effective Lagrangian relevant for π − → π + e − e −     Fig. (1a), are determined by O χ k in Eq. (4). The ππ (Fig. (1a)) and N N (Fig. (1b)) diagrams are the LO χEFT contributions to n 0 n 0 → p + p + e − e − .
(i.e., omitting N N and πN operators which do not contribute) is [37], Here, f π is the pion decay constant in the chiral limit, Λ 2 χ ≡ 8π 2 f 2 π is the scale of chiral symmetry breaking, and O χ k denote the leading χEFT operators corresponding to O k [39]. The χEFT LECs β k determine the ππ coupling, and are also essential input to study the nucleonic decay. The β k can be determined by evaluating the pion matrix elements of the O k in LQCD and matching them to the corresponding matrix elements of O χ k in Eq. (4).

B. Bare matrix elements
The pion matrix elements of each of the SMEFT operators in Eq. (2) are computed in LQCD using gauge-field ensembles with N f = 2 + 1 quark flavors generated by the RBC/UKQCD collaboration [42,43]. Each ensemble uses the Shamir kernel [44] for the domain-wall fermion action [45] and the Iwasaki action [46] for the gauge field. The parameters of each ensemble are detailed in Table I, and additional details regarding the ensemble generation can be found in Refs. [42,43,47]. The scale is set using the Wilson flow scale w 0 [40]. The pion mass, m π , the pion decay constant, f π , and the axial-vector renormalization constant, Z A , for each ensemble were determined in Ref. [34]. In the conventions used here, the physical pion decay constant [48] is f (phys) π = 130.2 MeV. The vector renormalization constant, Z V , for these ensembles was computed in the chiral limit in Refs. [40,41]. Because Z V ≈ Z A , the ensembles exhibit approximate chiral symmetry.
On each ensemble, the time-averaged two-point function and three-point functions where the pion interpolating operator χ π (x) = u(x)γ 5 d(x) has the quantum numbers of the π − and t + ≥ t x ≥ t − , are computed for each operator O k (x) in the BSM basis (Eq. (2)). Wall-source propagators are computed at each available time slice on each configuration, where "wall" denotes projection to vanishing three-momentum in the Coulomb gauge. Note that wall sources are not gauge-invariant, hence the need for gauge fixing. The two-point functions (Eq. (5)) are constructed using a wall source propagator at t − and a wall sink at t + t − , and the three-point functions (Eq. (6)) are constructed using wall source propagators at t − and t + and a point (local) sink at t x . The explicit Wick contractions are given in Appendix A.
The bare pion matrix elements in lattice units Subtracting 1 2 C 2pt (T /2)e mπ(2t−T /2) in the denominator of Eq. (8) isolates the backwards-propagating state in the two-point function, and in the 0 t T limit TABLE I. Parameters of the gauge field ensembles used in this study. Each ensemble was generated with two degenerate light quark flavors of mass m and one heavy quark flavor of mass ms. The lattice volumes are L 3 × T × Ls, with the fifth dimension having Ls sites. Derived quantities are computed in Ref. [34] (the pion mass mπ, the pion decay constant fπ, and the axial current renormalization ZA) and Refs. [40,41] (the vector current renormalization ZV and the inverse lattice spacing a −1 ).
O eff k (t) asymptotes to O k . The effective matrix elements are computed on between 33 and 53 gauge field configurations for each ensemble (details in Appendix B, Table III), resampled using a bootstrap procedure with n b = 50 bootstrap samples. The spectral decomposition of O eff k (t) up to and including the first excited state with energy m π + ∆, This function is used to model the temporal dependence as free parameters.
Fits of O eff k (t) to the model of Eq. (10) are performed using a correlated least-squares fit. Each fit is performed over a given range [t min , t max ], with the covariance matrix obtained from the bootstrapped sample covariance matrix via linear shrinkage with parameter λ [49,50]; the hyperparameters are varied, with t min ∈ [6,11], t max ∈ [30,32], and λ ∈ {0.1, 0.2, 0.3, 0.4}. Bayesian priors are placed on the model parameters, informed by the results of a two-state fit to C 2pt (t). The priors on the spectral coefficients are set to A (i) where µ ± σ denotes the normal distribution with mean µ and width σ. To enforce positivity, log-normal priors are chosen for the mass m (k) π and excited state gap ∆ (k) such that m (k) = m π ± δm π , where m π (δm π ) is the mean (standard deviation) of the pion mass (Table I), and ∆ (k) = 2m π ± m π . Statistically indistinguishable results are obtained for O k under variation of all hyperparameters within the ranges described above, and when widths of the priors are inflated by a factor of 2,

C. Renormalization
To make contact with phenomological calculations, lattice-regulated matrix elements must be renormalized in the MS scheme. In this calculation, the renormalization coefficients are computed non-perturbatively in the RI/sMOM-(γ µ , γ µ ) (abbreviated as RIγ) scheme [51,52] and perturbatively matched to MS. In terms of the operator basis {O k (x)} (Eq. (2)), the renormalized matrix elements can be expressed as where sums over repeated indices are implied. Here O (x; a) denotes the bare operator at lattice spacing a, and is the multiplicative matching coefficient from the RIγ to MS schemes, computed at one-loop in perturbation theory in the strong coupling α s (µ) [52,53]. Note that each renormalization coefficient is mass-independent and defined in the chiral limit.
The renormalization coefficients, Eq. (11), are conventionally computed in the Non-Perturbative Renormalization (NPR) operator basis, {Q n (x)}, which contains dif- ferent linear combinations of operators than the BSM basis of Eq. (2). Correlation functions involving the color-mixed operators O 1 (x), O 2 (x) may be rewritten with Fierz identities [54] as combinations of color-unmixed quark bilinears, which simplifies the calculation. The NPR basis is defined in terms of the quark bilinears: This basis is related to the positive-parity projection of the BSM basis, Eq. (2), as The space spanned by {Q n (x)} splits into three irreducible subspaces under chiral symmetry, with bases As both the MS and RIγ schemes obey chiral symmetry, the renormalization coefficients Z MS;Q nm (µ 2 ; a) and Z RIγ;Q nm (µ 2 ; a), which satisfy analogous equations to Eqs. (11) and (12), each factorize into a direct sum of three block diagonal matrices, each of which spans an irreducible subspace.
To renormalize the NPR basis operators, the four-point functions are computed on each ensemble, where V = L 3 × T is the lattice volume and q = p 2 − p 1 . Latin letters a, b, c, d denote color indices, while Greek letters α, β, γ, δ denote Dirac indices. All correlation functions used for the renormalization are computed in the Landau gauge with momentum sources [55] using 10 configurations for each ensemble, as the V 2 averaging from the momentum sources significantly reduces noise. The momenta are chosen subject to the symmetric constraint with the particular choice with q = p 2 − p 1 and j ∈ Z. The kinematic configuration corresponding to G n (q; a, m ) is depicted in Fig. (3). Note that with this choice of momentum, each value of q corresponds to a unique value of p 1 and p 2 , hence functions of (p 1 , p 2 , q) are labeled as functions of q for conciseness. The four-point functions are amputated, where is the Landau-gauge momentum-projected quark propagator. The ensemble dependence of Λ n (q), G n (q), and S(p) has been suppressed in Eq. (19) for clarity. Projectors (P n ) βαδγ badc are introduced to project (Λ m ) αβγδ abcd onto the NPR basis for RIγ [52] to yield a matrix of projected four-point functions with components F mn (q; a, m ) ≡ (P n ) βαδγ badc (Λ m ) αβγδ abcd (q; a, m ).
The remaining quantities which are computed nonperturbatively on each ensemble are the RIγ quark-field renormalization and the vector and axial-vector-renormalization coefficients, Z V (µ 2 ; a, m ) and Z A (µ 2 ; a, m ), whose computation is described in Appendix C.
is the vector three-point function, with V µ (x) = u(x)γ µ d(x) the vector-current. The quantities Z ∈ {Z RIγ q /Z V , F nm } display mild dependence on quark mass, and are extrapolated to the chiral limit via a joint fit over ensembles with different masses to the model Z(µ 2 ; a, m ) = Z(µ 2 ; a) +Z(µ 2 ; a)m (24) where Z(µ 2 ; a) andZ(µ 2 ; a) are fit coefficients, and Z(µ 2 ; a) is understood as the chiral limit of Z(µ 2 ; a, m ). Correlations between Z RIγ q /Z V and F nm on each ensemble are retained in the fits, and the covariance matrix is block-diagonal as data from different ensembles is uncorrelated. Fitted values of Z(µ 2 ; a) are statistically consistent when a constant model Z(µ 2 ; a, m ) = Z(µ 2 ; a) is used in place of the linear model of Eq. (24). The full set of extrapolations for (Z RIγ q /Z V )(µ 2 ; a) and F mn (q; a) for both the a = 0.11 fm and a = 0.08 fm ensembles is shown in Appendix D.
With the definitions above, the NPR-basis renormalization coefficients in the RIγ scheme can be computed as is the matrix of projections of the tree-level vertex function Λ (tree) n , and the notation | sym denotes evaluation at the symmetric kinematic point, Eq. (17). The renormalization coefficients Z RIγ;Q nm (µ 2 ; a)/Z 2 V are only computed non-perturbatively at scales µ j = 2π aL ||(j, j, 0, 0)|| corresponding to the lattice momenta given in Eq. (18), where || · || denotes the Euclidean norm of the lattice vector. However, the matching coefficients C MS←RIγ;Q nm (µ 2 , a) in Eq. (11) have been computed at µ = M ≡ 3 GeV [52,53], and therefore the renormalization coefficients must be perturbatively evolved from µ j to M . To minimize the artifacts from truncating the perturbative expansion of the matching coefficients, µ j must be chosen to lie in the Rome-Southampton window [56,57], with µ j taken to satisfy µ j ≤ M to minimize discretization artifacts. In practice, the scale µ 4 is used for renormalization at both a = 0.11 fm and a = 0.08 fm, as this is the nearest available scale to M satisfying these constraints. Numerically, these scales are µ 4 = 2.64 GeV for the a = 0.11 fm ensemble and µ 4 = 2.65 GeV for the a = 0.08 fm ensemble. Scale evolution from µ 4 to M is performed by integrating the evolution equation, where the NPR basis anomalous dimensions γ RIγ;Q nm (α s (µ)) have been computed at two-loop order in α s (µ) in Ref. [58].
The components corresponding to transitions between operators in different irreducible chiral representations are consistent with |Z MS;Q nm /Z 2 V | < 10 −5 and thus set to zero in Eq. (28). The renormalization coefficients have been computed for the NPR operator basis (Eq. (14)) in Ref. [52] using s quarks in place of d quarks. The results in Ref. [52] agree with Eq. (28) at the percent level, and deviations between the results are likely due to perturbative truncation errors, as Ref. [52] used non-pertubative step-scaling [56,57]. The NPR basis renormalization coefficients are converted to the BSM basis using the change of basis matrix, Eq. (15), and combined with the bare matrix elements to form renormalized matrix elements, O k (m π , f π , a, L) ≡ π + |O MS k (p = 0)|π − (m π , f π , a, L). (29) On a given ensemble, the renormalization coefficients and bare matrix elements are computed on different configurations, as the former are only computed on a subset of 10 of the configurations used to compute the matrix elements on each ensemble. As such, they are combined as an uncorrelated product and their errors are added in quadrature. The renormalized matrix elements are shown in Table II.

III. CHIRAL EXTRAPOLATION
The renormalized matrix elements O k (m π , f π , a, L), Eq. (29), computed on each ensemble, are extrapolated to the continuum and infinite volume limit and physical pion mass using χEFT at N 2 LO; the relevant expressions have been derived in Ref. [37] using the Lagrangian in Eq. (4). The chiral models F k for O k are given by where 2 π = m 2 π /Λ 2 χ is a power-counting parameter for χEFT, β k are the LO LECs defined in Eq. (4), and α k and c k are the additional NLO LECs. The matrix elements O 1 and O 2 have the same chiral behavior as O 1 and O 2 and are modeled by F 1 and F 2 , respectively, but with different LECs, α 1 , β 1 , c 1 and α 2 , β 2 , c 2 . The  4. Chiral extrapolation of renormalized matrix elements. The LQCD results are shown at 2 π = m 2 π /(8π 2 f 2 π ) calculated using the pion mass of each ensemble and the physical value of fπ, and the values of O k (mπ, fπ, a, L) have been shifted by −F k (mπ, fπ, a, L; α k , β k , c k ) + F k (mπ, f (phys) π , 0, ∞; α k , β k , c k ), where α k , β k , c k are the best-fit coefficients given in Table II. The physical pion mass is denoted by the dashed line.
are sums of modified Bessel functions K i (z) arising from one-loop, finite volume χEFT in the p-regime.
The models are fit to the data in Table II, using leastsquares minimization including the correlations between O k , m π , and f π on each ensemble. The final extrapolated results for the matrix elements and corresponding LECs are given in Table II. The resulting fits are shown in Fig. (4), where to isolate the pion-mass dependence of the matrix elements, 2 π has been rescaled by (f (lat) π /f (phys) π ) 2 and the values of O k (m π , f π , a, L) have been shifted by −F k (m π , f π , a, L; α k , β k , c k ) + F k (m π , f (phys) π , 0, ∞; α k , β k , c k ), where α k , β k , c k are the best-fit coefficients given in Table II. The extrapolation bands for each O k depict the functional form F k (m π , f (phys) π , 0, ∞; α k , β k , c k ). The results for π + |O MS k |π − obey the same hierarchy as the chiral SU (3) estimates [36], and are consistent with these results within two standard deviations.
The results for the renormalized, extrapolated, matrix elements are found to be in mild tension with the results of Ref. [37]. There are a number of differences between the two calculations which may account for the discrepancy. The present calculation was performed with the same domain-wall action for the valence and seaquarks and is thus unitary, while that of Ref. [37] used a mixed action where unitarity is only restored in the continuum limit. Using the domain-wall action for valence and sea-quarks yields matrix elements that have a mild dependence on the lattice spacing. In contrast, the mixed action results appear to have a larger dependence on the lattice spacing. However, the analysis of Ref. [37] was performed on nine ensembles with pion masses m π 310 MeV, including one ensemble with pion mass below the physical point, which allows for an interpolation to the physical point. Ref. [37] also uses three lattice spacings as opposed to the two used in this computation, which allows for higher control of discretization artifacts in the non-perturbative renormalization and the chiral and continuum extrapolation.

IV. CONCLUSION
This work presents a determination of the renormalized matrix elements and χEFT LECs for the short-distance operators that potentially arise from BSM physics at high scales and are relevant for the π − → π + e − e − transition. The present calculation is the first to use chiral fermions with the same valence and sea-quark actions. The domain-wall action yields a simple renormalization coefficient structure and straightforward extrapolation to the continuum and infinite volume limit and physical value of the light quark mass. With the results of Ref. [34], this completes the calculation of both the long and short-distance amplitudes for π − → π + e − e − on the same gauge-field ensembles.
One may compare the relative size of the decay amplitude of π − → π + e − e − induced by short-distance mechanisms, A SD , to that induced by long-distance mechanisms, A LD . In any model with a seesaw-type mechanism [59], for example the minimal left-right symmetric model [22], the effective Majorana neutrino mass m ββ scales as c/(G F Λ LNV ), where c is a Wilson coefficient.

This implies
where M 0ν is the long-distance nuclear matrix element for π − → π + e − e − . The final line of Eq. (32) arises by assuming that in a given BSM model, the dimensionless Wilson coefficients, c k and c, describing each amplitude are order 1, and by using dimensional arguments to approximate the matrix elements. In particular, the longdistance nuclear matrix element includes the convolution of a massless bosonic propagator with a bilocal QCD matrix element. The convolution picks out the dimensional scale 1/Λ 2 QCD , thereby enhancing the long-distance contribution compared to the short-distance one.
Since π + |O k |π − and M 0ν have now been computed consistently in LQCD, it is possible to compute the ratio of Eq. (32), quantitatively, given the Wilson coefficients c k and c from some model. For example, taking c k = c = 1, and using the LQCD results from this work and of Ref. [34] for the matrix elements yields ASD ALD = 6.1(2) × 10 −5 , consistent with expectations.
In addition to the pion-pion χEFT LECs, the other LECs contributing to nuclear 0νββ decay must be determined in future calculations in order to constrain models of new physics from experimental constraints on nuclear 0νββ decay rates. Knowledge of these LECs may be used as input for models of nuclear many-body physics, which may be used to estimate the half-lives of various nuclear 0νββ decay processes from short-distance mechanisms with increasing precision. The other LO LECs that are necessary for describing nuclear 0νββ decay are from the nucleon-nucleon interaction (Fig. (1b)), and may be determined with knowledge of the p + p + |O k (p = 0)|n 0 n 0 matrix elements [33]. Calculations of these matrix elements are ongoing and will provide the first direct LQCD probe of 0νββ decay in nuclear systems.
tion. These computations used the CPS [60], GLU [61], Grid [62], Chroma [63], QLUA [64], and QUDA [65] software packages. Least-squares fits were performed with the lsqfit software package [66]. Furthermore, the Tikz-Feynman package [67] was used to generate diagrams for this manuscript, and the RunDec3 package [68] was used to run the strong coupling in the perturbative match-ing for the renormalization. This work is supported in part by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics under grant Contract Number DE-SC0011090. WD is also supported by the SciDAC5 award DE-SC0023116. PES is additionally supported by the National Science Foundation under EAGER grant 2035015, and by the U.S. DOE Early Career Award DE-SC0021006.

A. THREE-POINT CONTRACTIONS
The correlation functions of Eq. (6) can be written in terms of the following contraction structures, where Γ 1 , Γ 2 are arbitrary Dirac matrices, Tr C (Tr D ) denotes a color (spin) trace, Tr = Tr C • Tr D denotes a full trace, and x = (x, t x ). Propagators S(t src → x) are computed with a zero three-momentum wall source at time t src ∈ {t − , t + } and a point sink at time t x , With the definitions of Eq. (33), the correlation functions are evaluated as where S = 1, P = γ 5 , V = γ µ , and A = γ µ γ 5 .

B. EFFECTIVE MATRIX ELEMENT FITS
Figs. (5)-(8) display the remaining fits to the effective matrix elements (Eq. (8)) that were not depicted in Fig. (2). The fit procedure is described in Section II B of the main text. The number of gauge field configurations per ensemble used in each matrix element extraction, n cfgs , and the corresponding bare matrix elements in lattice units, Eq. (7), are shown in Table III.

C. VECTOR AND AXIAL-VECTOR RENORMALIZATION COEFFICIENTS
Calculation of the scale and scheme-independent vector and axial-vector-current renormalization coefficients Z j (a), with j ∈ {V, A}, proceeds through the vector (Eq. (23)) and axial-vector three-point functions, where A µ (x) = u(x)γ µ γ 5 d(x). The momenta p 1 , p 2 , and q are subject to the symmetric constraint, Eq. (17), and parameterized identically to the modes used in the calculation of the four-quark operator renormalizations (Eq. (18)) with k ∈ {2, 3, 4, 5}. The lattice spacing dependence is made explicit in this section. The amputated three-point functions with j ∈ {V, A}, are used to compute the renormalization coefficients, wherep µ = 2 a sin( a 2 p µ ) is the lattice momentum. Note that the quark-field renormalization in Eq. (38) is defined in the RI/sMOM scheme [51], which differs from the RIγ scheme [52] of Eq. (22); Z V and Z A are scheme-independent, hence may be computed in any scheme. The chiral limits Z V (µ 2 ; a) and Z A (µ 2 ; a) of Z V (µ 2 ; a, m ) and Z A (µ 2 ; a, m ) are evaluated by a joint, correlated linear extrapolation of {Z RI/sMOM q , Z V , Z A } in m , identical to the procedure used in the am → 0 extrapolation of {Z RIγ q /Z V , F nm }, as described in Section II C of the text (Eqs. (21)-(24)). Although the renormalization coefficients Z V , Z A are scale-independent, the RI procedure introduces scaledependence from the kinematic setup (Eq. (17)). This scale-dependence is removed by fitting Z j (µ 2 ; a) to a power series in µ 2 and taking the µ 2 → 0 limit as described in Ref. [69], with fit model: Here Z j (a), c j (a), and c (2) j (a) are coefficients which are determined by correlated χ 2 minimization. The fits are shown in Fig. (9). The fits have χ 2 /dof ranging between 0.15 and 0.71. The best-fit value of Z j (a) is the value that is taken for the renormalization factor, and it is determined that Z V (0.11 fm) = 0.7119 (20) Z V (0.08 fm) = 0.7472(24) Z A (0.11 fm) = 0.7137 (19) Z A (0.08 fm) = 0.7462 (23).
The results show that Z V = Z A within statistical precision as expected. The determination presented in this work is consistent with the determination of Z V in Ref. [41] for the a = 0.08 fm and a = 0.11 fm ensembles, and with Z A in Ref. [34] for the a = 0.08 fm ensembles, although Z A differs from the a = 0.11 fm value in that work by about one standard deviation. This deviation may be due to discrepancies in the procedure used to extract Z A , as the fit model (Eq. (40)) does not capture all the discretization artifacts present in the data. /Z V and F nm , as described in Section II C of the text. Each renormalization coefficient is evaluated at q = 2π L (4, 4, 0, 0), which is the lattice momentum corresponding to the scale µ = µ 4 . In each of Figs. (10)- (17), the µ dependence of (Z RIγ q /Z V )(µ 2 ; a) and the q dependence of F nm (q; a) has been suppressed for clarity. The data is observed to have very mild dependence on am .  22), computed on the a = 0.11 fm ensembles at q = 2π aL (4, 4, 0, 0) and extrapolated to the chiral limit via a joint correlated linear extrapolation in am (Eq. (24)). The data is depicted in red, and the shaded band denotes the extrapolation.